Artur Pessoa

Ruslan Sadykov
email: ruslan.sadykov@inria.fr

Eduardo Uchoa
email: uchoa@producao.uff.br

François Vanderbeck

A Generic Exact Solver for Vehicle Routing and Related Problems *

Major advances were recently obtained in the exact solution of Vehicle Routing Problems (VRPs). Sophisticated Branch-Cut-and-Price (BCP) algorithms for some of the most classical VRP variants now solve many instances with up to a few hundreds of customers. However, adapting and reimplementing those successful algorithms for other variants can be a very demanding task. This work proposes a BCP solver for a generic model that encompasses a wide class of VRPs. It incorporates the key elements found in the best existing VRP algorithms: ng-path relaxation, rank-1 cuts with limited memory, path enumeration, and rounded capacity cuts; all generalized through the new concepts of "packing set" and "elementarity set". The concepts are also used to derive a branching rule based on accumulated resource consumption and to generalize the Ryan and Foster branching rule. Extensive experiments on several variants show that the generic solver has an excellent overall performance, in many problems being better than the best specific algorithms. Even some non-VRPs, like bin packing, vector packing and generalized assignment, can be modeled and effectively solved.

Introduction

Since its introduction by Dantzig and Ramser [START_REF] Dantzig | The truck dispatching problem[END_REF], the Vehicle Routing Problem (VRP) has been one of the most widely studied in combinatorial optimization. Google Scholar indicates that 691 works containing the exact string "vehicle routing" in the title were published only in 2018. VRP relevance stems from its direct use in the real systems that distribute goods and provide services, vital to the modern economies. Reflecting the large variety of conditions in those systems, the VRP literature is spread into dozens, perhaps hundreds, of variants. For example, there are variants that consider capacities, time windows, heterogeneous fleets, multiple depots, pickups and deliveries, optional customer visits, arc routing, etc.

In recent years, big advances in the exact solution of VRPs have been accomplished. A milestone was certainly the Branch-Cut-and-Price (BCP) algorithm of [START_REF] Pecin | Improved branch-cut-and-price for capacitated vehicle routing[END_REF][START_REF] Pecin | Improved branch-cut-andprice for capacitated vehicle routing[END_REF], that could solve Capacitated VRP (CVRP) instances with up to 360 customers, a large improvement upon the previous record of 150 customers. That algorithm exploits many elements introduced by several authors, combining and enhancing them. In particular, the new concept of limited memory cut proved to be pivotal. Improvements of the same magnitude were later obtained for a number of classical variants like VRP with Time Windows (VRPTW) [START_REF] Pecin | New enhancements for the exact solution of the vehicle routing problem with time windows[END_REF], Heterogeneous Fleet VRP (HFVRP) and Multi Depot VRP (MDVRP) [START_REF] Pessoa | Enhanced branch-cut-and-price algorithm for heterogeneous fleet vehicle routing problems[END_REF], and Capacitated Arc Routing (CARP) [START_REF] Pecin | Comparative analysis of capacitated arc routing formulations for designing a new branch-cut-and-price algorithm[END_REF]. For all those variants, instances with about 200 customers are now likely to be solved, perhaps in hours or even days. However, there is something even more interesting: many instances with about 100 customers, that a few years ago would take hours, are solved in less than 1 minute. This means that many more real world instances can now be tackled by exact algorithms in reasonable times.

Unhappily, designing and coding each of those complex and sophisticated BCPs has been a highly demanding task, measured on several work-months of a skilled team. In effect, this prevents the use of those algorithms in real world problems, that actually, seldom correspond exactly to one of the most classical variants. This work presents a framework that can handle most VRP variants found in the literature and can be used to model and solve many other new variants. In order to obtain state-of-the-art BCP performance, some key components found in the best specific VRP algorithms had to be generalized. The new concepts of packing set and elementarity set were instrumental for that.

The quest for general exact VRP algorithms can be traced back to Balinski and Quandt [START_REF] Balinski | On an integer program for a delivery problem[END_REF], where a set partitioning formulation valid for many variants was proposed. That formulation had only turned practical in the 1980's and 1990's, when the Branch-and-Price (BP) method was developed. At that time, it was recognized that the pricing subproblems could often be modeled as Resource Constrained Shortest Path (RCSP) problems and solved by labeling algorithms, leading to quite generic methods (for example, Desaulniers et al. [START_REF] Desaulniers | A unified framework for deterministic time constrained vehicle routing and crew scheduling problems[END_REF]). However, those BP algorithms only worked well on problems with "tightly constrained" routes, like VRPTW with narrow time windows. Many variants, including CVRP, were much better handled by Branch-and-Cut (BC) algorithms using problem-specific cuts (for example, Lysgaard et al. [START_REF] Lysgaard | A new branch-and-cut algorithm for the capacitated vehicle routing problem[END_REF]). In the late 2000's decade, after works like [START_REF] Fukasawa | Robust branch-and-cut-and-price for the capacitated vehicle routing problem[END_REF][START_REF] Baldacci | An exact algorithm for the vehicle routing problem based on the set partitioning formulation with additional cuts[END_REF][START_REF] Jepsen | Subset-row inequalities applied to the vehicle-routing problem with time windows[END_REF][START_REF] Desaulniers | Tabu search, partial elementarity, and generalized k-path inequalities for the vehicle routing problem with time windows[END_REF][START_REF] Ropke | Branch and cut and price for the pickup and delivery problem with time windows[END_REF][START_REF] Baldacci | An exact algorithm for the pickup and delivery problem with time windows[END_REF][START_REF] Bode | Cut-first branch-and-price-second for the capacitated arc-routing problem[END_REF], it became clear that the combination of cut and column generation performs better than pure BP or pure BC on almost all problems. Until today, BCP remains the dominant VRP approach. A first attempt of a generic algorithm based on cut and column generation was presented in Baldacci and Mingozzi [START_REF] Baldacci | A unified exact method for solving different classes of vehicle routing problems[END_REF], where 7 variants, all of them particular cases of the HFVRP, could be solved. Recently, [START_REF] Sadykov | A bucket graph based labeling algorithm with application to vehicle routing[END_REF] proposed a BCP for several particular cases of the HFVRP with time windows. The framework now proposed is far more generic than that.

The Basic Model

In this section, we provide a formal definition of a generic model that can be solved by a branchcut-and-price algorithm, where all pricing subproblems are modeled as RCSP problems. The main original feature of the basic model is the concept of mapping for linking the variables induced by the resource constrained paths to a user-defined generic objective function and constraints in the Master formulation.

Graphs for RCSP Subproblems

Let K be the set of subproblems. Define directed graphs G k = (V k , A k ), k ∈ K. Let V = ∪ k∈K V k and A = ∪ k∈K A k . The graphs are not necessarily simple and may even have loops. Vertices and arcs in all graphs are distinct. Each graph has special source and sink vertices: v k source and v k sink . The source and sink may be distinct vertices, but may also be the same vertex. Define set R k of resources. For each r ∈ R k and a ∈ A k , q a,r ∈ R is the consumption of resource r in arc a. There are finite accumulated resource consumption intervals [l a,r , u a,r ], for each r ∈ R k and a ∈ A k . Since in most applications these intervals are more naturally defined on vertices, we may define intervals [l v,r , u v,r ], v ∈ V k , meaning that [l a,r , u a,r ] = [l v,r , u v,r ] for every arc a ∈ δ -(v) (i.e., entering v). Resources without any negative consumption are called monotone, otherwise they are non-monotone. Set R k is divided into main resources R k M and secondary resources R k N . Main resources should be monotone, secondary resources may be monotone or non-monotone. As will be discussed in Section 5, the concept of main resource is directly related to key implementation issues in the labeling algorithm used for solving the pricing. For each cycle in G k , the net consumption of at least one resource should be strictly positive. Therefore, unless G k is acyclic, the existence of at least one resource is mandatory. Finally, resources are also classified as disposable or non-disposable. In the former case, it is possible to drop resources if this is needed order to satisfy the accumulated resource consumption intervals. Main resources are always disposable. By default, secondary resources are also assumed to be disposable. The existence of non-disposable resources should be explicitly indicated. A resource constrained path p = (v k source = v 0 , a 1 , v 1 , . . . , a n-1 , v n-1 , a n , v n = v k sink ) over a graph G k should have n ≥ 1 arcs, v j = v k source and v j = v k sink , 1 ≤ j ≤ n -1, and is feasible if:

• for every r ∈ R k that is disposable, the accumulated resource consumption S j,r at visit j, 0 ≤ j ≤ n, where S 0,r = 0 and S j,r = max{l aj ,r , S j-1,r + q aj ,r }, does not exceed u aj ,r ;

• for every r ∈ R k that is non-disposable, the accumulated resource consumption S j,r at visit j, 0 ≤ j ≤ n, where S 0,r = 0 and S j,r = S j-1,r + q aj ,r , lies in the interval [l aj ,r , u aj ,r ].

Note that some feasible paths may not be elementary, some vertices or arcs being visited more than once. For each k ∈ K, let P k denote the set of all feasible resource constrained paths in G k . Each set P k is finite, either because G k is acyclic or because the resources limit the number of times that each vertex or arc can be visited. Define P = ∪ k∈K P k . As vertices and arcs in different graphs are distinct, paths in different graphs are also distinct.

Master Formulation

For all a ∈ A and p ∈ P , let h p a indicate how many times arc a appears in path p. Given c ∈ R n1 , f ∈ R n2 , α ∈ R m×n1 , β ∈ R m×n2 , and d ∈ R m , the problem is formulated as follows.

Min

α ij x j + n2 s=1 β is y s ≥ d i , i = 1, . . . , m, (1b)
x j = k∈K p∈P k a∈M (xj ) h p a λ p , j = 1 . . . , n 1 , (1c)
L k ≤ p∈P k λ p ≤ U k , k ∈ K, (1d)
λ p ∈ Z + , p ∈ P, (1e)
x j ∈ Z, y s ∈ Z, j = 1, . . . , n1 , s = 1, . . . , n2 , (1f)
where x j , 1 ≤ j ≤ n 1 , y s , 1 ≤ s ≤ n 2 , and λ p , p ∈ P , are variables. The first n1 x variables and the first n2 y are defined to be integer; all λ variables are non-negative integer. Equation (1a) defines a general linear objective function, c and f are cost vectors. Inequalities (1b) define m general linear constraints over those variables, α and β are the coefficient matrices and d the right-hand side vector. Constraints (1b) may even contain exponentially large families of cuts, provided that suitable procedures are given for their separation. However, by simplicity, we continue the presentation as if all the m constraints are explicitly defined. For each variable x j , 1 ≤ j ≤ n 1 , M (x j ) ⊆ A defines its mapping into a non-empty subset of the arcs. We remark that mappings do not need to be disjoint, the same arc can mapped to more than one variable x j . Define M -1 (a) as {j | a ∈ M (x j )}. As not all arcs need to belong to some mapping, some M -1

sets may be empty. The relation between variables x and λ is given by (1c). For each k ∈ K, L k and U k are given lower and upper bounds on the number of paths from G k in a solution.

A feasible solution to Formulation (1) is composed of a set of paths, each path p ∈ P with multiplicity λ p in the solution, and additional decisions represented by the values assigned to variables y s , s = 1, . . . , n 2 . Hence, modelling a new problem as [START_REF] Lemon | Library for Efficient Modeling and Optimization in Networks[END_REF] requires it to contain one or more structures that can be cast into paths in one or more properly defined graphs. Then, resources should be created to model "intrapath" constraints while global "interpath" constraints, and the objective function, should be modelled as (1b), and (1a), respectively. Note that the values of x j , j = 1, . . . , n 1 , are completely defined as a function of the paths that compose the solution (and their multiplicities) through the mappings. Thus, these variables should only be created for allowing formulating (1a) and (1b).

Eliminating the x variables and relaxing the integrality constraints, the following LP is obtained:

Min k∈K p∈P k n1 j=1 c j a∈M (xj ) h p a λ p + n2 s=1 f s y s (2a) S.t. k∈K p∈P k n1 j=1 α ij a∈M (xj ) h p a λ p + n2 s=1 β is y s ≥ d i , i = 1, . . . , m, (2b)
L k ≤ p∈P k λ p ≤ U k , k ∈ K, (2c)
λ p ≥ 0, p ∈ P. (2d)
Master LP ( 2) is solved by column generation. Let π i , 1 ≤ i ≤ m, denote the dual variables of Constraints (2b), ν k + and ν k -, k ∈ K, are the dual variables of Constraints (2c). The reduced cost of an arc a ∈ A is defined as:

