Roberto Baldacci
email: r.baldacci@unibo.it

Sandra Ulrich Ngueveu
email: sandra.ulrich.ngueveu@laas.fr

Roberto Wolfler Calvo

R Wolfler Calvo

The Vehicle Routing Problem with Transhipment Facilities

Keywords: transhipment facilities, dual ascent heuristic, column-and-cut generation The VRPTF

This paper proposes an exact method for solving an optimization problem arising in several distribution networks, where customers can be served either directly, using vehicle routes from a central depot, or through transhipment facilities. The problem consists of optimizing the following interdependent decisions: selecting transhipment facilities, allocating customers to these facilities and designing vehicle routes emanating from a central depot to minimize the total distribution cost. This problem is called the Vehicle Routing Problem with Transhipment Facilities (vrptf). The paper describes two integer programming formulations for the vrptf, an edge-flow based formulation and a Set Partitioning (SP) based formulation. The LP-relaxation of the two formulations are further strengthened with the addition of different valid inequalities. Moreover, two new route relaxations that are used by dual ascent heuristics to find near-optimal dual solutions of the LP-relaxation of the SP model are described. The valid inequalities and the route relaxations are used in a branchand-cut-and-price approach to solve the problem to optimality. The proposed method is tested on a large family of instances, including real-world instances, and the computational results obtained indicate the effectiveness of the proposed method.

Introduction

In several distribution networks the shipment to a customer is performed either directly, using vehicle routes emanating from a central depot, or through intermediate depots or transhipment facilities. In the latter case, the shipment is first delivered to a transhipment facility by a vehicle route, and then it is successively delivered to the final customer. Transhipment facilities provide a way to consolidate

Literature review

The vrptf generalizes the well-known Capacitated Vehicle Routing Problem (cvrp). In the cvrp, a fleet of identical vehicles located at a central depot has to be optimally routed to supply a set of customers with known demands. Each vehicle performs at most one route, each customer must be visited exactly once, and the total demand of the customers visited by a route cannot exceed the vehicle capacity. The book edited by [START_REF] Toth | Vehicle Routing: Problems, Methods, and Applications[END_REF] provides a comprehensive overview of exact methods for the cvrp and other variants.

As far as the authors know, the vrptf has never been addressed in the literature. Closely related problems to the vrptf are the Capacitated m-Ring-Star Problem (cmrsp), the Multiple Vehicle Traveling Purchaser Problem (mvtpp), the Two-Echelon Capacitated Vehicle Routing Problem (2ecvrp), and the Location Routing Problem (lrp). The cmrsp, introduced in [START_REF] Baldacci | The Capacitated m-Ring-Star Problem[END_REF], arises in the design of urban optical telecommunication networks and it consists of designing a set of rings that pass through a telephone exchange and through some transition points (also called steiner nodes) and/or users. Each nonvisited user must be assigned to a visited point or to a user. The number of users visited and assigned to a ring is bounded by the capacity of the ring. The objective is to minimize the total routing cost plus the assignment costs. The special case of the cmrsp arising when the users can be assigned only to steiner nodes, can be solved as a vrptf with unit demands. The mvtpp described by Riera-Ledesma and Salazar-González (2012) models a family of routing problems combining stop selection and bus route generation. The problem consists of choosing a set of bus stops to which users are assigned, and simultaneously designing bus routes visiting such stops. The total number of users assigned to the stops of a route cannot exceed the seat capacity of a bus. The objective is to minimize the total length of all routes plus the total assignment cost. The undirected version of the mvtpp is equivalent to the vrptf with the additional constraint imposing that the customers can only be assigned to facilities (or bus stops) and cannot be visited by a route. Both [START_REF] Baldacci | The Capacitated m-Ring-Star Problem[END_REF] and [START_REF] Riera-Ledesma | Solving school bus routing using the multiple vehicle traveling purchaser problem: A branch-and-cut approach[END_REF] proposed branch-and-cut approaches for the solution of the cmrsp and mvtpp, respectively. Recently, Riera-Ledesma and González (2013) also proposed a branch-and-cut-and-price algorithm for the mvtpp. The 2e-cvrp is a two-level distribution system where the deliveries to customers from a depot are managed through intermediate capacitated depots, called satellites. The first level consists of vehicle routes visiting satellites only whereas the second level routes supply all customers. The main difference between the vrptf and the 2e-cvrp is that in the vrptf a customer can be either visited on a route or assigned to a facility, whereas in the 2e-cvrp each customer is visited once by exactly a second level route. The 2e-cvrp model is particularly useful when the facilities are part of the same company owing the main depot whereas in the vrptf model the facilities are generally owned by third-party contractors, which are in charge of delivering to the final customers the quantity consolidated at the facilities. Exact methods for the 2e-cvrp have been proposed by [START_REF] Jepsen | A branch-and-cut algorithm for the symmetric two-echelon capacitated vehicle routing problem[END_REF] and [START_REF] Baldacci | An Exact Algorithm for the Two-Echelon Capacitated Vehicle Routing Problem[END_REF]. The lrp is a special case of the 2e-cvrp and consists of opening a set of depots and designing a set of routes for each opened depot so that the total load of the routes operated from a depot does not exceed its capacity and each customer is visited by exactly one route. The objective is to minimize the sum of the fixed costs of the opened depots and the costs of the routes operated from the depots. A recent review of location routing problem variants and heuristic and exact algorithms can be found in [START_REF] Prodhon | A survey of recent research on location-routing problems[END_REF].

Another related problem to the vrptf is the Multi-Vehicle Covering Tour Problem (m-CTP) introduced by [START_REF] Hachicha | Heuristics for the multi-vehicle covering tour problem[END_REF]. In the m-CTP two sets of locations are given. The first set, consists of potential locations at which some vehicles may stop, and the second set are locations not actually on vehicle routes, but within an acceptable distance from a vehicle route. The m-CTP consists of determining a set of total minimum length vehicle routes on a subset of the first set of locations, subject to side constraints, such that every location of the second set is within a prespecified distance from a route. [START_REF] Há | An exact algorithm and a metaheuristic for the multi-vehicle covering tour problem with a constraint on the number of vertices[END_REF] proposed a branch-and-cut for the variant named the m-CTP-p where an upper bound on the number of vertices per route is given with a parameter p and the m number of vehicles used is a decision variable. The m-CTP-p is equivalent to the vrptf with the additional constraint imposing that the customers can only be assigned to facilities and cannot be visited by a route, and that all assignment costs are equal to zero.

The vrptf does not require any specific synchronization of incoming and outgoing vehicles at the facilities. In some practical applications, a correct synchronization can be required and in this case the facilities are generally referred as cross-docking facilities. For an overview of the cross-docking concept and extensive review of the existing literature the reader is referred to [START_REF] Belle | Cross-docking: State of the art[END_REF]. In this context, a generic class of vrps that has recently received attention in the literature is the class of vrps with Multiple Synchronization Constraints (vrpmss). vrpmss exhibit synchronization requirements between the vehicles, concerning spatial, temporal, and load aspects. A review of vrpms presenting a classification of different types of synchronizations and a discussion about heuristic and exact algorithms can be found in [START_REF] Drexl | Synchronization in vehicle routing -a survey of vrps with multiple synchronization constraints[END_REF].

Contributions of this paper

This paper addresses a new problem of practical relevance and proposes both heuristic and exact methods for its solution. More specifically, we introduce a two-index formulation (T I) and we describe different valid inequalities for it, both by adapting those already proposed for the cmrsp, and by introducing new ones specific for the vrptf. We also describe lower bounds derived from a set-partitioning based formulation (SP ) of the problem, and computed using two efficient dual ascent heuristics that use two new route relaxations, called q- * route and ng- * route, respectively. The proposed methods have been tested on a large family of instances, including both instances derived from the literature and real-world instances. The computational results show that real-world instances with up to 142 customers and 18 facilities were solved to optimality and that high quality solutions were computed for instances with up to 164 customers. In addition, tight lower bounds were computed, with average percentage deviations equal to 98.7% and 97.4% for real-word and literature-based instances, respectively. This paper is organized as follows. The next section formally introduces the vrptf and presents formulation T I for which different valid inequalities are described in Section 3. Section 4 presents formulation SP and lower bounds based on its LP-relaxation; some properties of the LP-relaxation of SP are also investigated in the section. A bounding method used to compute a lower bound on the vrptf is described in Section 5. Section 6 describes the exact method used to solve the vrptf to optimality together with two heuristic algorithms. Section 7 reports computational results, and concluding remarks are given in Section 8.

Problem description and Two-Index (TI) formulation

This section describes the vrptf and presents a edge-flow based formulation to model it.

The vrptf is defined on a mixed graph G = (V, E ∪ A), where V = {0} ∪ V is the node set, E = {{i, j} : i, j ∈ V, i = j} is the edge set, and A is the arc set. Node set V is partitioned into two subsets: V C = {1, . . . , n C } containing a node for each customer and V F = {n C + 1, . . . , n C + n F } containing a node for each transhipment facility. Node 0 represents a central depot. Each customer i ∈ V C requires a supply of q i units from the depot (we assume q i = 0, ∀i ∈ {0} ∪ V F ) that can be delivered either directly from a vehicle route emanating from the depot or through a facility selected from a set F i ⊆ V F of facilities to which customer i can be assigned. Set A represents the possible assignments between customers and facilities, i.e., A = {(i, j) : i ∈ V C , j ∈ F i }. Set E is the set of possible route edges, each edge e = {i, j} ∈ E is associated with a non-negative routing cost r e = r {i,j} , while each arc (i, j) ∈ A is associated with a non-negative assignment cost d ij . Henceforth, if e connects the two nodes i and j then {i, j} and e will be used interchangeably to denote the same edge.

A route is defined by a pair (R, A ) where R = (0, i 1 , . . . , i r , 0), r ≥ 1, is a simple cycle in G passing through the depot, visiting nodes V (R) = {i 1 , . . . , i r } ⊆ V , and A ⊆ A are assignments between customers of V C \ V (R) and nodes of V (R) ∩ V F . Notice that if r = 1 then route R represents the single-node route R = (0, i 1 , 0). We say that a customer i is assigned to a route R if it is either visited by the simple cycle (i.e., i ∈ V (R)) or it is connected to a node of the route representing a facility (i.e., a node j ∈ V (R) ∩ V F exists such that (i, j) ∈ A ). The total load of a route is computed as the sum of the demands of the customers assigned to the route. The route is feasible if its total load does not exceed the vehicle capacity Q. The cost of a route is equal to the sum of the routing costs of the edges forming the route plus the sum of the assignment costs of the arcs in A .

The aim of vrptf is to design a set of routes so that each customer is assigned to exactly one route, each intermediate facility is visited at most once and the sum of the route costs is minimized.

We will use the following notation throughout. For any S ⊆ V , let V C (S) = S ∩ V C and V F (S) = S ∩ V F denote the set of customers and of facilities in S, respectively. Let F i (S) = V F (S) ∩ F i denote the set of facilities in S associated with customer i ∈ V C . Also for any set S ⊆ V , define δ(S) = {{i, j} ∈ E : i ∈ S, j / ∈ S} (if S = {i}, we simply write δ(i) instead of δ({i})).

Let x e be an integer variable which takes value in {0, 1}, ∀e ∈ E \ {{0, j} : j ∈ V } and value in {0, 1, 2}, ∀e ∈ {{0, j} : j ∈ V }. Notice that x {0,j} = 2 when the single-node cycle R = (0, j, 0) is selected in the solution. For each arc (i, j) ∈ A, let z ij be a binary variable which is equal to 1 if and only if customer i is assigned to node j. Moreover, for each i ∈ V , let y i be a binary variable which is equal to 1 if and only if node i is on a route. Formulation T I is as follows:

(T I) min e∈E r e x e + (i,j)∈A d ij z ij (1) s.t. e∈δ(i)
x e = 2y i , ∀i ∈ V (2)

y i + j∈Fi z ij = 1, ∀i ∈ V C (3) e∈δ(S) x e ≥ 2 Q   i∈V C (S) q i y i + (i,j)∈A:j∈V F (S) q i z ij   , ∀S ⊆ V : S = ∅ (4) x e ∈ {0, 1}, ∀e ∈ E \ {{0, j} : j ∈ V } (5) x e ∈ {0, 1, 2}, ∀e ∈ {{0, j} : j ∈ V } (6) z ij ∈ {0, 1}, ∀(i, j) ∈ A (7) y i ∈ {0, 1}, ∀i ∈ V . (8)
Constraints (2) impose that the degree of each node i ∈ V is 2 if the node is on a route. Constraints (3) state that a customer i ∈ V C is either on a route or is assigned to one of its facilities. Inequalities (4) are the fractional route capacity inequalities (FrCC). These constraints, within the integrality of x, z and w variables, impose that for a given subset S of nodes, at least ( i∈S q i y i + (i,j)∈A:j∈S q i z ij )/Q routes are needed to visit the customers assigned to nodes in S.

Strengthening the LP-relaxation of formulation T I

A number of valid inequalities can be used to improve the quality of the lower bound obtained from the LP-relaxation of formulation T I. In this section, we first derive valid inequalities by extending the results proposed for the cmrsp by [START_REF] Baldacci | The Capacitated m-Ring-Star Problem[END_REF] to the vrptf. Then, a new class of valid inequalities specifically devised for the vrptf is introduced. The separation procedures for different valid inequalities are then described in Section 5.2.

Simple valid inequalities are the following:

(i) x {i,j} ≤ y j , i ∈ V C , j ∈ V C , i = j; (ii) x {i,j} ≤ y j , i ∈ V F , j ∈ V , i = j; (iii) x {i,j} + z ij ≤ y j , i ∈ V C , j ∈ F i , (iv) y j ≤ i∈V C :j∈Fi z ij , ∀j ∈ V F .
Further, the following inequalities are also valid. a) Connectivity inequalities (CI):

e∈δ(S) x e ≥ 2   y i + j∈V F (S)∩Fi z ij   , ∀S ⊆ V , ∀i ∈ V C (S), S = ∅. (9) b) Multistar inequalities (MI): e∈δ(S) x e ≥ 2 Q   i∈V C (S)
q i y i + (i,j)∈A:j∈V F (S)

