Nabil Absi
email: absi@emse.fr

Diego Cattaruzza
email: cattaruzza@centralelille.fr

Dominique Feillet
email: feillet@emse.fr

Maxime Ogier
email: ogier@centralelille.fr

Frédéric Semet
email: semet@centralelille.fr

A heuristic branch-cut

published or not. The documents may come

Introduction

In this paper, we propose a heuristic solution method for the Inventory Routing Problem introduced during the 2016 ROADEF/EURO challenge, which we coin as REC-IRP. Inventory routing has attracted researchers for many years due to both its practical and theoretical issues [START_REF] Coelho | Thirty years of inventory routing[END_REF]2013). Many different variants have been investigated in the literature. The REC-IRP proposed for the challenge is original and particularly complex for several reasons. Among them, we can cite the logistic ratio optimization objective, the hourly time-granularity for inventory constraints or the driver/trailer allocation management. Designing an exact solution approach is out of reach for large size instances as those proposed during the challenge. However, we decided to address the REC-IRP with a branch-cut-and-price framework: a cut-and-column generation procedure is developed, along with a heuristic pricing algorithm to generate new columns and a heuristic fixing procedure to generate integer solutions.

Several modeling opportunities were opened when deciding what should be a column in our approach. We decided that a column would be what we call a td-shift, namely, a trailer-driver-shift. A td-shift is an elementary piece of work for a driver (the driver is known), that starts with a specific trailer from the base (the trailer is known), visits a subset of customers and returns at the base, with possibly intermediate stops at one or several sources to refill the trailer. The td-shift is time-stamped, meaning that the timing of all operations is fixed. Conversely, the quantities delivered to the customers are unknown.

Let us define a trip as a sequence of customers starting from and ending at the base or a source. Because of the multi-trip dimension of shifts, one could have decomposed shifts in trips and defined columns as trips. These two options (shifts or trips) have been compared in [START_REF] Hernandez | Branch-and-price algorithms for the solution of the multi-trip vehicle routing problem with time windows[END_REF] for the Multi-trip Vehicle Routing Problem with Time Windows. The theoretical advantage of a column generation approach based on trips is to limit the combinatorial explosion at the pricing problem level (a trip contains fewer customers than a shift). The counterpart is a complex master problem where one has to select trips but also to combine them into shifts, and eventually a more complex implementation of the branchand-price framework. Experimental results in [START_REF] Hernandez | Branch-and-price algorithms for the solution of the multi-trip vehicle routing problem with time windows[END_REF] do not permit to conclude definitely on the relative efficiency of the two approaches. For this reason, and given the challenging problem in hand, we decided for the easiest one.

A key feature of the problem concerns the quantities delivered during the shifts. In this respect, the problem is very similar to the Split Delivery Vehicle Routing Problem (SDVRP), where partial deliveries to customers are feasible. An important decision here is to decide whether the quantities delivered should be fixed in the column definition (pricing problem level) or decided at the master problem level. The two alternatives have been considered in [START_REF] Feillet | Vehicle routing with time windows and split deliveries[END_REF] and [START_REF] Desaulniers | Branch-and-price-and-cut for the vehicle routing problem with time windows[END_REF], respectively. The clear conclusion in the context of the SDVRP with Time Windows is that deciding on the quantities at the pricing problem level is better. However, the interaction between quantities delivered in different routes and for different customers seems stronger in inventory routing, and it is not clear whether such a conclusion can be generalized to this context. In the inventory routing literature, both options have been analyzed. [START_REF] Desaulniers | A branch-price-and-cut algorithm for the inventory-routing problem[END_REF] extends the approach developed in [START_REF] Desaulniers | Branch-and-price-and-cut for the vehicle routing problem with time windows[END_REF]. [START_REF] Bard | A branch-and-price algorithm for an integrated production and inventory routing problem[END_REF] follows a similar approach while solving the pricing problem as an integer program. [START_REF] Le | A column generation-based heuristic algorithm for an inventory routing problem with perishable goods[END_REF] or [START_REF] André | An approach based on column generation for solving routing problems with fractional objective function[END_REF] decide of the quantities at the master problem level. We decided to follow this latest scheme. Note also that some authors address easier cases where the delivered quantities are limited to a few possible discrete values [START_REF] Christiansen | Robust inventory ship routing by column generation, Column generation[END_REF][START_REF] Christiansen | Robust inventory ship routing by column generation, Column generation[END_REF][START_REF] Michel | A column-generation based tactical planning method for inventory routing[END_REF][START_REF] Michel | A column-generation based tactical planning method for inventory routing[END_REF]. This could also have been an option for us in a heuristic approach.

A last important characteristic of the problem, that may a priori significantly influence the solution scheme is the logistic ratio objective. This objective is naturally written as a fractional function of the decision variables. We are only aware of a few implementa-tions of the branch-and-price methodology for mathematical formulations with fractional objectives. [START_REF] Garaix | Optimization of occupancy rate in dial-a-ride problems via linear fractional column generation[END_REF] implemented and compared two approaches, adapting to the column generation framework algorithms presented in [START_REF] Charnes | Programming with linear fractional functionals[END_REF] and [START_REF] Dinkelbach | On nonlinear fractional programming[END_REF]. Initially they developed their research works in the context of mathematical programming with fractional objective functions. In short, Charnes and Cooper's approach relies on changes of variables. Dinkelbach's algorithm expresses the objective as a weighted linear combination of the numerator and the denominator, and solves the resulting program repeatedly while updating the weights (see more details in Section 3.4). Computational experiments in [START_REF] Garaix | Optimization of occupancy rate in dial-a-ride problems via linear fractional column generation[END_REF] show a similar efficiency for both methods. The authors concluded that having a fractional objective is not particularly detrimental to the solution of the problem. In this paper, we decided to use Dinkelbach's approach.

Recently, in [START_REF] Archetti | Minimizing the logistic ratio in the inventory routing problem[END_REF] and [START_REF] André | An approach based on column generation for solving routing problems with fractional objective function[END_REF], the first column generation approaches for inventory routing problems with logistic ratio objectives were proposed. In both cases, the authors tackled a standard Inventory Routing Problem (apart from the objective) with the aim to develop an exact solution method. In particular, the additional difficulties described above are not present (even if [START_REF] André | An approach based on column generation for solving routing problems with fractional objective function[END_REF] progressively includes features of the REC-IRP). [START_REF] Archetti | Minimizing the logistic ratio in the inventory routing problem[END_REF] apply a Dinkelbach's scheme. At each step, the linear program is solved using the method developed in [START_REF] Desaulniers | A branch-price-and-cut algorithm for the inventory-routing problem[END_REF]. As mentioned above, delivered quantities are fixed in column definitions. Only small instances, including at most 15 customers and defined on a time horizon limited to at most 5 periods, are solved to optimality. Contrary to what could be observed in [START_REF] Garaix | Optimization of occupancy rate in dial-a-ride problems via linear fractional column generation[END_REF], the fractional objective highly complicates the problem. The reason identified by the authors is that the value of the objective function tends to be flat around the optimal solution. As a consequence, there exist many solutions with objective function values close to the optimal one. This creates difficulties in the solution algorithm as it is much harder to prune branch-and-bound nodes. Paradoxically, this might rather be helpful for a heuristic solution of the problem. [START_REF] André | An approach based on column generation for solving routing problems with fractional objective function[END_REF] use Charnes and Cooper's method. The column definition is very similar to ours. Their method allows finding good lower and upper bounds on the same instances as those used in [START_REF] Desaulniers | A branch-price-and-cut algorithm for the inventory-routing problem[END_REF] and on a new instance set including simplified Air Liquide instances of the same size.

The paper is organized as follows. Section 2 introduces the notation. In section 3, we present our mixed integer linear programming formulation. Section 4 explains why a simultaneous cut-and-column generation strategy is needed and how it is carried out. In Section 5, we propose what we call an implicit-relaxation based dynamic programming algorithm, that constructs a dual vector from an implicit linear relaxation of the master problem, and price columns according to this vector. Section 6 describes the increasing horizon heuristic, that we developed, based on the formulation we proposed. Computational results are presented in Section 7, before concluding with Section 8.

Problem description and notation

This paper being part of a special issue devoted to the REC-IRP, we omit the complete description of the problem. Readers are referred to the first paper of the special issue to find this description. However, the paper is organized so that it should be possible to read and understand it independently. In particular, the notation is entirely reintroduced, and all constraints are progressively introduced and explained.

Tables 1 and2 provide a first round of notation. Note that (contrary to the convention taken in the problem description) we introduce a customer in V 2 for each order placed by a call-in customer. In the following, we, however, call these orders customers. Especially, we call customer set the set V = V 1 ∪ V 2 . Note also that the time horizon is discretized in one-hour periods, that starts at period 0 and ends at period H: H = {0, . . . , H}. The term period in the paper will always refer to this time granularity of one hour. When data or variables are expressed in minutes, it will be emphasized.

H

set of one hour periods in the time horizon

V 1 set of VMI customers V 2 set of call-in customer orders V set of customers D set
i ∈ V 1 each time it is visited q min i accumulated minimal quantity delivered to customer i ∈ V 2 q max i accumulated maximal quantity delivered to customer i ∈ V 2 stock i 0 initial stock for customer i ∈ V 1 stock i max maximal stock for customer i ∈ V 1 q init t initial load in trailer t ∈ T d i h demand of customer i ∈ V 1 at period h ∈ H

Mathematical formulation

We propose to formulate the REC-IRP with an extended formulation based on the selection of td-shifts. We introduce R, the set of time-stamped shifts. By time-stamped shift, we mean that the physical itinerary followed during the shift and the time of all the events of the shift are known. Conversely, the quantities delivered to customers, the trailer used or the driver assigned to the trailer are not known. Note that the concept of tdshift, introduced in Section 1, will only be useful when describing the column generation procedure (see Section 4). However, it is noteworthy that a td-shift is obtained by assigning a trailer and a driver to a time-stamped shift. From now on, time-stamped shifts will be called shifts for short.

To limit the complexity of the model, we make some assumptions on set R:

1. each customer is visited at most once in a shift; 2. a shift lasts at least the duration of a period, i.e., it cannot start and end at the same period;

3. a trailer cannot have more than one visit to a source at a given period;

4. a shift is carried out at most once.

Assumption 1 makes sense in practice: it is certainly not desirable that a trailer would visit twice a customer in the same shift. This assumption should have a very small impact with respect to optimality on real-life instances. Assumption 2 and 3 are also expected to hold for most shifts, since a period is one-hour long. Actually, in the testbed used for the experiments, these assumptions are almost always satisfied (only three instances contain customers that can be delivered directly from the base within 54 minutes shift). Assumption 4 excludes the assignment of the same time-stamped shift to different drivers. This could be of interest in some specific situations, but, again, one can expect a very limited loss of optimality due to this assumption on real-life instances. We have to stress that Assumption 4 does not prevent from repeating the same tour at different times since a shift is time-stamped.

In the next two subsections, we introduce the notation associated with the shifts and the mathematical formulation.

