Daniel Chemla

AND Frédéric Meunier

Roberto Wolfler Calvo

BIKE SHARING SYSTEM: SOLVING THE STATIC REBALANCING PROBLEM

This paper deals with a new problem that is a generalization of the many to many pickup and delivery problem and which is motivated by operating self-service bike sharing systems. There is only one commodity, initially distributed among the vertices of a graph, and a capacitated single vehicle aims to redistribute the commodity in order to reach a target distribution. Each vertex can be visited several times and also can be used as a buer in which the commodity is stored for a later visit. This problem is NP-hard, since it contains several NP-hard problems as special case (the TSP being maybe the most obvious one). Even nding a tractable exact formulation remains problematic.

This paper presents ecient algorithms for solving instances of reasonable size, and contains several theoretical results related to these algorithms. A branch-and-cut algorithm is proposed for solving a relaxation of the problem. An upper bound of the optimal solution of the problem is obtained by a tabu search, which is based on some theoretical properties of the solution, once xed the sequence of the visited vertices. The possibility of using the information provided by the relaxation receives a special attention, both from a theoretical and practical point of view. It is proved that to build a feasible solution of the problem by using the one obtained by the relaxation is an NP-hard problem. Nevertheless, a tabu search initialized with the optimal solution of the relaxation often shows that it is the optimal one.

The algorithms have been tested on a set of instances coming from the literature, proving their eectiveness.

Introduction

Nowadays self-service bike sharing systems are ourishing all over the world. The exploitation of such transport systems implies various problems and one of them is to ensure people that they will be able to nd a bike or to park it at each station all the day long. Therefore, a regulation system is necessary for maintaining a prexed optimal number of bikes at each station in order to fulll at best the demand. It is done with a eet of vehicles which are driven around the city and move bikes from a place to another. This work restrains to the rebalancing problem with one vehicle and in the static case. Static means that users cannot act on the bikes during the rebalancing process. The static rebalancing problem can be classied as a Single Vehicle One-commodity Capacitated Pickup and Delivery Problem (SVOCPDP). This problem is the one faced by an operator during the night when the number of moving bikes is negligible and when the city is divided into districts. Each district is covered by a single truck that has to redistribute the bikes in order to respond to the morning peak at best. This is a realistic problem since it is during the night that the impact of the regulation is the most important, as explained by operators. There are even such systems (for instance the one in Lyon, France), which are closed during the night, to make this regulation task easier.

To get an idea of the size of such a system, the numerical features of the Vélib system in Paris are presented. This system oers more than 20 000 bikes deployed in about 1 200 stations that are in Paris and its border cities and a eet of twenty three trucks of capacity equal to 20 are used to move bikes during the day to match the demand. If the city is divided into areas on which only one truck operates, then each truck has to work on about 50 stations. This paper focus on this part of this operating problem: how to deal with a part of the city assigned to a single vehicle. The multiple vehicles case is currently under study (see [START_REF] Chemla | Balancing a bike-sharing system with multiple vehicles[END_REF]).

The problem can be formalized as follows. Let G = (V, A) be a complete directed graph where V = {0, . . . , n} is the vertex set composed by n + 1 vertices, the vertices in {1, . . . , n} representing the stations and the vertex 0 representing the depot, and where A is the set of arcs. For each arc (i, j) ∈ A, we denote by c (i,j) the cost of the arc (i, j). The cost is assumed to satisfy the triangular inequality (i.e., c (i,j) +c (j,k) ≥ c (i,k) for all i, j, k ∈ V ). Each vertex i has a capacity C i ∈ Z + . For each vertex i ∈ V , its initial state in bikes is dened by p i ∈ Z + and its target, or nal, state by q i ∈ Z + . A vertex is in excess (resp. in default ) if p i > q i (resp. p i < q i ). Some vertex can be initially balanced i.e. p i = q i . Moreover, throughout the paper the imbalance d i = p i -q i is used and the depot is always assumed to have no bike: C 0 = p 0 = q 0 = 0 and d 0 = 0. The vehicle has also a capacity Q.

A feasible solution for the SVOCPDP, also called a route, is a sequence of vertices, starting and nishing with the depot 0, together with bike displacements within the limits of capacity constraints, at the end of which the system is balanced: each vertex i has been brought from its initial state p i to its target state q i . In this case, the sequence of vertices is said to be induced by the route. The goal of the SVOCPDP is to nd the minimal cost route. Note that the convergence for each vertex from p i to q i is not required to be monotonous: drops and multiple visits of the vertices are allowed, i.e. bikes can be loaded from vertices in default or unloaded at vertices in excess and transfers can take place at initially balanced vertices. Figure 1 shows an example of instance with 9 vertices (the depot and 8 stations) including one initially balanced vertex. A pair of values representing (p i , q i ) is displayed next each vertex and the capacity of the vehicle is equal to 8. Complexity. The SVOCPDP is NP-hard since it contains NP-hard problems as special cases. It is easy to see it, but details are given for sake of completeness. The TSP is obviously one of them. Set p i = 0 and q i = 1 for all vertices i ∈ {1, . . . , n -1}, and set p n = n -1 and q n = 0. Add a depot at distance 0 from the vertex n. Set Q = n -2. The optimal solution of the SVOCPDP coincides with the optimal TSP solution.

The 2-partition is another special case. Let b 1 , . . . , b n be n non-negative integers such that i b i is even. Dene m = 1 2 n i=1 b i . The 2-partition problem is a decision problem asking whether there exists I ⊆ {1, . . . , n} such that i∈I b i = m. Take the complete graph with n + 3 vertices: n vertices with p i = b i and q i = 0; 2 vertices with p i = 0 and q i = m and the depot. Dene the c (i,j) to be 1 for each arc (i, j) and the capacity of the vehicle equal to m. The optimal solution of the SVOCPDP problem is equal to n + 3 if and only if there is a subset I ⊆ {1, . . . , n} such that i∈I b i = m.

Moreover, an optimal route of the SVOCPDP problem may not have an encoding that is polynomial in the size of the input. The simple case with three vertices the depot 0, one vertex 1 with p 1 = B and q 1 = 0 and one vertex 2 with p 2 = 0 and q 2 = B is enlightening. Assuming that Q = 1, the optimal sequence of vertices is 0, 1, 2, 1, 2, . . . , 0, with a number of terms = 2B + 2, although the input is in O(log 2 B).

1.2. Main contributions of this paper. This paper deals with a new problem (see [BBC + 11] for a rst reference on this problem with theoretical results such as approximation algorithms or special polynomial cases). We emphasize that we are not only interested by ecient algorithms for solving it, but also by theoretical results linked with these algorithms and giving a better understanding of the problem.

An exact model is presented, which does not seem to be tractable. Therefore, a relaxation, which turns out to be an integer program with an exponential number of constraints, is proposed and solved by a branchand-cut algorithm. The optimal solution of this relaxation generally provides a good lower bound of the optimal solution of the original problem (the quality of this lower bound is proven experimentally). For a xed sequence of vertices visited by the vehicle, we prove that it is possible to nd in polynomial time the operations that bring the system to the possible state nearest to the target state (Propositions 1 and 2).

In particular, we prove that it is possible to decide in polynomial time whether there is a feasible solution inducing precisely a xed sequence of visited vertices. It allows to derive a purely combinatorial tabu search providing upper bounds for the problem: the underlying neighborhood does not take care of the loading of the vehicle. The tabu search proves good behaviors. Moreover, a special attention is given to the way of using the information provided by the relaxation: both from a theoretical point of view and from a practical one. From a theoretical point of view, we are able to prove that even if the relaxation is quite natural and provides (experimentally) good lower bounds, it is an NP-complete problem to decide whether there is an optimal solution of same cost. More precisely, given the number of times each arc is traversed, and given the number of bikes carried on each visit on each arc, it is an NP-complete problem to decide whether there is an optimal solution realizing these numbers. This property is formalized in Propositions 58. From a practical point of view, it is experimentally shown that the solutions obtained with the relaxation are often the optimal ones. Indeed, the tabu search is in general able to nd a solution of almost the same cost, when initialized with the optimal solution of the relaxation. 1.3. Notations and basic notions. We dene for all subsets S ⊆ V :