q i z ij + i∈V C (S) j∈S q i x {i,j}   , ∀S ⊆ V , S = ∅. ( 10
)
where S = V \ S.

The VRPTF R. Baldacci, S. U. Ngueveu, R. Wolfler Calvo c) Rounded capacity constraints I (RCI):

e∈δ(S) x e ≥ 2     i∈S:Fi⊆S q i /Q     , ∀S ⊆ V , V C (S) = ∅. (11) d) Rounded capacity constraints II (RCII): e∈δ(S) x e ≥ 2             i∈V C (S) q i y i + (i,j)∈A: j∈V F (S) q i z ij      /Q        , ∀S ⊆ V , S = ∅. ( 12
)
Notice that CI inequalities are not dominated by MI inequalities whereas MI inequalities dominate FrCC inequalities. RCII inequalities ( 12) are clearly nonlinear. In the next section, we describe two ways of linearizing inequalities (12). The first linearization extends to the vrptf a similar linearization proposed for the cmrsp, whereas the second one is new and it is based on mixed integer optimization.

Linearized versions of inequalities RCII

A first family of valid inequalities can be obtained using the following lemma, proposed by [START_REF] Baldacci | The Capacitated m-Ring-Star Problem[END_REF].

Lemma 1 Let m, n and o be three non-negative integer values with m > o and mod(m, o) = 0:

m -n o ≥ m o - n mod(m, o) . ( 13
)
The term i∈V C (S) q i y i + (i,j)∈A:j∈V F (S) q i z ij of RCII inequalities (12) can be rewritten as:

q(V C ) -   i∈V C (S) q i y i + (i,j)∈A:j∈V F (S) q i z ij   (14)
and by using Lemma 1, from expression ( 14) we obtain the following inequality valid for any S ⊆ V , S = ∅:

e∈δ(S) 1 2 x e ≥ q(V C ) Q - 1 mod(q(V C ), Q)   i∈V C (S) q i y i + (i,j)∈A:j∈V F (S) q i z ij   , ( 15
)
hereafter called RCII-a inequalities.

The term i∈V C (S) q i y i + (i,j)∈A:j∈V F (S) q i z ij of RCII inequalities (12) can also be rewritten as:

q(V C (S)) -      (i,j)∈A: i∈V C (S),j∈V F (S) q i z ij - (i,j)∈A: i∈V C (S),j∈V F (S) q i z ij      , ( 16
)
The VRPTF R. Baldacci, S. U. Ngueveu, R. Wolfler Calvo and by using Lemma 1 and by disregarding the term (i,j)∈A: i∈V C (S),j∈V F (S) q i z ij from ( 16) we get:

e∈δ(S) 1 2 x e ≥ q(V C (S)) Q - 1 mod(q(V C (S)), Q) (i,j)∈A: i∈V C (S),j∈V F (S) q i z ij , ( 17
)
hereafter called RCII-b inequalities.

Proposition 1 of the e-companion to this paper shows that there is no dominance between inequalities RCII-a and RCII-b.

The following lemma is based on mixed integer optimization. For a number m ∈ R, define m = m -m to be its fractional part.

Lemma 2 Let o ∈ R with ô > 0 and T = {m ∈ R, n ∈ Z : m + n ≥ o, m ≥ 0}. The following inequality is valid for T : m + ô n ≥ ô o . ( 18
)
Proof. The proof is provided in the e-companion to this paper.

Based on the above lemma, a second family of valid inequalities for the vrptf can be obtained using the following theorem.

Theorem 1 Let α e ≥ 0, ∀e ∈ E, β i ≥ 0, ∀i ∈ V and γ ij ≥ 0, ∀(i, j) ∈ A and consider the following inequality valid for formulation T I:

e∈E α e x e + i∈V β i y i + (i,j)∈A γ ij z ij ≥ o ( 19
)
where o ∈ R and ô > 0. Then the following inequality:

e∈E ϕ o (α e )x e + i∈V ϕ o (β i )y i + (i,j)∈A ϕ o (γ ij )z ij ≥ o (20)
where

ϕ o (m) = m + min m ô , 1 , m ∈ R, n ∈ R, ô > 0, is also a valid inequality for formulation T I.
Proof. The proof is provided in the e-companion to this paper.

Notice that, inequality (19) can be scaled by a rational number t thus obtaining the following valid inequality for formulation T I:

e∈E ϕ to (tα e )x e + i∈V ϕ to (tβ i )y i + (i,j)∈A ϕ to (tγ ij )z ij ≥ to . ( 21
)
Starting from inequalities (4) and substituting the right-and side according to expressions ( 14) and ( 16) we get:

e∈δ(S) 1 2 x e + i∈V C (S) q i Q y i + (i,j)∈A:j∈V F (S) q i Q z ij ≥ q(V C ) Q , ( 22
)
and e∈δ(S)

1 2 x e + (i,j)∈A: i∈V C (S),j∈V F (S) q i Q z ij - (i,j)∈A: i∈V C (S),j∈V F (S) q i Q z ij ≥ q(V C (S)) Q . ( 23
)
First of all, notice that for m, n ∈ R, mod(m, n

) = n((m/n) -m/n ). Then, by setting o = q(V C ) Q and as ϕ o ( 1 2 ) = min Q 2mod(q(V C ),Q) , 1 and ϕ o ( qi Q ) = min qi mod(q(V C ),Q) , 1 , ∀i ∈ V C , from
Theorem 1 and inequality ( 22) we obtain the following valid inequality:

e∈δ(S) min Q 2 mod(q(V C ), Q) , 1 x e ≥ q(V C ) Q - i∈V C (S) min q i mod(q(V C ), Q) , 1 y i - (i,j)∈A:j∈V F (S) min q i mod(q(V C ), Q) , 1 z ij . ( 24
)
Also from Theorem 1, by disregarding the negative term of inequality ( 23) we obtain:

e∈δ(S) min Q 2 mod(q(V C (S)), Q) , 1 x e ≥ q(V C (S)) Q - (i,j)∈A: i∈V C (S),j∈V F (S) min q i mod(q(V C (S)), Q) , 1 z ij . ( 25
)
We call inequalities ( 24) and ( 25) RCII-c and RCII-d inequalities, respectively. Inequalities RCII-c and RCII-d are stronger than the pure integer rounding inequalities obtained from inequalities ( 22) and ( 23). In addition, notice that the coefficients of variables {x e } in both inequalities (24) and (25) are greater than 0.5 and less than or equal to 1. If q i = 1, ∀i ∈ V C , inequalities RCII-a and RCII-b dominate inequalities RCII-c and RCII-d. In general, no dominance relations exist among the four types of inequalities RCII-a, RCII-b, RCII-c and RCII-d.

Lower bounds based on a Set-Partitioning (SP) formulation

In this section, we first describe a Set-Partitioning (SP) based formulation for the vrptf. Then, we investigate lower bounds based on the LP-relaxation of formulation SP . We introduce a theorem that is used to derive two dual ascent heuristics to find near-optimal dual solutions of the LP-relaxation of the SP model. Then, we describe how the valid inequalities described for the T I formulation in the previous sections can be used for strengthening the value of the LP-relaxation of formulation SP . Finally, we derive some properties of the LP-relaxation of formulation SP .

Let R be the index set of all feasible routes. Given a route ∈ R, we denote with R the sequence (i 1 = 0, i 2 , . . . , i r = 0) of the nodes visited by the route and with V C (R ) and V F (R ) the sets V C ∩ V (R ) and V F ∩ V (R ), respectively. In addition, V A (R ) denotes the customers of the route assigned to facilities in V F (R ). Let a i be a (0-1) binary coefficient equal to 1 if node i ∈ V (R ), 0 otherwise. In addition, let b j i be a (0-1) binary coefficient equal to 1 if customer i ∈ V A (R ) is assigned to node j ∈ V F (R ), 0 otherwise. Given a route , we denote with c its routing cost computed as

|R | h=2 r {i h-1 ,i h } ,
and with p its assignment cost computed as j∈V F (R ) i∈V A (R ) b j i d ij . Let ξ , ∈ R, be a (0-1) binary variable equal to 1 if and only if route is in the optimal solution.

Formulation SP is as follows:

(SP ) min

∈R (c + p )ξ (26) s.t. ∈R a i ξ = 1, ∀i ∈ V C (27) ∈R a i ξ ≤ 1, ∀i ∈ V F (28) ξ ∈ {0, 1}, ∀ ∈ R, (29)
where

a i = a i + j∈V F (R ) b j i , i ∈ V C , ∈ R.
In the formulation, constraints ( 27) and ( 28) impose that each customer is assigned exactly once and each facility is visited at most once, respectively.

We denote by LSP the LP-relaxation of formulation SP and by DSP the dual of LSP . The variables of DSP are given by the vector u = {u 0 , u 1 , .

λ i ≤ 0, ∀i ∈ V F , with constraint (28). Let R i = { ∈ R : a i > 0}. For each i ∈ V C compute: ν i = q i min ∈Ri (c + p ) -j∈V C a j λ j -j∈V F a j λ j j∈V C a j q j . ( 30
)
A feasible DSP solution u of cost z(DSP (λ)) is given by the following expressions:

u 0 = 0 and u i = ν i + λ i , ∀i ∈ V C , and u i = λ i , ∀i ∈ V F . ( 31
)
Proof. The proof is provided in the e-companion to this paper.

The pricing problem associated with formulation SP is a strongly N P-hard problem, since it requires finding minimum cost elementary routes over a graph with both positive and negative edge and arc costs. In the special case where V F = ∅, the pricing problem consists of finding capacitated elementary cycles, a strongly N P-hard problem (see [START_REF] Poggi | Chapter 3: New exact algorithms for the capacitated vehicle routing problem[END_REF].

Therefore, in practice we enlarge the set of routes R to contain also nonnecessarily elementary routes, i.e., coefficients a i are general nonnegative integer, thus a node can be visited in a route more than once and/or a customer can be assigned more than once to facilities of the routes. Although non-elementary routes are infeasible, this relaxation has the advantage that the pricing subproblem becomes solvable in pseudo-polynomial time (by dynamic programming). Moreover, Theorem 2 remains valid if the set of routes R is enlarged to contain also nonnecessarily elementary routes.

In Section 5, we introduce two route relaxations called q- * route and ng- * route, used by two dual ascent heuristics based on Theorem 2 to find near-optimal solutions of problem DSP . q- * route and ng- * route relaxations are based on route relaxations already proposed for the cvrp and on the observation that given a route R = (i 1 = 0, i 2 , . . . , i r = 0), a lower bound on its cost c + p can be computed as

|R | h=2 r {i h-1 ,i h } + j∈V F (R ) lb j , where lb j ≤ i∈V A (R ) b j i d ij . Each value lb j , j ∈ V F (R )
, can be computed as the minimum of the costs of all possible assignments of facility j involving customers in {i :∈ V C : j ∈ F i } with a total load q = i∈V A (R ) b j i q i .

Formulation LSP can be strengthened by adding valid inequalities derived for the T I formulation as follows. For each ∈ R, let coefficients η e be defined as follows: if is a route covering node h only, then η {0,h} = 2 and η {i,j} = 0, ∀{i, j} ∈ E \ {0, h}; if is not a single-node route, then η {i,j} = 1 for each edge {i, j} traversed by route R , and η {i,j} = 0 otherwise. Any feasible solution ξ of SP can be transformed into a feasible T I solution (x, z, w) by setting:

x e = ∈R η e ξ , ∀e ∈ E, ( 32
)
z ij = ∈R b j i ξ , ∀(i, j) ∈ A, ( 33
)
y i = ∈R a i ξ = 1 - j∈Fi ∈R b j i ξ , ∀i ∈ V C , and (34)
y i = ∈R a i ξ , ∀i ∈ V F . ( 35
)
The following theorem shows that any feasible solution of formulation LSP already satisfies some valid inequalities derived from formulation T I.

Theorem 3

The LP-relaxation of the SP formulation satisfies both CI and FrCI inequalities, and a weak form of MI inequalities.

Proof. The proof is provided in the e-companion to this paper.

Bounding procedure

This section presents a method for computing a lower bound on the vrptf which combines in sequence two dual ascent heuristics (see Section 5.1), and a column-and-cut generation method (see Section 5.2), all based on formulation LSP .

Dual ascent heuristics

The dual ascent heuristics are based on Theorem 2 where the set of routes R is enlarged with set R > containing also nonnecessarily elementary routes (i.e., R > ⊇ R). In particular, two different route relaxations are used, called q- * route and ng- * route, to compute lower bounds LB 1 and LB 2 on the vrptf, respectively. The two dual ascent heuristics are based on a column generation-like method, called CG for solving the following problem:

LCG = max λ {z(DSP (λ))}. ( 36
)
CG executes a number of macro-iterations to compute a dual solution u of the master problem DSP , defined by the route subset R ⊆ R > , and then CG solves problem (36) with a predefined number M axit2 of subgradient iterations to modify the penalty vector λ.

Route relaxation q- * route

q- * routes are based on the q-path relaxation proposed by [START_REF] Christofides | Exact algorithms for the vehicle routing problem based on spanning tree and shortest path relaxation[END_REF]. We define a q- * path as a nonnecessarily elementary partial route in G from depot 0 to node i ∈ V with a load equal to q. In a q- * path a node i ∈ V can be visited more than once and a customer i ∈ V C can be assigned more than once. In the following, we describe a dynamic programming algorithm for computing q- * paths, with the restriction that a q- * path can not contain loops formed by three consecutive nodes. Let f (q, i) be the cost of the least cost q- * path from node 0 to node i and let π(q, i) be the node immediately before i in the least cost path of value f (q, i). Let g(q, i) be the cost of the least cost q- * path from node 0 to node i, such that γ(q, i) = π(q, i), where γ(q, i) is the node immediately before i in the least cost path corresponding to g(q, i). For a given value of q, let h(i, j) be the cost of the least cost q- * path from 0 to j, with i ∈ V just before j and without loops. In addition, for each facility k ∈ V F , let lb k (q) be a lower bound on the assignment cost of any assignment of load q of customers to the facility k. lb k (q), for each k ∈ V F and q min ≤ q ≤ Q, can be computed as the optimal solution cost of the following knapsack problem KP (q, k):

(KP (q, k)) lb k (q) = min i∈V C :k∈Fi d ik χ i (37) s.t. i∈V C :k∈Fi q i χ i = q (38) χ i ∈ {0, 1}, ∀i ∈ V C : k ∈ F i . ( 39
)
We assume that lb k (q) = ∞ if problem KP (q, k) does not admit a feasible solution for the given pair q and k. For each q = q min , . . . , Q and i, j ∈ V , i = j, compute:

h(i, j) =        f (q -q j , i) + r {i,j} , if π(q -q j , i) = j g(q -q j , i) + r {i,j} , otherwise. , j ∈ V C min q min ≤ w ≤ Q f (q -w, i) + r {i,j} + lb j (w), if π(q -w, i) = j g(q -w, i) + r {i,j} + lb j (w), otherwise. , j ∈ V F (40)
Then, compute:

f (q, j) = min i∈V \{j} {h(i, j)} π(q, j) = i ( 41
)
where i is the node producing the above minimum,

g(q, j) = min i∈V \{j,i } {h(i, j)} γ(q, j) = i ( 42
)
where i is the node producing the above minimum. The functions are initialized as follows:

• f (q j , j) = r 0j , π(q j , j) = 0, j ∈ V C ; • f (q, j) = ∞, π(q, j) = 0, q = 0, . . . , Q ,q = q j , j ∈ V C ; • f (q, j) = r 0j + lb j (q), π(q, j) = 0, q = 0, . . . , Q, j ∈ V F ; • g(q, j) = ∞, γ(q, j) = 0, q = 0, . . . , Q, j ∈ V .
A q- * route is obtained from a q- * path ending in i by adding arc (i, 0).

The VRPTF R. Baldacci, S. U. Ngueveu, R. Wolfler Calvo

Route relaxation ng- * route

ng- * routes are based on the route relaxations proposed by [START_REF] Baldacci | An Exact Method for the Capacitated Location-Routing Problem[END_REF] for the cvrp. Let N i ⊆ V be a set of selected nodes for node i ∈ V (according to some criterion) such that N i i and |N i | ≤ Γ, where Γ is a parameter (e.g., Γ = 5, ∀i ∈ V , and N i contains i and the four nearest nodes to i).

With a forward path P = (0, i 1 , . . . , i k ) , we associate a set Π(P ) ⊆ V defined as:

Π(P ) = {i r : i r ∈ k s=r+1 N is , r = 1, . . . , k -1} ∪ {i k }. (43)
A forward ng- * path (N G, q, i) is a non-necessarily elementary partial route P = (0, i 1 , . . . , i k-1 , i k = i) starting from the depot with a load equal to q, ending at customer i, and such that N G = Π(P ), and i / ∈ Π(P ), where P = (0, i 1 , . . . , i k-1 ). Let f (N G, q, i) be the cost of a least-cost forward ng- * path (N G, q, i). The dynamic programming (DP) recursion for computing functions f (N G, q, i) is defined on a state-space graph H = (E , Ψ) defined as:

E = {(N G, q, i) : q i ≤ q ≤ Q, ∀N G ⊆ N i s.t. N G i, ∀i ∈ V } Ψ = {((N G , q , j), (N G, q, i)) : ∀(N G , q , j) ∈ Ψ -1 (N G, q, i), ∀(N G, q, i) ∈ E }, (44)
where

Ψ -1 (N G, q, i) = {(N G , q -q i , j) : ∀N G ⊆ N j s.t. N G j and N G ∩ N i = N G \ {i}, ∀j ∈ V \ {i}}, if i ∈ V C , and Ψ -1 (N G, q, i) = {(N G , q , j) : 0 ≤ q ≤ q -min i∈V C {q i }, ∀N G ⊆ N j s.t. N G j and N G ∩ N i = N G \ {i}, ∀j ∈ V \ {i}}, if i ∈ V F .
The DP recursion for computing functions f (N G, q, i), for each state (N G, q, i) ∈ E is as follows:

i) i ∈ V F :f (N G, q, i) = min (N G ,q ,j)∈Ψ -1 (N G,q,i) {f (N G , q , j) + r {j,i} + lb i (q -q )}, ∀(N G, q, i) ∈ E , ii) i ∈ V C : f (N G, q, i) = min (N G ,q ,j)∈Ψ -1 (N G,q,i) {f (N G , q , j) + r {j,i} }, ∀(N G, q, i) ∈ E ,
where functions lb i (q) are computed as described in Section 5.1.1 and the initialization f ({0}, 0, 0) = 0 and f ({0}, q, 0) = ∞, ∀ 0 < q ≤ Q is required. We define a ng- * route as a route obtained by adding, to an ng- * path (N G, q, i), edge e = {0, i}; the cost of an ng- * route is equal to the cost of ng- * path (N G, q, i) plus r e .

