Ruslan Sadykov

Eduardo Uchoa

Artur Pessoa

A Bucket Graph Based Labeling Algorithm with Application to Vehicle Routing

We consider the Shortest Path Problem with Resource Constraints (SPPRC) arising as a subproblem in state-of-the-art Branch-Cut-and-Price algorithms for vehicle routing problems. We propose a variant of the bi-directional label correcting algorithm in which the labels are stored and extended according to the so-called bucket graph. Such organization of labels helps to decrease significantly the number of dominance checks and the running time of the algorithm. We also show how the forward/backward route symmetry can be exploited and how to eliminate arcs from the bucket graph using reduced costs. The proposed algorithm can be especially beneficial for vehicle routing instances with large vehicle capacity and/or with time window constraints. Computational experiments were performed on instances from the distance constrained vehicle routing problem, including multi-depot and site-dependent variants, on the vehicle routing problem with time windows, and on the "nightmare" instances of the heterogeneous fleet vehicle routing problem. Significant improvements over the best algorithms in the literature were achieved and many instances could be solved for the first time.

has several degrees of freedom. In particular, there are many possible ways of choosing the next label to be extended and how often and how extensively dominance checks are performed. Particular labeling algorithms are classified as being either label-setting or label-correcting. The defining property of a label-setting algorithm is that it only extends labels that can never be dominated by another label created after that extension.

 and Pugliese and Guerriero (2013) for surveys. However, somehow surprisingly, not so many papers discuss in depth issues related to label organization and dominance strategies. In fact, even though the pricing time is a bottleneck in many BCP algorithms for VRP, authors often skip the details of the labeling algorithm used. An exception is the recent algorithm for the Capacitated VRP (CVRP) in Pecin et al. (2017b). In that case, assuming that the consumptions of the single resource (capacity) are given by positive integer numbers, labels are organized in buckets corresponding to each possible consumption. Dominance is only performed for labels in the same bucket, that by definition have the same resource consumption. However, that approach cannot be applied in cases where resource consumptions are given by arbitrary real numbers. In fact, if the consumptions are given by larger integer numbers the algorithm becomes slow, since few dominance checks are performed and too many undetected dominated labels are expanded. The opposite approach is to perform full label dominance after each extension. In this case, the running time may suffer a lot from too many dominance checks with negative results.

. It is not known how Pulse would perform on larger instances and how it would handle the modifications in the SPPRC induced by cuts, essential for modern BCP algorithms for VRP.

The main original contributions of this paper are the following:

• In Section 3, a new variant of the labeling algorithm is proposed. The approach relies on a so-called bucket graph, consisting of buckets and bucket arcs. Labels for paths ending in the same vertex of the original graph and having similar resource consumption are grouped together in buckets. A bucket arc links two buckets if a label in the first one can be possibly extended to a label in the second. Labels within the same bucket are always pairwise undominated, i.e. there is no label which is dominated by another label in the same bucket. However, dominance checks between labels in different buckets are only performed before an extension, when it is checked whether a candidate label to be extended is dominated or not. Moreover, this inter-bucket dominance uses bounds on the costs of all labels in a bucket, avoiding many unnecessary checks. The key parameter of the algorithm is the step size, used for determining the resource consumption intervals that define each bucket. If the step size is sufficiently small, the bucket graph is acyclic and the algorithm becomes a label-setting algorithm. In general, the bucket graph has cycles and

Introduction

The best known exact approach for many classical variants of the Vehicle Routing Problem (VRP) is a Branch-Cut-and-Price (BCP) algorithm in which the master problem contains route variables, customer visiting constraints and additional cuts. Because of their exponential number, the route variables are generated dynamically by solving pricing subproblems modeled as Shortest Path Problems with Resource Constraints (SPPRC). In those problems, one looks for least cost paths joining a source vertex with a sink vertex such that the accumulated consumption of resources along the path respects given lower and upper limits. Dynamic Programming labeling algorithms are the most usual way of solving those problems.

Labeling algorithms differ from the traditional "table filling" dynamic programming algorithms because they only store (as labels) the reachable states, those representing feasible partial paths. More importantly, while traditional dynamic programming only performs dominance over identical states, labeling algorithms perform dominance between labels corresponding to non-identical states. Typically, a label L representing a path P L is dominated, and therefore eliminated, if there is another label L representing a path P L that ends in the same vertex and is not more costly and does not use more resources than P L . Mono-directional labeling algorithms start from a single label representing a null path at the source vertex. At each step, a non-extended label L is extended to additional labels corresponding to the possible ways of adding a single arc to P L . Dominance may be used to eliminate labels, avoiding future extensions. If there are no non-extended labels, the algorithm stops and optimal solutions are found among the labels representing paths ending in the sink vertex. The general labeling algorithm the algorithm is a label-correcting algorithm. A bi-directional version of the bucket graph based labeling algorithm is also proposed (since [START_REF] Righini | Symmetry helps: Bounded bi-directional dynamic programming for the elementary shortest path problem with resource constraints[END_REF] it is known that bi-directional labeling algorithms are often more efficient than their mono-directional counterparts). The concatenation of labels step is also accelerated by bounds on the costs. Bi-directional search also allows us to exploit the forward/backward symmetry of routes that exists on some VRP variants, for example, on the classical Capacitated VRP (CVRP). On those cases, the backward search is replaced by a "reversed copy" of the forward search, reducing the running time.

• In Section 4, a procedure for accelerating the labeling algorithm by removing arcs from the bucket graph using reduced cost arguments is proposed. We call it bucket arc elimination procedure. The concept of jump bucket arcs is introduced. The new fixing procedure generalizes both the procedures in [START_REF] Irnich | Path-reduced costs for eliminating arcs in routing and scheduling[END_REF] and in [START_REF] Pessoa | Exact algorithm over an arc-time-indexed formulation for parallel machine scheduling problems[END_REF]. After a fixing, some bucket arcs are eliminated, and thus the number of extensions and the running times are reduced in future calls to the labeling algorithm.

The new labeling algorithm was embedded as the pricing method in a modern BCP algorithm for VRPs. Our BCP also includes many ingredients found in other state-of-the-art BCP algorithms, like dynamic ng-path relaxation, rounded capacity cuts, limited arc memory rank-1 cuts, automatic dual price smoothing stabilization, a procedure for enumerating elementary paths, and multi-phase strong branching with pseudo-costs. The algorithm was computationally tested on several VRP variants with time window constraints and/or large vehicle capacity. Pretests have shown that the labeling algorithm could perform very well, but that performance is rather sensitive to the choice of the step size parameter. Moreover, the good step size values varied strongly from instance to instance. Therefore, a simple but effective scheme for an automatic dynamic adjustment of the step size was devised. Other specific experiments indicate that, on the hardest instances, the new bucket arc fixing procedure can indeed be better than existing schemes for fixing arcs by reduced cost. Lastly, extensive experiments show that the final BCP algorithm outperforms significantly other recent state-of-art algorithms for the VRPTW and the Multi-Depot VRP with Distance Constraints (MDVRPDC). Moreover, our algorithm is the first exact approach that can handle medium-sized instances of the classical Distance Constrained Vehicle Routing Problem (DCVRP), including the site-dependent variant. Our algorithm is also the first exact method successfully applied for a set of "nightmare" instances of the Heterogeneous Fleet VRP (HFVRP).

The remainder of this article has the following structure. Section 2 defines the exact SPPRC solved in this paper. It is explained how this SPPRC arises as a pricing subproblem for solving a Heterogeneous Fleet Vehicle Routing Problem with Time Windows (HFVRPTW) by a columnand-cut generation algorithm. Section 3 describes the bucket graph labeling algorithm proposed in this work. Section 4 presents the bucket arc elimination procedure and introduces the key concept of jump bucket arcs. Section 5 presents the complete BCP algorithm for HFVRPTW where the bucket graph labeling algorithm was embedded. Section 6 presents results of the computational experiments. There are experiments devised to assess the impact of the bucket step size and experiments for measuring the impact of the introduced bucket arc elimination procedure. Moreover, extensive experiments evaluate the overall BCP performance on several classical VRP variants that are particular cases of the HFVRPTW. Section 7 contains final remarks. Finally, very detailed computational results are presented in the online appendix.

Application

In order to give a precise definition of the SPPRC addressed in this paper, we first need to define its application; defining the family of VRPs that we ultimately want to solve and outlining the column-and-cut generation algorithm where that SPPRC arises as the pricing subproblem.

HFVRPTW Definition

The Heterogeneous Fleet Vehicle Routing Problem with Time Windows (HFVRPTW) [START_REF] Jiang | Vehicle routing problem with a heterogeneous fleet and time windows[END_REF]) is defined as follows. Let G = (V, A) be a directed graph with vertex set V = {0, . . . , n + 1}. For convenience, the depot vertex is split into the source vertex 0 and the sink vertex n + 1; vertices in V = V \ {0, n + 1} represent the n customers. The arc set is defined as (v, v ) = (0, n + 1)}. Each vertex v ∈ V has a demand w v , a service time (duration) s v , and a time window [l v , u v ] associated with it. Demands and service times are non-negative for customers v ∈ V and equal to zero for the depot vertices, l 0 is assumed to be equal to 0. The fleet is composed of a set M of different types of vehicles. For each m ∈ M , there are U m available vehicles, each with a capacity W m . Let W = max{W m : m ∈ M } be the largest capacity. Every vehicle type is associated with a fixed cost denoted by f m . For each arc a ∈ A and m ∈ M there is a cost c m a and a non-negative travel time t m a associated to the traversal of this arc by a vehicle of type m. We assume that G does not have any cycle where, for some vehicle type, all vertices have zero demands and zero service times and all arcs have zero time. A route P = (v 0 = 0, v 1 , . . . , v k , v k+1 = n + 1) for a vehicle of type m ∈ M is a walk where v 1 , . . . , v k are not necessarily distinct customers in V and is said to be feasible if 1) it satisfies vehicle capacity: k j=1 w vj ≤ W m ; and 2) the earliest start of service time S j at every visited vertex v j , 0 ≤ j ≤ k + 1, falls within the corresponding time window: l vj ≤ S j ≤ u vj , where S 0 = l 0 and S j = max{l vj , S j-1 + s vj-1 + t m (vj-1,vj ) }. A vehicle may arrive at a vertex before the beginning of its time window and wait, but it cannot arrive after the end of the time window. The cost of route P is calculated as f m + k+1 j=1 c m (vj-1,vj ) . A route is said to be elementary if the same customer is not visited more than once. The objective is to determine a set of feasible elementary routes with minimal total cost such that: (i) each customer is visited by exactly one route; (ii) the number of routes associated to a vehicle type does not exceed its availability. Clearly, some of most classical VRP variants, like CVRP, VRPTW and HFVRP, are particular cases of HFVRPTW. Other variants that are also particular cases of HFVRPTW are listed below.

A = {(v, v ) : v, v ∈ V, v = v , v = n + 1, v = 0,
In the Multi-Depot VRP (MDVRP), the vehicles themselves are identical, they only differ by being attached to different depots. It can be easily modeled as a HFVRP by associating each depot to a vehicle type and setting costs c m (0,v) and times t m (0,v) for leaving the depot, together with costs c m (v,n+1) and times t m (v,n+1) for entering the depot, that depend on m ∈ M . The Site Dependent VRP (SDVRP), a variant where vehicles differ only by the subset of the customers that they can visit, is also easily modeled as a HFVRP by setting infinite costs for c m (v,v ) if vehicle type m cannot visit either v or v .

The classical Distance Constrained Vehicle Routing Problem (DCVRP) can also be solved as a HFVRPTW. In this variant, besides the capacity constraint, the route length is also forbidden to be above a certain threshold D. Since service times are included in that "length", this is actually a time limit for returning to the depot. It can be modeled by setting time window [0, D] for every v ∈ V .

Set Partitioning Formulation and Cuts

Let Ω m be the set of all feasible elementary routes for vehicle type m ∈ M . A set partitioning formulation over those sets would lead to a hard pricing problem, probably intractable on large instances. A more tractable (but weaker) formulation can be defined using a relaxed set of routes that includes some non-elementary routes, i.e. routes in which some customers are visited more than once. In this work we assume that the ng-route relaxation [START_REF] Baldacci | New route relaxation and pricing strategies for the vehicle routing problem[END_REF], a route relaxation with a good tradeoff between formulation strength and pricing difficulty, is being used. Let N v ⊆ V be the neighborhood of v ∈ V , typically containing the closest customers to v. An ng-route can only revisit a customer v, forming a cycle, if the cycle contains another customer v with v / ∈ N v . In many cases, reasonably small neighborhoods (for example, with |N v | = 8) already provide bounds that are close to those that would be obtained by pricing elementary routes (Poggi andUchoa 2014, Contardo, Desaulniers, and[START_REF] Contardo | Reaching the elementary lower bound in the vehicle routing problem with time windows[END_REF].

Let Ω N m ⊇ Ω m be the set of ng-routes for vehicle type m ∈ M with respect to a given set of neighborhoods N = (N 1 , . . . , N n ). We denote by c P the cost of route P ∈ Ω N m , by x (v,v ) the number of times arc (v, v ) ∈ A appears in route P , and by y P v = a∈δ + ({v})∪δ -({v}) x P a /2 the number of times vertex v ∈ V is visited in route P (this "symmetric" definition of y P v is exploited in Section 3.6). Let λ P be a binary variable indicating whether route P is selected or not in the solution. The Set Partitioning Formulation (SPF) for the HFVRPTW considered in this paper is the following.

(SPF) min

m∈M P ∈Ω N m c P λ P (1) subject to m∈M P ∈Ω N m y P v λ P = 1 ∀v ∈ V , (2)
P ∈Ω N m λ P ≤ U m ∀m ∈ M, (3)
λ P ∈ {0, 1} ∀m ∈ M, P ∈ Ω N m . (4)
A column generation algorithm should be used to solve the linear relaxation of ( 1)-( 4). However, even using large neighborhoods (or even enforcing that all routes are elementary) the resulting bounds are often not good enough for building an effective branch-and-price algorithm. Therefore, the SPF should be reinforced by adding cuts.

Let C ⊆ V be a subset of the customers, define w(C) = v∈C w v as its total demand. The value w(C)/W is a valid lower bound on the number of vehicles that must visit C. Therefore, the following Rounded Capacity Cut (RCC) [START_REF] Laporte | A branch and bound algorithm for the capacitated vehicle routing problem[END_REF] is valid:

m∈M P ∈Ω N m   a∈δ + (C)∪δ -(C) x P a   λ P ≥ 2 w(C)/W . (5)
RCCs are known to provide a significant reinforcement of the SPF on the CVRP [START_REF] Fukasawa | Robust branch-and-cut-and-price for the capacitated vehicle routing problem[END_REF]). However, they do not work so well on many HFVRP and VRPTW instances. This happens because the expression w(C)/W , which depends only on the largest vehicle capacity W and disregards time windows, can be a poor bound on the actual number of vehicles visiting C. The coefficient of a variable λ P in a RCC is given by a linear expression over the values of x P a . This means that RCCs are robust cuts [START_REF] Pessoa | Robust branch-cut-and-price algorithms for vehicle routing problems[END_REF], such cuts do not have any impact in the difficulty of the pricing subproblem.

Using a Chvátal-Gomory rounding of a subset C ⊆ V of Constraints (2), relaxed to ≤, with multipliers p v (0 < p v < 1), v ∈ C, the following valid Rank-1 Cut (R1C) is obtained:

m∈M P ∈Ω N m v∈C p v y P v λ P ≤ v∈C p v . (6)
The Subset Row Cuts (SRCs) proposed in [START_REF] Jepsen | Subset-row inequalities applied to the vehiclerouting problem with time windows[END_REF] are the particular R1Cs with multipliers p v = 1/K, where K is an integer, for all v ∈ C. Recently, the optimal multiplier vectors for R1Cs with up to five rows have been determined by Pecin et al. (2017c) They are optimal in the sense that they generate R1Cs that are equivalent or dominate the R1Cs generated with any other multipliers. When non-elementary routes can be generated, one can also apply cuts ( 6) with C = {v} and p v = 1/2 for some v ∈ V .

R1Cs are known to be very effective. However, because of the rounding down operator, the coefficient of a variable λ P in a R1C is not given by a linear expression over the values of x P a , so those cuts are non-robust. The concept of limited-memory for reducing the negative impact of R1Cs in the pricing, first proposed in [START_REF] Pecin | Improved branch-cut-and-price for capacitated vehicle routing[END_REF]Pecin et al. ( , 2017b)), represented a breakthrough in the performance of BCP algorithms for VRP. In this paper, we assume that the generalized arc memory variant of the concept (Pecin et al. 2017a) is being used. Each R1C, indexed by , is associated to a customer subset C ⊆ V , a multiplier vector (p v ) v∈C , and an arc memory set AM with AM ⊆ A. The idea is that, when a route P ∈ Ω N m leaves the memory set, i.e. follows an arc not in AM , it "forgets" previous visits to nodes in C , yielding a coefficient for λ P in the cut that may be smaller than the original coefficient v∈C p v y P v . The memory set AM is defined during the cut separation procedure as a minimal set preserving the coefficients of the route variables λ P that have a positive value in the current linear relaxation solution. Given C ⊆ V , a multiplier vector p of dimension |C|, and a memory arc set AM ⊆ A, the limited-arc-memory R1C cut is defined as:

