1 Background and motivation

Sharp growth in e-commerce sales in developed countries worldwide over the last decade has led to substantial increases in urban freight goods transport, with associated negative impacts on road traffic, availability of kerbside space and air quality. In the UK, total measured national volumes in the parcels market increased by 10%, to 2.6 billion items, in 2018–19 (Ofcom 2019). This growth and worsening road traffic and parking conditions in city centres make parcel deliveries ever-more difficult to perform. Factors that increase vehicle use for parcel deliveries include the growing trend for same-day and ‘instant’ deliveries (within 2 h) leading to fragmentation of consignments and increasing the number and frequency of deliveries (Dablanc et al. 2017) and the substantial number of competing freight transport operators that result in much duplication of van activity (Browne et al. 2014). Last-mile parcel delivery operations, characterised by multi-player, multi-drop vehicle rounds where kerbside access is needed are at direct odds with an infrastructure designed and legislated in favour of passenger transportation (Allen and Browne 2014).

Many of our major cities and particularly London have seen a considerable shift to walking, cycling and the use of buses over the past 20 years together with a fall in car traffic (Transport for London 2016). As a result, road space is being increasingly reallocated in favour of dedicated cycle and bus lanes, as well as pavement widening programmes, which in central London has led to a 30% decrease in road network capacity for private motorised vehicles between 1993 and 2009 (Transport for London 2013). With declining kerbside stopping locations, carriers’ round efficiency declines due to the need for additional, unproductive driving whilst searching for parking locations and the potential increase in fines and general traffic disruption through illegal parking (Bates et al. 2017). With average speeds in cities also falling, for example, by 13% (from 16.6 to 14.5 mph) on local authority managed roads in London between 2015 and 2018 (Department for Transport 2020), there is growing interest in alternative methods for addressing the last-mile problem that aim to reduce reliance on road vehicles.

The study by Allen et al. (2018) on last-mile goods vehicle activity indicated that multi-drop parcel delivery drivers in central London typically walk 8 km (5 miles) while their vehicle is parked at the kerbside, taking up more than 60% of their time worked. Their observations suggest a two-echelon last-mile distribution model, with driving and walking as the higher and lower-level echelons, respectively. However, the design of such a last-mile distribution model gives rise to challenging optimisation problems, even for a single driver, in which decisions concerning the partitioning of customers that will be served either by driving or by walking, the sequence of locations to be visited in either mode, and the selection of parking points from where walking tours will start from and end at, need to be made simultaneously. In addition, the driving and walking has to be differentiated on the basis of travel time or cost between a pair of locations, which is likely to vary between these two modes.

Two-level distribution systems within the broader context of city logistics give rise to what is known as the class of two-echelon vehicle routing problems (2EVRPs) (see Cuda et al. 2015, for a survey). Whilst the definition of the 2EVRP allows for use of different types of vehicles at each echelon (e.g. Hemmelmayr et al. 2012), they are primarily differentiated by the capacities of the vehicles, but not by travel cost or time. The 2EVRP assumes the use of a set of ‘satellites’ at given locations, and where the first and the second echelon routes meet. The model described by Anderluh et al. (2017) that combines vans and cargo bikes defines a separate travel cost for each vehicle travelling between a pair of locations, but requires vans and cargo bikes to travel from separate depots, and assumes that the customers are already partitioned into those to be served by cargo bikes and vans. Our problem uniquely differs from the 2EVRP in that we associate different travel costs for the two different modes, and do not require the use of satellites.

There has recently been a stream of research dedicated to the combined use of walking and driving for home on-site services. We will not attempt to review this body of literature as it relates more to health care, but instead mention to the particular problem described by Fikar and Hirsch (2015) to serve as an example. Motivated by home health care services, but also applicable to others in home repair, maintenance and private tutoring, this problem concerns the use of vehicles to transport nurses, who can be dropped-off and picked-up at various locations. The nurses can choose to walk between a drop-off and pick-up location for site visits. Additional practical restrictions, such as time windows and mandatory working time and break regulations, also apply. The authors describe a heuristic algorithm and apply it on real-world data provided by the Austrian Red Cross.