• S = V \S • δ + (S) := {(i, j) ∈ A : i ∈ S; j ∈ S} • δ -(S) := {(i, j) ∈ A : i ∈ S; j ∈ S} • δ(0) := δ + ({0}) ∪ δ -({0}) • d(S) = j∈S d j • µ(S)
is equal to 1 whenever there is at least one initially non-balanced vertex in S, 0 otherwise.

Given z ∈ R A + , G[z]
is the directed graph obtained from G by deleting the arcs a with z a = 0. This graph is called the support graph of z.

A notion used several times in the paper is the one of a b-ow. A b-ow is an usual notion in combinatorial optimization (see for instance [Sch03, KV02]). Given a directed graph D = (U, A ), a value b ∈ R U , and

capacities l, u ∈ R A with l ≤ u, a b-ow is a map f : A → R such that l a ≤ f (a) ≤ u a for all a ∈ A and a∈δ + (v) f (a) = b v + a∈δ -(v) f (a) for all v ∈ U . If it
exists, a b-ow can be computed in strongly polynomial time. Moreover, when all the b v and the l a , u a are integral, if a b-ow exists, there is an integral one.

1.4. Plan. In Section 2, dierences between the SVOCPDP and problems in the literature are presented.

Section 3 presents the aforementioned proposition that enables to nd in polynomial time the operations that bring the system to the possible state nearest to the target state, given a xed sequence of vertices visited by the vehicle. An exact model of the problem is given in Section 4, and a relaxation is presented in Section 5. In Section 6, we prove that deciding whether a solution of the relaxation is a feasible solution for the SVOCPDP is NP-complete. The Section 7 contains the description of the branch-and-cut algorithm that is used to solve the relaxation. The proposition of Section 3 is used in Section 8 for deriving a tabu search. Finally, Section 9 presents the computational results on instances from the literature with a slight adaptation in order to t the SVOCPDP's features.

Literature review

In the literature, similar problems are the One-Commodity Pickup-and-Delivery TSP (1PDTSP) studied by [HPSG04a] and the Swapping Problem dened by [START_REF] Anily | The swapping problem[END_REF]. This latter has been solved by [START_REF] Bordenave | A branch-and-cut algorithm for the preemptive swapping problem[END_REF][START_REF] Bordenave | A branch-and-cut algorithm for the non-preemptive swapping problem[END_REF] in its preemptive and its non-preemptive versions. For solving the 1PDTSP, several papers have been published recently in the literature. They propose dierent exact and heuristic algorithms (see [HPSG04b],

[HPSG07], [START_REF] Hernandez Pérez | A hybrid grasp/vnd heuristic for the one-commodity pickup-and-delivery traveling salesman problem[END_REF] and [START_REF] Hana | A general variable neighborhood search for the one-commodity pickup-and-delivery travelling salesman problem[END_REF]). Since these two problems present several similarities with our problem, more details about them are given in the next paragraph. A recent conference paper by [RTF09] is motivated by the same application. They discuss dierent variants of the problem of repositioning the bikes which are modeled by mixed linear programs and solved using CPLEX. For some variants, they are able to solve several instances to optimality (up to 60 stations and 2 vehicles for their so-called Arc-Indexed variant). Our model is close to their Sequence-Indexed variant (see Section 3.4 of [RTF09]), but instead of computing a minimum cost route with xed target states, they try to nd the best repartition of bikes that can be achieved by a eet of one or several vehicles within a time limit. Moreover, in the Sequence Index formulation given by [RTF09], drops are not allowed, a thing that will be xed in an upcoming version by

[RT11].
The SVOCPDP gathers aspects from both the Swapping Problem and the 1PDTSP but diers by main features. In the Swapping Problem a single vehicle has to move unitary objects from a vertex to another.

There are m dierent types of objects and the vehicle has a unitary capacity: it can only move one object at the time. [START_REF] Anily | The swapping problem[END_REF] showed that a vertex could be visited at most three times in the optimal solution. This theorem was the starting point of the work of [START_REF] Bordenave | A branch-and-cut algorithm for the preemptive swapping problem[END_REF][START_REF] Bordenave | A branch-and-cut algorithm for the non-preemptive swapping problem[END_REF] for solving the Swapping Problem with a branch-and-cut algorithm. The SVOCPDP is dierent since there is only one type of object (bikes), but the supply and the demand are greater than one and the vehicle capacity is Q. Anily and Hassin's theorem does not hold anymore: for example in the trivial instance with one pickup vertex with p i = 5Q and q i = 0 and one delivery vertex with p i = 0 and q i = 5Q, the vehicle has to do ve round trips between the two vertices.

The main dierence with the 1PDTSP is that in the SVOCPDP a vertex can be visited several times, whereas in the 1PDTSP the solution is required to be an Hamiltonian cycle. In the 1PDTSP, the depot is assumed to have an unlimited number of objects, but this fact can be modelled in the SVOCPDP by a vertex at distance 0 from the depot with Q bikes available and a target state of Q bikes. A solution of the 1PDTSP is a feasible solution for the SVOCPDP, but the SVOCPDP can have a lower-cost optimal solution. Moreover, the 1PDTSP may have instances without feasible solutions (the instance of Figure 3, with Q = 3, is an example), something that can not happen for the SVOCPDP.

The fact that one can get better solutions in routing problems when customers are allowed to be visited several times has already be noticed in other works. Note for instance the work by [START_REF] Archetti | A Tabu Search Algorithm for the Split Delivery Routing Problem[END_REF], in which a Split Delivery Problem is solved through a tabu search. In the SVOCPDP, any vertex can be used as a buer where bikes can be indierently temporary loaded or unloaded before being moved to their nal destination.

In particular, initially balanced vertices are not required to be visited, but may act as temporary depot in some optimal solutions. The buers can improve the optimal solution in some instances such as shown in the example given in Figure 3 and in Figure 5. In Figure 3 the square with the 0 inside is the depot and its distance with vertex 1 is null. For any other vertex i, (p i ; q i ) is given. The capacity of the vehicle is Q = 3 and the distances between all the peripheral vertices to the central one is 1. We show here the optimal solution, whose value is equal to 10. In this solution a bike is temporary unloaded at the central vertex 6 and then reloaded on the vehicle. The optimal solution would be equal to 12, if drops would be forbidden, because one of the odd number vertices would have to be visited twice.

(2; 0) 0 (2; 0) (2; 0) Figure 4 shows the best solution of the 1PDTSP. The two rst vertices with initial and target states (Q; 0) and (0; Q) (Q being the capacity of the vehicle) are here to prevent a use of the unlimited number of bikes in an optimal solution, available in the depot for the 1PDTSP. Figure 5 represents the optimal solution for the SVOCPDP. In both gures, the square with 0 is the depot, as in Figure 3. In Figure 5, the vehicle takes a bike from the left-corner vertex of the square, increasing momentary the decit in bikes, before coming back with the two missing bikes. With Euclidean distances, the solution of the SVOCPDP is better than the one of the 1PDTSP.

(0; 3) (0; 3) (0; 0)

0

(0; 1)

(3; 0) 3 0 0 (Q; 0) (0; Q) Q (2; 3) (2; 3) 1 2 0 Figure 4.
An optimal solution for the 1PDTSP

Dealing with sequences and routes

The diculty of the SVOCPDP is that a feasible solution is identied by both a sequence of vertices and a set of numbers of bikes carried on arcs. The following proposition (and its companions Propositions 2 0 (0; 1)

(3; 0) 1 0 3 2 0 0 (Q; 0) (0; Q) Q 0
(2; 3) (2; 3)

Figure 5. An optimal solution for the SVOCPDP with the same input as for Figure 4, the non-monotonous convergence helps and 3) enables us to work only with sequences of vertices. It will be particularly useful in Section 8 when we will design a local search for the SVOCPDP.