m∈M P ∈Ω N m α(C, AM, p, P )λ P ≤ v∈C p v (7)
where the coefficient α(C, AM, p, P ) of the route variable λ P , P ∈ Ω N m , is computed by the following pseudocode:

Function α(C, AM , p, P = (v 0 = 0, v 1 , . . . , v k , v k+1 = n + 1)) α ← 0, S ← 0; for j = 1 to k do if (v j-1 , v j ) / ∈ AM then S ← 0; if v j ∈ C then S ← S + p vj ; if S ≥ 1 then S ← S -1, α ← α + 1;
return α;

This function can be explained as follows. Whenever a route P visits a vertex v ∈ C, the multiplier p v is added to the state variable S. When S ≥ 1, S is decremented and α is incremented. If AM = A, the function returns v∈C p v y P v and the limited-memory cut would be equivalent to the original cut. On the other hand, if AM ⊂ A, it may happen that P uses an arc not in AM when S > 0, causing the state S to be reset to zero and "forgetting" some previous visits to nodes in C. In this case, the returned coefficient may be less than the original coefficient. However, memory sets can be "corrected" in the next separation round. The final linear relaxation bound obtainable with limited memory R1Cs is exactly that potentially achieved with ordinary R1Cs.

Pricing Problem

We now define the pricing problem for the SPF over ng-routes and with additional RCCs and limited-memory R1Cs, which decomposes into subproblems, one for each vehicle type. In the subproblem for type m ∈ M , we search for a route P ∈ Ω N m with the minimum reduced cost. Let J be the current set of active RCCs (5), and L be the current set of active limited memory R1Cs (7). Each cut  ∈ J is determined by a set C  of customers. Each cut ∈ L is determined by the triple (C , AM , p ). Let (π, µ, ν, σ) be the current dual solution of the linear relaxation of (1)-( 4) with sets J and L of cuts added, where π corresponds to the partitioning constraints (2), µ corresponds to Constraints (3), ν corresponds to the set J of active RCCs, and σ corresponds to the set L of active limited memory R1Cs. Let cm (π, ν) be the vector of current arc reduced costs, where cm

(v,v ) (π, ν) = c m (v,v ) - π v 2 - π v 2 - ∈J : (v,v )∈δ + (C  )∪δ -(C  ) ν  , (8)
where π 0 and π n+1 are defined as zero. The pricing subproblem (PSP m ) for vehicle type m ∈ M is then formulated as min

P ∈Ω N m a∈P cm a (π, ν) - ∈L σ • α(C , AM , p , P ) -µ m .
3 Bucket Graph based Labeling Algorithm

SPPRC over Forward and Backward graphs

We now formulate the pricing problem defined in Section 2.3 as a SPPRC in a forward graph G = (V, A) that is identical to G, so A = A. In this section we consider the dual solution (π, µ, ν, σ) fixed and also drop index m for more clarity. Therefore, the reduced cost of an arc a ∈ A is denoted simply as ca . Define a set of two resources R = {1, 2} corresponding to vehicle capacity and time. The capacity resource consumption of each arc (v, v ) ∈ A equals to w v and the time resource consumption of the same arc equals to s v + t (v,v ) . For convenience, we use the same notation d a,r for the resource consumption of an arc a ∈ A for both resources r ∈ R. Each vertex v ∈ V has lower and upper bounds l v,r and u v,r on the accumulated resource consumption to v of each resource r ∈ R. For the second resource (time) these bounds correspond to original time windows. For the first (capacity) resource, we have l v,1 = 0 and u v,1 = W . To simplify the notation, we often drop the resource index r to represent a vector with the corresponding values for all resources in R, e.g. l v = (l v,1 , . . . , l v,|R| ) . Moreover, we use x y to denote the component-wise product of two vectors x and y. We always represent a forward path P in graph G as defined by its arcs: P = ( a P 1 , a P 2 , . . . , a P | P | ), with a P j = (v P j-1 , v P j ), v P 0 = 0, v P | P | = n + 1. For lightening the notation, we may drop the superscript P when it is clear from the context. Path P yields a vector of accumulated resource consumptions q P j ∈ IR |R| at vertices that can be computed as:

q P j = l 0 if j = 0; max q P j-1 + d aj , l vj , otherwise (9)
where the max denotes the component-wise maximum function of two vectors. A forward path P is feasible if it respects

• the resource consumption bounds:

q P j ≤ u vj , ∀j ∈ {0, . . . , | P |}; (10)
• and ng-path neighborhoods:

∀(i, j) ∈ {(i, j) : 0 ≤ i < j ≤ | P |, v i = v j = v} ∃h : i < h < j, v ∈ N v h . (11)
Condition ( 11) can be rewritten in the following way. Let F P j be the set of vertices defined as

F P 0 = ∅ and F P j = (F P j-1 ∩ N vj ) ∪ {v j } for all 1 ≤ j ≤ | P |. Then (11) is equivalent to v j ∈ F P j-1
for all 1 ≤ j ≤ | P |. Let c P (π, µ, ν, σ) be the total reduced cost of path P :

c P (π, µ, ν, σ) = a∈ P c a (π, ν) - ∈L σ • α(C , AM , p , P ) -µ.
Let also P be the set of all feasible forward paths. The pricing subproblem then can be reformulated as finding a path in P minimizing the total reduced cost: min P ∈ P c P (π, µ, ν, σ). In the remainder of this section, we use simplified notations cP and ca as the dual solution (π, µ, ν, σ) is fixed.

We now define the backward graph G = (V, A) with source n + 1 and sink 0. Set A contains one backward arc a = (v , v) for each forward arc a = (v, v ) ∈ A. Each backward arc a has the same reduced cost c a = c a and same resource consumption d a = d a as the corresponding forward arc.

A backward path P is denoted by P = ( a 1 , a 2 , . . . , a | P | ) in graph G, a j = (v j-1 , v j ), v 0 = n+1, v | P | = 0, possibly containing cycles, where we may add the superscript P if necessary as in the forward graph. The accumulated consumption vector q P j ∈ IR |R| in P is calculated recursively as

q P j = u n+1 , if j = 0; min q P j-1 -d aj , u vj , otherwise.
where the min denotes the component-wise minimum function of two vectors. Feasibility conditions for a backward path P are the same as for a forward path except that (10) is replaced by q P j ≥ l vj , for all j ∈ {0, . . . , | P |}. We denote as P the set of all feasible backward paths. Sets F P j , 0 ≤ j ≤ | P |, are defined in the same way as for a forward path.

It is easy to see that a forward path P belongs to P if and only if the corresponding backward path P belongs to P, where

| P | = | P |, a j = (v, v ) and a |P |+1-j = (v , v), 1 ≤ j ≤ | P |. Also, we have q P j ≤ q P | P |-j for 0 ≤ j ≤ | P |.
An example of a pair of forward and backward paths together with their consumption of one resource is given in Figure 1. Moreover, the total reduced costs of P and P are equal, as coefficients α(C, AM, p, P ) do not change if we traverse path P in the backward order in function α. Therefore the pricing problem defined in Section 2.3 can be also reformulated as finding a path in P minimizing the total reduced cost.

P ( )

P ( ) From now on, we use accent • for a forward entity, accent • for a backward entity, and no accent for an entity which can be both forward or backward (when it is applied). We also use accent •

q P j,1 v P 0 v P 6 v P 1 v P 5 v P 2 v P 4 v P 3 v P 3 v P 4 v P 2 v P 5 v P 1 v P 6 v P 0 0 = = n + 1 l v u v
• for an entity of the opposite sense.

Labels and Dominance Rule

We devise a labeling algorithm for the previously defined SPPRC where each label L = (c L , v L , q L , F L , S L ) corresponds to a partial forward or backward path P (partial means that we may have v

P |P | ∈ {0, n + 1}), where cL = cP , v L = v P |P | , q L = q P |P | , F L = F P |P |
, and S L ∈ IR |L| gives the current state (computed as in Function α) of each cut ∈ L for label L.

A label L dominates label L if v L = v L , q L ≤ q L (q L ≥ q L for backward labels), F L ⊆ F L , and cL -

∈L: S L >S L σ ≤ cL , (12)
as σ ≤ 0, ∈ L [START_REF] Jepsen | Subset-row inequalities applied to the vehiclerouting problem with time windows[END_REF]. Condition cL > cL is sufficient to verify that L does not dominate L .

Bucket Graph

A critical aspect of labeling algorithms that solve the previously described SPPRC is when to perform dominance checks. Given a set of labels that potentially dominate each other, dominance checks may be performed for all pairs of labels. However, it can be prohibitively time consuming. Alternatively, skipping too many dominance checks may cause a premature explosion on the number of maintained labels. One approach to address this issue is to partition the labels into buckets and to ensure that no pair of labels inside a bucket dominate each other. Pecin et al. (2017b) followed this approach by defining buckets for every possible resource consumption of labels and skipping inter-bucket dominance checks. In this paper we define buckets differently in order to avoid resource discretisation. Also, we perform inter-bucket dominance checks, but less frequently.

In the proposed labeling algorithm, labels are grouped into buckets based on their final vertices and on ranges defined for both accumulated resource consumption values. Moreover, we define bucket arcs connecting pairs of buckets through which the labels can be extended. Different bucket graphs are then defined for forward and backward labeling. Bucket graphs are useful because they help 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 the graph's strongly connected components, trying to minimize such reprocessing. Additionally, the bucket graph is used to avoid label extensions that are proved not to contribute to a solution that improves the current best one. This is achieved by removing arcs from the bucket graph based on a reduced cost argument. In what follows, we present all the notation required to formalize these ideas.

The set of buckets associated to a given vertex v ∈ V is determined by a step size dr for each resource r ∈ R. Let K v be the set of all |R|-dimensional integer vectors such that the r-th component belongs to {0, 1, . . . , (u v, r -l v, r )/ dr }, r ∈ R. A forward bucket b is defined by a pair of vertex and lower bound vector (ṽ b , l

b ) = (v, l v + κ d), v ∈ V , κ ∈ K v . A forward label L is contained in forward bucket b if v L = ṽ b and l b ≤ q L < l b + d.
Similarly, a backward bucket b is defined by a pair of vertex and upper bound vector (ṽ b , ũ Let B ( B) be the set of all forward (backward) buckets. We define as Γ ( Γ) the following set of directed forward (backward) bucket arcs.

b ) = (v, u v -κ d), v ∈ V , κ ∈ K v . A backward label L is contained in backward bucket b if v L = ṽ b and ũ b ≥ q L > ũ b -d. Let b L be the bucket containing label L.
Each forward/backward bucket arc γ ∈ Γ is defined by a pair (b γ , a γ ) of bucket and arc of the corresponding sense, b γ ∈ B, a γ ∈ A, and ṽbγ is the tail of a γ . The tail of a forward/backward bucket arc γ is bucket b Analogously, B Φ = ( B, φ( Γ) ∪ Γ Φ ) is the extended backward bucket graph. Those extended graphs are not used in the core of the labeling algorithm, they are needed only in the initialization to compute their sets of strongly connected components and appropriate topological orders for them. Let C b be the strongly connected component of a forward/backward bucket b.

