Robust_branch-and-cut-and-price_for_the_capacitated_vehicle_routing_problem

Robust Branch-and-Cut-and-Price for the Capacitated Vehicle Routing Problem Ricardo Fukasawa, Marcus Poggi de Araga˜o, Marcelo Reis, Eduardo Uchoa Relato´rios de Pesquisa em Engenharia de Produc¸a˜o RPEP Vol.3 no.8 (2003) Nitero´i, 10 de setembro de 2003 Robust Branch-and-Cut-and-Price for the Capacitated Vehicle Routing Problem Ricardo Fukasawa 1, Marcus Poggi de Araga˜o 2, Marcelo Reis 2, Eduardo Uchoa 3 1 Gapso Inc. R. Jardim Botaˆnico 674, sala 614, Rio de Janeiro, Brasil, 22461-000 rfukasawa@gapso.com.br 2 Departamento de Informa´tica, PUC-Rio R. Marqueˆs de Sa˜o Vicente, 225, Rio de Janeiro, Brasil, 22453-900 {mreis,poggi}@inf.puc-rio.br 3 Departamento de Engenharia de Produc¸a˜o, Universidade Federal Fluminense R. Passo da Pa´tria, 156, Bloco E, sala 440, Nitero´i, Brasil, 24210-240 uchoa@producao.uff.br (corresponding author) September 10th, 2003 Abstract During the eigthies and early nineties, the best exact algorithms for the Capacitated Vehicle Routing Problem (CVRP) utilized lower bounds obtained by Lagrangean relaxation or column generation. Next, the advances in the polyhedral description of the CVRP yielded branch-and- cut algorithms giving better results. However, several instances in the range of 50–80 vertices, some proposed more than 30 years ago, can not be solved with current known techniques. This paper presents an algorithm utilizing a lower bound obtained by minimizing over the intersec- tion of the polytopes associated to a traditional Lagrangean relaxation over q-routes and the one defined by bounds, degree and the capacity constraints. This is equivalent to a linear program with an exponential number of both variables and constraints. Computational experiments show the new lower bound to be superior to the previous ones, specially when the number of vehicles is large. The resulting branch-and-cut-and-price could solve to optimality almost all instances from the literature up to 100 vertices, nearly doubling the size of the instances that can be consistently solved. Further progress in this algorithm may be soon obtained by also using other known families of inequalities. 1 Introduction The Capacitated Vehicle Routing Problem (CVRP) consists of given an undirected graph G = (V,E) with vertices numbered as {0, 1, . . . , n} (vertex 0 represents the depot and the remaining vertices represent clients), client demands d(1), . . . , d(n), lengths `(e) associated to edges in E and K vehicles with capacity C; determine routes for each vehicle satisfying the following constraints: (i) each route starts and ends at the depot, (ii) each client is visited by a single vehicle, and (iii) the total demand of clients visited in a route is at most C. The objective is to minimize the sum of the routes length. This classical NP-hard problem was first proposed by Dantzig and Ramsey in 1959 [19] and have received a lot of attention from the optimization community since then, due to its widespread applications, and also because it is a natural generalization of the Travelling Salesperson Problem (TSP). A landmark exact algorithm for the CVRP was presented in 1981 by Christofides, Mingozzi and Toth [16] using a Lagrangean bound obtained by solving minimum q-route problems. A q-route starts at the depot and traverses a sequence of clients limiting the demand accumulated to at most C and returns to the depot. It is not necessarily a simple path, for some clients may be visited more than once. Therefore, the set of valid CVRP routes is contained in the set of q-routes. The resulting branch-and-bound could solve instances up to 25 vertices, a quite respectful size at that time. Several other branch-and-bound algorithms using Lagrangean bounds appear in the litera- ture. This same article [16] also describes a lower bound based on k-degree center trees, minimum spanning trees having degree K ≤ k ≤ 2K on the depot, plus 2K − k least cost edges. Other authors propose Lagrangean bounds based on K-trees, which are sets of n+K edges spanning G, like Fisher [21] and Martinhon, Lucena and Maculan [29]. There is also an algorithm based on minimum b-matchings having degree 2K at the depot and 2 on the remaining vertices by Miller [30]. The Lagrangean bounds can be improved by dualizing capacity inequalities [21, 30] and also combs and multistar inequalities [29]. Another kind of exact algorithms have its stem on the formulation of the CVRP as a set partition problem by Balinsky and Quandt [9]. In that formulation, a column covers a set of vertices S with total demand not exceeding C and have the cost of a minimum route over {0} ∪ S. Bramel and Simchi-Levi [13] proved that for certain natural classes of instances, the ratio between the lower bounds given by that formulation and the optimal solution values asymptotically approaches 1 as the number of clients grows. However, that formulation in itself is not practical because pricing over the exponential number of columns require the solution of capacitated prize-collecting TSPs, a problem almost as difficult as the CVRP itself. Agarwal, Marthur and Salkin [3] proposed a column generation algorithm on a modified set partition where column costs are given by a linear function over the vertices yielding a lower bound on the actual route cost. Columns with the modified cost can be priced by solving easy knapsack problems. Hadjconstantinou et al. [23] derive lower bounds from heuristic solutions to the dual of the set partitioning formulation. Those dual solutions are obtained by the so-called additive approach, combining the q-path with the K-shortest path relaxations. For further information and some comparative results on the above mentioned algorithms, we refer the reader to the survey by Toth and Vigo [37]. Starting in the 90’s, most of the research effort on the CVRP is now concentrated on the polyhedral description of convex hull of the edge incidence vectors that correspond to K fea- sible routes and on the development of effective separation algorithms for the families of valid inequalities identified (for instance, [4, 14, 18, 6, 7, 2, 31, 27]). In particular, Araque et al. [5], Augerat et al. [8], Blasum and Hochstattler [12], Ralphs et al. [36], Achutan, Caccetta and Hill [1] and Lysgard, Letchford and Eglese [28] describe complete branch-and-cut algorithms, some including sophisticated inequalities such as framed capacity, strengthened combs, multi- star, among others. (The taxonomy and nomenclature of valid inequalities for the CVRP is not uniform in the literature, we are following [28] in this article.) Although those are the best exact algorithms currently available for the CVRP, the lower bounds obtained at the root nodes, even after adding all those cuts, are not very tight for many instances with as little as 40 vertices. The quality of those bounds is specially problematic for larger values of K, say K ≥ 7. Many nodes in a branch-and-cut tree may have to be explored in order to close the resulting duality gaps. Even resorting to massive computational power (up to 80 processors running in parallel in a recent work by Ralphs [36, 35]) several instances with less than 80 vertices, including some proposed more than 30 years ago by Christofides and Eilon [15], can not be solved at all. In fact, it seems that branch-and-cut algorithms for the CVRP are experimenting a “diminishing returns” phenomenon, where substantial theoretical and implementation efforts leads to practical results that are only marginally better than those of previous works. This work presents a new exact algorithm for the CVRP that seems to break through this situation. The main idea is to combine the branch-and-cut approach with the old q-routes ap- proach (which we interpret as a column generation instead of the original Lagrangean relaxation) 3 to derive superior lower bounds. Since the resulting formulation has an exponential number of both columns and rows, this leads to a branch-and-cut-and-price algorithm. Computational experiments over the main instances from the literature show that this algorithm can consis- tently solve instances with up to 100 vertices. Seventeen open instances were solved for the first time. Those results were obtained using only bound, degree and capacity inequalities. An improved algorithm is soon expected by also separating framed capacity, strengthened combs, multistar, partial multistar and extended hypotour inequalities, but this was not done yet due to the complexity of implementing the corresponding separation heuristics. 1 The idea of combining column and cut generation in order to achieve improved lower bounds have existed since the eighties. The difficulty, pointed out for rarely performing that, was the fact that the new dual variables corresponding to separated cuts could change the structure of the pricing subproblem, leading to an intractable pricing (Barnhart et al.[11] and Wilhelm[39] comment on this matter). Recently, several researchers [38, 24, 25, 20, 10] have noted that cuts expressed in terms of variables from a suitable original formulation can be incorporated to the column generation without disturbing the pricing. We use the term “robust branch-and- cut-and-price” to refer to such algorithms: branch-and-bounds over linear programs with an exponential number of both columns and rows and where neither branching nor separation ever change the structure of the pricing subproblems. The article [33] is a detailed discussion on that matter. In particular, it proposes some reformulation techniques that extend the applicability of robust branch-and-cut-and-price algo- rithms to virtually any combinatorial optimization problem. Moreover, it is argued that such algorithms may lead to advances on a wide variety of classic problems, from TSP to graph coloring. The present article on the CVRP is part of a larger effort to support that claim. Very good results were also already obtained on the capacitated minimum spanning tree problem [22] and on the generalized assignment problem [32]. 2 The New Formulation A classical formulation for the CVRP [26] represents by xij the number of times a vehicle traverses edge (i, j) ∈ E. Let V+ be the set {1, . . . , n} of client vertices. Given a set S ⊆ V+, let d(S) be the sum of the demands of all vertices in S and δ(S) denote the cut-set defined by S. Let also k(S) = dd(S)/Ce. Consider the polytope in R|E| next defined: P1 =  ∑ e∈δ({i}) xe = 2 ∀ i ∈ V+ (1)∑ e∈δ({0}) xe = 2.K (2)∑ e∈δ(S) xe ≥ 2.k(S) ∀S ⊆ V+ (3) xe ≤ 1 ∀ e ∈ E \ δ({0}) (4) xe ≥ 0 ∀ e ∈ E Constraints (1) state that each non-depot vertex is visited once by a vehicle and constraints (2) that K vehicles must go out and in the depot. Constraints (3) are the capacity inequalities requiring that all subsets are served by enough vehicles. Constraints (4) ensure that each edge non adjacent to the depot can be traversed no more than once. Edges adjacent to the depot can be used twice, this is the case where a route serves only one client. The integer vectors x in P1 define all feasible solutions for the CVRP. Due to the exponential number of inequalities (3), the lower bound given by L1 = min{ ∑ e∈E `e.xe : x ∈ P1} 1J. Lysgard kindly promised to let us experiment with his implementation of those heuristics, described in [27, 28], in a near future. 4 has to be calculated by a cutting plane algorithm. In fact, as separating capacity inequalities is NP-hard, one usually resorts to separation heuristics. In that case, the actual bounds obtained may be a little worse than L1. Modern branch-and-cut algorithms for the CVRP, like [5, 8, 12, 36, 1, 28], may improve this bound by also separating several other known families of inequalities. Formulations with an exponential number of columns can be obtained by defining variables (columns) that correspond to valid CVRP routes, as first proposed by Balinski and Quandt [9]. Those formulations are not practical since pricing such columns amounts to solving a capacitated prize-collecting TSP, a strongly NP-hard problem. One way of dealing with this difficulty is by enlarging the set of the columns to correspond to q-routes. As already mentioned in the introduction, a q-route is a closed path visiting a sequence {0, v1, . . . , vr, 0} of vertices such that ∑r i=1 d(vi) ≤ C and vi 6= 0 for each i in {1, . . . , r}, r ≥ 1. Note that vertices can appear more than once in the q-route. The cost of a q-route is given by `(0, v1) + ∑r−1 i=1 `(vi, vi+1) + `(vr, 0). Finding a minimum cost q-route can be done in pseudo- polynomial time, O(n2.C). Christofides, Mingozzi and Toth [16] has shown that restricting the set of q-routes to those without 2-cycles (subpaths (i, j), (j, i)) does not change this complexity. Let Q be a m × p matrix where the columns are the edge incidence vectors of all p possible q-routes with no 2-cycles (except for those corresponding to single client routes). Denote by qej the coefficient associated to edge e in the jth column of Q and consider the following polytope in Rp+|E|: P2 =  p∑ j=1 qej .λj − xe = 0 ∀e ∈ E (5) p∑ j=1 λj = K (6)∑ e∈δ({i}) xe = 2 ∀ i ∈ V+ (1) xe ≥ 0 ∀ e ∈ E λj ≥ 0 ∀ j = 1, . . . , p Constraints (5) define the coupling between variables x and λ. Constraint (6) states the number of vehicles to utilized. Since the definition of the columns already imposes the capacity con- traints, it remains to add the degree constraints (1). The integer vectors in P2 also define all feasible solutions for the CVRP. This is because setting to 1 a variable λj corresponding to a q-route that is not a valid CVRP route, violates some degree constraints (1). Due to the exponential number of variables λ, the lower bound given by L2 = min{ ∑ e∈E `e.xe : x ∈ Projx(P2)} has to be calculated by a column generation algorithm, or equivalently, by Lagrangean re- laxation. A branch-and-bound using such relaxation [16] was the first really successful exact algorithm for the CVRP. The alternative way of describing the polyhedra associated to a column generation or to a Lagrangean relaxation in terms of two sets of variables, λ and x, used in the definition of P2, is called Explicit Master in [33]. The main contribution of this article is to propose a formulation that amounts to optimizing over the intersection of polytopes P1 and Projx(P2), which follows. The Explicit Master format makes easy to see that such formulation must be the following. 5 P1 ⋂ Projx(P2) = Projx  ∑ e∈δ({i}) xe = 2 ∀ i ∈ V+ (1)∑ e∈δ({0}) xe = 2.K ∀S ⊆ V+ (2)∑ e∈δ(S) xe ≥ 2.k(S) ∀S ⊆ V+ (3) xe ≤ 1 ∀ e ∈ E \ δ({0}) (4) p∑ j=1 qej .λj − xe = 0 ∀e ∈ E (5) p∑ j=1 λj = K (6) xe ≥ 0 ∀ e ∈ E λj ≥ 0 ∀ j = 1, . . . , p Remark that constraints (6) can be discarded. They are implied by constraints (2) and (5). The improved lower bound is, then, given by: L3 = min{ ∑ e∈E `e.xe : x ∈ P1 ⋂ Projx(P2)}. This bound can be calculated by solving a linear program with an exponential number of both variables and constraints. A more compact equivalent linear program can be obtained by elim- inating the x variables substituting the xe values given by (5) in constraints (1)-(4). We will refer to the resulting LP as the Dantzig-Wolfe Master problem (DWM). DWM =  L3 = min p∑ j=1 ∑ e∈E `e.q e j .λj (7) s.t. p∑ j=1 ∑ e∈δ(i) qej .λj = 2 ∀ i ∈ V+ (8) p∑ j=1 ∑ e∈δ({0}) qej .λj = 2.K (9) p∑ j=1 ∑ e∈δ(S) qej .λj ≥ 2.k(S) ∀S ⊆ V+ (10) p∑ j=1 qej .λj ≤ 1 ∀ e ∈ E \ δ({0}) (11) λj ≥ 0 ∀ j = 1, . . . , p Remark that bound constraints (11) are really necessary, because imposing upper bounds on the x variables can not be done only by upper bounds on λ variables. 3 The Branch-and-Cut-and-Price 3.1 Initialization We first try to tighten the upper bounds on the variables associated to edges incident to the depot. IfK−1 vehicles are not enough to serve the n−1 clients in the set V +\{i} (a bin-packing problem) then the upper bound of edge (0, i) can be reduced to 1. We do not actually solve the bin-packing, the simple lower bound dd(V + \ {i})/Ce is used instead. The DWM corresponding to the root node is initialized with a set of artificial variables of high cost covering constraints (8) and (9). Constraints (10) and (11) are omitted. 6 3.2 Column Generation The pricing subproblem amounts to determining whether there is a q-route corresponding to a column with a negative reduced cost. In particular, we want to find the q-route with the most negative reduced cost and, if there exists, several others negative reduced cost routes as different as possible from each other. We remind that only 2-cycle free q-routes are being considered. Let µ, ν, pi and ω be the dual variables associated to constraints (8), (9), (10) and (11), respectively. The reduced cost associated to edge e is given by c¯e =  `e − µi − µj − ∑ S|δ(S)3e piS − ωe e = {i, j} ∈ E \ δ({0}) `e − ν − µj − ∑ S|δ(S)3e piS e = {0, j} ∈ δ({0}) The reduced cost of a column (λ variable) is given by the sum of the reduced costs of the edges in the q-route associated to it. A minimum reduced cost q-route can be obtained by dynamic programming in O(n2.C) time. We resort the reader to [16, 17]. We generate columns to DWM after the root initialization, after some cut from (10) or (11) is added to DWM or, in the non-root nodes, after a branching. 3.3 Cut Generation Given a fractional solution λ¯ to DWM, the corresponding fractional solution x¯ is calculated by using equalities (5). All separation routines work over x¯. Bound constraints (11) are easily separated by inspection. Remark that having all those constraints in the current LP would add |E| rows and too many non-zeros to it. Capacity constraints (10) are heuristically separated using a set of six routines written by T. Ralphs as part of a demo implementation of a branch-and-cut for the CVRP on the Symphony framework [34]. Those routines are pretty fast, but they are not as good on finding violated capacity constraints as state-of-art separating routines, like those described on [8, 7, 12, 36, 28]. For that reason, we also use an exact separation routine using the MIP solver embedded in CPLEX. Let x¯ be a fractional solution. Define yi as a binary variable equal to 1 iff vertex i belongs to set S and wij be a variable that is equal to 1 if edge (i, j) belongs to δ(S). For each value of M in {1, . . . , dd(V +)/Ce − 1} we solve: z = min ∑ (i,j)∈E x¯ijwij s.t. wij ≥ yj − yi ∀ (i, j) ∈ E wij ≥ yi − yj ∀ (i, j) ∈ E∑ i∈V + d(i)yi ≥ M.C + 1 y0 = 0 yi ∈ {0, 1} ∀ i ∈ V + wij ≥ 0 ∀ (i, j) ∈ E If z ≤ 2(M + 1), then the capacity inequality over set S is violated. Solving those MIPs is computationally costly. In level 0 of the branch-and-bound tree (the root node), we separate capacity constraints exactly, calling the MIP solver whenever the heuris- tic separation fails, as many times as necessary to be sure that no capacity constraint is violated. In levels 1 up to 4 we perform only one round of exact separation. In deeper levels, only the heuristic separation is performed. 3.4 Branching We branched by calculating x¯ as above and selecting the variable with value closest to 0.5. Note that if for some edge (0, j), 1 < x¯0j < 2, then for some other edge (i, j), 0 < x¯ij < 1. In other words, we never need to branch over a fractional variable with value greater than 1. 7 Fixing a variable xe to 0 or 1 can be enforced in the DWM by a constraint similar to (11). Edges fixed to 0 can also be removed from the pricing subproblem. 4 Computational Experiments We present our computation results on 44 instances from the literature available at http://www.branchandcut.org/VRP , a site maintained by T. Ralphs. Those instances are: (i) the 12 instances from series A with no less than 50 vertices; (ii) the 13 instances from series B with no less than 50 vertice