ca = j∈M -1 (a) c j - m i=1 j∈M -1 (a) α ij π i . The reduced cost of a path p = (v 0 , a 1 , v 1 , . . . , a n-1 , v n-1 , a n , v n ) ∈ P k is: c(p) = n j=1 caj -ν k + -ν k -.
So, the pricing subproblems correspond to finding, for each k ∈ K, a path p ∈ P k with minimum reduced cost.

Example of a Basic Model

For the sake of illustration, we show how to model the Capacitated VRP (CVRP) in the proposed solver. The problem is formally described as: Goal: Find a minimum cost set of routes, starting and ending at the depot, visiting all customers and such that the sum of the demands of the customers in a route does not exceed vehicle capacity.

Data: Undirected graph G = (V, E), V = {0, . . . , n}, 0 is the depot, V + = {1, . . . ,
A possible model is:

Model: A single graph G = (V, A) (the index k is omitted), A = {(i, j), (j, i) : {i, j} ∈ E}, v source = v sink = 0; R = R M = {1}; q a,1 = (d i + d j )/2, a = (i, j) ∈ A (define d 0 = 0); l i,1 = 0, u i,1 = Q, i ∈ V . Integer variables x e , e ∈ E.
The formulation is:

Min e∈E c e x e (3a)
S.t.

e∈δ ( i) x e = 2, i ∈ V + ; (3b) M (x e ) = {(i, j), (j, i)}, e = {i, j} ∈ E; L = n i=1 d i /Q , U = n.
The variables x are explicitly defined by the user and indicate how many times each edge is traversed. The objective function (3a) directly defines (1a) and Constraints (3b), degree constraints over x variables, directly define (1b). Constraints (1c) are implicitly defined by the mapping M and Constraints (1d) are implicitly defined by the values of L and U . Variables λ are implicitly defined by the set of resource constrained paths in G;

-1-0, 0-1-2-0, 0-1-2-1-0, 0-1-3-0, 0-1-4-0, 0-2-0, 0-2-1-0, 0-2-3-0, 0-3-0, 0-3-1-0, 0-3-2-0, 0-4-
x 01 = 2λ 1 + λ 2 + 2λ 3 + λ 4 + λ 5 + λ 7 + λ 10 + λ 13 , x 02 = λ 2 + 2λ 6 + λ 7 + λ 8 + λ 11 , x 03 = λ 4 + λ 8 + 2λ 9 + λ 10 + λ 11 ,
x 04 = λ 5 + 2λ 12 + λ 13 ,

x 12 = λ 2 + 2λ 3 + λ 7 ,
x 13 = λ 4 + λ 10 , x 14 = λ 5 + λ 13 ,

x 23 = λ 8 + λ 11 , x 24 = x 34 = 0, 2 ≤ λ 1 + λ 2 + λ 3 + • • • + λ 13 ≤ 4,
x ∈ Z

2λ 12 + 2λ 13 = 2 2 ≤ λ 1 + λ 2 + λ 3 + • • • + λ 13 ≤ 4,
λ ≥ 0 Some remarks on the CVRP model:

• It can be seen that the mapping transforms degree constraints (3b) defined over the x variables into "set-partitioning" constraints over the λ variables.

• Set P contains pairs of equivalent paths, visiting the same customers in reversed order. Those pairs lead to identical λ variables in the LP. For example, variables λ 2 and λ 7 in the above LP are identical. This redundancy does not harm the column generation. If one of those variables is included in the restricted Master LP, the other will not have negative reduced cost. Actually, in Section 5.1.1 it will be seen that the pricing algorithm may even profit from that symmetry in the model.

• Set P contains non-elementary paths. In the example, path 0-1-2-1-0 corresponds to a route that visits customer 1 twice. Although variables associated to non-elementary paths never appear in integer solutions (note that λ 3 has coefficient 4 in the transformed degree constraint of customer 1), they may weaken the linear relaxation. Section 3.2 will present a mechanism that can be used for eliminating the most undesirable non-elementary paths in that model. 3 Generalizing State-of-the-Art Components Formulation (1) can be used to model most VRP variants and also some other non-VRP applications. It can be solved by a standard BP algorithm (or a standard robust BCP algorithm [START_REF] Poggi De Aragão | Integer program reformulation for robust branch-andcut-and-price[END_REF], if (1b) contains separated constraints), where the RCSP subproblems are handled by a labeling dynamic programming algorithm. However, its performance would be very poor when compared to the best existing specific algorithms. One of the main contributions of this work is a generalization of the key additional concepts found in those state-of-the-art algorithms, leading to the construction of a powerful and still quite generic BCP algorithm. In order to do that, the new concepts of packing sets and elementarity sets are introduced.

Packing Sets and Elementarity Sets

Let P ⊂ 2 A be a collection of mutually disjoint subsets of A such that the constraints:

p∈P a∈S h p a λ p ≤ 1, S ∈ P, (4)
are satisfied by at least one optimal solution (x * , y * , λ * ) of Formulation [START_REF] Lemon | Library for Efficient Modeling and Optimization in Networks[END_REF]. In other words, the arcs in each S ∈ P appear at most once in all paths that are part of some optimal solution. In those conditions, we say that each element of P is a packing set. Note that a packing set can contain arcs from different graphs and not all arcs in A need to belong to some packing set. The definition of a proper collection P is application specific and part of the modeling task. It does not follow automatically from the analysis of Formulation [START_REF] Lemon | Library for Efficient Modeling and Optimization in Networks[END_REF].

In many applications the packing sets are more naturally defined on vertices, so we also provide that modeling alternative. Let coefficient h p v indicate how many times vertex v appears in a path p. Let V = V \ {v k source , v k sink : k ∈ K} be the set of vertices that are neither source or sink of their graphs. Let P V ⊂ 2 V be a collection of mutually disjoint subsets of V such that the constraints:

p∈P v∈S h p v λ p ≤ 1, S ∈ P V , (5)
are satisfied by at least one optimal solution (x * , y * , λ * ) of Formulation [START_REF] Lemon | Library for Efficient Modeling and Optimization in Networks[END_REF]. In those conditions, we say that the elements of P V are packing sets on vertices. Actually, in Section 5.1.1 we show that in some symmetric problems there is a computational advantage in defining packing sets on vertices.

Let E ⊂ 2 A be a collection of mutually disjoint subsets of A such that the constraints:

a∈S h p a λ p ≤ 1, S ∈ E, p ∈ P, (6)
are satisfied by at least one optimal solution (x * , y * , λ * ) of Formulation [START_REF] Lemon | Library for Efficient Modeling and Optimization in Networks[END_REF]. In other words, the arcs in each S ∈ E appear at most once in each path that is part of some optimal solution. In those conditions, we say that each element of E is an elementarity set. Collections E V ⊂ 2 V under similar conditions define elementarity sets on vertices.

The definition of a proper collection E is also part of the modeling task. However, any collection of packing sets P is also a collection of elementarity sets. Similarly, a collection of packing sets on vertices P V is also a collection of elementarity sets on vertices. So, it is always possible to define E = P or E V = P V . In fact, this is done in most models, including 7 out the 8 model examples given in Section 4. However, there are some models where it is advantageous to define a distinct collection of elementarity sets. The current implementation of the solver assumes that if the packing sets are defined in arcs, the elementarity sets are also defined in arcs and E ⊇ P. Similarly, if the packing sets are being defined on vertices then E V ⊇ P V . So, the collection of elementary sets is either equal or is an extension of the collection of packing sets.

The following concepts -ng-paths, Limited Memory Rank-1 Cuts, path enumeration, accumulated consumption branching, and rounded capacity cuts -were originally proposed and used on the most classical VRP variants, often CVRP and VRPTW. In our proposed generalization, those problems will correspond to simple models where the sets in P V and E V are the singletons formed by each customer vertex. However, as will be seen in Section 4, there are models where those sets are significantly more complex.

ng-paths

When modeling classical VRPs, one of the weaknesses of linear relaxation (2) is often the existence of non-elementary paths in P that can not be part of any integer solution. In those cases, one would like to eliminate all those paths from the definition of P . However, this would make the pricing subproblems much harder, to the point of becoming intractable in many cases. A good compromise between formulation strength and pricing difficulty can be obtained by the so-called ng-paths, introduced in Baldacci et al. [START_REF] Baldacci | New route relaxation and pricing strategies for the vehicle routing problem[END_REF].

In our more general context, we say that a path is E-elementary if it does not use more than once arcs in the same elementarity set of E. Let P k elem be the subset of the paths in P k that are E-elementary, P elem = ∪ k∈K P k elem . For example, suppose that V = {0, 1, 2, 3} with v source = v sink = 0, and E = {{(1, 2), (2, 1)}, {(1, 3), (3, 1)}, {(0, 3), (2, 3), (3, 2)}}. In this case, the path 0-1-2-3-1-0 would be E-elementary, but the paths 0-1-2-1-0 and 0-1-2-3-1-2-0 would not.

Ideally, we would like to price only E-elementary paths. Instead, we settle for generalized E-ng-paths defined as follows. For each arc a ∈ A, let N G(a) ⊆ E denote the ng-set of a. A E-ngpath may use two arcs belonging to the same elementarity set S, but only if the subpath between those two arcs passes by an arc a such that S / ∈ N G(a). In the previous example, suppose

that N G((1, 2)) = N G((2, 1)) = N G((1, 3)) = N G((3, 1)) = {{(1, 2), (2, 1)}, {(1, 3), (3, 1)}}, N G((2, 3)) = N G((3, 2)) = {{(1,
3), (3, 1)}}, and all other arcs have empty ng-sets. In this case, the path 0-1-2-3-1-2-0 would be allowed because {(1, 2), (2, 1)} does not belong to N G((2, 3)).

The ng-sets may be determined a priori; but also dynamically, like in [START_REF] Roberti | Dynamic ng-path relaxation for the delivery man problem[END_REF] and [START_REF] Bulhoes | A branch-and-price algorithm for the minimum latency problem[END_REF].

If the elementarity sets are being defined on vertices, there is the similar concept of E Velementary path: a path that does not use more than once vertices in the same elementarity set of E V . We also denote by P k elem the subset of the paths in P k that are E V -elementary. In this context, for each vertex v ∈ V , let N G(v) ⊆ E V be the ng-set of v. A E V -ng-path may use two vertices belonging to the same elementarity set S, but only if the subpath between those two vertices passes by a vertex v such that S / ∈ N G(v).

When E or E V are clear from the context, we may still refer to E-ng-paths or E V -ng-paths simply as ng-paths.

Limited Memory Rank-1 Cuts

The Rank-1 Cuts (R1Cs) [START_REF] Petersen | Chvátal-gomory rank-1 cuts used in a dantzig-wolfe decomposition of the vehicle routing problem with time windows[END_REF][START_REF] Pecin | Limited memory rank-1 cuts for vehicle routing problems[END_REF][START_REF] Bulhões | On the complete set packing and set partitioning polytopes: Properties and rank 1 facets[END_REF] are a generalization of the Subset Row Cuts proposed by Jepsen et al. [START_REF] Jepsen | Subset-row inequalities applied to the vehicle-routing problem with time windows[END_REF]. Here, they are further generalized as follows. Consider a collection of packing sets P. A Chvátal-Gomory rounding of Constraints (4), using a non-negative multiplier ρ B for each packing set S ∈ P, yields:

p∈P S∈P ρ S a∈S h p a λ p ≤ S∈P ρ S . (7)
Those R1Cs are potentially very strong, but each added cut makes the pricing subproblems significantly harder.

The limited memory technique [START_REF] Pecin | Improved branch-cut-and-price for capacitated vehicle routing[END_REF] is essential for mitigating that negative impact. In this technique, a R1C, characterized by its vector of multipliers ρ, is associated to a memory arc-set A(ρ) ⊆ A. The limited-memory R1C (lm-R1C) is defined as:

p∈P α(ρ, A(ρ), p)λ p ≤ S∈P ρ S , (8)
where the coefficient α(ρ, A(ρ), p) is computed as in the pseudo-code that describes Function α.

Function α(ρ, A, p = (v 0 , a 1 , v 1 , . . . , a n-1 , v n-1 , a n , v n )) 1 α ← 0, s ← 0; 2 for j = 1 to n do 3 if a j / ∈ A(ρ) then 4 s ← 0; 5 for S ∈ P do 6 if a j ∈ S then 7 s ← s + ρ S ; 8 if s ≥ 1 then 9 s ← s -1, α ← α + 1;
10 return α;

If A(ρ) = A, constraints ( 7) and ( 8) are identical. Otherwise, variables λ p corresponding to paths p passing by arcs a / ∈ A(ρ) may have their coefficients decreased. However, if the memory sets are adjusted in such a way that variables λ p with positive values in the current linear relaxation have the same coefficients that they would have in [START_REF] Baldacci | An exact algorithm for the pickup and delivery problem with time windows[END_REF], the resulting lm-R1C is as effective as the original R1C. Yet, if the final A(ρ) is a small subset of A, as usually happens, the impact in the pricing is much reduced.

If the model defines its packing sets in vertices, the R1Cs are defined in a similar way. There is a non-negative multiplier ρ S for each S ∈ P V and the cut is:

p∈P S∈P V ρ S v∈S V h p v λ p ≤ S∈P V ρ S . (9)
Given a memory arc-set A(ρ) ⊆ A corresponding to the vector ρ, the lm-R1C is defined as:

p∈P α(ρ, A(ρ), p)λ p ≤ S∈P V ρ S , (10)
where Function α is the same, except that the line 5 is replaced by (for S ∈ P V ) and the condition of line 6 by (v j ∈ S)

Regardless of if the packing sets are being defined on arcs or on vertices, it is possible to use lm-R1Cs where the memories are defined by vertex-sets. In this case, a memory vertex-set V (ρ) ⊆ V should be assigned to the lm-R1C corresponding to vector ρ. Function α should receive V (ρ) instead of A(ρ) as parameter and the condition in line 3 should be changed to

(v j / ∈ V (ρ)).
As discussed in [START_REF] Pessoa | Enhanced branch-cut-and-price algorithm for heterogeneous fleet vehicle routing problems[END_REF] and [START_REF] Sadykov | A bucket graph based labeling algorithm with application to vehicle routing[END_REF], memory vertex-sets perform better for most instances of some classical VRPs. This happens because R1C memory adjustment converges in less iterations in that case. In the other hand, memory arc-sets may be better for some harder instances because they allow for a finer memory adjustment, leading to less impact in the pricing.

Path Enumeration

The path enumeration technique was proposed by Baldacci et. al. [START_REF] Baldacci | An exact algorithm for the vehicle routing problem based on the set partitioning formulation with additional cuts[END_REF], and later improved by Contardo and Martinelli [START_REF] Contardo | A new exact algorithm for the multi-depot vehicle routing problem under capacity and route length constraints[END_REF]. It consists in trying to enumerate into a pool all paths in a certain set P k that can possibly be part of an improving solution. After a successful enumeration, the corresponding pricing subproblem k can be solved by inspection, saving time. Moreover, standard fixing by reduced costs can be used to remove paths from the pools. If the enumeration has already succeeded for all k ∈ K and once the total number of paths in the tables is reduced to a reasonably small number (say, less than 10,000), the formulation restricted to those paths can be given and directly solved by a general MIP solver.

In our context, we try to enumerate all paths p ∈ P k elem such that c(p) < U B -LB, where U B is the best known integer solution cost, and LB the value of the current linear relaxation. Moreover, if two paths p and p in P k map to variables λ p and λ p with identical coefficients in the essential constraints in (2b), the one with a larger cost is dominated and can be dropped. The essential constraints are those that are required to make the formulation valid, constraints in (2b) added only to strengthen the linear relaxation are not essential.

However, the enumeration procedure would be highly inefficient if the dominance could only be checked for pairs of complete paths. Instead, it is necessary to perform dominance over the partial paths (E-elementary paths starting at the source vertex) that are being constructed along the procedure. Our procedure uses the following dominance rule: if p and p are partial paths ending at the same vertex and having already visited exactly the same elementarity sets in E (regardless of the visitation order), the one with larger cost (breaking ties arbitrarily) is considered dominated and dropped. No complete path in P k elem that is the completion of a dominated partial path will be produced. The following condition is sufficient to assure that such enumeration procedure is valid (i.e., no improving solution is ever missed):

Sufficient Condition for Enumeration. Every two feasible partial E-elementary paths starting in v k source that end in the same vertex and map to different coefficients in some essential constraint in (2b) should have visited different subsets of E.

We remark that the above condition can not be checked automatically. In fact, in general it is not even possible to automatically determine what are the essential constraints in (2b). It is up to the modeler to prove that the provided model satisfies the sufficient condition, so enumeration can be used. Happily, in many models (including all the examples in Section 4) it is easy to prove that the condition is satisfied. However, if it is not satisfied, it is up to the modeler to prove that the enumeration is valid for his model directly from the dominance rule. Otherwise, the enumeration should be turned off.

Branching

Branching over individual x and y variables (or over constraints defined over those variables) is simple and do not change the structure of the pricing subproblems. In many models this kind of branching is sufficient for correctness. However, there are models where Constraints (1e) need to be explicitly enforced. However, branching over individual λ variables should be avoided due to a big negative impact in the pricing and also due to highly unbalanced trees [START_REF] Vanderbeck | Reformulation and decomposition of integer programs[END_REF]. The model offers two ways of branching over sets of λ variables, both are based on packing sets:

• Choose distinct sets S and S in P. Let P (S, S ) ⊆ P be the subset of the paths that contain arcs in both S and S . The branching is over the value of p∈P (S,S ) λ p , either 0 or 1. This is a generalization of the Ryan and Foster branching rule [START_REF] Ryan | An integer programming approach to scheduling[END_REF]. It is still to be avoided if possible, because it makes the pricing harder. However, using that scheme leads to more balanced search trees.

• Choose S ∈ P, r ∈ R k M and a certain threshold value t * : in the left child make u a,r = t * , for all a ∈ S; in the right child make l a,r = t * , , for all a ∈ S. This branching over the accumulated consumption of a resource generalizes the strategy proposed by Gélinas et al. [START_REF] Gélinas | A new branching strategy for time constrained routing problems with application to backhauling[END_REF]. The branching is not likely to be complete, in the sense that some fractional λ solutions can not be eliminated by it. However, it does not increase the pricing difficulty and it may work well in practice, postponing (and even avoiding) the use of a branching that makes pricing harder.

Rounded Capacity Cut Separators

The Rounded Capacity Cuts (RCCs), first proposed for CVRP [START_REF] Laporte | A branch and bound algorithm for the capacitated vehicle routing problem[END_REF], are still useful on modern BCP algorithms for that problem and also for a number of other VRP variants. Moreover, a very good heuristic separation routine is available for it in CVRPSEP library [START_REF] Lysgaard | CVRPSEP: A package of separation routines for the capacitated vehicle routing problem[END_REF]. So, we decided to introduce the concept of RCC separator as a feature of our model.

The RCC separator can only be used if the packing sets are defined on vertices. For a vertex v ∈ V , define S(v) as the packing set of P V that contains v, S(v) = ∅ if v is not in any packing set. An RCC separator is defined by setting a capacity Q and a demand function d : P V ∪ ∅ → R + such that d(∅) = 0 and is valid if there exists an optimal solution (x * , y * , λ * ) of Formulation (1) such that:

1. n j=0 d(S(v j )) ≤ Q, for all p = (v 0 , a 1 , v 1 , . . . , a n-1 , v n-1 , a n , v n ) ∈ P with λ * p ≥ 1;
2. for all S ∈ P V such that d(S) > 0, the corresponding constraints in (5) should be satisfied with equality by (x * , y * , λ * ).

Again, it is up to the modeler to prove that the separator included in the model is valid.

Given a valid RCC separator, for S ⊆ P V , d(S) denotes S∈S d(S) and h p S is the number of times that an arc in path p ∈ P enters in S. We say that an arc (v j-1 , v j ) enters in S if S(v j-1 ) / ∈ S and S(v j ) ∈ S. A Rounded Capacity Cut is the following valid inequality:

p∈P h p S λ p ≥ d(S) Q . ( 11
)
Cuts in format [START_REF] Balinski | On an integer program for a delivery problem[END_REF] are robust. The dual variable of the cut corresponding to an S ⊆ P V is simply subtracted from the reduced cost of all arcs entering S.

It is possible to define multiple RCC separators in the same model, each one having its demand function and capacity. This can be useful for modeling VRPs where routes are constrained by multiple dimensions. Packing sets would have zero demand in the dimensions that they do not "participate".

Model Examples

We selected a number of classical problems to illustrate the modeling capabilities of our solver. For each formulation, we only present explicitly the objective function and the general constraints corresponding to (1a-1b). The constraints corresponding to (1c) and (1d) are implicitly given by the mappings and by the bounds L k and U k , respectively.

Generalized Assignment Problem (GAP)

Data: Set T of tasks; set K of machines; capacity Q k , k ∈ K; assignment cost c k t and machine load w k t , t ∈ T , k ∈ K.
Goal: Find an assignment of tasks to machines such that the total load in each machine does not exceed its capacity, with minimum total cost.

Model: RCSP generator graphs G k = (V k , A k ) for each k ∈ K: V k = {v k t : t = 0, . . . , |T |}, A k = {a k t+ = (v k t-1 , v k t ), a k t-= (v k t-1 , v k t ) : t = 1, . . . , |T |}, v k source = v k 0 , v k sink = v k |T | (see Fig. 2); R k = R k M = {r k }; q a k t+ ,r k = w k t , q a k t-,r k = 0, t ∈ T ; [l v k t ,r k , u v k t ,r k ] = [0, Q k ], t ∈ T ∪ {0}. Integer variables x k t , t ∈ T , k ∈ K. The formulation is: Min t∈T k∈K c k t x k t (12a) S.t. k∈K x k t = 1, t ∈ T ; (12b) M (x k t ) = {a k t+ }, t ∈ T, k ∈ K; L k = 0, U k = 1, k ∈ K. P = E = ∪ t∈T {{a k t+ : k ∈ K}}.
Branching is over the x variables. Enumeration is on.

Comments: Graphs G k , illustrated in Figure 2, are designed to model binary knapsack constraints: each path in P k corresponds to a possible assignment of a set of tasks to machine k. A collection of packing sets (that are also the elementarity sets) is provided, so the features described in Section 3, that extend the basic formulation, can be used. In this model, the validity of the chosen P is a clear consequence of Constraints (12b) and of the mapping. However, in other problems, the validity of the packing sets provided by the modeler may not be obvious. All constraints in (12b) are essential. It can be checked that the enumeration sufficient condition is satisfied, so the enumeration procedure can be used. Goal: Find a packing using the minimum number of bins, such that, for each dimension, the total weight of the items in a bin does not exceed its capacity.

k 0 v k 1 v k 2 v k 3 v k |T |-1 v k |T | a k 1+ a k 1- a k 2+ a k 2- a k 3+ a k 3- a k |T |+ a k |T |- . . . . . .
Model: A single graph G = (V, A) (we omit the index k in such cases): V = {v t : t = 0, . . . , |T |}, A = {a t+ = (v t-1 , v t ), a t-= (v t-1 , v t ) : t = 1, . . . , |T |}, v source = v 0 , v sink = v |T | . R = R M = D; q at+,d = w d t , q at-,d = 0, t ∈ T, d ∈ D; [l vt,d , u vt,d ] = [0, Q d ], t ∈ T ∪ {0}, d ∈ D. Continuous variables x t , t ∈ T ∪ {0}. The formulation is: Min x 0 (13a) S.t. x t = 1, t ∈ T ; (13b) M (x 0 ) = {a 1+ , a 1-}, M (x t ) = {a t+ }, t ∈ T ; L = 0, U = ∞. P = E = ∪ t∈T {{a t+ }}.
Branching over accumulated resource consumption and, if still needed, by Ryan and Foster rule. Enumeration is on.

Comments: Defining the x variables as integer would be useless, it would not be possible to branch over them (except, in very limited way, over x 0 ). Branching on λ variables is needed; Ryan and Foster rule suffices for correctness, however accumulated resource consumption branching (that may not suffice) should be performed first.

Capacitated Vehicle Routing Problem (CVRP)

The problem was already defined in Section 2.

Model:

A single graph G = (V, A), A = {(i, j), (j, i) : {i, j} ∈ E}, v source = v sink = 0; R = R M = {1}; q a,1 = (d i + d j )/2, a = (i, j) ∈ A (define d 0 = 0); l i,1 = 0, u i,1 = Q, i ∈ V .
Integer variables x e , e ∈ E. The formulation is:

Min e∈E c e x e (14a)
S.t. e∈δ ( i) x e = 2, i ∈ V + , (14b)
x e ≤ 1, e ∈ E \ δ(0); (14c) M (x e ) = {(i, j), (j, i)}, e = {i, j} ∈ E; L = n i=1 d i /Q , U = n. P V = E V = ∪ i∈V+ {{i}}. RCC separator on (∪ i∈V+ {({i}, d i )}, Q).
Branching on x variables. Enumeration is on.

Comments: Constraints (14c) are separated (by inspection) as user cuts. The packing sets are defined on vertices. In this problem, defining the resource consumption in a symmetric way (q (i,j),1 = q (j,i),1 ) improves the efficiency of the pricing, as will be discussed in Section 5.1.1. As constraints (14c) are not essential, the enumeration condition over (14b) is satisfied. The function d : P V ∪ ∅ → R + for the RCC separator is defined as the set of all pairs (S, d(S)) for which d(S) = 0.

Heterogeneous Fleet Vehicle Routing Problem (HFVRP)

Data: Undirected graph G = (V, E), V = {0, . . . , n}, 0 is the depot, V + = {1, . . . , n} are the customers; positive demand d i , i ∈ V + ; set of vehicle types K = {1, . . . , m}; number of available vehicles u k , k ∈ K; edge costs c k e , e ∈ E, k ∈ K (assume that fixed costs f k for using a vehicle of type k are included in the cost of the edges incident to the depot); vehicle type capacity Q k , k ∈ K.

Goal: Find a minimum cost set of routes, each route associated to a vehicle type and starting and ending at the depot, visiting all customers and such that the sum of the demands of the customers in a route does not exceed its vehicle type capacity. The number of routes for a vehicle type should not exceed its availability.

Model:

Graphs G k = (V k , A k ), V k = {v k 0 , . . . , v k n }, A k = {(v k i , v k j ), (v k j , v k i ) : {i, j} ∈ E}, v k source = v k sink = v k 0 , k ∈ K; R k = R k M = {r k }; q a,1 = (d i + d j )/2, a = (v k i , v k j ) ∈ A k , k ∈ K (define d 0 = 0); l v k i ,r k = 0, u v k i ,r k = Q k , v k i ∈ V k , k ∈ K. Integer variables x k e , e ∈ E, k ∈ K. The formulation is: Min k∈K e∈E c k e x k e (15a) S.t. k∈K e∈δ ( i) x k e = 2, i ∈ V + ; (15b) M (x k e ) = {(v k i , v k j ), (v k j , v k i )}, e = {i, j} ∈ E, k ∈ K; L k = 0, U k = u k . P V = E V = ∪ i∈V+ {{v k i : k ∈ K}}. RCC separator on (∪ i∈V+ {({v k i : k ∈ K}, d i )}, max k∈K Q k ).
Branching first on the aggregation of x variables corresponding to number of times that a vehicle of each type is used, then on the aggregation of x variables corresponding to the assignment of a customer to a vehicle type or on the aggregation of x variables corresponding to the edges of the original graph G . Enumeration is on.

Team Orienteering Problem (TOP)

Data: Directed graph G = (V, A), V = {0, . . . , n + 1}, 0 and n + 1 are the initial and final depots, respectively, V + = {1, . . . , n} are the customers; positive travel time t a , a ∈ A; profit p i , i ∈ V + ; maximum route duration T ; and fleet size F . Goal: Find a set of at most F routes, each one starting at 0, ending at n + 1 and not exceeding the maximum route duration, that visit each customer at most once and maximize the total profit of the visited customers. Goal: Find a minimum cost set of routes, starting and ending at the depot, performing all pickups and deliveries within the time windows (waiting is allowed) and such that, along a route, the total demand already collected but not yet delivered does not exceed vehicle capacity.

Model: A single graph G = (V, A), v source = 0, v sink = n + 1; R = R M = {1}; q a,1 = t a , a = (i, j) ∈ A; l i,1 = 0, u i,1 = T, i ∈ V . Integer variables x a , a ∈ A and binary variables y i , i ∈ V + . The formulation is: Min -i∈V+ p i y i (16a) S.t. i∈δ -(i) x a = y i , i ∈ V + ; (16b) M (x a ) = {a}, a ∈ A; L = 0, U = F . P V = E V = ∪ i∈V+ {{i}}.
Model: A single graph G = (V, A), v source = v sink = 0. R M = {n+2} (the time resource); R N = {1, . . . , n + 1}, the first n resources are non-disposable;

q (v,v ),v = 1, if v ∈ P , q (v,v ),v -n = -1, if v ∈ D, and q (v,v ),n+1 = d v , (v, v ) ∈ A; q a,n+2 = t a , a ∈ A; all other resource consumptions are zero; u v,r = 1, v ∈ P ∪ D, r = 1, . . . , n, u v,n+1 = Q, v ∈ V , and (l v,n+2 , u v,n+2 ) = (l v , u v ), v ∈
V ; all other resource bounds are zero. Integer variables x a , a ∈ A. The formulation is:

Min a∈A c a x a (17a) S.t. a∈δ -(v) x a = 1, v ∈ P ; (17b) M (x a ) = {a}, a ∈ A; L = 0, U = ∞. P V = E V = ∪ v∈P {{v}}.
Branching first on the aggregation of x variables corresponding to number of routes, then on individual x variables.

Comments: This example illustrates the use of non-monotone secondary resources, some of them being defined as non-disposable. The first n resources are used to enforce that if a route performs a pickup i it can not return to the depot before doing the corresponding delivery i + n.

If those resources were disposable, routes that visit i + n without having visited first i would be possible. Resource n + 1 controls the capacity along the route.

Capacitated Arc Routing (CARP)

Data: Undirected graph G = (V , E), V = {0, . . . , n}, 0 is the depot vertex; positive cost c e and non-negative demand d e , e ∈ E, set F of required edges:

F = {e ∈ E | d e > 0}; vehicle capacity Q.
Goal: Find a minimum cost set of routes, closed walks starting and ending at the depot, that serve the demands in all required edges. Edges in a route can be traversed either serving or deadheading (not serving). The sum of the demands of the served edges in a route can not exceed capacity.

Model: For i, j ∈ V , let D(i, j) ⊆ E be the set of edges in a chosen cheapest path from i to j, with cost C(i, j) = e∈D(i,j) c e . Define a dummy required edge r 0 = (0, 0 ) and F 0 = F ∪ {r 0 }. Define an auxiliary undirected complete graph G = (F 0 , E ). For each r = (w 1 , w 2 ) ∈ F 0 , define o(r, w 1 ) = w 2 and o(r, w 2 ) = w 1 .

The model has a single RCSP graph generator

G = (V, A), V = {v w r : r ∈ F 0 , w ∈ r}, A = {(v w1 r1 , v z1 r2 ), (v w1 r1 , v z2 r2 ), (v w2 r1 , v z1 r2 ), (v w2 r1 , v z2 r2 ) : r 1 = (w 1 , w 2 ), r 2 = (z 1 , z 2 ) ∈ F 0 , r 1 = r 2 }, v source = v 0 r0 , v sink = v 0 r0 ; R = R M = {1}; for a = (v w r1 , v z r2 ) ∈ A, q a,1 = d r2 ; l v,1 = 0, u v,1 = Q, v ∈ V . Integer variables x a , a ∈ A. For a = (v w r1 , v z r2 ) ∈ A, c a = C(w, o(r 2 , z)) + c r2 . The formulation is: Min a∈A c a x a (18a)
S.t.

a∈δ -({v

w 1 r ,v w 2 r }) x a = 1, r = (w 1 , w 2 ) ∈ F, (18b)
plus Lifted Odd-Cutsets [START_REF] Belenguer | The capacitated arc routing problem: Valid inequalities and facets[END_REF][START_REF] Bartolini | Improved lower bounds and exact algorithm for the capacitated arc routing problem[END_REF];

M (x a ) = {a}, a ∈ A; L = 0, U = ∞. P V = E V = ∪ r=(w1,w2)∈F {{v w1 r , v w2 r }}. RCC separator on (∪ r=(w1,w2)∈F {({v w1 r , v w2 r }, d r )}, Q).
Branching first on the aggregation of x variables corresponding to node degrees in original graph G or on the aggregation of x variables corresponding to edges of graph G .

Comments: A case where a more complex transformation (similar to [START_REF] Longo | Solving capacitated arc routing problems using a transformation to the cvrp[END_REF][START_REF] Bartolini | Improved lower bounds and exact algorithm for the capacitated arc routing problem[END_REF][START_REF] Benavent | A branch-price-and-cut algorithm for the min-max k-vehicle windy rural postman problem[END_REF]) is used to fit the problem into the model. The separation of Lifted Odd-Cutsets is performed using the Gomory-Hu tree algorithm.

Black-and-White Traveling Salesman Problem (BWTSP)

Data: Undirected graph G = (V = B ∪ W, E), where B is the set of black vertices and W is the set of white vertices; positive cost c e , e ∈ E; positive distance d e , e ∈ E; limits Q and D.

Goal: Find a minimum cost tour such that every subpath between two consecutive black vertices visits at most Q white vertices and has total distance at most D.

Model: A single graph G = (B ∪ W ∪ {0}, A), A = {(i, j), (j, i) : (i, j) ∈ E} ∪ {(0, i), (i, 0) : i ∈ B}, v source = v sink = 0; R = R M = {1, 2}.
Consumptions of the first (capacity) and second (distance) resources for arc a = (i, j) ∈ A are

q a,1 = Q, i ∈ B and j ∈ B, 0.5 • |{i, j} ∩ W |, otherwise,
q a,2 = l e , (i, j) = e or (j, i) = e, 0, otherwise, where |{i, j} ∩ W | is the number of white vertices in {i, j}. We also have

l i,1 = 0, u i,1 = Q, l i,2 = 0, u i,2 = D, i ∈ V . Integer variables x e , e ∈ E.
The formulation is:

Min e∈E c e x e (19a)
S.t. e∈δ(i) x e = 2, i ∈ V, ( 19b
) i∈V j ∈V x (i,j) ≥ 2, ∅ = V ⊂ V ; (19c) M (x e ) = {(i, j), (j, i)}, e = {i, j} ∈ E; L = U = |B|. P V = ∪ i∈W {{i}}, E V = ∪ i∈B {{i}} ∪ P V . RCC separator on (∪ i∈W {({i}, 1)}, Q).
Branching on x variables. Enumeration is on. Constraints (19c) (standard TSP subtour elimination constraints) are added dynamically as core cuts.

Comments: The graph G is used for generating the subpaths linking two consecutive black vertices. This model demonstrates the case when the collection of packing sets is not equal to the collection of elementarity sets. Each black vertex i ∈ B defines an elementarity set, but not a packing set. In fact, we cannot put a black vertex to any packing set as it appears in two paths in every feasible solution. However, a black vertex appears at most once in every subpath, therefore it can be put to an elementarity set. The model is illustrated in Figure 3.

Implementation

Algorithms used in the implementation of the solver are generalizations of already published algorithms. Thus, we only give an overview of these algorithms; together with references to the original papers. In fact, the main goal of this section is to explain how modeling decisions and solver parameters impact the implemented algorithms.

Labeling Algorithm for Pricing Problems

Pricing problems are solved by a bi-directional labeling dynamic programming algorithm [START_REF] Righini | Symmetry helps: Bounded bi-directional dynamic programming for the elementary shortest path problem with resource constraints[END_REF], using the bucket graph implementation proposed in [START_REF] Sadykov | A bucket graph based labeling algorithm with application to vehicle routing[END_REF]. A label is a data structure that represents either a forward partial path started in vertex v k source or a backward partial path started in vertex v k sink . Initially, only labels representing null paths at v k source or at v k sink exist. Labels are extended along arcs in A k , if the new extended path is feasible, then the corresponding new label is created. The key feature of a labeling algorithm is the use of dominance checks. Let p and p be two partial paths ending at the same vertex. If it can be proved that any complete path that is a completion of p will have a reduced cost greater than the reduced cost of the path obtained by the same completion over p, then p is dominated by p and the corresponding label can be dropped. Labels are grouped and stored in buckets. The idea is to perform frequent dominance checks only for labels in the same bucket, in order to avoid spending an excessive time on those operations. In our algorithm, labels with the same final vertex and within the same ranges of accumulated consumption values of main resources are put in the same bucket. Such bucket organization has advantages over a simpler organization based on resource discretization, used for example in [START_REF] Pecin | Improved branch-cut-andprice for capacitated vehicle routing[END_REF]: it improves significantly the algorithm performance in the cases where resource intervals are given by very large integer numbers and even allows the use of continuous resource consumptions.

Our implementation supports at most two main resources, as we believe that having three or more main resources would result in too many buckets and too few labels per bucket, not being computationally advantageous. Accumulated resource consumption ranges for buckets are specified using the step size dr defined for each main resource r ∈ R k M . Each bucket contains labels with a main resource consumption within two consecutive multiples of the step size for this resource. Initial step sizes are defined using parameter ψ buck which designates the maximum number of buckets per vertex. Let Qk r be the maximum spread of the resource consumption for resource

r ∈ R k M in graph G k : Qr = u v k sink ,r -l v k source ,r
. Then in the case with one main resource, d1 = Q1 /ψ buck . In the case two main resources dr = Qr / ψ buck for r = 1, 2. In the course of the algorithm, step sizes may be automatically reduced (i.e. number of buckets per vertex may be increased) if it determined that too many dominance checks are performed between labels in a same bucket. If this happens, the buckets are recreated. Binary parameter ψ reduc determines whether automatic step size reduction is activated.

The bucket graph defines one node per bucket and directed arcs connecting pairs of buckets through which the labels can be possibly extended. Distinct bucket graphs are then defined for forward and backward labeling. The concept of bucket graph is useful because:

1. It helps to determine an efficient order of treatment for the buckets. If the bucket graph is acyclic, it is desirable to process the buckets in its topological order because no further extension from a bucket is necessary after it has been processed. If the bucket graph contains cycles, then buckets are handled in the topological order of its strongly connected components, trying to minimize such reprocessing.

2. It is used to improve the efficiency of performing dominance checks between labels in different buckets.

3. Arcs can be removed from the bucket graph by reduced cost arguments, avoiding label extensions in future calls to the pricing. This fixing procedure is more powerful than the one in [START_REF] Irnich | Path-reduced costs for eliminating arcs in routing and scheduling[END_REF]. Binary parameter φ elimin determines whether this more sophisticated fixing procedure is applied.

Our labelling algorithm and the bucket arc elimination procedure are bi-directional. For the first main resource, we calculate a special resource consumption value q * equal to the average middle value of the resource consumption bounds in all vertices. In the forward (backward) labelling, we keep only labels with the first main resource consumption not larger (larger) than q * . Then the concatenation step is performed to obtain complete paths. Completion bounds are used to speed up concatenation. If the difference between the number of forward and backward non-dominated labels is large, the next call of the pricing will use an automatically adjusted value of q * . Parameter φ bidir defines whether bidirectional search is used. If φ bidir = 0, bidirectional search is never used. If φ bidir = 1, bidirectional search is always used. If φ bidir = 2, bidirectional search is used only for the exact labeling, and not used for the heuristic labeling presented in Section 5.1.3.

Benefiting from symmetry

We exploit forward-backward path symmetry when solving the pricing problem for graph G k , k ∈ K, if the following conditions are satisfied.

• Packing and elementarity sets are defined on vertices.

• For each arc a = (i, j) ∈ A k , k ∈ K, there exists arc a = (j, i) ∈ A k with the same resource consumptions and the same inverse mapping: q a,r = q a ,r , for all r ∈ R k , and M -1 (a) = M -1 (a ). For the sake of symmetry detection, vertices v k source and v k sink are treated as the same vertex, even if the user set them as being distinct vertices.

• Resource consumption bounds are the same for all vertices:

[l v,r , u v,r ] = [l v ,r , u v ,r ] = [0, Q r ]
, for all r ∈ R and for all pairs v, v ∈ V k .

Under these conditions, value q * is fixed to Q 1 /2. The backward labelling is not executed, this saves a significant time. The concatenation is performed just after the forward labelling, using symmetric copies of forward labels as backward labels, as described in [START_REF] Sadykov | A bucket graph based labeling algorithm with application to vehicle routing[END_REF].

Dominance on disposable and non-disposable resources

When performing dominance checks between labels corresponding to partial paths p and p ending at the same vertex, we need to compare their accumulated resource comsumptions. Path p can only dominate p if, for every disposable resource r ∈ R k , the accumulated consumption of p is not larger than the accumulated consumption of p . Note that a strictly smaller accumulated consumption does not prevent the dominance, because it is always possible to dispose the additional units of that resource in a completion of p. This is not true for non-disposable resources. Therefore, path p can only dominate p if, for every non-disposable resource r ∈ R k , the accumulated consumption of p is identical to the accumulated consumption of p . This means that models should only define resources as non-disposable if this is really necessary. Remark that monotone and non-monotone resources are not treated differently in the dominance.

A resource is binary if its accumulated consumption can only be 0 or 1. Sets of secondary binary resources of the same type (disposable or not disposable) are implemented in a special way, being represented as bitsets in the labels. This greatly decreases the time spent for dominance checks between labels.

Pricing heuristics

We use two labelling heuristics similar to [START_REF] Sadykov | A bucket graph based labeling algorithm with application to vehicle routing[END_REF]. In the lighter heuristic, only one label per bucket is kept, the one with the smallest reduced costs. In the heavier heuristic, only reduced cost and resource consumption is used to check dominance between labels. Label dimensions related to ng-paths and rank-1 cuts are ignored.

Column and cut generation

We use three-stage column generation. In the first and second stages, the lighter and the heavier labelling heuristics are used, and at most γ heur columns are generated per iteration and per k ∈ K. In the last stage, the exact labeling algorithm is used, and at most γ exact columns are generated per iteration and per k ∈ K. We use automatic dual price smoothing [START_REF] Pessoa | Automation and combination of linear-programming based stabilization techniques in column generation[END_REF] technique to improve the column generation convergence in each stage. Parameter σ stab specifies the minimum column generation stage in which stabilization is active. The column clean-up procedure is also used: each time the total number of columns exceeds 10,000, only 66% of the columns with the smallest reduced cost remain in the master problem while other columns are removed.

If enabled, the bucket arc elimination procedure for each pricing problem k ∈ K is performed after the first column generation convergence and each time the current primal-dual gap is reduced by more than 10% since the last call to it. Immediately after this procedure, the bidirectional path enumeration labeling algorithm for graph G k is executed. This algorithm is a generalization of the one used in [START_REF] Baldacci | An exact algorithm for the vehicle routing problem based on the set partitioning formulation with additional cuts[END_REF][START_REF] Pecin | Improved branch-cut-andprice for capacitated vehicle routing[END_REF]. It is aborted if the number of non-dominated labels (in this context they correspond to partial E-elementary paths) exceeds ω labels or the number of generated paths exceeds ω routes . If the enumeration algorithm for graph G k succeeds, i) the enumerated paths are stored in a pool ii) all columns corresponding to paths that are not Eelementary are deleted from the master, iii) in future column generation iterations the pricing is performed by inspection in the pool, and iv) the fixing by reduced costs also starts to be performed in the pool and removes paths from it. If the total number of paths for all graphs drops below ω MIP , all corresponding columns are added to the restricted master, and current node is solved by the MIP solver. Initial ng-sets may be explicitly given by the user. However, it is possible (and usually much more practical) to specify a distance matrix between elementarity sets and a parameter η init , which was always equal to 8 for the experiments reported in this paper. If elementarity sets are defined on arcs, N G(a) will include the η init closest elementarity sets to the elementarity set of a according to the distance matrix. If elementarity sets are defined on vertices, N G(a) will include the union of the η init closest elementarity sets to the elementarity sets of both extremities of a.