γ . Set Γ contains γ if l b γ + d a γ ≤ u v , where v is the head of a γ . The head of a forward bucket arc γ ∈ Γ, denoted by b head γ , is the bucket b ∈ B such that ṽ b = v and l b ≤ max l b γ + d a γ , l v < l b + d. Set Γ contains γ if ũ b γ -d a γ ≥ l v ,
Forward/backward label L, obtained by extension from label L along a bucket arc γ with b γ = b L , is not necessary contained in bucket b head γ . For forward labels, this happens when q L r + d aγ ,r ≥ lb head γ ,r + dr , for some r ∈ R. The backward case is analogous. So, we may have b head γ ≺ b L . Even in this case, we still consider that such an extension has been made along γ. The bucket arcs in Γ Φ ensure that there is always a path from b L to b L in the extended forward bucket graph. Therefore, C b L cannot be before C b L in the topological order of strongly connected components in Γ Φ .

source v = 1 v = 2 v = 3 v = 4 sink Figure 2: An example of Forward Graph
Figure 3 shows a small extended bucket graph that corresponds to the forward graph depicted in Figure 2. In this figure, each rectangle represents the space of possible resource consumption vectors for partial paths finishing in a corresponding vertex. The consumptions of the first and second resource determine the shifts in the horizontal and vertical directions, respectively. Thus, each square is associated with two ranges for the consumptions of both resources whose dimensions are determined by the bucket steps. Each node of the extended bucket graph is depicted inside the corresponding square. The bucket arcs that belong to Γ and Γ Φ correspond to arrows between rectangles and between squares of the same rectangle, respectively. Finally, an example of a strongly connected component of this extended bucket graph is bold printed in the figure.

= 1 v = 2 v = 3 v = 4 source sink l 2,1 u 2,1 r = 1 l 2,2 u 2,2 r = 2

Mono-directional labeling algorithm

In this subsection, we describe the mono-directional labeling algorithm, which can be run in either forward or backward sense to solve the pricing problem. An auxiliary function Extend(L , γ, L) presented below extends a label L to label L along a bucket arc γ ∈ Γ. It returns false if this extension is not possible. One can easily see that the steps performed by this function initialize the reduced cost and the final vertex of the new label, calculate its resource consumptions and check their bounds, obtain its ng-route information and check its feasibility, and finally update its reduced cost considering the active R1Cs.

The labeling algorithm is presented in Algorithm 1. In this algorithm, we maintain values cbest > cL , one can conclude that no label in b or in the component-wise smaller buckets can dominate L, and this search branch can be pruned. This function also receives the set B visited of visited buckets needed to avoid processing the same bucket twice. For our case with one or two resources it is easy to do a specific implementation to verify whether a bucket belongs to set B visited in a constant time. The same implementation can be used in Functions ConcatenateLabel and UpdateBucketsSet described later. Function Extend(L , γ, L)

cL ← cL + caγ , v L ← the head of aγ q L ← q L if L is a forward label then q L ← max{q L + da γ , l v L } if q L > u v L then return false if L is a backward label then q L ← min{q L -da γ , u v L } if q L < l v L then return false if v L ∈ F L then return false F L ← (F L ∩ N v L ) ∪ {v L } for ∈ L do if aγ ∈ AM then S L ← S L else S L ← 0 if v L ∈ C then S L ← S L + p v L if S L ≥ 1 then S L ← S L -1, cL ← cL -σ return true

Bi-directional labeling algorithm

We call a pair of forward and backward labels L and L ω-compatible if

v L = v L , q L + d (v L ,v L ) ≤ q L , F L ∩ F L = ∅, and c(P L ||P L ) < ω,
where P L ||P L is the path obtained by concatenation of partial paths P L and P L along arc (v L , v L ) ∈ A, and its total reduced cost c(P L ||P L ) is calculated as

c(P L ||P L ) = c L + c(v L ,v L ) + c L - ∈L: S L +S L ≥1 σ .
For a pair of forward and backward labels L and L such that

v L = v L , c L + c(v L ,v L ) + c L is a lower bound on value c(P L ||P L ).
For the bi-directional labeling algorithm, we need to choose a resource r * ∈ R and a threshold value q * ∈ [l 0,r * , u n+1,r * ] for it. The bi-directional labeling algorithm is presented in Algorithm 2. It uses the auxiliary procedure ConcatenateLabel( L, b, P best , B visited ), presented below, which tries to find a backward label L such that b L b, b L ∈ B visited , and such that pair ( L, L) is cP best -compatible, i.e. concatenation of partial paths P L and P L along arc (v L , v L ) improves on P best . We use values cbest as bounds to prune the search. From the reasoning of the previous paragraph, if c L + c(v L ,ṽ b ) + cbest b ≥ cP best then such a label does not exist. All buckets that may contain extensions for L that improve on P best are tried. In the literature, values cbest b are referred as completion bounds [START_REF] Christofides | State-space relaxation procedures for the computation of bounds to routing problems[END_REF].

Symmetric case

In this section, we suppose that 1) all time windows are the same; and 2) R1Cs arc memories are symmetric: (v, v ) ∈ AM if and only if (v , v) ∈ AM for all ∈ L. We now show that in this case the forward/backward route symmetry can be exploited. First, we redefine the resource consumption of arcs. We make the capacity resource consumption of each arc (v, v ) ∈ A equal to 1 2 w v + 1 2 w v , and we make the time resource consumption of this arc equal to 1 2 s v + t (v,v ) + 1 2 s v . Thus, the cost and the resource consumption of two arcs (v, v ) and (v , v), v, v ∈ V , become the same. Moreover, as for any route P coefficients x P (v,v ) and x P (v ,v) are the same in constraints (2) and cuts (5), reduced costs of these arcs are the same:

Algorithm 1: Mono-directional labeling algorithm if forward algorithm then L init ← (0, 0, l0, ∅, 0) if backward algorithm then L init ← (0, n + 1, un+1, ∅, 0) insert initial label L init to its bucket b L init and mark L init as non-extended foreach strongly connected component C in B Φ in a topological order do repeat foreach bucket b ∈ C do foreach non-extended label L : b L = b do if not DominatedInCompWiseSmallerBuckets(L , b L , ∅) then foreach bucket arc γ ∈ Γ such that bγ = b L do if Extend(L , γ, L) then if L is not dominated
c(v,v ) = c(v ,v) , ∀v, v ∈ V. ( 13
)
Consider now a backward label L and the corresponding partial path P = P L . Let P be the partial forward path such that v P j = v P j , 1 ≤ j ≤ | P |, and L be the forward label such that P L = P . From ( 13) and from the fact that the R1Cs arc memories are symmetric, we have c L = c L , F L = F L , and S L = S L . As time windows are the same, resource bounds are also the same: [l v , u v ] = [l, u] for all v ∈ V . Then, q L -l = u -q L . Moreover, for every backward bucket b, there exists a symmetric forward bucket b such that l b -l = u -ũ b and κ b = κ b . In this case (which we call symmetric), in the bi-directional labeling algorithm, we set q * = 1 2 l 0,r * + 1 2 u n+1,r * , skip the backward labeling step, and use symmetric forward buckets and labels instead of backward buckets and labels in the concatenation step.

Algorithm 2: Bi-directional labeling algorithm run forward labeling algorithm where only labels L, q L r * ≤ q * are kept run backward labeling algorithm where only labels L, q L r * > q * , are kept let P best be the best complete path obtained in the two algorithms above

foreach forward label L do foreach bucket arc γ ∈ Γ such that b γ = b L do if Extend( L, γ, L ) and q L r * > q * then let b ∈ B be the bucket such that ũ b ≥ q L > ũ b - d ConcatenateLabel( L, b, P best , ∅) return P best Procedure ConcatenateLabel( L, b, P best , B visited ) B visited ← B visited ∪ { b} if c L + c(v L ,ṽ b ) + cbest b ≥ cP best then return foreach label L : b L = b do if pair ( L, L) is cP best -compatible then P best ← P L , a = (v L , v L ), P L foreach b ∈ Φ b \ B visited do ConcatenateLabel( L, b , P best , B visited ) return

Bucket arc elimination

Arc elimination procedure by [START_REF] Irnich | Path-reduced costs for eliminating arcs in routing and scheduling[END_REF] can be employed to remove arcs from the original graph using the argument that any path using one of these arcs is not contributing to an optimal solution. In this section, we generalize this procedure to eliminate bucket arcs from the bucket graph using the same argument. However, this generalization is not straightforward because of the dominance between labels. Even labels from the same bucket may use different bucket arcs when extended through the same path in the original graph. Hence, the same path extension may be feasible for the dominated label and not feasible for the dominating label due to removal of some bucket arcs. We introduce below the new concept of jump bucket arcs to bring back compatibility between bucket arc elimination and dominance.

We call a path P ∈ Ω m , m ∈ M , improving if λ P participates in an integer solution of the original problem (SPF) with a cost smaller than the cost of the best solution found so far. We denote the latter cost as U B, as it is an upper bound on the optimal solution value of (SPF).

Let Ω I m be the set of improving paths in Ω m . We can exclude paths proved to be non-improving from the solution space of the pricing problem. Note that if there are no improving paths (i.e., the best known solution is optimal) any path can be correctly excluded, so we only need to be careful in the case where some Ω I m sets are not empty. So, we assume that case in the remainder of the section. Suppose that some variables λ P , P ∈ Ω I m , are fixed to zero in (SPF). Let (π, μ, ν, σ) be the dual solution of the current restricted linear relaxation of (SPF), corresponding to constraints (2), (3) and to constraints of type ( 5) and ( 7). Let z (RSPF) (π, μ, ν, σ) be the value of this solution. The value

LB Lagr = z (RSPF) (π, μ, ν, σ) + m∈M min 0, min P ∈Ω I m cP (π, μ, ν, σ) • U m
is a Lagrangian lower bound on the optimal solution value of (SP F ). It is obtained by dualizing constraints (2), ( 5) and ( 7) with multipliers π, ν, and σ. For any vector (z 1 , z 2 , . . . , z |M | ) such that z m ≤ min

P ∈Ω I m cP (π, μ, ν, σ), m ∈ M, (14)
we have

LB = z (RSPF) (π, μ, ν, σ) + m∈M min {0, z m } • U m ≤ LB Lagr . (15)
Thus, LB is a valid lower bound on the optimum solution of (SPF). If an algorithm for solving pricing sub-problem (PSP m ) returns value z m satisfying condition ( 14), we say that this algorithm satisfies ( 14). Such an algorithm can be used to calculate a valid lower bound for (SPF) even if it does not necessarily solve (PSP m ) to optimality. In the remainder of this section, we again drop index m and consider dual solution (π, μ, ν, σ) fixed for more clarity.

We say that a feasible forward/backward path P ∈ P passes through bucket arc γ ∈ Γ if there is an index j, 1 ≤ j ≤ |P |, such that a P j = a γ , and lbγ ≤ q P j-1 < lbγ + d for a forward path (ũ bγ ≥ q P j-1 > ũbγ -d for a backward path). We denote such a bucket arc as γ P j . We call a forward/backward γ ∈ Γ an improving bucket arc if there exists an improving path of the same sense passing through it. The bucket arc elimination procedure described in Section 4.2 reduces the current set Γ by removing some non-improving bucket arcs from it. However, the labeling algorithm proposed in Section 3, in which labels are only extended over a reduced set Γ of bucket arcs does not necessarily satisfy condition (14).

Figure 4 illustrates a case where the value z of the best solution found by the labeling algorithm does not satisfy ( 14). In this figure, the value of a single resource consumption r = 1 for each partial path is represented by the level of its end node in the vertical axis. P and P are improving and non-improving forward paths, respectively. Thus, for each traversed arc, the corresponding bucket arc that is passed through depends solely on the total resource consumed up to this point. Paths P and P traverse different arcs until vertex v P i-1 = v P i -1 and finish with the same subpaths from this point until the sink vertex (but with such subpaths passing through different bucket arcs).

Assume that the reduced bucket graph does not contain the non-improving bucket arc γ corresponding to (v P i -1 , v P i ). Suppose also that P has a reduced cost smaller than that of P , allowing the label L associated to P until vertex v P i -1 to dominate the label L associated to P until the same vertex. The reasoning for such domination is that one can always replace the partial path of L by that of L, causing that P is replaced by P which has a smaller reduced cost. However, absence of bucket arc γ in the reduced bucket graph prevents P from being found by the labeling algorithm. Thus, neither P nor P are obtained by the labeling algorithm, and the value z found by the algorithm may violate condition (14). The remainder of this section is divided in two parts. In Section 4.1, we show how to modify our labeling algorithm so that it satisfies condition ( 14) when executed over a superset of all improving bucket arcs. In Section 4.2, we present a procedure for finding non-improving bucket arcs, and show the correcteness of its combination with the modified labeling algorithm.

q P 1 s l0,1 l1,1 l2,1 l3,1 l4,1 v P i -1 v P i-1 v P i v P i v P i +1 v P i+1 v P |P | v P |P | L L domination P P ∈ Γ

Modified labeling algorithm

To guarantee the fulfillment of condition ( 14), labels in our forward (backward) labeling algorithm should not only be extended along the current reduced set Γ ( Γ) of bucket arcs, but also along the set Ψ ( Ψ) of jump bucket arcs obtained as shown below. Jump bucket arcs impose an additional resource consumption when traversed. , and a ψ = (ṽ, v ) ∈ A is its arc of the corresponding sense. In the labeling algorithm, each label L should be extended along a jump bucket arc ψ ∈ Ψ if b L = b base ψ . When extending label L along ψ ∈ Ψ, we first do the "jump", i.e. in the second line of function Extend(L ,ψ,L), instead of q L ← q L , we do q L ← max q L , lb jump ψ (q L ← min q L , ũb jump ψ for a backward label). Function ObtainJumpBucketArcs(Γ) presented below obtains the current set Ψ of jump bucket arcs given the current reduced set Γ of bucket arcs. This function can be used for both forward and backward sense. We denote as L P j the label corresponding to the partial path consisting of the first j arcs of path P . If P passes through γ, then there is an index j, 1 ≤ j ≤ |P |, such that label L P j is contained in bucket b γ . Let L(Γ, Ψ) and L(Γ, Ψ) be the sets of non-dominated labels generated by the modified forward and backward mono-directional labeling algorithms, in which each label L is extended along all bucket arcs γ ∈ Γ such that b γ = b L and along all jump bucket arcs