Proposition 1. Let i 1 , i 2 . . . , i k be a sequence of vertices, starting and nishing at the depot 0 = i 1 = i k .

There is a polynomial algorithm nding new initial and target states (p i , q i ) for each vertex i and a route inducing this sequence of vertices and satisfying these new states such that

• 0 ≤ p i ≤ p i and 0 ≤ q i ≤ q i for each vertex i

• i p i = i q i • the quantity i∈V p i is maximal.
In particular, it is possible to decide in polynomial time whether a sequence of vertices is induced by a feasible solution (in this later case, p i = p i for each vertex i). The proposition says roughly speaking that it is possible to nd the best bike displacements compatible with the sequence of vertices. The quantity i∈V (q iq i ) can be interpreted as a kind of degree of infeasibility of a sequence and the proposition shows how to minimize it polynomially.

Proof. Let us build a directed graph D = (U, A ) as follows. U has k + 2 elements: each vertex i j (make as many copies of a vertex i of G as there are occurrences of i in the sequence) and two more vertices s and t. The arcs in A are of four types:

(1) one arc between s and the rst occurrence of each vertex i in the sequence, with capacity p i , (2) one arc (i j , i j+1 ) for each j = 1, . . . , k -1, with capacity Q, (3) one arc between the rth occurrence of each vertex i and its (r + 1)th occurrence, if there is one, with capacity C i , (4) one arc between the last occurrence of each vertex i in the sequence and t, with capacity q i . See Figure 6 for an illustration with the sequence 0 → 1 → 2 → 3 → 1 → 3 → 1 → 2 → 0. Computing a maximum s-t ow in this graph leads to the proposition. Indeed, any s-t ow on D encodes possible bike displacements compatible with the given sequence of vertices. The numbers of bikes to be moved while going from i j to i j+1 are given by the ow on arc (i j , i j+1 ) (arcs dened in 2.); the number of bikes remaining in a vertex i after the rth visit of the vehicle is given by the ow on arc between the rth occurrence of vertex i and its (r + 1)th occurrence (arcs dened in 3.); the initial and nal numbers of bikes in a vertex i are given respectively by the ows on arcs dened in 1. and 4. And conversely, any bike displacements compatible with the given sequence of vertices induce an s-t ow.

p i is then the value of the ow in the arc between s and the rst occurence of i, and q i the value of the ow in the arc between the last occurence of i and t. If a vertex i is not present in the sequence, then we set p i = q i = min(p i , q i ).

In the case that all C i are suciently large, Proposition 1 has a nice and maybe quite unexpected corollary, formalized by the following proposition. If we are not allowed to change the initial states, but only the nal ones, the solution given by Proposition 1provides a route with new nal states closest to the original ones.

Proposition 2. Assume that we have C i = +∞ for all vertices i. Let i 1 , i 2 . . . , i k be a sequence of vertices, starting and nishing at the depot 0 = i 1 = i k . Consider the problem of nding bike displacements along this sequence leading to the nal state q closest to q when starting from the initial state p, where closest means the state that minimizes the L 1 norm ||q -q|| 1 = i∈V |q iq i |. The bike displacements given by Proposition 1 for this sequence are precisely the solution of this problem.

Proof. See Appendix.

As a by-product of the proof of Proposition 1, we get Proposition 3. Assume that we have a feasible solution to the SVOCPDP with fractional numbers of bikes carried during some moves. Then there is also a feasible solution of the same cost with only integral numbers of bikes carried during the moves.

Thanks to this proposition, we can forget the integrality constraint without changing the cost of the optimal solution.

An exact model

The exact model is based on the assumption that a constant β i is given for each vertex i, which is an upper bound of the number of times the vehicle has to visit i in any optimal solution. For instance, given a feasible solution of total cost O, we can set

β i := O/ min j =i c (i,j) ([Min10]
). The model is intractable with a direct approach, since the β i computed in this way can be quite huge and since the model involves variables with four indices. It would be interesting to compute a tighter β i , which would reduce the size of the exact model, but whether it is possible from a priori considerations remains an open question.

We introduce the following variables. x i,t takes the value 1 only if vertex i is visited at least t times. z i,t,i ,t takes the value 1 only if the vehicle visits vertex i for the t th time just after having visited i for the tth time, and y i,t,i ,t is then the number of bikes that is carried by the vehicle during this move. The variables y are not required to be integral since, according to Proposition 3, it does not change the optimal value of the linear program. u i,t,i ,t is a variable which takes the value 1 if the t th visit of vertex i comes in the route after the tth time of vertex i (but not necessarily right after).

(P ) z(P ) = min (i,i )∈A β i t=1 β i t =1 c (i,i ) z i,t,i ,t s.t. xi,t ≤ xi,t-1 ∀i ∈ V, ∀t ∈ {2, . . . , βi} (1)
i ∈V \{i}

β i t =1 z i,t,i ,t = xi,t ∀i ∈ V, ∀t ∈ {1, . . . , βi} except for (i, t) = (0, 2) (2) i ∈V \{i} β i t =1 z i ,t ,i,t = xi,t ∀i ∈ V, ∀t ∈ {1, . . . , βi} except for (i, t) = (0, 1) (3) u i,t,i ,t ≥ u i,t,i ,t + z i ,t ,i ,t -1 ∀i, i , i ∈ V, ∀t ∈ {1, . . . , βi}, ∀t ∈ {1, . . . , β i }, ∀t ∈ {1, . . . , β i } (4) u i,t,i ,t ≥ z i,t,i ,t ∀i, i ∈ V, ∀t ∈ {1, . . . , βi}, ∀t ∈ {1, . . . , β i } (5) ui,t,i,τ = 0 ∀i ∈ V, ∀t, τ ∈ {1, . . . , βi} with τ ≤ t (6) y i,t,i ,t ≤ Qz i,t,i ,t ∀i, i ∈ V, ∀t ∈ {1, . . . , βi}, ∀t ∈ {1, . . . , β i } (7) 0 ≤ i ∈V \{i} t τ =1 β i τ =1 (y i ,τ ,i,τ -y i,τ,i ,τ ) + pi ≤ Ci ∀i ∈ V, ∀t ∈ {1, . . . , βi} (8)
i ∈V \{i}

β i τ =1 β i τ =1 (y i ,τ ,i,τ -y i,τ,i ,τ ) + pi = qi ∀i ∈ V (9) z i,t,i ,t , u i,t,i ,t ∈ {0, 1} ∀i, i ∈ V, ∀t ∈ {1, . . . , βi}, ∀t ∈ {1, . . . , β i } (10) y i,t,i ,t ∈ R+ ∀i, i ∈ V, ∀t ∈ {1, . . . , βi}, ∀t ∈ {1, . . . , β i } (11) xi,t ∈ {0, 1} ∀i ∈ V, ∀t ∈ {1, . . . , βi} (12)
A way to see the exactness of this formulation consists in introducing a new graph, G β , obtained by making β i copies of each vertex i of G, and whose arcs link any pair of copies of distinct vertices. The idea is similar to the one for dynamic ows, for which one introduces usually time-expanded graph. A vertex of G β is of the form (i, t) with i ∈ V and t ∈ {1, . . . , β i } and any arc of the form ((i, t), (i , t )), with i = i . A vertex (i, t) in G β models the tth passage of the vehicle on vertex i. A feasible solution of the SVOCPDP is equivalent to an elementary directed path in G β , starting from (0, 1), nishing at (0, 2) and visiting (i, t) only if all (i, τ ) with τ < t have been visited.

(P) models this equivalent version on G β . Constraints (1) are logical constraints implied by the denition of x i,t . Constraints (2) and (3) ensure that we get in G β an elementary path starting at (0, 1) and nishing at (0, 2). At rst glance, there might be circuits as well. Constraints (4), (5), and (6) ensure that the solution is coherent with the t index: the tth visit of vertex i comes after the τ th visit, for any τ < t. As a direct by-product, we get also that there is actually no circuit. Constraints (7) limit the number of bikes transshipped on an arc at each move to the capacity of the vehicle. Constraints (8) bound the number of bikes parked at a vertex i between 0 and its maximal capacity C i at any time step. Constraints (10), ( 11) and (12) ensure that the vertices have all reached their target state at the end of the route.