We have implemented a dynamic extension of ng-sets [START_REF] Roberti | Dynamic ng-path relaxation for the delivery man problem[END_REF][START_REF] Bulhoes | A branch-and-price algorithm for the minimum latency problem[END_REF], which is useful for some problems. Given a fractional solution λ obtained after column generation convergence, we augment ng-sets based on "elementarity set cycles" in paths p ∈ P k such that λp > 0. An elementarity set cycle is a partial path started and finished by arcs belonging to a same elementarity set S. For each considered path, we find all cycles of size at most 5 or, if there are no such cycles, a minimum size cycle. Then we add S to N G(a) for every arc in these cycles. At most 100 paths p with cycles are considered, those with the largest values λp . The size of ng-sets is limited by parameter η max . All columns corresponding to paths which are not feasible with respect to new ng-sets are deleted from the restricted master.

In each cut round, we generate at most 100 rounded capacity cuts (if RCC separators are defined). We call a rank-1 cut l-row if the corresponding vector ρ of multipliers has l positive components. In each round, we generate at most θ num l-row rank-1 cuts for each l ∈ {1, 4, 5, 6, 7}, and at most 1.5 • θ num 3-row rank-1 cuts (which are subset-row cuts [START_REF] Jepsen | Subset-row inequalities applied to the vehicle-routing problem with time windows[END_REF]). Parameter θ rows determines the maximum number of rows in rank-1 cuts, i.e. if l > θ rows , l-row rank-1 cuts are not separated. l-row rank-1 cuts are separated by exhaustive enumeration for l = {1, 3, 4}. For separating other rank-1 cuts, a local search heuristic is used. It is based on the distance matrix between packing sets (contained in the distance matrix between elementarity sets) given by the user. The local search heuristic separates only cuts corresponding to multipliers ρ such that every packing set S, ρ S > 0, is among 16 closest to any other packing set S , ρ S > 0.