Function ObtainJumpBucketArcs(Γ) Ψ ← ∅ foreach bucket b ∈ B do foreach a = (v, v ) ∈ A such that v = ṽb do if ∃γ ∈ Γ such that bγ = b, aγ = a then B ←
ψ ∈ Ψ such that b base ψ = b L .
The next proposition is auxiliary but very important. It shows that for every improving path P , the modified mono-directional labeling algorithm either generates P or a path dominating P . Moreover, this is true for any prefix of P . This proposition applies for both forward and backward sense. Its proof is illustrated in Figure 4. Proposition 1. Let set Γ be a superset of all improving bucket arcs, and Ψ be the set of jump bucket arcs computed using procedure ObtainJumpBucketArcs(Γ). Then, for each improving path P ∈ P and for each index j, 0 ≤ j ≤ |P |, either L P j ∈ L or there exists a label L ∈ L dominating L P j , where L = L(Γ, Ψ).

Proof. Consider an improving path P ∈ P. We prove this proposition by induction on the value of the index j, 0 ≤ j ≤ |P |. For j = 0 the proposition is true since there is a single label L init in the bucket b L init that corresponds to vertex v P 0 . Now, assume that the proposition is true for

j = i -1, 1 ≤ i ≤ |P |.
We must prove that it is also valid for j = i. By the inductive hypothesis, either L P i-1 ∈ L or there exists a label L ∈ L dominating L P i-1 . Let us use the notation L defined for the latter case to also denote L P i-1 in the former case. In all cases, b L b L P i-1 . Let γ and γ be the bucket arcs such that b

γ = b L , b γ = b L P i-1
, and a γ = a γ = a P i . If γ ∈ Γ, the label L obtained by extending L through a γ dominates L P i . In this case, either L ∈ L or L is dominated by some label from L, which should also dominate L P i . In both situations, the proof is finished. Otherwise, since γ ∈ Γ, γ is an improving arc, and b L ≺ b γ , there exists a jump bucket arc ψ ∈ Ψ such that a ψ = a . Existence of such a jump bucket arcs is guaranteed by the definition of procedure ObtainJumpBucketArcs(Γ). Thus, extending L through a ψ leads to a label L that dominates L P i . As before, either L ∈ L or L is dominated by some label from L, which dominates L P i , finishing the proof.

We have just shown essentially that the modified mono-directional labeling algorithm satisfies condition ( 14). In Proposition 3, we show that the same condition is satisfied by the modified bidirectional labeling algorithm. First observe that any complete path can be divided into forward and backward parts by removing any arc it passes through. Moreover, the total reduced cost of the concatenation of these parts is equal to the reduced cost of the complete path.

Observation 1. For a pair of corresponding forward/backward path P and the opposite sense path P P and each value j, 0 ≤ j < |P |, pair of labels L P j and L L P P |P |-j-1 is ω-compatible for all ω > cP .

Secondly, we will need the following auxiliary proposition. Proposition 2. Suppose that i) forward and backward labels L and L are ω-compatible, ii) either L = L or label L dominates L , iii) either L = L or label L dominates L . Then labels L and L are also ω-compatible.

Proof. Let a = (v L , v L ) = (v L , v L ).
From the definition of dominance in Section 3.2 and the definition of ω-compatible labels in Section 3.5, we have 12) and the definition of c(P L ||P L ) it follows

v L = v L = v L = v L , q L + d a ≤ q L + d a ≤ q L ≤ q L . We also have F L ⊆ F L , F L ⊆ F L , and F L ∩ F L = ∅. It follows then that F L ∩ F L = ∅. From (
c(P L ||P L ) = c L + ca + c L - ∈L: S L +S L ≥1 σ ≤ c L + ∈L: S L >S L σ ≥c L as L dominates L +c a + c L + ∈L: S L >S L σ ≥c L as L dominates L - ∈L: S L +S L ≥1 σ = c L + ca + c L - ∈L: S L +S L ≥1 σ =c(P L ||P L )<ω + ∈L: S L +S L ≥1, S L +S L <1 σ ≤0 as σ≤0 + ∈L: S L >S L σ + ∈L: S L >S L σ - ∈L: S L +S L <1, S L +S L ≥1 σ ≤0 as σ≤0
It follows then that c(P L ||P L ) < ω. Hence, all conditions for L and L to be ω-compatible are met.

We define as z( Γ, Γ, Ψ, Ψ) the reduced cost of the best path obtained by the modified bidirectional labeling algorithm in which forward labels are extended along bucket arcs in Γ and Ψ, and backward labels are extended along bucket arcs in Γ and Ψ.

Proposition 3. Let sets Γ and Γ be supersets of all improving forward and backward bucket arcs, and Ψ and Ψ be the sets of forward and backward jump bucket arcs computed using procedure ObtainJumpBucketArcs(Γ). Then z ≤ c P for any improving path P ∈ P, where z = z( Γ, Γ, Ψ, Ψ).

Proof. Consider an arbitrary improving path P ∈ P. Let j, 0 ≤ j ≤ | P | be the largest index such that q L P j r * ≤ q * . We divide the proof into two cases: j = | P | and j < | P |.

For j = | P |, let L = L P j . We have v L = n + 1. As P is improving, by Proposition 1, either L ∈ L( Γ, Ψ) or there is a label L ∈ L( Γ, Ψ) dominating L. As v L = n + 1, path P L ∈ P and the total reduced cost of P L is ≤ c P , thus z ≤ c P . Suppose now that j < | P |. Let P be the corresponding backward path for P . By Observation 1, the pair of labels L = L P j and L = L P | P |-j-1 is ω-compatible for all ω > c P . As both P and P are improving, by Proposition 1 either L ∈ L( Γ, Ψ) or there exists a label L ∈ L( Γ, Ψ) dominating L . Also, by Proposition 1 either L ∈ L( Γ, Ψ) or there exists a label L ∈ L( Γ, Ψ) dominating L . Thus, by Proposition 2 there exists a pair of ω-compatible labels generated by the modified labeling algorithm for all ω > c P . Then c( P L || P L ) ≤ c P , thus z ≤ c P .

We have just shown that our labelling algorithm can be modified in order to satisfy condition ( 14) by adding jump bucket arcs. Now it can be executed over any superset Γ of improving bucket arcs. In the next section, we show how we can find non-improving bucket arcs, which then can be safely removed from the current set Γ.

A procedure for finding non-improving bucket arcs

In this paragraph, we temporarily consider again all subproblems (PSP m ), m ∈ M . Let (π, μ, ν, σ) be the current dual solution of the restricted linear relaxation of (SPF) together with constraints ( 5) and ( 7). Suppose that the current sets Γ m , m ∈ M , include all improving buckets arcs, and the modified labelling algorithms have returned values z m for all m ∈ M . By Proposition 3, the modified labeling algorithm satisfies condition ( 14). Thus the value LB given by ( 15) is a valid lower bound for (SPF). A necessary condition for a path P ∈ Ω m to be improving is LB -min {0, z m } + cP (π, μ, ν, σ) < U B. We denote

θ m = U B -LB + min {0, z m } . (16)
Then, a sufficient condition for a forward/backward path P ∈ Ω m to be non-improving is cP (π, μ, ν, σ) ≥ θ m . The next proposition gives a necessary condition for a bucket arc to be improving. The proposition applies for both forward and backward sense. In the following, we drop again index m. Proof. If γ is an improving bucket arc, then there exists an improving path P and an index j, 0 ≤ j < |P |, such that b γ = b L P j and a γ = a P j+1 . We have cP < θ, as it is a necessary condition for P to be improving. Then by Observation 1, labels L Proposition 4 is the basis for the procedure to find non-improving bucket arcs. To prove that a bucket arc γ is non-improving, we need to verify that there does not exist a pair of θcompatible labels L ∈ L and L L ∈ L L (one at the tail and one at the head of a γ ) such that b L b γ and b b L L b b arr γ . A direct approach is to exhaustively enumerate all such pairs of labels for every bucket arc and verify whether at least one of them is θ-compatible. This approach is slow as the same pair of labels may be verified for θ-compatibility several times.

To guarantee that no pair of labels is verified for θ-compatibility more than once, we compute

B B visited ← B B visited ∪ { b b} if cL + cv L ,ṽ b b + cbest b b ≥ θ then return if b b ∈ B B then foreach label L L ∈ L L, b b L L = b b do if pair (L, L L) is θ-compatible then B B ← B B ∪ { b b} foreach b b ∈ Φ b b \ B B visited do UpdateBucketsSet(θ, L, B B, b b , L L, B B visited )
Finally, we present formally the main procedure BucketArcElimination(θ, Γ, Ψ, L, L L) which can be used for both forward and backward sense. The procedure removes some non-improving bucket arcs from Γ, given the set L of non-dominated labels of the same sense as Γ and the set L L of non-dominated labels of the opposite sense. In this procedure, in order to do the correct initialization of sets B B a,b , buckets b are considered in a lexicographic order of κ b . When considering bucket b, procedure UpdateBucketsSet is called to augment set B B a,b only if there exists γ ∈ Γ such that b = b γ or there exists ψ ∈ Ψ such that b = b base ψ . If such a (jump) bucket arc does not exist, there is no bucket arc γ ∈ Γ such that a γ = a and b γ b, and thus set B B a,b is not used in the reminder of the algorithm.

The next proposition proves formally that the combination of the modified bi-directional labeling algorithm and the bucket arc elimination procedure is correct, i.e. it allows us to compute a valid lower bound for the original problem.

Proposition 5. Let z be the reduced cost of the best path obtained by the modified bi-directional labeling algorithm after performing j calls to procedure BucketArcElimination, for some integer j ≥ 0. Then z ≤ c P for any improving path P ∈ P.

Proof. It is enough to prove that the current forward (backward) bucket graph obtained after j calls of Procedure BucketArcElimination is a superset of all improving forward (backward) bucket arcs, since this proposition is a direct consequence of this statement and Proposition 3. We prove it by induction on the value of j. For j = 0, this is true because the initial forward and backward bucket graphs contain all possible forward and backward bucket arcs. Now, suppose by inductive hypothesis that this is valid for j = i, i.e. all improving forward Procedure BucketArcElimination(θ, Γ, Ψ, L, L L) Since this pair of labels prevents γ to be removed in the next call to procedure BucketArcElimination, the statement under consideration is also true for j = i + 1, finishing the proof.

foreach bucket b ∈ B in a lexicographic order of κ b do if Φ b = ∅ then foreach arc a = (ṽ b , v ) ∈ A do B B a,b ← ∪ b ∈Φ b B B a,b else foreach arc a = (ṽ b , v ) ∈ A do B B a,b ← ∅ foreach jump bucket arc ψ ∈ Ψ such that b base ψ = b do foreach L ∈ L such that b L = b do UpdateBucketsSet(θ, L, B B a ψ ,b , b b arr ψ , L L, ∅) foreach bucket arc γ ∈ Γ such that bγ = b do foreach L ∈ L such that b L = b do UpdateBucketsSet(θ, L, B B aγ ,b , b b arr γ , L L, ∅) if ∃ b b ∈ B B aγ
Let us now give some final remarks on our bucket arc elimination procedure.

• Our procedure can be viewed as a "resource-value dependent" arc elimination and thus a generalization of the arc elimination procedure by [START_REF] Irnich | Path-reduced costs for eliminating arcs in routing and scheduling[END_REF]. The latter can readily be extended to remove extensions of forward labels L along an arc if q L exceeds a certain threshold (and extensions of backward labels L if q L is lower than a threshold).

The opposite case ("less than" for the forward sense and "greater than" for the backward sense) is more complicated as shown in our analysis and necessitates introduction of jump bucket arcs. Our procedure is more generic than both, and may result in elimination of disjoint "resource-value" intervals for an arc.

• Necessity of adding jump bucket arcs essentially comes from the fact that the dominating label may be contained in a different bucket than the dominated one. Restricting dominance checks only to ones within buckets does not make jump bucket arcs redundant. This is because two labels in the same bucket with distinct resource consumption may be extended (over the same arc) to labels in different buckets. Only restricting the dominance checks to labels with the same resource consumption allows one to skip the generation of jump bucket arcs.

• "Resource-value dependent" arc elimination has already been used in [START_REF] Pessoa | Exact algorithm over an arc-time-indexed formulation for parallel machine scheduling problems[END_REF] and in Pecin et al. (2017b). However, it was done at the cost of resource discretization, limited dominance checks, and dropping the exploitation of the route symmetry. Our procedure allows one to overcome these drawbacks.

5 Branch-Cut-and-Price algorithm

We implemented a BCP algorithm for the HFVRPTW using the new bucket graph labeling algorithm in its pricing and also the bucket arc elimination procedure. Other algorithmic elements are similar to the ones presented in [START_REF] Pecin | Improved branch-cut-and-price for capacitated vehicle routing[END_REF], Pecin et al. (2017b) and Pecin et al. (2017a), except that a different dynamic ng-relaxation by [START_REF] Bulhoes | A branch-and-price algorithm for the minimum latency problem[END_REF] and automatic dual price smoothing stabilization by [START_REF] Pessoa | Automation and combination of linear-programming based stabilization techniques in column generation[END_REF] are used. However, to make this paper self-contained, we give a succinct algorithm description in what follows.

For each node of the branch-and-bound tree, the lower bound on the optimal cost is computed by solving the linear relaxation of SPF enhanced with RCCs (5) and limited-memory R1Cs (7). Before each cut generation round, the current SPF relaxation is solved through column generation, using automatic (parameterless) dual pricing smoothing stabilization described in [START_REF] Pessoa | Automation and combination of linear-programming based stabilization techniques in column generation[END_REF]. To further speed up the convergence, three-stage column generation is implemented. At the first stage, the "light" pricing heuristic is applied in which only one label per bucket (with the smallest reduced cost) is kept. In the second stage, a more expensive pricing heuristic is used in which states F and S are not taken into account in the dominance checks. In the last stage, the exact labeling algorithm is used. We use parameters χ heur and χ exact for the maximum number of columns generated at each iteration of column generation at heuristic and exact stages. When the number of columns in the restricted master exceeds χ max , we clean them up and leave only χ perc % columns with the least reduced cost. Basic columns are never removed, however, non-basic columns with zero reduced cost may be removed.

The threshold value q * for the bi-directional labeling algorithm is an important parameter in the non-symmetric case. Pecin et al. (2017b) and [START_REF] Tilk | Asymmetry matters: Dynamic half-way points in bidirectional labeling for solving shortest path problems with resource constraints faster[END_REF] proposed to modify the labeling algorithm to automatically adjust this value. Here we adopt a simpler but still effective approach. Value q * is initialized as the average middle value for all resource consumption intervals. After each exact pricing we compare the number of forward and backward nondominated labels. If one number exceeds another by more than 20%, we adjust value q * by moving it by 5% towards the maximum or minimum possible value.

The limited-arc-memory R1Cs are separated by complete enumeration for |C| ≤ 3 and by a local search heuristic for |C| ∈ {4, 5}. This heuristic is launched for each possible optimal vector p. Given p, it tries to find a subset C of customers such that the corresponding cut is violated with a full memory. If such C exists, an arc-memory AM of minimal size such that the cut violation remains the same is identified. Then AM is completed with all arcs between vertices in C. When there already exists a cut defined for the same subset C and vector p but for a different arc-memory AM , this cut is simply enhanced by enlarging its memory to AM ∪ AM . Otherwise, a new cut is added for C, p and AM . Also, memory sets are kept symmetric. This means that if arc (v, v ) is added to AM , then (v , v) must also be added. Parameter vector β is used in which β k , k ∈ {0, 1, 3, 4, 5} is the maximum number of R1C with |C| = k separated per round (β 0 is the maximum number of RCCs generated per round). Other parameters are δ perc % and δ num which are used to stop the cut generation by tailing-off. If the number of times the dual-primal gap is decreased by less than δ perc % reaches δ num , then branching is performed. At the root node, R1Cs are separated only if RCCs are exhausted or tail off. Before each cut round, non-active cuts are removed from the restricted master.

In our algorithm, we use the dynamic variant of the ng-relaxation from Bulhoes, Sadykov, and Uchoa (2018) which is inspired by [START_REF] Roberti | Dynamic ng-path relaxation for the delivery man problem[END_REF]. We start by generating ng-neighborhood for each customer by including their η init closest customers (including itself). After the convergence of column generation, we augment neighborhoods in order to forbid the cycles in the columns forming the current fractional solution. For each such column we either forbid all cycles of size at most η 1 or, if there is no one, the minimum size cycle. Cycles in at most η 2 columns are forbidden in each round. After augmentation, the neighborhood of any customer cannot exceed size η max . The same tailing-off parameters are applied here as for cutting planes. We start the root node by augmenting the ng-neighborhoods. Then, after tailing-off, RCCs are separated. After another tailing-off, R1Cs are separated.

Each additional active inequality (7) or ng-neighborhood increase make the pricing problem harder, since they weaken the dominance conditions. Because of that, separation is stopped when the average exact pricing time since the last round of cuts has exceeded a given threshold τ . In some extreme cases, R1Cs may be even removed by the roll-back procedure activated if a single pricing time exceeds τ , as described in Pecin et al. (2017b).

The bucket arc elimination procedure is called after the initial convergence, and also each time the current primal-dual gap decreased by at least υ% since the last call to the bucket arc elimination procedure. A bi-directional enumeration procedure is called after each reduced cost fixing to try to generate all improving elementary routes, using the algorithm proposed in Pecin et al. (2017b). The enumeration procedure is interrupted if the number of labels exceeds ω labels or if the number of enumerated routes exceeds ω routes . After a successful enumeration for a pricing subproblem corresponding to a vehicle type, this subproblem passes to the enumerated state and the (SPF) relaxation is updated by excluding non-elementary columns coming from this subproblem. The pricing subproblem in the enumerated state is solved by inspection. Elimination procedure based on reduced costs is also performed for enumerated pricing subproblems, but considering each route individually. If after some point, the number of enumerated routes for each subproblem falls below ω MIP , the residual problem is handled by a standard MIP solver, and the node is considered as processed.

Recall that to perform bucket arc elimination, the forward and backward passes of the labeling algorithm must be performed completely. Thus the running time may increase substantially in comparison with the "standard" exact labeling in which the labels beyond the bi-directional threshold value q * are not stored. To accelerate the bucket arc elimination procedure, we use what we call exact completion bounds technique: for every label L just created by extension beyond q * , we verify exhaustively whether there exists an opposite sense label L L such that L and L L are θ-compatible. Here θ is as defined in Section 4. We keep label L only if such a label L L exists. A similar approach is used in the bi-directional enumeration procedure. Using completion bounds for labels before the bi-directional threshold value q * , as it is done in Pecin et al. (2017b), was not computationally beneficial for our labeling algorithm.

If tailing-off condition is attained, or at least one non-enumerated pricing subproblem becomes too time consuming to solve, branching must be performed. This is done by adding constraints to the master LP that correspond to tightening lower and upper bounds on the value of the following aggregated variables corresponding to different branching strategies:

• g m = P ∈Ωm v ∈V 1 2 x P (0,v ) + x P (v ,n+1) λ P , m ∈ M ; • g m v = P ∈Ωm y P v λ P , v ∈ V , m ∈ M ; • g {v,v } = m∈M P ∈Ωm (x P (v,v ) + x P (v ,v) )λ P , v, v ∈ V , v = v .
g m corresponds to the number of used vehicles of type m ∈ M ; g m v characterizes whether customer v ∈ V is served by a vehicle of type m; g {v,v } characterizes whether edge (v, v ) is used in the solution. Given a fractional branching variable g with fractional value f , two child nodes are created by adding the constraints g ≤ f and g ≥ f to the master LP. A generic branching constraint m∈M P ∈Ωm a∈A α m a x P a λ P ≥ β added to the master LP with associated dual variable ι can be considered by the pricing algorithm by subtracting ια m a from the expression of cm a (π, ν) given by ( 8), for each m ∈ M and a ∈ A, thus, not making the subproblem harder. Moreover, the symmetry is preserved, as for any route P , coefficients x P (v ,v) and x P (v,v ) are the same in all branching constraints we use. The selection of the branching variable at each node is done using a sophisticated hierarchical evaluation strategy similar to the one proposed in Pecin et al. (2017b). The idea is to spend more time evaluating branching variables in the lowest levels of the branch-and-bound 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: 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 the three strategies. Within the same strategy, the candidates are chosen based on the distance from its fractional value to the closest integer and for branching on edges, distance to the closest depot (the smaller the better) is also taken into account.

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 heuristic column generation, but not cut generation. The best candidate is also selected by the Product Rule.

The parameters ζ 0 and ζ 1 are used to specify the maximum number of candidates chosen in phases 0 and 1.

Computational experiments

The BCP algorithm described in the previous section was coded in C++ and compiled with GCC 5.3.0. The BaPCod package by [START_REF] Vanderbeck | BaPCod -a generic Branch-And-Price Code[END_REF] was used to handle the BCP framework. IBM CPLEX Optimizer version 12.6.0 was used as the LP solver in column generation and as the IP solver for the set partitioning problem with enumerated columns. 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 eight copies of the algorithm running in parallel on the server. Each instance is solved by one copy of the algorithm using a single thread.

Instances

VRPTW instances For the vehicle routing problem with time windows, 14 most difficult (according to Pecin et al. (2017a)) instances with 100 customers by [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF], as well as 120 instances by [START_REF] Gehring | Parallelization of a two-phase metaheuristic for routing problems with time windows[END_REF] with 200 and 400 customers were considered. Solomon instances used are C203, C204, RC204, RC207, RC208, R202, R203, R204, R206, R207, R208, R209, R210, and R211. The instances are divided into classes of clustered (C), random-clustered(RC), and random (R) location of customers on the plane. Note that the capacity resource in Solomon instances, as well as Gehring and Homberger instances of classes C2, RC2, and R2, is not tight. Therefore, the subproblem only considers the time resource in those classes. The capacity resource is imposed through RCCs, which are treated as "core" cuts, i.e. they are separated for every fractional but also for every integer solution found during the search. The initial upper bounds were obtained by the heuristic in [START_REF] Vidal | A hybrid genetic algorithm with adaptive diversity management for a large class of vehicle routing problems with time-windows[END_REF]. The same values were used in Pecin et al. (2017a).

DCVRP instances For the distance constrained vehicle routing problem, we used seven classic instances from [START_REF] Christofides | The vehicle routing problem[END_REF] with 50-200 customers, named CMT6, CMT7, CMT8, CMT9, CMT10, CMT13, and CMT14. In spite of the fact that those instances are ubiquitous in the heuristic literature for the CVRP, no exact methods could solve them. The best known solutions for these instances are taken from [START_REF] Nagata | Edge assembly-based memetic algorithm for the capacitated vehicle routing problem[END_REF].

MDVRP instances For the multi-depot vehicle routing problem, we used 22 standard distance constrained instances by Cordeau, Gendreau, and Laporte (1997) with 80-360 customers, named p08-p11, p13, p14, p15, p17, p19, p20, p22, p23, and pr01-pr10. The best known solutions are taken from [START_REF] Vidal | A hybrid genetic algorithm for multidepot and periodic vehicle routing problems[END_REF].

SDVRP instances For the site-dependent vehicle routing problem, we used 10 standard distance constrained instances by [START_REF] Cordeau | A tabu search algorithm for the site dependent vehicle routing problem with time windows[END_REF] with 48-288 customers, named pr01-pr10. The best known solutions are taken from [START_REF] Cordeau | A parallel iterated tabu search heuristic for vehicle routing problems[END_REF]. et al. (2019). Following the convention used in those papers, we impose that each customer should be visited exactly once. As the triangle inequality is not always satisfied in those instances, allowing multiple visits could possibly reduce costs.

Algorithm parameterization

The column management parameter values are χ heur = 30, χ exact = 150, χ max = 10 000, χ perc = 66%. The cut generation parameter values are The pricing time thresholds (in seconds) depend on instance class, instance size, and on the number of resources (dimension of buckets). These thresholds as well as which resource is chosen as r * in bidirectional labeling are shown in Table 1.

β 0 = β 1 = β 4 = β 5 =
The bucket step size in the labeling algorithm is adjusted dynamically as explained below.

Impact of bucket step size and improved arc elimination

As two main contributions of the paper are using the bucket graph in the labeling algorithm and improved procedure for bucket (or resource value dependent) arc elimination, we first estimate their impact on the solution time of the BCP algorithm.

These tests were done on the set of 74 VRPTW instances with 100 and 200 customers. Only the root node is considered. For instances which cannot be solved at the root, we establish the "target" lower bounds. These values are determined in such a way that once they are achieved by a good parameterisation of the algorithm in the root node, it is reasonable to branch because either i) cuts are tailing off, or ii) pricing becomes too expensive. In Table 2, we give the instance name, its optimum or Best Known Solution (BKS) value, and the Target Lower Bound (TLB) value we use. For instances not shown in this table, the target lower bound equals the optimum solution value.