Remark. The coherence on the t indices implied by the variables u i,t,i ,t and the constraints (4), ( 5) and ( 6) in (P ) can alternatively be obtained by the introduction of variables h i,t ∈ {1, . . . , n i=1 β i } encoding the instant of tth passage on vertex i and the use of big-M constraints, in a similar spirit as for the TSP with time windows (see for instance [START_REF] Desrochers | Improvements and extensions to the miller-tucker-zemlin subtour elimination constraints[END_REF]).

Relaxations

This section presents two equivalent mixed integer programs. They represent a relaxation of the original problem, since even solved to optimality they produce a lower bound of the optimal solution. 5.1. A rst relaxation. Using the variables of the exact model (P ) of Section 4, a relaxation for the SVOCPDP can be obtained by dening z (i,i ) = βi t=1 β i t =1 z i,t,i ,t , which represents the number of times the arc (i, i ) ∈ A is traversed in the solution, and by dening y (i,i ) = βi t=1 β i t =1 y i,t,i ,t , which represents the total number of bikes transported on it.

The relaxation is

(RP 1) z(RP 1) = min (i,j)∈A c (i,j) z (i,j) (13) s.t. j∈V z (i,j) = j∈V z (j,i) ∀i ∈ V (14) i∈V \{0} z (0,i) = 1 (15) i∈V \{0} y (0,i) = 0 (16) p i + j∈V \{i} y (j,i) = q i + j∈V \{i} y (i,j) ∀i ∈ V (17) (i,j)∈δ + (S) z (i,j) ≥ µ(S) S ⊆ V \ {0} (18) 0 ≤ y (i,j) ≤ Qz (i,j) ∀(i, j) ∈ A (19) z (i,j) ∈ Z + , ∀(i, j) ∈ A (20)
Again, requiring that the y (i,j) are integral does not improve the value of the relaxation (RP1). Indeed, whenever the values of the z (i,j) are xed, the variables y (i,j) are solutions of a b-ow problem, and hence, since the p i , q i and Q are integer numbers, there is an optimal solution with integral values for the y (i,j) . Constraints ( 14) and (15) come from the constraints (2), and (3). Constraints ( 16), ( 17) and ( 19) come from the constraints (7) and (9). Finally, constraints (18) ensure the connectivity of the solution, while initially balanced vertices might be skipped in an optimal solution. The proposed relaxation does not take into account the evolution on the number of bikes on each vertex at each time-step: all the moves are considered simultaneously and so the sequential dimension of the problem disappears. Therefore, although any solution of the SVOCPDP provides a solution of (RP1), a solution of model (RP1) might not be a solution of the SVOCPDP, as shown by the example given in Figure 7. In this example the values (p i ; q i ) are displayed next to each vertex, Q ≥ 2 and the optimal solution of (RP1) is represented: each arc a in the gure corresponds to z a = 1, the others z a being equal to 0 and the numbers near each arc are the y a . The optimal solution satises model (RP1), but violates SVOCPDP: the vehicle would have to take one bike from the rst vertex that is not yet arrived there.

However, the solution of (RP1) is often a solution to the SVOCPDP as illustrated in the computational results section (Section 9). 5.2. A second relaxation. The former mixed integer program requires two families of variables z (i,j) and y (i,j) . Linear program solvers are sensitive to the number of variables, and even if the y (i,j) can be assumed to be real, (RP1) is too big to be solved by standard solvers. A way to get a tractable formulation consists in reducing the number of variables. Therefore, we now dene an integer program which contains only the z (i,j) variables and we show that solving the new integer program is equivalent to solving the former one.

(RP 2) z(RP 2) = min (i,j)∈A c (i,j) z (i,j) (21) s.t. j∈V z (i,j) = j∈V z (j,i) ∀i ∈ V (22) 0 0 0 1 (0; 1) 0 (2; 0) (0; 2) 2 0 (1; 0) 1 (0; 0) Figure 7. A solution of (RP1) that is not a solution of the SVOCPDP i∈V \{0} z (0,i) = 1 (23)
(i,j)∈δ + (S)

z (i,j) ≥ µ(S), ∀S ⊆ V \ {0} (24)
(i,j)∈δ + (S)\δ(0)

z (i,j) ≥ d(S) Q ∀S ⊆ V (25) z (i,j) ∈ Z + , ∀(i, j) ∈ A (26)
Any feasible solution of the SVOCPDP induces z (i,j) that satisfy the constraints of program (RP2). Variables y (i,j) disappear in (RP2). However, for any feasible solution z (i,j) of (RP2), there is a feasible solution to (RP1) with the same z (i,j) (and hence the same value of the objective function), and conversely. This fact is summarized in the following proposition.

Proposition 4. Let y, z be a feasible solution of (RP1). Then z is a feasible solution of (RP2). Conversely, let z be a solution of (RP2). Then there exists y such that y, z is a feasible solution of (RP1).

This property is also proven by [START_REF] Hernandez Pérez | The one-commodity pickup-and-delivery travelling salesman problem[END_REF] when z encodes an Hamiltonian circuit. Here the proposition is proven in a more general case, and the proof is maybe slightly simpler since it does not involve Bender's decomposition but the cut condition for b-ows.

Proof. Take a feasible solution y, z of (RP1). We show that the z (i,j) satisfy (RP2). The only thing that has to be checked is that the z (i,j) satisfy constraints (25). Let S ⊆ V . Using constraints ( 16) and (19), aggregated over (i, j) ∈ S we get (i,j)∈δ + (S)\δ(0)

z (i,j) ≥ 1 Q (i,j)∈δ + (S)
y (i,j) .

Hence (i,j)∈δ + (S)\δ(0)

z (i,j) ≥ 1 Q   (i,j)∈δ + (S)
y (i,j) -

(i,j)∈δ -(S)
y (i,j)   .

Since (i,j)∈δ + (S)

y (i,j) - (i,j)∈δ -(S)
y (i,j) = i∈S   j∈V \{i} y (i,j) - j∈V \{i} y (j,i)   ,
we get with the help of constraints ( 17) that (i,j)∈δ + (S)\δ(0)

z (i,j) ≥ 1 Q d(S).
Conversely, take a feasible solution z of (RP2). The only thing that has to be checked is that there exists a non-negative b-ow y with b i = p iq i for all i ∈ V with a capacity = Qz (i,j) on each arc (i, j). But constraints (25) are precisely the well-known cut condition for ows on networks (see for instance [START_REF] Gallai | Maximum-minimum theorems for networks[END_REF]).

The integrality is a consequence of the integrality of b i .

Note that the two programs (RP1) and (RP2) are NP-hard, since they obviously contain the TSP as a special case. Moreover, the continuous relaxation of (RP2) provides a better lower bound than the continuous relaxation of (RP1) thanks to • in constraints (25).

Relaxation vs original problem

The hardness of SVOCPDP, already emphasized in Subsection 1.1, appears also through the following four propositions, which show that even if we have a feasible solution or an optimal solution of (RP1) or (RP2), it is an NP-complete problem to decide whether this solution is induced by a feasible solution for the SVOCPDP.

Proposition 5. Let y, z be a feasible solution of (RP1). Deciding whether there is a feasible solution of the SVOCPDP inducing y, z is NP-complete.

Proof. Let b 1 , . . . , b r be r non-negative integers with l b l even. Dene m = 1 2 r l=1 b l . Consider then the graph of Figure 8. It encodes a feasible solution y, z of program (RP1): each arc (i, j) has z (i,j) = 1 and each number next to it is the corresponding value for y (i,j) . Assume that the capacity of the vehicle is m + 1. We are going to prove now that there is a feasible solution for the SVOCPDP inducing such y (i,j) , z (i,j) if and only if there is a subset I ⊆ {1, . . . , r} such that l∈I b l = m, whence showing that deciding whether such a route exists is NP-complete. Note that there are exactly 2m + 1 bikes in the network.