Procedure CG

Let R ⊆ R > be a subset of routes satisfying a given route relaxation. Moreover, given a route , we denote with q(R ) = i∈V C (R) q i + i∈V A (R ) q i its load. Procedure CG works as follows.

Step 1. Initialization. Generate a route set R to initialize the master problem which corresponds to LSP , where R is replaced with R. We assume that R contains at least one route containing each customer i ∈ V C . Set LCG = 0 and iter = 1.

Step 2. Find a master dual solution u of cost z. Initializes z = 0 and performs M axit2 iterations of the following operations.

(i) Compute a dual solution u of the master of cost z by means of expressions ( 30) and ( 31), where R is replaced with R and by using the current vector λ. Let R be the index set of routes producing ν i , i ∈ V C , in expressions (30), and let (i) be the index of the route in R associated with

ν i , i ∈ V C . Define a non-necessarily feasible solution ξ of LSP as ξ = i∈V C a i qi q(R ) ζ i , ∈ R, by setting ζ i (i) = 1 and ζ i = 0, ∀ ∈ R \ { (i)}, ∀i ∈ V C . If z > z, update z = z, ξ = ξ, u = u. (i) Update the penalty vectors λ as follows. Compute α i = ∈ R a i ξ , i ∈ V C , and α i = ∈ R a i ξ , i ∈ V F . Then, vector λ is modified as follows: λ i = λ i -γ(α i -1), i ∈ V C , and λ i = min{0, λ i -γ(α i -1)}, i ∈ V F .
where is a positive constant and γ = 0.2z

i∈V (αi-1) 2 . Step 3. Check if u is a feasible DSP solution. Generate the largest subset N ⊆ R > of routes
having negative reduced cost with respect to the current dual master solution u and such that |N | ≤ ∆ (∆ is an a priori defined parameter). If N = ∅ and z is greater than LCG,

then LCG = z, u * = u, ξ * = ξ and λ * = λ; otherwise, R = R ∪ N is updated.
Step 4. Termination criterion. Set iter = iter + 1. If iter = M axit1, stop.

Computing lower bound LB 1 Lower bound LB 1 corresponds to lower bound LCG computed by procedure CG using q- * route relaxation. The initial route set R of the master problem contains a feasible solution generated with the heuristic algorithm described in 6.1. We initialize λ = 0.

Define the modified routing cost r {i,j} = r {i,j} -(1/2)(u i +u j ), ∀{i, j} ∈ E (we assume u 0 = 0), and the modified assignment cost

d ij = d ij -u i , ∀(i, j) ∈ A, with respect to the current dual solution u.
The set N is computed as follows. We compute functions lb k (q), f (q, i) and g(q, i) using the modified routing and assignment costs r {i,j} and d ij instead r {i,j} and d ij . Let h(i) = min qi≤q≤Q {f (q, i)+r {0,i} }, if ∀i ∈ V C , and h(i) = min q min ≤q≤Q {f (q, i) + r {0,i} }, ∀i ∈ V F . The set N contains any q- * route corresponding to h(i) < 0, i ∈ V . Set u 1 = u * , λ 1 = λ * , and LB 1 = LCG.

Computing lower bound LB 2 Lower bound LB 2 corresponds to lower bound LCG computed by procedure CG using ng- * route relaxation.

We initialize λ

= λ 1 , define r 1 {i,j} = r {i,j} -(1/2)(u 1 i + u 1 j ), ∀{i, j} ∈ E (we assume u 1 0 = 0), d 1 ij = d ij -u 1 i , ∀(i, j) ∈ A,
and compute N i to be the Γ nearest nodes to i according to r 1 {i,j} . We compute functions f (N G, q, i) and lb k (q) using r 1 {i,j} and d 1 ij instead of r {i,j} and d ij , respectively, and the costs

h(i) = min (N G,q,i)∈E {f (N G, q, i) + r 1
{0,i} }, of the least cost ng- * route visiting i immediately before arriving at the depot. The initial route set R contains the ng- * routes corresponding to h(i) < 0, i ∈ V . At each iteration of procedure CG, to generate the set N , we compute functions f (N G, q, i) and lb k (q) with the modified routing cost r {i,j} = r {i,j} -(1/2)(u i + u j ), ∀{i, j} ∈ E, and the modified assignment cost

d ij = d ij -u i , ∀(i, j) ∈ A, with respect to the current solution u. N contains every ng- * route corresponding to h(i) = min (N G,q,i)∈E {f (N G, q, i) + r {0,i} } < 0, i ∈ V . Set LB 2 = LCG.

Column-and-cut generation method

In this section, we describe a bounding procedure that computes a lower bound on the vrptf as the cost of an optimal solution of problem LSP obtained from formulation LSP by substituting the route set R with the set R > of ng- * route and by adding valid inequalities derived from a family F of valid inequalities described for formulation T I.

Any valid inequality t ∈ F can be expressed in general form as

e∈E α t e x e + i∈V β t i y i + (i,j)∈A γ t ij z ij ≥ ω t , ( 45
)
and can be transformed into the following valid inequality for formulation SP using equations ( 32)-( 35), where R is substituted by R > :

∈R > (ϕ t + φ t + ψ t )ξ ≥ ω t , ( 46
)
where

ϕ t = e∈E α t e η e , φ t = i∈V β t i a i and, ψ t = (i,j)∈A γ t ij b j i .
The bounding procedure solves problem LSP by using column and cut generation. The initial master problem is obtained from the computation of lower bound LB 2 by replacing the route set R > with the route set R generated by procedure CG during the computation of LB 2 . The initial set of valid inequalities F is set to the empty set. At each iteration (say k), the procedure performs the following steps.

1. Solve problem LSP . Let ξ and (u, v) be the optimal primal and dual solutions, respectively.

Vector u is given by u 3. Given the solution vector ξ, compute the corresponding solution vector (x, z, w) by means of equations ( 32)-( 35) where R is substituted by R. Solve the separation problems associated with the set of valid inequalities F (see below) and add, if any, violated inequalities to set F .

= {u 0 , u 1 , . . . , u |V C | , u |V C |+1 , . . . , u |V | },
It can be easily shown that the complexity of the pricing algorithm solved at Step 2 of the above procedure is not sensitive to the addition of the valid inequalities in F , since the values of the corresponding dual variables can be translated into subproblem costs. Indeed, at each iteration of the procedure, to generate the set N , we compute the ng- * route functions f (N G, q, i) and lb k (q) with the modified routing cost

r {i,j} = r {i,j} -(1/2)(u i + t∈F β t i v t )-(1/2)(u j + t∈F β t j v t )-t∈F α t {i,j} v t , ∀{i, j} ∈ E, and the modified assignment cost d ij = d ij -u i -t∈F γ t ij v t , ∀(i, j) ∈ A, with respect to the current dual solution (u, v). N contains every ng- * route corresponding to h(i) < 0, i ∈ V .
We conducted preliminary experiments to identify a good separation strategy to be used at Step 3. As a result of our experimentation, we decided to use the following inequalities to define the family set F : CI, MI, RCI, RCII-a, RCII-b, RCII-c, and RCII-d inequalities. For a given solution (x, z, w), we identified (as far as possible) violated inequalities of above seven types by applying the corresponding separation procedures as described below.

Separation procedures

The separation problems of CI, RCII-a and RCII-c inequalities can be reduced to max-flow/min-cut problems using a standard construction, and therefore solved in polynomial time; we omit the details for sake of brevity (see [START_REF] Baldacci | The Capacitated m-Ring-Star Problem[END_REF]). Concerning MI inequalities, the following theorem holds.

Theorem 4 Let (x, z, y) be a solution of the LP-relaxation of formulation T I and assume that q i ≤ Q, ∀i ∈ V C , and that x e = 0, e = {i, j} ∈ E \ {{0, h} : h ∈ V }, if q i + q j > Q. The separation problem for MI inequalities (10) is solvable in polynomial time.

Proof. The proof is provided in the e-companion to this paper. RCI, RCII-b and RCII-d inequalities are separated using a heuristic separation procedure. The procedure is a Multistart Local Search that, at each iteration, generates a starting point and evolves it through a Local Search procedure. We start by generating a set S of 10(n -1) subsets of V as follows. For the RCI inequalities the first |V C | sets of S are obtained by inserting in each set, for i = 1, . . . , |V C |, the nodes in F i . The remaining sets are generated by first computing a random number m drawn from a uniform distribution in [1, . . . , n -1], and then by randomly selecting m different nodes of V , again using a uniform distribution. For the RCII-b and RCII-c inequalities all the sets are randomly generated as above. Each set S ∈ S is then iteratively expanded by adding one node at each iteration until S = V . For a given set S, let θ(S) denote the difference between the left-hand side and the right-hand side value of the considered inequality (i.e., the inequality can be rewritten as θ(S) ≥ 0 and the separation problem corresponds to compute arg min S⊆V {θ(S)}). Each set S is expanded by choosing the node i ∈ V \ S such that θ(S ∪ {i}) is minimized

Solving the vrptf to Optimality

In this section, we describe the method implemented for solving the vrptf to optimality. We start by describing two heuristic algorithms that compute primal bounds used to initialize the exact method. The exact method is a branch-and-cut-and-price (BCP) solution method based on the SCIP (see [START_REF] Achterberg | Scip: Solving constraint integer programs[END_REF]) BCP solution framework.