Notation for shifts

Let us note R d ⊆ R the subset of shifts compatible with driver d ∈ D. Shift r is compatible with driver d if: (1) shift r is included in one of the time windows associated with driver d;

(2) the driving duration of shift r does not exceed the maximal driving duration allowed for driver d; (3) if r is a layover shift, r respects the maximum layover duration associated with driver d.

We note R t ⊆ R the subset of shifts compatible with trailer t ∈ T . A shift r is compatible with trailer t if: (1) all customers visited in shift r are compatible with trailer t;

(2) all sources visited in shift r are compatible with trailer t; (3) for each trip in r, the sum of the minimal delivery quantities of the VMI customers served in the trip does not exceed the capacity of trailer t.

We note R dt the subset of shifts that are simultaneously compatible with driver d ∈ D and with trailer t ∈ T , i.e.,

R dt = R d ∩ R t .
We note c rdt the cost of shift r ∈ R dt when operated by driver d and trailer t. We note st r the starting time of shift r, and et r its ending time (in minutes). st r and et r are included in the time interval covered by the trip.

Notation i ∈ r indicates that customer i ∈ V is visited by shift r. Given i ∈ r, we note st ri the one-hour period during which the delivery of the commodity to customer i starts. Last, we define δ ihr = 1 if i ∈ r and st ri = h, 0 otherwise. We introduce T rip r the set of trips performed during in shift r. Recall that a trip is defined as a sequence of customer visits starting from/ending at a source and/or the base. We note i ∈ s if customer i ∈ V is visited in trip s. We note V s the subset of customers visited in trip s.

Notation introduced for shifts are reported in Table 3. To ease notation, we introduce additional subset notation reported in

V 1 r subset of customer set V 1 visited in shift r V 2 r subset of customer set V 2 visited in shift r V 1 s subset of customer set V 1 visited in trip s V 2 s subset of customer set V 2 visited in trip s R i subset of shifts that visit customer i

Mathematical model

The REC-IRPcan be formulated as an integer linear program as follows. Let us first introduce four sets of decision variables:

       x rdt ∈ {0, 1} (d ∈ D, t ∈ T , r ∈ R dt ) q irt ≥ 0 (r ∈ R, t ∈ T r , i ∈ V r ) q th ≥ 0 (t ∈ T , h ∈ H) S i h ≥ 0 (i ∈ V 1 , h ∈ H)
x rdt is a binary variable equal to 1 if shift r is assigned to driver d and trailer t, 0 otherwise. Variable q irt indicates the quantity delivered to customer i in shift r using trailer t. Variable q th is the load in trailer t at the end of period h. Variable S i h represents the inventory level at VMI customer i ∈ V 1 at the end of period h.

We propose formulation (1)-( 20), with new notation detailed next: minimize d∈D t∈T r∈R dt c rdt x rdt r∈R t∈Tr i∈Vr q irt

(1) subject to compatibility constraints:

r∈K d m t∈T rd x rdt ≤ 1 (d ∈ D, m ∈ St(R d )), π d m ( 2
)
r∈K t m d∈Drt x rdt ≤ 1 (t ∈ T , m ∈ St(R t )), π t m (3) d∈Dr t∈T rd x rdt ≤ 1 (r ∈ R), π r (4)
trailer capacity constraints:

q th ≤ q t,h-1 - r∈R i∈Vr δ ihr q irt + Q t d∈D r∈R dt δ src hr x rdt (t ∈ T , h ∈ H \ {0}), γ th (5) q t0 ≤ q init t - r∈R i∈Vr δ i0r q irt + Q t d∈D r∈R dt δ src 0r x rdt (t ∈ T ), γ t0 (6)
q th ≤ Q t - r∈R i∈Vr δ src+ ihr q irt - r∈R i∈Vr δ start hr (1 -δ src hr )δ ihr q irt (t ∈ T , h ∈ H), γ up th (7) q t,h-1 ≥ r∈R i∈Vr δ src- ihr q irt + r∈R i∈Vr δ end hr (1 -δ src hr )δ ihr q irt (t ∈ T , h ∈ H \ {0}), γ down th ( 8
)
q init t ≥ r∈R i∈Vr δ src- i0r q irt + r∈R i∈Vr δ end 0r (1 -δ src 0r )δ i0r q irt (t ∈ T ), γ down t0 (9) i∈Vs q irt ≤ Q t d∈Drt x rdt (t ∈ T , r ∈ R t , s ∈ T rip r ), γ trs (10)
delivery constraints:

t∈Tr q irt ≥ q min i d∈Dr t∈T rd x rdt (r ∈ R, i ∈ V 1 r ), β ri (11)
r∈R i t∈Tr q irt ≥ q min i (i ∈ V 2 ), β min i ( 12
)
r∈R i t∈Tr q irt ≤ q max i (i ∈ V 2 ), β max i ( 13
)
customer inventory constraints:

S i h-1 + r∈R i t∈Tr δ ihr q irt -d i h = S i h (i ∈ V 1 , h ∈ H \ {0}), σ ih (14)
stock i 0 + r∈R i t∈Tr δ i0r q irt -d i 0 = S i 0 (i ∈ V 1 ), σ i0 (15)
S i h ≤ stock i max (i ∈ V 1 , h ∈ H), σ max ih ( 16
)
and with decision variable definitions:

x rdt ∈ {0, 1} (d ∈ D, t ∈ T , r ∈ R dt ), ( 17
) q irt ≥ 0 (r ∈ R, t ∈ T r , i ∈ V r ), (18) q th ≥ 0 (t ∈ T , h ∈ H), ( 19
) S i h ≥ 0 (i ∈ V 1 , h ∈ H). (20)
We first describe equations that can quickly be explained. Other constraints are explained in details after, as well as new notation introduced for the model. On the extreme right of the equations, we introduce, in red, notation for dual variables. These will be used later in the branch-cut-and-price algorithm.

Fractional objective function (1) is the logistics ratio. Constraints (4) state that every shift is selected at most once in a solution. These constraints stem from Assumption 4 and are useful for inventory constraints (11) below. Constraints (11) ensure that minimal quantities are delivered to VMI customers at each visit. Constraints ( 12) and ( 13) respectively impose a minimal and a maximal accumulated quantities delivered to each call-in customer. Constraints ( 14) and ( 15) avoid stockouts at VMI customers. These constraints remain easy to formulate thanks to Assumption 1. Constraints ( 16) limit the stock level at VMI customers. In these three sets of constraints, the inventory level is checked at the end of the period, once the deliveries have been made and the demand is satisfied.