Assume that such a feasible solution exists. Consider again Figure 8 and the following instant: the vehicle takes arc (u, u ). At this time, the vehicle carries m + 1 bikes. Therefore, it has already taken the arc (s, u) (in order to have m bikes) and exactly one of the arcs (v, v ) (in order to get the remaining bike). It cannot have already taken the other arc (v, v ), otherwise it should go back to the depot. It has not taken the arc (u , s) yet. Thus, m + 1 bikes among the 2m + 1 are on the vehicle, and m bikes are on vertex v : there is no bike left on vertex s. When the vehicle goes back to vertex s, using arc (u , s), it carries exactly m bikes. These m bikes have to be carried to vertex v, using arcs (s, v). Hence, the feasible solution selects (s, v) arcs whose b l sum exactly to m.

Conversely, assume that there is a subset I such that l∈I b l = m. The sequence starting with the depot, then going through s, u, s in this order, then going back and forth between s and v, using the arcs indexed by I, to carry m bikes to v, then going through v , u, u , s in this order, then using the remaining arcs between s and v, and nally going through v, v in order to nish again at the depot is a feasible solution for the SVOCPDP. Those two propositions deal with feasible solutions of programs (RP1) and (RP2). Now, even if we have a special solution, especially the optimal one, the question whether it is induced by an optimal one of the SVOCPDP is NP-complete. We have indeed the following propositions.

Proposition 7. Let z * be an optimal solution of (RP2). Deciding whether there is a feasible solution of the SVOCPDP inducing z * is NP-complete.

Proof. See Appendix. Proposition 8. Let y * , z * be an optimal solution of (RP1). Deciding whether there is a feasible solution of the SVOCPDP inducing y * , z * is NP-complete.

Proof. See Appendix.

Lower Bound

The integer program (RP2) has an exponential number of capacity and connectivity constraints. The problem is solved through a branch-and-cut algorithm. At each node of the branch-and-cut tree, the continuous relaxation of the problem (RP2) is solved, but only a subset of constraints (24) and (25) are activated.

The continuous optimal solution z is checked by dierent routines that try to determinate violated constraints. If a violated constraint is found, it is added into the linear relaxation which is again optimally solved. If no violated constraint is found, but z is still fractional, branching goes on.

If the solution is integer, then it is a feasible solution of (RP2). Its value is kept and can be used as an upper bound on the optimal value of (RP2). When the algorithm reaches the end, updating the upper bound during the branch-and-cut process gives at the end the optimal solution of (RP2). If it stops before (for instance if there is a time limit), we get at the end the best feasible solution encountered during its exploration and a lower bound on the optimal solution value, used for calculating the gap. 7.1. Separation of connectivity constraints. Given a solution z, computed at a node of the branch-andcut tree, we check that constraints (24) are satised as follows. We consider the support graph G[z] together with a capacity equal to z(i,j) on each arc (i, j). For each unbalanced vertex i, we compute the minimum cut δ + (S) separating i from the depot 0 (with the algorithm by [START_REF] Edmonds | Theoretical improvements in algorithmic eciency for network ow problems[END_REF]). If this minimum cut has a value strictly lower than 1, then we add the corresponding constraint to the linear program. 7.2. Separation of capacity constraints. As already noted when we have dened (RP2), there is no known polynomial algorithm that checks the capacity constraints (25).

The following constraints are less tight but can be separated in polynomial time.

(i,j)∈δ + (S)\δ(0)

z (i,j) ≥ d(S) Q , (∀S ⊆ V ) (27)
These constraints are the relaxation of the constraints (25) in (RP2). They are called relaxed capacity constraints. A method for nding violated relaxed capacity constraints for the CVRP is proposed by

[ABB + 99].
In the SVOCPDP, the dierence is that a vertex can either be in excess or in default. Therefore their method has been adapted to t the specic characteristics of the SVOCPDP. The two arcs incident to the vertex 0 are deleted and two new vertices s and t are added. Vertex s is linked to all vertices in excess with the capacity κ (s,i) = di Q whereas all vertices in default are linked to vertex t with the capacity κ (i,t) = -di Q . The capacity is equal to z(i,j) on the original arcs (i, j). We compute an s-t min-cut (again with the Edmonds-Karp algorithm). Let X be the subset of vertices such that δ + (X) is the s-t min-cut. s is in X and we dene S := X \ {s}. We still dene S as V \ S, where V is the vertex set of the original graph, and contains neither s nor t.

capacity of δ + (X) = i∈S\{0},j∈ S\{0} z(i,j) + i∈ S\{0} κ (s,i) + i∈S\{0} κ (i,t) (28) = (i,j)∈δ + (S)\δ(0) z(i,j) - d(S) Q + i∈V : d(i)>0 d(i) Q . ( 29
)
If the capacity of δ + (X) minus i∈V :

d(i)>0 d(i)
Q (which is a xed value) is less than 0, it means that a relaxed capacity constraint is violated and the corresponding constraint are added to the linear relaxation.

Otherwise, there is no violated relaxed capacity constraint as shown in Figure 9, which represents a possible solution of the linear relaxation when solving the instance presented in Figure 2, but with a capacity of the vehicle set to Q = 10. The fractional value of each arc is given next to it.

When the relaxed capacity constraints are respected a second procedure runs looking for violated capacity constraints. The interest to look for these constraints in addition to the former is to increase the speed of the general algorithm. Indeed, we get in this case tighter constraints. The used procedure is again an heuristic presented in [ABB + 99], since there is no known polynomial algorithm. It is a tabu search that tries to nd a subset S violating the constraints. A short explanation of the tabu search method is done later in Section 8. The main idea is starting from a subset S of vertices, vertices are either added to S or removed from S in order to nd a new subset where capacity constraints are violated. Roughly speaking, the vertex that is added to or removed from the subset is kept in a tabu list and cannot be used for a given number of iterations. If the procedure nds a violated capacity constraint, it is added to the linear program.

7.3. Initial relaxation, separation strategy and branching rules. The initial relaxation solved is (RP2) without (24) and with (25) only for S = {i} for all i ∈ V . The separation routines are called with respect to the following order:

• connectivity constraints (24)

• relaxed capacity constraints (25)

• capacity constraints (27)

The branch-and-bound tree framework SCIP has been used for the computational results. Two dierent branching strategies have been tested. The rst one is the so-called reliability branching rule dened in their paper (see in [START_REF] Achterberg | Branching rules revisited[END_REF] for more details), which is the one that seems the most eective one for solving SVOCPDP after several tests that have been done. For the node selector, the rule is the best estimate search rule, which is the default rule of SCIP.

The second one consists in branching rst on j∈V z (j,i) . If there is a vertex i for which this quantity is equal to a non-integral u, we add two constraints: j∈V z (j,i) ≥ u and j∈V z (j,i) ≤ u .

The vertex for which uu is the nearest from 0.5 is selected in priority. When there is no such vertex, the branching uses the former strategy. This second strategy (called in the following degree branching rule) is interesting since it based on the idea of branching on constraints instead of branching simply on variables.

For the node selector, the rule is the same of the previous branching. 7.4. Adding valid inequalities. The clique cluster inequalities dened in Section 3.1 of [HPSG07] can be straightforwardly adaptated. We have implemented them but the fact that the variables are not restricted to binary values have lead to huge computation times. It comes partially from the fact that their ecient version (7) of that paper is no longer valid. The still valid version is separated in a more dicult way. The results were dominated by the version without them. We have not included the results in the present paper.

It seems that adding cuts in ecient way for the SVOCPDP is an open problem. Other inequalities could be adaptated, such as the comb inequalities. We have not tried them.

Upper Bound

The algorithm chosen for computing an upper bound of the optimal value is a tabu search. Several problems like: job shop scheduling, graph coloring, traveling salesman problem and other vehicle routing have been successfully tackled by the tabu search. For an introduction to this method, see for instance [START_REF] Glover | Tabu search[END_REF].

While a greedy local search stops when no improving move is found in the neighborhood N (s) of the current solution s, the basic principle of a tabu search is to continue the search by allowing non-improving moves.