The first parameter we test is the "maximum number of buckets per vertex" which we denote as ξ. For instances with one time resource, bucket step size d is determined as (u n+1 -l 0 )/ξ. For instances with two resources, bucket step size d1 for the first (capacity) resource is determined as W/ √ ξ, and bucket step size d2 for the second (time) resource is determined as (u n+1 -l 0 )/ √ ξ. All values dr are rounded off to one decimal place.

The second binary parameter to test is whether the standard arc elimination procedure (as in [START_REF] Irnich | Path-reduced costs for eliminating arcs in routing and scheduling[END_REF]) or improved bucket arc elimination procedure is applied. In the former case, we fix arc a ∈ A if all bucket arcs γ ∈ Γ such that a γ = a can be fixed. When fixing original arcs, no jump bucket arcs are needed. The procedure for fixing original arcs runs faster, as once 2010) is valid only for the single resource case with positive integer consumption. This procedure can be viewed as a fixing in a bucket graph with unitary step size ( d = 1) under the condition that the dominance checks are performed only between labels with the same resource consumption.

The base variant of the algorithm which we denote as (O) is parameterized with ξ = 1, which corresponds to the basic variant of the label correcting algorithm, and uses the standard arc elimination procedure. The base variant is compared with variant (A) which is parameterized with a fixed value ξ ∈ Ξ = {1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000} and uses the standard arc elimination procedure. Then variant (A) is compared with variant (A + ) which uses a procedure described below to automatically adjust value ξ and the standard arc elimination procedure. Finally variant (A + ) is compared with variant (A + B) which uses dynamically adjusted value ξ and the improved procedure for bucket arc elimination.

Different bucket step sizes used in the labeling algorithm usually change the set of routes returned by the "light" pricing heuristic to the master. This happens because a given nondominated label may be extended if one bucket step size is used and not extended if another bucket step size is used. Thus, for each value ξ, the sequence of primal-dual master solutions in the column generation may follow a different trajectory and converge to distinct alternative optima. In turn, this changes the set of cuts returned by the separation procedures, affecting the pricing difficulty for the next iterations. This "random noise" may make the experimental results less clear. To avoid this, for each instance i ∈ I we choose a value ξi ∈ Ξ. Then, two copies of the labeling algorithm are created, the first one with value ξi , and the second one with the value being tested ξ ∈ Ξ. Both copies are executed in each iteration of columns generation, and the arc elimination procedure is also performed twice. Only the set of columns produced by the first copy of the algorithm is sent to the master. Therefore, the sequence of primal-dual solutions of the master problem is always the same regardless of the bucket step size parameterisation.

For each run, i.e. for each pair (i, ξ ), we collect the following statistics. The first one is the ξ-subproblem time which includes the time to create the bucket graph and a topological order for it, the total pricing time, the total bucket arc elimination time, and the total enumeration time. The second statistics is the master time which is the difference between the overall running time and the sum of ξi -subproblem time and ξ -subproblem time. The total time is the sum of the master time and the ξ -subproblem time. For each instance i ∈ I and each value ξ ∈ Ξ, we calculate ratios tot i,ξ and sp i,ξ between the corresponding total time (or ξ -subproblem time) and the minimum total time (minimum ξ-subproblem time) among all ξ ∈ Ξ for this instance. Ratio tot i,1 shows the speed-up of variant (A) of the algorithm over variant (O) for instance i ∈ I.