In Constraints ( 2

(R d ) = {st r , r ∈ R d }.
Note that the time scale for these constraints is the minute. Notation introduced in these constraints, as well as in the next constraints, are recalled in Table 5. 5: Additional notation introduced for shift compatibility management Constraints (3) are equivalent for trailers. We note St(R t ) the set of starting times of shifts in R t . We construct a set K t m for each m ∈ St(R t ) such that K t m = {r ∈ R t : st r ≤ m, et r ≥ m}. Then constraints (3) enforce that two shifts requiring the same trailer cannot overlap. Again, the time scale for these constraints is the minute. Constraints (5)-( 10) are by far the less easily readable. They are introduced to manage trailer capacity. One important aspect of the REC-IRP is that the remaining load in a trailer at the end of a trip is available for the next trip. Due to this specific feature, we cannot model the trailer capacity constraint by simply bounding the quantity delivered on each trip. We proposed to trace the trailer loading levels with variables q th . These constraints also imply introducing notation given in Table 6 and are at the origin of assumptions 2 and 3.

Constraints (5) are active when no source is attained by trailer t during period h. Then, the loading level of the trailer at the end of the period cannot exceed its level at δ src hr 1 if shift r reaches a source at period h, 0 otherwise δ src+ ihr 1 if, during period h, shift r both visits, in this order, a source and customer i, 0 otherwise δ src- ihr 1 if, during period h, shift r both visits, in this order, customer i and a source, 0 otherwise δ start hr 1 if shift r starts at period h, 0 otherwise δ end hr 1 if shift r ends at period h, 0 otherwise Table 6: Additional notation introduced for the trailer capacity management the beginning of the period minus the quantity delivered during the period. Note that when the source is visited, Constraints (5) are equivalent to the sum of Constraints ( 7) and Constraints (8). Constraints (6) consider the initial period h = 0.

When a source is reached by a trailer during period h, Constraints (7) limit the load of the trailer at the end of the period. The value subtracted to Q t on the right-hand side of the constraint evaluates the quantity of the commodity delivered after the source in period h, both by the shift that visits the source and possibly successive shifts. We emphasize that a shift is assumed to last more than one period and that the source is assumed to be visited at most once in a given period. Therefore, at most two shifts can be associated with trailer t at a given period. In such a case, the first shift may include a visit to the source while the second shift includes visits to customers only. The latter starts in this period. If no source is reached, these constraints are dominated by Constraints (5).

When a source is reached by a trailer during period h, Constraints (8) guarantee that the total quantity delivered does not exceed the quantity available before reaching a source. More precisely, the two terms on the right-hand side evaluate the quantity delivered before the source in period h, either by the shift that visits the source or by a shift scheduled before, i.e., that ends at period h. If no source is visited, these constraints are dominated by Constraints (5). Constraints (9) play the same role for the initial period (h = 0). Constraints (10) are only active when shift r is not assigned to a trailer and a driver. It prevents from delivering commodities with shift r in this case. Otherwise, the maximal amount of commodity delivered is restricted by the previous constraints. Note that these constraints are tight.

Modified objective function

In the algorithm presented in Section 6, the objective function is actually slightly modified to favor early deliveries and identify feasible solutions more rapidly. We proceed as follows. We introduce a delivery profit vector c h that indicates the profit raised when delivering a commodity unit to a customer at a specified time h ∈ H (see Section 6.2 for details on the values given to c h ). Note that the actual objective function of the REC-IRP is retrieved if

Linearized objective function

In the algorithm presented in Section 6, the objective function is linearized by following the principle of Dinkelbach's algorithm. Dinkelbach's method solves nonlinear problems with a fractional objective function, a concave numerator and a (non-zero) convex denominator.

Since the linear relaxation of our model fits these conditions, its application to our problem is justified.

To briefly explain how Dinkelbach's algorithm works, let us consider a linear program with a linear fractional objective:

(P ) min x∈X N (x) D(x) ,
with N (x) and D(x) linear and positive functions for all x ∈ X. We note f (Z) the value of the optimal solution of the parametric linear program:

min x∈X (N (x) -Z × D(x)).
One can easily prove that f has a unique root Z * (i.e., f (Z * ) = 0) and that

Z * = min x∈X N (x) D(x) .
To solve the problem (P ), Dinkelbach's algorithm computes Z * applying Newton's algorithm: it starts from a reasonable value Z 0 , and it iterates as follows:

Z k+1 = Z k - f (Z k ) f (Z k ) = N (x * k ) D(x * k )
, where x * k is the optimal solution of the parametric linear program for parameter value Z k . In our case, the objective function of the parametric linear program, for a given parameter value Z, is then:

minimize d∈D t∈T r∈R dt c rdt x rdt -Z r∈R t∈Tr i∈Vr ( h∈H δ ihr c h )q irt ,
To simplify notation, we introduce Z ir = Z h∈H δ ihr c h . Then, the objective function becomes:

minimize d∈D t∈T r∈R dt c rdt x rdt - r∈R t∈Tr i∈Vr Z ir q irt ,
In Section 6.2, we precise how parameter Z varies in the course of the algorithm.

Valid inequalities

The linear relaxation of the model can be enriched by the introduction of the following three classes of valid inequalities related to visits to customers.

For call-in customers, at least one visit is needed:

r∈R i d∈Dr t∈T rd x rdt ≥ 1 (i ∈ V 2 ). α i (21)
For VMI customers, the minimal number of visits can be deduced from the accumulated demand of the customer on the horizon, the inventory capacity of the customer and the maximal capacity of a trailer. Let call n i this value. For h, h ∈ H with h ≤ h , we introduce

D i hh = h l=h d i l the cumulated demand of customer i ∈ V 1
between the beginning of period h and the end of period h . Then,

n i = D i 0H -stock i 0 min{stock max i , max r∈R i ,t∈Tr Q t } .
and the following constraints can be imposed:

r∈R i d∈Dr t∈T rd x rdt ≥ n i (i ∈ V 1 ). α i (22)
Additional valid inequalities can be added to better manage the visit times at VMI customers. For these inequalities, we introduce functions π i : H → H ∪ {-1}, i ∈ V 1 defined as follows:

• if D i hH > stock max i , π i (h) is the smallest period h ≥ h with D i hh > stock max i ; • if D i hH ≤ stock max i , π i (h) = -1.
Then, the valid inequalities are as follows:

r∈R i d∈Dr t∈T rd l=π i (h) l=h δ ilr x rdt ≥ 1 (i ∈ V 1 , h ∈ H : π i (h) > 0). α ih (23)

Artificial variables

In a branch-cut-and-price algorithm, the model presented in Section 3.2 will be solved over a restricted number of td-shifts. As a drawback, the restricted model may not admit a feasible solution even if the problem does. To prevent this situation we introduce some so-called artificial variables. In the following we detail the modifications on the objective function and on some constraints of the model due to the inclusion of these variables. We introduce four types of non-negative artificial variables. Basically, all these variables enable some deliveries to be skipped:

(1) y i , to deactivate Valid Inequalities ( 21) or ( 22) for customers i ∈ V;

(2) w i , to deactivate Valid Inequalities (23) for customers i ∈ V 1 ;

(3) oos ih , to allow customer i ∈ V 1 to get out of stock at period h;

(4) drop i , to allow not to visit customer i ∈ V 2 .

Artificial variables are included in the objective function multiplied by large values M i (i ∈ V).

The objective function becomes:

minimize d∈D t∈T r∈R dt c rdt x rdt - r∈R t∈Tr i∈Vr Z ir q irt + i∈V M i y i + i∈V 1 M i w i + i∈V 1 h∈H M i oos ih + i∈V 2 M i drop i .
Constraints ( 12), ( 14), ( 21), ( 22), ( 23) become:

r∈R i t∈Tr q irt ≥ q min i -drop i (i ∈ V 2 ), S i h-1 + r∈R i t∈Tr δ ihr q irt -d i h = S i h -oos ih (i ∈ V 1 , h ∈ H), r∈R i d∈Dr t∈T rd x rdt ≥ 1 -y i (i ∈ V 2 ), r∈R i d∈Dr t∈T rd x rdt ≥ n i (1 -y i ) (i ∈ V 1 ), r∈R i d∈Dr t∈T rd π i (h) l=h δ ilr x rdt ≥ 1 -w i (i ∈ V 1 , h ∈ H : π i (h) > 0).
In the following of the paper, we will say that a mathematical program is feasible when it admits a solution with all artificial variables equal to zero. Otherwise we will say that the mathematical program is infeasible.

Master problem

In the rest of the paper we will call Master Problem the linear relaxation of the mathematical model defined by ( 1)-(20) with the modifications described in the preceding subsections: modified and linearized objective function, valid inequalities and artificial variables added. The Master Problem will be denoted M P .

Column generation framework

Because of its number of variables, M P cannot be solved directly with a linear programming solver. In this section we describe how a column generation mechanism can be applied to progressively enrich a restricted version of M P .

Let us introduce some notation that will be used in this section. We note T DR the complete set of triplets (r, d, t), with r ∈ R, d ∈ D r and t ∈ T rd . T DR is the set of all the td-shifts. A variable x rdt is present in M P for each triplet in T DR. Similarly, we note T R the set of pairs (r, t), with r ∈ R and t ∈ T r , and DR the set of pairs (r, d), with r ∈ R and d ∈ D r .

We call RM P k the restricted M P at a step k of the algorithm and we explain how RM P k is constructed from a subset T DR k ⊆ T DR of td-shifts. Together, we explain how sets T R k ⊆ T R, DR k ⊆ DR and R k ⊆ R are constructed. To this aim, we equivalently explain how we proceed when a td-shift is added to T DR k .

When a new td-shift (r, d, t) is added to T DR k :

• variable x rdt is added to the model;

• if (r, t) is not in T R k , (r, t) is added to T R k , variables q irt are
added to the model for all i ∈ r, and Constraints (10) are added for pair (r, t) for every trip s ∈ T rip r . Furthermore, if the starting time of r is not in St(R t k ), Constraint (3) is added for pair (t, st r ); 2) is added for pair (d, st r ); 4) is added to the model for r, as well as Constraint (11) for every i ∈ V 1 r .