To avoid cycling on the same solution, the last moves are kept in the tabu list that is updated at each iteration. Its size is a parameter.

Various stopping criteria can be used, such for example: the predetermined number of iterations with no improvement, the maximum number of iterations, the maximum amount of time spent in the method.

To dene a tabu search properly, we have to describe how to compute the cost of the current solution, the neighborhood, the way to encode the tabu list, a heuristic to build the rst current solution and the stopping criterion.

8.1. Cost of the current solution. Proposition 1 allows to encode the solution as a sequence of vertices, starting and ending with the depot 0, without specifying bike displacements. The score of a sequence i 1 → i 2 . . . → i k is the cost of the traveled distance

k-1 j=1 c (ij ,ij+1
) to which we add a penalty when the sequence is infeasible, i.e. when there is no feasible solution for the SVOCPDP inducing this sequence. The feasible sequence space may be disconnected, whence to explore eciently the used neighborhood is enlarged allowing the tabu in visiting infeasible solutions. Then, the infeasibility is penalized in the objective function.

Thanks to Proposition 1 we are able to evaluate the infeasibility of any sequence indicated in the following with s. It is the degree of infeasibility we have dened in Section 3. The relaxed constraints are (8) and ( 9) and the resulting cost function f (s) dened in (30) is inspired by the one proposed by [START_REF] Gendreau | A tabu search heuristic for the vehicle routing problem[END_REF] for the VRP and is the following:

(30)

f (s) = k-1 j=1 c (ij ,ij+1) + γ i∈V (p i -p i )
where γ is a positive constant and the p i are computed through the max-ow algorithm of Proposition 1.

Recall that according to that proposition, p i ≤ p i .

In our experiment we have dened γ as 10 times the mean distance of a direct trip from the depot to a vertex. This choice is a way to convert at the right scale in term of cost the number of bikes that remain to be displaced (and was experimentally tuned).

Initial solution.

We propose two distinct methods to compute the initial sequence needed for the tabu search.

8.2.1. Greedy heuristics. The method rst tries to close vertices. Closing a vertex means bringing a vertex from its original state to its target state q i in one move. The rst criterion used, for deciding which vertex visit rst, consists in ranking the unbalanced vertices with respect to their distance from the vehicle position and then in choosing the nearest vertex it can close. If it is impossible to close any vertex, then for each vertex the number of bikes the vehicle can load or unload is computed and the vehicle is driven toward the vertex with the highest absolute number of bikes that can be exchanged (i.e., either loaded or unloaded). In case of equality between several vertices, it goes to the nearest one. With this method we are sure to start from a feasible sequence.

The solution of (RP2).

A solution z of relaxation (RP2) of Section 5 is such that G[z] is an Eulerian graph. A closed Eulerian path in the supported graph G[z] provides a sequence of vertices that can be used as an initial solution of the tabu search. To compute such a closed Eulerian path various classical methods are available (see for instance page 31 of the book by [KV02]). Since the relaxation is quite good, this sequence is often not too far from an optimal solution. Note that if the sequence is feasible for the SVOCPDP, it is an optimal solution. Unfortunately, even if we have an instance for which the relaxation is tight, there is no polynomial algorithm for nding the optimal solution of the original problem (Proposition 7 above, Section 6).

8.3. Neighborhood description. At each iteration the whole neighborhood is explored. The denition of neighborhood is essential, since dierent neighborhoods oer dierent ways to explore the solution space but require dierent computational eorts. For solving the addressed problem, four dierent moves have been dened and used in the tabu search framework. Note that any time a vertex appears two consecutive times in a sequence one of the occurrences can be removed.

Let k denote the length of the sequence and note that k is in general larger than n. We therefore evaluate the complexity for computing the neighborhood in terms of k instead of n. To illustrate the dierent moves the solution displayed on Figure 2 is used each time as initial solution.

0 → 7 → 6 → 4 → 5 → 2 → 8 → 1 → 0
Following this sequence, the vehicle leaves the depot to go to vertex 7. Then it goes on with vertex 6 and so on until vertex 1 from where it goes back to the depot. 8.3.1. 2-OPT. This move is classical for routing problems. A pair of non consecutive arcs are removed from the sequence and two new arcs are inserted. They link the two tails and the two heads of the removed arcs and therefore the travelling direction of the subset of nodes between the two removed arcs is reversed to obtain a new sequence. An example of 2-OPT is obtained by removing arcs (0,7) and (5,2), by inserting arcs (7,2) and (0,5) and reversing the sequence 5, 4, 6, 7, therefore the following sequence is obtained:

0 → 5 → 4 → 6 → 7 → 2 → 8 → 1 → 0
This move is tested for each pair of arcs in the sequence. It takes O(k 2 ) to test all the moves of this type. 8.3.2. Suppression. This move tries to delete a vertex in the sequence. If we delete vertex 4 in the example of Figure 2, we obtain

0 → 7 → 6 → 5 → 2 → 8 → 1 → 0
This move is tested for each vertex in the sequence. It takes O(k) to test all moves of this type. 8.3.3. Add unbalanced vertex. This move is active when the current sequence is infeasible. All the vertices are checked and both the most unbalanced vertex i in excess and the most unbalanced vertex j in default are selected (vertices for which |(q iq i ) -(p ip i )| are maximal). One neighbor is obtained by adding, just after i in the sequence, the moves → j → i →. The other neighbor is obtained by adding the moves → i → j → in the sequence just after j. If neither i nor j are present in the sequence, then there is only the neighbor obtained by adding → i → j at the end of the sequence. Testing all moves of this type takes O(k 2 ).

Once more with the example reported in Figure 2, if the capacity of the vehicle is lower than 8, the vertices are not balanced at the end of the sequence. Let Q = 6, the vehicle cannot bring vertices 1, 6 and 8 to their target state: 2 bikes need to be loaded from 6 while 1 bike is requested by both node 1 and node 8. The moves → 8 → 6 are tried to be added and we obtain the following sequence

0 → 7 → 6 → 8 → 6 → 4 → 5 → 2 → 8 → 1 → 0
This sequence is a feasible solution. 8.3.4. Add buer vertex. This move tries to add a second time a vertex in the sequence. The second copy of a vertex can be used as buer vertex, since in the buer vertex bikes are either delivered or picked up. Thus, if a vertex is used as a buer vertex, it has to be visited at least twice. Every vertex can be used as a buer vertex. A copy of each vertex already present in the sequence is tried to be inserted at each possible position. If a vertex is not present in the sequence (for example the initially balanced vertices), then it is inserted twice. This move is very time consuming and testing all moves of this type takes O(k 3 ). Thus it is not called at each iteration but in 20% of the cases. 8.4. The tabu list. The fact that arcs could be used more than once leads to a management of the tabu list slightly dierent from the classical one used, for example, for solving classical routing problems. When a move is made, some arcs leave the solution and others enter the solution. Each exiting arc is kept in the tabu list along with the position of its occurrence in the last solution. This arc is forbidden for iterations, but it is allowed to be reinserted in a dierent position.

The tabu list is kept using the multimap structure in the STL of C++, which ensures a quick access.

8.5. The tabu search. The components previously described are integrated in the general algorithm summarized in Algorithm 1. The following notations are used in the pseudo-code: s is the current solution; s * is the best feasible solution encountered from the beginning of the tabu search. f (s) is the value of the objective function for the solution s; N bIterM ax is the maximum number of iterations done before the tabu search stops. i is the current iteration; sX is the best solution encountered using the move X; s # is the best solution in the whole neighborhood. a is a random variable uniformly distributed between 0 and 1. ) and tested on a PC AMD Athlon 5600+ clocked at 2.8 GHz, with 16 GB RAM. To the best of our knowledge, this paper is the rst solving instances of the SVOCPDP, therefore the following Subsection 9.1 presents how the instances have been created, the results obtained are reported in Subsection 9.2 while Subsection 9.3 reports a discussion on the results. 9.1. Instances. The instances are taken from [HPSG04a] which denes, for various values of n, 10 instances named from A to J, with an algebraic demand di between -10 and 10 for each vertex i. We work on these instances for n ∈ {20, 40, 60, 100}. In order to get an initial state p i and a target state q i for each vertex i, as required for the SVOCPDP, the demands have been modied. We do it twofold. A rst set of instances has been obtained by setting p i = 10, q i = 10 + di and C i = 20 for each vertex i. We refer to these instances by the indication pi = qi = 10 (average state). A second set of instances is obtained by setting p i = 30, q i = 30 + 3 × di and C i = 60 for each vertex i. We refer to these instances by the indication pi = qi = 30

