- Email: [email protected]
Improved lower bounds for the Split Delivery Vehicle Routing Problem
Operations Research Letters 38 (2010) 302–306
Contents lists available at ScienceDirect
Operations Research Letters journal homepage: www.elsevier.com/locate/orl
Improved lower bounds for the Split Delivery Vehicle Routing Problem Lorenza Moreno a,∗ , Marcus Poggi de Aragão a , Eduardo Uchoa b a
Departamento de Informática, PUC-Rio, Brazil
b
Departamento de Engenharia de Produção, Universidade Federal Fluminense, Brazil
article
info
Article history: Received 31 July 2008 Accepted 22 February 2010 Available online 20 April 2010 Keywords: Split Delivery Extended formulation Column generation Cutting-plane algorithm
abstract This paper presents an algorithm to obtain lower bounds for the Split Delivery Vehicle Routing Problem. An extended formulation over a large set of variables is provided and valid inequalities are identified. The algorithm combined column and cut generation and improved the best known lower bounds for all instances from the literature. Some reasonably sized instances are solved to optimality for the first time. © 2010 Elsevier B.V. All rights reserved.
1. Introduction The Split Delivery Vehicle Routing Problem (SDVRP) is a relaxation of the classical Capacitated Vehicle Routing Problem (CVRP) where a client can be attended by more than one vehicle. The problem was introduced to the literature by Dror and Trudeau [6], who empirically showed that allowing split deliveries can significantly reduce the total distance travelled. Archetti et al. [1] prove that this reduction can achieve up to 50%. Belenguer et al. [2] provide good lower bounds for SDVRP using a cutting-plane algorithm. Jin et al. [9] improve the lower bounds for several instances with a column generation algorithm. For the same set of instances, the average gap is respectively 10.5% and 6.5%. The few exact algorithms for the SDVRP are able to solve only small instances. Dror et al. [5] present an integer programming formulation and valid inequalities for the SDVRP. They developed a cutting-plane algorithm and solved to optimality an instance with 10 clients. Lee et al. [10] use a dynamic programming model to solve instances with 9 clients and 6 vehicles. A two-stage algorithm with clustering and routing phases is proposed by Jin et al. [8]. It can successfully solve instances with up to 23 clients. Bigger instances, even those that are easily solved when splits are not allowed (by CVRP algorithms), could not be solved by any known SDVRP algorithm. The main difficulty lies in the fact that in the SDVRP the amount of goods delivered by a vehicle in each visit to a client
∗ Corresponding address: R. Marquês de São Vicente, 225 RDC, Rio de Janeiro RJ, 22453-900, Brazil. E-mail address: [email protected] (L. Moreno). 0167-6377/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.orl.2010.04.008
becomes a decision variable, leading to MIP formulations with significantly larger integrality gaps. This work uses column generation and cutting planes to obtain stronger lower bounds to the SDVRP. In Section 2, we present formulations and valid inequalities for the problem. The proposed algorithm is discussed in Section 3. Experimental results, reported in Section 4, show that the lower bounds of all instances tested were improved. The algorithm yielded optimal integer solutions in 11 cases. 2. Formulation and valid inequalities Let G = (V , A) be a directed graph. Let V = {0, . . . , n} be the vertex set, where 0 represents a depot and the remaining vertices are client locations. For each client i, let di denote its demand, a positive quantity of goods to be delivered at location i. The arc set of G is defined by A = {(i, j) : i, j ∈ V , i 6= j} and each arc a = (i, j) has an associated distance ca . There are m identical vehicles available, their capacity is denoted by Q . The SDVRP consists of finding m vehicle routes and corresponding deliveries, such that: (1) each vehicle route starts and ends at the depot; (2) vehicles are only loaded at the depot, respecting the capacities Q ; and (3) each client i receives deliveries, from one or more vehicles, totaling di units. The goal is to minimize the total distance travelled. We present an extended formulation for the SDVRP. Some extra notation is necessary. Assume that Q and all demands di have integer values and that d0 = 0. Let Di (q) = min{q, di } denote the maximum amount that can be delivered to client i when the total load available on the vehicle is q. Let GQ = (V , AQ ) be a directed multigraph. For each arc a = (i, j) ∈ A \ δ − (0), there are arcs aqd = (i, j)qd ∈ AQ , one for each possible load size q = {1, . . . , Q } and each possible delivery quantity d = {1, . . . , Dj (q)}. For each
L. Moreno et al. / Operations Research Letters 38 (2010) 302–306
arc arriving at the depot a = (i, 0) ∈ δ − (0), there is one arc a00 ∈ AQ representing an empty vehicle (q = d = 0). It is assumed that δQ− (S ) and δQ+ (S ) are the subsets of arcs in AQ with any indices q
qd
and the elimination of the xa variables, by replacing those variables in (10) and (11) using (9).
qd
and d entering and leaving S respectively. Let xa be an integer variable indicating the number of vehicles that traverse arc a = (i, j) with exactly q units of goods after delivering d units to client i. Let V+ = {1, . . . , n}. Minimize
X
ca .xqd a
(1)
X
Minimize
r ∈R
aqd ∈A
s.t.
d.
= d i ∀ i ∈ V+
s.t.
X
a∈δ − (i)
Dj (q)
x(aq+d)d =
X
X
xqd a
a=(i,j)∈δ + (i) d=1
d=1
∀ i ∈ V+ , q = 1 , . . . , Q − 1 i (Q ) X DX
a∈δ − (i) d=1
X
x(aq+d)d = x(i,0) q0
∀i ∈ V+ , q = 0
x00 (i,0) = m
(3) (4)
∀aqd ∈ AQ .
(5) (6)
Eqs. (3), (4), (5) and the redundant equation
X
qd
x(0,i) = m
(7)
(0,i)qd ∈δQ+ (0)
λr = m
(16)
λr ≥ 0 ∀r ∈ R.
(17)
The DWM has an exponential number of variables and is solved by column generation. The pricing subproblem aims at finding minimum cost paths. Let (18) represent a generic cut (or original qd constraint) i over the extended variables xa ,
X
αaiqd .xqd a ≥ bi ,
(18)
qd
where αai is the coefficient of arc aqd in the cut and bi is its righthand side coefficient. Let t be the current number of constraints and cuts in the DWM and βi be the dual variable associated to the ith constraint. The contribution of each extended arc aqd to a path qd reduced cost (denoted by c¯a ) is: c¯aqd = ca −
t X
αaiqd βi .
(19)
X
Hence, the reduced cost c¯ of a path is:
X
c¯ =
X
2.1. Valid inequalities 2.1.1. Homogeneous extended capacity cuts For every set S ⊆ V+ the following balance equality is valid:
X
(8)
qd qd pqd ∈ AQ ar .λr − xa = 0 ∀a
− aqd ∈δQ (S )
yq =
(21)
X
X
xqd a
= d i ∀ i ∈ V+
qd
q = 1, . . . , Q ,
(22)
qd
q = 1, . . . , Q .
(23)
x(i,j) ,
Dj (q)
zq =
X
X
(i,j)∈δ + (S ) d=1
(10)
x(i,j) ,
The balance equation over S with aggregate variables becomes:
aqd ∈δ − (i)
X
q.xqd a = d(S ).
+ aqd ∈δQ (S )
Let yq and zq be aggregated variables defined as follows:
(9)
r ∈R
x00 (i,0) = m
(11)
(i,0)∈δ − (0)
λr ≥ 0 ∀r ∈ R xqd a
X
q.xqd a −
(i,j)∈δ − (S ) d=1
d.
(20)
Dj (q)
ca .xqd a
X
pqd caqd . ar .¯
aqd ∈AQ
aqd ∈AQ
s.t.
(15)
i=1
can be viewed as defining a network flow of m units over an acyclic layered graph GN = (VN , AN ). For each vertex i ∈ V+ there are Q vertices in VN , one for each possible load size q = 0, . . . , Q − 1. The depot corresponds to two additional vertices in VN . Each arc aqd ∈ AQ corresponds to an arc in AN . As the flow of m units in GN has just one source and one destination (the vertices associated to the depot) any integer solution can be decomposed into a set of m paths, each one corresponding to a route and a delivery pattern. Constraints (2) impose that each client receives exactly its demand. The direct use of this formulation is not practical even for small sized instances, so we provide an equivalent reformulation suitable to be solved by column generation. Let R denote the set of origin– destination paths in GN , each one corresponding to a possible route and delivery pattern. For each r ∈ R, let λr be a continuous variqd able and let par indicate whether a vehicle executing r traverses arc a = (i, j) with q units of goods and delivers exactly d units to j. Minimize
∀ i ∈ V+
aqd ∈δ − (i)
aqd ∈AQ
(i,0)∈δ − (0)
xqd a ∈ Z+
.λr = di d.pqd ar
r ∈R
aqd ∈δQ (i) D i ( Q −q )
X
r ∈R
(2)
(14)
Q
−
X
.λr ca .pqd ar
X X
xqd a
X
aqd ∈AQ
X
303
∈ Z ∀a ∈ AQ . qd
(12) (13)
It can be seen that Equalities (3) and (4) are implied by the definition of R and by (9). The so-called Dantzig–Wolfe Master problem (DWM) is obtained with the relaxation of integrality constraints
Q X q =1
q.yq −
Q −1 X
q.zq = d(S ).
(24)
q =1
For each possible values of Q and D = d(S ), define P (Q , D) as the polyhedron induced by the integer solutions of (24). The Homogeneous Extended Capacity Cuts (HECCS) [13] are all inequalities valid for those polyhedra. The HECCs, already used on other problems [11,12], are also valid for the SDVRP.
304
L. Moreno et al. / Operations Research Letters 38 (2010) 302–306
2.1.2. Edge upper bound
3.1. Column generation
Dror and Trudeau [6] stated the following theorem: Theorem 1. If distance matrix satisfies triangular inequality, there is an optimal solution to the SDVRP where any two paths have at most one client in common. In other words, this theorem assures that there is an optimal solution where two clients appear on the same path at most once. A weaker property also valid for SDVRP is: Corollary 1. If distance matrix satisfies triangular inequality, there is an optimal solution where, for any two vertices i, j ∈ V + , i 6= j, the edge (i, j) of the undirected graph is traversed at most once. Thus, Eq. (25) is a valid inequality for SDVRP. Dj (q) Q X X q =1
d=1
Di (q) qd x(i,j)
+
X
! qd x(j,i)
≤ 1,
∀i, j ∈ V+ , i 6= j.
(25)
d=1
The requirement of a distance matrix satisfying the triangle inequality in the above propositions is actually not restrictive. A preprocessing procedure can always be used to obtain such a distance matrix. If the direct distance between clients i and j, c(i,j) , is greater than the distance from i to another client k and then to j, the value of c(i,j) can be replaced by c(i,k) + c(k,j) . This transformation is valid because the arc (i, j) can now be interpreted as a vehicle travelling from i to k (and delivering nothing) and then travelling from k to j. We applied this transformation to all instances in our tests, so inequality (25) can always be used.
2.1.3. Split deliveries cuts Eqs. (2) impose that the sum of all deliveries at vertex i must be equal to the amount requested by i. Valid inequalities of the polyhedron induced by the integer solutions of (2) are referred to as Split Deliveries Cuts (SDCs). These inequalities are valid for SDVRP.
The generation of minimum reduced cost paths is a NP-hard problem that can be solved in pseudo-polynomial time. Note that the SDVRP definition allows paths that visit the same client more than once, i.e., cycling is permitted. The generation of paths can be executed by a dynamic programming algorithm (DP) with O(|V |.Q ) states and pseudo-polynomial complexity O(|A|.Q 2 ). When triangular inequality holds, there is always an optimal solution in which no path has a cycle. Thus, we eliminate paths with 2-cycles (subpaths i → j → i, i, j ∈ V +) in order to strengthen the formulation. As already noted by [4], this elimination of 2-cycles does not change the pricing complexity. In order to accelerate the pricing, two dynamic programming based heuristics (suggested by Fukasawa et al. [7]) are also used. The scaling heuristic consists of performing the dynamic programming over a simplified instance obtained by dividing the capacity by an integer g > 1 (and rounding down) and also dividing all demands by the same g (but rounding up). A path with negative reduced cost obtained in this simplified instance corresponds to a path with negative reduced cost in the original instance. The sparsification heuristic restricts the set of arcs used in the dynamic programming, for each client i only the k cheaper arcs incident to it are considered. Of course, the dynamic programming based heuristics can fail to find a minimum reduced cost path. Hence, column generation stops only when the pricing without any acceleration technique is used and no path with negative reduced cost is found. 3.2. Cutting At each cut generation iteration, the current fractional solution qd of formulation (14)–(17) is converted into xa variables using Eqs. (9). Violated cuts are then separated and introduced into the DWM LP in format (18). Inequalities (25) are separated by inspection, while Homogeneous Extended Capacity Cuts are separated by the same procedure described in [11]. The separation of Split Delivery Cuts begins by discarding all clients that do not have split deliveries in the current solution. Then, we look for violated inequalities of the form:
X
drdeaqd = drdi e,
(26)
−
aqd ∈δQ (S )
3. Lower bounding algorithm The lower bounding algorithm combines column and cut generation over the DWM relaxation presented in Section 2. That linear program is initialized with a set of variables (paths) constructed as follows:
P −1 • The first path visits clients 1, 2, . . . , j such that ij= 1 di < Pj Q and i=1 di ≥ Q . It delivers all the demand of clients Pj−1 1, 2, . . . , (j − 1) and Q − i=1 di units of client j. The next path transports the remaining units of j and all demands of the following clients until the vehicle capacity is achieved again, and so on.
• There is also |V+ | additional paths, each one delivering all demand of a single client. At the beginning, the DWM only contains the constraints from the original formulation, (15)–(17). Since the linear relaxation of a set covering problem is easier to solve by column generation than a set partitioning problem, Eqs. (15) are relaxed to ≥. This change has no impact on the lower bounds obtained.
where 0 < r < 1. The example below shows a fractional solution and a violated inequality. Suppose that the demand of a client i is 7 and that its demand is covered by 3 distinct paths in a fractional solution. Each path is associated to a variable λr with value 0.5. The first and second paths are cycle-free and deliver 6 and 4 units to client i respectively. The third path visits the client twice and delivers 2 units at each stop. Relaxing equation (26) to ≥ and using r = 0.5, a violated inequality is found: 4w7 + 3w6 + 3w5 + 2w4 + 2w3 + w2 + w1 ≥ 4.
(27)
4. Experiments The tests were executed in a machine with processor Intel Core 2 Duo 2 GHz and with 3 GB of RAM. CPLEX 11.0 was used to solve the linear programs. This algorithm was tested for 3 distinct sets of instances found in the SDVRP Pnliterature. In all those sets the number of vehicles m is set as d i=1 di /Q e. The first set was generated by Belenguer et al. [2] and has 14 instances with 50, 75 or 100 clients and vehicle capacity of 160 units. Instances with name ending in ‘‘d1’’ have demands ranging between 1% and 10% of vehicle capacity. Instances with names ending
L. Moreno et al. / Operations Research Letters 38 (2010) 302–306 Table 1 Comparison with [2] on their set of instances.
305
Table 3 Results on the set of instances by [3].
Instance
UB
LB [2]
Our LB
Our time
Instance
n
UB [3]
Our LB
Our time
s51d1 s51d2 s51d3 s51d4 s51d5 s51d6 s76d1 s76d2 s76d3 s76d4 s101d1 s101d2 s101d3 s101d5
458 726 972 1677 1440 2327 594 1147 1474 2257 716 1393 1975 2915
454.0 676.6 905.2 1520.7 1272.9 2113.0 584.9 1020.3 1346.3 2011.6 700.6 1271.0 1739.7 2630.4
454.4 694.2 930.7 1539.0 1313.4 2141.7 586.6 1061.1 1395.9 2046.1 704.9 1337.1 1832.2 2737.1
88 123 200 609 369 525 125 393 557 1438 753 1025 1683 1168
8.6%
5.9%
sd1 sd2 sd3 sd4 sd5 sd6 sd7 sd8 sd9 sd10 sd11 sd12 sd13 sd14 sd15 sd16 sd17 sd18 sd19 sd20 sd21
8 16 16 24 32 32 40 48 48 64 80 80 96 120 144 144 160 160 192 240 288
228.3 714.4 430.6 631.1 1 408.1 831.2 3 714.4 5 200.0 2 059.8 2 749.1 13 612.1 7 399.1 10 367.1 11 023.0 15 271.8 3 449.1 26 665.8 14 546.6 20 559.2 40 408.2 11 491.7
228.3 708.3 430.6 631.0 1 390.6 831.2 3 640.0 5 068.3 2 044.2 2 684.9 13 275.0 7 175.8 10 053.6 10 588.2 14 908.5 3 381.3 26 317.2 14 029.2 19 707.2 39 252.8 11 271.0
0 1 0 2 6 0 56 92 11 70 431 748 2045 165 2703 6 2192 3060 723 3785 205
Avg. gap
Table 2 Comparison with [9] on the set of instances by [2]. Instance s51d1 s51d2 s51d3 s51d4 s51d5 s51d6 s76d1 s76d2 s76d3 s76d4 s101d1 s101d2 s101d3 s101d5 Avg. gap
UB
LB [9]
Our LB
723.4 968.9 1586.5 1355.5 2197.8
695.0 922.7 1505.4 1297.5 2108.6
1185.7 1504.9 2136.4
1066.2 1397.4 2019.9
1474.5 2012.9 2846.2
1349.8 1837.3 2725.5
457.1 700.4 938.5 1549.7 1318.9 2154.7 592.6 1071.3 1407.5 2059.8 716.8 1358.9 1853.1 2767.6
6.5%
5.0%
Time [9] (s) 3555 323 14 2 <1 7511 549 188 16 507 2698 10
Our time 184 131 225 631 320 529 146 525 771 1554 671 1326 1068 1670
in ‘‘d2’’ to ‘‘d6’’ have demands ranging between the following fractions of the vehicle capacity: (10%–30%), (10%–50%), (10%–90%), (30%–70%) and (70%–90%). Tables 1 and 2 compare our results on this first set of instances with those obtained by Belenguer et al. [2] using a cutting-plane algorithm and with those obtained by Jin et al. [9] using column generation. While [9] compute euclidean distances without rounding, as double precision values, [2] round them to the nearest integer using the TSPLIB convention. Therefore, we show separate results considering rounded distances and fractional distances. It can be seem that improved bounds were obtained in all those instances. All times presented in the tables are in seconds. The times of [9] were obtained in a machine with a processor Pentium 4, 2 GHz and 2 GB memory. Although their algorithm can be significantly faster than ours on instances with large demands, it can also be significantly slower on instances with small demands (they could not run instances of type ‘‘d1’’). The newly proposed lower bounding algorithm is much more stable in terms of running time. Article [2] did not provide times for comparison. The second set has 21 instances created by Chen et al. [3] (Table 3). In all problems, vehicle capacity is 100 and client demand is either 60 or 90. These instances are highly symmetric, with clients placed in concentric circles around the depot. The euclidean distances are not rounded. Those authors used a heuristic algorithm to provide upper bounds. The average integrality gap for these instances is 1.6%. Lower bound values in bold are optimal integer solutions, they were obtained for 9 instances of this set. The largest solved instance has 288 clients. The third set of instances is composed by some instances from the TSPLIB (Table 4). The results are compared with the lower bounds given by Belenguer et al. [2]. The first 4 instances of this
Avg. gap
1.6%
Table 4 Comparison with [2] on TSPLIB instances. Instance
n
Q
Q0
UB
LB [2]
eil22 eil23 eil30 eil33 eil51 eila76 eilb76 eilc76 eild76 eila101 eilb101
21 22 29 32 50 75 75 75 75 100 100
6000 4500 4500 8000 160 140 100 180 220 200 112
60 4500 180 800 160 140 100 180 220 200 112
375 569 510 835 521 832 1023 735 683 817 1077
375.0 569.0 508.0 833.0 511.6 782.7 937.5 706.0 659.4 793.5 1005.9
375.0
3
510.0 834.7 517.8 807.6 981.4 717.8 666.1 799.8 1040.6
157 1065 92 101 80 124 167 632 249
3.2%
1.6%
Avg. gap
Our LB
Our time
set have very large vehicle capacities Q . Since the complexity of our pricing subproblem is O(|A|.Q 2 ), we compute the greatest common divisor g between the vehicle capacity and each demand value. Then, each demand value di can be replaced by d0i = di /g and the vehicle capacity Q is reduced to Q 0 = Q /g. We could not run our algorithm in instance eil23 as Q 0 is still 4500. Optimal integer solutions were obtained for instances eil22 and eil30. Although the solution of instance eil33 is fractional, the lower bound of 834.7 is enough to prove the optimality of a previously known solution with value 835. The algorithm proposed is quite different from the previous algorithms to the SDVRP. The extended formulation provides a more detailed view of the problem and a more convenient structure for the identification of effective inequalities. In particular, by using extended variables, the knapsack aspect of the problem could be better explored. The algorithm is able to achieve sharper lower bounds for almost all SDVRP instances tested. However, it has a limitation: it is not practical for problems with very large vehicle capacity values. In these cases, the algorithm becomes quite slow and can even fail to complete due to memory requirements. Acknowledgements LM and MPA were partially financed by CNPq Grants 140715/ 2004-5 and 311997/06-6. EU received support from CNPq Grant 304533/02-5 and Faperj Grant E-26/100.504/2007. The authors also thank the technical support by GAPSO Inc., Brazil.
306
L. Moreno et al. / Operations Research Letters 38 (2010) 302–306
References [1] C. Archetti, M.W.P. Savelsbergh, M.G. Speranza, Worst-case analysis for Split Delivery Vehicle Routing Problems, Transportation Science 40 (2) (2006) 226–234. [2] J.M. Belenguer, M.C. Martinez, E. Mota, A lower bound for the Split Delivery Vehicle Routing Problem, Operations Research 48 (5) (2000) 801–810. [3] S. Chen, B.L. Golden, E.A. Wasil, The Split Delivery Vehicle Routing Problem: applications, algorithms, test problems, and computational results, Networks 49 (4) (2007) 318–329. [4] N. Christofides, A. Mingozzi, P. Toth, Exact algorithms for the vehicle routing problem based on spanning tree and shortest path relaxations, Mathematical Programming 20 (1981) 255–282. [5] M. Dror, G. Laporte, P. Trudeau, Vehicle routing with Split Deliveries, Discrete Applied Mathematics 50 (3) (1994) 239–254. [6] M. Dror, P. Trudeau, Savings by Split Delivery Routing, Transportation Science 23 (1989) 141–145. [7] R. Fukasawa, H. Longo, J. Lysgaard, M. Poggi de Aragão, M. Reis, E. Uchoa, R. Werneck, Robust branch-and-cut-and-price for the capacitated vehicle routing problem, Mathematical Programming 106 (3) (2006) 491–511.
[8] M. Jin, K. Liu, R.O. Bowden, A two-stage algorithm with valid inequalities for the Split Delivery Vehicle Routing Problem, International Journal of Production Economics 105 (1) (2007) 228–242. [9] M. Jin, K. Liu, B. Eksioglu, A column generation approach for the Split Delivery Routing Problem, Operations Research Letters 36 (2) (2008) 265–270. [10] C. Lee, M.A. Epelman, C.C. White III, Y.A. Bozer, A shortest path approach to the multiple-vehicle routing with split pick-ups, Transportation Research Part B: Methodological 40 (4) (2006) 265–284. [11] A. Pessoa, E. Uchoa, M. Poggi de Aragão, A robust branch-cut-and-price algorithm for the heterogeneous fleet vehicle routing problem, Networks 54 (2009) 167–177. [12] A. Pessoa, E. Uchoa, M. Poggi de Aragão, R. Rodrigues, Algorithms over arc-time indexed formulations for single and parallel machine scheduling problems, Technical Report RPEP, vol. 8, No. 8, Universidade Federal Fluminense, Engenharia de Produção, Niterói, Brazil, 2008. [13] E. Uchoa, R. Fukasawa, J. Lysgaard, A. Pessoa, M. Poggi de Aragão, D. Andrade, Robust branch-cut-and-price for the capacitated minimum spanning tree problem over a large extended formulation, Mathematical Programming 112 (2) (2007) 443–472.
Contents lists available at ScienceDirect
Operations Research Letters journal homepage: www.elsevier.com/locate/orl
Improved lower bounds for the Split Delivery Vehicle Routing Problem Lorenza Moreno a,∗ , Marcus Poggi de Aragão a , Eduardo Uchoa b a
Departamento de Informática, PUC-Rio, Brazil
b
Departamento de Engenharia de Produção, Universidade Federal Fluminense, Brazil
article
info
Article history: Received 31 July 2008 Accepted 22 February 2010 Available online 20 April 2010 Keywords: Split Delivery Extended formulation Column generation Cutting-plane algorithm
abstract This paper presents an algorithm to obtain lower bounds for the Split Delivery Vehicle Routing Problem. An extended formulation over a large set of variables is provided and valid inequalities are identified. The algorithm combined column and cut generation and improved the best known lower bounds for all instances from the literature. Some reasonably sized instances are solved to optimality for the first time. © 2010 Elsevier B.V. All rights reserved.
1. Introduction The Split Delivery Vehicle Routing Problem (SDVRP) is a relaxation of the classical Capacitated Vehicle Routing Problem (CVRP) where a client can be attended by more than one vehicle. The problem was introduced to the literature by Dror and Trudeau [6], who empirically showed that allowing split deliveries can significantly reduce the total distance travelled. Archetti et al. [1] prove that this reduction can achieve up to 50%. Belenguer et al. [2] provide good lower bounds for SDVRP using a cutting-plane algorithm. Jin et al. [9] improve the lower bounds for several instances with a column generation algorithm. For the same set of instances, the average gap is respectively 10.5% and 6.5%. The few exact algorithms for the SDVRP are able to solve only small instances. Dror et al. [5] present an integer programming formulation and valid inequalities for the SDVRP. They developed a cutting-plane algorithm and solved to optimality an instance with 10 clients. Lee et al. [10] use a dynamic programming model to solve instances with 9 clients and 6 vehicles. A two-stage algorithm with clustering and routing phases is proposed by Jin et al. [8]. It can successfully solve instances with up to 23 clients. Bigger instances, even those that are easily solved when splits are not allowed (by CVRP algorithms), could not be solved by any known SDVRP algorithm. The main difficulty lies in the fact that in the SDVRP the amount of goods delivered by a vehicle in each visit to a client
∗ Corresponding address: R. Marquês de São Vicente, 225 RDC, Rio de Janeiro RJ, 22453-900, Brazil. E-mail address: [email protected] (L. Moreno). 0167-6377/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.orl.2010.04.008
becomes a decision variable, leading to MIP formulations with significantly larger integrality gaps. This work uses column generation and cutting planes to obtain stronger lower bounds to the SDVRP. In Section 2, we present formulations and valid inequalities for the problem. The proposed algorithm is discussed in Section 3. Experimental results, reported in Section 4, show that the lower bounds of all instances tested were improved. The algorithm yielded optimal integer solutions in 11 cases. 2. Formulation and valid inequalities Let G = (V , A) be a directed graph. Let V = {0, . . . , n} be the vertex set, where 0 represents a depot and the remaining vertices are client locations. For each client i, let di denote its demand, a positive quantity of goods to be delivered at location i. The arc set of G is defined by A = {(i, j) : i, j ∈ V , i 6= j} and each arc a = (i, j) has an associated distance ca . There are m identical vehicles available, their capacity is denoted by Q . The SDVRP consists of finding m vehicle routes and corresponding deliveries, such that: (1) each vehicle route starts and ends at the depot; (2) vehicles are only loaded at the depot, respecting the capacities Q ; and (3) each client i receives deliveries, from one or more vehicles, totaling di units. The goal is to minimize the total distance travelled. We present an extended formulation for the SDVRP. Some extra notation is necessary. Assume that Q and all demands di have integer values and that d0 = 0. Let Di (q) = min{q, di } denote the maximum amount that can be delivered to client i when the total load available on the vehicle is q. Let GQ = (V , AQ ) be a directed multigraph. For each arc a = (i, j) ∈ A \ δ − (0), there are arcs aqd = (i, j)qd ∈ AQ , one for each possible load size q = {1, . . . , Q } and each possible delivery quantity d = {1, . . . , Dj (q)}. For each
L. Moreno et al. / Operations Research Letters 38 (2010) 302–306
arc arriving at the depot a = (i, 0) ∈ δ − (0), there is one arc a00 ∈ AQ representing an empty vehicle (q = d = 0). It is assumed that δQ− (S ) and δQ+ (S ) are the subsets of arcs in AQ with any indices q
qd
and the elimination of the xa variables, by replacing those variables in (10) and (11) using (9).
qd
and d entering and leaving S respectively. Let xa be an integer variable indicating the number of vehicles that traverse arc a = (i, j) with exactly q units of goods after delivering d units to client i. Let V+ = {1, . . . , n}. Minimize
X
ca .xqd a
(1)
X
Minimize
r ∈R
aqd ∈A
s.t.
d.
= d i ∀ i ∈ V+
s.t.
X
a∈δ − (i)
Dj (q)
x(aq+d)d =
X
X
xqd a
a=(i,j)∈δ + (i) d=1
d=1
∀ i ∈ V+ , q = 1 , . . . , Q − 1 i (Q ) X DX
a∈δ − (i) d=1
X
x(aq+d)d = x(i,0) q0
∀i ∈ V+ , q = 0
x00 (i,0) = m
(3) (4)
∀aqd ∈ AQ .
(5) (6)
Eqs. (3), (4), (5) and the redundant equation
X
qd
x(0,i) = m
(7)
(0,i)qd ∈δQ+ (0)
λr = m
(16)
λr ≥ 0 ∀r ∈ R.
(17)
The DWM has an exponential number of variables and is solved by column generation. The pricing subproblem aims at finding minimum cost paths. Let (18) represent a generic cut (or original qd constraint) i over the extended variables xa ,
X
αaiqd .xqd a ≥ bi ,
(18)
qd
where αai is the coefficient of arc aqd in the cut and bi is its righthand side coefficient. Let t be the current number of constraints and cuts in the DWM and βi be the dual variable associated to the ith constraint. The contribution of each extended arc aqd to a path qd reduced cost (denoted by c¯a ) is: c¯aqd = ca −
t X
αaiqd βi .
(19)
X
Hence, the reduced cost c¯ of a path is:
X
c¯ =
X
2.1. Valid inequalities 2.1.1. Homogeneous extended capacity cuts For every set S ⊆ V+ the following balance equality is valid:
X
(8)
qd qd pqd ∈ AQ ar .λr − xa = 0 ∀a
− aqd ∈δQ (S )
yq =
(21)
X
X
xqd a
= d i ∀ i ∈ V+
qd
q = 1, . . . , Q ,
(22)
qd
q = 1, . . . , Q .
(23)
x(i,j) ,
Dj (q)
zq =
X
X
(i,j)∈δ + (S ) d=1
(10)
x(i,j) ,
The balance equation over S with aggregate variables becomes:
aqd ∈δ − (i)
X
q.xqd a = d(S ).
+ aqd ∈δQ (S )
Let yq and zq be aggregated variables defined as follows:
(9)
r ∈R
x00 (i,0) = m
(11)
(i,0)∈δ − (0)
λr ≥ 0 ∀r ∈ R xqd a
X
q.xqd a −
(i,j)∈δ − (S ) d=1
d.
(20)
Dj (q)
ca .xqd a
X
pqd caqd . ar .¯
aqd ∈AQ
aqd ∈AQ
s.t.
(15)
i=1
can be viewed as defining a network flow of m units over an acyclic layered graph GN = (VN , AN ). For each vertex i ∈ V+ there are Q vertices in VN , one for each possible load size q = 0, . . . , Q − 1. The depot corresponds to two additional vertices in VN . Each arc aqd ∈ AQ corresponds to an arc in AN . As the flow of m units in GN has just one source and one destination (the vertices associated to the depot) any integer solution can be decomposed into a set of m paths, each one corresponding to a route and a delivery pattern. Constraints (2) impose that each client receives exactly its demand. The direct use of this formulation is not practical even for small sized instances, so we provide an equivalent reformulation suitable to be solved by column generation. Let R denote the set of origin– destination paths in GN , each one corresponding to a possible route and delivery pattern. For each r ∈ R, let λr be a continuous variqd able and let par indicate whether a vehicle executing r traverses arc a = (i, j) with q units of goods and delivers exactly d units to j. Minimize
∀ i ∈ V+
aqd ∈δ − (i)
aqd ∈AQ
(i,0)∈δ − (0)
xqd a ∈ Z+
.λr = di d.pqd ar
r ∈R
aqd ∈δQ (i) D i ( Q −q )
X
r ∈R
(2)
(14)
Q
−
X
.λr ca .pqd ar
X X
xqd a
X
aqd ∈AQ
X
303
∈ Z ∀a ∈ AQ . qd
(12) (13)
It can be seen that Equalities (3) and (4) are implied by the definition of R and by (9). The so-called Dantzig–Wolfe Master problem (DWM) is obtained with the relaxation of integrality constraints
Q X q =1
q.yq −
Q −1 X
q.zq = d(S ).
(24)
q =1
For each possible values of Q and D = d(S ), define P (Q , D) as the polyhedron induced by the integer solutions of (24). The Homogeneous Extended Capacity Cuts (HECCS) [13] are all inequalities valid for those polyhedra. The HECCs, already used on other problems [11,12], are also valid for the SDVRP.
304
L. Moreno et al. / Operations Research Letters 38 (2010) 302–306
2.1.2. Edge upper bound
3.1. Column generation
Dror and Trudeau [6] stated the following theorem: Theorem 1. If distance matrix satisfies triangular inequality, there is an optimal solution to the SDVRP where any two paths have at most one client in common. In other words, this theorem assures that there is an optimal solution where two clients appear on the same path at most once. A weaker property also valid for SDVRP is: Corollary 1. If distance matrix satisfies triangular inequality, there is an optimal solution where, for any two vertices i, j ∈ V + , i 6= j, the edge (i, j) of the undirected graph is traversed at most once. Thus, Eq. (25) is a valid inequality for SDVRP. Dj (q) Q X X q =1
d=1
Di (q) qd x(i,j)
+
X
! qd x(j,i)
≤ 1,
∀i, j ∈ V+ , i 6= j.
(25)
d=1
The requirement of a distance matrix satisfying the triangle inequality in the above propositions is actually not restrictive. A preprocessing procedure can always be used to obtain such a distance matrix. If the direct distance between clients i and j, c(i,j) , is greater than the distance from i to another client k and then to j, the value of c(i,j) can be replaced by c(i,k) + c(k,j) . This transformation is valid because the arc (i, j) can now be interpreted as a vehicle travelling from i to k (and delivering nothing) and then travelling from k to j. We applied this transformation to all instances in our tests, so inequality (25) can always be used.
2.1.3. Split deliveries cuts Eqs. (2) impose that the sum of all deliveries at vertex i must be equal to the amount requested by i. Valid inequalities of the polyhedron induced by the integer solutions of (2) are referred to as Split Deliveries Cuts (SDCs). These inequalities are valid for SDVRP.
The generation of minimum reduced cost paths is a NP-hard problem that can be solved in pseudo-polynomial time. Note that the SDVRP definition allows paths that visit the same client more than once, i.e., cycling is permitted. The generation of paths can be executed by a dynamic programming algorithm (DP) with O(|V |.Q ) states and pseudo-polynomial complexity O(|A|.Q 2 ). When triangular inequality holds, there is always an optimal solution in which no path has a cycle. Thus, we eliminate paths with 2-cycles (subpaths i → j → i, i, j ∈ V +) in order to strengthen the formulation. As already noted by [4], this elimination of 2-cycles does not change the pricing complexity. In order to accelerate the pricing, two dynamic programming based heuristics (suggested by Fukasawa et al. [7]) are also used. The scaling heuristic consists of performing the dynamic programming over a simplified instance obtained by dividing the capacity by an integer g > 1 (and rounding down) and also dividing all demands by the same g (but rounding up). A path with negative reduced cost obtained in this simplified instance corresponds to a path with negative reduced cost in the original instance. The sparsification heuristic restricts the set of arcs used in the dynamic programming, for each client i only the k cheaper arcs incident to it are considered. Of course, the dynamic programming based heuristics can fail to find a minimum reduced cost path. Hence, column generation stops only when the pricing without any acceleration technique is used and no path with negative reduced cost is found. 3.2. Cutting At each cut generation iteration, the current fractional solution qd of formulation (14)–(17) is converted into xa variables using Eqs. (9). Violated cuts are then separated and introduced into the DWM LP in format (18). Inequalities (25) are separated by inspection, while Homogeneous Extended Capacity Cuts are separated by the same procedure described in [11]. The separation of Split Delivery Cuts begins by discarding all clients that do not have split deliveries in the current solution. Then, we look for violated inequalities of the form:
X
drdeaqd = drdi e,
(26)
−
aqd ∈δQ (S )
3. Lower bounding algorithm The lower bounding algorithm combines column and cut generation over the DWM relaxation presented in Section 2. That linear program is initialized with a set of variables (paths) constructed as follows:
P −1 • The first path visits clients 1, 2, . . . , j such that ij= 1 di < Pj Q and i=1 di ≥ Q . It delivers all the demand of clients Pj−1 1, 2, . . . , (j − 1) and Q − i=1 di units of client j. The next path transports the remaining units of j and all demands of the following clients until the vehicle capacity is achieved again, and so on.
• There is also |V+ | additional paths, each one delivering all demand of a single client. At the beginning, the DWM only contains the constraints from the original formulation, (15)–(17). Since the linear relaxation of a set covering problem is easier to solve by column generation than a set partitioning problem, Eqs. (15) are relaxed to ≥. This change has no impact on the lower bounds obtained.
where 0 < r < 1. The example below shows a fractional solution and a violated inequality. Suppose that the demand of a client i is 7 and that its demand is covered by 3 distinct paths in a fractional solution. Each path is associated to a variable λr with value 0.5. The first and second paths are cycle-free and deliver 6 and 4 units to client i respectively. The third path visits the client twice and delivers 2 units at each stop. Relaxing equation (26) to ≥ and using r = 0.5, a violated inequality is found: 4w7 + 3w6 + 3w5 + 2w4 + 2w3 + w2 + w1 ≥ 4.
(27)
4. Experiments The tests were executed in a machine with processor Intel Core 2 Duo 2 GHz and with 3 GB of RAM. CPLEX 11.0 was used to solve the linear programs. This algorithm was tested for 3 distinct sets of instances found in the SDVRP Pnliterature. In all those sets the number of vehicles m is set as d i=1 di /Q e. The first set was generated by Belenguer et al. [2] and has 14 instances with 50, 75 or 100 clients and vehicle capacity of 160 units. Instances with name ending in ‘‘d1’’ have demands ranging between 1% and 10% of vehicle capacity. Instances with names ending
L. Moreno et al. / Operations Research Letters 38 (2010) 302–306 Table 1 Comparison with [2] on their set of instances.
305
Table 3 Results on the set of instances by [3].
Instance
UB
LB [2]
Our LB
Our time
Instance
n
UB [3]
Our LB
Our time
s51d1 s51d2 s51d3 s51d4 s51d5 s51d6 s76d1 s76d2 s76d3 s76d4 s101d1 s101d2 s101d3 s101d5
458 726 972 1677 1440 2327 594 1147 1474 2257 716 1393 1975 2915
454.0 676.6 905.2 1520.7 1272.9 2113.0 584.9 1020.3 1346.3 2011.6 700.6 1271.0 1739.7 2630.4
454.4 694.2 930.7 1539.0 1313.4 2141.7 586.6 1061.1 1395.9 2046.1 704.9 1337.1 1832.2 2737.1
88 123 200 609 369 525 125 393 557 1438 753 1025 1683 1168
8.6%
5.9%
sd1 sd2 sd3 sd4 sd5 sd6 sd7 sd8 sd9 sd10 sd11 sd12 sd13 sd14 sd15 sd16 sd17 sd18 sd19 sd20 sd21
8 16 16 24 32 32 40 48 48 64 80 80 96 120 144 144 160 160 192 240 288
228.3 714.4 430.6 631.1 1 408.1 831.2 3 714.4 5 200.0 2 059.8 2 749.1 13 612.1 7 399.1 10 367.1 11 023.0 15 271.8 3 449.1 26 665.8 14 546.6 20 559.2 40 408.2 11 491.7
228.3 708.3 430.6 631.0 1 390.6 831.2 3 640.0 5 068.3 2 044.2 2 684.9 13 275.0 7 175.8 10 053.6 10 588.2 14 908.5 3 381.3 26 317.2 14 029.2 19 707.2 39 252.8 11 271.0
0 1 0 2 6 0 56 92 11 70 431 748 2045 165 2703 6 2192 3060 723 3785 205
Avg. gap
Table 2 Comparison with [9] on the set of instances by [2]. Instance s51d1 s51d2 s51d3 s51d4 s51d5 s51d6 s76d1 s76d2 s76d3 s76d4 s101d1 s101d2 s101d3 s101d5 Avg. gap
UB
LB [9]
Our LB
723.4 968.9 1586.5 1355.5 2197.8
695.0 922.7 1505.4 1297.5 2108.6
1185.7 1504.9 2136.4
1066.2 1397.4 2019.9
1474.5 2012.9 2846.2
1349.8 1837.3 2725.5
457.1 700.4 938.5 1549.7 1318.9 2154.7 592.6 1071.3 1407.5 2059.8 716.8 1358.9 1853.1 2767.6
6.5%
5.0%
Time [9] (s) 3555 323 14 2 <1 7511 549 188 16 507 2698 10
Our time 184 131 225 631 320 529 146 525 771 1554 671 1326 1068 1670
in ‘‘d2’’ to ‘‘d6’’ have demands ranging between the following fractions of the vehicle capacity: (10%–30%), (10%–50%), (10%–90%), (30%–70%) and (70%–90%). Tables 1 and 2 compare our results on this first set of instances with those obtained by Belenguer et al. [2] using a cutting-plane algorithm and with those obtained by Jin et al. [9] using column generation. While [9] compute euclidean distances without rounding, as double precision values, [2] round them to the nearest integer using the TSPLIB convention. Therefore, we show separate results considering rounded distances and fractional distances. It can be seem that improved bounds were obtained in all those instances. All times presented in the tables are in seconds. The times of [9] were obtained in a machine with a processor Pentium 4, 2 GHz and 2 GB memory. Although their algorithm can be significantly faster than ours on instances with large demands, it can also be significantly slower on instances with small demands (they could not run instances of type ‘‘d1’’). The newly proposed lower bounding algorithm is much more stable in terms of running time. Article [2] did not provide times for comparison. The second set has 21 instances created by Chen et al. [3] (Table 3). In all problems, vehicle capacity is 100 and client demand is either 60 or 90. These instances are highly symmetric, with clients placed in concentric circles around the depot. The euclidean distances are not rounded. Those authors used a heuristic algorithm to provide upper bounds. The average integrality gap for these instances is 1.6%. Lower bound values in bold are optimal integer solutions, they were obtained for 9 instances of this set. The largest solved instance has 288 clients. The third set of instances is composed by some instances from the TSPLIB (Table 4). The results are compared with the lower bounds given by Belenguer et al. [2]. The first 4 instances of this
Avg. gap
1.6%
Table 4 Comparison with [2] on TSPLIB instances. Instance
n
Q
Q0
UB
LB [2]
eil22 eil23 eil30 eil33 eil51 eila76 eilb76 eilc76 eild76 eila101 eilb101
21 22 29 32 50 75 75 75 75 100 100
6000 4500 4500 8000 160 140 100 180 220 200 112
60 4500 180 800 160 140 100 180 220 200 112
375 569 510 835 521 832 1023 735 683 817 1077
375.0 569.0 508.0 833.0 511.6 782.7 937.5 706.0 659.4 793.5 1005.9
375.0
3
510.0 834.7 517.8 807.6 981.4 717.8 666.1 799.8 1040.6
157 1065 92 101 80 124 167 632 249
3.2%
1.6%
Avg. gap
Our LB
Our time
set have very large vehicle capacities Q . Since the complexity of our pricing subproblem is O(|A|.Q 2 ), we compute the greatest common divisor g between the vehicle capacity and each demand value. Then, each demand value di can be replaced by d0i = di /g and the vehicle capacity Q is reduced to Q 0 = Q /g. We could not run our algorithm in instance eil23 as Q 0 is still 4500. Optimal integer solutions were obtained for instances eil22 and eil30. Although the solution of instance eil33 is fractional, the lower bound of 834.7 is enough to prove the optimality of a previously known solution with value 835. The algorithm proposed is quite different from the previous algorithms to the SDVRP. The extended formulation provides a more detailed view of the problem and a more convenient structure for the identification of effective inequalities. In particular, by using extended variables, the knapsack aspect of the problem could be better explored. The algorithm is able to achieve sharper lower bounds for almost all SDVRP instances tested. However, it has a limitation: it is not practical for problems with very large vehicle capacity values. In these cases, the algorithm becomes quite slow and can even fail to complete due to memory requirements. Acknowledgements LM and MPA were partially financed by CNPq Grants 140715/ 2004-5 and 311997/06-6. EU received support from CNPq Grant 304533/02-5 and Faperj Grant E-26/100.504/2007. The authors also thank the technical support by GAPSO Inc., Brazil.
306
L. Moreno et al. / Operations Research Letters 38 (2010) 302–306
References [1] C. Archetti, M.W.P. Savelsbergh, M.G. Speranza, Worst-case analysis for Split Delivery Vehicle Routing Problems, Transportation Science 40 (2) (2006) 226–234. [2] J.M. Belenguer, M.C. Martinez, E. Mota, A lower bound for the Split Delivery Vehicle Routing Problem, Operations Research 48 (5) (2000) 801–810. [3] S. Chen, B.L. Golden, E.A. Wasil, The Split Delivery Vehicle Routing Problem: applications, algorithms, test problems, and computational results, Networks 49 (4) (2007) 318–329. [4] N. Christofides, A. Mingozzi, P. Toth, Exact algorithms for the vehicle routing problem based on spanning tree and shortest path relaxations, Mathematical Programming 20 (1981) 255–282. [5] M. Dror, G. Laporte, P. Trudeau, Vehicle routing with Split Deliveries, Discrete Applied Mathematics 50 (3) (1994) 239–254. [6] M. Dror, P. Trudeau, Savings by Split Delivery Routing, Transportation Science 23 (1989) 141–145. [7] R. Fukasawa, H. Longo, J. Lysgaard, M. Poggi de Aragão, M. Reis, E. Uchoa, R. Werneck, Robust branch-and-cut-and-price for the capacitated vehicle routing problem, Mathematical Programming 106 (3) (2006) 491–511.
[8] M. Jin, K. Liu, R.O. Bowden, A two-stage algorithm with valid inequalities for the Split Delivery Vehicle Routing Problem, International Journal of Production Economics 105 (1) (2007) 228–242. [9] M. Jin, K. Liu, B. Eksioglu, A column generation approach for the Split Delivery Routing Problem, Operations Research Letters 36 (2) (2008) 265–270. [10] C. Lee, M.A. Epelman, C.C. White III, Y.A. Bozer, A shortest path approach to the multiple-vehicle routing with split pick-ups, Transportation Research Part B: Methodological 40 (4) (2006) 265–284. [11] A. Pessoa, E. Uchoa, M. Poggi de Aragão, A robust branch-cut-and-price algorithm for the heterogeneous fleet vehicle routing problem, Networks 54 (2009) 167–177. [12] A. Pessoa, E. Uchoa, M. Poggi de Aragão, R. Rodrigues, Algorithms over arc-time indexed formulations for single and parallel machine scheduling problems, Technical Report RPEP, vol. 8, No. 8, Universidade Federal Fluminense, Engenharia de Produção, Niterói, Brazil, 2008. [13] E. Uchoa, R. Fukasawa, J. Lysgaard, A. Pessoa, M. Poggi de Aragão, D. Andrade, Robust branch-cut-and-price for the capacitated minimum spanning tree problem over a large extended formulation, Mathematical Programming 112 (2) (2007) 443–472.