• if (r, d) is not in DR k , (r, d) is added to DR k and, if the starting time of r is not in St(R d k ), Constraint (
• if r is not in R k , r is added to R k , Constraint (
In our algorithm, see Section 6, we start with an initial set of promising td-shifts T DR 0 and regularly insert new td-shifts to drive the search towards feasible/better solutions. Since the modifications described above impact RM P k when new td-shifts are added, standard column generation cannot be applied. We recall in the next subsection how to adapt column generation for the simultaneous generation of columns and cuts. Then, we explain how it is applied to the REC-IRP.

About simultaneous cut and column generation

We note D the dual of M P and D k the dual of RM P k . When the addition of columns to a restricted master problem also implies the addition of cuts, it raises new difficulties that invalidate the standard column generation procedure (see [START_REF] Feillet | A note on branchand-cut-and-price[END_REF]). Indeed, the usual optimality criterion used in column generation relies on the fact that when no column of negative reduced cost exists, the current dual solution (obtained from the solution of the restricted master problem) is feasible for the dual of the master problem. In case of simultaneous cut and column generation it does not hold, since the cuts, that are missing in the restricted master problem, translate into missing variables into the dual of the master problem. Then, no guarantee of feasibility for the current dual solution with regards to the dual of the master problem is possible, without defining a value for these missing variables. [START_REF] Feillet | A note on branchand-cut-and-price[END_REF] propose to tackle this difficulty with a 3-step methodology:

1. from optimal solution X * k of RM P k , construct a feasible solution X of M P , with cost

z M P = z * RM P k ; 2. from optimal solution Π * k of D k , construct a (not necessarily feasible) solution Π of D, with cost z D = z * D k , i.e.
, z D = z M P ; 3. if Π is feasible, the optimality criterion is met: stop the algorithm; otherwise, at least one constraint of D is violated; add to T DR k one td-shift from T DR \ T DR k for which a constraint is violated and pass to the next iteration.

Unfortunately, no simple general rule seems to exist to construct solutions X and Π. In the next section, we see how this scheme can be applied to the REC-IRP.

Application to the REC-IRP

Let us start with a solution X * k of RM P k and the associated dual solution Π * k of D k . The first step of the methodology consists in constructing a feasible solution X of M P , having the same cost as X * k . We simply fix to 0 all variables x rdt for (r, d, t) ∈ T DR \ T DR k and q irt for (r, t) ∈ T R \ T R k , and keep other variables to their value in X * k . X is feasible for M P and z M P = z * RM P k . The next step is then to construct a not necessarily feasible solution Π for the dual D, having the same cost as Π * k . For this step, we need to consider D carefully. D is composed of eight series of constraints, one for each type of primal variables (x rdt , q irt , q th , S i h ), plus one for each type of artificial variables (y i , w i , oos ih , drop i ). The differences between D k and D are the followings:

• the constraints derived from primal variables x rdt with (r, d, t) ∈ T DR \ T DR k and primal variables q irt with (r, t) ∈ T R \ T R k are missing in D k

• some variables π d m , π t m , π r , γ trs and β ri , respectively derived from primal constraints (2), (3), ( 4), ( 10) and ( 11) are missing in D k .

We define Π as follows. First, all variables already defined in D k keep the same value as in Π * k . Then, variables π d m , π t m and π r not in D k are set to 0. The remaining variables to be fixed are variables γ trs for (r, t) ∈ T R \ T R k and s ∈ T rip r , and

β ri for r ∈ R \ R k and i ∈ V 1 r .
Before defining these variables, let us underline that, independently of the values to be assigned, several important conditions already hold. First, z D = z * D k . Indeed, variables γ trs and β ri derive from primal Constraints (10) and ( 11) that have no fixed term, and thus do not appear in the dual objective function. Second, Π satisfy all the dual constraints already present in D k . Indeed, none of the remaining variables appear in these constraints. Thus, we now have to focus on finding values for the remaining variables and evaluating the feasibility of the missing constraints.

We write the missing constraints below. In order to facilitate the understanding of the constraints, we define by convention all dual variables nonnegative, except variables σ ih (i ∈ V 1 , h ∈ H) associated with Constraints ( 14) and ( 15) that are unsigned.

For (r, d, t) ∈ T DR \ T DR k , the dual constraint in D derived from primal variable x rdt is:

c rdt + i∈V 1 r q min i β ri + {m∈St(R d ):r∈K d m } π d m + {m∈St(R t ):r∈K t m } π t m + π r - i∈Vr α i - i∈V 1 r {h∈H:h≤st ri ≤π i (h)} α ih -Q t h∈H δ src hr γ th -Q t s∈T ripr γ trs ≥ 0. ( 24
) For (r, t) ∈ T R \ T R k , i ∈ V 1 r
, and s the trip of customer i in td-shift r, the dual constraint derived from primal variable q irt is:

h∈H (δ src+ ihr + δ start hr (1 -δ src hr )δ ihr )γ up th + h∈H (δ src- ihr + δ end hr (1 -δ src hr )δ ihr )γ down th + h∈H δ ihr γ th + γ trs -β ri - h∈H δ ihr σ ih ≥ Z ir .
Last, for (r, t) ∈ T R \ T R k , i ∈ V 2 r , and s the trip of customer i in td-shift r, the dual constraint derived from primal variable q irt is:

h∈H (δ src+ ihr + δ start hr (1 -δ src hr )δ ihr )γ up th + h∈H (δ src- ihr + δ end hr (1 -δ src hr )δ ihr )γ down th + h∈H δ ihr γ th + γ trs + β max i -β min i ≥ Z ir .
For the sake of simplifying these constraints, we introduce :

δ up ihr = δ src+ ihr + δ start hr (1 -δ src hr )δ ihr , and
δ down ihr = δ src- ihr + δ end hr (1 -δ src hr )δ ihr .
Note that δ up ihr = 1 in two special cases: a source and customer i are visited in this order at period h in shift r, or shift r starts at period h and visits customer i at the same period but does not visit a source. Similarly, δ down ihr = 1 if customer i and a source are visited in this order at period h in shift r, or shift r ends at period h and visits customer i at the same period but does not visit a source.

Given (r, t) ∈ T R \ T R k , i ∈ V 1
r , and s the trip of customer i in td-shift r, the dual constraint derived from primal variable q irt is then:

h∈H δ up ihr γ up th + h∈H δ down ihr γ down th + h∈H δ ihr γ th + γ trs -β ri - h∈H δ ihr σ ih ≥ Z ir , (25)
and given (r, t) ∈ T R \ T R k , i ∈ V 2 r , and s the trip of customer i in td-shift r, the dual constraint derived from primal variable q irt is:

h∈H δ up ihr γ up th + h∈H δ down ihr γ down th + h∈H δ ihr γ th + γ trs + β max i -β min i ≥ Z ir . ( 26
)
When defining the remaining dual variables in Π, we thus only have to consider the feasibility of Constraints ( 24), ( 25) and ( 26). Before providing a formal proof, we give some intuition on how we proceed. Focusing on Constraints ( 25) and ( 26), it is easy to see that choosing large values for variables γ trs and small values for variables β ri allows satisfying these constraints. However, having large values for variables γ trs and small values for variables β ri decreases the chances that Constraints (24) are satisfied. We propose to select values that ensure the feasibility of Constraints ( 25) and ( 26), but with variables γ trs as small as possible and variables β ri as large as possible, so that the chances, that Constraints (24) are feasible, are increased. However, clearly, looking at Constraints (25) a compromise has to be done between the values of variables γ trs and β ri . We state and solve this compromise with Property 1.

Property 1. Let consider r ∈ R \ R k , and assume that variables γ trs are set for all t ∈ T r , s ∈ T rip r and that variables β ri are set for all i ∈ V 1 r . Let γ 1 be the minimal value of all these variables. Let γ 2 be the minimal gap between the left-hand side and the right-hand side for all Constraints (26) associated with r. Let finally assume that γ = min(γ 1 , γ 2 ) > 0. Then, subtracting γ to all the variables mentioned above increases the chances that dual solution Π is feasible.

Proof. First, the only constraints impacted by the values of the variables mentioned in the property are the constraints (24), ( 25) and ( 26) associated with r. It is easy to see that the subtraction of γ does not impact Constraints (25) and that Constraints (26) remain feasible. Given a td-shift (r, d, t), the impact on Constraint (24) is the addition to the left-hand side of term -( i∈V 1 r q min i )γ + |T rip r |Q t γ.

i∈V 1 r q min i is a lower bound on the quantities delivered on shift r, while |T rip r |Q t is an upper bound. It proves that the term added to Constraint (24) is nonnegative, which concludes the proof.

We now detail how values are fixed for γ trs for (r, t) ∈ T R \ T R k and s ∈ T rip r , and

β ri for r ∈ R \ R k and i ∈ V 1
r . Two cases are considered. First, r ∈ R k . In this case, variables β ri are already fixed and we need to define γ trs for all t such that (r, t) ∈ T R \ T R k and for all s ∈ T rip r . Clearly, the smallest γ trs is the best choice with respect to Constraint (24). So, we fix γ trs to the minimal value leaving Constraints ( 25) and ( 26) feasible:

γ 1 = max i∈V 1 s Z ir - h∈H δ up ihr γ up th - h∈H δ down ihr γ down th - h∈H δ ihr γ th + β ri + h∈H δ ihr σ ih , γ 2 = max i∈V 2 s Z ir - h∈H δ up ihr γ up th - h∈H δ down ihr γ down th - h∈H δ ihr γ th -β max i + β min i , γ trs = max(γ 1 , γ 2 , 0).
Second, r / ∈ R k . In this case, none of the variables γ trs for t ∈ T r and q irt for i ∈ V 1 r are fixed. In view of Property 1, we favor small values for these variables:

γ 1 = max i∈V 1 s Z ir - h∈H δ up ihr γ up th - h∈H δ down ihr γ down th - h∈H δ ihr γ th + h∈H δ ihr σ ih , γ 2 = max i∈V 2 s Z ir - h∈H δ up ihr γ up th - h∈H δ down ihr γ down th - h∈H δ ihr γ th -β max i + β min i , γ trs = max(γ 1 , γ 2 , 0).
Note that γ trs is computed assuming β ri variables set to 0. Note that with variables γ trs fixed, variables β ri are as large as possible.

β ri = min
With this definition, the dual solution Π has the same cost as the dual solution Π * k of the restricted dual program and all constraints are feasible except possibly Constraints (24) for some td-shifts (r, t, d) ∈ T DR \ T DR k . The goal of the pricing problem is to search for td-shifts (r, d, t) such that this constraint is not satisfied, i.e.:

c rdt + i∈V 1 r q min i β ri + {m∈St(R d k ):r∈K d m } π d m + {m∈St(R t k ):r∈K t m } π t m - i∈Vr α i - i∈V 1 r {h∈H:h≤st ri ≤π i (h)} α ih -Q t h∈H δ src hr γ th -Q t s∈T ripr γ trs < 0. ( 27
)
In the following, we will call the left-hand side of ( 27) the reduced cost of the td-shift.

In our algorithm, presented in Section 6, we developed two procedures to insert new td-shifts with negative reduced cost in RM P k . The first procedure considers existing tdshifts and tries to generate new td-shifts by changing the trailer, the driver and/or the time-stamp of these td-shifts. In this procedure, the reduced cost is evaluated exactly. This procedure is detailed in Section 6.3.3.

The second procedure aims at generating completely new td-shifts. It is classically based on the solution of a shortest path problem with resource constraints by dynamic programming. However, because of the difficulty of the pricing, reduced costs are only evaluated heuristically in what we call an implicit-relaxation-based dynamic programming procedure. This procedure is described in the next section.

Implicit-relaxation-based dynamic programming procedure

The implicit-relaxation-based dynamic programming procedure is called with a few parameters that define the subset of td-shifts that is explored: the td-shifts start at time m 0 (in minutes), have to be finished at time m 1 (in minutes), use driver d 0 and trailer t 0 . For short, the notation introduced below is not indexed by these parameters. In Section 6, we detail how these parameters are fixed and when the procedure is called.

The dynamic programming procedure is named implicit-relaxation-based because it is implicitly applied on a relaxation of M P . By implicitly, we mean that this relaxation is not implemented. Instead, when RM P k is solved, a dual vector consistent with the dual space of this relaxation is constructed. Using this dual information, we derive a (heuristic) reduced-cost condition that is checked for the generation of new td-shifts. We call pricing problem, the problem of identifying td-shifts that do not satisfy this condition. This condition presents the advantage to be easier to consider in a dynamic programming scheme than condition (27). In particular, as it will be seen, the relaxation is constructed so that waiting is of no interest when td-shifts are generated in the pricing problem. However, contrary to standard column generation, there is no guarantee that the actual reduced cost of the td-shift in M P is negative. Hence, this reduced cost has to be computed exactly after the pricing problem is solved, which is done using the formula presented in Section 4.

We now describe the relaxation, the construction of the dual vector, the reduced-cost condition and the dynamic programming algorithm developed to identify td-shifts that do not satisfy this condition.

Implicit relaxation and dual vector

In this section, we consider the different families of constraints that are relaxed and explain how they are relaxed and with which consequences for the dual vector.

(23): These constraints enforce visiting at least once a customer during a given period if its consumption during this period is greater than its inventory capacity. Let consider i ∈ V1 . These valid inequalities are initially written:

r∈R i d∈Dr t∈T rd l=π i (h) l=h δ ilr x rdt ≥ 1 -w i (h ∈ H : π i (h) > 0)
We introduce set

H i = {h ∈ H : [m 0 , m 1 ] ∩ [h, π i (h)] = ∅ 1
} and we focus on the constraints associated with periods h ∈ H i . In these constraints, we change the term

l=π i (h) l=h δ ilr x rdt to l∈[m 0 ,m 1 ]∪[h,π i (h)] δ ilr x rdt for all triplets (r, d, t) such that d = d 0 , t = t 0 . This new formulation is a relaxation because [h, π i (h)] ⊆ [m 0 , m 1 ] ∪ [h, π i (h)].
Formally, the constraints become:

r∈R i d∈Dr t∈T rd l=π i (h) l=h δ ilr x rdt ≥ 1 -w i (h ∈ H \ H i : π i (h) > 0) r∈R i (d,t)∈Dr×T rd \{(d 0 ,t 0 )} l=π i (h) l=h δ ilr x rdt + r∈R i l∈[m 0 ,m 1 ]∪[h,π i (h)] δ ilr x rd 0 t 0 ≥ 1 -w i (h ∈ H i : π i (h) > 0)
In the pricing problem, the contribution of these constraints to the reduced cost of td-shifts that visit i is h∈H i α ih . The interest of this relaxation relies on the fact that this term does not depend on the visiting time of i. For short, we note: h∈H i α ih = αi .

Constraints (5) and (6):

These constraints are related to the trailer capacities. For these constraints, we apply a surrogate relaxation. We call h 0 and h 1 the periods (in hour) that contain m 0 and m 1 , respectively. For t = t 0 , we replace the constraints defined for h ∈ {h 0 , . . . , h 1 } with their sum:

q t 0 h 1 ≤ q t 0 ,h 0 -1 - h∈{h 0 ,...,h 1 } r∈R i∈Vr δ ihr q irt 0 + Q t 0 d∈D r∈R dt 0 n src r x rdt 0
with n src r equal to the number of visits to a source in r on interval [h 0 , h 1 ]. If h 0 = 0, q t,h 0 -1 is replaced by q init t in the above expression.

We note γt 0 the dual associated with this constraint. In the pricing problem, the contribution of this constraint to the reduced cost of a td-shift is γt 0 each time the td-shift visits a source. As the constraint is not in M P , we need to define a value for γt 0 . We arbitrarily set:

γt 0 = 1 h 1 + 1 -h 0 h=h 1 h=h 0 γ t 0 h .
The rationale is that every time the source is visited, the td-shift will receive a dual price equal to the average value of dual variables γ t 0 h instead of the value of the γ t 0 h variable associated with the specific period of the visit. Again, the main advantage is that the new dual variable γt 0 is not indexed by time.

Constraints (7) to (9):

These constraints also concern trailer capacities. In these constraints, we set to 0 the coefficient of all q irt 0 variables. It defines a relaxation because the coefficients where initially all nonnegative on the left-hand side and nonpositive on the right-hand side of the constraints. The interest here is to get rid of the burden of carrying dual variables γ up t 0 h and γ down t 0 h , that we suspect not being very influential.

Constraints ( 14) and ( 15):

These constraints prevent from stockouts at VMI customers. Again, these constraints are simplified using a surrogate relaxation. We reuse notation h 0 and h 1 . Given i ∈ V 1 , we replace the constraints defined for h ∈ {h 0 , . . . , h 1 } with their sum:

S i h 0 -1 + r∈R i t∈Tr q irt - h∈{h 0 ,...,h 1 } d i h = S i h 1 - h∈{h 0 ,...,h 1 } oos ih
If h 0 = 0, S i h 0 -1 is replaced by stock i 0 in the above expression. We note σi the dual variable associated with this constraint. As the constraint is not in M P , we need to define a value for this variable. Again, we define σi as the average value of variables σ ih on the interval:

σi = 1 h 1 + 1 -h 0 h=h 1 h=h 0 σ ih .
Intuitively, the dual price associated with a visit of customer i is averaged on the different periods of interval {h 0 , . . . , h 1 }.

Evaluation of the reduced cost

We now explain how we evaluate the reduced cost of the td-shifts r searched in the pricing problem. Given i ∈ V 1 visited by r and s ∈ T rip r the trip of r in which i is visited, the constraint derived from primal variable q irt 0 in the relaxed formulation is: γt 0 + γ t 0 rs -β ri -σi ≥ Z ir and given i ∈ V 2 visited by r, the dual constraint derived from primal variable q irt 0 is:

γt 0 + γ t 0 rs + β max i -β min i ≥ Z ir
Because we are only interested in finding td-shifts not in R k , we assume that variable γ t 0 rs is not known. Following the approach presented in Section 4.2, we set:

γ 1 = max i∈V 1 r Z ir -γt 0 + σi γ 2 = max i∈V 2 r Z ir -γt 0 -β max i + β min i γ t 0 rs = max(γ 1 , γ 2 , 0)
Then, to avoid managing values γ trs for t = t 0 , we set:

β ri = 0
Note that following the scheme developed in Section 4.2, we would have set β ri = min t∈T γt + γ trs -σ i -Z, but 0 also ensures the validity of dual constraints associated with dual variables q irt 0 .

The goal of the pricing problem is to find td-shifts (r, d 0 , t 0 ), starting at time m 0 and ending not after time m 1 , such that:

c rd 0 t 0 + {m∈St(R d 0 k ):r∈K d 0 m } π d 0 m + {m∈St(R t 0 k ):r∈K t 0 m } π t 0 m - i∈Vr α i - i∈V 1 r αi -Q t 0 n src r γt 0 -Q t 0 s∈T ripr γ t 0 rs < 0 (28)
It is noteworthy that this implicit relaxation can also be interpreted as a stabilization method. Indeed, when solving RM P k with the simplex method, dual variables receive extreme values. An important drawback with dual variables indexed by time is that a dual variable can receive a large value while the dual variables in adjacent periods receive zero. Then, with the column generation mechanism, the dual value can be restricted for the specific period that received a large value, but this value can (repeatedly) be moved to a close period. Then many iterations are potentially required to attain a good evaluation of the dual prices. One advantage of our implicit relaxation scheme is to average the value of these dual variables.

Dynamic programming algorithm

The pricing is addressed through dynamic programming as an elementary shortest path problem with resource constraints. The algorithm follows the approach described in [START_REF] Feillet | An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems[END_REF], except for the label definition, resource extension functions and dominance rule. Before detailing these features of the algorithm, we introduce additional notation (see Table 7) that completes the notation already reported in Tables 1 and2

i ∈ V [a ik , b ik ] k th time window for customer i ∈ V (1 ≤ k ≤ tw i )
Table 7: Notation (other data) Labels are defined with ten attributes plus a vector of size |V|, L = (ending vertex, cost, time, driving time, minimal quantity, gamma, sources, trips, layover, layover customer,visited[]), that can be interpreted as shown in Table 8. Denoting base the index of the base in the location set, the algorithm is initialized with label:

L = (base, m 0 +rest d 0 -1 m=m 0 π d 0 m , m 0 , T d 0 , 0, 0, 0, 0, 0, 0, [0, . . . , 0]) .
Resource extension functions are defined in Tables 9 and10. In these tables, label L = (i, c, t, dt, minq, gamma, nbs, nbt, lay, laycust, visited[]) is extended to vertex j ∈ V to form label L = (i , c , t , dt , minq , gamma , nbs , nbt , lay , laycust , visited []). The last column in Table 9 indicates if the formula applies to the source (s), the base (b) or both (s+b). The last column in Table 10 is empty when the formula always applies. (lay.) indicates that it applies when a layover is carried out, and (no lay.) stands for the case with no layover.

The inclusion of a layover is somewhat complex. Three cases are considered. When the waiting time before the opening of a time window exceeds the layover duration, a layover is imposed (formally when t ≤ b jk ), it is inserted before, otherwise it will be inserted when L is extended, after j.

+ s i + t ij + t layover d 0 ≤ a jk with k = min{l : 1 ≤ l ≤ tw j and t + s i + t ij ≤ b jl },
The extension of label L to j is also subject to some few feasibility tests: i) If j is a customer, the extension is not allowed if visited[j] = 1; ii) the extension with a layover is not allowed when lay = 1 (note that only the first case above is concerned); iii) a time window k computed as explained in Table 10 must exist (t+s i +t ij ≤ b jtw j ). If j is the base or a source, condition t ≤ m 2 is required. If j is the base the combination lay = 1 and laycust = 0 is not allowed. In all cases (customer, source or base), the remaining driving time should be nonnegative: dt ≥ 0, and the minimal delivered quantity mq should not exceed Q t 0 .

