An Exact Algorithm Based on Cut-and-Column Generation for the Capacitated Location-Routing Problem

This article was downloaded by: [128.122.253.212] On: 15 October 2014, At: 04:19 Publisher: Institute for Operations Research and the Management Sciences (INFORMS) INFORMS is located in Maryland, USA INFORMS Journal on Computing Publication details, including instructions for authors and subscription information: http://pubsonline.informs.org An Exact Algorithm Based on Cut-and-Column Generation for the Capacitated Location-Routing Problem Claudio Contardo, Jean-François Cordeau, Bernard Gendron To cite this article: Claudio Contardo, Jean-François Cordeau, Bernard Gendron (2014) An Exact Algorithm Based on Cut-and-Column Generation for the Capacitated Location-Routing Problem. INFORMS Journal on Computing 26(1):88-102. http://dx.doi.org/10.1287/ ijoc.2013.0549 Full terms and conditions of use: http://pubsonline.informs.org/page/terms-and-conditions This article may be used only for the purposes of research, teaching, and/or private study. Commercial use or systematic downloading (by robots or other automatic processes) is prohibited without explicit Publisher approval, unless otherwise noted. For more information, contact [email protected]. The Publisher does not warrant or guarantee the article’s accuracy, completeness, merchantability, fitness for a particular purpose, or non-infringement. Descriptions of, or references to, products or publications, or inclusion of an advertisement in this article, neither constitutes nor implies a guarantee, endorsement, or support of claims made of that product, publication, or service. Copyright © 2014, INFORMS Please scroll down for article—it is on subsequent pages INFORMS is the largest professional society in the world for professionals in the fields of operations research, management science, and analytics. For more information on INFORMS, its publications, membership, or meetings visit http://www.informs.org http://pubsonline.informs.org http://dx.doi.org/10.1287/ijoc.2013.0549 http://dx.doi.org/10.1287/ijoc.2013.0549 http://pubsonline.informs.org/page/terms-and-conditions http://www.informs.org INFORMS Journal on Computing Vol. 26, No. 1, Winter 2014, pp. 88–102 ISSN 1091-9856 (print) — ISSN 1526-5528 (online) http://dx.doi.org/10.1287/ijoc.2013.0549 © 2014 INFORMS An Exact Algorithm Based on Cut-and-Column Generation for the Capacitated Location-Routing Problem Claudio Contardo Département de management et technologie, ESG UQÀM, Montréal (Québec), Canada H2X 3X2, [email protected] Jean-François Cordeau Canada Research Chair in Logistics and Transportation, HEC Montréal and CIRRELT, Montréal (Québec), Canada H3T 2A7, [email protected] Bernard Gendron Département d’informatique et de recherche opérationnelle, Université de Montréal and CIRRELT, Montréal (Québec) Canada H3C 3J7, [email protected] In this paper we present an exact algorithm for the capacitated location-routing problem (CLRP) based oncut-and-column generation. The CLRP is formulated as a set-partitioning problem that also inherits all of the known valid inequalities for the flow formulations of the CLRP. We introduce five new families of inequalities that are shown to dominate some of the cuts from the two-index formulation. The problem is solved by column generation, where the subproblem consists in finding a shortest path of minimum reduced cost under capacity constraints. We first use the two-index formulation for enumerating all of the possible subsets of depot locations that could lead to an optimal solution of cost less than or equal to a given upper bound. For each of these subsets, the corresponding multiple depot vehicle routing problem is then solved by means of column generation. The results show that we can improve the bounds found in the literature, solve to optimality some previously open instances, and improve the upper bounds on some other instances. Key words : location routing; vehicle routing; branch-and-cut-and-price; column generation History : Accepted by Karen Aardal, Area Editor for Design and Analysis of Algorithms; received July 2011; revised April 2012; accepted January 2013. Published online in Articles in Advance June 14, 2013. 1. Introduction In the capacitated location-routing problem (CLRP) we are given a set I of potential facilities and a set J of customers that represent source and destina- tion points, respectively, of a unique commodity. With every facility i ∈ I are associated a fixed operating cost fi and a capacity bi representing the maximum amount of commodity that can be delivered to the customers using the vehicles assigned to that facility. With every customer j ∈ J is associated a demand dj for the commodity. The problem can be defined on an undirected graph G = 4V 1E5, where V = I ∪ J is the vertex set and E is the edge set. With every edge e = 8i1 j9 we associate a routing cost cij , and these costs are assumed to satisfy the triangle inequality. The fleet of vehicles is assumed to be homogeneous, each vehicle having a capacity Q representing the maxi- mum amount of commodity that can be delivered to the customers. The problem is to choose a subset of facilities to operate and to construct (possibly several) vehicle routes around these facilities so that every customer is visited exactly once, the amount of com- modity delivered by each vehicle does not exceed its capacity, the amount of commodity delivered from each facility respects its capacity, and the total cost (the sum of the fixed costs and the routing costs) is minimized. The CLRP arises in several real-world applica- tions. Labbé and Laporte (1986) solve the problem of locating postal boxes while minimizing a linear combination of routing costs (those of the mail col- lecting trucks) and customer inconvenience produced by their distance to the nearest postal box. Billionet et al. (2005) consider a location problem arising in mobile networks. The problem consists in locating radio-communication stations, designing rings, and building antennae inside these rings at minimum cost. Gunnarsson et al. (2006) solve a location-routing problem arising in the pulp distribution industry in Scandinavia. The CLRP can be formulated as a three-index mixed-integer program (Perl and Daskin 1985). 88 D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS 89 In such a formulation, asymmetries in the distance matrix and heterogeneities in the vehicle capacities can easily be taken into account. However, because of the large number of variables and its poor lin- ear programming relaxation it has no practical use within an enumeration method such as branch and bound. In the context of exact algorithms for solv- ing the CLRP, Belenguer et al. (2011) developed a two-index formulation and proposed several fam- ilies of valid inequalities, such as y-capacity cuts (y-CC), path constraints (PC), facility degree con- straints (FDC), imparity constraints (IC), and facility capacity inequalities (FCI). They solve the problem by means of branch and cut and their algorithm suc- ceeds in solving small- and medium-size instances with up to 50 customers. Contardo et al. (2011) intro- duced three new formulations of the CLRP based on vehicle flows and commodity flows. They intro- duced strengthenings of the FCI as well as location- routing comb inequalities (LR-CI), location-routing generalized large multistar inequalities (LRGLM), and y-generalized large multistar inequalities (y-GLM), exploiting the fact that facilities have limited capac- ities. Their algorithms were able to solve instances containing up to 100 customers, the largest for branch- and-cut methods (Contardo et al. 2011). Akca et al. (2009) developed a set-partitioning formulation based on a Dantzig-Wolfe decomposition of the three-index model. They solve the problem by means of branch and price, where the subproblem is a shortest path problem under capacity constraints (SPPRC). Their formulation provides reasonably good bounds at the root node of the search tree but does not appear to be effective for closing the gap using branch- ing. Baldacci et al. (2011) also formulate the CLRP as a set-partitioning problem. They use three differ- ent relaxations of the formulation that are applied sequentially in an additive manner. In the last step, they solve a small number of multiple depot vehicle- routing problems (MDVRP) by means of a cut-and- price-and-branch method, in which the root node is solved by column generation, and then enumerate all of the remaining columns whose reduced cost is smaller than a given gap. The resulting integer pro- gram is then solved by means of a general-purpose integer programming solver. Baldacci et al. (2011) use a strengthened version of the capacity cuts (CC) as well as clique inequalities. The bounds provided by their model are very tight and able to solve instances with up to 199 customers and 14 facilities. The CLRP is NP-hard as it generalizes both the capacitated VRP (CVRP) and the capacitated facil- ity location problem (CFLP). Moreover, the presence of capacities for both the vehicles and the facilities makes it particularly hard. Because of this, solu- tion approaches for solving medium- and large-size instances have mainly focused on the development of heuristics. These heuristics, in most cases, use some decomposition scheme to divide the prob- lem into a design subproblem for the location deci- sions and an operational subproblem for the routing part (Perl and Daskin 1985, Hansen et al. 1994, Lin et al. 2002, Wu et al. 2002, Liu and Lin 2005). More recently, Prins et al. (2006a, b; 2007) proposed sev- eral metaheuristics that include memetic algorithms, cooperative Lagrangean relaxation with tabu search, and greedy randomized adaptive search procedure (GRASP). Computational experience shows that the second approach is the most effective one for tackling large instances of the CLRP. The contributions of this paper can be summarized as follows. (i) We adapt the set-partitioning formulation of Akca et al. (2009) so that all of the cuts valid for the two-index formulation of the CLRP (Belenguer et al. 2011, Contardo et al. 2011) can be easily incorporated. (ii) We introduce three procedures that are applied sequentially and that allow, in most cases, reducing the CLRP to a series of multiple depot VRP, as in Baldacci et al. (2011). Our computational results show that our procedures can be stronger than those used by Baldacci et al. (2011) for some instances. (iii) We introduce several new families of cuts that are effective for closing the optimality gap. More- over, our computational experience shows that using state-space relaxation in the pricing problem suffices to get bounds close to those obtained by pricing on elementary routes (routes that do not contain cycles). (iv) We introduce a new fathoming rule that accel- erates the solution of the pricing subproblems. As a result, our algorithm is able to solve all instances that are also solved by the exact method of Baldacci et al. (2011) as well as four previously open instances. Additionally, we improve the best known feasible solution for three other instances. Moreover, for the instances that remain unsolved we improve the best known lower bounds. The paper is organized as follows: In §2 we present some formulations of the CLRP, namely, the two-index vehicle-flow formulation of Belenguer et al. (2011) as well as the set-partitioning formulation of Akca et al. (2009). In §3 we describe the valid inequalities used through this paper. It includes some known valid inequalities from the two-index formulation and the set-partitioning problem as well as new valid inequal- ities that are shown to be valid for the set-partitioning formulation of the CLRP. In §4 we describe the exact algorithm used to solve the CLRP to optimality. We first give an overview of the complete algorithm, which is done in three stages that are applied sequen- tially. We also describe the separation algorithms used to find violated valid inequalities. We then describe D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP 90 INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS the different procedures used at each of the three stages as well as the pricing algorithms used to solve the corresponding set-partitioning problems. Finally, we discuss some computational issues that are mostly implementation specific and that have an important impact on the performance of the algorithm. In §5 we present our computational results and compare against the state-of-the-art solvers for solving the CLRP. We conclude in §6 with a summary of the pro- posed methodology and discuss possible avenues of future research. 2. CLRP Formulations In this section we first present the two-index vehicle- flow formulation of the CLRP of Belenguer et al. (2011) and the set-partitioning formulation introduced by Akca et al. (2009). We also show that any inequal- ity valid for the two-index formulation can easily be extended to the set-partitioning formulation. 2.1. Two-Index Vehicle-Flow Formulation Belenguer et al. (2011) proposed the following two- index vehicle-flow formulation for the CLRP: For every vertex set U , let „4U5 be the edge subset con- taining all edges with exactly one endpoint in U . For two disjoint vertex sets T , U , let 4T 2 U5 be the edge subset containing all edges with one endpoint in T and the other in U . For every facility i ∈ I , let zi be a binary variable equal to 1 iff facility i is selected for opening. For every edge e ∈ E, let xe be a binary variable equal to 1 iff edge e is traversed once by some vehicle. Finally, for every edge e ∈ „4I5 let ye be a binary variable equal to 1 iff edge e is used twice by some vehicle. For a given edge set F ⊆ E let x4F 5= ∑ e∈F xe, y4F 5 = ∑ e∈F ye. For a given customer subset S ⊆ J , let d4S5 = ∑ j∈S dj , r4S5 = ‘d4S5/Q’ (which is a lower bound on the number of vehicles needed to serve the customers in S), and S̄ = J\S. The formula- tion is the following: min { ∑ i∈I fizi + ∑ e∈E cexe + 2 ∑ e∈„4I5 ceye } (TIF) subject to x4„4j55+ 2y4I2 8j95= 2 j ∈ J 1 (1) x4„4S55+ 2y4I2 S5≥ 2r4S5 S ⊆ J 1 —S— ≥ 21 (2) xij + yij ≤ zi i ∈ I1 j ∈ J 1 (3) x4I2 8j95+ y4I2 8j95≤ 1 j ∈ J 1 (4) x ( 4I\8i95∪ S̄2 S ) + 2y4I\8i92 S5≥ 2 i ∈ I1 S ⊆ J 1 d4S5 > bi1 (5) x4„4S55≥ 2 ( x48h92 I ′5+ x48j92 I\I ′5 ) S ⊆ J 1 —S— ≥ 21 h1 j ∈ S1 I ′ ⊂ I1 (6) zi ∈ 80119 i ∈ I1 (7) xe ∈ 80119 e ∈ E1 (8) ye ∈ 80119 e ∈ „4I50 (9) Constraints (1) are the degree constraints at cus- tomer nodes. Constraints (2) are capacity cuts (CC) whose role is to forbid at the same time proper tours disconnected from facilities and tours serving a demand larger than Q. Constraints (3) ensure that there is no outgoing flow leaving from closed facili- ties. Constraints (4) are the path constraints for single customers. They forbid routes of the form i1 → j → i2, i11 i2 ∈ I , i1 6= i2, j ∈ J . Constraints (5) are the facil- ity capacity inequalities. They forbid the existence of routes leaving from a same facility i and serving a demand larger than bi. Constraints (6) are the path constraints for multiple customers. Their role is to prevent the route of a single vehicle from joining two different facilities. 2.2. Set-Partitioning Formulation We now describe the set-partitioning formulation introduced by Akca et al. (2009) and used by Baldacci et al. (2011), and link it to the two-index vehicle-flow formulation so that all of the known cuts for the CLRP are also valid. Let us denote by ìi the set of all routes (possibly containing cycles) starting and ending at facility i ∈ I and servicing a subset of at least two cus- tomers with total demand of Q or less, and let ì = ⋃ i∈I ìi be the set of all possible routes servicing two or more customers with total accumulated demand of Q or less. For every l ∈ ì let us associate a binary variable ‹l equal to 1 if l appears in the optimal solu- tion of the CLRP and 0 otherwise, and a cost cl for using this route. For every edge e ∈ E and route l ∈ì let qel be the number of times that edge e appears in route l. If ì is restricted to contain only elementary routes then qel is a binary constant, otherwise it can be a general integer. On the other hand, let us define binary variables yij for 8i1 j9 ∈ „4I5 equal to 1 iff cus- tomer j is served from facility i by a single-customer route. For a subset of facilities H ⊆ I and a subset of customers S ⊆ J we let y4H2 S5 = ∑ i∈H ∑ j∈S yij . Note that given that distances satisfy the triangle inequal- ity, the optimal solution of this problem will only con- tain elementary paths even if ì is enlarged to contain routes with cycles. In fact, in this case it is always possible to build from a solution with cycles another solution with elementary routes at lower cost. Let us extend the demands to facility nodes by letting dv = 0 for every v ∈ I . A valid formulation for the CLRP is min { ∑ i∈I fizi + ∑ l∈ì cl‹l + 2 ∑ e∈„4I5 ceye } (SPF) D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS 91 subject to ∑ l∈ì ∑ e∈„48j95 qel ‹l + 2y4I2 8j95= 2 j ∈ J 1 (10) ∑ l∈ìi ∑ e=8h1j9∈E 4dh+dj5q e l ‹l+2 ∑ e∈„48i95 djyij ≤2bizi i∈ I , (11) ‹l ∈ 80119 l ∈ì1 (12) ye ∈ 80119 e ∈ „4I51 (13) zi ∈ 80119 i ∈ I 0 (14) In this formulation, constraints (10) ensure that each customer is served exactly once. Constraints (11) are the facility capacity inequalities. They ensure that the demand served from any facility i will not exceed its capacity bi. The distinction between single- and multiple-customer routes naturally defines a relation- ship between vehicle-flow variables x from the two- index formulation and ‹, as follows: ∑ l∈ì qel ‹l − xe = 0 e ∈ E0 (15) In this way, all of the valid inequalities from the two-index formulation of the CLRP can be trans- lated into the set-partitioning formulation by using identities (15). 3. Valid Inequalities In this section we describe the valid inequalities that can be applied to formulation (SPF) and that strengthen the LP relaxation. First, we describe some of the valid inequalities that have been developed in the context of the two-index and three-index for- mulations by Belenguer et al. (2011) and Contardo et al. (2011). We then describe new families of valid inequalities that are shown to dominate several of the former and that effectively strengthen formula- tion (SPF). 3.1. Valid Inequalities for the Two-Index Formulation The valid inequalities for formulation (TIF) include several different families. After a series of preliminary tests, we have decided to keep only a subset of them, namely, strengthened comb inequalities (SCI), framed capacity inequalities (FrCI), effective strengthened facility capacity inequalities (ESFCI), facility degree constraints (FDC), and location-routing comb inequal- ities (LRCI). To include any of these constraints into formulation (SPF) we use identity (15). For details on the inequalities we refer to Lysgaard et al. (2004), Belenguer et al. (2011), and Contardo et al. (2011). 3.2. Valid Inequalities for the Set-Partitioning Formulation The valid inequalities for the set-partitioning formu- lation include a strengthening of the SCC introduced by Baldacci et al. (2008) for solving the CVRP and also used by Baldacci et al. (2011) for the CLRP. These constraints can also be seen as a strengthening of the y-CC introduced by Belenguer et al. (2011), which are valid for formulation (TIF). We also introduce strengthenings of the degree constraints (10), of the SFCI, ESFCI, and FrCI. We complement this with the addition of subset-row inequalities (SRI). Any constraint in the two-index space can be trans- lated into a constraint in the route space by using identity (15). However, the constraints translated this way will not take into account the fact that a route can cross more than once a given subset of vertices. Let us start by defining some notation. For a given subset of customers S ⊆ J and a route l ∈ì, let ŽlS be a binary constant equal to 1 iff route l visits at least one customer in S. We also denote, for a given route vec- tor ‹ and a route subset R ⊆ ì, Ž‹4R1 S5= ∑ l∈R Ž l S‹l. Finally, for a given subset of facilities H ⊆ I we define ìH = ⋃ i∈H ìi. We now describe some new families of valid inequalities that strengthen the set-partitioning formulation. 3.2.1. y-Strengthened CC (y-SCC). Let us con- sider, for a given customer set S ⊆ J and subset S ′ ⊂ S such that r4S\S ′5 = r4S5, the corresponding y-CC as described in Belenguer et al. (2011) and Contardo et al. (2011) in terms of the route variables after using identities (15): ∑ l∈ì ∑ e∈„4S5 qel ‹l + 2y4I2 S\S ′5≥ 2r4S50 (16) In these constraints, the quantity ∑ e∈„4S5 q e l is greater than or equal to 2 for routes visiting set S. Baldacci et al. (2008) noted that the CC (2) can be strengthened by setting the coefficient of a given path variable ‹l to be 0 if l does not serve a customer in S and 2 oth- erwise, rather than counting the number of edges of l that are also in „4S5. Let us call these constraints SCC, and for details we refer to Baldacci et al. (2011). For the y-CC we can apply the same reasoning, as stated in the following proposition. Proposition 1. Let S ⊆ J be a subset of customers, and S ′ ⊂ S such that r4S\S ′5 = r4S5. The following constraint is valid for the CLRP and dominates the y-CC used in Belenguer et al. (2011) and Contardo et al. (2011) and the SCC used in Baldacci et al. (2011): Ž‹4ì1S5+ y4I2 S\S ′5≥ r4S50 (17) We call this constraint the y-strengthened CC (y-SCC). Proof. See the supplemental material (available as supplemental material at http://dx.doi.org/10.1287/ ijoc.2013.0549). 3.2.2. Strengthened Degree Constraints (SDEG). Degree constraints (10) count the number of times that a certain node is traversed. This is normally equal to 1 if a route visits that node or 0 otherwise. When D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP 92 INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS the route set ì is enlarged to contain nonelemen- tary routes, a node can be visited several times in a route. In that case, a stronger version of the degree constraint is Ž‹4ì1 8j95+ y4I2 8j95≥ 1 j ∈ J 0 (18) In our algorithm we have found that the addition of these constraints when pricing on nonelementary routes is an effective method to get bounds close to the ones obtained by pricing on elementary routes. Indeed, the problem of finding an appropriate balance between speed and lower bound quality for differ- ent variants of the SPPRC has already been studied and is a key aspect in the performance of column generation-based algorithms for vehicle routing prob- lems (see, e.g., Boland et al. 2006, Desaulniers et al. 2008, Righini and Salani 2008). This intuition is sup- ported by the following proposition and corollary. Proposition 2. Assume that ì is enlarged to contain routes with cycles. After the addition of a strengthened degree constraint for a customer j ∈ J , the solution of the linear relaxation of problem (SPF) is such that ‹l = 0 for every l ∈ì visiting j more than once. Proof. See the online supplement. Corollary 1. After the addition of all strengthened degree constraints, the value of the linear relaxation of problem (SPF) is the same as if ì was restricted to contain only elementary routes. 3.2.3. Set-Partitioning SFCI (SP-SFCI). The strengthened facility capacity inequalities (SFCI) were introduced by Belenguer et al. (2011) and Contardo et al. (2011) and are a generalization of constraints (5). Let us define for a customer set S ⊆ J and a facility subset H ⊂ I the quantity r4S1H5 = max801 ‘4d4S5 − b4H55/Q’}, with b4H5 = ∑ i∈H bi. Here, r4S1H5 is a lower bound on the number of vehicles that are needed to serve customers in S from facilities other than those in H . Let S ⊆ J , S ′ ⊂ S, and H ⊆ I be such that r4S1H5 = r4S\S ′2 H5. The following inequality, called set-partitioning SFCI (SP-SFCI), is valid for the CLRP and dominates the SFCI: Ž‹4ì\ìH1 S5+ y4I\H2 S\S ′5≥ r4S1H50 (19) Proposition 3. Constraints (19) are valid for the CLRP and dominate the SFCI. Proof. See the online supplement. 3.2.4. Set-Partitioning ESFCI (SP-ESFCI). The effective SFCI were introduced by Belenguer et al. (2011) and Contardo et al. (2011) and are valid for the two-index formulation. They can be seen as a strengthening of the SFCI by noticing that the right- hand side of such a constraint can in fact be lifted whenever zi = 0 for some i ∈ H . Indeed, let S ⊆ J , S ′ ⊂ S, and H ⊆ I be such that r4S\S ′1H5 = r4S1H5 and r4S\S ′1H\8i95 = r4S1H\8i95, and let i ∈ H be any facility in H . The following constraint, named set- partitioning ESFCI (SP-ESFCI), is valid for the CLRP: Ž‹4ì\ìH1 S5+ y4I\H2 S\S ′5 ≥ r4S1H\8i95+ zi ( r4S1H5− r4S1H\8i95 ) 0 (20) The validity proof follows from the validity of the SP-SFCI for the two cases zi = 1 and zi = 0. Indeed, in the former case when facility i is open constraint (20) corresponds to the SP-SFCI. In the latter case, when facility i is closed it can be ignored from the SP-SFCI. Because r4S1H\8i95 ≥ r4S1H5 the case zi = 0 corre- sponds to a lifting of the original SP-SFCI. 3.2.5. Strengthened Framed Capacity Inequali- ties (SFrCI). The framed capacity inequalities were developed by Augerat (1995) for the CVRP and later succesfully used by other authors in the develop- ment of algorithms based on cutting planes and col- umn generation (Lysgaard et al. 2004, Fukasawa et al. 2006). Given a customer set S, that we call the frame, and a partition of it 4Si5ti=1, the related FrCI seen in formulation (TIF) is x4„4S55+ 2y4I2 S5+ t ∑ i=1 ( x4„4Si55+ 2y4I2 Si5 ) ≥ 2 ( BPP4S — 4Si5 t i=15+ t ∑ i=1 r4Si5 ) 1 (21) where BPP4S — 4Si5ti=15 represents the solution of the following bin-packing problem. For every i = 11 0 0 0 1 t consider ‘d4Si5/Q’ items of size Q except for the last item that will have size d4Si5− 4‘d4Si5/Q’−15Q. Also, set the bins to have size Q. In addition to using iden- tity (15) to adapt this constraint to formulation (SPF), the same observation as done for the y-SCC, SDEG, SP-SFCI, and SP-ESFCI can be applied. The following constraint, called strengthened FrCI (SFrCI), is valid for the CLRP and also dominates the FrCI: Ž‹4ì1S5+ t ∑ i=1 Ž‹4ì1Si5+ 2y4I2 S5 ≥ BPP4S — 4Si5 t i=15+ t ∑ i=1 r4Si50 (22) Before proving the validity of the SFrCI we need the following lemma. Lemma 1 (Augerat 1995). Let S ⊆ J and 4Si5 t i=1 be a partition of S. If ‘d4S1 ∪ S25/Q’ = ‘d4S15/Q’ + ‘d4S25/Q’ then BPP4S — S11 S21 0 0 0 1 St5 ≥ BPP4S — S1 ∪ S21 S31 0 0 0 1 St5. Otherwise BPP4S — S11 S21 0 0 0 1 St5 + 1 ≥ BPP4S — S1 ∪ S21 S31 0 0 0 1 St5. D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS 93 Proof. See Augerat (1995). Proposition 4. Constraints (22) are valid for (SPF). Proof. See the online supplement. 3.2.6. Subset-Row Inequalities. The subset-row inequalities (SRI) are a special case of the clique inequalities (Balas and Padberg 1976) and are valid for the set-partitioning formulation of the CLRP. They were proposed by Jepsen et al. (2008) for the vehicle routing problem with time windows (VRPTW) with the objective of capturing the strength of the clique inequalities without compromising the performance of the pricing algorithm. Let us consider the conflict graph H‹ constructed as follows: The vertices of H‹ are the routes l ∈ ì such that ‹l > 0. Two vertices in V 4H‹5 are linked by an edge if they share at least one customer. A clique in H‹ is a maximal complete induced subgraph of H‹. For every clique C⊆H‹, the following clique inequality is valid for the CLRP: ∑ v∈V 4C5 ‹v ≤ 10 (23) The addition of clique inequalities into the master problem SPF has, however, an important drawback: they make the pricing problem of finding routes (with or without cycles) of negative reduced cost much more difficult. Indeed, during the pricing problem it must be checked whether a partial path participates in a clique. This is equivalent to checking if a par- tial path intersects every column already in a clique in at least one customer node, which in practice is difficult to do. A subset-row inequality C3 of three customers is a special case of a clique inequality to which we associate a subset of customers •4C35⊆ J of cardinality 3, and such that every column in C3 inter- sects •4C35 in at least two customers. Now, the pricing problem can be accelerated as only three comparisons are needed to check if a given route participates in a subset-row inequality. The results obtained by Jepsen et al. (2008) show that the subset-row inequalities of three customers provide the best balance between the quality of the lower bound and the computational effort needed to achieve it. 4. Solution Methodology In this section we describe the exact algorithm that solves the CLRP to optimality. We first provide an overview of the proposed algorithm, and we comple- ment it with a high-level pseudocode in which we describe the interactions between the different stages of our procedure, which will be detailed later in the text. Then, we describe the separation algorithms used to find violated inequalities. Then, we describe each of the three different stages of our algorithm in detail, which are applied sequentially. The first stage is based on the solution of formulation (TIF) by branch and cut, and the two following stages are based on the solu- tion of problem (SPF) by column-and-cut generation. Finally, we describe the computational issues in the implementation of the proposed algorithm. 4.1. Overview of the Proposed Algorithm The proposed algorithm is performed in three stages that are applied sequentially. In the first stage, we solve formulation (TIF) by branch and cut. For this, we have used the code of Contardo et al. (2011) to include many of the inequalities used in that arti- cle. The main difference with respect to the original method introduced in Contardo et al. (2011) is the branching rule that now forces only integrality on the location variables z. The objective is to enumer- ate all possible configurations of open facilities that may lead to an improvement of the current upper bound. For this, the best known solutions reported by Baldacci et al. (2011) are used as cutoff values dur- ing the enumeration process. In the second stage, for each of the possible configurations of open facilities, we apply column-and-cut generation on formulation (SPF) subject to the facilities that are open, and solve the linear relaxation. This procedure provides a lower bound on the solution of the corresponding multiple depot CVRP. If that lower bound is strictly smaller than the best known solution, in a third stage we apply column generation to enumerate all columns whose reduced costs are less than or equal to the remaining gap. We then impose integrality on the route variables and solve the resulting integer pro- gram with a general-purpose solver such as CPLEX. The following pseudocode describes more precisely the interactions between the three stages. Algorithm 1 (CLRP) 1. Let zBKS be the value of the current best known solution. 2. Solve problem (TIF) by partial branch and bound, let I⊆P4I5 be the set of configurations H whose associated lower bound v14H5 is strictly smaller than zBKS. 3. Sort I by nondecreasing value of v14H5 and let gap14H5← zBKS − v14H5 for all H . 4. for all H ∈I such that v14H5 < zBKS do. 5. Solve the corresponding MDVRP by column- and-cut generation and make three attempts to generate all columns whose reduced costs are less than or equal to the current gap. 6. Let v24H5 be the final lower bound associated and let gap24H5← zBKS − v24H5 be the associated gap. 7. if gap24H5 < 0 then 8. Discard H 9. end if D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP 94 INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS 10. end for 11. for all H ∈I such that v24H5 < zBKS do 12. If columns were not generated in the previous stage for H , enumerate all columns whose gap is less than or equal to gap24H5 and include them in the associated problem (SPF). 13. Impose integrality over all variables and solve the resulting IP with a general-purpose solver. 14. Update zBKS if the value of the new solution is better than the incumbent. 15. end for. 4.2. Separation Algorithms We now describe the separation algorithms used to separate the different families of valid inequalities used in our algorithm. Our separation strategy is as follows: We first try to generate cuts translated from the two-index formulation. If no such cuts can be found, we try to generate cuts SDEG, y-SCC, SP-SFCI, SP-ESFCI, and SFrCI. If it fails, we try to generate cuts SRI. This strategy allows us to keep the number of strong constraints small as their inclusion in the pric- ing algorithm makes it harder. 4.2.1. Inequalities Translated from TIF. For the valid inequalities translated from the two-index for- mulation using identity (15), such as y-CC, SFCI, ESFCI, SCI, LR-CI, or FrCI, we use the separation algo- rithms introduced by Lysgaard et al. (2004), Belenguer et al. (2011), and Contardo et al. (2011). 4.2.2. SDEG, y-SCC, SP-SFCI, SP-ESFCI, and SFrCI. Although there is a polynomial number of SDEG constraints, we do not add them all at the beginning of the algorithm, rather we check if for a certain weak degree constraint, its related strong con- straint is violated, and add it to the problem. For the remaining constraints, we use the same principle. In fact, we check if, for any previously found weak con- straint y-CC, SFCI, ESFCI, or FrCI, its related strong constraint is violated and in this case we add it to formulation SPF. 4.2.3. Subset-Row Inequalities. The separation of the subset-row inequalities is done by enumeration just as in Jepsen et al. (2008). We check for every triplet 4i1 j1 k5 ∈ J 3, i < j < k if the corresponding SRI is violated. If it is the case, it is added to the master problem. 4.3. Enumeration of Candidate Subsets by Branch and Cut In this procedure, an enumeration method based on a branch-and-cut algorithm (Contardo et al. 2011) is applied to problem (TIF) after dropping the integral- ity constraints on the edge variables x and y. This procedure is used to obtain candidate subsets H ⊆ I of facilities such that the problem restricted to open- ing these facilities could lead to a feasible solution with cost less than or equal to a given upper bound. We denote the set that contains the subsets H by I. For finding the subsets in I, a good upper bound is needed to prune nodes in the branching tree. In our method, we have used the best feasible solutions found in the literature. For large instances, however, the computation of the whole branching tree can be prohibitive. In this case, the branch-and-bound algo- rithm is terminated earlier and the uninspected nodes are also added to I. Now, the facilities in a given subset H ∈ I are not only those that are open but also those that could not be fixed in the current node. Throughout the process, different families of valid inequalities are added to strengthen the formulation. However, we only add cuts in nodes whose depth is less than or equal to 5. For each candidate set H ⊆ I generated by the algorithm we proceed as follows. (i) Based on reduced costs, perform variable fix- ing on the location variables z, in case set H contains facilities that remained unfixed. (ii) Based on reduced costs, perform variable fixing on the edge variables x. The edge variables x fixed to zero are removed from the shortest path problems used in the column generation in §§4.4 and 4.5. (iii) Compute the optimal dual variables associated to the degree constraints (1). These values are used to initialize the stabilization centers during the column generation algorithm, as explained in §4.6.2. (iv) Compute Km4H5 as an upper bound on the max- imum number of routes that serve two or more custom- ers, namely, Km4H5 = max8 1 2x4„4H552 4x1y1 z5 ∈ A9, where A stands for the set of constraints (1)–(9) plus the generated cuts and after dropping the integral- ity conditions. This quantity is used to compute a lower bound of the solution of the current multiple depot vehicle routing problem associated to facilities H at each iteration of the column generation algo- rithm during the second stage, as explained in §4.4.4. Note that although the method developed by Baldacci et al. (2011) also decomposes the problem into three stages that are applied sequentially, in the first stage their procedure computes a global lower bound by solving a relaxation of the set-partitioning problem. This bound is then used to discard non- promising subsets H ⊆ I . In §5 we present computa- tional results comparing the first stage of our method with the one suggested by Baldacci et al. (2011). 4.4. Column-and-Cut Generation In this procedure, the following state-space relaxation of formulation (SPF) is solved by means of column- and-cut generation for a fixed candidate subset of facilities H ∈ I. Instead of considering elementary D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS 95 routes (i.e., routes without cycles), we allow routes that contain cycles of length three or more, i.e., for nodes i 6= j 6= k 6= i the subpaths i → i, i → j → i are forbidden, but the sequence i → j → k → i is permitted. The pricing problem consists in finding routes without cycles of length one or two and such that the reduced costs are minimized. This prob- lem is known in the literature as the 2-cyc-SPPRC (Desrochers et al. 1992). This is an important differ- ence with respect to the method of Baldacci et al. (2011) in which the resolution of the subproblem is restricted to elementary routes. During the computa- tion, we add the cuts described in §3. The violation threshold for the strong cuts is initially set to 003. When no more columns of negative reduced cost or violated cuts can be detected, the current objective function value is in fact a valid lower bound for the problem. Let us call this lower bound z∗. We run algo- rithm ENUM-ESPPRC (described in the online sup- plement) to price out the remaining columns l ∈ ì such that c̄l ≤ zUB − z∗. We have set two hard limits to algorithm ENUM-ESPPRC: the number of labels cannot exceed at any time a maximum ”max = 106, and the total number of generated columns cannot exceed ãmax = 107. In case of success of this proce- dure, the columns generated are stored in a column pool P and the violation threshold for strong con- straints is lowered to 0001. Otherwise, we lower the violation threshold (thus generating more cuts) and continue with the process. This is done at most three times before finishing the column generation process. For instance, for the case of constraints SDEG, the sequence of violation thresholds is 400310025100210015. Whenever the column enumeration ENUM-ESPPRC is done with success, at every following iteration of the column generation method, we do not solve the pricing problem 2-cyc-SPPRC but rather check the reduced costs of columns in P. Note that the size of set P can be huge and computing the reduced cost of every column in it can be very cumbersome. To deal with this issue, at every iteration after the cre- ation of P in which no columns of negative reduced cost were found, we also delete from the pool all the columns l such that c̄l > zUB − z∗. At the very end of the column-and-cut generation, we either prune the current node if the final lower bound is greater than or equal to zUB, or otherwise solve the integer problem with the columns generated so far, with the hope of improving the upper bound. In what follows, we first describe the decomposition of the reduced costs for the constraints translated from formulation (TIF), namely, all of the constraints in (SPF) plus the cuts that are valid for this formulation. We then show how to incorporate the set-partitioning constraints, such as y-SCC, SDEG, SP-SFCI, SP-ESFCI, SFrCI, and SRI into the computation of the reduced costs. We then describe the pricing problem 2-cyc-SPPRC that suits our problem with the additional cuts. We end by describing how we compute lower bounds out of the result of the pricing problem. 4.4.1. Decomposition of the Reduced Costs Edge by Edge. Let us first suppose that only constraints (10) and (11) (with duals  and ‚, respectively) have been added to the problem. For every i ∈ H , define the reduced cost of an edge e ∈ E4J 5∪ „48i95 as c̄e =      ce−4h+j5 −4dh+dj5‚i if e=8h1j9∈E4J 5, ce−j −dj‚i if e=8i1j9∈„48i950 (24) Let us write a route l ∈ ìi as a sequence of edges in E, that is l = 4et5 p t=1 (in the case where cycles are permitted, edges may appear more than once in the sequence). Thus, the reduced cost of such a route is given by the following expression: c̄l = p ∑ t=1 c̄et 0 (25) It follows that in this case a column of minimum reduced cost can be computed as the solution of —H — shortest path problems with resource constraints. Moreover, the addition of any cut of the general form ∑ i∈H ’izi + ∑ e∈E ∑ l∈ì qel”e‹l + ∑ e∈„4H5 œeye ≤ (26) produces a contribution to the computation of the reduced cost of the columns that can still be decom- posed by edge, thus without breaking the shortest path structure of the pricing. This is the case for all of the cuts valid for the two-index formulation of the CLRP after being translated to formulation (SPF) using identity (15). 4.4.2. Addition of the Strong Constraints and Effect on the Reduced Costs. When a constraint can- not be decomposed edge by edge, as for constraints (17)–(20), (22), or (23), the contribution to the reduced cost cannot be decomposed edge by edge, and thus the original structure of the SPPRC is broken. Indeed, consider a SRI for a clique C3 ∈ C3 such that for •4C35= 8i1 j1 k9 with dual variable ‘ ≤ 0. The reduced cost c̄l of a route l ∈ ì that crosses at least two of those three customers must be augmented by −‘ units. For the other strong constraints SDEG, y-SCC, SP- SFCI, SP-ESFCI, or SFrCI, the contribution to the reduced cost is related to the simple intersection of path l with the sets describing the constraints. For instance, if we consider a SDEG constraint associated to a customer j and with dual variable ‘ ≥ 0, then the D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP 96 INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS reduced cost of a route l will be reduced by ‘ units if l passes through node j . Now, consider a constraint y-SCC for given sets S ⊆ J and S ′ ⊂ S as in (17) with dual value ‘ ≥ 0. The contribution to the reduced cost will reduce it by ‘ units if l intersects set S. For a SP-SFCI or SP-ESFCI associated to sets S ⊂ J , S ′ ⊂ S, H with dual variable ‘ ≥ 0, the contribution to the reduced cost will reduce it by ‘ units if route l crosses set S but is not linked to a facility in H . Finally, for the SFrCI associated to set S and partition Si, i = 11 0 0 0 1 t and with dual variable ‘ ≥ 0, the reduced cost must be reduced by ‘ units once for each time that route l intersects either S or any of its subsets. 4.4.3. The Pricing Problem. The pricing problem corresponds to solving —H — 2-cyc-SPPRC, one for each facility in H . For each i ∈ H , the objective is to find routes (eventually containing cycles) with negative reduced cost starting and ending at i such that the total demand served is less than or equal to Q. The resources associated to each label during the recur- sion are (1) vehicle load; (2) binary resources related to constraints SDEG, y-SCC, SP-SFCI, SP-ESFCI, and SFrCI; and (3) resources for taking into account the SRI. The algorithm used to solve these problems is based on dynamic programming (DP), as was done by several authors (Desrochers et al. 1992, Feillet et al. 2007, Baldacci et al. 2008, Jepsen et al. 2008, Righini and Salani 2008). Moreover, it is also possible to solve it by means of bidirectional DP (BDP). In classical uni- directional DP, paths are extended until reaching the depot node while ensuring that loads do not exceed capacity. In BDP, however, paths are extended until reaching half of the capacity for later joining paths pairwise. In this section we describe the 2-cyc-SPPRC algorithm used in the context of the CLRP. For general use of the dynamic programming method for solving the SPPRC we refer to the papers cited earlier. 4.4.3.1. Resources Description. As said before, three different types of resources are considered in the problem: vehicle load resource; resources associated to constraints SDEG, y-SCC, SP-SFCI, SP-ESFCI, and SFrCI; and resources associated to SRI. Vehicle load. The vehicle load is defined by an inte- ger variable q that keeps track of the load of the cur- rent path. It is updated every time a path is extended to a customer node. Resources associated to SDEG, y-SCC, SP-SFCI, SP- ESFCI, and SFrCI. For each of the constraints SDEG, y-SCC, SP-SFCI, and SP-ESFCI, the associated resource is defined by a single Boolean variable that takes the value true if the path intersects the proper set, as described earlier. We designate those sets as critical sets, and denote them by S4C5 for every constraint C. For each constraint SFrCI, there will be not one, but as many Boolean variables as the size of the partition, plus one for the frame. Each of these variables will take the value true if the path crosses the proper set. Now, we do not have one but several critical sets that we denote by S4C1k5. Any time that one of these Boolean variables passes from false to true, the reduced cost of the current path is reduced according to the value of the dual variable. Resources associated to SRI. For every clique C3 with •4C35 = 8i1 j1 k9 we associate three binary variables rC34k51 k = 11213 that are initialized to 0 until the path crosses one of the customers, in which case the proper variable is set to 1, and the reduced cost of a path will be updated whenever rC3415 + rC3425 + rC3435 reaches the value 2. 4.4.3.2. The 2-cyc-SPPRC Algorithm. In this section we describe the bidirectional dynamic programming algorithm used to find routes of negative reduced cost. This algorithm is run independently for each facility i ∈ H , and so we assume that i is fixed throughout this section. We start by giving the defi- nition of a label in the dynamic programming algo- rithm. Intuitively, a label corresponds to a partial path starting at i and visiting a subset of customers, some of which are possibly visited more than once. We also describe the dominance and fathoming rules used to discard labels. Next, we describe the path-joining pro- cedure to construct feasible paths from a given pair of labels. At the end, we describe the skeleton of the algorithm by means of a pseudocode. Label definition. A label L is defined by the following: (i) A node v4L5 that is the end node of the partial path represented by label L. (ii) The reduced cost c̄4L5 of the partial path represented by label L. (iii) The accumulated load q4L5 of the partial path represented by label L. (iv) Resources resC4L5 associated to the binding constraints SDEG, y-SCC, SP-SFCI, SP-ESFCI, SFrCI, and SRI. For constraints SFrCI and SRI we write resC4L1k5 for the different subresources associated to these constraints. (v) A pointer to the predecessor label pred4L5 of L. (vi) An integer variable vdom4L5 initially set to −1 and updated whenever L is found to be dominated by a label L′, in which case we set vdom4L5= v4pred4L′55. (vii) A Boolean variable proc4L5 initialized to false and updated to true whenever the algorithm processes the label and inspects its neighbors. (viii) A list succ4L5 of pointers to the successors of L. Dominance rule. Let L, L′ be two labels. We de- note SRILL′ = 8C3 ∈ SRI: ∑ k resC34L1k5 ≤ 1 and 6 ∑ k resC34L ′1 k5≥ 2 or ∃k s.t. resC34L1k5 < resC34L ′1 k579, nC1L1L′ = —8k: resC4L1k5 < resC4L′1 k59— and OTHL1L′ = D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS 97 8C ∈ SDEG ∪ y-SCC ∪ SP-SFCI ∪ SP-ESFCI: resC4L5 < resC4L′59. We will say that L is dominated by L′ if (i) v4L5= v4L′5; (ii) q4L5≥ q4L′5; (iii) c̄4L5≥ c̄4L′5− ∑ C ∈SRILL′ ‘C + ∑ C ∈SFrCI nC1L1L′‘C + ∑ C ∈OTHLL′ ‘C . The dominance rule is a direct application of the one used by Archetti et al. (2011) for the inclusion of SRI and k-path inequalities in the context of the VRP with split deliveries and time windows (VRPSDTW). A label L that is dominated by another label L′ cannot be directly eliminated unless v4pred4L55 = v4pred4L′55 or if vdom4L5 y 8−11v4pred4L′559. In that case, label L is removed and recursively we also remove all of its successors in succ4L5. Otherwise, vdom4L5 is set to v4pred4L′55. Note that the inclusion of SDEG constraints allows the dominance rule to weaken with respect to a traditional elementarity con- straint, in which the condition for dominance would be resC4L5≥ resC4L′5 for each C ∈ SDEG. Fathoming rule. In addition to the dominance crite- rion, a fathoming rule can be applied if a lower bound on the cost of extending a path can be computed. For- mally, let L be a label and let LB4L5 be a lower bound on the reduced cost that can be obtained by extend- ing L, computed as follows. First of all, discard SRI as their dual variables are negative. For every binding strong constraint C ∈SC= SDEG∪y-SCC∪SP-SFCI∪ SP-ESFCI ∪ SFrCI, with dual variables 4‘C5C∈SC, and for every edge e crossing the critical sets related to these constraints, we decrease the reduced cost of that edge by ‘C/2 units. We refer to this procedure as underestimation of constraint C. Because a route that crosses a customer set S must have at least two edges in „4S5 then the reduced cost of a path computed in this way will in fact be a lower bound on the real reduced cost. We then solve the related 2-cyc-SPPRC with no resources associated to strong constraints, and compute functions f , g, and  as follows: f 4p1 j5= min { c̄4L52 v4L5= j1 q4L5≤Q− p+ dj } 1 (27) 4p1 j5= v ( pred4arg min8f 4p1 j595 ) 1 (28) g4p1 j5= min { c̄4L52 v4L5= j1 q4L5≤Q− p+ dj1 v4pred4L55 6=4p1 j5 } 0 (29) The meaning of functions f , , and g is the fol- lowing: The quantity f 4p1 j5 represents the mini- mum reduced cost among all partial paths visiting a demand less than or equal to Q− p+ dj and having j as end node. This is used as the completion bound for a path represented by a label L with load q4L5= p, end node v4L5 = j , and whose penultimate node is differ- ent from 4p1 j5. The quantity 4p1 j5 represents the penultimate node of the partial path reaching f 4p1 j5. Finally, function g represents the minimum reduced cost among the partial paths serving a demand of Q − p + dj , having j as end node and a node differ- ent from 4p1 j5 as penultimate node. This quantity is used as the completion bound for paths represented by labels L such that q4L5 = p, v4L5 = j and such that the penultimate node is 4p1 j5. Therefore, we define h4L5 as the completion bound for a label L, as follows: h4L5= { f 4q4L51v4L55 if 4q4L51v4L55 6= v4pred4L55, g4q4L51v4L55 otherwise0 Let us define nC1 j = —8k2 j ∈ S4C1k59— for a SFrCI and a customer j ∈ J . A lower bound on the reduced cost reachable by extending label L can be computed as LB4L5 = c̄4L5+h4L5+ 1 2 ∑ C ∈SC\SFrCI v4L5∈S4C5 ‘C + 1 2 ∑ C ∈SFrCI nC1v4L5‘C 0 (30) The two sums aim to compensate the fact that the contribution of the underestimated constraints C ∈SC is being considered at least 1.5 times in c̄4L5 and h4L5 whenever v4L5 ∈ S4C5 or nC1v4L5 > 0, thus tightening LB4L5. If a label L is such that LB4L5 > 0, then L can be discarded. Similar fathoming rules have been implemented by Christofides et al. (1981) and Baldacci et al. (2008, 2010, 2011), for instance. Note that we have used unidirectional DP for computing functions f , g, . From an implementation point of view, it only differs from the BDP in the fact that now all labels are inspected for extension and not only those whose load is less than or equal to Q/2, so at the end the joining of paths is not necessary. Note also that this fathoming procedure can be generalized (and also strengthened) by keeping as resources, thus without underestimating, the k constraints C ∈ SC with the largest duals, where k is a parameter defined a priori. After doing a series of experiments, we let k = min8201 —SC—/59. For these constraints, the coeffi- cients in the sums in (30) can now be lifted to 1, as the contribution to the reduced cost of customer v4L5 such that v4L5 ∈ S4C5 or nC1v4L5 > 0 is being counted twice. Path joining. Because the labeling algorithm is bidi- rectional, the labels must be joined to construct feasi- ble paths. Given two labels L, L′ such that v4L5= v4L′5 and q4L5 + q4L′5 ≤ Q + dv4L5, they will produce a fea- sible path (one that satisfies capacity constraints and such that its reduced cost is negative) if (i) min8q4L51 q4L′59≥ 4q4L5+ q4L′5− dv4L55/2; (ii) max8q4L51 q4L′59≤ 4q4L5+ q4L′5+ dv4L55/2; (iii) v4pred4L55 < v4pred4L′55; (iv) the reduced cost of the concatenated path P = 4L1L′5 is negative. D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP 98 INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS The first two conditions are the median conditions (Baldacci et al. 2008) that ensure that labels L and L′ are the closest possible to half of the load. The third condition ensures that if path P = 4L1L′5 is kept, then path P ′ = 4L′1L5 will be discarded. This way, symmet- ric or repeated paths will not be added to the master problem. The dynamic programming algorithm. Let us describe the labeling algorithm by means of a pseudocode. Let L0 be the label representing an empty path starting at the facility, such that all of the resources are set at their default values. Also, let us note that labels will be stored in buckets, and let B4q1v5 be the bucket storing labels L whose loads are q4L5 = q and such that v4L5= v. Algorithm 2 (2-cyc-SPPRC) 1. Compute functions f 1 g1  using DP. 2. B40105← 8L091V← 8i91R← ™. 3. repeat 4. Take node v from V and set V←V\8v9. 5. for q = 0 to Q/2 do 6. for all L ∈ B4q1v5 such that proc4L5= false do 7. Set proc4L5← true. 8. for all w ∈ Neighbors of v, w 6= 0 and q4L5+ dw ≤Q and pred4L5 6=w do 9. Create L′ such that v4L′5=w and pred4L′5= L. Update resources accordingly. 10. Apply fathoming rule and eventually discard L′. 11. Apply dominance rule and eventually discard L′. 12. if L′ has not been discarded then 13. Make B4q4L′51w5← B4q4L′51w5∪ 8L′9. 14. Apply dominance rule and eventually delete other labels in B4q4L′51w5. 15. Make V←V∪ 8w9. 16. end if 17. end for 18. end for 19. end for 20. until V= ™ 21. R← 84L1L′52 v4L5= v4L′51 q4L5+ q4L′5≤Q+ dv4L5, L and L′ satisfy (i)–(iv)}. 22. return R. 4.4.4. Computing Lower Bounds. When pricing problems are solved to optimality, it is possible to obtain a lower bound on the problem. This lower bound can then be used for fathoming the current node as well as for early termination criteria. The fol- lowing proposition provides a way of computing a lower bound on the CLRP. Proposition 5. Let c̄min be the minimum reduced cost at the current iteration for columns in ì, and let z̄ be the value of the master problem at the current iteration. Also, let Kmax be an upper bound on the number of vehicles that serve two or more customers. A valid lower bound for the CLRP is given by zLB = z̄+Kmaxc̄min0 (31) Proof. See the online supplement. For every candidate set H we use Kmax =Km4H5 as described in the §4.3. 4.5. Enumeration of Remaining Columns For each subset of facilities H generated during the first stage and not discarded after the second stage procedure, we let zLB and ‘ be the lower bound at the end of the column-and-cut generation and the dual variables associated to such lower bound. If pro- cedure ENUM-ESPPRC was successful to generate the column set P, we simply compute the reduced cost of columns in P and add to the master prob- lem those columns l such that c̄l < zUB − zLB. We then solve the resulting integer problem using a general- purpose solver such as CPLEX. If, however, we were not able to obtain set P, we first check whether the upper bound zUB improved during the column-and- cut generation procedure after the consideration of set H . In this case, we run algorithm ENUM-ESPPRC again but now with the updated upper bound, as the performance of algorithm ENUM-ESPPRC strongly depends on the gap zUB −zLB. Otherwise, we start the following procedure with the hope of getting a better upper bound (if any), and in the worst case it gives us a method for tightening the gap: (i) Let ã← 4zUB − zLB5/10. Set k ← 1. (ii) Let z′UB ← zLB +kã and try to generate all of the columns whose reduced costs are less than or equal to kã. If more than ãmax = 106 columns are found or if we run out of memory, we exit. Otherwise we go to step (iii). (iii) Solve the resulting integer problem to optimal- ity. If a new upper bound was found with value z∗ < zUB, set zUB ← z∗ and zLB ← min8zUB1 z′UB9. If zLB = z ′ UB then exit. Otherwise, if either z′UB < z ∗ or the problem was solved to optimality but no integer solution was found with value less than zUB, set zLB = z′UB. If k < 10 do k ← k+ 1 and go back to (ii). This method generalizes the one proposed by Baldacci et al. (2011) by artificially lowering the opti- mality gap and iteratively increasing it, thus reduc- ing the negative impact of an initial upper bound of poor quality. In the online supplement we describe the algorithm for solving the column enumeration problem. This algorithm is a variation of the elemen- tary SPPRC (ESPPRC) and we call it ENUM-ESPPRC. 4.6. Computational Issues We now make some observations that can help to accelerate the algorithm. D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS 99 4.6.1. Initial Set of Columns. An initial set of columns is required in column generation algorithms. Indeed, at every iteration of the CG, a feasible solu- tion of the master problem is needed for running the pricing algorithms. In our algorithm, we let the ini- tial set of columns contain only the single-customer variables y. Additionally, we also add slack and arti- ficial variables to the formulation so the problem will always have a feasible solution. 4.6.2. Stabilization of the Column Generation. With the aim of reducing the oscillation of the dual variables during the first iterations of the column gen- eration process, we use a box-pen method (du Merle et al. 1999) for stabilizing the duals of the degree constraints (10). For every set H ∈ I, the centers are initially set to the optimal dual variables of the degree constraints (1) obtained after the solution of formulation (TIF) by partial branch and cut. 4.6.3. Column Pool Management. For some in- stances, the quantity of columns added can be huge and, moreover, most of them will be useless. In fact, it is known that at the beginning of the column gener- ation process, many columns are generated that soon will become nonbasic for the rest of the algorithm. We keep a pool of columns and keep track of the number of consecutive iterations that columns have been nonbasic. At every 30 iterations we check and delete all columns having been inactive for more than 30 iterations. Note that after the creation of set P during the column-and-cut generation procedure, the columns deleted from the problem must be inserted back into P. 4.6.4. Memory Management. The dynamic pro- gramming algorithms can be very demanding in terms of memory. In fact, every new created label needs to be allocated in memory. In this context, the new and delete operators of C++ (or malloc and free operators in the case of C) can be very inefficient. We have decided to manage our own memory pool, in which dynamic memory is allocated in chunks of 400 MB. The newly created labels are thus allocated inside the previously allocated memory. 5. Computational Experiments We have run our method on an Intel Xeon E5462, 3.0 GHz processor with 16 GB of memory. The code was compiled with the Intel C++ compiler v11.0 and executed on Linux, kernel 2.6. Linear and integer programs were solved by CPLEX 12.2. The pricing algorithms 2-cyc-SPPRC and ENUM-ESPPRC have been coded in C++ using the same compiler as before. The algorithm has been tested over five sets of instances from the literature, containing in total 71 instances. The first family (F1) has been adapted by Barreto (2004) from other vehicle routing prob- lems in the literature and contains 16 instances with capacitated vehicles and facilities. The second set of instances (F2) has been developed by Belenguer et al. (2011) and contains 30 instances with capacitated vehicles and facilities. The third set of instances (F3) has been introduced by Akca et al. (2009) and contains 12 instances with capacitated vehicles and facilities. The fourth set of instances (F4) has been intro- duced by Tuzun and Burke (1999) and contains nine instances with capacitated vehicles and unca- pacitated facilities. The fifth and last set of instances (F5) has been introduced by Baldacci et al. (2011) and contains four instances with capacitated vehicles and uncapacitated facilities. The dimensions of the instances vary from very small instances with 12 cus- tomers and two facilities up to very large instances with 199 customers and 14 facilities. We compare our results against those obtained by other exact algo- rithms, namely, the methods of Belenguer et al. (2011), Contardo et al. (2011), and Baldacci et al. (2011). For every instance we use as upper bound the best solu- tion available in the literature. In Table 1 we report the aggregate results obtained by our algorithm on each family of instances and for each of the three stages in our method. The columns in this table are labeled as follows: (i) Family: family of instances. (ii) gap1, t1: gap obtained (in %) and CPU time (in minutes) taken by the first stage procedure, namely, the enumeration of candidate subsets by par- tial branch and cut. The gap is computed as 4z∗ −zLB15/ z∗ × 100. (iii) gap2, t2: gap obtained (in %) and CPU time (in minutes) taken by the second-stage procedure, namely, the column-and-cut generation procedure to solve the linear relaxation of problem (SPF). (iv) gap31 t3: gap obtained (in %) and CPU time (in minutes) taken by the third-stage procedure, namely, the column enumeration algorithm. (v) t: overall CPU time (in minutes). (vi) opt: number of instances that were solved to optimality on the total number of instances in the family. As shown by Table 1, our algorithm is capable of solving 58 of the 71 instances considered in our study. Moreover, all instances of families F1, F3, and F5 (32 in total) are solved to optimality. For families F2 and F4, we report separately the aver- age results over all instances and over those that were solved to optimality or not. By looking more carefully at the detailed results (see the online supple- ment), one can see that instances Chr-75×10 ba, ppw- 50×5-2b, ppw-100×5-2b, and ppw-200×10-3a were solved to optimality for the first time, and that we have improved the best feasible solution for three D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP 100 INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS Table 1 Aggregate Results of the Proposed Method Family gap1 t1 gap2 t2 gap3 t3 t opt F1 2012 2076 0003 13027 0000 0001 16004 16/16 F2 all 3060 10000 0050 214067 0037 180044 405011 20/30 F2 solved 2087 3064 0013 19053 0000 6097 30013 F2 unsolved 5004 22074 1025 604095 1010 527038 11155007 F3 2065 0003 0002 0006 0000 0000 0010 12/12 F4 all 4024 52085 0027 11237038 0020 229063 11519086 6/9 F4 solved 3026 15021 0000 115013 0000 0000 130034 F4 unsolved 6021 128013 0080 31481086 0060 688090 41298089 F5 9085 860056 0013 91696040 0000 0019 101557013 4/4 58/71 more instances (ppw-100x5-0b, P113112, and P131112). Finally, our method is able to solve all instances with 85 customers or less. In Table 2 we compare our average results against those of the branch-and-cut methods of Belenguer et al. (2011) (column labeled BBPPW) and Contardo et al. (2011) (column labeled CCG-BC). A detailed comparison is presented in the online supplement. In the case of method CCG-BC we consider the branch-and-cut method on the two-index vehicle-flow formulation only, and with a maximum computing time of two hours. The results for the method intro- duced in this paper appear in the columns labeled CCG-CG. In Table 2, columns labeled gaplr and tlr stand for the gap (in %) and CPU time (in min- utes), respectively, obtained after solving the linear relaxation of formulation (TIF) by branch and cut. Columns labeled gap and t stand for the final gap and CPU time, respectively. Columns labeled opt rep- resent the number of instances solved to optimality by each method. Finally, columns labeled gap1 and t1 are as before. Results in bold characters indicate that a method has produced a better average gap than the other two, or that it solved more instances. Note that the instances considered for method BBPPW are a subset of those considered by CCG-BC, which are themselves a subset of the instances considered in this study. To make a fair comparison, two different Table 2 Average Comparison with Branch-and-Cut Methods BBPPW CCG-BC CCG-CG Family gaplr tlr gap t opt gaplr tlr gap t opt gap1 t1 gap t opt F†1 5032 0008 0050 5094 8/9 4091 0005 0036 13079 8/9 1052 0092 0000 8028 9/9 F‡1 5005 0011 0081 46067 10/16 2012 2076 0000 16004 16/16 F†2 5039 0006 0070 56028 6/12 4085 0013 0055 60084 6/12 2063 0007 0000 1063 12/12 F‡2 4091 5005 1022 89086 7/24 3009 2049 0027 176095 19/24 F3 6045 0002 0000 6011 12/12 6005 0002 0000 4082 12/12 2065 0003 0000 0010 12/12 26/33 26/33 33/33 29/52 47/52 †Restricted to instances reported by methods BBPPW and CCG-BC. ‡Restricted to instances reported by method CCG-BC. averages are reported. The first is the one restricted to the instances considered by method BBPPW that are also considered by CCG-BC and this method. The second one is restricted to the instances con- sidered by CCG-BC that are also included in this study. As shown in this aggregate table, our partial branch-and-cut produces better bounds than the flow- based algorithms at the root node. This is not surpris- ing because this procedure uses the code of CCG-BC for doing a partial branch and bound on the loca- tion variables. However, note that methods BBPPW and CCG-BC provide tighter bounds at the linear relaxation than our method on a few instances (ppw- 50x5-2a, ppw-50x5-2b, ppw-100x5-0a). This counterin- tuitive result is because our partial branch and cut does not include all valid inequalities used by the other two methods, notably framed capacity inequali- ties, multistar inequalities, hypotour inequalities, and location-routing comb inequalities, which have not been included to speed up the algorithm. When looking at the final gaps and CPU times, one may observe that our method is able to produce tighter gaps than the other two. Although a direct com- parison of computing times can be difficult (because each algorithm was run on a different machine), it is worth noting that our method was much faster on the instances of family F3 and on the subset of instances of family F2 considered by BBPPW. Moreover, we can solve all instances considered by method BBPPW and 47 out of 52 instances when including the ones considered by method CCG-BC, 18 more than the latter method. Finally, we compare the proposed methodology against the column generation method of Baldacci et al. (2011). In Table 3 we report the average results of the method of Baldacci et al. (2011) under the columns labeled BMW, but we retain the label CCG-CG for our method. Again, detailed results can be found in the online supplement. Note that, although the second and third stages in both methods are similar, the D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS 101 Table 3 Average Comparison with Column Generation Method BMW CCG-CG Family gap1 t1 gap2 t2 gap3 t3 t opt gap1 t1 gap2 t2 gap3 t3 t opt F†1 4065 6026 0007 14097 0000 1079 23002 15/15 2012 2076 0003 13027 0000 0001 16004 15/15 F†2 all 3057 5034 0076 2061 0047 85040 93035 17/24 3013 8024 0042 44033 0028 101004 153062 19/24 F†2 solved 3056 3008 0041 1071 0000 15077 20057 2071 1067 0014 2077 0000 6034 10078 F†2 unsolved 3044 6087 1064 4037 1064 283062 294086 4066 33037 1053 151053 1035 445047 630038 F3 4060 1098 0022 0077 0000 0001 2076 12/12 2065 0003 0002 0006 0000 0000 0010 12/12 F4 all 6080 65069 0048 129022 0040 153076 348068 6/9 4024 52085 0027 11237038 0020 229063 11519086 6/9 F4 solved 6010 78057 0012 106055 0000 1056 186068 3026 15021 0000 115013 0000 0000 130034 F4 unsolved 8020 39094 1020 174057 1020 458018 672069 6021 128013 0080 31481086 0060 688090 41298089 F5 7020 101059 0025 21439011 0000 26067 21567037 4/4 9085 860056 0013 91696040 0000 0019 101557015 4/4 54/64 56/64 †Restricted to instances reported by method BMW. procedure used by Baldacci et al. (2011) in the first stage of their algorithm is based on the solution of a relaxation of the set-partitioning formulation, whereas in our case it is based on the two-index vehicle-flow formulation of the problem that is solved by par- tial branch and cut. Results in bold characters sig- nify that a gap dominates the other, or a method is able to solve more instances than its counterpart. As seen in Table 3, our method is able to produce tighter bounds than that of Baldacci et al. (2011) at every stage for most families of instances considered in our study. Our partial branch-and-cut procedure is quite effective whenever branching decisions on the loca- tion variables have a significant impact on either the bounds or the feasibility of the problem, obtaining smaller gaps than those of Baldacci et al. (2011) on most instances. Regarding the column-and-cut gener- ation procedure, again our algorithm obtains gaps that are at least as good as those obtained by method BMW in most instances considered. This shows the strength of the set-partitioning formulation with the additional cuts. The third stage of the algorithm, although it can be very time consuming, is shown to be effective for solving instances Chr-75×10ba and ppw-200×10- 3a in which the initial upper bounds are significantly improved during this procedure. Overall, our algo- rithm is able to solve four open instances (two more than method BMW when restricted to the instances reported in their article), and improves the best known feasible solution on three other instances. For the instances that remain unsolved, our method always produces stronger lower bounds. However, for the instances in family F5 our method is outperformed by that of Baldacci et al. (2011). The overall results sug- gest that our method is competitive against that of Baldacci et al. (2011). This is the result of several refine- ments with respect to their method, namely, the use of the new cuts, as well as the use of efficient pric- ing algorithms that properly handle these new cuts. This includes the use of stronger fathoming proce- dures based on the solution of a 2-cyc-SPPRC with resources. 6. Concluding Remarks In this paper, we have presented an exact method for solving the CLRP. The methodology consists in formulating the CLRP as a set-partitioning problem that is solved in three stages: In a first stage we consider the two-index formulation and branch on the location variables. This strategy works well for instances where branching decisions on the location variables have a significant impact on the feasibility or the bound at the resulting nodes in the branch- ing tree. The remaining gap is then closed in the two following stages. In the second stage, the linear relax- ation of the set-partitioning formulation is solved by means of column-and-cut generation, which usually provides tight gaps. In the third and final stage, we try to enumerate all columns that could potentially lead to an improvement of the best known solution and solve the resulting integer program by branch and bound. The algorithm proposed in this paper is able to produce the tightest gaps on a large num- ber of instances. In addition, it has solved to opti- mality four previously open instances and improved the best known feasible solution for three additional ones. The methodology can be easily adapted to solve other routing problems. For instance, it would be interesting to measure the impact of y-SCC and SFrCI cuts on solving hard instances of the CVRP. With respect to the pricing algorithm introduced in this paper, the consideration of SDEG cuts allows lower bounds that are comparable to those obtained when pricing on elementary routes in a fraction of the computational effort. Indeed, in most cases only a fraction of SDEG cuts need to be added to the mas- ter problem to obtain significant improvements in the lower bound. Moreover, we show how to take advan- tage of this pricing problem in the computation of D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d. Contardo, Cordeau, and Gendron: An Algorithm Based on Cut-and-Column Generation for CLRP 102 INFORMS Journal on Computing 26(1), pp. 88–102, © 2014 INFORMS tight fathoming rules that speed up the whole algo- rithm. Further research related to the methodology introduced in this paper should address the devel- opment of new cutting planes for the set-partitioning formulation and to adapt some of them to other routing problems. Supplemental Material Supplemental material to this paper is available at http://dx .doi.org/10.1287/ijoc.2013.0549. Acknowledgments The authors thank Simon Spoorendonk, Enrico Bartolini, and two anonymous referees for their valuable comments. The authors gratefully acknowledge the financial support from the Natural Sciences and Engineering Research Coun- cil of Canada (NSERC) and Le fonds québécois de la recherche sur la nature et les technologies (FQRNT) [Grants 184122 and 227837-09]. References Akca Z, Berger RT, Ralphs TK (2009) Modeling and solving loca- tion, routing, and scheduling problems. Chinneck JW, ed. Proc. Eleventh INFORMS Comput. Soc. Meeting (Charleston, SC), 309–330. Archetti C, Bouchard M, Desaulniers G (2011) Enhanced branch and price and cut for vehicle routing with split deliveries and time windows. Transportation Sci. 45:285–298. Augerat P (1995) Approche polyhédrale du problème de tournées de véhicules. Ph.D. thesis, Institut National Polytechnique de Grenoble, France. Balas E, Padberg MW (1976) Set-partitioning: A survey. SIAM Rev. 18:710–760. Baldacci R, Christofides N, Mingozzi A (2008) An exact algo- rithm for the vehicle routing problem based on the set par- titioning formulation with additional cuts. Math. Programming 115:351–385. Baldacci R, Mingozzi A, Wolfler Calvo R (2011) An exact method for the capacitated location-routing problem. Oper. Res. 59: 1284–1296. Baldacci R, Bartolini E, Mingozzi A, Roberti R (2010) An exact solu- tion framework for a broad class of vehicle routing problems. Comput. Management Sci. 7:229–268. Barreto S (2004) Análise e modelização de problemas de localização-distribuição. Ph.D. thesis, University of Aveiro, Campus Universitário de Santiago, Aveiro, Portugal. Belenguer JM, Benavent E, Prins C, Prodhon C, Wolfler Calvo R (2011) A branch-and-cut algorithm for the capacitated location routing problem. Comput. Oper. Res. 38:931–941. Billionet A, Elloumi S, Djerbi LG (2005) Designing radio-mobile access networks based on synchronous digital hierarchy rings. Comput. Oper. Res. 32:379–394. Boland B, Dethridge J, Dumitrescu I (2006) Accelerated label set- ting algorithms for the elementary resource constrained short- est path problem. Oper. Res. Lett. 34:58–68. Christofides N, Mingozzi A, Toth P (1981) Exact algorithms for the vehicle routing problem, based on spanning tree and shortest path relaxations. Math. Programming 20:255–282. Contardo C, Cordeau J-F, Gendron B (2011) A computational com- parison of flow formulations for the capacitated location- routing problem. Technical report CIRRELT-2011-47, Université de Montréal, Canada. Desaulniers G, Lessard F, Hadjar A (2008) Tabu search, partial elementarity, and generalized k-path inequalities for the vehi- cle routing problem with time windows. Transportation Sci. 42:387–404. Desrochers M, Desrosiers J, Solomon M (1992) A new optimization algorithm for the vehicle routing problem with time windows. Oper. Res. 40:342–354. du Merle O, Villeneuve D, Desrosiers J, Hansen P (1999) Stabilized column generation. Discrete Math. 194:229–237. Feillet D, Gendreau M, Rousseau L-M (2007) New refinements for the solution of vehicle routing problems with branch and price. INFOR 45:239–256. Fukasawa R, Longo H, Lysgaard J, Poggi de Aragão M, Reis M, Uchoa E, Werneck RF (2006) Robust branch-and-cut-and-price for the capacitated vehicle routing problem. Math. Programming Ser. A 106:491–511. Gunnarsson H, Rönnqvist M, Carlsson D (2006) A combined ter- minal location and ship routing problem. J. Oper. Res. Soc. 57:928–938. Hansen PH, Hegedahl B, Hjortkjaer S, Obel B (1994) A heuris- tic solution to the warehouse location-routing problem. Eur. J. Oper. Res. 76:111–127. Jepsen M, Petersen B, Spoorendonk S, Pisinger D (2008) Subset-row inequalities applied to the vehicle-routing problem with time windows. Oper. Res. 56:497–511. Labbé M, Laporte G (1986) Maximizing user convenience and postal service efficiency in post box location. Belgian J. Oper. Res., Statist. Comput. Sci. 26:21–35. Lin CKY, Chow CK, Chen A (2002) A location-routing-loading problem for bill delivery services. Comput. Oper. Res. 43:5–25. Liu SC, Lin CC (2005) A heuristic method for the combined location routing and inventory problem. Internat. J. Advanced Manufac- turing Tech. 26:372–381. Lysgaard J, Letchford AN, Eglese RW (2004) A new branch-and-cut algorithm for the capacitated vehicle routing problem. Math. Programming 100:423–445. Perl J, Daskin MS (1985) A warehouse location-routing problem. Transportation Res. B 19:381–396. Prins C, Prodhon C, Wolfler Calvo R (2006a) A memetic algorithm with population management (MA — PM) for the capacitated location-routing problem. Gottlieb J, Raidl GR, eds. EvoCOP 2006, Lecture Notes in Computer Science, Vol. 3906 (Springer- Verlag, Berlin), 183–194. Prins C, Prodhon C, Wolfler Calvo R (2006b) Solving the capaci- tated location-routing problem by a GRASP complemented by a learning process and path relinking. 4OR 4:221–238. Prins C, Prodhon C, Ruiz A, Soriano P, Wolfler Calvo R (2007) Solving the capacitated location-routing problem by a coop- erative Lagrangean relaxation-granular tabu search heuristic. Transportation Sci. 41:470–483. Righini G, Salani M (2008) New dynamic programming algorithms for the resource constrained elementary shortest path problem. Networks 51:155–170. Tuzun D, Burke LI (1999) A two-phase tabu search approach to the location routing problem. Eur. J. Oper. Res. 116:87–99. Wu T-H, Low C, Bai J-W (2002) Heuristic solutions to multi-depot location-routing problems. Comput. Oper. Res. 29:1393–1415.D ow nl oa de d fr om in fo rm s. or g by [ 12 8. 12 2. 25 3. 21 2] o n 15 O ct ob er 2 01 4, a t 0 4: 19 . Fo r pe rs on al u se o nl y, a ll ri gh ts r es er ve d.