More relevant to our setting in freight distribution is the truck-and-trailer routing problem (TTRP) (Derigs et al. 2013; Villegas et al. 2010) where the ‘truck-and-trailer’ route in the TTRP corresponds to the driving route in our setting, while the ‘truck only’ route corresponds to that of walking. Variants of the TTRP that take into account the time requirements on the deliveries as well as the limited carrying capacity of both modes of delivery have also been studied (Lin et al. 2011), which give rise to difficult optimization problems that have only recently been optimally solved with up to 100 customers (Parragh and Cordeau 2017). Such problems however implicitly assume the use of a single mode of transport in both layers. One study that looks at such a two-echelon distribution model where walking and driving is treated separately is by Nguyên et al. (2019) where time window constraints are also present. A similar distribution problem involving delivery to clusters of customers by a fleet of trucks operating from a single depot appears in De Grancy and Reimann (2016), where each customer is served on foot from the parking location within the cluster it belongs to, and where there are time window restrictions on the deliveries made. However the problems described by Nguyên et al. (2019) and De Grancy and Reimann (2016) both assume that delivery locations are clustered in advance and is provided as input data.

The paper by Lin (2011) studies a similar problem to the one presented here, namely a pick-up and delivery problem with two modes of delivery (e.g., a van and a foot courier) and time windows, in which the two modes are treated separately. In this problem, the upper-echelon mode of transport has capacity restrictions whereas the lower-echelon does not, which is different to our setting. A similar problem arising in the routing and scheduling of technicians of an energy provider is described by Coindreau et al. (2019), where workers are allowed to either share a vehicle or walk to get from one site to another. Each vehicle has a driver and can accommodate a worker, who can both choose to walk to a site. After any walking, the driver must return to the vehicle. If a worker is dropped-off at a particular location by one driver, they can be picked-up at another location by another driver, where the travel between the two locations is undertaken by foot by that worker. For this problem, the authors describe a formulation and a heuristic algorithm that is applied to a set of instances.

The delivery model we study here also bears resemblance to those that involve simultaneous routing of ground vehicles and drones, such as the one described in Luo et al. (2017) in which drones are used for deliveries and are carried on the ground vehicles. However, as the carriage capacity of drones is often limited to one parcel (Murray and Chu 2015), such problems assume that drones serve only one customer at a time, and return to the ground vehicle carrying it to pick up a new parcel for delivery to another customer. In addition, drones are also limited by their flight range. For this reason, the problem we study here can be regarded as a more general version of those involving drones and vehicles.

This paper contributes to the literature by (1) introducing a two-echelon last-mile delivery system that explicitly considers driving and walking decisions in an integrated manner, (2) describing an associated optimisation model, valid inequalities and a branch-and-cut algorithm to solve the model, and (3) presenting results on real instances derived from last-mile delivery operations in London. The rest of the paper is structured as follows. We present a formal description of the problem in Sect. 2. We then describe an optimisation model in Sect. 3, where we also present valid inequalities and describe the branch-and-cut algorithm. Computational results are given in Sect. 4. The paper concludes in Sect. 5.

2 Planning a combined driving and walking distribution model

The proposed distribution system requires explicit differentiation between walking and driving. In particular, the time required for travelling between two locations will differ between walking and driving, sometimes by a significant margin, particularly in urban areas. Furthermore, the dependency between the two echelons, one corresponding to driving and the other to walking, suggest that designing the two sets of routes separately may result in a sub-optimal solution. For a given area to be serviced by a single driver, the main decisions that need to be made in the design of a combined distribution system include:

  • The locations that are to be visited by either mode of transport,

  • The selection of locations where the vehicle will be parked,

  • The ordering of delivery locations in both the driving route and the walking routes, bearing in mind the capacity restrictions in the latter (e.g., associated with using a backpack, a wheeled bag or a trolley).

The distribution design problem concerns a single vehicle that is pre-loaded with parcels and destined to serve a given set of customers with known locations and delivery requirements specified in terms of number and size of items. The aim is to (1) group the customers into clusters (where singleton clusters are allowed), (2) find a driving route across the clusters such that one node within each cluster is designated as the parking node, and (3) find a walking route within each cluster that starts and ends at the parking node. The objective is to minimise the total time of all travel, under the following constraints: (1) the demands of all customers must be satisfied, (2) the total weight and volume of all deliveries made within each cluster is limited to the total capacity that the driver can carry, and (3) the driving starts from and ends at a depot of a known location. Figure 1left shows an instance to the problem with the depot shown by the dark square and 10 customers shown by the small circles. A feasible solution to the problem is shown on the right of Fig. 1 which indicates the driving by the solid lines and the walking by the dashed lines.