In Figure 5, we show the graphs for average and maximum values for both tot i,ξ and sp i,ξ , as well as graphs for these ratios for some representative individual instances. The averages are shown for two groups of instances: with one and two resources. The time limit was reached for ξ = 1, 2, 5, 10, 20, 50, 100 for 6, 6, 4, 4, 3, 1, 1 instances accordingly. Therefore average and maximum values given in Figure 5 for these values of ξ are in fact their lower bounds. It can be seen from the figure that parameter ξ has a large impact on the total pricing time, as well as on the overall time of the root node. The best value for this parameter is instance-dependent. Two extreme values ξ = 1 and ξ = 2000 are bad parameterisations on average. The variant (A) with parameterisation ξ = 200 is at least 2.85 times faster on average than variant (O). The maximum speed-up reaches at least the factor of 20 (it is expected to be much higher for instances stopped by the time limit). In Figure 5, the graph for instance RC 1 2 5 also shows the interest to have smaller buckets than it is required to make the extended bucket graph acyclic. For this instance, setting parameter ξ to 500 makes the extended bucket graph acyclic. However, the subproblem time and the total time of the algorithm are, respectively, 14% and 12% smaller for parameterisation ξ = 1000. Most of the time, however, the best parameterisation results in an extended bucket graph with cycles.

In variants (A + ) and (A + B) we use a scheme for dynamic adjustment of parameter ξ which outperforms the best parameterisation (ξ = 200) for variant (A). In this scheme, we start with ξ = 25. After each column generation convergence, we calculate the average ratio (among all calls to exact pricing) of the number of dominance checks between labels in a same bucket and the total number of non-dominated labels. If this average ratio is larger than 500 and the average number of non-fixed bucket arcs (including jump bucket arcs) per vertex is less than 10'000, we multiply ξ by 2. After each such multiplication, the bucket graph is regenerated.

In Figure 6(a), we show the performance profile which compares the total time for variant (A) with ξ = 200 and for variant (A + ). The geometric mean of the speed-up ratios between these two variants is 7%, and the maximum speed-up reaches 83%.

In Figure 6(b) we show the performance profile which compares variants (A + ) and (A + B) of the algorithm. Although the geometric mean of the speed-up ratios due to improved bucket arc elimination procedure is only 6%, for two instances the speed-up ratio reaches the factor of two. These are instances RC2 2 8 and R2 2 4. In general, we noticed that one significantly benefits from the improved bucket arc elimination procedure on hard instances with small primal-dual gap. Such instances have usually long routes.

(A + B) (A + ) (b)
Total time for a variant is at most X times larger than the best Figure 6: Performance profiles for variants (A), (A + ), and (A + B) of the algorithm Concluding this part, we summarize the results in Table 3. In it, we show the geometric means of total solution time ratios for different variants of the algorithm. There, variant (A) uses parameterisation ξ = 200. We give these ratios for all 74 instances as well as for subclasses of instances. Looking at Table 3, we can say that the most computational improvement is achieved when passing from variant (O) to variant (A) of the algorithm provided that the bucket step size is well chosen. Additional speed-ups of variants (A + ) and (A + B) are relatively small on average but can be significant for certain instances.

Instances (O)→(A) (A)→(A + ) (A + )→(A + B)

Impact of two-dimensional buckets

As proposed in Section 3.3, labels are grouped into buckets based on the accumulated resource consumption values. Thus, if there are two resources, we can say that buckets are "twodimensional". However, buckets can also be formed based only on the accumulated consumption of only one of two resources. Implementation of the labeling algorithm and bucket arc elimination procedure with such "one-dimensional" buckets is easier as sets Φ b , b ∈ B, contain at most one element. In this section we computationally verify whether an additional effort to implement two-dimensional buckets pays off, i.e. it results in improved algorithm efficiency.

In this experiment, we use variant (A + B) of the algorithm to solve the root node of instances with two resources. As the number of VRPTW instances with two resources is only 30 (these are instances of classes C1, RC1, and R1), we add in this experiment MDVRP, SDVRP, and DCVRP instances. In total, we have 69 instances. As before, for instances which cannot be solved at the root, we set the "target" lower bounds, given in Table 4. Again, for instances not shown in the table, the target lower bound equals the optimum solution value. As before, when using two-dimensional buckets, given parameter ξ, bucket step size d1 is determined as W/ √ ξ and d2 is calculated as (u n+1 -l 0 )/ √ ξ. Naturally, when using one-dimensional buckets, "bucketization" is performed according to the most critical resource. This is the capacity resource for VRPTW instances ( d1 = W/ξ) and the time resource for the remaining ones ( d1 = (u n+1 -l 0 )/ξ). The most critical resource can easily be determined by running a reasonable heuristic for the problem and analyzing the solutions obtained. As in the previous section, two copies of the labeling algorithm are created to obtain the same sequence of primal-dual solutions of the master problem for both tested configurations.

In Figure 7 we show the performance profile which compares variants with one-dimensional and two-dimensional buckets. The geometric means of the speed-up ratios between the two variants is 21%. On individual instances the maximum speed-up factor is 6.2. From detailed results (not given here) it can be observed that it is important to use two-dimensional buckets when both resources are similarly constraining.

one-dimensional buckets two-dimensional buckets

Total time for a variant is at most X times larger than the best Figure 7: Performance profiles for variants with one-dimensional and two-dimensional buckets

In Table 5 we give the geometric means of total solution time ratios for five variants of the algorithm applied for instances with two resources. First three variants (O), (A), (A + ), (A + B) use one-dimensional buckets, and the last variant (A + B) uses two-dimensional buckets. Again, variant (A) uses parameterisation ξ = 200. We time ratios for all 69 instances as well as for instances of different problems separately. Results in Table 5 suggest the most computational improvement is again achieved when passing from variant (O) to variant (A). Significant speed-ups achieved when passing from (A) to (A + ) and from one-dimensional buckets to two-dimensional. Impact of the bucket arc elimination procedure is very small here.

Results for the BCP algorithm

In this section, we present the performance of the proposed BCP algorithm on instances from the literature. We start with 74 VRPTW instances having 100 and 200 customers. We compare our algorithm with the state-of-the-art BCP algorithm for this problem by Pecin et al. (2017a). The One-dimensional buckets The computer used in Pecin et al. (2017a) and our computer have approximately the same speed, according to www.cpubenchmark.net/singleThread.html. In Table 6, we present the comparison of the two algorithms. Our algorithm outperforms the one by Pecin at al. for all classes of VRPTW instances, except class R1 for which the results are similar. Our algorithm proved optimality for five out of nine open instances with 200 customers for the first time. On the set of instances solved by the both algorithms, our approach is three times faster if we compare average solution times, and four times faster if we compare the geometric mean of solution times. Speed-up factor on individual instances reaches 47.5 as demonstrated in the detailed results presented in the online appendix. Interestingly, the largest improvement factors are reached on instances for which the improved arc elimination procedure had the most impact.