Heuristic algorithms

Primal bounds for the vrptf are computed by means of two different types of heuristic algorithms: a constructive heuristic and a Lagrangean heuristic.

The basis of the constructive algorithm is a heuristic to solve the cvrp. Given an instance of vrptf, we define a complete graph G = (V , E) where the node set V = {0} ∪ V C contains the depot and the customer nodes. Each edge e ∈ E has a cost given by r e . Each customer i ∈ V C has a demand equal to q i and the capacity of the vehicles is set to Q. Roughly speaking, we solve a problem obtained from vrptf by disregarding the facility nodes (set V F ) and the connection arcs (set A). The cvrp instance is solved through an iterative multistart procedure based on a cluster-first, route-second heuristic procedure. Each iteration consists of three phases: (i) determine a partition of the customers into a number of subsets each one satisfying the capacity constraint; (ii) for each set, find the route of a single vehicle that serves all the customers in the set (i.e. we solve an instance of a Traveling Salesman Problem (tsp)); (iii) locally optimize the solution obtained at step (ii). The cvrp solution so far obtained, is then locally optimized by iteratively applying two post-optimization procedures specifically devised for the vrptf.

The Lagrangean heuristic is based on procedure CG described in Section 5.1.3. Procedure CG is interwoven with an algorithm that produces a feasible vrptf solution using the route set R (see Step 2 of procedure CG). The route set R is first modified to contain only customers visited at most once. Then, unrouted customers are inserted in order to obtain a feasible solution. The solution obtained is further optimized by applying the same post-optimization procedures used by the constructive algorithm.

A step-by-step description of the heuristics are given in the e-companion to this paper.

Details of the BCP method

The lower bound at the root node of the enumeration tree is first computed by using the bounding procedure described in Section 5, then by using the column-and-cut generation method described in Section 5.2. The master problem at a generic node except the root node is initialized with the set of valid inequalities F and the set of routes R of the parent node, where set R is further modified by extracting the largest set of routes satisfying the branching conditions.

To choose a node-selection rule, we first performed some preliminary experiments with different rules and, based on these results, we decided to adopt the best-first strategy for all the computations of Section 7. We did not implement primal heuristics but the algorithm was initialized with the best primal solution found by the two heuristic algorithms described in the previous section that are executed at the root node. We used the default branching scheme of the SCIP framework, namely the hybrid branching scheme (see [START_REF] Achterberg | Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems: 6th International Conference[END_REF], that combines ideas from pseudocost branching [START_REF] Benichou | Experiments in mixed-integer linear programming[END_REF]) and strong branching [START_REF] Applegate | The Traveling Salesman Problem: A Computational Study[END_REF].

Computational Results

This section reports on the computational results of the exact method described in this paper and analyses the effectiveness of the dual ascent heuristics and of the different types of inequalities on the bounding procedure procedure described in Section 5.

The algorithms were coded in C++ and linked with the SCIP 3.1.1 BCP solution framework (see Achterberg 2009) using the IBM Cplex 12.6.1 linear programming solver (see IBM CPLEX 2014). The experiments were performed on an Intel Core 2 Duo at 2.66 GHz personal computer equipped with 4 Gb of RAM.

The exact method has been tested on real-world instances and on instances derived from lrp instances already proposed in the literature, used to further evaluate the performance of our algorithms. The same instances have been also used to generate 2e-cvrp instances. The following sections 7.1 and 7.2 briefly describe the real-world and lrp based instances, respectively, and report on the results obtained by the different algorithms. The complete details of the instances are provided in the e-companion to this paper.

Based on the results of preliminary experiments to identify good parameter settings for our method, we decided to use the following settings for our bounding procedure (see Section 5):

• in computing lower bound LB 1 : M axit1 = 50, M axit2 = 50, = 1.5 and ∆ = 50;

• in computing lower bound LB 2 : Γ = 12, M axit1 = 100, M axit2 = 50, = 2.0 and ∆ = 50;

• in the column-and-cut-generation method: ∆ = 100 at the root node of the BCP whereas ∆ = 50 for the remaining nodes.

Results on real-world instances

The data of this set of instances were provided by a major Italian transportation company that distributes non-perishable products over the whole Italian peninsula. The company operates through three main distribution areas (North, Centre and South) using three main central depots located in the provinces of Milan, Rome and Naples.

The three distribution areas operate independently in the corresponding areas to serve customer orders using an existing set of intermediate facilities. The customer orders are placed into Europallet and distributed either to the final customers or the intermediate facilities by means of a fleet of identical capacitated vehicles which are stationed at the different central depots and whose capacity is expressed in terms of pallets. All the facilities are owned by third-party contractors, that are in charge of delivering to the final customers the orders consolidated at the facilities.

The company was interested in analyzing different distribution scenarios associated with the three distribution areas. A total number of 18 instances were provided by the company, six instances per each area or depot. The following naming convention was adopted to identify the different instances. The instance name is a string area a×b Qc, where area represents the area (i.e., North, Centre, South), a represents the number of customers, b corresponds to the number of facilities, and c is the vehicle capacity.

In Table 1, we report the results obtained by the heuristic algorithms, the bounding procedure and the BCP method. The columns of the table report the instance name (Name), the cost of the best solution found by the heuristics and BCP algorithms (z * ), the percentage deviation of the upper bound computed by the constructive heuristic (%U B 1 ), the percentage deviation of the upper bound computed by the lagrangean heuristic (%U B 2 ), the percentage deviation of lower bound LB 1 (%LB 1 ), the percentage deviation of lower bound LB 2 (%LB 2 ), the total computing time of lower bounds LB 1 and LB 2 that also includes the time spent for computing U B 2 (t DA ), the percentage deviation of the lower bound LB computed at the root-node of the BCP algorithm and the corresponding computing time (%LB, t LB ), the cardinality of the sets F and R associated with lower bound LB (#cuts and #cols), the total number of nodes of the exact algorithms (#N), the percentage deviation of the best lower bound achieved by the exact method (%Opt), and the total computing time in seconds spent by the exact method (t T OT ), that also include the time spent for computing upper bound U B 1 . The percentage deviation of value x is computed as 100 × x/z * . In order to evaluate the quality of the different lower bounds, we also computed, for each instance, the value of the lower bound obtained by solving the LP-relaxation of formulation T I strengthened with the different valid inequalities (using the separation strategy described in Section 5.2). In the table, column %LB C reports the percentage deviation of the final lower bound obtained whereas column t C displays the corresponding computing time.

Table 1: Results on real-world instances

Name z * #r #f #c %U B 1 %U B 2 %LB 1 %LB 2 t DA %LB C t C %LB t LB #cuts #cols #N %Opt t
For each instance, Table 1 also reports the following details about the best solution found: the number of routes in the solution (#r), the number of facilities visited (#f ) and the number of customers assigned to a facility (#c).

For these set of instances, a time limit of 7,200 seconds was imposed to the SCIP framework.

The last row of the table reports averages computed over the different columns. The average reported under column t T OT is computed over the instances solved to optimality within the imposed time limit. If a value of 100.0 is reported for column %Opt, then the algorithm terminated with an optimal solution.

Table 1 shows that 8 out of 18 instances were solved to optimality and that the final lower bound LB is on average quite tight, being equal to 98.7%. The largest instance solved to optimality involves 142 customers and 18 facilities. On these set of instances, lower bounds LB 1 and LB 2 have the same quality and are on average superior to lower bound LB CP , thus showing the effectiveness of our q- * route and ng- * route relaxations. Moreover, the different valid inequalities can substantially increases the lower bound, as shown by the improvements on instances north-68x7-Q24 and south-54x4-Q34.

The table shows that upper bound U B 2 is always better than upper bound U B 1 and that the BCP algorithm can further improve the upper bounds in almost all instances, thus producing high quality primal solutions also whenever the algorithm terminates without proving the optimality of the solution found.

It is worth mentioning that the time spent for computing upper bound U B 1 is on average equal to 187.4 seconds and that the time spent by the procedure used to compute upper bound U B 2 (called during the computation of lower bound LB 2 ) is on average equal to 226.8 seconds. Therefore, both U B 1 and U B 2 can be computed efficiently in practice.

Results on lrp based instances

This set of instances was derived from 75 lrp instances used in [START_REF] Baldacci | An Exact Method for the Capacitated Location-Routing Problem[END_REF] and [START_REF] Contardo | An Exact Algorithm Based on Cut-and-Column Generation for the Capacitated Location-Routing Problem[END_REF] for solving the lrp and proposed by different authors. We derived two classes of test instances (A and B) having the same topology of the underlying graph, but with different cost structures.

We generated a total number of 150 instances, 75 instances per class. The dimensions of the instances vary from very small instances with 12 customers and two facilities up to large instances with 150 customers and 20 facilities. The instance name is a string name<a×b>, where name represents the instance name, a represents the number of customers and b corresponds to the number of facilities.

For sake of presentation, the instances were grouped into the following three groups accordingly to the original lrp source: ii) [START_REF] Prins | Nouveaux algorithmes pour le problème de localisation et routage sous contraintes de capacité[END_REF]: 24 instances involving 20, 50, and 100 customers, 5 or 10 facilities;

iii) Different authors: 39 instances, involving up to 150 customers and 20 facilities.

For this set of instances, a time limit of 3,600 seconds was imposed to the SCIP framework.

Tables 2 and3 summarize the results obtained on both classes A and B. In the tables, column #Opt reports for each group of instances the total number of instances solved to optimality within the imposed time limit.

The meaning of the remaining columns is the same described in the previous section, but in the tables their values are relative to averages computed over the instances composing the three groups. The values reported under column t T OT are computed over the instances solved to optimality within the imposed time limit.

Tables 2 and3 show that 30 and 24 out of 75 instances were solved to optimality within the imposed time limit for classes A and B, respectively. For these instances, lower bound LB 2 is on average superior with respect lower bound LB 1 . As the feasible solutions associated with these instances are characterized (on average) by a larger number of customers per route, the ng- * route relaxation performs in practice better than q- * route relaxation. Also for these instances, the different valid inequalities can substantially increase the final lower bound (see column %LB). Instances of Class B are more difficult with respect to the corresponding instances of class A. This is due to the different cost structure of class B instances and it is testified by the worse quality of lower bounds LB CP and of the final lower bound LB. Nonetheless, lower bounds LB 1 and LB 2 show the same quality of class A instances.

Concerning the upper bounds, the tables show that both the two upper bounding procedures can compute good quality solutions. The average computing time of upper bound U B 1 (U B 2 ) is equal to 70.8 and 72.9 seconds (148.0 and 89.5 seconds) for classes A and B, respectively. Therefore, the computation of LB 2 requires a higher computing time with respect to the real-world instances and this is due to the larger vehicle capacity that characterizes most of the instances in classes A and B.

The detailed results reported in the e-companion to this paper show that instances with up to 100

Effectiveness of the dual ascent heuristics and valid inequalities

Table 4 reports an analysis of the effectiveness of the dual ascent heuristics when used to initialize the master problem of problem LSP (see Section 5.2). In order to assess the quality of lower bounds LB 1 and LB 2 , we solved problem LSP without adding valid inequalities, i.e., we computed the optimal solution cost LB SP of formulation LSP and the LP-relaxation of formulation SP with ng- * route.

In addition, the Lagrangean heuristic has been disabled during the computation of LB 1 and LB 2 .

The table reports the average percentage deviations of lower bounds LB 1 , LB 2 , and LB SP under columns %LB 1 , %LB 2 and %LB SP , respectively. The table then reports, under heading t LB SP , the average total computing times spent in computing lower bound LB SP under the following options: (a) without computing lower bounds LB 1 and LB 2 (b) by computing lower bound LB 1 (c) by computing lower bound LB 2 , and (d) by computing both lower bounds LB 1 and LB 2 . In case (a), the master problem of LSP is initialized with single-customer routes whereas in case (b), the master problem is initialized using the dual solution provided by lower bound LB 1 , that is used to generate an initial set of ng- * route. In cases (c) and (d), the master problem is initialized with the route set generated by procedure CG during the computation of LB 2 (as described in 5.2). Moreover, in case (c) the master problem associated with the computation of LB 2 , is initialized as for LB 1 , i.e., using the solution provided by the constructive heuristic described in Section 6.1.

All values in the table are relative to averages computed over the instances composing the three groups of classes A and B, and over the real-world instances. The last row of the table reports averages computed over all instances.

The table shows that the bounding procedure based on the use of both lower bounds LB 1 and LB 2 (case (d)) is about five times faster than the standard column generation method (case (a)). Generally speaking, standard column generation methods are time-consuming as the LP-relaxation of the master problem is usually highly degenerate and degeneracy implies alternative optimal dual solutions. Consequently, the generation of new columns and their associated variables may not change the value of the objective function of the master problem, the master problem may become large, and the overall method may become slow computationally. In case (d), the bounding procedure starts from a near-optimal dual solution of the LP-relaxation of SP with ng- * route provided by lower bound LB 2 , as shown by the percentage deviations of lower bounds LB 2 and LB SP . This allows us to generate an initial master problem containing the routes having a very small reduced cost that are likely to be in the optimal LSP solution.

The analysis of cases (b) and (c) shows that it is also computationally convenient to compute LB 1 or LB 2 . In particular, computing LB 1 before the computation of LB 2 speedup the computation of LB 2 as procedure CG used to compute LB 2 takes advantage from the master initialization provided by the dual solution corresponding to LB 1 . Table 4 also reports the computational results obtained when calculating the lower bound LB 1 under the following two ways of initializing the corresponding master problem: (i) by using the heuristic solution provided by the constructive heuristic (case (e)) (ii) by using single-customer routes (case (f)). The table shows that on average, the difference is slightly marginal. Nevertheless, as in our implementation the constructive heuristic is executed before computing LB 1 , it is worthwhile to initialize the master of LB 1 with the solution found by the heuristic.

Table 5 analyses the impact of the valid inequalities on the column-and-cut bounding procedure described in Section 5.2 at the root node of the BCP method.