Fig. 1
figure 1

A feasible solution to an instance with one depot and 10 customers, where walking is shown by the dashed lines and driving by the solid lines

Full size image

The design of the distribution system that includes driving and walking therefore hinges on the trade-off between the time taken by the two modes, which requires the two modes to be treated separately. The model we describe in the following section explicitly addresses this particular aspect of the problem.

Time window constraints for parcel delivery are not explicitly considered here. Our analysis of a data set from a major carrier operating in London covering the period 4–9 June 2018 revealed that, of about 13,000 consignments that were delivered, only 1.6% required delivery by 9am, 2.3% by 10am and 7.2% by 12pm, whereas the rest (88.9%) could be delivered anytime until 6pm on a given day. Given the relatively small proportion of time-sensitive parcels, we assume that these could be delivered separately to ensure the delivery deadlines are met, and therefore focus our attention on parcels that do not have any time delivery constraints attached to them.

3 Mathematical model and branch-and-cut algorithm

This section presents a formulation of the problem described above along with valid inequalities. We also describe a branch-and cut algorithm to solve the formulation and in which the new valid inequalities are used.

3.1 Mathematical modelling

Let \(G=(V,A)\) be a complete graph with \(V = V'\cup \{0\}\) as the set of nodes, where node 0 corresponds to the depot, \(V'\) is the set of customers, and A is the arc set. Two different (non-negative) travel times are associated with each arc \((i,j) \in A\) that represent travel from node \(i \in V\) to \(j \in V\), namely the time \(c^d_{ij}\) of driving and the time \(c^w_{ij}\), of walking. We allow the costs to be asymmetric. We denote by \(v_i\) (\(w_i\)) the volume (weight) of the parcels to be delivered to customer \(i\in V'\).

There are no capacity or volume constraints on the vehicle as it is pre-loaded at the time of departure, but each walking route is limited to carrying parcels that total no more than \(Q_v\) units in volume and \(Q_w\) units in weight. In addition, any individual customer with demand larger than \(Q_v\) in volume or heavier than \(Q_w\) in weight must be visited and served by the vehicle. Such a customer cannot appear on a walking route, and we assume that they also cannot be a parking node. Let r(S) denote the number of walking routes needed to meet the demand in a given subset S of \(V'\), calculated as \(r(S) = \max \{\lceil v(S)/Q_v \rceil , \lceil w(S) / Q_w \rceil \}\) where \(v(S) = \sum _{i \in S} v_i\) is the total volume of parcels to be delivered to S, and \(w(S) = \sum _{i \in S} w_i\) is the total weight of parcels to be delivered to S.

The mathematical model is defined with respect to a linear multiplier \(\alpha \ge 0\) to the driving time, pre-defined as an input to the optimisation problem. It inversely reflects the relative importance given to walking over driving. In particular, increasing \(\alpha \) places less emphasis on minimising the total walking time, and puts more weight on reducing the overall driving time. This may be required if there are uncertainties on the driving times, or for when the overall vehicle mileage (and the associated fuel consumption) of the vehicle needs to be reduced.

The model uses a binary variable \(x^{d}_{ij}\), which takes the value 1 if driving takes place from node \(i \in V\) to node \(j \in V \setminus \{i\}\), and 0 otherwise. Similarly, a binary variable \(x^{w}_{ij}\) takes the value 1 if walking takes place from node \(i \in V\) to node \(j \in V \setminus \{i\}\), and 0 otherwise. Finally, a binary variable \(z_i\) takes the value 1 if the driver parks the vehicle at the location of customer \(i\in V'\) (hereafter called the parking node), and 0 otherwise.

Table 1 lists the parameters and the variables used in the model.

Table 1 A list of the parameters and the variables used in the model
Full size table

The complete formulation is given as follows:

$$\begin{aligned} \hbox {Minimise} \alpha \sum _{(i,j) \in A} c^d_{ij} x^{d}_{ij} + (1-\alpha ) \sum _{(i,j) \in A} c^w_{ij} x^{w}_{ij} \end{aligned}$$
(1)

subject to

$$\begin{aligned} \sum _{j \in V} x^d_{ij} + \sum _{j \in V'} x^w_{ij}&= 1+z_i&\quad \forall i \in V' \end{aligned}$$
(2)
$$\begin{aligned} \sum _{i \in V} x^{d}_{0i}&= 1&\quad \end{aligned}$$
(3)
$$\begin{aligned} \sum _{i \in V} x^{d}_{i0}&= 1&\quad \end{aligned}$$
(4)
$$\begin{aligned} \sum _{i,j \in S} x^d_{ij}&\le |S|-1&\quad \forall S \subset V' \end{aligned}$$
(5)
$$\begin{aligned} \sum _{j \in S} x^d_{ij}&\ge z_i&\quad \forall i \in V' \end{aligned}$$
(6)
$$\begin{aligned} \sum _{j \in V} x^d_{ij}&= \sum _{j \in V} x^d_{ji}&\quad \forall i \in V' \end{aligned}$$
(7)
$$\begin{aligned} \sum _{j \in S} x^w_{ij}&\ge z_i&\quad \forall i \in V' \end{aligned}$$
(8)
$$\begin{aligned} \sum _{j \in V'} x^w_{ij}&= \sum _{j \in V'} x^w_{ji}&\quad \forall i \in V' \end{aligned}$$
(9)
$$\begin{aligned} x^w_{ij} = x^w_{ji}&= 0&\quad \forall i \in V' \end{aligned}$$
(10)
$$\begin{aligned}&\quad v_i>Q_v \text{ or } w_i>Q_w \nonumber \\ \sum _{i,j \in S} x^w_{ij}&\le |S|-r(S) +\sum _{i\in S} z_{i}&\quad \forall S \subset V' \end{aligned}$$
(11)
$$\begin{aligned} x^w_{ij} + z_{i} + z_{j}&\le 2&\quad \forall i,j \in V' \end{aligned}$$
(12)
$$\begin{aligned} \sum _{k \in S} x^w_{ik} + \sum _{k,l \in S} x^w_{kl} + \sum _{k \in S} x^w_{kj}&\le |S|+2 - z_i-z_j&\quad \forall S \subset V',\ S\ne \emptyset \end{aligned}$$
(13)
$$\begin{aligned} \qquad&i\in V' \setminus S, j\in V'\setminus S, i \ne j \nonumber \\ x^d_{ij}, x^w_{ij}&\in \{0,1\}&\quad \forall (i,j) \in A \end{aligned}$$
(14)
$$\begin{aligned} z_{i}&\in \{0,1\}&\quad \forall i \in V'. \end{aligned}$$
(15)

The objective function (1) represents the total travel time parameterised with respect to \(\alpha \). For example, setting \(\alpha = 0.5\) minimises the total travel time by giving equal importance to both driving and walking, which would effectively assume that the unit cost of driving is the same as the unit cost of walking. Constraints (2) impose degree restrictions on the nodes, with the added constraint that the parking nodes must be visited in both driving and walking routes. As a consequence of these constraints, we limit to one the number of walking routes based on the same parking node. Constraints (3) and (4) are relevant to the driving route ensuring that the vehicle leaves from and returns to the depot, where subtours are prevented by constraints (5). If a node is designated as a parking node, constraints (6) ensure that it is visited by the driver, for which route continuity is ensured by (7). Constraints (8) and (9) play a similar role for the walking routes, where customers with excess volume or weight are excluded by constraints (10). For a given subset S of customers, constraints (11) ensure that there are a sufficient number of parking nodes if they are served by a walking route. Finally, (12) and (13) are path elimination constraints that eliminate walking paths that do not start from and end at the same parking node.

Assuming that unit costs of walking and driving per unit time are available, the model above can be used to minimise the total costs (as opposed to time). In this case, the model should be run with different values of \(\alpha \), and after evaluating the cost of each resulting solution, the one with the minimum total cost should be chosen. This would entail solving the above formulation for as many times as the number of different values of \(\alpha \) used. We provide numerical results on the way that the walking and driving times change with different values of \(\alpha \) in Sect. 4.

3.2 Valid inequalities

This section presents two new sets of valid inequalities to enhance the algorithm that will be described later in this section. Both inequalities serve to improve the connectivity of two graphs, one induced by the driving and the other by the walking routes. To shed some light into the connectivity just mentioned, we present a fractional solution in Fig. 2 for an instance with seven nodes with two disconnected components, where driving (solid lines) and walking (dashed lines) takes place in each component. The solution is obtained at the root node of the branch-and-cut tree whilst solving the formulation above, where the integrality constraints (14) are relaxed, but where a number of subtour (5), walking capacity (11) and path elimination (12) and (13) constraints are added. The support graph in Fig. 2 shows that the component defined on nodes 1, 3, 4 and 5 is disconnected from the other component.

Fig. 2
figure 2

A fractional solution where the numbers on the arcs correspond to the value of the \(x_{ij}^d\) variables (solid lines) and \(x_{ij}^w\) variables (dashed lines) in the fractional solution. In this solution, the subset \(S^d=\{4,5\}\) violates an extended subtour elimination constraint for the driving route and the subset \(S^w=\{1,3,4,5\}\) violates an extended subtour elimination constraint for the walking routes. Only node 2 is used as a rendezvous node, i.e., \(z_2=1\)

Full size image

We now describe the new inequalities that are used to connect disconnected solutions such as the one shown in Fig. 2. Let \(S \subset V'\) be a subset of customers with \(|S| \ge 2\) and let \(\gamma (S)\) denote the set of arcs with both endpoints in set S.

Proposition 1

The following extended subtour inequalities for the driving routes are valid for the problem:

$$\begin{aligned} \sum _{i\in S}\sum _{j\notin S}x^d_{ij} \ge \frac{x^d(\gamma (S))}{|S|-1}. \end{aligned}$$
(16)

Proof

The inequality is trivially valid if \(x^d(\gamma (S)) = 0\). If not, then \(x^d(\gamma (S)) = \sum _{(i,j) \in S} x^d_{ij} \ge 1\), which means that the driver visits any customer in S from another customer in S, and which implies that there must be at least one arc from S to \(V\setminus S\) used by the driver due to constraint (5). Furthermore, \(x^d(\gamma (S))\le |S|-1\) is always satisfied in any feasible solution where an arc with both endpoints in S is used by the driver, so the right hand side of (16) is always less than or equal to 1. \(\square \)

Similar constraints on the arcs arriving to subset S can also be written as follows.

$$\begin{aligned} \sum _{i\notin S}\sum _{j\in S}x^d_{ij} \ge \frac{x^d(\gamma (S))}{|S|-1}. \end{aligned}$$
(17)

The extended subtour inequalities can also be written for the walking routes. In this case the z variables related to the parking nodes are included in the inequality since the number of arcs used in a walking route either leaving or arriving into subset \(S\in V'\) will also depend on the number of nodes in S used as parking nodes.

Proposition 2

The following extended subtour inequalities for the walking routes are valid for the problem:

$$\begin{aligned}&\sum _{i\in S}\sum _{j\notin S}x^w_{ij} + \sum _{i \in S} z_i \ge \frac{x^w(\gamma (S))}{|S|-r(S)+1}, \end{aligned}$$
(18)
$$\begin{aligned}&\sum _{i\notin S}\sum _{j\in S}x^w_{ij} + \sum _{j \in S} z_j \ge \frac{x^w(\gamma (S))}{|S|-r(S)+1}. \end{aligned}$$
(19)

Proof

We differentiate between the following cases:

  • If all customers in S are served by the driver, i.e., no walking is required, then the right hand side of either inequality is equal to 0 and they are trivially satisfied.

  • If there is a feasible closed walking path in S then \(x^w(\gamma (S))=|S|\) and \(r(S)=1\), so the right hand side of either inequality would take the value 1, which would in turn require the left hand side to attain the value of at least 1. This implies that at least one node of S is a parking node, or at least one arc either leaves S through (18) or arrives into S through (19), or both.

  • If \(r(S)>1\), this means that at least two walking routes are required to serve the demand in S. In this case, the maximum value \(\frac{|S|}{|S|+1-|S|/2}\) of the right hand side is given when \(r(S)=|S|/2\). Then, the right hand side of either inequality could increase by at most r(S), which would necessitate the use of an additional parking node, or an additional arc either leaving S through (18) or arriving into S through (19). In this case, the right hand side of either inequality is less than |S|/2 and the left hand side of either inequality should be greater than or equal to |S|/2 since that would be the minimum number of parking nodes needed.

\(\square \)

It is worth highlighting that the denominator of the right hand side of constraints (18) and (19) is one unit greater than the denominator used in right hand side of (16) and (17) due to the fact that a closed route could be performed by walking in S, but not by driving. The driving must always start from and end at the depot.

Revisiting the example in Fig. 2, it can now be seen that the inequalities (16) and (17) are violated for \(S=\{4,5\}\), and the inequalities (18) and (19) are violated for \(S=\{1,3,4,5\}\).

3.3 Branch-and-cut algorithm

The branch-and-cut algorithm is designed to optimally solve the integer linear programming formulation presented in Sect. 3.1 and starts by solving the linear programming (LP) relaxation defined by the objective function (1) and constraints (2), (3), (4), (6), (7), (8), (9), (10). Constraints (5), (11), (13) and the integrality constraints (14) are relaxed at the root node; instead the \(x^w\) and \(x^d\) variables are restricted to be within the interval [0, 1].

Each node of the branching tree consists of a linear programme (LP) with a given set of constraints solved to optimality, and a set of separation routines which aim to find violated constraints. If a violated constraint is identified at a given node, it is added to the LP of that node and the LP is reoptimised. This process is repeated until the separation routines are not able to find any new violated inequality, following which the branch-and-cut algorithm will branch on a variable which violates the integrality constraints (14).

For a fractional solution \({\bar{x}}\) of the LP, we construct a support graph \(\overline{G^w} = (V,\overline{A^w})\) for the walking routes where \((i,j) \in \overline{A^w}\) iff \({\bar{x}}^w_{ij}>0\). Similarly, let \(\overline{G^d} = (V,\overline{A^d})\) be the support graph for the driver where \((i,j) \in \overline{A^d}\) iff \({\bar{x}}^d_{ij}>0\). At any node of the branching tree, one of the two cases below may arise, depending on whether the integrality conditions (14) are satisfied:

  1. 1.

    If the solution of an LP at a given node is integer, then we check if constraints (5), (11) and (13) are also satisfied. This can be done in linear time by shrinking the nodes in the support graphs \(\overline{G^d}\) and \(\overline{G^w}\). If no violated inequalities are identified, then the solution at that node would be feasible and no further branching would be needed. Otherwise, the LP is reoptimised after appending the violated inequality. We note that the addition of violated inequalities to the LP may once again make the solution fractional.

  2. 2.

    If the solution to the LP is fractional we use separation routines for the following inequalities:

    1. (a)

      Subtour/capacity inequalities for the driving route (5).

    2. (b)

      Capacity inequalities for the walking routes (11).

    3. (c)

      Extended subtour inequalities for the driving and walking routes (18), (19), (16) and (17) as described in Sect. 3.2.

    4. (d)

      Path elimination constraints (13).

The separation procedures used for the extended subtour elimination constraints are presented in Algorithms 1 and 2 for the driving and walking routes, respectively. In both cases we first compute the connected components of the support graph (line 2 in Algorithms 1 and 2) and check all the nodes of each connected component (lines 4–5 in Algorithms 1 and 2). Then we apply a greedy heuristic to build promising subsets S in order to identify new cuts (5), (16) and (17) for the driving route (lines 8–20 in Algorithms 1) and (11), (18) and (19) for the walking routes (lines 13–25 in Algorithm 2). Path elimination constraints for the walking routes (13) are checked only when the corresponding connected components are integer.

figure a
figure b

The separation routines shown in Algorithms 1 and 2 differ from the classical separation procedures in that one separation routine is used to find inequalities from different families. The main reason behind this idea is that we use several inequalities that use a subset S, which would need to satisfy similar conditions across the different inequalities. A detailed computational analysis of how the cutting plane approach used at each node of the branching tree enhances the branch-and-cut algorithm is presented in Sect. 4.

4 Computational experiments

The branch-and-cut algorithm was tested on a set of instances, extracted from a real data set of a courier operating in the Southwark area in London. The original data set contained information on 117 deliveries, including location of delivery, and volume and weight of each item delivered. The weights of the items ranged between 2 and 5 kg (with an average of 1.41 kg) and the volumes ranged between 20 and 82 L (with an average of 21.8 L), with \(Q_v\) = 200 L and \(Q_w\) = 25 kg. From this data set, we create 26 instances by randomly sampling between \(|V'| = 5, 6, \ldots , 30\) unique locations, include the depot, and repeat this five times for each value of \(|V'|\) to result in a total of \(26 \times 5 = 130\) instances. Walking and driving times between all pairs of points were computed off-line using the Google Maps Distance Matrix API. The experiments were run using a computer with an Intel(R) Core (TM) i9-9980XE CPU at 3.00 GHz running Windows 10 with 64 bits and with 64 GB of RAM. The algorithm was coded in MVS2017 in which CPLEX 12.9 was used as the optimiser, and part of which included the use of the lazy and cut callback functions for Algorithms 1 and 2 . For comparison purposes, each instance was subject to a computational time limit of 3600 s. Testing was performed using the four following Branch-and-Cut (BC) solution strategies:

  • BC1 is a branch-and-bound algorithm that uses no cutting planes, and where the path elimination constraints (13) are separated only on integer nodes.

  • BC2 is the same as BC1 in which an initial solution (and a corresponding upper bound) is provided to the solver. The initial solution assumes that all deliveries are made by driving, obtained by the solution of the corresponding Travelling Salesman Problem, which can be solved efficiently using Concorde (http://www.math.uwaterloo.ca/tsp/concorde.html).

  • BC3 is a branch-and-cut algorithm, initialised with an upper bound as in BC2, and where Algorithms 1 and 2 are used at each node of the branching tree only for separating inequalities (5) and (11).

  • BC4 is similar to BC3 which additionally uses the extended subtour inequalities (16), (17), (18) and (19).

Table 2 shows the computational results obtained with the four different strategies, all using the setting \(\alpha =0.9\), where each line is an average of five instances. The column ‘Best UB’ shows, the smallest average objective value obtained by any of the four strategies. The table also presents, for each strategy, the average computational time (in seconds) required to solve the corresponding instance under column ‘Time’, and the average final optimality gap (in percent) under column ‘Gap’. An optimality gap less than \(0.01\%\) for a given strategy indicates that all five instances with \(|V'|\) customers were solved to optimality within the time limit by that strategy. The smallest average solution times in each line are indicated in bold. Further details on each strategy, such as the number of cuts added, the number of nodes in the branch and bound tree, and the lower bounds obtained at the root node are presented in the appendix, for each individual instance.

Table 2 Comparison results between the four different BC strategies showing average values of five instances per line with \(\alpha = 0.9\)
Full size table
Table 3 Computational results with BC4 and \(\alpha \in \{0.6, 0.8, 0.9\}\)
Full size table

The results shown in Table 2 suggest that, for relatively small-sized instances (i.e., those with \(|V'| \le 14\)), BC1 and BC2 provide optimal solutions quicker than the other variants, but BC3, and more so BC4 present a superior performance as the instance size increases, either in terms of the computational time or the final optimality gap. Whilst a great majority of instances have been optimally solved, the difficulty of solving others of up to 30 nodes suggests that the complexity of this problem is well beyond that of the VRP. Overall, BC4 yields the smallest average optimality gap and the fastest solution time) for all instances, for which reason it is used in the remainder of the experiments.

Next, we investigate the impact of changing the multiplier \(\alpha \) on the total driving and walking times in the resulting solutions. To this end, we test the values in the set \(\alpha \in \{0.6, 0.8, 0.9\}\). The results are given in Table 3, where each row shows the average results for the group of five instances with \(|V'|\) customers. The column titled ‘Gap/Time’ shows either the average solution time (in seconds), or the average optimality gap (appended by ‘%’) if optimality was not attained within the one hour time limit. For gaps that are strictly positive, the best solutions obtained were used to analyse the results. The column ‘D/W’ shows the percentage split between the total time spent driving (D) and that of walking (W) in the resulting solutions. To provide a meaningful figures, the stem driving time, defined as the time taken to drive (and back) from the depot to the closest customer and equal to 1897 s in our data set, has been excluded in the calculation of the driving time in all instances.

As Table 3 shows, all instances using the setting \(\alpha = 0.6\) have been solved to optimality. In this case, the majority of instances yielded optimal solutions that only make use of driving, which corresponds to the initial solutions used in the BC algorithm for upper bounds. The amount of driving is reduced using the settings \(\alpha =0.8\) and \(\alpha = 0.9\), but this increases the difficulty of solving the model to optimality. In the case of \(\alpha = 0.9\), for example, the algorithm has failed to find optimal solutions for some instances with 21, 27, 28 and 30 customers, although the optimality gap for these instances remains within 2%, and generally below 1%. The results also indicate that, using \(\alpha =0.6\), an average of \(99.3\%\) of the deliveries are done by driving, which is reduced to \(51\%\) when \(\alpha =0.8\) and further to \(2.7\%\) when \(\alpha =0.9\), showing the significant effect that the parameter \(\alpha \) has on the resulting solutions.

For illustrative purposes, we provide a sample visualisation of solutions for an instance with 28 addresses in Fig. 3 using different values for \(\alpha \), namely 0.6, 0.8 and 0.9. It should be kept in mind that the routes shown in the figure are designed with respect to the capacity and volume restrictions of the items and based on street-level distance matrices.

Fig. 3
figure 3

Visualisation of the best solutions found for an instance with 28 customers, using the settings \(\alpha = 0.6\) (left), \(\alpha = 0.8\) (middle) and \(\alpha = 0.9\) (right). The shaded dots represent the customer locations and the filled dots indicate the parking nodes. Red lines correspond to driving and black lines show the walking routes. The depot, located further east, is not depicted in the figures. (Leaflet, ©OpenStreetMap contributors, ©CartoDB)

Full size image

5 Conclusions

This paper proposes that a combined driving and walking delivery system could be a viable alternative to last-mile delivery problems in urban areas to mitigate the adverse effects of more traditional means of distribution, such as traffic congestion and use of kerbside space. Our primary contribution has been to describe such a system, to propose a mathematical model along with valid inequalities, and an optimisation algorithm to prescribe solutions for such a system to operate in practice. Computational experience suggests that, for instances of up to 30 customers, the algorithm is able to either provide optimal solutions, or those that are within a reasonable gap to the optimal value. The results also suggest that the complexity of the problem increases not only with an increased number of deliveries but also when the relative weighting of one mode over the other changes. In particular, the problem becomes more difficult to solve with increased amount of walking, which may be explained by the additional constraints, such as weight and volume restrictions, required to model the walking paths.

5.1 Practical considerations

One practical consideration that may arise in the implementation of the proposed model relates to the application of the model to larger areas. In reality, carriers divide up geographical areas into sub-rounds which consist of postcode patches that are linked together by the road network. Once this segmentation is undertaken, each patch and their postcode constituent parts would be identified to create rounds using the proposed system for a single driver operating on that patch. The specific elements of the round would lend to walking routes where drop-off points for the van would be identified, and from where the walking routes would start and finish.

One other practical aspect concerns the walking routes. Practical work undertaken by the authors (Allen et al. 2018) suggested that many drivers preferred walking as opposed to continually moving the vehicle very small distances while trying to find parking spaces. This crucially depends on the average size of parcel being delivered as the number of trips needed back and forth to the van will increase as package size increases. Our previous surveys undertaken in London showed drivers walking up to 10 km on a round. As deliveries were all within a relatively small area, it would be feasible to walk the entire round. For this reason, the amount of walking undertaken was not constrained for the application context considered here, particularly as the objective of the problem, in part, serves to minimise the length of the walking tours. Also, for larger areas (and less dense deliveries), driving would naturally take precedence as being a faster mode of delivery over longer distances. However, it is possible to introduce constraints into the model, if the walking distances were considered to be too great.

5.2 Future research directions

The nature of day-to-day parcel delivery operations within urban areas present further challenges that may form the basis of future work. First, and further to the discussion in Sect. 2, it is possible to introduce time window constraints into the model described in this paper to account for time-sensitive parcels. Such constraints, however, may result in non-linearities for particular types of formulations (Nguyên et al. 2019), and whilst they can be linearised, the additional constraints needed can significantly increase the complexity of solving the model to optimality. Alternative formulations or solution methods that can effectively deal with time window constraints is one possible avenue for further research. Second, traffic conditions within urban areas present a degree of uncertainty which impacts on driving times, whereas there is very little uncertainty associated with walking times. Incorporating such partial stochasticity falls within the class of stochastic vehicle routing problems (Gendreau et al. 2016), which are considerably more difficult to solve compared to their deterministic counterparts. Third, whilst a great majority of the deliveries are known and planned for ahead of the actual operations, there are occasional requests (e.g. collections) that dynamically arise. Accounting for such requests would result in a dynamic routing problem (Psaraftis et al. 2016), which is also significantly harder to solve than the static variants. Fourth, whereas capacity limitations have been explicitly represented in our model, this has assumed a two-dimensional treatment of the parcels. Practical requirements may dictate that both the weight and the volume of the parcels would need to be explicitly represented, and that capacity is modelled as a three-dimensional constraint. Finally, the problem as described in this paper can be extended to allow multiple vehicles for applications in larger areas.