There is parameter θ mem to determine how rank-1 cuts limited memory is computed. If θ mem = 1, the arc memory is used as described in Section 3.3. If θ mem = 2, the vertex memory is used, which is "projected" to arcs between vertices in the memory before executing the labeling algorithm. Vertex memory allows the algorithm to converge faster, as rank-1 cuts are stronger. However, the impact on the running time of the labeling algorithm is larger. Thus, the obtained dual bounds may be weaker due to the time thresholds τ soft and τ hard defined below. If θ mem = 0, the root node is solved two times, first time with arc memory for rank-1 cuts, and the second time with vertex memory. After solving the root, it is automatically determined based on the latest pricing time and the dual bound obtained which memory is used for the remaining of the branch-and-bound tree.

Increasing of ng-sets can also be considered as cut generation, as this procedure improves dual bound given by the master problem. In the root node, ng-set augmentation has the highest priority, robust cuts (RCCs and user cuts) have medium priority, and rank-1 cuts have the smallest priority. We separate cuts with lower priority only if the tailing-off condition is satisfied for cuts of higher priority. The tailing condition is fulfilled when the dual gap is reduced by less than δ gap % after δ num round of cuts. In the other nodes of the branch-and-bound tree, priority of all cuts is the same, i.e. they all are separated in every round.