The table reports average percentage deviations of the lower bounds obtained by the bounding procedure under the following cases: (i) without adding valid inequalities (under column heading "no cuts") (ii) by adding CI, MI and RCI inequalities ("+ CI + MI + RCI") (iii) by adding CI, MI, RCI, RCII-a, and RCII-b inequalities ("+ RCII-a + RCII-b"), and (iv) by adding CI, MI, RCI, RCII-a, RCII-b, RCII-c, RCII-d inequalities ("+ RCII-c + RCII-d"). The last case corresponds to the final procedure we adopted in our computational results and, as mentioned in Section 5.2, the sequence of separation procedures was defined after conducting preliminary computational experiments performed to identify a good separation strategy.

For each group, the table reports the average percentage deviations of the lower bounds obtained and the corresponding average computing times (%LB, t LB ), and the average cardinalities of the sets F associated with the lower bound computation (#cuts). As for Table 4, the Lagrangean heuristic has been disabled during the computation of LB 1 and LB 2 . In addition, the time t LB also includes the time spent for computing LB 1 and LB 2 .

As for Table 4, all values in the table are relative to averages computed over the instances composing the three groups of classes A and B, and over the real-world instances. The last row of the table reports averages computed over all instances.

The table shows that the average percentage gaps left by considering in turn the different three groups of valid inequalities are equal to 3.5, 2.9 and 2.8, respectively. With respect to the "no cuts" case, a final gap reduction of about 1.4% has been achieved. The contribution given by inequalities RCII-c and RCII-d is on average equal to 0.1% as shown by the table. During preliminary computational experiments, we observed that their addition generally results in separating additional RCI and RCII-b inequalities, which separation procedures are heuristics.

Conclusions

In this paper, we considered a vehicle routing problem with transhipment facilities, called the Vehicle Routing Problem with Transhipment Facilities (vrptf), that was motivated by a real-world application of interest to an Italian company operating in the production and distribution of nonperishable products. The vrptf consists of selecting transhipment facilities, allocating customers to these facilities and designing vehicle routes emanating from a central depot to minimize the total distribution cost. A feature of the problem is that a customer can be either served on a vehicle route emanating from the central depot or through an intermediate facility, where the demand is first delivered by a vehicle route, and then it is successively delivered to the final customer.

We proposed two integer programming formulations for the vrptf, a two-index formulation (T I) and a set-partitioning based formulation (SP ). The formulations were used to derive a bounding method based on two dual ascent heuristics and a column-and-cut generation procedure. In particular, we proposed valid inequalities to strengthen the linear relaxations of the two formulations and two different route relaxations, called q- * route and ng- * path, that have the advantage that the pricing subproblem associated with the linear relaxation of formulation SP can be efficiently solved (by dynamic programming).

All our findings have been used to develop branch-and-cut-and-price algorithm that has been tested on a large family of instances, including both real-world instances and instances derived from the literature.

The implementation solved to optimality different instances from our real-world instances involving up to 142 customers and 18 facilities. The implementation was also tested on literature-based instances to better evaluate the limits of the algorithms, and the new approaches can find optimal solutions on some difficult instances with up to 100 customers and 10 facilities.

e-companion 9 Proofs of statements

Proposition 1 There are no dominance relations between inequalities RCII-a and RCII-b.

Proof. Consider a vrptf instance with |V

C | = 6, |V F | = 1, with V C = {1, 2, 3, 4, 5, 6} and V F = {7}.
In addition, let Q = 10 and q 1 = 4, q 2 = 3, q 3 = 5, q 4 = q 5 = q 6 = 2. First select a set S = {1, 2, 4, 7} and such that the nodes {1, 2, 4, 7} are visited on a route while customers {3, 5, 6} are assigned to facility node 7 of S. The right-hand side of inequality ( 15) becomes 18/10 -0 = 2, while the right-hand side of ( 17) has value 9/10 -0 = 1 and ( 15) is stronger than ( 17). Now consider a set S = {1, 2, 3, 4, 5} containing five customers all visited on a route. Customer 6 is associated with a node in V \S. The right-hand side of (15) takes value 18/10 -1/4 = 7/4, while the right-hand side of ( 17) becomes 16/10 -0 = 2 and the second inequality is stronger than the first one. Moreover, note that the FrCC inequalities (4) are dominated by ( 15), for the first example, and by both ( 15) and ( 17) for the second one.

Lemma 3 (Lemma 2.) Let o ∈ R with ô > 0 and T = {m ∈ R, n ∈ Z : m + n ≥ o, m ≥ 0}. The inequality m + ôn ≥ ô o (18)
is valid for T .

Proof. We have two cases:

(i) n ≥ o . As m ≥ 0, we have m ô ≥ 0, hence m ô + n ≥ o ; (ii) n ≤ o . As 0 < ô < 1 we have that o -n ≥ ô( o -n). ( 47
)
Since o = o + ô and using inequality (47), inequality m + n ≥ o can be rewritten as:

m ≥ ô + ô( o -n). ( 48
)
The right-hand side of inequality ( 48), can be rewritten as:

ô(1 + o ) -ôn = ô o -ôn (49)
thus obtaining m ≥ ô o -ôn.

Theorem 5 Let α e ≥ 0, ∀e ∈ E, β i ≥ 0, ∀i ∈ V and γ ij ≥ 0, ∀(i, j) ∈ A and consider the following inequality valid for formulation T I:

e∈E α e x e + i∈V β i y i + (i,j)∈A γ ij z ij ≥ o ( 19
)
where o ∈ R and ô > 0. Then the following inequality:

e∈E ϕ o (α e )x e + i∈V ϕ o (β i )y i + (i,j)∈A ϕ o (γ ij )z ij ≥ o (20)
where 19) we can round up the coefficients in E 2 , V 2 and A 2 to obtain:

ϕ o (m) = m + min m ô , 1 , m ∈ R, o ∈ R, ô > 0, is also a valid inequality for formulation T I. Proof. Let E 1 ⊆ E, E 2 = E \ E 1 , V 1 ⊆ V , V 2 = V \ V 1 and A 1 ⊆ A, A 2 = A \ A 1 . Starting from inequality (
e∈E 1 α e x e + i∈V 1 β i y i + (i,j)∈A 1 γ ij z ij + e∈E 2 α e x e + i∈V 2 β i y i + (i,j)∈A 2 γ ij z ij ≥ o. ( 50
)
Writing

α e = α e + αe , ∀e ∈ E 1 , β i = β i + βi , ∀i ∈ V 1 , and γ ij = γ ij + γij , ∀(i, j) ∈ A 1
and re-arranging terms, we get:

(

e∈E 1 αe x e + i∈V 1 βi y i + (i,j)∈A 1 γij z ij )+ ( e∈E 1 α e x e + e∈E 2 α e x e + i∈V 1 β i y i + i∈V 2 β i y i + (i,j)∈A 1 γ ij z ij + (i,j)∈A 2 γ ij z ij ) ≥ o. ( 51
)
The first part of inequality ( 51) is non-negative, and the second part is integral for all x, w and z integral. Applying Lemma 2 we get:

1 ô ( e∈E 1 αe x e + i∈V 1 βi y i + (i,j)∈A 1 γij z ij )+ ( e∈E 1 α e x e + e∈E 2 α e x e + i∈V 1 β i y i + i∈V 2 β i y i + (i,j)∈A 1 γ ij z ij + (i,j)∈A 2 γ ij z ij ) ≥ o . ( 52
)
The coefficients of variables {x e } in (52) are α e + αe ô if e ∈ E 1 and α e if e ∈ E 2 . Similarly, the coefficients of variables {y i } in (52) are

β i + βi ô if i ∈ V 1 and β i if i ∈ V 2 and the coefficients of variables {z ij } are γ ij + γij ô if (i, j) ∈ A 1 and γ ij if (i, j) ∈ A 2 .
The best choices of coefficients for the sets

E 1 , V 1 and A 1 are E 1 = {e ∈ E : αe ≤ ô}, V 1 = {i ∈ V : βi ≤ ô} and A 1 = {(i, j) ∈ A : γij ≤ ô}, respectively. By defining ϕ o (m) = m + min m ô , 1 , m ∈ R, o ∈ R, ô > 0, inequality (52)
becomes inequality (20).

Theorem 6 Let us associate penalties λ

i ∈ R, ∀i ∈ V C , with constraints (27), and λ i ≤ 0, ∀i ∈ V F , with constraint (28). For each i ∈ V C , define a i = a i + j∈V F (R ) b j i , and let R i = { ∈ R : a i > 0}. For each i ∈ V C compute: ν i = q i min ∈Ri (c + p ) -j∈V C a j λ j -j∈V F a j λ j j∈V C a j q j (30)
The VRPTF R. Baldacci, S. U. Ngueveu, R. Wolfler Calvo A feasible DSP solution u of cost z(DSP (λ)) is given by the following expressions:

u i = ν i + λ i , ∀i ∈ V C , and u i = λ i , ∀i ∈ V F . (31) Proof. Consider a route ∈ R. Since ∈ R i , ∀i ∈ V C (R )
, from expression (30) we derive:

ν i ≤ q i (c + p ) -j∈V C a j λ j -j∈V F a j λ j j∈V C a j q j , ∀i ∈ V C (R ). ( 53
)
Given a route ∈ R, from expression (31) we obtain:

i∈V C a i u i + i∈V F a i u i ≤ i∈V C a i q i (c + p ) -j∈V C a j λ j -j∈V F a j λ j j∈V C a j q j + i∈V C a i λ i + i∈V F a i λ i . ( 54
)
Inequality ( 54) can be written as:

i∈V C a i u i + i∈V F a i u i ≤ c + p , (55)
that corresponds to the constraint of problem DSP for route .

Let E(S) denote the set of edges in G with both end-nodes in S and, given two disjoint vertex sets S 1 , S 2 , let E(S1 : S2) denote the set of edges crossing from S 1 to S 2 (i.e., E(S 1 : S 2 ) = δ(S 1 ) ∩ δ(S 2 )) (if S 1 = {i}, we simply write E(i : S 2 ) instead of E({i} : S 2 )).

Theorem 7

The LP-relaxation of the SP formulation satisfies both CI and FrCI inequalities, and a weak form of MI inequalities.

Proof. Consider a set S ⊆ V with V C (S) = ∅ and let T = V C (S) be the set of customers contained in S. Define the surrogate constraint obtained by adding partitioning constraints (27) corresponding to customers in T after having multiplied the equation associated with i ∈ T by

q i : ∈R q (T )ξ = q(T ), (56)
where q(T ) = i∈S q i and q (T ) = i∈S q i a i . Since q (T ) ≤ min[Q, q(T )], we have

∈R(T ) ξ ≥ max[1, q(T )/Q], (57)
where R(T ) = { ∈ R : a i = 1 for some i ∈ T }. Given a route ∈ R(T ), define q (T ) as the total demand of the customers not in T assigned to the route, i.e. q (T ) = q (T ), where

T = (V C \ T ) ∩ V C (R ).
As q (T ) + q (T ) ≤ Q we have:

∈R(T ) Qξ ≥ ∈R q (T )ξ + ∈R q (T )ξ . ( 58
)
The VRPTF R. Baldacci, S. U. Ngueveu, R. Wolfler Calvo

From equations (56) and inequalities (58) we derive:

∈R(T ) ξ ≥ max    1, 1 Q   q(T ) + ∈R q (T )ξ      . ( 59
)
Note that any route ∈ R(T ) contains at least two edges, one having an ending node in S and the other in S. Therefore, we have:

{i,j}∈δ(S)

η ij ξ ≥ 2ξ . ( 60
)
Adding inequality (60) for all ∈ R(T ) we obtain:

∈R(T ) {i,j}∈δ(S) η ij ξ ≥ 2 ∈R(T ) ξ . ( 61
)
Thus inequalities (59) become

∈R(T ) ρ (S)ξ ≥ 2 max    1, 1 Q   q(T ) + ∈R q (T )ξ      , ( 62
)
where ρ (S) = {i,j}∈δ(S) η ij . Since

q (T ) ≥ j∈T k∈V F (S) b k j q j + 1 2 j∈T q j {i,h}∈E(S:{j}) η ih (63)
from ( 63) we obtain:

∈R(T ) ρ (S)ξ ≥ 2 max{1, 1 Q (q(T ) + ∈R ( j∈T k∈V F (S) b k j q j + 1 2 j∈T q j {i,h}∈E(S:{j}) η ih )ξ )}. ( 64
)
We have:

i) q(T ) = ∈R q (T )ξ = i∈T q i ( ∈R a i ξ ) + i∈T q i ( ∈R k∈V F (R ) b k i ξ ); ii) ∈R j∈T k∈V F (S) b k j q j ξ = j∈T q j ( ∈R k∈V F (S) b k j ξ ).
Using the equations ( 32)-( 35) linking variables variables ξ with (x, z, w), we derive:

i) i∈T q i ( ∈R a i ξ ) = i∈V C (S) q i y i ; ii) i∈T q i ( ∈R k∈V F (R ) b k i ξ ) + j∈T q j ( ∈R k∈V F (S) b k j ξ ) ≥ (i,j)∈A: j∈V F (S) q i z ij ; iii) 1 2 j∈T q j ( ∈R ) {i,h}∈E(S:{j}) η ih ξ = 1 2 i∈V C (S)
j∈S q i x {i,j} .

Using the above equations, from (64) we obtain:

e∈δ(S)
x e ≥ 2 max

   1, 1 Q   i∈V C (S)
q i y i + (i,j)∈A:j∈V F (S)

q i z ij + 1 2 i∈V C (S) j∈S q i x {i,j}      (65)
The VRPTF R. Baldacci, S. U. Ngueveu, R. Wolfler Calvo Theorem 8 Let (x, z, y) be a solution of the LP-relaxation of formulation T I and assume that q i ≤ Q, ∀i ∈ V C , and that x e = 0, e = {i, j} ∈ E \ {{0, h} : h ∈ V }, if q i + q j > Q. The separation problem for MI inequalities (10) is solvable in polynomial time.

Proof. Consider the MI inequality for a given set S ⊆ V , S = ∅:

e∈δ(S) x e ≥ 2 Q   i∈V C (S)
q i y i + (i,j)∈A:j∈V F (S)