Regarding the dominance rule, we propose very simplified (non-exact) conditions. Given two labels L 1 and L 2 ,

L 1 = (i 1 , c 1 , t 1 , dt 1 , minq 1 , gamma 1 , nbs 1 , nbt 1 , lay 1 , laycust 1 , visited 1 []) dominates L 2 = (i 2 , c 2 , t 2 , dt 2 , minq 2 , gamma 2 , nbs 2 , nbt 2 , lay 2 , laycust 2 , visited 2 []) if: 1. i 1 = i 2 2. c 1 ≤ c 2 3. t 1 ≤ t 2 4. dt 1 ≥ dt 2 end. ver. i = j (s+b) c = c+dc d 0 ×(s i +t ij )+tc t 0 ×d ij + t -1 m=t π t 0 m + t +rest d 0 -1 m=t+rest d 0 π d 0 m - Q t 0 γt 0 (s) c = c+dc d 0 ×(s i +t ij )+tc t 0 ×d ij + t m=t π t 0 m + t +rest d 0 m=t+rest d 0 π d 0 m (b) time t = t + s i + t ij (s+b) dr. time dt = dt -t ij (s+b) min. qtty mq = 0 (s+b) gamma gamma = 0 (s+b) sources nbs = nbs + 1 (s) nbs = nbs (b) trips nbt = nbt (s+b) layover lay = lay (s+b) lay. cust. laycust = laycust (s+b) visited[] visited [u] = visited[u] for all u ∈ V
Table 9: Resource extension functions to a source or to the base

end. ver. i = j cost c = c+dc d 0 ×(t -t)+tc t 0 ×d ij + t -1 m=t π t 0 m + t +rest d 0 -1 m=t+rest d 0 π d 0 m - α j -αj + gamma × Q t 0 -gamma × Q t 0 (no lay.) c = c+dc d 0 ×(t -t)+tc t 0 ×d ij + t -1 m=t π t 0 m + t +rest d 0 -1 m=t+rest d 0 π d 0 m - α j -αj +gamma×Q t 0 -gamma ×Q t 0 -dc d 0 ×t layover d 0 +c layover d 0 (lay.) time t = max(t+s i +t ij , a jk ) with k = min{l : 1 ≤ l ≤ tw j and t+ s i + t ij ≤ b jl } (no lay.) t = max(t + s i + t ij + t layover d 0 , a jk ) with k = min{l : 1 ≤ l ≤ tw j and t + s i + t ij + t layover d 0 ≤ b jl } (lay.) dr. time dt = dt -t ij (no lay.) dt = min(T d 0 , T d 0 + dt -t ij ) (lay.) min. qtty mq = mq + q min j gamma gamma = max(gamma, Z -γt 0 + σj ) if j ∈ V 1 , gamma = max(gamma, Z -γt 0 -β max j + β min j ) if j ∈ V 2 sources nbs = nbs trips nbt = nbt + 1 if i / ∈ V, nbt = nbt otherwise layover lay = lay (no lay.) lay = 1 (lay.) lay. cust. laycust = laycust if j is not a layover customer, laycust = 1 otherwise visited[] visited [u] = visited[u] for all u ∈ V \ {j}, visited [j] = 1 Table 10: Resource extension functions to a customer j ∈ V c 1 + (gamma 1 -gamma 2 ) × Q t 0 ≤ c 2
The last condition can be explained as follows. When gamma 1 ≤ gamma 2 , the condition is necessary satisfied from the second one. When gamma 1 > gamma 2 the reduced cost c 1 benefits from this difference compared to c 2 . Indeed, the term Q t 0 s∈T ripr γ t 0 rs is subtracted in the reduced cost evaluation. The maximum future benefit that could be obtained from this term for label L 2 and not be obtained to L 1 is (gamma 1 -gamma 2 )×Q t 0 .

Increasing horizon heuristic

In order to deal with the large time horizon of the instances provided for the ROADEF/EURO challenge, we designed an iterative increasing horizon heuristic. The principle of this heuristic is to progressively construct a solution from the beginning to the end of the planning horizon.

We divide the time horizon in three parts. The first part goes from the beginning of the horizon to, let say, period H 1 . This is the active part where the algorithm takes decisions, i.e., determines the td-shifts that will be part of the solution. In order not to make blind choices, constraints (consumption of VMI customers, call-in orders, time windows, etc.) are taken into account until period H 2 > H 1 . Finally, to reduce the complexity of the problem, consumption of customers from H 2 + 1 to H is not considered.

When the algorithm has computed the solution in [0, H 1 ], the values of H 1 and H 2 are increased and the procedure is repeated. The active part of the horizon is still [0, H 1 ], but now it is wider. That is why we call this approach a increasing horizon heuristic instead of a more classic but misleading rolling horizon algorithm.

The column generation framework presented in the previous section should be integrated in a branching scheme in order to be able to solve the REC-IRP. However, the large size of instances proposed for the ROADEF/EURO challenge makes this option irrelevant. As a consequence, we select td-shifts,that will be part of the solution, using a variable fixing procedure. This heuristic will be described in Section 6.4.