One → two Instances Number (O)→(A) (A)→(A + ) (A + )→(A + B) (A + B
We have also tested our algorithm on VRPTW instances with 400 customers. In the literature only one of these instances has been solved to optimality, by [START_REF] Kallehauge | Lagrangian duality applied to the vehicle routing problem with time windows[END_REF]. We managed to solve to optimality 24 instances out of 60. We observed three main reasons (in decreasing order of the number of impacted instances) why algorithm could not solve more 400-customer instances: (i) the labeling algorithm is not fast enough, so only relatively few R1Cs can be added before the pricing time limit is exceeded, leaving a gap that is too large (ii) R1Cs with up to five rows can be added almost until convergence, but they still leave a gap that is too large, and (iii) column generation convergence is very slow.

In the next experiment, we compare our algorithm with the exact algorithm by [START_REF] Contardo | A new exact algorithm for the multi-depot vehicle routing problem under capacity and route length constraints[END_REF] on the set of distance-constrained MDVRP instances. The main features of that previous algorithm that differ from ours are: • Absence of column generation stabilization;

• Separation of different families of valid inequalities, particularly, only a subset of R1Cs is used (Subset Row Cuts), and the limited memory technique is not applied;

• Exclusive use of route enumeration for finishing the instance, without any branching.

The speed of the computer used in [START_REF] Contardo | A new exact algorithm for the multi-depot vehicle routing problem under capacity and route length constraints[END_REF] is approximately 1.6 times slower than ours, according to www.cpubenchmark.net/singleThread.html. Therefore, in the following comparison, as well in the detailed results in the online appendix, the solutions times of [START_REF] Contardo | A new exact algorithm for the multi-depot vehicle routing problem under capacity and route length constraints[END_REF] are divided by 1.6. In Table 7, we present the comparison of the two algorithms. Our algorithm solves all instances, most of them for the first time. When compared on the set of instances solved by the both algorithms, our approach is 110 times faster in terms of the average time, and 16.5 times faster in terms of the geometric mean time (33 times faster if the smallest instance "pr01" is not taken into account). Speed-up factor on individual instances reaches three orders of magnitude as demonstrated in the detailed results presented in the online appendix. We believe that all the differences pointed above, not only the bucket graph pricing, contribute to this big improvement. For the instance "pr10", we improved the best known solution. For other instances, the values obtained by [START_REF] Vidal | A hybrid genetic algorithm for multidepot and periodic vehicle routing problems[END_REF] are optimal.

In the following Table 8 we give the summary results for remaining instances. To our knowledge, these instances are treated by an exact approach for the first time. We solved to optimality six out of seven classic distance-constrained CVRP instances by [START_REF] Christofides | The vehicle routing problem[END_REF]. All best known solution values for these instances given in [START_REF] Nagata | Edge assembly-based memetic algorithm for the capacitated vehicle routing problem[END_REF] are optimal. The only unsolved instance CMT13 is clustered and has many long arcs in the solution. The gap for this instance is large even after adding many R1Cs. This gap can be reduced by adding even more cuts, but exact pricing becomes too expensive.

We solved to optimality seven out of ten of distance-constrained site-dependent instances. These instances happen to be harder than similar multi-depot instances in class "pr". Best known solution values from [START_REF] Cordeau | A parallel iterated tabu search heuristic for vehicle routing problems[END_REF] were improved for four instances.

Finally, we have considered HFVRP instances from [START_REF] Duhamel | Efficient frameworks for greedy split and new depth first search split procedures for routing problems[END_REF] who called these instances "nightmare". Particularities of these instances are large capacities of vehicles, very heterogeneous customer demands, and many -up to eight -very heterogeneous vehicle types. We confirm that these instances are particularly difficult. We could not solve several instances which have slightly more than 100 customers. Nevertheless, we have solved to optimality all instances with up to 100 customers, 29 from 38 instances with 101-150 customers, 13 from 30 instances with 151-200 customers. Improved solutions were found for 43 instances.

The detailed results for all tested instances are given in the online appendix. The new algorithms were embedded in an advanced BCP algorithm, that includes several other features picked from the literature. The computational results were highly positive. On the classical VRPTW, the new BCP is at least four times faster on average than a very recent BCP algorithm (Pecin et al. 2017a). Moreover, the new BCP solved five additional VRPTW instances to optimality. Since most other elements in those two BCP algorithms, like ng-paths and rank-1 cuts, are similar, the significant improvement is clearly due to the techniques introduced in the paper. The experimental investigation suggests that the largest improvement comes from the organisation of labels in buckets which decreases the number of dominance checks in the pricing. In fact, this proposed labeling implementation can better handle the large numbers that define the time resource in the benchmark VRPTW instances. Other components as dynamic adjustment of the bucket step size and an improved arc elimination procedure bring less speed-up on average but can be important for some instances.

Very good results were also obtained on instances from other important VRP variants, including classical DCVRP instances proposed almost 40 years ago, but only now solved to optimality.

Although the bucket implementation presented in this paper considered up to two resources, it is not difficult to extend it for handling the SPPRCs with several resources that arise in other VRP variants. We believe that those harder problems may also benefit from the approach. However, substantial additional experimental work would be needed for determining how many and which resources should be used for defining the buckets, the remaining resources being treated as in the standard labeling. helped to improve the quality of the paper.

(ratio of the running times of the literature algorithm and ours). In Tables 9101112131415, "Time PCDU" refers to the running time of the algorithm by Pecin et al. (2017a). In Tables 2223, "Time CM" refers to the running time of the algorithm by [START_REF] Contardo | A new exact algorithm for the multi-depot vehicle routing problem under capacity and route length constraints[END_REF]. In the tables, * near the best solution value indicates that its optimality was proven for the first time. If the best solution value is underlined, it improves on the best known solution value. All obtained solutions were proven to be optimal except when the total time is given with the ">" sign.

Figure 1 :

 1 Figure 1: An example of a forward and the corresponding backward paths in the case of one resource

  For a forward (backward) bucket b, let κ b be the vector κ ∈ K ṽb such that lb = l ṽb + κ d (ũ b = u ṽb -κ d). We say that bucket b is component-wise smaller than bucket b (denoted as b ≺ b) if buckets are of the same sense, ṽb = ṽb , and κ b ≤ κ b .

  where v is the head of a γ . The head of a backward bucket arc γ, denoted by b head γ , is the bucket b ∈ B such that ṽ b = v and ũ b ≥ min ũ b γ -d a γ , u v > ũ b -d. Let B = ( B, Γ) and B = ( B, Γ) be the directed forward and backward bucket graphs, respectively. For each forward/backward bucket b ∈ B, we define the set Φ b of "adjacent" component-wise smaller buckets: Φ b = b : ṽb = ṽb , ∃r ∈ R, κ b r = κ b r -1, κ b r = κ b r , ∀r ∈ R \ {r } . As |R| = 2, we have |Φ b | ≤ 2. By definition, we have b ≺ b for all b ∈ Φ b . Let Γ Φ ( Γ Φ ) be the set of directed forward (backward) bucket arcs (b , b) such that b ∈ B, b ∈ Φ b . Let B Φ = ( B, φ( Γ) ∪ Γ Φ ) be the extended forward bucket graph, where φ( Γ) = {(b γ , b head γ ) : γ ∈ Γ}.

  v

Figure 3 :

 3 Figure 3: An example of Extended Bucket Graph for the forward graph of Figure 2

b

  equal to the smallest reduced cost of labels L such that b L b. It also uses the auxiliary function DominatedInCompWiseSmallerBuckets(L, b, B visited ) presented below, which checks whether a label L is dominated by a label contained in a bucket b such that b b, b ∈ B visited . To avoid checking all component-wise smaller buckets, it assumes that the values of cbest b are updated (for all buckets preceding b L in the topological order) and uses it to prune the search. If the search reaches a bucket b such that cbest b

  Two buckets b and b such that b ∈ Φ b may belong to the same strongly connected component. This may happen, for example, when b and b are forward buckets and there is another forward bucket b with ṽb = ṽb = ṽb such that lb ,r < lb ,r < lb,r , for some r ∈ R. In that case, the extended bucket graph may contain the directed cycle ((b , b), (b, b ), (b , b )). Therefore values cbest b for buckets in the same strongly connected component are calculated in a lexicographic order of κ b to ensure that cbest b is set before cbest b .

  by a label in b L then mark L as non-extended and insert in b L remove labels dominated by L from b L mark L as extended until all labels in all buckets b ∈ C are extended foreach bucket b ∈ C in a lexicographic order of κ b do cbest b ← min min L: b L =b {c L }, min b ∈Φ b {c best b } return P best = argmin P L :v L =v cL , where v is the sink in G Function DominatedInCompWiseSmallerBuckets(L, b, B visited ) B visited ← B visited ∪ {b} if C b precedes C b L in the topological order used and cL < cbest b then return false if b = b L and L is dominated by a label in bucket b then return true for b ∈ Φ b \ B visited do if DominatedInCompWiseSmallerBuckets(L, b , B visited ) then return true return false

Figure 4 :

 4 Figure 4: Illustration of the proof of Proposition 1

  Each forward/backward jump bucket arc ψ ∈ Ψ is characterized by a triple (b base ψ , b jump ψ , a ψ ), where b base ψ is its base bucket, b jump ψ is its jump bucket such that ṽb base ψ = ṽb jump ψ = ṽ, b base ψ ≺ b jump ψ

  {b b : ∃γ ∈ Γ : bγ = b , aγ = a} remove from B all buckets which are not component-wise minimal foreach b ∈ B do Ψ ← Ψ ∪ {ψ}, where b base ψ = b, b jump ψ = b , a ψ = a return Ψ

  P i , b base ψ = b L , and b jump ψ b L P i

.

  For a (jump) bucket arc γ ∈ Γ (ψ ∈ Ψ), we denote as b b arr γ ( b b arr ψ ) its arrival bucket of the opposite sense. For a forward bucket arc γ ∈ Γ, its arrival bucket b arr γ is the backward bucket b such that ũ b ≥ max l b γ + d a γ , l v > ũ b -d, where v is the head node of a γ . For a backward bucket arc γ ∈ Γ, its arrival bucket b arr γ is the forward bucket b such that l b ≤ min ũ b γ -d a γ , u v < l b + d, where v is the head node of a γ . For a jump bucket arc ψ ∈ Ψ, the definition of its forward/backward arrival buckets is identical, just replacing b γ by b jump ψ Proposition 4. Let set Γ be a superset of all improving bucket arcs, Ψ be the set of jump bucket arcs computed using procedure ObtainJumpBucketArcs(Γ), and value θ be calculated according to (16). If γ ∈ Γ is an improving bucket arc then there exists a pair of θ-compatible labels L ∈ L and L L ∈ L L such that b L b γ and b b L L b b arr γ , where L = L(Γ, Ψ), and L L = L L(Γ, Ψ).

  P j and L L P P |P |-j-1 are θ-compatible. If label L P j ∈ L, let L = L P j , otherwise by Proposition 1 there exists a label L ∈ L dominating L P j such that b L b L P j = b γ . If label L L P P |P |-j-1 ∈ L L, let L L = L L P P |P |-j-1 , otherwise by Proposition 1 there exists a label L L ∈ L L dominating L L P P |P |-j-1 such that b b L L b b L L P P |P |-j-1 b b arr γ . Finally, by Proposition 2, L and L L are θ-compatible.

  sets B B a,b of buckets of opposite sense with respect to bucket b at the head of arc a. B B a,b contains every bucket b b containing a label L L ∈ L L which is θ-compatible with a label L ∈ L such that b L b. Thus, a bucket arc γ is non-improving if B B aγ ,bγ does not contain any bucket b b b b arr γ . An advantage of using sets B B a,b is that they can be constructed incrementally using adjacent buckets: every set B B a,b is initialized as the union of sets B B a,b , b ∈ Φ b , and then augmented using the recursive procedure UpdateBucketsSet which tries to find opposite sense buckets containing a label which is θ-compatible with a label in b. We now present formally the auxiliary procedure UpdateBucketsSet(θ, L, B B, b b, L L, B B visited ), which adds to set B B all opposite sense buckets b b , b b b b, b b ∈ B B visited , containing at least one label L L ∈ L L such that pair (L, L L) of labels is θ-compatible. Procedure UpdateBucketsSet(θ, L, B B, b b, L L, B B visited )

  ,b : b b b b arr γ then Γ ← Γ \ {γ} (backward) bucket arcs are present in the current forward (backward) bucket graph before the (i + 1)th call to procedure BucketArcElimination. Then, by Proposition 4, for each improving (forward or backward) bucket arc γ ∈ Γ, the bi-directional labeling algorithm generates a pair of θ-compatible labels L and L L such that b L b γ and b b L L b b arr γ .

  100, β 3 = 150. The tailing off parameter values are δ perc = 1.5%, δ num = 3. Dynamic ng-relaxation parameter values are η init = 8, η 1 = 5, η 2 = 100, η max = 16. Parameter value for the bucket graph elimination procedure is υ = 10%. Enumeration parameter values are ω labels = 1 000 000, ω routes = 5 000 000, ω MIP = 10 000. The strong branching parameter values are ζ 0 = 50, ζ 1 = 3.

Figure 5 :

 5 Figure 5: Impact of parameter ξ on the subproblem and total time

Table 1 :

 1 Parameterisation of pricing time thresholds HFVRP instances For the heterogeneous fleet vehicle routing problem, we used 96 instances by Duhamel, Lacomme, and Prodhon (2011) with 20-256 customers. Each instance represents one of 96 French departments. The best known solutions are taken from Penna

	Problem	Size	Resources	τ	τ	r *
	VRPTW	100-200 1	10 30 time
	VRPTW	200	2	5 10 capacity
	VRPTW	400	1	20 60 time
	VRPTW	400	2	10 20 capacity
	HFVRP	any	1	5 10 capacity
	DCVRP	any	2	5 10 time
	MDVRP, SDVRP any	2	10 20 time

Table 2 :

 2 Target lower bounds for VRPTW instances we determine that a bucket arc γ ∈ Γ cannot be fixed, we do not try to fix bucket arcs γ ∈ Γ such that a γ = a γ . The resource value dependent arc elimination procedure proposed inPessoa et al. (

	Instance	BKS	TLB	Instance	BKS	TLB	Instance	BKS	TLB
	R208	701.0	699.0	R2 2 4	1928.5 1925.0	RC1 2 7	3177.8 3150.0
	C2 2 3	1763.4 1753.1	R2 2 6	2675.4 2674.6	RC1 2 8	3060.0 3044.8
	C2 2 4	1695.0 1663.5	R2 2 8	1842.4 1831.6	RC1 2 9	3074.8 3035.7
	R1 2 3	3373.9 3352.2	R2 210	2549.4 2545.3	RC1 210 2990.5 2966.0
	R1 2 4	3047.6 3035.3	RC1 2 1	3516.9 3499.1	RC2 2 2	2481.6 2477.5
	R1 2 5	4053.2 4044.1	RC1 2 2	3221.6 3204.4	RC2 2 4	1854.8 1848.0
	R1 2 6	3559.2 3544.2	RC1 2 3	3001.4 2980.7	RC2 2 7	2287.7 2279.7
	R1 2 8	2938.4 2933.3	RC1 2 4	2845.2 2832.8	RC2 2 8	2151.2 2150.3
	R1 2 9	3734.7 3726.0	RC1 2 5	3325.6 3310.0	RC2 2 9	2086.6 2072.0
	R1 210	3293.1 3277.3	RC1 2 6	3300.7 3289.0	RC2 210 1989.2 1962.0

Table 3 :

 3 Comparison of different variants of the algorithm on VRPTW instances

	100 customers	1.54	1.18	1.06
	200 customers	2.01	1.04	1.07
	Series 1	2.10	1.11	1.05
	Series 2	1.78	1.03	1.07
	Series C	1.20	1.21	1.03
	Series R	1.62	0.99	1.08
	Series RC	4.15	0.97	1.09
	All	1.91	1.07	1.06

Table 4 :

 4 Target lower bounds for MDVRP, SDVRP, and DCVRP instances

	Instance	BKS	TLB	Instance	BKS	TLB	Instance	BKS	TLB
	CMT8	865.94	857.8	mdvrp p10	3631.11 3626.2	sdvrp pr04 3449.84 3432.3
	CMT9	1162.55 1153.3	mdvrp p11	3546.06 3540.4	sdvrp pr05 4377.35 4332.6
	CMT10	1395.85 1387.3	mdvrp pr05 2331.20 2318.7	sdvrp pr06 4422.02 4329.2
	CMT13	1541.30 1475.2	mdvrp pr06 2676.30 2662.9	sdvrp pr08 2971.01 2962.6
	mdvrp p08 4372.78 4347.6	mdvrp pr10 2868.26 2860.6	sdvrp pr09 3536.20 3527.7
	mdvrp p09 3858.66 3849.9	sdvrp pr03	2575.36 2565.6	sdvrp pr10 4639.62 4659.4

Table 5 :

 5 Comparison of different variants of the algorithm on instances with two resources

	)

Table 6 :

 6 Comparison of our BCP algorithm with the one inPecin et al. (2017a) on the VRPTW instances difference between the two algorithms mainly resides in the different labeling routines included. Another difference is a better dual price smoothing stabilization with the automatically adjusted parameter[START_REF] Pessoa | Automation and combination of linear-programming based stabilization techniques in column generation[END_REF]) used by us. Also, we use only standard RCCs and limited-memory R1Cs with up to five rows, whereas algorithm inPecin et al. (2017a) uses also elementary inequalities (a family of R1Cs) with six and seven rows and clique cuts (after enumeration). Finally, different mechanisms to dynamically adjust ng-neighborhoods are employed. Other components used in both algorithms are similar.

Table 7 :

 7 Comparison of our BCP algorithm with the one in[START_REF] Contardo | A new exact algorithm for the multi-depot vehicle routing problem under capacity and route length constraints[END_REF] on the distance-constrained MDVRP instances• Use of standard mono-directional label correcting algorithm for solving the pricing problem (improved by decremental state-space relaxation);

			Our algorithm		Algorithm in Contardo and Martinelli (2014)
	Class	Solved	Average Geomean time time	Solved	Average time	Geomean time
	"p"	12/12	4538	201	6/6	938	342
	"pr"	10/10	13773	937	4/7	39058	903
	Overall	22/22	7881	392	10/13	16186	504
	Solved by both		147	30		16186	504

Table 8 :

 8 Results for instances never attacked by exact methods before7 ConclusionsThis article proposed a new bucket-graph based variant of the labeling algorithm for the SPPRC. It can handle the particularities of the SPPRCs that arise as pricing suproblems in the most effective existing exact algorithms for the VRP: the partial route elementarity defined by ng-sets and the numerous additional dimensions induced by limited-memory R1Cs. The algorithm is especially better than existing variants when the consumptions, in at least one of the resources, are given either by real numbers or by large integer numbers. Besides, a new procedure that generalizes the principle of arc elimination by reduced costs to the bucket graph context was given.

Table 15 :

 15 Detailed results for 200-customer Gehring and Homberger VRPTW instances of class RC2

	Instance	Initial Robust U B bound	bound	Root time	Nodes #R1C number	Total time	Final U B	Time PCDU	Improv. factor
	C203	588.8	588.7	588.7	27s	0	1	27s	588.7		6m05s	13.3
	C204	588.2	588.1	588.1	40s	0	1	40s	588.1 27m28s	40.9
	R202	1029.7	1022.2 1029.6	2m25s	72	1 2m25s 1029.6		5m16s	2.2
	R203	870.9	866.9	870.8	1m16s	17	1 1m16s	870.8		9m25s	7.4
	R204	731.4	724.9	731.3	3m04s	56	1 3m04s	731.3 24m55s	8.1
	R206	876.0	866.8	875.9	2m30s	51	1 2m30s	875.9		9m06s	3.6
	R207	794.1	790.1	794.0	1m43s	27	1 1m43s	794.0 12m08s	7.0
	R208	701.1	691.6	697.7 14m57s	138	9 1h14m	701.0 17h48m	14.3
	R209	854.9	841.4	854.8	3m45s	105	1 3m45s	854.8		9m59s	2.7
	R210	900.6	889.4	900.5	4m34s	107	1 4m34s	900.5 12m09s	2.7
	R211	746.8	734.6	746.7	3m42s	110	1 3m42s	746.7 22m00s	5.9
	RC204	783.6	779.3	783.5	2m11s	15	1 2m11s	783.5 19m27s	8.9
	RC207	963.0	947.3	962.9	3m21s	91	1 3m21s	962.9		5m47s	1.7
	RC208	776.2	766.6	776.1	2m02s	36	1 2m02s	776.1		9m59s	4.9
			Table 9: Detailed results for the hardest 100-customers Solomon
			VRPTW instances
	Instance	Initial Robust U B bound	bound	Root time	Nodes #R1C number	Total time	Final U B	Time PCDU	Improv. factor
	C1 2 1	2698.7	2698.6 2698.6	8s	0	1	8s 2698.6		20s	2.2
	C1 2 2	2694.4	2694.3 2694.3	20s	0	1	20s 2694.3		44s	2.1
	C1 2 3	2675.9	2659.7 2675.8 4m39s	66	1 4m39s 2675.8 10m26s	2.2
	C1 2 4	2625.7	2619.8 2625.6 6m51s	42	1 6m51s 2625.6	5m40s	0.8
	C1 2 5	2695.0	2694.9 2694.9	11s	0	1	11s 2694.9		24s	2.1
	C1 2 6	2695.0	2694.9 2694.9	15s	0	1	15s 2694.9		26s	1.7
	C1 2 7	2695.0	2694.9 2694.9	18s	0	1	18s 2694.9		22s	1.2
	C1 2 8	2684.1	2683.7 2684.0	19s	0	1	19s 2684.0		28s	1.4
	C1 2 9	2639.7	2637.2 2639.6	54s	0	1	54s 2639.6	1m24s	1.5
	C1 2 10	2624.8	2620.2 2624.7 2m40s	34	1 2m40s 2624.7	1m53s	0.7
		Table 10: Detailed results for 200-customer Gehring and Homberger
		VRPTW instances of class C1
	Instance	Initial Robust U B bound	bound	Root time	Nodes #R1C number	Total time	Final U B		Time PCDU	Improv. factor
	C2 2 1	1922.2	1915.0 1922.1	3m06s	70	1	3m06s	1922.1 14m20s	4.6
	C2 2 2	1851.5	1839.4 1851.4	6m57s	80	1	6m57s	1851.4	1h29m	12.9
	C2 2 3	1763.5	1737.0 1752.1 32m16s	287	37 11h10m 1763.4 *	>96h	>8.6
	C2 2 4	1695.1	1652.3 1663.2 28m59s	187	355	>60h	1695.1	>96h	-
	C2 2 5	1869.7	1858.7 1869.6	3m57s	52	1	3m57s	1869.6 23m48s	6.0
	C2 2 6	1844.9	1836.6 1844.8	3m22s	23	1	3m22s	1844.8 34m03s	10.1
	C2 2 7	1842.3	1834.9 1842.2	3m03s	21	1	3m03s	1842.2 44m33s	14.6
	C2 2 8	1813.8	1801.3 1813.7	3m11s	12	1	3m11s	1813.7 42m22s	13.3
	C2 2 9	1815.1	1795.0 1815.0	9m17s	101	1	9m17s	1815.0	2h01m	13.1
	C2 2 10	1791.3	1765.7 1791.2 13m04s	96	1 13m04s	1791.2	2h56m	13.5
			Table 11: Detailed results for 200-customer Gehring and Homberger
			VRPTW instances of class C2