There are two thresholds τ soft and τ hard for the running time of the exact labeling algorithm. If the running time for at least one pricing problem exceeds τ soft , after the column generation convergence, cuts are not separated and branching is performed. If the running time exceeds τ hard , the labelling algorithm is interrupted and the rollback procedure [START_REF] Pecin | Improved branch-cut-andprice for capacitated vehicle routing[END_REF] is executed: the cuts added in the last separation round are deleted from the master, and branching is performed.

Safe dual bounds

An optimal dual solution returned by the restricted master problem is actually an approximation of a real one, as competitive LP solvers use floating-point arithmetic. Thus the sign of the minimum reduced cost found by the pricing problem can hardly be decided. This can lead to premature termination or to endless loops [START_REF] Held | Maximum-weight stable sets and safe lower bounds for graph coloring[END_REF]. For most applications, no practical difficulty occurs. However for some problems, wrong dual bound can be calculated. Vertex coloring and bin packing are examples of such problems.

Our solver has an option to compute a numerically safe dual bound similarly to the method proposed in [START_REF] Held | Maximum-weight stable sets and safe lower bounds for graph coloring[END_REF]. This approach can be applied in our solver if i) variables y are absent (i.e. n 2 = 0), and ii) all coefficients α are non-negative.

We use the property that a valid (Lagrangian) dual bound can be computed for any vector of dual values. Given as a parameter a large integer constant K, we call the pricing problem with modified dual values Kπ i , i = 1, . . . m, Kν k + , Kν k -, k ∈ K. The dual values corresponding to the cuts are modified in the same way. The modified cost of variable x j in the pricing problems is equal to Kc j , j = 1, . . . , n 1 . This rounding procedure allows us to obtain a feasible dual solution at the end of column generation convergence under two conditions given in the previous paragraph. As the objective function of any modified pricing problem is integer, its solution value can be computed exactly. The (Lagrangian) dual bound is then calculated in the usual way and then multiplied by K-1 . When the option to use safe dual bounds is activated ( K > 0), the reduced cost tolerance of the LP solver is set to the minimum possible value (10 -9 for CPLEX).

Branching

The user has the possibility to set the priority for each branching strategy he uses for his application. The selection among branching candidates with the same priority is done using a sophisticated hierarchical evaluation strategy similar to the one proposed in [START_REF] Røpke | Branching decisions in branch-and-cut-and-price algorithms for vehicle routing problems[END_REF][START_REF] Pecin | Improved branch-cut-andprice for capacitated vehicle routing[END_REF]. The idea is to spend more time evaluating branching candidates in the lowest levels of the branch-andbound tree where each selection has a greater impact on the overall time, and spend less time as the level increases, taking advantage of the history of previous evaluations. The following three evaluation phases are used: Phase 0: Up to a half of the candidates are chosen from history using pseudo-costs (if history is not empty). The remaining candidates are chosen in a balanced way between all branching strategies of the current priority. Within the same strategy, the candidates are chosen based on the distance from its fractional value to the closest integer, the larger distance is better.