(average state). . An optimal solution already found by the branch-and-cut for the instance n20G, with Q = 10 and pi = qi = 10 9.2. Computational results. The proposed algorithms have been tested on these instances for all Q ∈ {10, 30, 45, 1000}. Note that even with huge capacity for the vehicle our problem is not a TSP with pickup and delivery. Initially balanced vertices can be skipped, nevertheless visiting a vertex several times can improve the cost of the solution (see for example the solution reported in Figures 4 and5). The algorithms tested are the tabu search initialized with the greedy heuristics from Subsection 8.2.1, the branch-and-cut of Section 7 which provides a lower bound for the problem, and the tabu search initialized with the solution of the branch-and-cut as explained in Subsection 8.2.2.

The results achieved by all our algorithms for solving the generated instances are reported respectively in Tables 1, 3, 5 and 7 for the reliability branching rule and the tabu search, and Tables 9, 11, 13 and 15 for the degree branching rule. For each pair n, Q, we give the results for two instances the one achieving the smallest gap and the one achieving the largest one and also the average performance on the ten instances.

We have treated the results for n = 100 in distinct tables since the results obtained for such dimensions are less good than those obtained with n < 100. The following webpage gathers the best results we have obtained so far for the complete set of instances. http://cermics.enpc.fr/∼meuniefr/SVOCPDP.html

The tabu search initialized with the greedy heuristic (indicated in the following with TS1) has the following set of parameters: maximum number of iterations = 1 000, the size of the tabu list (i.e. ) = 30 and the maximum amount of time = 1 000 seconds. The method also stops when no improvement has been done for 80 consecutive iterations. The tabu search initialized with the solution obtained by the branch-and-cut LB is a lower bound on the optimal solution of the relaxation. LBr is the lower bound that has been found at the root node. Gap % is the percentage gap between UB-BC and LB.

Time is the computational time. # Nodes is the number of nodes of the Branch-and-Cut tree. Figure 10 gives the solution for the instance n20G with Q = 10, and an average number of bikes per station equal to 10. The depot appears as a square with a 0 and (p i ; q i ) is given next to each vertex. The solution displayed is the one obtained at the end of the following sequence of algorithms: the tabu search rst, then the branch-and-cut with the result of the tabu search as upper bound, and, nally, the tabu search which starts from the solution given by the branch-and-cut. In this example the route displayed is exactly the solution of the branch-and-cut. It is the optimal solution even if a vertex (vertex 18) is visited twice.

Compare with Figure 11. It is the same instance, but this time with an average number of bikes per station equal to 30. Again, the branch-and-cut is able to nd the optimal solution. Interestingly, we can see that, with an increasing average number of bikes but with the same vehicle capacity, optimal solutions may visit vertices several times. 9.3. Discussion of the performances. As expected, the computational time and the number of nodes of the branch-and-cut algorithm increase with the number of vertices and decrease with the capacity of the vehicle. In Tables 2, 4, 6 and 8 (and 10, 12, 14 and 16), the number of nodes used in the branch-and-cut could seem really huge. However, since variables are integers (not binary), several branchings must be done on the same variable before obtaining the optimal solution. Except for the large instances with 100 vertices, the results reported in Tables 1 and5 show that the distance between the best upper bound found and the lower bound is really small. The gap is on average less that 5%. The local search is very ecient on small and medium instances (up to 60 vertices) and becomes less eective for larger instances. This is probably a consequence of the size of the neighborhood, which can be quite huge when the vehicle has to make several visits at some vertices. For n = 100, the second tabu search makes sometimes only one or two iterations.

Note also that the smaller is the capacity, the harder is the problem. It is in accordance with the intuition, since for instances with small vehicle capacity, the mean number of visits by vertex increases. Instances n20B for Q = 10 and Q = 45 and n40E for Q = 10 and Q = 30 in Table 1 are illustrations for such a phenomenon.

We have also such a phenomenon for n20B and n40E in Table 5 for Q = 10 and Q = 30.

The reliability branching rule seems to be slightly better than the degree branching one. However, this latter improves sometimes the lower bound: for instance, n60I, with Q = 30 and 10 bikes per station in average (Tables 1 and9), or n60D, with Q = 45 and 30 bikes per station in average (Tables 5 and13).

If we go back to the size of the Vélib system in Paris (one truck of capacity 20 for 50 vertices of capacity 30), the instances that are close to these numerical features are the instances with n = 60, Q = 30, and 10 bikes in average per station (Tables 1 and9) and the instances with n = 40, Q = 30 and 30 bikes in average per station (Tables 5 and13). We get solutions that have most of the time a gap smaller than 2%. It would certainly be possible to improve the time consumption of the ow algorithm in the tabu search. Stopping the branch-and-cut after a xed time would also be a reasonable solution to reduce the total time to get a good solution. However, we see that the instances of this kind can already be solved within an average optimality gap of less than 2% in less than 40 minutes. solution q obtained with the same bike displacements except for these bikes, which are not allowed to leave their initial vertices. Denoting by pi the number of bikes on vertex i that are not allowed to move, we have i pi = i: qi<q i q iq i and qi = q i + pi for each vertex i such that q i ≥ q i q i + pi for each vertex i such that q i < q i Note that pi ≥ qiq i . The distance to q is equal to

||q -q|| 1 = i: qi≥q i |q i -q i | + i: qi<q i |q i -q i | ≤ i: qi≥q i |q i -q i | + i: qi≥q i pi + i: qi<q i pi = i |q i -q i | = ||q -q|| 1 .
We can therefore choose the target state q ∈ Q minimizing ||q -q|| 1 such that at least max(0, qiq i ) bikes have not been moved for each vertex i. Let us now dene, for the route with these target states, p M i the number of bikes among the p i initial ones that have left vertex i and q M i the number of bikes that have been brought to vertex i. Note that according to our choice for q, we have q M i ≤ q i (and of course p M i ≤ p i ).

We have

(31) ||q -q|| 1 = i∈V |q i -q i | = i∈V |p i -p M i + q M i -q i |.
Let δ i := min(p ip M i , q iq M i ) and dene p i := p M i + δ i and q i := q M i + δ i for each vertex i. Note that the initial and target states p i and q i are feasible in the sense of Proposition 1. According to Equation (31), we have ||q -q|| 1 = i∈V |p ip i + q iq i |. Using the fact that for each i, at least one of p ip i and q iq i is equal to 0, we get that ||q -q|| 1 = 2 i∈V (p ip i ).

(b) Conversely, assume that we have loading and unloading operations compatible with initial and target states p i and q i , with

• 0 ≤ p i ≤ p i and 0 ≤ q i ≤ q i for each vertex i • i p i = i q i • the quantity i∈V p i is maximal.

We show that these operations are also compatible with the initial state p and a nal state q ∈ Q such that ||q -q|| 1 = 2 i∈V (p ip i ).

Since C = +∞, if p i < p i and q i < q i , we could have taken into account one more bike on i (increasing by one p i ), which would have not left vertex i. Therefore, for each i, we have at least one of p ip i and q iq i that is equal to 0 and thus, we have 2 i∈V (p ip i ) = i∈V |p ip i + q iq i |. Let us now consider what we have as a nal state when we start instead with p i , while making the same loading and unloading operations: we get qi := p ip i + q i bikes on each vertex i, for which we have ||q -q|| 1 = 2 i∈V (p ip i ).