In order to reduce the computation times, we introduced in our procedure what we called the catalog CAT . CAT is a long-term memory which stores promising templates identified so far by the algorithm. Templates are sequences of loading and delivery operations such that the time between each pair of operations is determined, but the starting time has not been set yet. More precisely, when we assign a starting time to a template it becomes a shift, while assigning a driver and a trailer to a shift makes it a td-shift.

The determination of a td-shift from a template and the exact evaluation of its reduced cost can be done very efficiently and in a much quicker way than calling the implicitrelaxation-based dynamic programming procedure presented in Section 5. As soon as a td-shift is generated (with constructive heuristics explained in Section 6.3.1, or with the dynamic programming pricing procedure presented Section 4), the corresponding template is added to CAT . It will be considered to generate new td-shifts by assigning to the template an effective starting time, a trailer and a driver. Note than we can obtain different shifts from the same template considering different starting times. Similarly, we can obtain different td-shifts from the same shift taking into account different trailer-driver pairs.

Algorithm Solve RM P and apply dynamic programming (Section 5) to add new templates to CAT 6:

Solve RM P and, if needed, recover feasibility (Section 6.5)

7: repeat 8:
Solve RM P and add new td-shifts to RM P from CAT (Section 6.3.4)

9:

until P cat nbSearch repetitions has been performed 10:

Over the horizon [0, H 1 ] fix x rdt variables to 1 in RM P (Section 6.4)

11:

Solve RM P and, if needed, recover feasibility (Section 6.5)

12:

Reduce RM P if it contains too many variables (Section 6.6)

13:

Increment H 1 and H 2 14:

Repopulate CAT 15: end while 16: Improve solution The framework of the solution procedure is presented in Algorithm 1. First, CAT is initialized with a constructive heuristic presented in Section 6.3.1 (Line 1) and the values of H 1 and H 2 are set to P rh short and P rh long respectively (Line 2). The first RM P is initialized with a set of td-shifts (Line 3). Modifications generated to RM P in order to deal with the increasing horizon are described in Section 6.1.

The goal of the main loop of Algorithm 1 (Lines 4 to 15) is to generate td-shifts over horizon [0, H 2 ] and to fix td-shifts during the current active part of the horizon, namely [0, H 1 ]. New templates are added to CAT with the dynamic programming pricing procedure presented in Section 5 (Line 5). This is detailed in Section 6.3.2. Whenever RM P is not feasible (i.e., based on definition introduced in Section 3.6, the optimal solution of RM P has non-zero artificial variables), the procedure presented in Section 6.5, is applied to recover feasibility, namely an attempt is made to obtain a solution with all artificial variables set to zero. Then, RM P is solved P cat nbSearch times. At each step, new td-shifts are generated from CAT using the procedure presented in Section 6.3.4 to improve the quality of the td-shifts considered in RM P .

To generate integer solutions we apply the variable-fixing procedure presented in Section It sets x rdt variables to 1 and determine the td-shifts that are included the solution. More precisely, x rdt variables of the current RM P associated with td-shifts starting in [0, H 1 ] are considered as binary. A mixed integer version of RM P is then solved. Then, variables x rdt that equal to 1 are fixed and the corresponding td-shifts are part of the final solution.

It may happen that once variables have been fixed, RM P becomes infeasible. In such a case, the procedure described in Section 6.5 is applied again to recover feasibility.

Then, the horizon is increased and values H 1 and H 2 are incremented by a value P rh + (or set to H if the increment would exceed H). At this point we evaluate the size of RM P , i.e., the td-shifts considered by the current RM P . If it is too large, some variables are excluded, as detailed in Section 6.6. Finally, the catalog CAT is repopulated calling the implicit-relaxation-based dynamic programming procedure.

Once the whole horizon has been considered, the algorithm provides with a solution for the REC-IRP. A final step, presented in Section 6.7, is then applied to try to improve its quality.

Restricted master problem on a limited time horizon

Since we consider a restricted time horizon, the RM P has to be defined only on a part of the whole horizon. Here we present the modifications we apply to RM P .

• for all i ∈ V 1 , for all periods h > H 2 :

variable S i h is unbounded; -Constraint ( 16) is relaxed; if π i (h) ≥ H 2 , the corresponding Valid inequality (23) is relaxed;

• for all i ∈ V 1 :
the value of n i involved in Valid inequalities (22), is computed over [0, H 2 ];

• for all i ∈ V 2 such that earliest time for delivery is in [H 2 + 1, H]:

-Constraints ( 12) and ( 13) are relaxed; -Valid inequality ( 21) is relaxed.

Management of the objective function

One of the specific feature of the REC-IRP is the fractional objective function that we manage using the Dinkelbach's algorithm (Section 3.4). In this section we explain the implementation details related to our procedure and how we modify the objective function to help find feasible solutions.

Application to the increasing horizon

During the resolution of the REC-IRP with the increasing horizon heuristic we initially set the value of Z, the parameter used to linearize the objective function, to Z 0 = 0.01.

Then, the value of Z is updated at each iteration of the algorithm, before calling the fixing procedure described in Section 6.4. The ratio of the objective function is updated P ratio iterM ax times according to the Dinkelbach's rule (3.4) and, after each update, the RM P is solved again. Whenever the current value of Z -N (x) D(x) is lower than a parameter P ratio eps , the procedure stops.

Note that this approach is consistent with the implementation proposed in [START_REF] Garaix | Optimization of occupancy rate in dial-a-ride problems via linear fractional column generation[END_REF], while [START_REF] Archetti | Minimizing the logistic ratio in the inventory routing problem[END_REF] only update Z when the parametric linear program is completely solved (by column generation).

Early delivery

During our preliminary experiments, we noticed that, on some instances, the algorithm worked correctly until a certain time period in the horizon. Then it encounters difficulties to find td-shifts to be part of the solution and to build a solution over the whole horizon. We guessed that this problem was due to a lack of anticipation in serving customers that would eventually run out of stock. Once such a problem appears, the procedure tries to insert a particular td-shift in the horizon unsuccessfully since it is blocked by other already fixed td-shifts.

To circumvent this difficulty, we slightly modify the objective function introducing the delivery profits c h . c h indicates the profit raised when delivering a unit of goods to a customer at a specified time h ∈ H. The expected benefit was to favor early delivery, and so, help finding feasible solutions. The values of c h are as follows:

• c h = 1 for h ≤ H 1 ; • c h = 0 for h ≥ H 2 ;
• c h decreases linearly between H 1 and H 2 .

Management of the catalog CAT

Initialization of the catalog CAT

CAT is first initialized with basic templates, i.e., templates including a visit to a source and a visit to only one customer. Obviously, only feasible template are kept into CAT . By feasible, we mean a template from which feasible td-shifts can be obtained. two templates generation procedures, one based on driver's time windows and one based on customer requirements are executed. These procedures are detailed thereafter. The procedure based on customer requirements is first executed once, then the one based on driver's time windows is executed at most P initCat nbIter times, and stops when CAT contains more than P initCat sizeM ax templates.

Initialization procedure based on driver's time windows

This procedure is based on the exploration of the set of driver time windows. It may happen that a driver time window is as wide as the whole horizon. This represents the availability of the driver during the entire planning horizon. Since two shifts assigned to the same driver d must be separated by at least rest d units of time, the driver time window is divided into smaller time windows in order to: (1) allow to deliver a reasonable number of customers together with a layover (more precisely the time window length is set to the average maximum driving time, plus the layover duration, plus 10 times the setup time);

(2) satisfy the minimal time that has to take place between two consecutive time windows.

For each time window, considered in chronological order, and all trailers compatible with the driver, the following steps are carried out:

1. identify critical customers, i.e., a set of customers compatible with the trailer that must be visited in the time window in order to avoid stockouts;

2. estimate delivery quantities for these customers based on their consumption; 3. generate a set of td-shifts that cover all the selected customers, by means of a bestinsertion-like procedure as follows:

• first, the customer is tried to be inserted into a trip of an existing shift;

• if this is not possible, a new trip with a visit to the source and then to the customer is created;

4. update inventory levels accordingly. As a consequence, the set of critical customers may change;

5. for each td-shift the corresponding template is obtained by removing the association to the driver and to the trailer and by unfixing the starting time. The template is then inserted into CAT .

In order to generate several templates with the same procedure, the following steps are randomized:

• for each time window, the order of the trailers compatible with the driver is shuffled;

• the order of the customers the best-insertion-like procedure is shuffled;

• in the selection of the customers, it is possible to randomly add customers that do not need to be visited to avoid stockouts;

• when the set of shifts has been generated, it is possible to try to add other customers using the best-insertion-like procedure.

Initialization based on customer requirements

The procedure described here considers only the customers that must be delivered during the horizon in order to avoid stockouts. For each of them, we consider a deliverytime window. It starts when the minimal quantity can be delivered and ends when the customer runs out of stock. The customers are then sorted based on the closing date of the delivery time windows.

The best-insertion-like procedure is applied on this set of customers in order to obtain a set of td-shifts. Among all td-shifts that have been produced by the procedure, only one is kept. Delivery time windows are then re-computed for the customers served by this td-shift. If some customers do not require additional deliveries, they are not anymore considered. The procedure is repeated until sufficient quantities have been delivered to all customers to avoid stockouts over the time horizon. For each of the td-shifts that are obtained the corresponding template is inserted into CAT .

Insertion of new templates with the impicit-relaxation-based dynamic programming procedure

In order to populate CAT with new templates, the dynamic programming pricing procedure is called. This procedure depends on a few parameters that define the subset of td-shifts that is explored: a driver d 0 , a trailer t 0 , a starting time m 0 and a maximum ending time m 1 . To diversify the templates generated, the procedure is applied for all the feasible trailer-driver pairs. For each pair, several values of m 0 and m 1 are considered in order to cover the horizon [0, H 2 ]. We consider values of m 0 and m 1 such that m 1 -m 0 is large enough to make a layover shift. As a consequence, we set m 1 -m 0 to twice the maximum driving time, plus the layover duration, plus 10 times the setup time to allow enough time for the operation. Moreover, we take into account the td-shifts that have already been fixed. For each trailer-driver pair, we compute the availability time windows, i.e., where they do not operate fixed td-shifts. Then, we ensure that m 0 and m 1 values are inside these time windows. Algorithm 2 details the search for negative reduced-cost shifts using the implicit-relaxation-based dynamic programming procedure.

2 Generate shifts with negative RC 1: S ← ∅ 2: for all d ∈ D do end for 18: end for 19: return S

Pricing of the columns in catalog CAT

For each template r stored in CAT , we evaluate the potential benefit of including it into the model, namely, to create a td-shift r based on template r. We proceed as follows. A shift r is obtained by postponing by an integer number of periods over [0, H 2 ], the template r. The scheduling of the operations is not modified. The procedure considers all time-shiftings that satisfy time constraints, i.e., for every customer the delivery occurs during one of its time window).

Then, for each of these shifts r, compatible trailer-driver pairs are determined in order to generate a td-shift r. To this aim, we take into account the current partial solution. When a variable x r dt is fixed to 1, it means that trailer t is no longer available between st r and et r . Similarly, driver d is no longer available between st r -rest d and et r + rest d . Hence, driver d and trailer t are compatible with shift r if:

• d and t are compatible (t ∈ T d );

• t can operate shift r;

• d can operate shift r; • t is available during [st r, et r]; • d is available during [st r, et r].
Then, for each feasible td-shift r derived from template r, the corresponding exact reduced cost is computed using Expression (27). Only the td-shift with the lowest reduced cost is considered for inclusion into RM P . Note that certain terms of Expression ( 27) depend only on the sequence of the customers and need to be computed only once.

A parameter P eval impr is used to determine how much improvement in the reduced cost is needed to consider updating the best td-shift generated from template r (the td-shift considered for inclusion into RM P ). Since starting times are ordered, this avoid to postpone the starting time of a shift when the impact on the reduced cost is limited.

Selection of td-shifts to enter RM P

To construct shifts from templates in CAT , initially a template is randomly selected. The reduced cost of associated td-shifts is calculated as explained in Section 6.3.3. If the resulting reduced cost is lower than a negative value P eval rc , the td-shift is memorized and will be considered to enter RM P . The procedure terminates when P cat nbT emplates have been memorized or CAT has been completely scanned. Finally, the P cat nbAddBySearch td-shifts with the lowest reduced cost are inserted into RM P .

Variable-fixing heuristic

The solution of the current RM P is by nature fractional. In order to obtain a solution for the REC-IRP, we need to retrieve an integer solution with respect to variables x rdt . To this aim, we propose a variable-fixing heuristic. The purpose of this procedure is to select the td-shifts included in the REC-IRP solution by fixing the corresponding x rdt variables to 1. This variable-fixing heuristic solves a mixed integer programming problem M IX-RM P derived from the current RM P . To do so, all the x rdt associated with td-shifts that start before H 1 are considered as binary variables. The other variables remains continuous. After solving M IX-RM P , binary variables x rdt with value 1 are fixed to 1 in RM . Note that the mixed integer program is solved with a time limit of P M IP tl seconds, and a relative gap of P M IP gap . It is noteworthy that each time the variable-fixing heuristic is called, the previous fixing decisions can be updated. Hence a variable previously fixed can be unfixed. Moreover, the reduction of the size of RM P presented in Section 6.6 permits to remove the x rdt with value 0. This allows to maintain a reasonable size of x rdt variables over the horizon [0, H 1 ], hence to maintain a reasonable number of integer variables when solving the mixed integer program.

After fixing the x rdt variables over the horizon [0, H 1 ], it could happen, even after applying the recovery procedure presented in Section 6.5, that RM P does not admit a solution without non-zero artificial variables. In this case we detect an incorrect fixing strategy and another mixed integer program M IX-RM P obtained is solved. M IX-RM P is obtained from RM P by considering as binary variables all variables in X f ixed , where X f ixed contains all the x rdt variables currently fixed. Then, the two following constraints are added:

x∈X f ixed x ≤ |X f ixed | -1;
(29)

x∈X f ixed x ≥ |X f ixed | -P M IP uf . (30)
Constraint (29) forces to unfix at least one variable while Constraint (30) limits the number of unfixed variables. Note that if RM P turns to be still unfeasible, the horizon is increased by incrementing the values of H 1 and H 2 .

Recovery procedure

It may occur that the current RM P only admits a feasible solution where at least one artificial variable is non-zero. That means that the td-shifts together with the variables that have been fixed (see Section 6.4) in the current RM P cannot provide a feasible solution for the REC-IRP. In this case, we propose the following recovery procedure.

The procedure detects the source of by generating the set U of customers associated with a non-zero artificial variable (we recall from Section 3.6 that all the artificial variables are indexed by a customer).

Then CAT is first scanned focusing only on templates that contain at least a customer in U. The associated P cat nbAddBySearch td-shifts with the most negative reduced cost are then added to RM P .

A second set of P cat nbAddBySearch td-shifts is added to RM P using the procedure detailed in Section 6.3.4.

Those two searches into CAT are repeated until RM P is feasible or it has been executed P recover nbIter times. In case RM P is still unfeasible after the P recover nbIter iterations:

• CAT is repopulated calling the implicit-relaxation-based dynamic programming procedure;

• the whole procedure is applied using the updated catalog.

If RM P is still unfeasible after P recover nbU pdateCat updates of the CAT , at most P M IP uf

x rdt variables are unfixed as described in Section 6.4.

Reduction of the size of the RMP

Each time the current values of H 1 and H 2 are incremented we check if it is possible to reduce the size of RM P . The procedure simply checks if the current solution of RM P admits more than P reduc maxN ul variables x rdt set to zero or more than P reduc maxU nf ix unfixed variables. In this case, the procedure removes from RM P all variables x rdt set to zero. For each variable x rdt removed from RM P the associated variables q irt and the constraints related to td-shift r are also removed.

Solution improvement

Once a feasible solution has been found, we try to improve the current best solution. We create a set S containing each td-shift included in the current solution. Then, for each td-shift in S, we add to S additional td-shifts built as follows:

(1) customers receiving less than 1% of the total quantity delivered in the td-shifts are removed;

(2) visits, for which only the minimal quantities are delivered, are removed from the td-shift;

(3) the possible initial stop at the source is removed; than 200 customers with a planning horizon of 10 days. Namely these are instances 13, 14, 15, 16, 19, 24, 25, 26 and X2. The remaining instances are large instances, namely 12, 17, 18, 20, 21, 22, 23, X1, X3, X4 and X5. The CPU time limits were 3 hours for mid-size instances, 6 hours for large instances with less than 20 periods or with less than 10 drivers, and 12 hours for instances with more than 20 periods and more than 10 drivers. In Tables 12 and 13, we report the results obtained with our algorithm. The column Instance contains the name of the instance, while the column Ratio contains the corresponding value of the logistic ratio (multiplied by 100). Since tests were performed with two different values of P rh + , the column P rh + gives the value of this parameter that allows to obtain the corresponding ratio.

The branch-cut-and-price heuristic provides solutions on all small and medium-size instances, and on 9 out of 11 large instances. The solution quality is reasonable for most instances with respect to the best known solution values. However, this is at the expense of large computation times.

Conclusions and perspectives

This paper deals with the problem proposed in the context of the ROADEF/EURO challenge. The problem is a rich variant of the so-called Inventory Routing Problem encountered by Air-Liquide.

Even if the designing of exact solution approaches would not be suitable for large size instances as those proposed by the organizers of the ROADEF/EURO challenge, we proposed a branch-cut-and-price framework to tackle the REC-IRP. To deal with real-life instances the cut-and-column generation procedure is developed along with a heuristic pricing algorithm. To produce new columns, the pricing procedure relies on a surrogate relaxation of the constraints. Moreover, a heuristic fixing procedure has been integrated into the framework to generate feasible solutions.

The proposed method provides solutions of reasonable quality on all but two instances. However, this is at the expense of large computation times.

  St(R d ) set of starting times of shifts in R d K d m set of shifts incompatible at time m for driver d St(R t ) set of starting times of shifts in R t K t m set of shifts incompatible at time m for trailer t Table

  dt c rdt x rdt r∈R t∈Tr i∈Vr ( h∈H δ ihr c h )q irt

δ

  ihr γ th + γ trs -h∈H δ ihr σ ih -Z ir .

  t ∈ T d do 4: maxDuration ← Twice the maximum driving time, plus the layover duration, plus 10 setup times 5: T W ← time windows where d and t are available (based on the current td-shifts fixed into the RMP) min {m 0 + maxDuration, end(tw)} 9: while m 0 < H 2 and m 0 < m 1 do 10: S DY N ← td-shifts from dynamic programming (m 0 , m 1 , d, t) 11: S ← S ∪ S DY N 12: end ← max (r,d,t)∈S DY N {et r } 13: m 0 ← end + gap(d) 14:m 1 ← min {m 0 + maxDuration; end(tw)}

Table 1 :

 1 Notation (sets)

			of drivers
		T	set of trailers
	T H last period in set H
	rest d	minimal resting time between two shifts for driver d ∈ D
	Q t	vehicle capacity for trailer t ∈ T
	q min i	minimal quantity delivered to customer

d set of trailers that can be driven by driver d ∈ D D t set of drivers that can drive trailer t ∈ T T W d set of allowed time windows for driver d ∈ D

Table 2 :

 2 Notation (additional data)

Table 4 .

 4

	R d	subset of shifts compatible with driver d
	R t	subset of shifts compatible with trailer t
	c rdt	cost of shift r when operated by driver d and trailer t
	st r	starting time of shift r (in minutes)
	et r	ending time of shift r (in minutes)
	V r	subset of customers visited in shift r
	st ri	starting time of service for customer i in shift r (in hours)
	δ ihr	1 if i ∈ r and st ri = h, 0 otherwise
	T rip r set of trips in shift r
	V s	subset of customers visited in trip s

Table 3 :

 3 Notation introduced for shifts (part 1) R dt subset of shifts compatible with driver d and trailer t D r subset of drivers compatible with shift r T r subset of trailers compatible with shift r D rt subset of drivers compatible with shift r and trailer t T rd subset of trailers compatible with shift r and driver d

Table 4 :

 4 Notation introduced for shifts (part 2)

  ) we introduce notation St(R d ) and K d m . Given a time instant (minute) m and a driver d, K d m = {r ∈ R d : st r ≤ m, et r + rest d ≥ m}. K d m represents a set of shifts that are not compatible at time m either because they overlap or because they would not satisfy the minimal duration rest d imposed between two shifts assigned to the same driver. Constraints (2) impose these incompatibilities. They are only needed at time instants (minutes) where a new shift comes into play, i.e., at the starting times of the shifts. St(R d ) is the set of these starting times: St

  .

	dc d	working cost (per working time unit) for driver d ∈ D
	tc t	distance cost (per distance unit) for trailer t ∈ T
	s i	setup time for location i (0 if i is the base)
	t ij	travel time from location i to location j
	d ij	distance from location i to location j
	t layover d c layover d	duration of a layover for driver d ∈ D cost of a layover for driver d ∈ D
	T d	maximum driving duration for driver d ∈ D
	tw i	number of time windows for customer

Table 8 :

 8 that is, k indicates the earliest time window allowing to reach j from L). Inserting a layover in this case is imposed by the REC-IRP definition. The two

	vertex	last visited vertex
	cost	evaluation of the reduced cost
	time	starting time of the service at ending vertex if it is a source or a
		customer, arrival time otherwise (base), expressed in minutes
	driving time	remaining allowed driving time, expressed in minutes
	minimal quantity minimal quantity to be delivered in the current trip, according to
		values q min i
	gamma	current value for dual variable γ t 0 rs , accordingly with the vertices
		already visited in the current trip
	sources	number of visits to a source
	trips	number of trips (started or finished)
	layover	boolean value indicating if a layover has been planned
	layover customer boolean indicating if a layover customer has been visited
	visited[u]	boolean indicating if customer u ∈ V has been visited or not