Phase 1: Evaluate the selected candidates from phase 0 by solving the current restricted master LP modified for each created child node, without generating columns. Select the variables with the maximum value of ∆LB 1 ×∆LB 2 , where ∆LB i denotes the increase in the current lower bound obtained for the ith child node, for i = 1, 2 (Product Rule, [START_REF] Achterberg | Constraint Integer Programming[END_REF]).

Phase 2: Evaluate the selected candidates from phase 1 by solving the relaxation associated to each created child node, including column generation with heuristic labelling for solving the pricing problem. Cut generation is not applied in this phase. The best candidate is also selected by the Product Rule.

The maximum number of candidates evaluated in phase ρ = 1, 2 in branch-and-bound node bbn is equal to min{ζ num ρ , T S(n bb ) • ζ estim ρ }, where T S(n bb ) is an estimation of the size of the subtree rooted at the parent of node bbn (it is equal to infinity for the root). Calculation of estimation T S(bbn) follows that of [START_REF] Kullmann | Handbook of Satisfiability, chapter Fundaments of branching heuristics[END_REF]. See also [START_REF] Le | An abstract model for branching and its application to mixed integer programming[END_REF] for a related work.

Primal heuristics

Restricted master and diving heuristics [START_REF] Sadykov | Primal heuristics for branch-and-price: the assets of diving methods[END_REF] are built-in to the solver to improve the primal solution during the search. When used, these heuristics are executed at each branch-and-bound node before branching. The first heuristic uses a MIP solver to solve the current restricted master problem as a mixed-integer program with the time limit χ rm . The second heuristic is the diving heuristic with Limited Discrepancy Search (LDS). It is applied in branch-and-bound nodes of depth at most χ div . The diving heuristic uses two parameters χ disc and χ depth which correspond to parameters maxDiscrepancy and maxDepth in [START_REF] Sadykov | Primal heuristics for branch-and-price: the assets of diving methods[END_REF]. In each dive of the diving heuristic, we iteratively fix a column with largest value in the fractional solution to one and then we solve the resulting residual master problem with column generation. In the diving heuristic with LDS, one does several dives. Thus, a search tree is formed in which backtracking is allowed. In each node n div of the diving search tree, we keep a tabu list of columns forbidden to be fixed. Tabu list of a child node is initialized with the current tabu list of the parent node (the initial tabu list at the root is empty). After a backtrack to a node n div , we insert in its tabu list the column which was fixed by the previous fixing decision in n div . An additional fixing decision in n div is allowed (i.e. other child nodes of n div can be created) if the size of its tabu list does not exceed χ disc and the depth of n div is not larger than χ depth .

Parameterization

Based on the description above, we list in Table 1 the solver parameters available to the user. In addition the user may provide the following information to improve the solver performance.

• Designation of the first main resource which will be used for the bi-directional labelling.

• Priorities for branching strategies.

• Distance matrix between elementarity sets for defining initial ng-sets and for local search heuristic separation of l-row rank-1 cuts with l ≥ 5.

We now give recommendations for tuning the parameters. The time thresholds τ soft and τ hard may be reduced if these thresholds are attained and the branch-and-bound tree is large. Parameter ψ buck may be increased if the pricing time is large during the first column generation convergence. Bucket arc elimination may be switched off (φ elim = 0) if it takes a large part of the solution time. Number of generated columns per iteration γ exact , γ heur may be reduced if set K of graphs is relatively large. Enumeration parameters ω labels , ω routes may be reduced if enumeration takes a large part of solution time (this often happens for easy instances). They may be increased for hard instances with small primal-dual gap. Parameter ω MIP may be decreased if enumerated MIP takes a lot of time to be solved by the MIP solver (and increased in the opposite case). Sometimes ω MIP may be set to zero if separation of lazy constraints is used: in Description Notation Default value(s) Time thresholds for the labeling algorithm τ soft , τ hard 10 sec., 20 sec. Calculation of step sizes for buckets ψ buck , ψ reduc 25, 1 (on) Bi-directional search and bucket arc elimination φ bidir , φ elim 2, 1 (on) Max. # of generated columns per iteration γ exact , γ heur 150, 30 Minimum stage for stabilization σ stab 0 (everywhere) Max. # of labels and paths in the enumeration ω labels , ω routes 10 6 , 10 χ div , χ depth , χ disc -1 (off), 0, 0 Table 1: Parameters of the solver available to the user and their default values this case the enumerated MIP may be resolved many times. Parameters η init and η max should be set to zero if the elementarity constraints are always satisfied, f.e. when the graphs are acyclic. Dynamic ng-relaxation should be used (η init should be set to zero) if the distance matrix between elementarity sets is not specified. Parameters η init and η max may be changed according to the difficulty of the pricing problem.

Limited-memory rank-1 cuts parameters may be quite important for the performance. For instances with long paths, arc memory should be used (θ mem should be set to 1). For a mix of instances with short and long paths, automatic memory should be used (θ mem should be set to 0). For instances with relatively short routes and a large branch-and-bound tree, parameter θ rows may be increased to 7. Number of generated rank-1 cuts θ num may be reduced if the pricing time increases significantly after one cut generation round.

Cut generation tailing-off parameters δ gap , δ num may be changed depending on the size of the branch-and-bound tree. In the case of a large tree or a large cut separation time, tailing-off may be reduced (δ gap increased and δ num decreased). It may be augmented in the opposite case.

Numerically safe dual bound should be used ( K should be set to a large value) for problems for which the objective function value is small and integer. This is often true for problems similar to bin packing and vertex coloring, where one needs to minimize the number of columns in the solution. Strong branching parameters should be increased if good branching candidates are rare or the quality of branching candidates is very heterogeneous.

Restricted master and diving heuristics should be used if a good initial solution is not available. The diving heuristic may be inefficient if the pricing problem generates many paths which are not E-elementary.

Computational Experiments

The generic BCP solver optimization algorithms were coded in C++. The interface to the solver is implemented in Julia 1.2 language using JuMP [START_REF] Dunning | JuMP: A modeling language for mathematical optimization[END_REF] and LightGraphs packages. We also used:

• BaPCod C++ library [START_REF] Vanderbeck | BaPCod -a generic Branch-And-Price Code[END_REF] which implements the BCP framework;

• the C++ code developed by Sadykov et al. [START_REF] Sadykov | A bucket graph based labeling algorithm with application to vehicle routing[END_REF] implementing the bucket graph based labeling algorithm, bucket arc elimination procedure, path enumeration, and the separation of limited-memory rank-1 cuts;

• CVRPSEP C++ library [START_REF] Lysgaard | CVRPSEP: A package of separation routines for the capacitated vehicle routing problem[END_REF] which implements heuristic separation of rounded capacity cuts;

• IBM CPLEX Optimizer version 12.8.0 as the LP solver in column generation and as the solver for the enumerated MIPs.

• Boost C++ library version 1.55 [START_REF]Boost C++ libraries[END_REF] • LEMON C++ library version 1.31 [START_REF] Lemon | Library for Efficient Modeling and Optimization in Networks[END_REF] for the bucket graph creation for the labeling algorithm.

The experiments were run on a 2 Deca-core Ivy-Bridge Haswell Intel Xeon E5-2680 v3 server running at 2.50 GHz. The 128 GB of available RAM was shared between 8 copies of the algorithm running in parallel on the server. Each instance is solved by one copy of the algorithm using a single thread.

In Table 2, we show computational results for 14 applications. The first column is the problem acronym, second column refers to data sets, the third indicates the number of instances. Next is the time limit per instance. The last three columns show the results obtained by our generic solver, as well as by two other algorithms, those with the best (to our knowledge) published results for the data set. For each algorithm, we give the number of instances solved within the time limit, the average time in brackets (geometric mean time if the time limit is 10 hours or more), and its reference. For instances not solved, the time limit is considered as the solution time. Best results are marked in bold. The performance of our solver depends significantly on initial primal bounds given by the user. In the experiments, we always used the same primal bounds as in the works we compare with.

For each problem below we give details concerning the models and the parameterization employed, instances considered, initial primal bounds used, as well as analysis of computational results.

Capacitated Vehicle Routing Problem (CVRP)

The model is given in Section 4. The considered E-M instances are the 12 hardest ones, those considered in [START_REF] Pecin | Improved branch-cut-andprice for capacitated vehicle routing[END_REF]. The considered X instances are those with less than 400 customers. The number of vehicles for E-M instances is fixed as it is usual in the literature. The number of vehicles is unbounded for X instances. The same initial upper bounds are used as in the literature [START_REF] Pecin | Improved branch-cut-andprice for capacitated vehicle routing[END_REF][START_REF] Uchoa | New benchmark instances for the capacitated vehicle routing problem[END_REF].

The results show that our solver outperforms noticeably state-of-the-art on X instances. Results for E-M instances are comparable to the state-of-the-art. Note that the initial upper bound for the hardest instance M-n200-k16 is different by 4 units from the optimum solution. This difference introduce randomness to the running time of algorithms.

According to CVRPLIB (http://vrp.atd-lab.inf.puc-rio.br), in November 2018 there were 55 open CVRP instances in the X set [START_REF] Uchoa | New benchmark instances for the capacitated vehicle routing problem[END_REF]. We started long runs of the generic solver on the most promising ones, using a specially calibrated parameterization. We could solve 6 instances to optimality for the first time, as indicated in Table 3. Improved best known solutions are underlined.

Vehicle Routing Problem with Time Windows (VRPTW)

We use the same model as for the CVRP except that an additional main resource is defined which represents the time. Different bounds on the consumption of the time resource prevent us sec., φ bidir = 1, η max = 16, θ mem = 1, ω routes = 5

• 10 6 , ζ num 1 = 50, ζ estim 1 = ∞, ζ estim 2 = ∞, δ gap = 1.5%.
The results in Table 2 show that our solver outperformed significantly the algorithm proposed in [START_REF] Pecin | New enhancements for the exact solution of the vehicle routing problem with time windows[END_REF]. It is more than three times faster on average and could solve 7 more instances within the time limit. All Solomon instances with 100 customers are solved in less than 5 minutes each, except the hardest instance R208 which is solved in 37 minutes. The latter was solved in about 17 hours by the algorithm in [START_REF] Pecin | New enhancements for the exact solution of the vehicle routing problem with time windows[END_REF].

Heterogeneous Fleet Vehicle Routing Problem (HFVRP)

The model is given in Section 4. We use the instances with 50-100 clients proposed in [START_REF] Baldacci | A unified exact method for solving different classes of vehicle routing problems[END_REF]. They include instances with limited and unlimited fleet, with and without fixed vehicle cost, and with vehicle dependent and independent routing costs. The same initial upper bounds are used as in [START_REF] Pessoa | Enhanced branch-cut-and-price algorithm for heterogeneous fleet vehicle routing problems[END_REF]. Our solver significantly outperforms the algorithm from [START_REF] Baldacci | A unified exact method for solving different classes of vehicle routing problems[END_REF]. It also has a better performance than a recent branch-cut-and-price algorithm specific for that problem [START_REF] Pessoa | Enhanced branch-cut-and-price algorithm for heterogeneous fleet vehicle routing problems[END_REF].

Multi-Depot Vehicle Routing Problem (MDVRP)

We use a model similar to the one used on HFVRP, defining a graph G k for every depot k ∈ K. Variables associated with edges incident to a depot k are mapped to the corresponding arcs in graph G k . Variables associated with edges between customers are mapped to corresponding arcs in all graphs. The same parameterization and the same priorities for branching as for the HFVRP are used.

We use instances proposed in [START_REF] Cordeau | A tabu search heuristic for periodic and multi-depot vehicle routing problems[END_REF]. Only instances without distance constraints are considered, so there is a single capacity resource. The same initial upper bounds are used as in [START_REF] Pessoa | Enhanced branch-cut-and-price algorithm for heterogeneous fleet vehicle routing problems[END_REF]. The solver obtained a performance similar to [START_REF] Pessoa | Enhanced branch-cut-and-price algorithm for heterogeneous fleet vehicle routing problems[END_REF] and significantly better than [START_REF] Contardo | A new exact algorithm for the multi-depot vehicle routing problem under capacity and route length constraints[END_REF].

Pickup and Delivery Problem with Time Windows (PDPTW)

The model is given in Section 4.6. The parameterization of the solver is the following: τ soft = 5 sec., τ hard = 10 sec., ψ buck = 200, ω labels = 5 • 10 5 , ω routes = 2.5 • 10 6 , ω MIP = 7000, δ num = 2,

ζ estim 1 = 30, ζ estim 1 = 1.0, ζ num 2 = 2. Diving heuristic is used with parameters χ div = 10, χ depth = ∞, χ disc = 1.
The convention in PDPTW literature is to first minimize the number of routes, then minimize the transportation cost. So, we add the constant 10,000 to the cost of arcs leaving the depot. The same initial upper bounds are used as in [START_REF] Baldacci | An exact algorithm for the pickup and delivery problem with time windows[END_REF]. Results for this problem are mixed. Worse performance for Li& Lim instances can be explained by the fact that our solver does not incorporate some labeling algorithm acceleration techniques specific to the PDPTW. For the Ropke & Cordeau instances however, generic state-of-the-art BCP components mitigate the effect of lacking ad-hoc enhancements. Employing diving heuristic is important here as initial upper bounds are not tight for some instances.

Team Orienteering Problem (TOP)

The model is given in Section 4.5. The parameterization of the solver is the following: τ soft = 5 sec., τ hard = 10 sec., ψ buck = 200, φ bidir = 1, ω labels = 5 • 10 5 , ω routes = 2.5 • 10 6 , θ mem = 1,

δ num = 2, ζ estim 1 = 50, ζ estim 1 = 1.0.
Both restricted master and diving heuristics are used with parameters χ rm = 1 sec., χ div = 0, χ depth = ∞, χ disc = 1.

No initial upper bound is defined, so using heuristics is important. Instances of class 4, the most difficult one according to [START_REF] Bianchessi | A branch-and-cut algorithm for the team orienteering problem[END_REF], are considered for TOP. Our solver clearly outperforms the state-of-the-art on this problem. Its advantage can be explained by the fact that some important recent improvements in BCP algorithms for variants like CVRP and VRPTW were not yet adapted and used in TOP, possibly due to the complexity of their implementation.

Capacitated Team Orienteering Problem (CTOP)

The same model as for the TOP is used, except that an additional capacity resource is considered. The first main resource for the bi-directional labelling is the time resource. Pricing problem can be difficult for the instances considered. Therefore, we used a special parameterization for the number of buckets to be able to solve them in a reasonable time. Bucket arc elimination procedure and thus route enumeration are not used due to performance issues. Rank-1 cuts are also not separated due to the difficulty of the pricing problem.

The parameterization of the solver is the following: τ soft = 5 sec., τ hard = 10 sec., ψ buck = 1000, φ bidir = 1, φ elim = 0 (off),

σ stab = 1, η max = 30, θ rows = -1 (off), δ gap = 1.5%, ζ estim 1 = 50, ζ estim 1 = 1.0. Diving heuristic is used with parameters χ div = 0, χ depth = ∞, χ disc = 1.
No initial upper bound is defined. We used all instances in Set 1 from [START_REF] Archetti | The capacitated team orienteering and profitable tour problems[END_REF]. We also used only open instances in Set 2 from [START_REF] Archetti | The capacitated team orienteering and profitable tour problems[END_REF] (which were not solved by Archetti et al. [START_REF] Archetti | Optimal solutions for routing problems with profits[END_REF]), as Set 2 contains easier instances (with a reduced vehicle capacity and a reduced time limit). We have not solved only one instance: "p09". All other instances were solved to optimality in less than 15 minutes each, seven of them for the first time.

Capacitated Profitable Tour Problem (CPTP)

The same model as for the CTOP, except that i) the only resource is the vehicle capacity, ii) the objective is the difference between the total profit and the transportation cost. The parameterization of the solver is the following: τ soft = 5 sec., τ hard = 10 sec., ψ buck = 200, φ bidir = 1, σ stab = 1, ω labels = 5 • 10 5 , ω routes = 2.5 • 10 6 , η max = 30, θ mem = 1, δ gap = 1.5%,