q i z ij + i∈V C (S) j∈S q i x {i,j}   . ( 66
)
We have:

e∈δ(S)
x e = e∈E(0:S)

x e + e∈E(S:S)

x e , (67)
and for each i ∈ V (see equation ( 2)):

e∈δ(i) x e = 2y i = x {0,i} + e∈E(i:S) x e + e∈E(i:S) x e . ( 68
)
From equation ( 68), the term i∈V C (S) q i y i of inequality (66) can be rewritten as follows:

i∈V C (S)

q i y i = i∈V C (S) q i 2   x {0,i} + e∈E(i:S) x e + e∈E(i:S) x e   . ( 69
)
The MI inequality (66) can be rewritten as:

e∈E(0:S)
x e + e∈E(S:S)

x e ≥ i∈V C (S)

q i Q x {0,i} + i∈V C (S) e∈E(i:S) q i Q x e +
i∈V C (S) e∈E(i:S)

q i Q x e + 2 Q (i,j)∈A:
j∈V F (S)

q i z ij + i∈V C (S) j∈S q i Q x {i,j} + i∈V C (S) j∈S q i Q x {i,j} . ( 70
)
We also have:

e∈E(0:S)
x e = e∈E(0:V C (S))

x e + e∈E(0:S\V C (S))

x e , (

and e∈E(S:S)

x e = j∈S {i,j}∈E(j:V C (S))

x {i,j} + j∈S {i,j}∈E(j:(V \V C (S))\S)

x {i,j} .

(72)

Notice that S \ V C (S) = V F (S) and that (V \ V C (S)) \ S = V F (S)
. Then, inequality (70) can be rewritten as:

e∈E(0:V C (S))
x e + e∈E(0:V F (S))

x e + j∈S {i,j}∈E(j:V C (S))

x {i,j} + j∈S {i,j}∈E(j:V F (S))

x {i,j} ≥ i∈V C (S)

q i Q x {0,i} + i∈V C (S) e∈E(i:S) q i Q x e + i∈V C (S) e∈E(i:S) q i Q x e + 2 Q (i,j)∈A:
j∈V F (S)

q i z ij + i∈V C (S) j∈S q i Q x {i,j} + i∈V C (S) j∈S q i Q x {i,j} . (73)
Notice that as q i = 0, ∀i ∈ V F , we have:

i∈V C (S) e∈E(i:S) q i Q x e =
i∈S e∈E(i:S)

q i Q x e = i∈S e∈E(i:V C (S)) q i Q x e + i∈S e∈E(i:V F (S)) q i Q x e , (74)
The VRPTF R. Baldacci, S. U. Ngueveu, R. Wolfler Calvo i∈V C (S) e∈E(i:S)

q i Q x e = i∈S e∈E(i:S) q i Q x e , (75)
and

i∈V C (S) j∈S q i Q x {i,j} = i∈S j∈S q i Q x {i,j} . (76)
Inequality ( 73) can be rewritten as:

{0,i}∈E(0:V C (S)) (1 -q i /Q)x {0,i} + e∈E(0:V F (S))
x e + j∈S {i,j}∈E(j:V C (S))

(1 -

(q i + q j )/Q)x {i,j} + j∈S {i,j}∈E(j:V F (S)) (1 -q j /Q)x {i,j} ≥ i∈S e∈E(i:S) q i Q x e + 2 Q (i,j)∈A:
j∈V F (S)

q i z ij + i∈S j∈S q i Q x {i,j} . (77)
Since i∈S e∈E(i:S)

q i Q x e + i∈S j∈S q i Q x {i,j} = j∈S     {j,i}∈δ(j): j<i q i Q x {j,i} + {i,j}∈δ(j): i<j q i Q x {i,j}     = j∈V     {j,i}∈δ(j): j<i q i Q x {j,i} + {i,j}∈δ(j): i<j q i Q x {i,j}     - j∈S     {j,i}∈δ(j): j<i q Q x {j,i} + {i,j}∈δ(j): i<j q i Q x {i,j}     , (78) and
(i,j)∈A: j∈V F (S)

q i z ij = q(V C ) - i∈V C (S) q i y i - i∈V C (S)
q i y i -(i,j)∈A: j∈V F (S)

q i z ij = q(V C ) - i∈V C q i y i - (i,j)∈A: j∈V F (S) q i z ij , (79)
inequality (77) can be rewritten as:

{0,i}∈E(0:V C (S)) (1 -q i /Q)x {0,i} + e∈E(0:V F (S))
x e + j∈S {i,j}∈E(j:V C (S))

(1 -

(q i + q j )/Q)x {i,j} + j∈S     {j,i}∈δ(j): j<i q i Q x {j,i} + {i,j}∈δ(j): i<j q i Q x {i,j}     + j∈S {i,j}∈E({j}:V F (S)) (1 -q j /Q)x {i,j} + 2 Q (i,j)∈A: j∈V F (S) q i z ij ≥ j∈V     {j,i}∈δ(j): j<i q i Q x {j,i} + {i,j}∈δ(j): i<j q i Q x {i,j}     + 2 Q (q(V C ) - i∈V C q i y i ).
(80) Notice that, as

q i ≤ Q, ∀i ∈ V C , x e = 0, e = {i, j} ∈ E \ {{0, h} : h ∈ V }, if q i + q j > Q, all
the variable coefficients of the above inequality are nonnegative whereas the right-hand-side of the inequality does not depend on the set S.

The most violated constraint (80) can now be found by computing a minimum s-t cut on an directed capacitated graph G = (V , A) with V = V ∪ {s, t} and A = {(i, j), (j, i) : ∀{i, j} ∈ E \ {{0, j} : j ∈ V }} ∪ {(s, i) : {0, i} ∈ E} ∪ {(i, t) : i ∈ V }. The additional nodes s and t represent source and sink node, respectively. The arcs capacities are defined as follows:

• Every arc (s, i), i ∈ V C is associated with a capacity (1 -q i /Q)x {0,i} ;

• Every arc (s, i), i ∈ V F is associated with a capacity x {0,i} ;

• Every arc (i, j), i ∈ V C , j ∈ V is associated with a capacity (1 -(q i + q j )/Q)x {i,j} ;

• Every arc (i, j), i ∈ V F , j ∈ V is associated with a capacity (1 -q j )x {i,j} ;

• Every arc (j, t), j ∈ V C , is associated with a capacity {j,i}∈δ(j):

j<i qi Q x {j,i} + {i,j}∈δ(j): i<j qi Q x {i,j} ;
• Every arc (j, t), j ∈ V F , is associated with a capacity 2 Q (i,j)∈A: j∈V F (S)

q i z ij + ( {j,i}∈δ(j): j<i qi Q x {j,i} + {i,j}∈δ(j): i<j qi Q x {i,j} ).
Let (S, V \ S) be the minimum s-t cut of G and assume that t ∈ S. One can see that if the cut capacity is strictly smaller than right-hand-side of inequality (80) then node set S = S \ {t} defines the most violated inequality (80). No violated inequality exists if the cut has a capacity greater than or equal to the value of the right-hand-side of (80).

Details of the heuristic algorithms 10.1 A constructive heuristic

Given an instance of vrptf, we define a complete graph G = (V , E) where the node set V = {0}∪V C contains the depot and the customer nodes. Each edge e ∈ E has a cost given by r e . Each customer i ∈ V C has a demand equal to q i and the capacity of the vehicles is set to Q. Let m = i∈Vc q i /Q be a lower bound on the minimum number of routes required. The details of our implementation of the three phases are as following.

(i) m "seed" customers are randomly selected to initialize the m routes of the emerging cvrp solution. The remaining customers are partitioned into m subsets by heuristically solving a Generalized Assignment Problem (gap) where each bin k is associated with the k-th customer selected to initialize a route. The assignment cost a ik for allocating customer i to bin k is r 0i + αr {i,k} -βr 0k , where α, β are nonnegative parameters. The gap is solved heuristically.

If, for a given m, the gap solution is infeasible, we set m = m + 1 and we repeat the above procedure.

(ii) The route for a subset of customers is determined by solving a tsp on the subgraph induced by the subset. We apply a 3-opt procedure to a starting tour obtained by generating a random sequence of the customers.

(iii) The solution obtained at step (ii) is locally optimized using a classical multiroute improvement procedure consisting of two types of operations: (a) movement of a customer from one route to another; (b) exchange of two customers belonging to different routes. We try all possible such operations until no improvement can be obtained. Each route of the modified solution is then re-optimized with the 3-opt procedure.

Since the initial partitioning and the tsp solutions are based on random choices, we can obtain different solutions executing the three phases several times. In our implementation we run them 2000 times: for the first 1000 runs we set α = 1.1 and β = 0.7 (see step (i)), while in the last 1000 runs we set β = 0, leaving α unchanged.

The cvrp solution so far obtained, is optimized by iteratively applying two re-optimization procedures: procedure Squeeze used by [START_REF] Baldacci | The Capacitated m-Ring-Star Problem[END_REF] for the cmrsp, and procedure LS-multiple. The two procedures are repeated in sequence until the current solution can be improved. Procedure Squeeze tries to re-optimize the routing (i.e. the set E) by allowing a few changes in the customer connections (set A). Procedure LS-multiple is a multiroute improvement procedure based on customer exchanges among the routes of the current solution.

A Lagrangean heuristic

Procedure CG is interwoven with a heuristic algorithm that produces a feasible vrptf solution of cost ẑ using the route sets R (see Step 2 of procedure CG). Given the current DSP (λ) solution, define vector ξ as follows:

ξ = i∈V C a i q i q(R ) ζ i , ∈ R, ( 81
) by setting ζ i (i) = 1 and ζ i = 0, ∀ ∈ R \ { (i)}, ∀i ∈ V C . Define C( ) = V C (R ) ∪ V A (R ), i.e. C(
) is the set of customers either visited on the route or assigned to facilities in V F (R ). The heuristic algorithm performs the following steps.

1. Initialization. Initialize ẑ = 0, SOL = ∅ and δ(i) = 0, and ∀i ∈ V .

2. Extract a subset of routes SOL ⊆ R. Let * be the route of R where ξ * = max{ ξ : ∈ R}.

Remove * from R. If δ(i) = 0, for some i ∈ C( * ), then update SOL = SOL ∪ { * }, δ(i) = δ(i) + a i , ∀i ∈ C( * ), and δ(i) = δ(i) + a i , ∀i ∈ V F (R * ). Repeat step 2 until R = ∅.

3. Modify the route set SOL so that δ(i) ≤ 1, ∀i ∈ V . a) Remove from SOL any route ∈ SOL such that δ(i) > 1, ∀i ∈ C( * ), and update δ(i), ∀i ∈ V , accordingly. For each ∈ SOL, compute the savings that can be achieved by removing from route every customer i ∈ C( * ) having δ(i) > 1. Let * ∈ SOL be the route of maximum saving. Remove from route * every customer i ∈ C( * ) with δ(i) > 1, and update δ(i). Repeat step 3.a until δ(i) ≤ 1, for each i ∈ V C . b) For each ∈ SOL, compute the total number α( ) of customers assigned to every facility i ∈ V F (R * ) having δ(i) > 1. Let * be the route having the minimum α( ) value. Remove from route * every facility i ∈ V F (R * ) with δ(i) > 1, update δ(i), and δ(j), ∀j ∈ V C , accordingly. Repeat step 3.b until δ(i) ≤ 1, for each i ∈ V F . c) For each ∈ SOL, remove any facility i ∈ V F (R ) with δ(i) = 1 and without customers assigned to it, and update δ(i) = δ(i) -1.

4. Insert unrouted customers. For each unrouted customer i (i.e., δ(i) = 0) perform the following operations. Compute the minimum extra-cost exc(i, ) for inserting i in route ∈ SOL without considering assignment of i to facilities in V F (R ). We set exc(i, ) = ∞ if the total load of the resulting route exceeds the vehicle capacity Q. Let * be such that exc(i, * ) = min ∈SOL [exc(i, )].

If exc(i, * ) = ∞, then set ẑ = ∞ and stop; otherwise, insert customer i in route * in the position of cost exc(i, * ) and set δ(i) = 1.

5. Define the vrptf solution ξ. Define ξ = 1, for each ∈ SOL, and ξ = 0, for each ∈ R \ SOL.

6. Local optimization. Locally optimize solution ξ by iteratively applying the two re-optimization procedures Squeeze and LS-multiple used also for the constructive heuristic. • The set of customers is selected from the customers that are currently served on a daily basis by using different criteria. The number of customers varies from a minimum of 54 up top a maximum of 164. The customer demands were computed based on historical data;

• The set of facilities corresponds to the existing set of facilities and can also include new facilities that the company want evaluate in order to revise the current distribution network. Instances with 4, 6, 7, 9, 12, 13, and 18 facilities were generated;

• Two types of fleet of vehicles were considered: vehicles with capacity equal to 24 pallets (singleunit 3 axes type of trucks) and vehicles with capacity equal to 36 pallets (single-trailer 3 axes type of trucks), respectively;

• The routing cost of a pair of nodes i and j of the network were computed as r {i,j} = c dist ij where dist ij represents the distance in kilometer between nodes i and j computed using a digital map of the Italian territory, and c is the routing cost per kilometers (currency in expressed in Euro (e)) associated with the type of vehicle (either single-unit or single-trailer);

• Set {F i } of facilities to which the customers can be assigned are defined directly by the company using different criteria. These criteria take into account the customer demand, required level of service, a priori agreements between the customers and the company, and the distance matrix [d ij ] used to compute the routing costs;

• The distribution from the facilities to the customers is performed by means of a fleet of singleunit 2 axes type of trucks with a vehicle capacity ranging from 6 to 8 pallets. The distribution cost from the facilities depends on the type of contract that has been defined between the company and the third-party contractor and vary from facility to facility. The distribution cost is a function of the number of pallets associated with the order and the distance between the customer location and the facility. Therefore the assignment cost matrix is defined by the company using the current distribution tariff agreed with the third-party companies.

lrp based instances

From each lrp instance we derived a vrptf instance as follows.

i) The set V F of facilities, the set V C of customers (and the associated demands), correspond to the set of depots and customers of the original lrp instance;