Label definition other cases result from a simplification. When j is a layover customer and L contains no layover (lay = f alse), we decided to insert a layover before (second case) or after (third case) j. First, the earliest time window k = min{l : 1 ≤ l ≤ tw j and t + s i + t ij ≤ b jl } in j is computed. If the layover can be inserted before the closing time of the time window (t + s i + t ij + t layover d 0

1

  Increasing horizon heuristic 1: Initialize catalog CAT 2: H 1 ← P rh short , H 2 ← P rh long 3: Initialize RM P with parameters (CAT , H 1 , H 2 ) 4: while H 1 ≤ H do

	5:

Table 12 :

 12 Best known Instance |V| |H| (in |D| |T | P rh Results on mid-size instances.

	+	Ratio (×100) ratio (×100)

Table 13 :

 13 known Instance |V| |H| (in days) |D| |T | P rh Results on large instances.

	+	Ratio (×100) ratio (×100)

note that these two intervals are not expressed with the same time unit: one has to understand that the two time intervals overlap

(4) a final stop at the source is added (before going back to the base);

(5) new visits are considered by applying a best-insertion like procedure (this step is applied as well to td-shifts generated in previous points).

When S contains more than P improve maxSize td-shifts, the process is stopped. The formulation proposed in Section 3 restricted to S is then solved in the hope of improving the solution quality. In case of improvement, the procedure is repeated.

Computational experiments

The algorithm presented in Section 6 was run on a Intel(R) Core(TM) i7-5600U CPU @ 2.60GHz. Cplex 12.6 was used to solve the (mix-integer) linear programs considered in the procedure and described in the previous sections. The parameters used in the algorithm are set as presented in Table 11.

A maximum time is allocated to solve each instance. We distinguish two sets of instances based on their size. Mid-size instances include less than 100 customers, or less


Achraf W Aguaf

Romain Alvernhe

Ludivine Fadel

Marjorie Grzeskowiak

Energy Harvesting with 2.45 GHz Rectenna for urban application

Keywords: WI-FI, RF energy harvesting, Rectenna, Wireless Sensor Networks

This paper describes the conception of a 2.45 GHz rectenna for energy harvesting application. Electromagnetic simulations have been carried out using Computer Simulation Technology software. The enhancement of the simulated output DC voltage is obtained with resonant circuit and Cockroft Walton boost (containing four zero bias diodes) for a sinusoidal input voltage. The simulated proposed rectenna with single zerobias diode and RF-DC boosting circuit have been realized and measured inside an anechoic chamber. The behavior of measured plots is in accordance with the simulated ones. The simulated RF-DC boosting rectenna enhanced the output DC voltage up to 140 mV for 1 µW/cm 2 power density.

I. INTRODUCTION

In recent years, Wireless Sensor Networks (WSNs) push harvester development to provide almost infinite lifetime to sensor from environmental energy. Further-more, the use of wireless electronic devices has become relevant in many fields (military, medical) and especially in unsafe places where common power supply remains restrictive. Among different renewable ambient power sources such as solar, vibration, or electromagnetic waves, radio frequency (RF) energy presents ubiquitous availability with low power density (around 1 µW/cm 2 at 2.45 GHz indoors) in comparison with the energy sources (up to 15 mW/cm 2 outdoors and 10100 µW/cm 2 indoors for solar energy) [START_REF] Kim | Georgiadis Ambient RF Energy-Harvesting Technologies for Self-Sustainable Standalone[END_REF], [START_REF] Anfr | Etude de lexposition du public aux ondes radiolectriques[END_REF].

This work addresses the issue of energy harvesting with the development of COTS based RF harvesters for urban areas with the ultimate goal of using the large transparent surfaces of the urban environment as a support for implementation. The RF energy harvester, commonly named rectenna, is dedicated to transform electromagnetic waves into an electrical signal by combining an antenna and an RF-DC conversion circuit, [START_REF] Takhedmit | Transparent Rectenna and Rectenna Array for RF[END_REF]. The rectifier is usually realized with Schottky diodes because of their low threshold voltage in RF harvesting applications.

While a RF-DC rectifier provides 0.3 V with 20 dBm [START_REF] Gudan | Ultralow power 2.4GHz RF energy harvesting and storage system with -25dBm sensitivity[END_REF], the voltage requirement for transistor switching operation is at least 0.5 V. To overcome low rectified voltage, DC-DC converter or RF-DC boost have been studied [START_REF] Ermeey | Indoor 2.45 GHz Wi-Fi Energy Harvester with Bridgeless Converter[END_REF][START_REF] Adami | A Flexible 2.45-GHz Power Harvesting Wristband With Net System Output From 24.3 dBm of RF Power[END_REF][START_REF] Kanaya | Energy Harvesting Circuit on a One-Sided Directional Flexible Antenna[END_REF][START_REF] Masotti | Genetic-based design of a tetra-band high-efficiency radio-frequency energy harvesting system[END_REF][START_REF] Vera | Design of a 2.45GHz rectenna for electromagnetic (EM) energy scavenging[END_REF][START_REF] Olgun | Investigation of rectenna array configurations for enhanced RF power harvesting[END_REF][START_REF] Bergs | A dual Band 915MHz/2.44GHz RF Energy Harvester[END_REF]. Finally, to improve the input voltage of the rectifier, and so increase the efficiency of the RF-DC conversion, a resonator is added between the antenna and the rectifier [START_REF] Kanaya | Energy Harvesting Circuit on a One-Sided Directional Flexible Antenna[END_REF]. This paper focuses on the design and implementation of a rectenna in the unlicensed ISM Bands at 2.45 GHz on FR4 printed circuit boards (PCB) with off-the-shelf components. Each part of the schematic bloc is studied; in section 2, with electromagnetic software for the antenna de-sign and with electrical software for boost structure in harmonic balance to take into account the non-linearity of the diodes. In section 3, co-simulation of the antenna, matching network, resonant circuit and Cockroft Walton circuit have been carried out under CST (Computer Simulation Technology) and compared with measurement in anechoic room.

II. RECTENNA BUILDING BLOCKS

The basic architecture of an RF energy harvester is presented in Fig. 1. It is composed of antenna to receive the electric field, L-matching circuit to optimize impedance matching between the antenna and the boosting circuit, a resonant circuit to improve the RF signal and a RF-DC boost to convert the RF signal in DC output voltage. The power density of the electromagnetic source varies from 1 to 30 W/cm 2 corresponding to ambient (E ≈1.94 V/m) and transferred energy (E ≈10.63 V/m).

Fig. 1: Block diagram of rectenna

The different blocks have been designed and realised on 1.6 mm thickness FR4 substrate to obtain a cheap prototype by microetching.

A. Antenna

The first bloc of the rectenna is a strip-loop antenna designed on a 1.6mm thickness FR4 substrate [START_REF] Takhedmit | Transparent Rectenna and Rectenna Array for RF[END_REF] whose electrical characteristics are ( r = 4.4 and tan δ = 0.025). The design of the antenna is illustrated in Fig. 2. A capacitor of 100 pF is placed in the gap on the side W 1 . The antenna impedance has been simulated in the measurement plan mentioned in Fig. 2 and reported in Fig. 3 with the value of Zin = 343 Ω at 2.72 GHz. The proposed rectenna have been realized by mi-croetching and measured inside an anechoic chamber. The measurement setup contains a 9 dB transmitting antenna, RF generator (Agilent 864D) and a gain power amplifier (Amplifier research 15S1G3), Fig. 9. The final dimensions of the rectenna are 7×7.6 cm 2 .

Fig. 9: Rectenna measurement setup

The distance between the emitting and receiving antennas is 2 m to be in far field conditions. A voltage meter (Voltcraft VC150) is connected in parallel with the output to measure DC output voltage. A RF probe (RF Survey Meter, EMR-300 Broadband) is placed next to the receiving antenna to measure the power density. In simulation a wave plane linearly polarized is applied on the rectenna (antenna + circuit). DC output voltage is obtained on the resistive load. On Fig. 10, DC output voltage is reported versus frequency from 1 to 4 GHz. The behavior of the measured and simulated plots presents three voltage peaks. The discrepancies between simulation and measurement are twofold. First, the measured working frequency is lower than the simulated ones. Secondly, the level of DC voltage is weaker experimentally. DC output voltage is shown when the power density varies from 0.5 µW/cm 2 to 21 µW/cm 2 . While the measured DC output voltage is lower for some power density values in comparison with the simulation, the evolution versus power density is such as expected. Each bloc of the rectenna is measured experimentally and confirms the matching between the antenna with LC matching circuit and the Cockroft Walton associated to the resonant circuit. A frequency shift is observed for the antenna resonance while the reflection parameter is less than 15 dB from 1 to 4 GHz for the RF-DC boosting circuit: the antenna with LC matching circuit is matched at 1.81 GHz (S11 = -15 dB). The discrepancy between measured and simulated results is attributed to the weak measured power efficiency of the RF-DC boosting circuit. The RF-DC boosting circuit has to be realized again and the system shifted to the Wifi working frequency, 2.45 GHz.

IV. CONCLUSION

A compact RF-DC boosting rectenna has been developed with components (capacitors, inductors and diodes) inside the receiving part of the structure. Co-simulation and measurement have been compared: DC output level is observed versus frequency and power density, while RF-DC boost has been highlighted in transient simulation. The rectenna prototype has to be improved experimentally, especially the power efficiency of the RF-DC boosting circuit: to discuss the RF-DC boosting circuit, the power conversion efficiency η (%) and DC voltage -RF power transfer function (mV/mW) could be reported for variable input power. In the following study, the FR4 dielectric could be replaced by Plexiglas substrate, not only to reduce the dielectric losses, but also to realize prototypes on transparent substrates, such as building glasses with invisible conductors. Others stages of voltage doublers can be added to obtain the required DC voltage, or array of rectennas to improve the RF power.

Fig. 2 :

 2 Fig. 2: Antenna dimensions in mm : L 1 = 37, L 2 = 100, W 1 = 49, W 2 = 100, W 3 = 2, S = 1

Fig. 3 :

 3 Fig. 3: Simulated antenna impedance

Fig. 4 :

 4 Fig. 4: L-matching circuit

Fig. 5 :

 5 Fig. 5: Cockroft Walton voltage multiplier circuit

Fig. 6 :Fig. 7 :

 67 Fig. 6: Simulated S11 parameter of resonant circuit

Fig. 8 :

 8 Fig. 8: Transient simulation ADS

Fig. 10 :

 10 Fig. 10: Simulated and measured DC-output voltage

Fig. 11 :

 11 Fig. 11: Simulated and measured DC-output voltage

ACKNOWLEDGMENT

The circuits have been designed thanks to Computer Simulation Technology software (CST) that is available in Electrical Engineering School -ENSEIRB-MATMECA.

The authors would like to thanks CST for their soft-ware licence, Jean-Marc Micouleau, David Bedenes and Jerome Le Gall, ENSEIRB-MATMECA INP Bordeaux, France for the circuits realization and their technical assistance, especially with anechoic chamber