ζ estim 1 = 50, ζ estim 1 = 1.0.
Only open instances from [START_REF] Archetti | The capacitated team orienteering and profitable tour problems[END_REF] are considered, which could not be solved by Bulhoes et al. [START_REF] Bulhoes | The vehicle routing problem with service level constraints[END_REF]. The same initial upper bounds were used as for the branch-and-price algorithm in [START_REF] Bulhoes | The vehicle routing problem with service level constraints[END_REF]. Note that only one of them was improved: for instance "p13-4-200" the optimum value is 304.15, whereas the best known solution is 303.18. So the heuristic suggested in [START_REF] Bulhoes | The vehicle routing problem with service level constraints[END_REF] is of a very good quality. Among 28 open instances of the CPTP, we solved to optimality 24 within 1 hour. Two more instances "p13-4-200" and "p10-20-200" were solved within 2 hours each. The only remaining open instances are now "p15-15-200", and "p16-20-200".

VRP with Service Level constraints (VRPSL)

VRPSL is a generalisation of CVRP in which a service weight is defined for each customer. For each predefined group of customers, total service weight of visited customers should not be below a threshold. The model contains edge variables x and visiting variables y, as for the TOP. For each group, a knapsack constraint over y variables is defined. Branching is both on y and x variables. The parameterization of the solver is the following: τ soft = 3 sec., τ hard = 6 sec., ψ buck = 200, ω labels = 2 • 10 5 , η max = 30, θ rows = 4,

θ mem = 1, δ gap = 1.5%, ζ num 1 = 50, ζ estim 1 = 1.0.
The instances proposed in [START_REF] Bulhoes | The vehicle routing problem with service level constraints[END_REF] are considered. The same initial upper bounds were used as in [START_REF] Bulhoes | The vehicle routing problem with service level constraints[END_REF]. None of them was improved. Our solver outperformed largely the branch-and-price algorithm by Bulhoes et al. [START_REF] Bulhoes | The vehicle routing problem with service level constraints[END_REF]. The reason is that the latter does not use many state-of-theart techniques for routing problems. Performance of our solver can probably be improved by separating valid inequalities for the knapsack constraints.

Generalized Assignment Problem (GAP)

The model is given in Section 4.1. We used 6 classic OR-Library instances of the most difficult type D with up to 20 machines and 200 tasks. Also we used instances by Nauss [START_REF] Nauss | Solving the generalized assignment problem: An optimizing and heuristic approach[END_REF] with |T | = 90, 100 and |K| = 25, 30. The parameterisation of the solver is the following: ψ buck = 100, ψ reduc = 0 (off), φ bidir = 1, γ exact = 10, γ heur = 10, ω labels = 10 5 , ω routes = 2 • 10 5 , ω MIP = 4000, η init = 0, η max = 0, θ rows = 4. Diving heuristic is used with parameters χ div = 0, χ depth = 3, χ disc = 2, as no initial upper bound is defined for both types of instances.

Performance of our solver is slightly worse than the state-of-the-art. Both "competitors" [START_REF] Posta | An exact method with variable fixing for solving the generalized assignment problem[END_REF] and [START_REF] Avella | A computational study of exact knapsack separation for the generalized assignment problem[END_REF] solved one more instance to optimality. In Table 2, times of [START_REF] Posta | An exact method with variable fixing for solving the generalized assignment problem[END_REF] and [START_REF] Avella | A computational study of exact knapsack separation for the generalized assignment problem[END_REF] are adjusted according to computer speeds. Note that our solver is the first algorithm in the literature which solves the pricing problem as the resource constrained shortest path problem. Although it is slower than specialised knapsack solvers, it supports rank-1 cuts and enumeration.

Our solver is much more efficient than the algorithm by Nauss [START_REF] Nauss | Solving the generalized assignment problem: An optimizing and heuristic approach[END_REF] and the MIP formulation for the GAP solved by Gurobi. Note that our solver obtains particularly good results for instances with relatively small number of tasks per machine, as in this case path enumeration procedure is very efficient. For instances with large number of tasks per machine our solver is less efficient. In particular, a more advanced stabilisation technique is required for such instances, as discussed in [START_REF] Pessoa | Automation and combination of linear-programming based stabilization techniques in column generation[END_REF].

Bin Packing Problem (BPP)

The model is given in Section 4.2. The parameterization of the solver is the following: ψ buck = 200, ψ reduc = 0 (off), φ bidir = 1, γ exact = 100, γ heur = 100, σ stab = 2, ω routes = 2 • 10 6 , ω MIP = 10 5 , η init = 0, η max = 0, θ rows = 4, δ num = 5, K = 10 10 , ζ num 1 = 20, ζ num 2 = 1. Diving heuristic with parameters χ div = 10, χ depth = 0, χ disc = 0 is used for all instances except set ANI.

We give smaller priority for the Ryan and Foster branching and larger priority for the branching over the accumulated resource consumption. The latter showed to be so effective that Ryan and Foster rule was never used in our experiments. For each instance, we use initial primal bound which is equal to the rounded up value of the column generation dual bound plus 1 unit (there is a long-standing conjecture that the optimal solution of a BPP instance is never larger than this). Solutions with these objective values are easily obtainable by simple heuristics.

Our solver obtains the best results for the two most difficult instance classes AI and ANI. For other less difficult instances, our algorithm is slower than the state-of-the-art. However, it can solve all instances to optimality in a relatively small time. Note that the algorithm of Clautiaux et al. [START_REF] Clautiaux | Iterative aggregation and disaggregation algorithm for pseudo-polynomial network flow models with side constraints[END_REF] showed better results for the class AI. However, the authors communicated to us that they discovered an issue with their code. Therefore, their results are not included in Table 2.

The bottleneck of our algorithm for solving instances of classes AI and ANI with 600 items or more is the LP solver. For such instances, very often when the pricing problem finds columns with negative reduced cost, LP solver does not include them in the basis. This happens because the absolute value of the reduced cost is smaller than the minimum possible reduced cost tolerance of the LP solver. Thus, the final dual solution of column generation is not feasible, and the Lagrangian bound obtained by the safe procedure is weaker than it can potentially be.

Vector Packing Problem (VPP)

The model is given in Section 4.2. The parameterization of the solver is the following: ψ buck = 2000, ψ reduc = 0 (off), φ bidir = 1, γ exact = 100, γ heur = 100, σ stab = 1, ω labels = 10 5 , ω MIP = 10 5 , η init = 0, η max = 0, θ rows = -1 (off), ζ num 2 = 1. We use diving heuristic with parameters χ div = 0, χ depth = 2, χ disc = 3.

We do not use initial upper bounds. We consider only largest instances from the literature with 200 items and only with 2 resources. We use instances of the most difficult, according to [START_REF] Heßler | Stabilized branch-and-price algorithms for vector packing problems[END_REF], classes 1, 4, 5, and 9. Other instances are significantly easier, all of them are solved in the literature. Our solver outperforms largely the algorithms in the literature.

Capacitated Arc Routing Problem (CARP)

The model is defined in Section 4.7. The parameterization of the solver is the following: τ soft = 6 sec., τ hard = 15 sec., ψ buck = 50, γ exact = 50, γ heur = 300, σ stab = 1, ω labels = 5 • 10 5 , θ num = 50, δ gap = 1%, δ num = 5, ζ estim 1 = 1.0, ζ num 2 = 5. In the distance matrix for ng-sets, the distance between two required edges is defined as the sum of costs of four paths between midpoints of the edges. Each path is the shortest path starting from a given vertex incident to the first edge and ending at a given vertex incident to the second edge.

The branching is done on the aggregation of x variables corresponding to node degrees in the original graph or on the aggregation of x variables corresponding to whether required two edges are served immediately one after another by the same route or not. The same priority is used for both branching strategies. The Eglese dataset [START_REF] Eglese | Efficient routeing for winter gritting[END_REF] is standard in the literature and it is used in all recent works on the CARP. We have used the same initial upper bounds as in [START_REF] Pecin | Comparative analysis of capacitated arc routing formulations for designing a new branch-cut-and-price algorithm[END_REF]. The performance of our solver is similar to the most recent exact algorithm [START_REF] Pecin | Comparative analysis of capacitated arc routing formulations for designing a new branch-cut-and-price algorithm[END_REF] for the problem. Other algorithms in the literature are significantly less efficient. The generality of our solver opens a way to quickly obtain excellent computational results for many variants of the arc routing problems.

Black-and-White Traveling Salesman Problem (BWTSP)

The model is defined in Section 4.8. The parameterization of the solver is the following: τ hard = 25 sec., ω labels = 5 • 10 5 ω routes = 5 • 10 6 , ω MIP = 5000, θ mem = 1, ζ estim One third of 720 instances do not have the constraint on the maximum distance between visiting two black vertices. For these instances, only one (capacity) resource is defined. For other instances, the first main resource is the distance resource.

Our solver is much more efficient for the GLR instances than the best approach in the literature [START_REF] Gouveia | Extended formulations and branchand-cut algorithms for the black-and-white traveling salesman problem[END_REF], which is a branch-and-cut algorithm. The average solution time is 72 times smaller for our solver.

Results with the default parameterisation

In Table 4 we present results for 7 VRP variants obtained using the default parameterisation of the solver. The built-in primal heuristics are used, with default parameters, if and only if no initial primal bound is specified. It can be seen that the performance is deteriorated, but not much. This indicates that finding a good problem-specific parameterization is important but not critical.