ii) The vehicle capacity Q is equal to the vehicle capacity of the original lrp instance;

iii) The depot coordinates were defined as follows. Let x min and x max be the minimum and maximum x-coordinates among the customers and the facilities x-coordinate, respectively; similarly define

  where u 0 = 0 and u 1 , . . . , u |V C | are associated with constraints (27), and u |V C |+1 , . . . , u |V | , with constraints (28). Vector v = {v 1 , . . . , v |F | } is associated with the family of valid inequalities F . 2. Generate the largest subset N ⊆ R > of ng- * route having negative reduced cost with respect to the current dual master solution (u, v) and such that |N | ≤ ∆ (∆ is an a priori defined parameter). If N = ∅, the procedure terminates; otherwise a new iteration is made. At iteration k + 1, the procedure solves a new master problem LSP by replacing R with R ∪ N and the valid inequalities of F violated by the LSP solution ξ achieved by iteration k.

  i)[START_REF] Akca | A branch-and-price algorithm for combined location and routing problems under capacity restrictions[END_REF]: 12 instances involving 5 facilities, and 30 or 40 customers;The VRPTF R. Baldacci, S. U. Ngueveu, R. Wolfler Calvo

Figure 1 :

 1 Figure 1: Real-world instances: North area

Figure 2 :

 2 Figure 2: Real-world instances: Centre area

  . . , u |V C | , u |V C |+1 , . . . , u |V | }, where u 0 = 0 for the depot, u 1 , . . . , u |V C | are associated with constraints (27), and u |V C |+1 , . . . , u |V | , with constraints (28). In the following, we denote with q min = min i∈V C {q i }. The following theorem holds. Let us associate penalties λ i ∈ R, ∀i ∈ V C , with constraints (27), and

	Theorem 2

Table 2 :

 2 Summary results on Class A instances%U B 1 %U B 2 %LB 1 %LB 2 t DA %LB CP t CP %LB t LB #Opt t T OT

	Akca et al. (2009) 100.3 100.6	94.3	96.8	9.9	96.1	4.6 98.5	2.1 10/12 145.3
	Prins et al. (2004) 100.2 100.4	93.7	96.0 48.6	94.0 78.6 97.8 14.7 10/24 221.9
	Different authors 100.2 101.0	91.9	94.3 297.5	93.7 339.8 96.6 121.3 10/39 213.0

Table 3 :

 3 Summary results on Class B instances%U B 1 %U B 2 %LB 1 %LB 2 t DA %LB CP t CP %LB t LB #Opt t T OT

	Akca et al. (2009) 102.6 101.3	94.6	96.9	5.2	95.6	4.1 98.1	3.7 9/12 274.2
	Prins et al. (2004) 101.4 101.0	94.2	95.9 52.8	93.0 76.1 97.2 14.2 7/24 184.0
	Different authors 101.1 102.5	92.4	94.2 187.6	93.4 324.1 96.1 141.2 8/39 284.3

Table 4 :

 4 Effectiveness of the dual ascent heuristics

	t LB SP	t LB 1

(e) route set R initialized with single-customer route (f) route set R initialized with the solution provided by the constructive heuristic customers and 10 facilities were solved to optimality.

Table 5 :

 5 Effectiveness of the different type of inequalities on column-and-cut generation procedure no cuts + CI + MI + RCI + RCII-a + RCII-b + RCII-c + RCII-d %LB t LB %LB t LB #cuts %LB t LB #cuts %LB t LB #cuts

	A	Akca et al. (2009) 96.9 2.3 97.9	3.1	7.3 98.5	4.0 200.6 98.5	4.0 169.9
		Prins et al. (2004) 96.1 5.3 97.1 10.0	8.3 97.8 17.2 521.5 97.8 18.2 539.8
		Different authors 94.7 53.5 95.7 92.1	23.8 96.5 139.3 1083.2 96.6 164.4 1181.3
	B	Akca et al. (2009) 96.9 2.7 97.4	3.6	13.3 98.0	5.5 655.6 98.1	5.6 763.0
		Prins et al. (2004) 96.0 5.8 96.3 10.7	14.6 97.1 16.2 623.8 97.2 17.7 737.0
		Different authors 94.8 56.9 95.6 109.0	61.4 96.0 173.3 1385.3 96.1 184.3 1566.1
	Real-word	98.6 4.2 98.7	6.6	4.1 98.7	6.9 106.0 98.7 12.9 121.1
		95.8 28.0 96.5 50.8	24.9 97.1 78.8 809.2 97.2 88.1 905.2

Table 6 :

 6 Results on Class A:[START_REF] Akca | A branch-and-price algorithm for combined location and routing problems under capacity restrictions[END_REF] lrp based instances Name z * #r #f #c %U B1 %U B2 %LB1 %LB2 t DA %LB C t C %LB t LB #cuts #cols #N %Opt t T OT

	cr30x5a-1 621	5	0	0 100.0 100.0	99.4 100.0 2.6	98.8 2.7 100.0 0.1	0	303	1 100.0	8
	cr30x5a-2 665	5	1	2 100.2 102.1	90.6	94.0 2.8	97.5 2.5 99.1 3.2	275 2185	23 100.0	20
	cr30x5a-3 575	5	1	2 100.7 100.0	90.8	96.5 2.5	98.5 2.2 98.5 1.3	252 1160	93 100.0	29
	cr30x5b-1 727	5	1	1 100.0 100.0	96.4	98.5 2.9	99.7 2.9 100.0 1.3	479	959	1 100.0	10
	cr30x5b-2 826	6	0	0 100.0 100.6	92.5	93.6 1.5	93.0 3.7 97.6 2.1	78 1609 367 100.0	79
	cr30x5b-3 788	7	0	0 100.1 100.0	94.6	96.1 32.5	94.2 2.8 97.6 1.9	94 1096 1335 100.0 1061
	cr40x5a-1 738	7	3	8 100.0 100.0	95.5	97.5 19.5	94.9 4.2 97.9 2.5	98 1530 151 100.0	82
	cr40x5a-2 786	6	1	1 100.5 101.0	92.7	96.9 5.2	95.8 8.3 98.2 3.8	357 1627 986	99.3 3615
	cr40x5a-3 807	6	4	4 101.1 101.5	94.8	98.0 21.4	95.2 5.3 98.5 2.8	405 1212 562	98.8 3631
	cr40x5b-1 964	8	0	0 100.0 101.0	95.9	98.1 11.1	95.6 6.0 98.9 1.3	238	646 187 100.0	49
	cr40x5b-2 901	8	2	3 100.0 100.1	93.1	95.7 3.3	95.5 7.9 98.1 3.4	372	910	41 100.0	34
	cr40x5b-3 887	8	2	5 100.5 100.5	95.2	96.7 13.5	95.0 6.4 98.1 1.5	295	876 378 100.0	81

Table 7 :

 7 Results on Class A: Prins et al. (2004) lrp based instances

	Name	z * #r #f #c %U B1 %U B2 %LB1 %LB2 t DA %LB C	t C %LB t LB #cuts #cols	#N %Opt t T OT
	ppw-20-5-0-a 253	5	1	1 100.0 100.0	95.2	96.7	1.5	98.5	1.2 100.0 0.2	11	394	1 100.0	5
	ppw-20-5-0-b 211	3	1	2 100.0 100.0	85.3	94.4	7.8	100.0	0.5 97.8 0.6	120	803	17 100.0	12
	ppw-20-5-2-a 247	5	1	3 100.0 100.0	90.8	94.5	3.2	95.4	1.3 97.4 0.4	157	711	197 100.0	14
	ppw-20-5-2-b 189	3	1	2 100.0 100.0	85.7	93.1	1.6	100.0	0.3 100.0 1.1	181 2212	1 100.0	5
	ppw-50-5-0-a 616 12	1	1 100.2 100.2	97.0	98.9	1.6	93.9 14.8 98.9 2.3	224	794	332 100.0	87
	ppw-50-5-0-b 400	6	1	2 100.0 102.0	91.8	95.6	6.3	96.8	9.8 96.7 11.9	362 2964	15	96.9 3621
	ppw-50-5-2'-a 653 12	0	0 100.5 100.0	95.9	96.1	1.6	95.6 13.4 99.5 2.6	323	827	85 100.0	48
	ppw-50-5-2'-b 351	6	0	0 100.3 100.0	91.1	92.9	5.3	99.0 10.6 99.4 6.9	973 2391	204	99.7 3619
	ppw-50-5-2-a 587 12	1	2 100.0 100.0	96.7	98.3	1.5	93.6 11.1 98.9 2.1	828	880	949 100.0	126
	ppw-50-5-2-b 357	6	0	0 100.0 100.3	92.3	93.4 11.3	96.7 13.3 96.7 6.5	514 2174	61	97.1 3625
	ppw-50-5-3-a 586 12	1	3 100.2 100.0	94.8	95.4 20.5	92.2 12.4 95.9 2.0	152	728 20734	97.9 3630
	ppw-50-5-3-b 381	6	0	0 100.0 100.0	91.2	94.4 28.3	95.9	8.4 96.8 6.4	493 3335	395	97.4 3644
	ppw-100-5-0-a 1158 25	2	2 100.9 100.9	97.7	98.8	3.1	93.8 216.9 99.3 10.6	12 1417 1312 100.0	507
	ppw-100-5-0-b 679 11	3	3 100.0 100.0	91.2	95.3 146.6	91.6 115.2 96.7 42.6	2019 5451	7	96.8 3819
	ppw-100-5-2-a 1010 24	1	1 100.6 100.4	96.8	97.4	3.1	91.7 193.0 97.5 14.0	640 1672 5639	97.8 3735
	ppw-100-5-2-b 569 12	1	1 100.0 100.2	94.8	95.6 99.5	90.0 110.1 96.0 30.8	1040 3019	302	96.2 3783
	ppw-100-5-3-a 1068 24	1	1 100.7 100.7	97.1	98.2 29.1	93.6 146.8 98.6 8.8	561 1312 8235	98.9 3686
	ppw-100-5-3-b 612 11	0	0 100.0 100.7	92.7	95.6 55.7	93.8 131.5 97.1 27.9	1105 5207	183	97.4 3733
	ppw-100-10-0-a 1215 24	1	1 100.1 100.7	97.9	98.7 47.9	89.9 154.4 98.9 11.2	372 1040 2420 100.0	991
	ppw-100-10-0-b 693 11	1	1 100.0 101.9	94.2	96.6 208.8	93.2 159.1 97.4 38.3	1006 4232	84	97.6 3869
	ppw-100-10-2-a 1030 24	0	0 100.8 100.0	96.2	97.7 134.6	90.6 163.8 97.7 9.9	85 1182 6064	98.2 3853
	ppw-100-10-2-b 582 11	0	0 100.0 100.9	93.4	95.7 127.1	91.2 133.3 96.4 54.6	829 4798	30	96.5 3809
	ppw-100-10-3-a 1055 24	4	6 100.0 100.3	97.6	98.3 121.6	88.7 142.1 99.3 14.8	30 1327	245 100.0	424
	ppw-100-10-3-b 608 11	1	2 100.0 101.6	91.2	92.9 98.7	89.7 123.0 94.6 46.8	919 6930	24	94.8 3773