Table 20 :

 20 Detailed results for 400-customer Gehring and Homberger VRPTW instances of class RC1

	Instance	Initial Robust U B bound	bound	Root time	Nodes #R1C number	Total time	Final U B
	RC2 4 1	6147.4	6120.1 6146.6	3h18m	238	3	3h37m 6147.3 *
	RC2 4 2	5407.6	5389.5 5407.5	4h06m	458	1	4h06m 5407.5 *
	RC2 4 3	4573.1	4549.6 4573.0 35h25m	1294	1 35h25m 4573.0 *
	RC2 4 4	3598.0	3535.9 3535.9	8h25m	0	5	>60h	3598.0
	RC2 4 5	5396.7	5363.0 5388.6	4h47m	466	7	>60h	5396.7
	RC2 4 6	5332.1	5270.9 5308.5	9h04m	738	3	>60h	5332.1
	RC2 4 7	4987.9	4944.9 4968.6	6h08m	453	19	>60h	4987.9
	RC2 4 8	4703.2	4602.0 4618.7	4h53m	216	31	>60h	4703.2
	RC2 4 9	4519.1	4448.9 4457.0	5h35m	278	21	>60h	4519.1
	RC2 4 10 4252.4	4157.0 4157.0	2h11m	0	83	>60h	4252.4

Table 21 :

 21 Detailed results for 400-customer Gehring and Homberger VRPTW instances of class RC2

	Inst.	n	m	Initial U B	Robust bound	bound	Root time	Nod. #R1C num.	Total time	Final U B	Time CM	Impr. factor
	p08	249	2 4372.90 4294.01 4347.67 49m01s	567	29	5h43m 4372.78 * ,a	-	-
	p09	249	3 3858.80 3798.47 3849.96	1h34m	650	19	2h25m	3858.66 *	-	-
	p10	249	4 3631.30 3582.48 3626.26 37m26s	522	3 41m46s 3631.11 * ,a	-	-
	p11	249	5 3546.20 3488.64 3540.43 44m15s	597	3 49m46s 3546.06 * ,a	-	-
	p13	80	2 1319.10 1318.95 1318.95	1s	0	1	1s	1318.95	31s	17.5
	p14	80	2 1360.30 1360.12 1360.12	1s	0	1	1s	1360.12	38s	20.5
	p16	160	4 2572.40 2572.23 2572.23	12s	0	1	12s	2572.23	8m27s	39.7

Table 26 :

 26 Detailed results for Duhamel et al. HFVRP instances with up to 100 customers

	Instance	n	m	Initial U B	Robust bound	bound	Root time	Nodes #R1C number	Total time	Final U B
	DLP 90	102	4	2347	2249.36	2257.75	5m05s	136	655 40h29m	2278.05 *
	DLP 17	105	3	5363	5278.74	5326.71	2m50s	131	7	3m30s	5362.83 *
	DLP 84	105	4	7228	7117.02	7142.28	2m02s	136	19 21m01s	7227.88 *
	DLP 81	106	4 10687 10483.95 10521.23	3m06s	182	313	5h36m 10583.50 *
	DLP 2B 107	6	8463	8300.92	8332.22	4m17s	37	1639	>60h	8463.00
	DLP 07	108	4	8090	8009.35	8049.38	3m51s	104	19 12m07s	8071.97 *
	DLP 87	108	4	3754	3658.33	3692.55	6m32s	236	7 29m51s	3753.87 *
	DLP 47	111	5 16157 15971.69 16123.03	3m54s	250	5	5m49s 16156.12 *
	DLP 48	111	5 21310 21041.94 21228.34 19m51s	321	7 28m25s 21257.38 *
	DLP 61	111	3	7293	6997.61	7027.77	3m15s	33	1435	>60h	7293.00
	DLP 12	112	4	3544	3410.14	3410.14	1m53s	0	25	1h39m	3543.99 *
	DLP 30	112	3	6314	6185.97	6207.60	4m11s	186	1489 31h29m	6278.62 *
	DLP 2A 113	6	7794	7708.25	7793.16	1m42s	77	1	1m42s	7793.16 *
	DLP 53	115	3	6435	6281.51	6376.47 12m26s	209	31	2h14m	6434.83 *
	DLP 05	116	5 10877 10764.90 10869.04	29s	30	1	29s 10869.04 *
	DLP 13	119	5	6697	6618.19	6649.75	5m20s	160	1261	>60h	6696.43
	DLP 06	121	8 11689 11618.62 11682.98	2m37s	53	1	2m37s 11682.98 *
	DLP 03	124	4 10710 10635.69 10669.77	55s	106	5	1m13s 10709.66 *
	DLP 83	124	4 10020	9871.11	9899.78	2m33s	127	59 51m57s 10001.80 *
	DLP 68	125	4	8971	8799.96	8852.89	5m39s	192	89	3h00m	8889.03 *
	DLP 74	125	5 11587 11543.40 11563.64	41s	78	7	4m18s 11586.58 *
	DLP 21	126	3	5140	5079.13	5110.52	8m11s	289	11 25m05s	5139.84 *
	DLP 26	126	5	6434	6352.66	6382.52	6m00s	230	453 10h10m	6406.16 *
	DLP 88	127	5 12389 12328.90 12385.74	29s	25	1	29s 12385.74 *
	DLP 16	129	6	4157	4000.74	4006.42	4m33s	16	1573	>60h	4157.00
	DLP 51	129	3	7722	7628.81	7666.22 10m11s	123	3 29m17s	7721.47 *
	DLP 31	131	8	4092	4045.03	4083.80 15m20s	284	5 18m24s	4091.52 *
	DLP 40	132	5 11119 10974.62 11015.52	2m32s	160	317	2h59m 11056.13 *
	DLP 89	134	5	7087	6982.53	7023.43	4m11s	239	1027	>60h	7087.00
	DLP 41	135	7	7598	7444.36	7444.36	4m48s	0	3305	>60h	7598.00
	DLP 34	136	6	5748	5700.07	5727.33	7m35s	306	27	1h10m	5739.02 *
	DLP 60	137	4 17037 16768.33 16934.97	5m40s	285	37 44m38s 17012.42 *
	DLP 73	137	5 10196 10138.46 10174.15	57s	88	3	2m07s 10195.33 *
	DLP 28	141	5	5531	5478.74	5499.59	4m24s	208	4705	>60h	5525.90
	DLP 25	143	6	7207	7030.07	7042.25	5m46s	100	959	>60h	7207.00
	DLP 85	146	4	8774	8690.09	8726.51	7m27s	295	111	2h40m	8763.90 *

Table 27 :

 27 Detailed results for Duhamel et al. HFVRP instances with 101-150 customers

	Instance	n	m	Initial U B	Robust bound	bound	Root time	Nodes #R1C number	Total time	Final U B
	DLP 69	152	4	9148	9003.97	9086.90	7m52s	335	147	3h36m	9127.16 *
	DLP 76	152	8 11995 11900.23 11933.72	5m28s	212	1195	>60h	11994.22
	DLP 56	153	4 31031 30571.34 30764.30	8m26s	250	1591 58h50m 30905.95 *
	DLP 86	153	5	9028	8935.26	8957.90	8m23s	274	249	8h45m	9020.63 *
	DLP 37	161	5	6851	6751.15	6787.92	7m57s	354	2017	>60h	6838.72
	DLP 64	161	3 17136 16552.72 16552.72	6m38s	25	995	>60h	17135.16
	DLP 24	163	4	9102	8990.07	8990.07	3m26s	0	1907	>60h	9102.00
	DLP 57	163	4 44782 44302.32 44409.19	8m53s	194	1565	>60h	44782.00
	DLP 29	164	4	9143	9044.54	9078.84	9m59s	182	43	2h18m	9132.03 *
	DLP 09	167	5	7604	7543.04	7563.84 12m52s	274	135	2h16m	7599.72 *
	DLP 35	168	6	9556	9492.98	9521.16	2m55s	114	3	3m17s	9522.45 *
	DLP 45	170	3 10477 10294.91 10343.43	8m34s	124	929	>60h	10477.00
	DLP 80	171	3	6817	6712.75	6716.88	8m39s	29	219 13h07m	6816.89 *
	DLP 44	172	3 12192 11895.60 11895.60 12m31s	19	1187	>60h	12192.00
	DLP 54	172	4 10352 10089.37 10089.37	7m45s	21	1323	>60h	10352.00
	DLP 67	172	5 10885 10604.14 10662.94	9m53s	233	325 23h44m 10772.81 *
	DLP 63	174	5 19952 19670.12 19799.92	8m54s	108	97	5h24m 19890.65 *
	DLP 14	176	4	5645	5511.85	5511.85	4m37s	0	1023	>60h	5639.98
	DLP 42	178	7 10818 10615.08 10615.08	7m00s	0	889	>60h	10805.94
	DLP 02	181	4 11676 11603.82 11632.01	4m27s	104	25	9m34s 11649.81 *
	DLP 04	183	4 10749 10582.45 10673.39	9m41s	277	985 13h53m 10714.84 *
	DLP 95	183	2	6176	6100.98	6119.06	6m53s	290	1709	>60h	6170.20
	DLP 71	186	3	9835	9711.77	9771.86 10m30s	292	87	4h08m	9798.06 *
	DLP 72	186	4	5904	5780.39	5831.66	9m35s	322	365 31h31m	5870.43 *
	DLP 50	187	6 12371 12035.29 12035.29	9m56s	0	1521	>60h	12371.00
	DLP 15	188	7	8221	8083.44	8083.44	3m27s	0	873	>60h	8221.00
	DLP 33	189	7	9411	9248.07	9283.05	9m51s	135	1129	>60h	9333.39
	DLP 77	190	3	6917	6821.36	6826.04	7m49s	10	1215	>60h	6917.00
	DLP 78	190	4	7036	6899.96	6899.96	6m40s	0	1099	>60h	7036.00
	DLP 59	193	6 14283 13786.96 13786.96 23m42s	0	737	>60h	14283.00
	DLP 91	196	4	6375	6228.90	6244.50 10m26s	157	1167	>60h	6348.36

Table 28 :

 28 Detailed results for Duhamel et al. HFVRP instances with 151-200 customers

	Instance	n	m	Initial U B	Robust bound	bound	Root time	Nodes #R1C number	Total time	Final U B
	DLP 23	203	4	7742	7590.28	7603.76	8m46s	70	1043 >60h	7742.00
	DLP 38	205	5 11195 11048.66 11092.13 14m28s	375	869 >60h 11195.00
	DLP 58	220	6 23371 22800.97 22800.97 36m42s	10	969 >60h 23371.00
	DLP 27	220	5	8423	8315.95	8343.41 10m15s	174	687 >60h	8423.00
	DLP 65	223	3 13044 12648.64 12648.64 41m17s	24	657 >60h 13044.00
	DLP 19	224	5 11687 11503.00 11526.84	8m25s	98	1053 >60h 11687.00
	DLP 62	225	5 23011 22594.36 22645.25 13m36s	47	1047 >60h 22987.90
	DLP 22	239	2 13069 12923.02 12956.02 13m12s	357	1193 >60h 13050.00
	DLP 32	244	8	9383	9205.76	9235.38 14m29s	108	695 >60h	9383.00
	DLP 49	246	8 16182 15681.80 15681.80 15m27s	0	735 >60h 16182.00
	DLP 46	250	5 24567 23856.71 23856.71 14m14s	21	913 >60h 24567.00
	DLP 18	256	5	9653	9491.55	9547.80 22m29s	519	1159 >60h	9649.05

Table 29 :

 29 Detailed results for Duhamel et al. HFVRP instances with more than 200 customers

Acknowledgements

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" (see www.plafrim.fr/en/home).

AP and EU visited

Bordeaux with funding from project FAPERJ E26/110.075/2014 (International Cooperation FAPERJ/INRIA).

The authors are thankful to Thibaut Vidal for running the hybrid genetic algorithm and communicating initial upper bounds for VRPTW instances with 200 and 400 customers.

The authors are grateful to two anonymous referees and the associate editor whose comments

A Detailed results

In this appendix, we present the detailed results for our BCP algorithm on all tested instances.

In the following tables, we show the instance name, n -the number of customers, m -the number of different vehicle types or the number of depots, when applicable, the initial upper bound, the "robust" dual bound (before adding non-robust cuts), the dual bound obtained in the root node, the root node time, the number of R1Cs at the end of the root node, the number of nodes, the total time, and the best solution value obtained by the algorithm. For instances solved before by another algorithm in the literature, we give the corresponding time (which takes into account the difference between the speed of the computers used), and the improvement factor