According to (b), starting from initial state p, the same loading and unloading operations as for an optimal solution in the sense of Proposition 1 lead to a nal state q such that ||q -q|| 1 = 2 i∈V (p ip i ). (a) shows that this q minimizes ||q -q|| 1 over Q, as required.

Proof of Proposition 6. The proof is very similar to the one of Proposition 5. Again, we work with a reduction from the 2-partition and with the graph of Figure 8. In order to adapt the proof for (RP2), we simulate the y (i,j) by subdividing each arc a = (i, j) in ve arcs, introducing four new vertices u a , u a , v a and v a , with p ua = 0, q ua = y (i,j) , p va = m + 1, q va = 0, p v a = 0, q v a = m + 1, p u a = y (i,j) , q u a = 0. Proof of Proposition 7. We can force the solution built in the proof of Proposition 6 to be optimal, and then the same proof does the job. To force the solution to be optimal, we consider now the graph G obtained in the proof of Proposition 6 by subdividing each arc in ve new arcs. We assume that the cost to traverse each of these arcs is 1. Now, we build from this graph a complete graph for which each arc (i, j) gets for cost the cost of a shortest path between i and j. This complete graph provides an instance of the SVOCPDP. The arc of G encodes a solution z for (RP2), which is optimal since each new vertex u a , u a , v a , v a has to be visited.

[Proof of Proposition 8. There is only one solution for y once z has been dened according to the proof of Proposition 6, whence the same proof as for Proposition 7 works.

(enpc) CERMICS, Université Paris est (Paris13) LIPN (UMR 7030), 99, av. Jean-Baptiste Clément, 93430 Villetaneuse, France.

FigureFigure 1 .

 1 Figure 1. Example of an instance.

Figure 2 .

 2 Figure 2. Example of a feasible solution (i.e. a route)

Figure 3 .

 3 Figure 3. Drops can help optimality

Figure 6 .

 6 Figure 6. Example of the graph used in the algorithm of Proposition 1

  Constraints (22), (23) and (24) are similar to those of the rst relaxation (RP1). Constraints (25) are the capacity constraints requiring that, for any subset S ⊆ V , the vehicle goes at least |d(S)| Q times into S. The absolute value | • | is useless in (RP2). Indeed, if d(S) ≥ 0, it is not necessary. If d(S) < 0, (25) are redundant and in anyway, when (25) is written with S, we obtain precisely what would have been obtained with the |.| for S. Constraints (25) are well-known in the CVRP context and it is still an open question whether a polynomial algorithm separating these constraints exists (see for example [ABB +99]). In Subsection (7.2), we explain how to deal with them.

Figure 8 .

 8 Figure 8. Proof of NP-completeness

Figure 9 .

 9 Figure 9. Construction of the new graph from the graph of Figure 2 to check the relaxed capacity constraints

  Figure10. An optimal solution already found by the branch-and-cut for the instance n20G,

Figure 12 .

 12 Figure 12. Construction of G in the proofs of Propositions 6, 7, 8

Table 1 .

 1 Performance of the overall method for pi = qi = 10, less than 60 vertices, and a vertex capacity equal to 20 reliability branching

	Instance name n Q UB1 Time Iter. UB2 T.t. n20A 20 10 4764 9 132 4702 10 4702.00 LB Gap % 0.00 n20B 20 10 4776 7 124 4769 8 4769.00 0.00 20 10 7 126 9 0.00 n20D 20 30 4156 4 93 4089 5 4089.00 0.00 n20E 20 30 4556 5 88 4556 5 4299.00 5.97 20 30 5 105 6 1.13 n20B 20 45 3792 5 102 3792 5 3792.00 0.00 n20C 20 45 4045 8 88 4045 8 3891.00 3.96 20 45 6 129 7 0.83 n20B 20 1000 3792 7 90 3792 7 3792.00 0.00 n20C 20 1000 4042 6 138 4042 6 3891.00 3.88 20 1000 7 138 8 0.83 n40B 40 10 6377 198 121 5949 468 5949.00 0.00 n40F 40 10 7373 391 216 7138 10473 6651.57 7.31 40 10 367 220 4266 1.54 n40B 40 30 5297 142 91 5110 145 5110.00 0.00 n40C 40 30 4806 261 179 4692 262 4644.00 1.03 40 30 225 157 228 0.10 n40B 40 45 5297 141 136 5110 144 5110.00 0.00 n40C 40 45 4941 185 113 4656 187 4522.00 2.96 40 45 227 110 230 0.53 n40B 40 1000 5297 152 91 5110 157 5110.00 0.00 n40H 40 1000 5060 183 128 5006 187 4772.00 4.90 40 1000 273 130 279 0.96 n60H 60 10 9100 1007 76 8123 11043 7747.60 4.84 n60F 60 10 9234 1004 54 8374 11034 7387.74 13.35 60 10 1003 69 11036 10.48 n60J 60 30 6543 1026 85 6389 1177 6389.00 0.00 n60I 60 30 6918 995 77 6480 11028 6153.33 5.31 60 30 1013 73 2218 1.05 n60J 60 45 6602 988 60 6374 1047 6374.00 0.00 n60H 60 45 6122 1037 51 6122 1055 5825.00 5.09 60 45 998 67 1046 2.25 n60D 60 1000 6345 1022 57 6223 1946 6223.00 0.00 n60I 60 1000 6387 988 65 6355 1469 5914.00 7.45 60 1000 1007 67 1255 2.19

Table 5 .

 5 Performance of the overall method for pi = qi = 30, less than 60 vertices and a vertex capacity equal to 60 reliability branching Instance name n Q UB1 Time Iter. UB2 T.t.

	LB Gap %

(indicated in the following with TS2) has the following parameters: the tabu list size is set to 3, the maximum number of iterations to 50 and the maximum amount of time to 300 seconds. The search is stopped whenever 15 consecutive iterations did not improve the best solution or if a solution is found with a cost equal to the lower bound obtained by the branch-and-cut. The abbreviations used in these tables are the following: n is the number of vertices. Q is the vehicle capacity. UB1 is the value obtained by TS1. Time is the cpu time (in seconds) used by the Tabu Search TS1. Iter. is the number of iterations. UB2 is the best solution obtained by the tabu searches TS1 and TS2. T.t. is the cpu time (in seconds) elapsed for executing TS1+B&C+TS2. LB is the best lower bound found by the branch-and-cut on (RP2). Gap % is the percentage gap between the previous LB and UB2.

The results achieved by the branch-and-cut alone are reported respectively in Tables 2, 4, 6 and 8 for the reliability branching rule, and Tables 10, 12, 14 and 16 for the degree branching rule. The maximum time was set to 10 000 seconds. The abbreviations used in these tables are the following: n is the number of vertices. Q is the vehicle capacity. UB-BC is the best solution found by the branch-and-cut. , A branch-and-cut algorithm for a traveling salesman problem with pickup and delivery, Discrete Applied Mathematics 145 (2004), 126 139.

[HPSG04b]

, Heuristics for the one-commodity pickup-and-delivery traveling salesman problem, Transportation Science 38 (2004), 245255.

[HPSG07] , The one-commodity pickup-and-delivery travelling salesman problem: inequalities and algorithms, Networks 4 (2007), 258272.

Appendix

Proof of Proposition 2. Denote by Q the set of all target states qi such that there exists a feasible solution for initial states p i and target states qi inducing the sequence. We are interested in the solution of min q∈Q ||q -q|| 1 .

(a) We rst show that for any solution of min q∈Q ||q -q|| 1 , there exist loading and unloading operations compatible with initial and target states p i and q i , such that • 0 ≤ p i ≤ p i and 0 ≤ q i ≤ q i for each vertex i • i p i = i q i and such that ||q -q|| 1 = 2 i∈V (p ip i ).

We have the following central observation: let q ∈ Q; there exists q ∈ Q such that ||q -q|| 1 ≤ ||q -q|| 1 and such that, on each vertex i, there are at least max(0, qiq i ) bikes that have not been moved. Indeed, for each vertex i, choose max(0, q iq i ) bikes among the q i bikes at the end of the route. Consider now the