Table 8 :

 8 Results on Class A: different authors lrp based instances

	Name	z * #r #f #c %U B1 %U B2 %LB1 %LB2	t DA %LB C	t C %LB t LB #cuts #cols #N %Opt t T OT
	Christ-50x5 514	5	2	2 100.6 103.1	95.7	97.3	58.1	98.2	8.3 99.2	7.9	1086 3231	54 100.0	130
	Christ-50x5 B 533	5	1	2 100.2 100.4	92.5	95.9	84.8	96.4	7.9 96.8	7.7	592 4284 175	97.2 3698
	Christ-75x10 783	9	3	6 100.0 100.0	93.4	94.6	40.0	93.6	58.2 97.1 28.6	665 6206	22	97.3 3676
	Christ-75x10 B 814	9	3	5 100.0 101.8	94.0	95.2 112.9	94.2	46.0 97.5 33.3	1096 6417 235	97.8 3750
	Christ-100x10 831	8	0	0 100.0 101.4	92.9	93.8 321.8	94.8	63.6 96.8 126.1	1335 17779	5	96.9 3987
	Gaskell-21x5 371	4	1	2 100.0 100.0	97.1	98.7	1.3	98.1	1.3 100.0	0.3	157	534	1 100.0	5
	Gaskell-22x5 554	3	3	4 102.0 100.0	82.3	88.6	83.7	97.3	1.0 99.8 46.0	301 3410	5 100.0	145
	Gaskell-29x5 503	4	1	1 102.2 100.0	88.3	94.0 119.1	93.7	1.6 98.0 67.1	223 1936	67 100.0	683
	Gaskell-32x5-2 427	3	0	0 100.0 100.0	92.4	98.9 459.5	100.0	1.4 100.0 97.3	6		1 100.0	567
	Gaskell-32x5 479	4	1	1 100.0 100.0	91.7	95.7 224.3	98.7	2.4 100.0 49.8	330 2149	1 100.0	280
	Gaskell-36x5 411	4	1	1 100.2 100.2	96.4	96.7	8.8	99.0	3.7 100.0	1.3	140 1875	1 100.0	17
	Min-27x5 3083	4	1	1 100.0 100.0	89.4	95.2	17.1	99.2	1.7 100.0	1.5	364 1357	1 100.0	24
	Perl83-12x2 100	2	0	0 100.0 100.0	92.1	99.3	0.8	100.0	0.1 100.0	0.1	25	417	1 100.0	2
	Perl83-55x15 453 10	3	3 101.5 101.3	96.8	97.9	38.2	94.6	17.5 99.3	3.5	324 2422 467 100.0	278
	Perl83-85x7 618 11	1	1 100.2 101.1	96.9	97.7	36.1	92.7	59.5 98.2 13.0	580 3516 1835	98.8 3736
	P111112-100x10 1346 11	0	0 100.0 100.6	92.1	94.5 163.2	92.2 112.2 95.4 46.4	976 7651 196	95.6 3846
	P111122-100x20 1252 11	1	2 100.0 102.6	93.8	96.1 453.1	93.3 148.3 98.7 58.2	70 5927 625	98.8 4138
	P111212-100x10 1266 10	0	0 100.0 100.6	92.8	95.8	46.8	93.1 118.8 96.8 57.4	1144 4786	5	96.9 3718
	P111222-100x20 1338 11	1	1 100.0 100.4	91.5	94.0 379.0	90.8 208.6 96.1 74.5	1542 5446	19	96.1 4053
	P112112-100x10 1236 11	3	3 100.0 100.0	89.6	93.2 196.6	93.1 173.0 96.7 134.8	2344 9884	15	96.8 3889
	P112122-100x20 1047 10	3	3 100.0 100.0	84.6	86.8 485.6	92.4 278.6 94.2 227.5	2348 16311	3	94.2 4177
	P112212-100x10 892 11	2	2 100.4 100.0	89.0	90.6 220.3	88.6 107.2 92.1 83.5	1045 9920	15	92.1 3918
	P112222-100x20 1006 10	1	1 100.0 103.0	93.4	94.4	93.4	94.7 170.8 95.7 176.1	1039 7578	8	95.7 3747
	P113112-100x10 1158 11	0	0 100.0 102.9	89.5	91.8 319.1	93.8 176.5 94.0 58.5	1295 6782	8	94.2 4007
	P113122-100x20 1190 11	4	6 100.0 102.3	87.8	90.3 227.5	92.9 221.0 96.1 163.0	2200 10465	13	96.1 3914
	P113212-100x10 1154 10	1	1 100.0 104.9	92.9	93.6	48.9	93.1 125.8 95.2 63.8	927 4778	97	95.3 3717
	P113222-100x20 1078 11	0	0 100.0 100.0	90.6	91.8	73.2	94.4 235.0 94.3 101.7	1164 5276	68	94.5 3748
	P131112-150x10 1833 16	1	1 100.8 100.0	93.7	95.1 171.8	90.2 669.4 95.6 134.3	1112 7283	37	95.7 3946
	P131122-150x20 1769 16	1	1 100.0 100.8	92.5	95.5 579.3	89.8 872.2 95.9 137.6	34 9383	56	95.9 4411
	P131212-150x10 1802 16	1	2 100.0 101.2	93.5	96.2 318.2	91.1 500.7 97.1 171.2	2076 9637	18	97.2 4147
	P131222-150x20 1802 15	2	2 100.0 100.5	93.2	95.2 600.7	89.8 892.7 95.6 119.4	870 7002	87	95.7 4371
	P132112-150x10 1783 16	3	3 100.0 101.0	91.8	95.0 969.1	93.3 1147.1 96.5 286.1	3288 11779	60	96.7 4815
	P132122-150x20 1541 15	1	1 100.0 100.1	88.8	90.4 731.2	91.0 1240.3 93.8 429.3	3329 23858	2	93.8 4508
	P132212-150x10 1251 16	0	0 100.0 101.0	91.5	92.6 240.4	90.8 766.6 94.3 180.1	2126 11721	13	94.3 4073
	P132222-150x20 1184 16	1	1 100.0 100.0	92.7	93.7 513.2	89.9 985.7 95.8 514.0	2487 11655	3	95.8 4368
	P133112-150x10 1899 16	2	2 100.0 103.2	92.4	93.8 617.6	93.0 1182.5 94.7 409.5	2545 12407	6	94.8 4434
	P133122-150x20 1498 16	2	2 100.0 100.8	92.5	93.4 1061.3	90.8 755.5 94.5 190.7	1966 9534	88	94.6 4908
	P133212-150x10 1245 16	1	1 100.0 103.9	93.1	94.1 649.9	93.4 818.1 95.6 108.3	1602 7306	39	95.7 4510
	P133222-150x20 1551 16	0	0 100.0 100.1	90.3	91.0 836.0	89.3 1069.9 91.5 321.5	1317 10764 222	91.6 4688

Table 9 :

 9 Results on Class B:[START_REF] Akca | A branch-and-price algorithm for combined location and routing problems under capacity restrictions[END_REF] lrp based instances Name z * #r #f #c %U B1 %U B2 %LB1 %LB2 t DA %LB C t C %LB t LB #cuts #cols #N %Opt t T OT

	cr30x5a-1 4176	5	2	8 100.0 100.2	98.8	99.2 3.4	97.0 2.6 99.5 1.9	426	770	13 100.0	13
	cr30x5a-2 4428	5	1	6 100.6 103.3	89.9	92.8 3.2	96.3 3.2 98.2 5.2	338 1431	50	98.3 3609
	cr30x5a-3 3655	5	2 11 101.8 100.5	91.5	96.7 3.2	98.7 2.4 98.2 3.3	633 1610	34 100.0	24
	cr30x5b-1 4844	5	1	5 102.4 100.0	94.7	97.1 3.0	97.3 3.0 98.4 5.9	893 1141	22 100.0	26
	cr30x5b-2 4931	6	4 13 103.3 102.0	95.9	97.2 4.1	95.2 3.1 98.6 0.9	236	656 100 100.0	28
	cr30x5b-3 4626	6	3 15 108.3 100.6	95.4	96.5 1.7	94.5 2.1 96.6 0.4	1	707 3330	98.9 3605
	cr40x5a-1 4221	6	4 21 100.0 100.5	95.4	96.8 3.5	94.6 5.1 97.4 3.0	15 2112	90 100.0	50
	cr40x5a-2 4804	6	4 17 102.0 101.4	91.6	96.9 4.3	94.6 7.1 97.3 5.4	2647	919 282	98.3 3616
	cr40x5a-3 4577	6	4 25 106.9 103.9	94.9	97.4 7.4	96.7 3.7 98.8 10.4	1546 1605 213 100.0 2085
	cr40x5b-1 6334	9	3 13 102.5 100.5	94.1	96.3 15.2	92.0 5.4 97.0 2.4	592	993 175 100.0	74
	cr40x5b-2 5933	8	3 14 100.6 102.1	95.3	96.8 7.5	93.4 6.3 97.2 2.4	219	930 362 100.0	136
	cr40x5b-3 5279	8	3 14 102.7 100.9	97.3	98.6 5.6	97.1 5.1 99.4	1610	955	45 100.0	33

Table 10 :

 10 Results on Class B:[START_REF] Prins | Nouveaux algorithmes pour le problème de localisation et routage sous contraintes de capacité[END_REF] lrp based instances Name z * #r #f #c %U B1 %U B2 %LB1 %LB2 t DA %LB C t C %LB t LB #cuts #cols #N %Opt t T OT

	ppw-20-5-0-a 1535	5	2	7 101.2 100.0	95.3	97.0	0.8	95.6	1.4 99.4 0.4	29	342	7 100.0	4
	ppw-20-5-0-b 1315	3	1	8 103.3 100.0	93.0	94.4	1.3	99.2	0.6 99.6 7.6	632 3212	5 100.0	13
	ppw-20-5-2-a 1488	5	3 11 102.1 100.0	94.8	95.6	2.1	92.4	1.4 97.3 0.4	300	324	133 100.0	11
	ppw-20-5-2-b 1085	3	2 11 103.0 100.0	94.7	99.2	1.9	100.0	0.4 100.0 0.2	1	539	1 100.0	4
	ppw-50-5-0-a 4159 12	3 12 102.3 100.7	97.5	99.3	1.7	93.5 15.7 99.3 0.9	7	650	645 100.0	96
	ppw-50-5-0-b 2638	6	3 15 100.0 104.2	91.0	96.5 17.6	97.0 13.1 97.1 8.5	981 2751	283	97.3 3633
	ppw-50-5-2'-a 4522 12	2	9 100.0 100.5	95.2	95.6 19.3	94.7 16.4 98.1 2.0	316	849 1057 100.0	214
	ppw-50-5-2'-b 2439	6	1	7 100.9 100.3	90.8	92.5	7.7	97.7 12.1 98.9 9.8	630 2811	847	99.1 3633
	ppw-50-5-2-a 3910 12	2	8 101.8 101.0	97.0	98.3 18.5	92.8 11.3 98.3 0.7	342	719 8849 100.0	946
	ppw-50-5-2-b 2389	6	2 10 101.1 100.0	90.7	91.8	7.6	93.6 11.2 94.9 12.8	1414 2733	60	95.2 3623
	ppw-50-5-3-a 3649 12	4 14 103.8 101.3	95.3	95.7 39.6	93.7 12.8 95.9 0.8	2	673 19638	98.3 3668
	ppw-50-5-3-b 2421	6	3 25 100.0 100.6	90.0	93.8	4.1	92.9	5.7 95.7 9.6	1596 3006	199	96.2 3620
	ppw-100-5-0-a 8009 25	3 10 101.5 101.1	97.2	98.5	3.1	93.3 210.1 98.7 7.2	11 1262 12954	99.6 3689
	ppw-100-5-0-b 4629 11	4 13 100.0 103.7	93.0	95.7 219.0	92.4 116.9 97.3 50.9	3561 4394	38	97.4 3891
	ppw-100-5-2-a 6838 24	3 13 104.1 100.0	96.6	97.2 105.8	91.6 169.9 97.3 6.5	645 1672 8072	97.7 3838
	ppw-100-5-2-b 3925 11	3 20 100.4 100.0	93.6	94.8 110.1	89.8 133.6 94.9 17.5	657 2893	326	95.1 3793
	ppw-100-5-3-a 7184 24	4 15 103.7 102.7	97.0	97.9 18.2	92.0 124.9 98.3 10.2	1291 1617 5523	98.7 3785
	ppw-100-5-3-b 4141 11	4 21 100.0 101.5	93.3	95.3 27.3	92.0 133.4 96.1 44.6	3258 5104	50	96.3 3706
	ppw-100-10-0-a 7960 24	8 28 101.2 100.0	95.5	96.8 166.0	88.3 109.0 96.8 5.9	2 1030 6269	97.4 3824
	ppw-100-10-0-b 4698 26	6 35 100.0 103.6	93.4	95.3 61.8	91.7 164.9 95.6 24.0	274 3132	55	95.8 3721
	ppw-100-10-2-a 6883 23	6 24 101.8 100.0	94.9	96.4 115.1	90.8 178.6 96.4 6.5	1 1988 2439	96.6 3831
	ppw-100-10-2-b 3984 11	5 22 100.0 100.8	92.5	94.4 95.3	90.1 120.4 94.5 35.3	37 2982	71	94.6 3780
	ppw-100-10-3-a 7060 24	9 32 102.3 100.0	95.6	96.4 147.3	87.4 124.5 97.1 10.5	5 1191 3973	97.5 3802
	ppw-100-10-3-b 4081 11	6 22 100.0 102.0	92.2	94.1 75.2	90.5 138.9 95.2 68.3	1697 5956	92	95.3 3750

Table 11 :

 11 Results on Class B: different authors lrp based instancesName z * #r #f #c %U B1 %U B2 %LB1 %LB2 t DA %LB C t C %LB t LB #cuts #cols #N %Opt t T OT

	Christ-50x5 3157	5	4 24 100.0 105.3	97.5	98.8	7.0	99.1	10.4 99.6	8.8	1410 2341 27 100.0	125
	Christ-50x5 B 3238	5	5 21 103.9 103.0	96.5	98.5	6.7	99.0	8.8 99.3 11.1	1137 2467 63	99.7 3619
	Christ-75x10 5241	9	5 17 100.0 101.4	92.1	92.9 60.4	90.8	45.9 94.4 42.1	1498 4130 83	94.6 3697
	Christ-75x10 B 5394	9	5 18 100.0 102.5	93.5	94.6 50.8	92.2	37.8 96.0 55.8	1375 4226	3	96.0 3687
	Christ-100x10 5377	8	5 34 100.0 103.9	92.3	93.8 70.2	93.3	63.8 94.9 91.9	813 9094 19	94.9 3735
	Gaskell-21x5 2057	4	3 14 100.6 104.6	96.8	98.4	1.8	98.3	1.0 99.4	0.4	285	555	5 100.0	5
	Gaskell-22x5 3097	3	3 14 105.5 100.0	78.4	75.9 86.7	96.4	0.8 99.6 147.4	560 3709	7 100.0	268
	Gaskell-29x5 3035	3	2	5 103.9 105.0	91.4	96.2 103.1	93.7	1.6 99.4 177.0	69 3273 24	99.4 3708
	Gaskell-32x5-2 2537	3	1	8 101.7 101.7	89.1	97.2 488.0	98.8	2.1 99.4 262.4	305 2208 39 100.0 1320
	Gaskell-32x5 2691	4	2	9 105.1 100.0	97.1	99.2 206.7	99.3	1.7 99.5 55.2	24 2038	7 100.0	314
	Gaskell-36x5 2355	4	3 24 109.8 100.0	98.7	98.8	1.3	99.8	3.8 100.0	6.8	1699 2661	5 100.0	19
	Min-27x5 20229	4	2 10 103.9 102.2	82.6	93.8	4.0	95.1	1.4 96.5	4.9	368 1738 100	96.9 3848
	Perl83-12x2	519	2	2 12 101.2 100.0	99.3	99.4	0.5	100.0	1.0 100.0	0.5	32	151	1 100.0	1

Acknowledgements

The authors thank two anonymous referees for making several suggestions that improved the presentation of the paper.

The VRPTF R. Baldacci, S. U. Ngueveu, R. Wolfler Calvo y min and y max . The coordinate (x, y) of the central depot are defined as follows:

and y = y min + (y max -y min )/2 . (82)

The routing and connection costs were generated as follows.

Class A. Routing and assignment costs of a pair of nodes i, j are equal to the Euclidean distance e ij , computed according to the TSPLIB EUC 2D standard.

Class B. For each pair of nodes i, j, the routing cost is r {i,j} = αe ij , while the assignment cost is

, where α = 7.0.

For all the instances, every customer can be assigned to every facility, i.e.

We generated a total number of 150 instances, 75 instances per class. The dimensions of the instances vary from very small instances with 12 customers and two facilities up to very large instances with 150 customers and 20 facilities.

12 Details about the computational results on lrp based instances This section reports the complete details about the computational results on lrp based instances.