Results without defining packing sets

The concept of packing set (and the related concept of elementarity set) is crucial to our models, since it is used in all the advanced BCP components described in Section 3. In some models, like those for CVRP and VRPTW, the packing sets are very simple: each customer vertex defines a packing set. However, there are models where the packing sets are more complex. For examples:

(1) on GAP each packing set is composed by |K| arcs, but only one of the two arcs entering each of the |K| vertices representing a task belongs to the packing set; (2) on CARP each packing set is composed by a pair of vertices. In order to justify the need for such a general concept of packing sets, we test the impact of not using packing sets on GAP and CARP.

In those tests, the solver parameterisation was the same as for the main experimental results. Of course, when packing sets are not defined, the solver cannot use ng-path relaxation, rank-1 and capacity cuts, enumeration, and specific branching. For the GAP, if packing sets are not defined, only two OR-Lib instances are solved in 2 hours (one in 16 minutes, another in 53 minutes) and none of Nauss instances has been solved in 1 hour. For the CARP, only one Eglese instance has been solved within the time limit of 30 hours (it was solved in 36 minutes). Comparison with results in Table 2 makes clear that defining packing sets is absolutely necessary for having a reasonable performance on those problems.

A deeper conclusion of the experiment is the following. If we had defined the advanced BCP components using a simpler concept, like "the set of customers", we could have built a general solver with a good performance on a number of classic VRP variants. However, that solver would not work well on other important problems. In order to have a good performance over a wider range of problems, the more general concept of packing set, including its arc and vertex variants, is indeed needed. that their use may become standard, at least for the purpose of having baseline results to be compared with results of proposed specialized algorithms. Furthermore, we believe that there is plenty of room for "creative modeling", where users may find original ways of fitting new problems into the proposed model. In fact, as already demonstrated on generalized assignment problem and on bin/vector packing problems, not only VRP variants can be efficiently treated. It may be also possible to model problems from scheduling, from network design and from other discrete optimization subareas. Several emerging VRP variants include complex features like time-dependent travel times, stochastic components or synchronization. There are two recent works where creative modeling was used for handling some of those complex variants with the proposed VRP solver: (1) on twoechelon VRP, [START_REF] Marques | An improved branch-cut-and-price algorithm for the two-echelon capacitated vehicle routing problem[END_REF] proposes an exponentially-sized number of Master constraints for modeling the synchronization between the two echelons; (2) on robust CVRP under knapsack uncertainty, [START_REF] Pessoa | Branch-and-cutand-price for the robust capacitated vehicle routing problem with knapsack uncertainty[END_REF] proposes a transformation for representing uncertain demands by a set of deterministic RCSP problems. Yet, most of those emerging variants are out of the modeling capabilities of the VRP solver.

We believe that the most promising course for future improvements in the modeling power of the solver is to add the possibility of using other types of resources. This may include resources with arc consumption and/or cost dependent on its own accumulated consumption or even dependent on accumulated consumption of other resources, resources with soft or multiple interval limits, non-linear and stochastic resources, and others, as discussed in [START_REF] Irnich | Resource extension functions: Properties, inversion, and generalization to segments[END_REF] and in [START_REF] Parmentier | Algorithms for non-linear and stochastic resource constrained shortest path[END_REF]. However, devising and implementing algorithms that support any of those more complex resources, in an efficient way and preserving the compatibility with all the existing features of the solver, will be a major challenge.

  Consider the instance over a complete graph G with n = 4, d 1 = 3, d 2 = 4, d 3 = 5, d 4 = 7, and Q = 10. Edge costs are c 01 = 5, c 02 = 9, c 03 = 7, c 04 = 3, c 12 = 7, c 13 = 8, c 14 = 7, c 23 = 4, c 24 = 9, and c 34 = 5. The resulting directed graph G is depicted in Figure 1, arc consumptions and mappings are indicated. Graph G generates a set P having 13 resource constrained paths: 0

Figure 1 :

 1 Figure 1: CVRP model graph and mapping. The resource consumptions (in bold) and mappings (the same for each pair of opposite arcs) are depicted. The upper bound on accumulated resource consumption is 10 for all vertices.

4. 2

 2 Vector Packing (VPP) / Bin Packing (BPP) Data: Set T of items; set D of dimensions; bin capacities Q d , d ∈ D; item weight w d t , t ∈ T , d ∈ D. (Bin packing is the case where |D| = 1).

  v

Figure 2 :

 2 Figure 2: GAP model graph, RCSPs correspond to binary knapsack solutions.

4. 6

 6 Pickup and Delivery VRPTW (PDPTW) Data: Directed graph G = (V, A), where V = {0} ∪ P ∪ D, P = {1, . . . , n} is the set of pickup vertices and D = {n + 1, . . . , 2n} the set of corresponding deliveries (a pickup at i correspond to a delivery at i + n); vehicle capacities Q; traveling cost c a and time (including service time) t a , a ∈ A; positive demands d v , v ∈ P (d v = -d v-n , v ∈ D); and time windows [l v , u v ], v ∈ V .

1 Figure 3 :

 13 Figure 3: BWTSP model graph and mapping. Consumption of the capacity resource is shown in bold. Dashed arcs have zero consumption for both resources.

3 . 1 =

 31 The parameterization of the solver is the following (values different from defaults): τ hard = 25 sec., ω routes = 5 • 10 6 , θ mem = 0, ζ estim 02, δ gap = 1.5%. The distance matrix between elementarity sets corresponds to the distance matrix between clients. The same matrix is used in all routing problems below, except CARP.

4 . 1 = 0 . 2 , ζ num 2 = 5 , ζ estim 2 =

 4102252 The parameterization of the solver is the following: τ soft = 6 sec., τ hard = 15 sec., ω routes = 5 • 10 6 , θ mem = 1, ζ estim 0.02. The branching on the number of paths for each vehicle type has the highest priority. Other branching strategies have the same priority.

1 =

 1 02, δ gap = 1.5%. We use restricted master heuristic with parameter χ rm = 10 sec.

  0, and 0-4-1-0. The complete formulation (corresponding to Formulation (1)) for that instance is: Min 5x 01 + 9x 02 + 7x 03 + 3x 04 + 7x 12 + 8x 13 + 7x 14 + 4x 23 + 9x 24 + 5x 34

S.t.

x 01 + x 12 + x 13 + x 14 = 2,

x 02 + x 12 + x 23 + x 24 = 2,

x 03 + x 13 + x 23 + x 34 = 2,

x 04 + x 14 + x 24 + x 34 = 2,

  Min 10λ 1 + 21λ 2 + 24λ 3 + 20λ 4 + 15λ 5 + 18λ 6 + 21λ 7 + 20λ 8 + 14λ 9 + 20λ 10 + 20λ 11 + 6λ 12 + 15λ 13 S.t. 2λ 1 + 2λ 2 + 4λ 3 + 2λ 4 + 2λ 5 + 2λ 7 + 2λ 10 + 2λ 13 = 2 2λ 2 + 2λ 3 + 2λ 6 + 2λ 7 + 2λ 8 + 2λ 11 = 2 2λ 4 + 2λ 8 + 2λ 9 + 2λ 10 + 2λ 11 = 2

10

+ , λ ∈ Z 13 + Eliminating the x variables and relaxing the integrality, the Master LP (corresponding to (2)) that should be solved by column generation is obtained:

  Branching on x or on y variables. Enumeration is on.Comments: The y variables, that indicate which customers are visited, are not mapped to any arc. So, they are kept in the Master LP.

Table 4 :

 4 Problem-dependent parameterisation vs default parameterisation

	Problem Data set	# T.L. Specific params Default params
	CVRP E-M [25, 26]	12 10h	12 (61s)	11 (65s)
		X [82]	58 60h	36 (147m)	35 (142m)
	VRPTW Solomon Hardest [81] 14 1h	14 (5m)	13 (10m)
		Homberger 200 [40] 60 30h	56 (21m)	52 (31m)
	HFVRP BaldacciMingozzi [9] 40 1h	40 (144s)	40 (191s)
	MDVRP Cordeau [29]	11 1h	11 (6m)	11 (7m)
	PDPTW RopkeCordeau [76] 40 1h	40 (5m)	38 (8m)
	TOP	Chao class 4 [24]	60 1h	55 (8m)	54 (9m)
	CARP Eglese [36]	24 30h	22 (36m)	21 (41m)

Conclusions

We proposed a new generic way of modeling VRPs and related problems, so that they can be solved by a BCP algorithm that already includes most state-of-the-art components introduced for the most classical VRP variants. It combines existing modeling concepts, like the use of RCSPs for defining the valid routes, with new ones: (1) the concept of mapping for linking the route variables with a generic set of constraints in the Master formulation; (2) the concepts of packing and elementarity sets for allowing the use of some advanced BCP components over a wide range of problems. The experiments show that the generic solver has a performance either comparable or better than the specific algorithms for all VRP variants tested. The cases where the performance was much better can be explained by the fact that previous works on some problems did not use some of those advanced components, possibly due to the complexity of their implementation. However, if generic VRP solvers become publicly available, we believe

Acknowledgements

We would like to thank Teobaldo Bulhoes and Guillaume Marques for a large part of the implementation of the Julia-JuMP interface to the solver; Teobaldo Bulhoes, Guillaume Marques and Eduardo Queiroga for implementing, over that interface, the models corresponding to the examples of this paper; and Laurent Facq for a general support of the computing environment. Experiments presented in this paper were carried out using the PlaFRIM (Federative Platform for Research in Computer Science and Mathematics), created under the Inria PlaFRIM development action with support from Bordeaux INP, LABRI and IMB and other entities: Conseil Régional d'Aquitaine, Université de Bordeaux, CNRS and ANR in accordance to the "Programme d'Investissements d'Avenir". This study was financed in part by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), grant 313601/2018-6 (Produtividade 1B), and by the Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ), grant E-26/202.887/2017 (Cientista do Estado).

Problem Data set # T.L. Gen. BCP Best Publ. 2nd Best CVRP E-M [START_REF] Christofides | An algorithm for the vehicle-dispatching problem[END_REF][START_REF] Christofides | Combinatorial Optimization, chapter The vehicle routing problem[END_REF] 12 10h 12 (61s) 12 (49s) [START_REF] Pecin | Improved branch-cut-andprice for capacitated vehicle routing[END_REF] 10 (432s) [START_REF] Contardo | A new exact algorithm for the multi-depot vehicle routing problem under capacity and route length constraints[END_REF] X [START_REF] Uchoa | New benchmark instances for the capacitated vehicle routing problem[END_REF] 58 60h 36 (147m) 34 (209m) [START_REF] Uchoa | New benchmark instances for the capacitated vehicle routing problem[END_REF] -VRPTW Solomon Hardest [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF] 14 1h 14 (5m) 13 (17m) [START_REF] Pecin | New enhancements for the exact solution of the vehicle routing problem with time windows[END_REF] 9 (39m) [START_REF] Baldacci | New route relaxation and pricing strategies for the vehicle routing problem[END_REF] Homberger 200 [START_REF] Gehring | Parallelization of a two-phase metaheuristic for routing problems with time windows[END_REF] 60 30h 56 (21m) 50 (70m) [START_REF] Pecin | New enhancements for the exact solution of the vehicle routing problem with time windows[END_REF] 7 (-) [START_REF] Kallehauge | Lagrangian duality applied to the vehicle routing problem with time windows[END_REF] HFVRP BaldacciMingozzi [START_REF] Baldacci | A unified exact method for solving different classes of vehicle routing problems[END_REF] from exploiting forward-backward path symmetry in the pricing.

The considered Solomon instances (all with 100 customers) are the hardest ones according to [START_REF] Pecin | New enhancements for the exact solution of the vehicle routing problem with time windows[END_REF]. As it is standard in the literature, we divide the instances in two groups. First group contains instances of type C1, RC1, and R1. For this group of instances, the first main resource for the bi-directional labelling is the capacity resource. Second group contains instances of type C2, RC2, R2. For the second group of instances, capacity resource is not used, and the feasibility related to the vehicle capacity is guaranteed by the rounded capacity cuts, separated both for fractional and integer solutions of the master problem. Therefore, those cuts become essential constraints and the sufficient condition for enumeration given in Section 3.4 is not verified (paths visiting the same customers in distinct orders may have different coefficients in some RCC). However, it is still possible to show that the enumeration procedure is valid. The same initial upper bounds are used as for the algorithm in [START_REF] Pecin | New enhancements for the exact solution of the vehicle routing problem with time windows[END_REF].

The parameterization of the solver for the first group of instances is τ soft = 5 sec., τ hard = 10 sec., φ bidir = 1, η max = 16, θ rows = 7, θ mem = 1, ω routes = 10