- Email: [email protected]
A column generation approach for the split delivery vehicle routing problem
Operations Research Letters 36 (2008) 265 – 270
Operations Research Letters www.elsevier.com/locate/orl
A column generation approach for the split delivery vehicle routing problem Mingzhou Jin ∗ , Kai Liu, Burak Eksioglu Department of Industrial and Systems Engineering, Mississippi State University, P.O. Box 9542, MS 39762, USA Received 22 December 2005; accepted 18 May 2007 Available online 10 July 2007
Abstract A column generation approach is presented for the split delivery vehicle routing problem with large demand. Columns include route and delivery amount information. Pricing sub-problems are solved by a limited-search-with-bound algorithm. Feasible solutions are obtained iteratively by fixing one route once. Numerical experiments show better solutions than in the literature. © 2007 Elsevier B.V. All rights reserved. Keywords: Split delivery vehicle routing; Column generation; Integer programming
1. Introduction The split delivery vehicle routing problem (SDVRP) was introduced to the vehicle routing problem (VRP) literature by [4] in the late 1980s. Unlike the capacitated vehicle routing problem (CVRP), the SDVRP allows the delivery to a demand point to be split between two or more vehicles. In many cases, allowing split deliveries results in significant savings both in the total distance traveled and the number of required vehicles compared to the corresponding CVRP solution. An example in [9] shows that the ratio of the optimal value of the SDVRP over that of the corresponding CVRP approaches 1/2 when the number of demand points goes to +∞. The savings are, in general, more significant when average customer demand is large (e.g., more than 10% of the vehicle capacity) [3]. There are a few exact algorithms in the literature to solve small-size SDVRP problems. A dynamic programming (DP) model with infinite state and action spaces is developed in [12]. The largest instance solved by [12] has nine demand points and six vehicles. A two-stage algorithm with valid inequalities is developed by [11] to successfully solve instances with up to 21 points. A cutting-plane algorithm incorporating the valid inequalities in [5] is used to solve small instances with up to 10 demand points optimally. In [2] a polyhedral study of the SDVRP is given and their cutting-plane algorithm with several valid inequalities yields good lower bounds and upper bounds ∗ Corresponding author.
E-mail address: [email protected] (M. Jin). 0167-6377/$ - see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.orl.2007.05.012
for large instances with up to 100 demand points. The gaps are around 10% for instances with 50, 75, or 100 demand points. The numerical experiments presented in Section 6 use the instances provided by [2] as benchmarks to evaluate the performance of the proposed column generation (CG) approach. A Tabu search algorithm is developed in [1] for the SDVRP. Various practical problems are modeled as SDVRPs in the literature, and they are generally solved using heuristic algorithms. Ref. [7] and [8] study the SDVRP and its variation with time windows on a grid network and propose constructive heuristics. A heuristic algorithm is used in [13] for a livestock feed distribution problem encountered on a cattle ranch in Arizona. In [14] an algorithm based on CG techniques is provided along with other heuristics to determine a flight schedule for helicopters to off-shore platform locations for exchanging work crews in the North Sea. The flight schedule does not restrict the number of total trips but has a constraint on the number of total stops in a trip. A CG approach is developed for solving the SDVRP with time windows in [9]. The CG approach proposed in this paper for the general SDVRP is similar to the algorithm, developed for one specific application, presented in [14]. However, our approach is different from that of [14] in the following ways: (1) our approach requires the solutions to use a fixed number of vehicles but does not have a restriction on the longest distance of each trip; (2) the decision variables in our master problem are defined as binary variables rather than general integer variables in order to improve the lower bound; (3) the algorithm to obtain the feasible
266
M. Jin et al. / Operations Research Letters 36 (2008) 265 – 270
solution fixes one route with the largest product of variable value and its total delivery amount and guarantees feasibility over iterations; and (4) we conduct numerical experiments on the instances published in [2] and show that the CG approach provides good results for instances with large average demand. 2. The SDVRP definition and model The SDVRP can be represented by a graph G = (V , E), where V = {0, 1, . . . , N) is the set of vertices and E = {(i, j ) : i = j } is the set of edges. Vertex 0 denotes the depot at which a fleet of identical vehicles with capacity Q are located. The other N vertices of V correspond to demand points. A positive demand di is associated with each demand point i (d0 = 0). A positive weight cij is associated with each edge (i, j ) ∈ E and represents the distance between vertices i and j. The distance from each vertex to itself is assumed to be a large positive number (cii = T , i = 0, . . . , N and T = a large number). The distances are assumed to be symmetric (cij = cj i ) and satisfy the triangular inequality (cij +cj k cik ). The SDVRP is to find routes such that the total distance traveled is minimized and the following conditions are satisfied: • The demand di of demand point i is delivered by one or more vehicles. • All vehicle routes start and end at the depot. • The total demand serviced by any vehicle does not exceed its capacity. In [3] and [12] mixed integer programming (MIP) formulations of the SDVRP are provided by fixing the total number of vehicles at K = N i=1 di /Q. This paper also requires to use K vehicles. This requirement implies that the delivery amount of each vehicle should be at least N i=1 di − (K − 1)Q. 3. The CG formulation and the column generation algorithm to obtain a lower bound (CGALB) The proposed CG approach is composed of a master problem M and a pricing sub-problem P, which are given below following the notations. The notations for the master problem M: air :
The amount delivered to demand point i by route r (in other words, a.r represents route r with delivery quantities); : The set of columns (routes); r : The minimum distance traveled to visit demand points with positive demand (air > 0) on route r; A binary variable that represents whether or not route br : r is chosen in the solution. The master problem: r br (1) M: min Z M = s.t.
r∈
air br di ,
i = 1, . . . , N,
(2)
r∈
r∈
A column r, which has N + 1 entries, is defined by a.r plus a “1”. Note that the definition of a column in the master problem M is different compared to the column definitions in CG approaches developed for routing problems without split delivery because a column in M includes delivery information. Thus, M is not the traditional set-partitioning or set-covering problem. The objective of M (1) is to minimize the total distance traveled. Constraints (2) ensure that the demand is satisfied at all points. Rather than “=” constraints in the model provided by [14], we use “ ” constraints because numerical experiments show that “ ” constraints improve the speed of the proposed CG approach. Constraint (3) guarantees that exactly K columns are used. Unlike in [14], br are defined as binary variables rather than integer variables because of the following theorem regarding the optimal solution: Theorem 1 (Dror and Trudeau [5]). Given cij is the distance between point i and point j. If the {cij } matrix satisfies the triangular inequality, then no two routes in the optimal solution to M can have more than one split demand point in common. The optimal solution should not have two identical columns that visit two or more than two points or have total delivery amount less than Q. Let be the set of all possible columns (i.e., routes with delivery quantities) that are different from each other and di /Q duplicates of column i for any demand point i for which di Q. In column i the ith entry is Q, the (N + 1)th entry is “1”, and all other entries are zero. To guarantee that are used all columns in should exactly K vehicles N satisfy N i=1 air i=1 di − (K − 1)Q. With this definition of , defining br as binary variables rather than general integer variables can improve the lower bound provided by the linear programming (LP) relaxation of M. When = , M yields an optimal solution. It is computationally intractable to enumerate all possible routes with delivery amounts. Instead, we use a CG approach to create routes using information obtained from the LP relaxation of M. Let i (i = 1, . . . , N ) be the dual price of the ith constraint for the current set of columns in the LP relaxation of M featured by . Let be the dual price of constraint (3) (the (N + 1)th constraint) in the LP relaxation of M featured by a given . The following pricing sub-problem P uses the dual prices to generate a new column: P: min
Z = P
br ∈ {0, 1}.
(3)
N
cij xij −
i=0 j =0
s.t.
N
x0j =
j =1 N j =1
i ai −
(4)
i=1
N
xj 0 = 1,
(5)
j =1
xij =
N
xj i = yi ,
i = 1, . . . , N,
(6)
j =1
ai di yi , i = 1, . . . , N, N N di − (K − 1)Q ai Q, i=1
br = K,
N N
i=1
ui − uj + (N + 1)xij N, xij , yi ∈ {0, 1}, ai , ui 0.
(7) (8)
i, j = 1, . . . , N, (9)
M. Jin et al. / Operations Research Letters 36 (2008) 265 – 270
267
bound (ICGAUB) is developed in Section 5 by fixing one column (route with delivery) iteratively. 4. LSB algorithm for the pricing sub-problem
Fig. 1. The CG algorithm for obtaining a lower bound (CGALB).
The pricing sub-problem P generates a route by minimizing the total cost (i.e., total distance traveled minus the benefit from delivery). The route starts from and ends at the depot (Constraints (5)). The model has two sets of binary variables: xij = 1 indicates that vehicle (route) r visits point j right after point i, and yi = 1 indicates that the vehicle visits demand point i (Constraints (6)). At each demand point i, the delivery amount ai should not exceed the demand (Constraints (7)). The total delivery of the route cannot exceed the vehicle capacity but should be at least N i=1 di − (K − 1)Q to make the column feasible for the SDVRP (Constraint (8)). Constraints (9) are sub-tour elimination constraints. If the solution to the pricing sub-problem P results in Z P < 0, then {ai , i = 1, . . . , N} for the first N constraints and a “1” for the (N + 1)th constraint are added to M as a new column in the next iteration. On the other hand, if the solution to P results in Z P 0, the current optimal value of the LP relaxation of M featured by is the same as the optimal value of the LP relaxation of M featured by , which is a lower bound to the SDVRP. The CG algorithm for obtaining the lower bound (CGALB) is illustrated in Fig. 1. The initial set of columns created by Step 0 of the CGALB should have a feasible solution to continue. There are exactly K columns in the initial set. The first vehicle (column) visits points 1 −1 1 i = 1, . . . , t1 , such that ti=1 di < Q and ti=1 di Q and the delivery amounts are ai1 = di for i = 1, . . . , t1 − 1 and at1 ,1 = 1 −1 di . Obviously, a solution of {xi = 1, i = 1, . . . , K} is Q − ti=1 feasible to the model M featured by the initial set of columns. In Step 2, the LP relaxation of M is solved using ILOG CPLEX 9.0 [10]. In Step 3, the pricing sub-problem P is solved, which is similar to the traveling salesman problems (TSPs) with profits in [6]. However, the TSPs with profits have a fixed profit from visiting one demand point while the benefit in the pricing sub-problem P depends on the delivery amount. Furthermore, the pricing sub-problem P has the vehicle capacity and minimal delivery constraints. The pricing sub-problem P, like the TSPs with profits, cannot be solved in a reasonable amount of time as an integer program by existing commercial solvers. The gap is still larger than 50% after 2 h when ILOG CPLEX 9.0 is used to solve a pricing sub-problem P with 50 demand points (e.g., instance S51D2 in the numerical experiments). Therefore, we develop an algorithm called limited-search-with-bounds (LSB) in Section 4 to solve the pricing sub-problem. Numerical experiments show that M featured by the set of columns that are created by the CGALB cannot be solved as an IP directly because of memory overflow and long computational time. Therefore, the iterative Column Generation approach to obtain an upper
The pricing sub-problem P determines the route and the delivery amounts for one or more vehicles. The benefit at demand point i is proportional to the delivery amount (ai ) and the dual price (i ). Thus, the delivery amount in an optimal solution satisfies the following Lemma. Lemma 2. In an optimal solution to the pricing sub-problem, each visited demand point receives its demand (ai = di ) except for one visited point with the smallest dual price that may receive partial (less-than demand) delivery. Lemma 2 is straightforward because ai appears only in Constraints (7) and (8). Without loss of generality, assume that 1 2 · · · n . Let TSP(I) denote the minimum distance traveled to visit all demand points in set I ⊆ {1, 2, . . . , N} and the depot. Let e(I ) be the largest demand point index in I. Lemma 3. If the route visits a set of demand points I plus any other demand points whose indices are larger than e(I), a lower bound to the pricing sub-problem is e(I )+J TSP(I ) − i∈I i di − i=e(I )+1 i di − e(I )+J +1 (Q − e(I )+J i∈I di − i∈I di . Here J is i=e(I )+1 i di ) − , when Q > e(I )+J a nonegative integer such that i∈I di + i=e(I )+1 i di Q e(I )+J +1 and i∈I di + i=e(I )+1 i di > Q. All terms with an index larger than N will not be included in this summation. Lemma 3 is based on three facts: (1) visiting more demand points increases the total distance traveled because of the triangular inequality assumption; (2) the dual prices i are in a descending order; and (3) the route has at most one demand point that is partially served. The following LSB algorithm is developed based on Lemmas 2, 3, and the capacity constraint (Fig. 2). Note that the average demand at each demand point is typically larger in an SDVRP compared to a CVRP. Therefore, the TSP in Step 1 of the LSB algorithm can be solved relatively easily by using CPLEX 9.0. Set I is fathomed and does not need any further exploration when (1) there is no more capacity; (2) the set includes the last demand point that has the smallest dual price; or (3) L(I ) U , where L(I ) is a lower bound for the solution including the demand points in I and other demand points with indices larger than e(I ) (i.e., smaller or the same dual price). If set I is not fathomed then it is explored by adding the next demand point e(I )+1 to I in Step 4. If set I is fathomed, then backtracking is performed in Step 5 until e(I ) = N and I = ∅. If e(I ) = N , then e(I ) is replaced with e(I ) + 1. When I =∅, the algorithm stops. With the LSB algorithm, all columns created by P have at most one demand point whose delivery is split. In other words, the upper bound obtained by the proposed CG is bounded from below by the master problem M with = {all columns with at most one split delivery point}.
268
M. Jin et al. / Operations Research Letters 36 (2008) 265 – 270
Fig. 2. The LSB algorithm.
5. The iterative column generation approach to obtain an upper bound (ICGAUB) When no more columns are created, the solution to the linear relaxation of M featured by the created column set yields a lower bound for the SDVRP. Usually the solution has fractional numbers so it is not feasible to the SDVRP. Solving the original M with br defined as binary variables can provide a feasible solution. However, this is very time consuming. Based on numerical experiments, when there are more than 50 demand points, using optimization solvers to solve the original M featured by created by CGALB is not practical. In fact, after 24 h, ILOG CPLEX 9.0 cannot even find a feasible solution though a feasible solution definitely exists because of the initial set of columns created by Step 0 in the CGALB. We propose the following ICGAUB in order to find a feasible solution by fixing one route with delivery amounts at each iteration. In the first iteration, after solving the LP relaxation of M for the original SDVRP, the algorithm fixes the variable in the optimal solution with the largest product of xr N i=1 air to 1. The rule of selecting which column (route with delivery amounts) should be fixed is based on the following two conjectures: 1. The columns with higher variable values in the solution to the LP relaxation of M are “good” routes. 2. Less demand left for the reduced SDVRP obtained after fixing the column can improve the solution to the reduced SDVRP. We conducted numerical experiments using other rules, but the above given rule was better than others. Fixing the column
means that the points whose demand are fully satisfied (i.e., air = br ) should be removed. All columns that visit any of these fully satisfied points should also be removed from the column set. The demand of the point that is partially served (i.e., 0 < air < br ) is reduced by the served amount. These operations (Step 3 in Fig. 3) create a new and smaller SDVRP with |{i : air∗ = di }| fewer demand points and 1 less vehicle. The new SDVRP is solved by using the CG algorithm in Fig. 1 again. Since the minimal delivery requirement for each route of the new SDVRP may be smaller than that of the SDVRP in the previous iteration, some of the old columns may violate the new minimal delivery requirement and should be removed. Note that Step 0 in the CGALB is still necessary to make the first LP for the reduced SDVRP feasible. Because all columns at each iteration meet the minimal delivery requirement for the reduced SDVRP, the reduced SDVRP at iteration l always has a feasible solution by using K −l +1 vehicles. When the reduced SDVRP is small enough (e.g., the product of the number of remaining points and the number of required vehicles is less than a threshold value Th = 80), the exact solution approach introduced in [11] is used to obtain an optimal solution. An upper bound for the original SDVRP is obtained after adding back the costs of all fixed columns. 6. Numerical experiments A desktop computer with a Pentium 4 processor, 2.8 GHz CPU and 2 G memory is used to conduct numerical experiments. The results are compared to the numerical results of [2]. In [2] 25 randomly generated instances are solved. However,
M. Jin et al. / Operations Research Letters 36 (2008) 265 – 270
269
Fig. 3. The iterative column generation approach to obtain an upper bound (ICGAUB).
Table 1 Comparison of the CG Algorithm with the cutting-plane algorithm Instance name
S51D2 S51D3 S51D4 S51D5 S51D6 S76D2 S76D3 S76D4 S101D2 S101D3 S101D5 eil76B
K
9 15 27 23 41 15 23 37 20 31 48 14
Cutting plane
Column generation
UB
LB
Gap %
UB
LB
Gap %
T1
T2
726 972 1677 1440 2327 1147 1474 2257 1393 1975 2915 1023
676.63 905.22 1520.67 1272.86 2113.03 1020.32 1346.29 2011.64 1270.97 1739.66 2630.43 937.47
6.80 6.87 9.32 11.61 9.20 11.04 8.66 10.87 8.76 11.92 9.76 8.36
723.37 968.85 1657.61 1439.92 2300.21 1185.72 1504.94 2219.07 1474.51 2012.86 2954.96 1063.75
694.98 922.72 1505.35 1297.46 2108.59 1066.17 1397.43 2019.91 1349.77 1837.33 2725.5 981.14
3.92 4.76 9.19 9.89 8.33 10.08 7.14 8.97 8.46 8.72 7.77 7.77
3555 323 14 2 <1 7511 549 188 16 507 2698 10 30 222
5978 607 260 46 243 12 806 2030 1813 47 658 7959 847 45 084
The underlined numbers show the results where the CG approach was worse compared to the cutting-plane approach.
we only include 12 of those 25 instances in our analysis. The included instances, listed in Table 1, are the ones in which average demand at a point is about 15% larger than the vehicle capacity. This is done because the CVRP rather than the SDVRP is a more appropriate model for instances with small average demand. We note that the proposed CG approach is not efficient for solving instances with small average demand due to the resulting large
tree in the LSB algorithm and the large computational time for solving the embedded TSPs. The full set of 25 instances is available at www.uv.es/∼belengue/SDVRP/sdvrplib.html. In order to obtain a better solution and speed up the algorithm, rather than adding only the best column, other “good” columns out of the LSB algorithm are also added to the master problem M. If the objective function value corresponding to set
270
M. Jin et al. / Operations Research Letters 36 (2008) 265 – 270
I is negative and lower than a threshold value, then the column corresponding to set I is added to M. In Table 1, column K shows the minimum number of vehicles required. The results of the cutting-plane method are provided by [2]. In the results of the CG approach, T 1 shows the CPU time (in seconds) spent by the CGALB to obtain the lower bound and T 2 shows the total time (in seconds) spent by the ICGAUB to obtain the upper bound. The CG algorithm outperforms the cutting-plane method in all instances with respect to the error gap. The CG algorithm finds better upper bounds in six instances and better lower bounds in ten instances out of the twelve instances. In general, the CG approach takes much less time when the average demand at each point is larger. Similar phenomenon is observed in the CG approach developed for the SDVRP with time windows in [9]. This result is interesting because most of other algorithms (e.g., [2,12]) developed for the SDVRP require more time as average demand increases for problems with the same number of demand points, which is perhaps caused by the increase in the number of binary variables. In our CG algorithm, larger average demand means more efficient implementation of the LSB algorithm because of the capacity constraint. Furthermore, larger demand instances have fewer column candidates that have at most one split delivery (i.e., fewer iterations of the CG algorithm). 7. Conclusion To improve the lower bounds, cuts can be added to the problem. For example, the number of vehicles visiting a subset of demand points should be enough for serving the total demand of the points in the subset. Though it is easy to generate the cuts, more studies are necessary to continue the CG approach
by modifying the algorithm to solve the pricing sub-problem with the cuts. References [1] C. Archetti, A. Hertz, M.G. Speranza, A Tabu search algorithm for the split delivery vehicle routing problem, Transp. Sci. 40 (2006) 64–73. [2] J.M. Belenguer, M.C. Martinez, E. Mota, A lower bound for the split delivery vehicle routing problem, Oper. Res. 48 (2000) 801–810. [3] M. Dror, G. Laporte, P. Trudeau, Vehicle routing with split deliveries, Discrete Appl. Math. 50 (1994) 239–254. [4] M. Dror, P. Trudeau, Savings by split delivery routing, Transp. Sci. 23 (1989) 141–145. [5] M. Dror, P. Trudeau, Split delivery routing, Naval Res. Logist. 37 (1990) 383–402. [6] D. Feillet, P. Dejax, M. Gendreau, Traveling salesman problems with profits, Transp. Sci. 39 (2005) 188–205. [7] P.W. Frizzell, J.W. Giffin, The bounded split delivery vehicle routing problem with grid networks distances, Asia Pacific J. Oper. Res. 9 (1992) 101–116. [8] P.W. Frizzell, J.W. Giffin, The split delivery vehicle scheduling problem with time windows and grid network distance, Comput. Oper. Res. 22 (1995) 655–667. [9] M. Gendreau, D. Feillet, P. Dejax, C. Gueguen, Vehicle routing with time windows and split deliveries, Technical Paper 2006-851, Laboratoire d’Informatique d’Avignon, 2006. [10] ILOG, ILOG CPLEX 9.0 User’s Manual, 2003. [11] M. Jin, K. Liu, R.O. Bowden, A two-stage algorithm with valid inequalities for the split delivery vehicle routing problem, Int. J. Production Econ. 105 (2007) 228–242. [12] C.G. Lee, M. Epelman, C.C. White, A shortest path approach to the multiple-vehicle routing with split picks-up, Transp. Res.—Part B 40 (2006) 265–284. [13] A. Mullaseril, M. Dror, J. Leung, Split-delivery routing heuristics in livestock feed distribution, J. Oper. Res. Soc. 48 (1997) 107–116. [14] G. Sieksma, G.A. Tijssen, Routing helicopters for crew exchanges on off-shore locations, Ann. Oper. Res. 76 (1998) 261–286.
Operations Research Letters www.elsevier.com/locate/orl
A column generation approach for the split delivery vehicle routing problem Mingzhou Jin ∗ , Kai Liu, Burak Eksioglu Department of Industrial and Systems Engineering, Mississippi State University, P.O. Box 9542, MS 39762, USA Received 22 December 2005; accepted 18 May 2007 Available online 10 July 2007
Abstract A column generation approach is presented for the split delivery vehicle routing problem with large demand. Columns include route and delivery amount information. Pricing sub-problems are solved by a limited-search-with-bound algorithm. Feasible solutions are obtained iteratively by fixing one route once. Numerical experiments show better solutions than in the literature. © 2007 Elsevier B.V. All rights reserved. Keywords: Split delivery vehicle routing; Column generation; Integer programming
1. Introduction The split delivery vehicle routing problem (SDVRP) was introduced to the vehicle routing problem (VRP) literature by [4] in the late 1980s. Unlike the capacitated vehicle routing problem (CVRP), the SDVRP allows the delivery to a demand point to be split between two or more vehicles. In many cases, allowing split deliveries results in significant savings both in the total distance traveled and the number of required vehicles compared to the corresponding CVRP solution. An example in [9] shows that the ratio of the optimal value of the SDVRP over that of the corresponding CVRP approaches 1/2 when the number of demand points goes to +∞. The savings are, in general, more significant when average customer demand is large (e.g., more than 10% of the vehicle capacity) [3]. There are a few exact algorithms in the literature to solve small-size SDVRP problems. A dynamic programming (DP) model with infinite state and action spaces is developed in [12]. The largest instance solved by [12] has nine demand points and six vehicles. A two-stage algorithm with valid inequalities is developed by [11] to successfully solve instances with up to 21 points. A cutting-plane algorithm incorporating the valid inequalities in [5] is used to solve small instances with up to 10 demand points optimally. In [2] a polyhedral study of the SDVRP is given and their cutting-plane algorithm with several valid inequalities yields good lower bounds and upper bounds ∗ Corresponding author.
E-mail address: [email protected] (M. Jin). 0167-6377/$ - see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.orl.2007.05.012
for large instances with up to 100 demand points. The gaps are around 10% for instances with 50, 75, or 100 demand points. The numerical experiments presented in Section 6 use the instances provided by [2] as benchmarks to evaluate the performance of the proposed column generation (CG) approach. A Tabu search algorithm is developed in [1] for the SDVRP. Various practical problems are modeled as SDVRPs in the literature, and they are generally solved using heuristic algorithms. Ref. [7] and [8] study the SDVRP and its variation with time windows on a grid network and propose constructive heuristics. A heuristic algorithm is used in [13] for a livestock feed distribution problem encountered on a cattle ranch in Arizona. In [14] an algorithm based on CG techniques is provided along with other heuristics to determine a flight schedule for helicopters to off-shore platform locations for exchanging work crews in the North Sea. The flight schedule does not restrict the number of total trips but has a constraint on the number of total stops in a trip. A CG approach is developed for solving the SDVRP with time windows in [9]. The CG approach proposed in this paper for the general SDVRP is similar to the algorithm, developed for one specific application, presented in [14]. However, our approach is different from that of [14] in the following ways: (1) our approach requires the solutions to use a fixed number of vehicles but does not have a restriction on the longest distance of each trip; (2) the decision variables in our master problem are defined as binary variables rather than general integer variables in order to improve the lower bound; (3) the algorithm to obtain the feasible
266
M. Jin et al. / Operations Research Letters 36 (2008) 265 – 270
solution fixes one route with the largest product of variable value and its total delivery amount and guarantees feasibility over iterations; and (4) we conduct numerical experiments on the instances published in [2] and show that the CG approach provides good results for instances with large average demand. 2. The SDVRP definition and model The SDVRP can be represented by a graph G = (V , E), where V = {0, 1, . . . , N) is the set of vertices and E = {(i, j ) : i = j } is the set of edges. Vertex 0 denotes the depot at which a fleet of identical vehicles with capacity Q are located. The other N vertices of V correspond to demand points. A positive demand di is associated with each demand point i (d0 = 0). A positive weight cij is associated with each edge (i, j ) ∈ E and represents the distance between vertices i and j. The distance from each vertex to itself is assumed to be a large positive number (cii = T , i = 0, . . . , N and T = a large number). The distances are assumed to be symmetric (cij = cj i ) and satisfy the triangular inequality (cij +cj k cik ). The SDVRP is to find routes such that the total distance traveled is minimized and the following conditions are satisfied: • The demand di of demand point i is delivered by one or more vehicles. • All vehicle routes start and end at the depot. • The total demand serviced by any vehicle does not exceed its capacity. In [3] and [12] mixed integer programming (MIP) formulations of the SDVRP are provided by fixing the total number of vehicles at K = N i=1 di /Q. This paper also requires to use K vehicles. This requirement implies that the delivery amount of each vehicle should be at least N i=1 di − (K − 1)Q. 3. The CG formulation and the column generation algorithm to obtain a lower bound (CGALB) The proposed CG approach is composed of a master problem M and a pricing sub-problem P, which are given below following the notations. The notations for the master problem M: air :
The amount delivered to demand point i by route r (in other words, a.r represents route r with delivery quantities); : The set of columns (routes); r : The minimum distance traveled to visit demand points with positive demand (air > 0) on route r; A binary variable that represents whether or not route br : r is chosen in the solution. The master problem: r br (1) M: min Z M = s.t.
r∈
air br di ,
i = 1, . . . , N,
(2)
r∈
r∈
A column r, which has N + 1 entries, is defined by a.r plus a “1”. Note that the definition of a column in the master problem M is different compared to the column definitions in CG approaches developed for routing problems without split delivery because a column in M includes delivery information. Thus, M is not the traditional set-partitioning or set-covering problem. The objective of M (1) is to minimize the total distance traveled. Constraints (2) ensure that the demand is satisfied at all points. Rather than “=” constraints in the model provided by [14], we use “ ” constraints because numerical experiments show that “ ” constraints improve the speed of the proposed CG approach. Constraint (3) guarantees that exactly K columns are used. Unlike in [14], br are defined as binary variables rather than integer variables because of the following theorem regarding the optimal solution: Theorem 1 (Dror and Trudeau [5]). Given cij is the distance between point i and point j. If the {cij } matrix satisfies the triangular inequality, then no two routes in the optimal solution to M can have more than one split demand point in common. The optimal solution should not have two identical columns that visit two or more than two points or have total delivery amount less than Q. Let be the set of all possible columns (i.e., routes with delivery quantities) that are different from each other and di /Q duplicates of column i for any demand point i for which di Q. In column i the ith entry is Q, the (N + 1)th entry is “1”, and all other entries are zero. To guarantee that are used all columns in should exactly K vehicles N satisfy N i=1 air i=1 di − (K − 1)Q. With this definition of , defining br as binary variables rather than general integer variables can improve the lower bound provided by the linear programming (LP) relaxation of M. When = , M yields an optimal solution. It is computationally intractable to enumerate all possible routes with delivery amounts. Instead, we use a CG approach to create routes using information obtained from the LP relaxation of M. Let i (i = 1, . . . , N ) be the dual price of the ith constraint for the current set of columns in the LP relaxation of M featured by . Let be the dual price of constraint (3) (the (N + 1)th constraint) in the LP relaxation of M featured by a given . The following pricing sub-problem P uses the dual prices to generate a new column: P: min
Z = P
br ∈ {0, 1}.
(3)
N
cij xij −
i=0 j =0
s.t.
N
x0j =
j =1 N j =1
i ai −
(4)
i=1
N
xj 0 = 1,
(5)
j =1
xij =
N
xj i = yi ,
i = 1, . . . , N,
(6)
j =1
ai di yi , i = 1, . . . , N, N N di − (K − 1)Q ai Q, i=1
br = K,
N N
i=1
ui − uj + (N + 1)xij N, xij , yi ∈ {0, 1}, ai , ui 0.
(7) (8)
i, j = 1, . . . , N, (9)
M. Jin et al. / Operations Research Letters 36 (2008) 265 – 270
267
bound (ICGAUB) is developed in Section 5 by fixing one column (route with delivery) iteratively. 4. LSB algorithm for the pricing sub-problem
Fig. 1. The CG algorithm for obtaining a lower bound (CGALB).
The pricing sub-problem P generates a route by minimizing the total cost (i.e., total distance traveled minus the benefit from delivery). The route starts from and ends at the depot (Constraints (5)). The model has two sets of binary variables: xij = 1 indicates that vehicle (route) r visits point j right after point i, and yi = 1 indicates that the vehicle visits demand point i (Constraints (6)). At each demand point i, the delivery amount ai should not exceed the demand (Constraints (7)). The total delivery of the route cannot exceed the vehicle capacity but should be at least N i=1 di − (K − 1)Q to make the column feasible for the SDVRP (Constraint (8)). Constraints (9) are sub-tour elimination constraints. If the solution to the pricing sub-problem P results in Z P < 0, then {ai , i = 1, . . . , N} for the first N constraints and a “1” for the (N + 1)th constraint are added to M as a new column in the next iteration. On the other hand, if the solution to P results in Z P 0, the current optimal value of the LP relaxation of M featured by is the same as the optimal value of the LP relaxation of M featured by , which is a lower bound to the SDVRP. The CG algorithm for obtaining the lower bound (CGALB) is illustrated in Fig. 1. The initial set of columns created by Step 0 of the CGALB should have a feasible solution to continue. There are exactly K columns in the initial set. The first vehicle (column) visits points 1 −1 1 i = 1, . . . , t1 , such that ti=1 di < Q and ti=1 di Q and the delivery amounts are ai1 = di for i = 1, . . . , t1 − 1 and at1 ,1 = 1 −1 di . Obviously, a solution of {xi = 1, i = 1, . . . , K} is Q − ti=1 feasible to the model M featured by the initial set of columns. In Step 2, the LP relaxation of M is solved using ILOG CPLEX 9.0 [10]. In Step 3, the pricing sub-problem P is solved, which is similar to the traveling salesman problems (TSPs) with profits in [6]. However, the TSPs with profits have a fixed profit from visiting one demand point while the benefit in the pricing sub-problem P depends on the delivery amount. Furthermore, the pricing sub-problem P has the vehicle capacity and minimal delivery constraints. The pricing sub-problem P, like the TSPs with profits, cannot be solved in a reasonable amount of time as an integer program by existing commercial solvers. The gap is still larger than 50% after 2 h when ILOG CPLEX 9.0 is used to solve a pricing sub-problem P with 50 demand points (e.g., instance S51D2 in the numerical experiments). Therefore, we develop an algorithm called limited-search-with-bounds (LSB) in Section 4 to solve the pricing sub-problem. Numerical experiments show that M featured by the set of columns that are created by the CGALB cannot be solved as an IP directly because of memory overflow and long computational time. Therefore, the iterative Column Generation approach to obtain an upper
The pricing sub-problem P determines the route and the delivery amounts for one or more vehicles. The benefit at demand point i is proportional to the delivery amount (ai ) and the dual price (i ). Thus, the delivery amount in an optimal solution satisfies the following Lemma. Lemma 2. In an optimal solution to the pricing sub-problem, each visited demand point receives its demand (ai = di ) except for one visited point with the smallest dual price that may receive partial (less-than demand) delivery. Lemma 2 is straightforward because ai appears only in Constraints (7) and (8). Without loss of generality, assume that 1 2 · · · n . Let TSP(I) denote the minimum distance traveled to visit all demand points in set I ⊆ {1, 2, . . . , N} and the depot. Let e(I ) be the largest demand point index in I. Lemma 3. If the route visits a set of demand points I plus any other demand points whose indices are larger than e(I), a lower bound to the pricing sub-problem is e(I )+J TSP(I ) − i∈I i di − i=e(I )+1 i di − e(I )+J +1 (Q − e(I )+J i∈I di − i∈I di . Here J is i=e(I )+1 i di ) − , when Q > e(I )+J a nonegative integer such that i∈I di + i=e(I )+1 i di Q e(I )+J +1 and i∈I di + i=e(I )+1 i di > Q. All terms with an index larger than N will not be included in this summation. Lemma 3 is based on three facts: (1) visiting more demand points increases the total distance traveled because of the triangular inequality assumption; (2) the dual prices i are in a descending order; and (3) the route has at most one demand point that is partially served. The following LSB algorithm is developed based on Lemmas 2, 3, and the capacity constraint (Fig. 2). Note that the average demand at each demand point is typically larger in an SDVRP compared to a CVRP. Therefore, the TSP in Step 1 of the LSB algorithm can be solved relatively easily by using CPLEX 9.0. Set I is fathomed and does not need any further exploration when (1) there is no more capacity; (2) the set includes the last demand point that has the smallest dual price; or (3) L(I ) U , where L(I ) is a lower bound for the solution including the demand points in I and other demand points with indices larger than e(I ) (i.e., smaller or the same dual price). If set I is not fathomed then it is explored by adding the next demand point e(I )+1 to I in Step 4. If set I is fathomed, then backtracking is performed in Step 5 until e(I ) = N and I = ∅. If e(I ) = N , then e(I ) is replaced with e(I ) + 1. When I =∅, the algorithm stops. With the LSB algorithm, all columns created by P have at most one demand point whose delivery is split. In other words, the upper bound obtained by the proposed CG is bounded from below by the master problem M with = {all columns with at most one split delivery point}.
268
M. Jin et al. / Operations Research Letters 36 (2008) 265 – 270
Fig. 2. The LSB algorithm.
5. The iterative column generation approach to obtain an upper bound (ICGAUB) When no more columns are created, the solution to the linear relaxation of M featured by the created column set yields a lower bound for the SDVRP. Usually the solution has fractional numbers so it is not feasible to the SDVRP. Solving the original M with br defined as binary variables can provide a feasible solution. However, this is very time consuming. Based on numerical experiments, when there are more than 50 demand points, using optimization solvers to solve the original M featured by created by CGALB is not practical. In fact, after 24 h, ILOG CPLEX 9.0 cannot even find a feasible solution though a feasible solution definitely exists because of the initial set of columns created by Step 0 in the CGALB. We propose the following ICGAUB in order to find a feasible solution by fixing one route with delivery amounts at each iteration. In the first iteration, after solving the LP relaxation of M for the original SDVRP, the algorithm fixes the variable in the optimal solution with the largest product of xr N i=1 air to 1. The rule of selecting which column (route with delivery amounts) should be fixed is based on the following two conjectures: 1. The columns with higher variable values in the solution to the LP relaxation of M are “good” routes. 2. Less demand left for the reduced SDVRP obtained after fixing the column can improve the solution to the reduced SDVRP. We conducted numerical experiments using other rules, but the above given rule was better than others. Fixing the column
means that the points whose demand are fully satisfied (i.e., air = br ) should be removed. All columns that visit any of these fully satisfied points should also be removed from the column set. The demand of the point that is partially served (i.e., 0 < air < br ) is reduced by the served amount. These operations (Step 3 in Fig. 3) create a new and smaller SDVRP with |{i : air∗ = di }| fewer demand points and 1 less vehicle. The new SDVRP is solved by using the CG algorithm in Fig. 1 again. Since the minimal delivery requirement for each route of the new SDVRP may be smaller than that of the SDVRP in the previous iteration, some of the old columns may violate the new minimal delivery requirement and should be removed. Note that Step 0 in the CGALB is still necessary to make the first LP for the reduced SDVRP feasible. Because all columns at each iteration meet the minimal delivery requirement for the reduced SDVRP, the reduced SDVRP at iteration l always has a feasible solution by using K −l +1 vehicles. When the reduced SDVRP is small enough (e.g., the product of the number of remaining points and the number of required vehicles is less than a threshold value Th = 80), the exact solution approach introduced in [11] is used to obtain an optimal solution. An upper bound for the original SDVRP is obtained after adding back the costs of all fixed columns. 6. Numerical experiments A desktop computer with a Pentium 4 processor, 2.8 GHz CPU and 2 G memory is used to conduct numerical experiments. The results are compared to the numerical results of [2]. In [2] 25 randomly generated instances are solved. However,
M. Jin et al. / Operations Research Letters 36 (2008) 265 – 270
269
Fig. 3. The iterative column generation approach to obtain an upper bound (ICGAUB).
Table 1 Comparison of the CG Algorithm with the cutting-plane algorithm Instance name
S51D2 S51D3 S51D4 S51D5 S51D6 S76D2 S76D3 S76D4 S101D2 S101D3 S101D5 eil76B
K
9 15 27 23 41 15 23 37 20 31 48 14
Cutting plane
Column generation
UB
LB
Gap %
UB
LB
Gap %
T1
T2
726 972 1677 1440 2327 1147 1474 2257 1393 1975 2915 1023
676.63 905.22 1520.67 1272.86 2113.03 1020.32 1346.29 2011.64 1270.97 1739.66 2630.43 937.47
6.80 6.87 9.32 11.61 9.20 11.04 8.66 10.87 8.76 11.92 9.76 8.36
723.37 968.85 1657.61 1439.92 2300.21 1185.72 1504.94 2219.07 1474.51 2012.86 2954.96 1063.75
694.98 922.72 1505.35 1297.46 2108.59 1066.17 1397.43 2019.91 1349.77 1837.33 2725.5 981.14
3.92 4.76 9.19 9.89 8.33 10.08 7.14 8.97 8.46 8.72 7.77 7.77
3555 323 14 2 <1 7511 549 188 16 507 2698 10 30 222
5978 607 260 46 243 12 806 2030 1813 47 658 7959 847 45 084
The underlined numbers show the results where the CG approach was worse compared to the cutting-plane approach.
we only include 12 of those 25 instances in our analysis. The included instances, listed in Table 1, are the ones in which average demand at a point is about 15% larger than the vehicle capacity. This is done because the CVRP rather than the SDVRP is a more appropriate model for instances with small average demand. We note that the proposed CG approach is not efficient for solving instances with small average demand due to the resulting large
tree in the LSB algorithm and the large computational time for solving the embedded TSPs. The full set of 25 instances is available at www.uv.es/∼belengue/SDVRP/sdvrplib.html. In order to obtain a better solution and speed up the algorithm, rather than adding only the best column, other “good” columns out of the LSB algorithm are also added to the master problem M. If the objective function value corresponding to set
270
M. Jin et al. / Operations Research Letters 36 (2008) 265 – 270
I is negative and lower than a threshold value, then the column corresponding to set I is added to M. In Table 1, column K shows the minimum number of vehicles required. The results of the cutting-plane method are provided by [2]. In the results of the CG approach, T 1 shows the CPU time (in seconds) spent by the CGALB to obtain the lower bound and T 2 shows the total time (in seconds) spent by the ICGAUB to obtain the upper bound. The CG algorithm outperforms the cutting-plane method in all instances with respect to the error gap. The CG algorithm finds better upper bounds in six instances and better lower bounds in ten instances out of the twelve instances. In general, the CG approach takes much less time when the average demand at each point is larger. Similar phenomenon is observed in the CG approach developed for the SDVRP with time windows in [9]. This result is interesting because most of other algorithms (e.g., [2,12]) developed for the SDVRP require more time as average demand increases for problems with the same number of demand points, which is perhaps caused by the increase in the number of binary variables. In our CG algorithm, larger average demand means more efficient implementation of the LSB algorithm because of the capacity constraint. Furthermore, larger demand instances have fewer column candidates that have at most one split delivery (i.e., fewer iterations of the CG algorithm). 7. Conclusion To improve the lower bounds, cuts can be added to the problem. For example, the number of vehicles visiting a subset of demand points should be enough for serving the total demand of the points in the subset. Though it is easy to generate the cuts, more studies are necessary to continue the CG approach
by modifying the algorithm to solve the pricing sub-problem with the cuts. References [1] C. Archetti, A. Hertz, M.G. Speranza, A Tabu search algorithm for the split delivery vehicle routing problem, Transp. Sci. 40 (2006) 64–73. [2] J.M. Belenguer, M.C. Martinez, E. Mota, A lower bound for the split delivery vehicle routing problem, Oper. Res. 48 (2000) 801–810. [3] M. Dror, G. Laporte, P. Trudeau, Vehicle routing with split deliveries, Discrete Appl. Math. 50 (1994) 239–254. [4] M. Dror, P. Trudeau, Savings by split delivery routing, Transp. Sci. 23 (1989) 141–145. [5] M. Dror, P. Trudeau, Split delivery routing, Naval Res. Logist. 37 (1990) 383–402. [6] D. Feillet, P. Dejax, M. Gendreau, Traveling salesman problems with profits, Transp. Sci. 39 (2005) 188–205. [7] P.W. Frizzell, J.W. Giffin, The bounded split delivery vehicle routing problem with grid networks distances, Asia Pacific J. Oper. Res. 9 (1992) 101–116. [8] P.W. Frizzell, J.W. Giffin, The split delivery vehicle scheduling problem with time windows and grid network distance, Comput. Oper. Res. 22 (1995) 655–667. [9] M. Gendreau, D. Feillet, P. Dejax, C. Gueguen, Vehicle routing with time windows and split deliveries, Technical Paper 2006-851, Laboratoire d’Informatique d’Avignon, 2006. [10] ILOG, ILOG CPLEX 9.0 User’s Manual, 2003. [11] M. Jin, K. Liu, R.O. Bowden, A two-stage algorithm with valid inequalities for the split delivery vehicle routing problem, Int. J. Production Econ. 105 (2007) 228–242. [12] C.G. Lee, M. Epelman, C.C. White, A shortest path approach to the multiple-vehicle routing with split picks-up, Transp. Res.—Part B 40 (2006) 265–284. [13] A. Mullaseril, M. Dror, J. Leung, Split-delivery routing heuristics in livestock feed distribution, J. Oper. Res. Soc. 48 (1997) 107–116. [14] G. Sieksma, G.A. Tijssen, Routing helicopters for crew exchanges on off-shore locations, Ann. Oper. Res. 76 (1998) 261–286.