H L M Kerivin

M Lacroix

A Quilliot

H Toussaint

Heuristic for the preemptive asymmetric stacker crane problem

Keywords: preemptive asymmetric stacker crane problem, reloads, routing, local search, heuristic design

In this paper, we deal with the preemptive asymmetric stacker crane problem in an heuristic way. We first present some theoretical results which allow us to turn this problem into a specific tree design problem. We next derive from this new representation a simple, efficient local search heuristic, as well as an original LIP model. We conclude by presenting experimental results which aim at both testing the efficiency of our heuristic and at evaluating the impact of the preemption hypothesis

Introduction

Pickup and delivery problems, which consist in scheduling the transportation of sets of goods/passengers from origin nodes to destination nodes while using a given set of vehicles, have been intensively studied for decades from both theoretical and practical points of view. Many variants have been considered and a lot of methods have been designed in order to improve the resolution of such problems. (One can refer to [START_REF] Berbeglia | Static pickup and delivery problems: a classification scheme and survey[END_REF][START_REF] Parragh | A survey on pickup and delivery problems: Part II: Transportation between pickups and delivery locations[END_REF][START_REF] Savelsberg | The general pickup and delivery problem[END_REF] for surveys on these problems and methods.) Among all the pick-up-and-delivery like problems which have been addressed by searchers, the Stacker Crane Problem is characterized by the fact that only one vehicle is involved, which can deal with only one demand unit at the same time. In this paper, we give a heuristic to handle what we call the Preemptive Stacker Crane Problem, that means the case when some demands can be dropped anywhere in the network and reloaded afterwards in order to gain time.

A rough description of the Stacker Crane Problem (SCP) can come as follows: G being some transit network whose oriented links or arcs are endowed with lengths or costs and which is provided with some specific Depot node, we are required to schedule the route of a single vehicle V, which is required to address a Demand set K, each demand k K being defined by some origin node o(k) and by some destination node d(k). Namely, addressing the demand k means transporting some unique load unit L(k) from o(k) to d(k) while using the vehicle V, whose capacity is such that it cannot contain more than one load unit L(k), k K, at a given time. Thus, scheduling V means designing a route inside the network G, which is going to start and end in Depot and to make possible for V to handle every demand k in K, and solving the Stacker Crane Problem will mean computing this route in such a way that this route is the shortest possible. Two versions of the SCP may be distinguished. In the first one, called non-preemptive Stacker Crane Problem (NPSCP), every demand has to be directly carried from its origin to its destination. In the second version, which is called Preemptive Stacker Crane Problem (PSCP), any load unit L(k) related to demand k may be dropped (unloaded) at any node x of the transit network G, before being reloaded a little further and this unload/reload process, which we call reload process, may be performed several times before the load unit L(k) reaches the destination node d(k). In case the cost or length function, which to any arc (x, y), make correspond some length DIST(x, y), is symmetric, we talk about Symmetric SCP (symmetric NPSCP or symmetric PSCP), and in the case the converse is true, we talk about asymmetric SCP (asymmetric NPSCP or asymmetric PSCP).

The Stacker Crane problem was first introduced by Frederickson et al. in [START_REF] Frederickson | Approximation Algorithms for Some Routing Problems[END_REF], under its non preemptive symmetric form. They proved its NP-hardness by using a reduction from the TSP. They also got a 9/5-approximation scheme for this problem. Moreover, they proposed a natural extension of this problem to n identical vehicles, and obtained for this extension a (1+ -1/k)-approximation scheme, where  corresponds to the bound of the approximation for the problem with only one vehicle.

Atallah and Kosaraju [START_REF] Atallah | Efficient Solutions to Some Transportation Problems with Applications to Minimizing Robot Arm Travel[END_REF] were the first to consider the preemptive version of the symmetric SCP. They studied both non-preemptive and preemptive versions of the symmetric SCP in the case when the underlying graph is an elementary path or an elementary cycle. They proved that in such a case, both versions are polynomial-time solvable. Frederickson and Guan [START_REF] Frederickson | Preemptive Ensemble Motion Planning on a Tree[END_REF][START_REF] Frederickson | Nonpreemptive ensemble motion planning on a tree[END_REF] studied both preemptive and non-preemptive versions of the symmetric SCP in the case when the underlying graph is a tree. They proved that the preemptive version is polynomial-time solvable and yielded two exact algorithms. However, the non-preemptive version was shown to be NP-hard and severalapproximations were provided. Kerivin et al. [START_REF] Kerivin | The Eulerian closed walk with precedence path constraints problem[END_REF] first proved that the optimal solutions of the preemptive stacker crane problem can be determined by the simple knowledge of the arc sets related to the vehicle route and to the demand paths. Using this result, they introduced, to the best of our knowledge, the first integer linear model for both symmetric and asymmetric versions of the preemptive SCP [START_REF] Kerivin | Models for the single-vehicle preemptive pickup and delivery problem[END_REF]. This formulation has a polynomial number of variables and an exponential number of constraints. However, the authors showed that the linear relaxation of the formulation can be solved in polynomial time.

Several variants of the pickup and delivery problems closely related to the SCP have been also studied. We should mention the Pickup and Delivery Traveling Salesman Problem (PDTSP) which corresponds to the non-preemptive stacker crane problem where no capacity constraint is taken into account (the vehicle V can contain as many object as required) and where every node of the network G is the origin or the destination of exactly one demand. Rodin and Ruland [START_REF] Ruland | The pickup and delivery problem: Faces and branch-and-cut algorithm[END_REF] presented an integer linear formulation for the PDTSP and used a branch-and-cut algorithm to solve it exactly. A polyhedral study of this formulation and of several other valid constraints was then made by Dumitrescu in [START_REF] Dumitrescu | Polyhedral results for the pickup and delivery travelling salesman problem[END_REF]. The PDTSP has also been well studied from a heuristic point of view and many local search algorithms [START_REF] Healy | A new extension of local search applied to the Dial-A-Ride Problem[END_REF][START_REF] Psaraftis | k-interchange procedures for local search in a precedence constrained routing problem[END_REF][START_REF] Renaud | Perturbation heuristics for the pickup and delivery traveling salesman problem[END_REF][START_REF] Renaud | A heuristic for the pickup and delivery traveling salesman problem[END_REF] have been tested on the PDTSP, which involved a local transformation procedure defined as extensions of the k-interchange procedure defined by Lin [START_REF] Lin | Computer solutions to the traveling salesman problem[END_REF] and Lin et Kernighan [START_REF] Lin | An effective heuristic algorithm for the traveling salesman problem[END_REF] for the TSP. The asymmetric version of the PDTSP was considered by Kalantari et al. [START_REF] Kalantari | An algorithm for the traveling salesman problem with pickup and delivery customers[END_REF] who developed a branch-and-bound algorithm based on Little et al. scheme for the asymmetric TSP [START_REF] Little | An algorithm for the traveling salesman problem[END_REF]. Furthermore, if every node may be incident several demands while the vehicle route is imposed to define a Hamiltonian circuit, one can check that the asymmetric PDTSP is nothing but the Precedence Constrained Asymmetric TSP (PCATSP) also called the Sequential Ordering Problem (SOP). Polyhedral approaches and branch-and-cut algorithms [START_REF] Ascheuer | A cutting plane approach to the sequential ordering problem (with applications to job scheduling in manufacturing)[END_REF][START_REF] Ascheuer | A Branch & Cut algorithm for the Asymmetric Traveling Salesman Problem with Precedence Constraints[END_REF][START_REF] Balas | The precedence constrained asymmetric traveling salesman problem[END_REF][START_REF] Gouveia | On extended formulations for the precedence constrained asymmetric traveling salesman problem[END_REF] as well as heuristics [START_REF] Chen | Commonality and genetic algorithms[END_REF][START_REF] Gambardella | An ant colony system hybridized with a new local search for the sequential ordering problem[END_REF][START_REF] Montemanni | A heuristic manipulation technique for the sequential ordering problem[END_REF] have been devised to solve this problem.

Variants of the pickup and delivery problem with a single vehicle and capacity constraints have also been considered. For instance, Hernández-Pérez and Salazar-González [START_REF] Hernández-Pérez | The multicommodity one-to-one pickup-anddelivery traveling salesman problem[END_REF] considered the SOP with capacity constraints, which they called the multi-commodity one-to-one pickup-and-delivery Traveling Salesman Problem (m-PDTSP). They gave mixed-integer linear formulations which they solved through branch-and-cut algorithms. Furthermore, the method given by Kalantari et al. in [START_REF] Kalantari | An algorithm for the traveling salesman problem with pickup and delivery customers[END_REF] may also be applied for the asymmetric PDTSP with capacity constraints. Kerivin et al. [START_REF] Kerivin | Models for the single-vehicle preemptive pickup and delivery problem[END_REF] extended their model for the preemptive asymmetric stacker crane problem when the vehicle (respectively demands) has a capacity (respectively volume). Some of the previously mentioned problems have also been studied with additional constraints such as time windows [START_REF] Mitrovic-Minic | Pickup and delivery problem with time windows: A survey[END_REF][START_REF] Parragh | A survey on pickup and delivery problems: Part II: Transportation between pickups and delivery locations[END_REF][START_REF] Savelsberg | The general pickup and delivery problem[END_REF], precedence constraints imposed to demand processes [START_REF] Fiala Timlin | Precedence constrained routing and helicopter scheduling[END_REF][START_REF] Fiala Timlin | Precedence constrained routing and helicopter scheduling: heuristic design[END_REF] or LIFO loading policy [START_REF] Cordeau | A branch-and-cut algorithm for the pickup and delivery traveling salesman problem with LIFO loading[END_REF]. Moreover, when the transportation involve human beings, the objective may not only to minimize the total cost of the vehicle route but may also put at stake people dissatisfaction [START_REF] Cordeau | The dial-a-ride problem: models and algorithms[END_REF][START_REF] Psaraftis | A Dynamic Programming Solution to the Single-Vehicle Manyto-Many Immediate Request Dial-a-Ride Problem[END_REF] (which can be expressed using people riding time or difference between desired time and arrival one).

The stacker crane problem can be extended to the case where every demand has several origins and destinations. This extension, named the Swapping Problem, belongs to the class of many-to-many pickup and delivery problems (See [START_REF] Berbeglia | Static pickup and delivery problems: a classification scheme and survey[END_REF] for a classification of routing problems). Usually, symmetric costs are considered. Moreover, it is assumed that every node is the origin of one demand and the destination of another one (a null demand may be considered if necessary). Swapping problems may be preemptive, non-preemptive, symmetric or non symmetric. They were first introduced by Anily and Hassin [START_REF] Anily | The swapping problem[END_REF], which exhibited a 2.5-approximation scheme for some mixed version of the problem. Anily et al. [START_REF] Anily | The preemptive swapping problem on a tree[END_REF] also considered the preemptive swapping problem in the case when the network is a tree. They proved that the problem remains NP-hard and gave in this case a 1.5-approximation. Recently, Bordenave et al. [START_REF] Bordenave | A branch-and-cut algorithm for the preemptive swapping problem[END_REF] presented a branch-and-cut algorithm for the preemptive swapping problem which allowed them to deal with instances with 100 nodes and 8 demands. They also designed a two-phase heuristic [START_REF] Bordenave | Heuristics for the mixed swapping problem[END_REF] for the asymmetric mixed swapping problem and applied to instances with up to 10000 nodes with an average optimality gap which did not exceed 1%.

Though many variants of single-vehicle pickup and delivery problems have been studied, we still may notice that few of them take into account the possibility of reloads in the transportation of the demands. In fact, what is more often considered is the case when reloads correspond to transshipments, that is, when a demand is unloaded from a vehicle in order to be reloaded into another vehicle. Such possibility may appear to pickup and delivery problems involving several vehicles. Despite the fact that these problems are quite different from the SCP, one must mention some of them in order to complete our bibliographic overview: Pickup and Delivery Problem with Transfers (PDPT) [START_REF] Cortés | The Pickup and Delivery Problem with Transfers: Formulation and Solution Approaches[END_REF], the Pickup and Delivery Problem with Time-Windows and Transshipments (PDPTWT) [START_REF] Mitrovic-Minic | The pickup and delivery problem with time windows and transshipment[END_REF] and the Pickup and Delivery Problem with Reloads (RPDP) [START_REF] Oertel | Routing with Reloads[END_REF]. When dealing with these problems, which sometimes involve account time-windows and capacity constraints, authors often prefer the terms transfer or transshipment to the words reloads or preemption.

Mitrovic-Minic and Laporte [START_REF] Mitrovic-Minic | The pickup and delivery problem with time windows and transshipment[END_REF] gave a two-phase heuristic to approximately solve the PDPTWT. They first construct an initial solution using multi-start cheapest insertion procedure, and next improve this solution by successively removing and reinserting every demand, with one or no reload. The experimental results they obtained show that allowing transshipment may be very useful to reduce total travel distance.

In order to get a model for the RPDP, Oertel [START_REF] Oertel | Routing with Reloads[END_REF] defined an auxiliary graph by considering two copies of every origin/destination node. Using this new graph, he first proposed a mixed-integer formulation for the problem, and next designed a tabu search insertion based algorithm, which could efficiently deal with instances with more that seventy demands. Cortés et al. [START_REF] Cortés | The Pickup and Delivery Problem with Transfers: Formulation and Solution Approaches[END_REF] used the same kind of trick and handled their model Benders decomposition.

The focus of this paper will be on the preemptive asymmetric stacker crane problem, which we shall denote by APSCP. We are first going (Section II) by setting our problem in a formal way, while assuming the triangle inequality for the cost function and while doing in such a way that origin/destination pairs of nodes become pairwise disjoint, that network be complete and that the asymmetric costs satisfy the triangle inequalities. Next (Section III) we shall prove some structural results which will allow us to turn our problem into a non constrained tree design problem. Thus reformulation of the problem will lead us to design (Section IV) in a natural way local search heuristic scheme together with a linear integer programming model, which will be implemented and tested in Section V, providing us with rather satisfactory numerical results.

II. A Formal Description of the APSC Problem.

Notations.

For any sequence = {x 1 , .., x n } and any object x = x i in , we denote by Succ( ,x) (Pred( ,x)), the successor (predecessor) x i+1 (x i-1 )of x in , and by Rank( ,x) the rank i of x = x i in . A sequence with only one element x is denoted by {x} and the empty sequence is denoted by Nil. We call subsequence of any sequence ' which may be written {x i1 , .., x ip } with i 1 < i 2 < .. < i p . The first (last) element x 1 (x n ) of is denoted by First( ) (Last( )). The number n of element of is denoted by Length( ). We denote by the concatenation of operator, which takes two sequences = {x 1 , .., x n } and ' = {y 1 , .., y m } and concatenates them into a unique sequence ' = {x 1 , .., x n , y 1 , .., y m }. We denote by * the operator which construct a sequence from its first element x 1 and from its tail subsequence Tail( ) = {x 2 , .., x n }: = First( )* Tail( ) = {x 1 , .., x n }. If x = x i and y = x j are two elements of such that i = Rank( , x) ≤ j = Rank( ,y), then we denote by I( , x, y) the subsequence {x i , .., x j } of  which is defined by all z such that Rank( , x) ≤ Rank( ,z) ≤ Rank( ,y). Any subsequence ' of which may be written I( , x, y) is also called a segment of . We call cut of any decomposition c = ( ', ") of as a concatenation ' " where both 'and " are segments of .

Modelling the Asymmetric Pre-emptive Stacker Crane Problem (APSCP);

At it was told in the introduction, the Asymmetric Pre-emptive Stacker Crane Problem can be roughly described as follows:

-a single vehicle V is required in order to address a set K of transportation demands, while performing some tour inside a given transit network G. Any demand k K is expressed as a pair (o(k), d(k)) of nodes of G, according to the following semantics: o(k) is the origin node of k; d(k) is the destination node of k;

V must transport exactly one load unit L(k) from o(k) to d(k). -the load capacity of V is equal to 1; -V is allowed to address the demands of K in a pre-emptive way: it may, while carrying the load L(k), stop at some node x, unload L(k), deal with other demands and next come back to x, load again L(k), and keep on with the handling of demand k. Such an intermediate load is then called a reload node for the demand k. -V starts and ends its tour in a "Depot" node, and try to do it as fast as possible, in the sense of a cost (length) function which is supposed to be defined on the arcs of the network G.

In order to set a formal model for this problem, we make copies of the original nodes of the network G in such a way that the nodes Depot, o(k), d(k), k K, and the eventual reload nodes becomes all distinct. That means that we consider a node set X which may be written according to the following partition:

X = {Depot} X O X D X R in such a way that: X O = {o(k), k K); X D = {d(k), k K};
X R contains a copy of every element in {Depot} X O X D together with a set of other eventual reload nodes. Then the original cost (length) function which was defined on the arc set of the network G gives rise, through a shortest path computing process, to a X.X indexed (shortest path) distance matrix DIST, such that if x is in {Depot} X O X D and if x' is the copy of x in X R then DIST(x, x') = 0. We suppose of course that DIST satisfies the Triangle Inequality Property, but that it does not need to be symmetric.

Then we define labelled link (L.L), as being any triple r = (x, y, k), where x and y are nodes of X and k is a label in the set {0} K: x (y) is called the starting node (ending node) of the labelled link r and is denoted by Start(r) (End(r)); k is called the label of r and is denoted by Label(r).

It comes that a tour defined on X is going to be a sequence of labelled links. For such a tour , and for any label k in {0} K, we denote by (k) the labelled link sequence which derives in a natural way from by only considering the labelled link r such that Label(r) = k, and we call it the deriving subsequence of G related to the label k. Of course, we understand that the meaning of the label k = Label(r) of some labelled link which appears in a tour G related to the activity of the vehicle V, is that k = 0 means that V is empty when running from x = Start(r) to y = End(r) and that k ≠ 0 means that V is then carrying the unit load L(k).

The cost of is then defined in a natural way as the quantity:

Cost( ) =  r DIST(Start(r), End(r)).
Clearly, not any tour may be viewed as reflecting the activity of a vehicle V which conveniently handles every demand k K. We see that in order to get it, we need to impose G to be valid, which will mean that:

-

III. Some Structural Results.

We are now going to state and prove some results which will be the basis for the design of the heuristics which will be described in Section IV.

III. 1. A Theorem.

Let some valid tour. For any labelled link r = (x, o(k), 0) in , we denote by ( , r) the unique labelled link (d(k), y, 0) which is also in . By the same way, if x is some node in X R , such that a triple a labelled link r = (y, x, k), k ≥ 1, is in , then we also denote by ( , r) the unique triple r' = (x, z, k) which exists in G and which is such that: -Rank( , r') > Rank( , r); -Rank( , r') is minimal with this property. We say that two Labelled links r and r' in are overlapping if we have: Rank( , ( , r')) > Rank( , ( , r)) > Rank( , r') > Rank( , r).

Then we may state:

Theorem 1.
Let be some optimal tour for the APSCP Problem, which we suppose chosen in such a way that: -(A): Length( ) is the smallest possible; -(B): the number of labelled links r in which are such that Label(r) ≠ 0 is the smallest possible, (A) being supposed to be satisfied; Then, the following assertions must be true:

-(S1): does not contain two occurrences of the same labelled link r = (x, y, k), with k ≠ 0; -(S2): does not contain two consecutive labelled links r and r' such that Label(r) = Label(r'); -(S3): does not contain two overlapping labelled links r and r'; -(S4): does not contain two labelled links r and r' such that End(r) = End(r').

Proof.

We assume that is given, which is an optimal solution of APSCP and which is such that (A) and (B) are true. Part (S1). If r = (x, y, k), k > 0 appears twice in , with respectively rank s and s', then x and y are both reload nodes, and we may replace, in any labelled link r" which is such that: -s ≤ Rank( , r") < s'; -Label( r") = k; the label value Label(r") by 0. While doing it, we keep on with a valid tour which is an optimal solution of APSCP and we get a contradiction on the (B) hypothesis.

Part (S2).

If r = (x, y, k) and r' = (y, z, k) were two consecutive labelled links of such that Label(r) = Label(r'), then we would be able to remove both r and r' from , and replace them by a unique labelled link (x, z, k). While doing this, we would also keep, because of the Triangle property on the DIST matrix, an optimal solution of APSCP, and this solution would contradict the (A) hypothesis.

Part (S3).

Let us suppose that there exists two labelled links r = (y, x, k) and r' = (y', x', k') which are overlapping in . We may choose them in such a way that Rank( , ( , r')) -Rank( , r) is the smallest possible.

(E1) Because of (E1), we see that x and x' must be reload nodes: if, for instance, x were an origin node o(h), then we would have k' ≠ h and there would exist an labelled link r" such that:

-Label(r") = h; -Rank(G, r) < Rank(G, r") < Rank(G, r') < Rank(G, ( , r")) < Rank(G, ( , r)) < Rank(G, ( , r')). Then we might deduce a new overlapping pair (r", r') which would induce a contradiction on the minimality assumption (E1). By the same way, we may check that x' cannot be an origin node o(h). Thus x and x' are both reloads nodes, and we clearly have: k ≠ k' ≠ 0. It comes that we may write:

-r = (y, x, k), r' = (y', x', k') ; -( , r) = (x, z, k), ( , r') = (x', z', k'). So, we set : 5 . Of course, the lengths of and are Aux equal, as well as their respective costs. So we state:

Lemma 1.

Aux defined above is a valid tour.

Proof-Lemma.

We only need to check that switching 2 , 3 and 4 does not break any sequence (k), or, in other words, that for any k ≠ 0, we have (k) = Aux (k). If the converse were true, we would be able to find k" ≠ k', k, k" ≠ 0, as well as two labelled links r" and ( , r"), with label k" or with the ending node of r" equal to o(k"), such that one of the three following relations would be true:

-r" 2 and ( , r") 3 ;

(E2); -r" 2 and ( , r") 4 ;

(E3); -r" 3 and ( , r") 4 ;

(E4); In case (E2) or (E3) were true, r" and r' would be overlapping, and would contradict the (E1) hypothesis, related to the minimality of Rank( , ( , r')) -Rank( , r) In case (E4) were true, r and r" would be overlapping, and would contradict the (E1) hypothesis, related to the minimality of Rank( , ( , r')) -Rank( , r). In any case, we become able to conclude. END-LEMMA.

The above lemma allows us to conclude the proof of (S3) by noticing that r and ( , r) become consecutive in the valid tour Aux , which implies (proof of statement (S2)) that r and ( , r) may be replaced in Aux by a unique labelled link (y, z, k) in such a way that Cost( Aux ) does not increase and that Length( Aux ) decreases, inducing a contradiction on the (A) hypothesis.

Part (S4).

Let us suppose that contains two labelled links r and r' such that End(r) = End(r') and such that Rank( , r) < Rank( , r'). Since the starting node x of r cannot be in {Depot} X D , it must be a reload node in X R which is used twice as a reload node. So, r, ( , r) and r' and ( , r') may be written:

-r = (y, x, k), k ≠ 0; -r' = (y', x, k'), k' ≠ 0, k; -( , r) = (x, z, k); -( , r') = (x, z', k'). Because of (S3) we must have:

-Rank( , r) < Rank( , ( , r)) < Rank( , r') < Rank( , ( , r')) (E5) or -Rank( , r) < Rank( , r') < Rank( , ( , r')) < Rank( , ( , r)).

(E6) Let us first suppose that (E5) holds. Then we set: 5 . Of course, the lengths of and are Aux equal, as well as their respective costs, and we proceed as in the proof of (S3) in order to prove that Aux must be a valid tour. But we also notice, as in the proof of (S3), that Aux can be shortened by replacing the consecutive labelled links r and r* by a unique labelled link (y, z, k), in such a way that Cost( Aux ) does not increase and that Length( Aux ) decreases, inducing a contradiction on the (A) hypothesis.

We apply exactly the same kind of reasoning in case (E6) holds. END-THEOREM.

III. 2. A Tree Representation of the APSCP Problem.

Theorem 1 leads us to introduce the following definition:

Strongly Valid Tour: a valid tour is a strongly valid tour if it satisfies the (S1)…(S4) properties which are listed into the statement of Theorem 1.

The following figure 2 shows us how the valid tour of Figure 1 may be turned into a strongly valid tour whose cost ' is no more than the cost of : Clearly, solving the APSCP Problem means finding a strongly valid tour with minimal cost value. Now, we are going to see that any strongly valid tour may be represented as some kind of tree, and this will provide us with the basis (section IV) for the algorithms which we are going to design in order to deal with APSCP.

Bipartite Ordered Trees: we say that a tree T is a bipartite ordered tree if: its nodes can be split into two classes A and B in such way that nodes in class A have their sons in class B and conversely;

for every node x in T which is not a terminal node (leaf) the son set (T, x) associated with x is linearly ordered: thus (T, x) is described as a sequence.

We say that a bipartite ordered tree T is consistent from the APSCP instance defined by the demand set K and by the node set X if:

-the nodes in T can be identified with the demands k K (we shall then talk about demand nodes) or with nodes in {Depot} X R , (and then we talk about reload nodes) and any possible demand node k K appears in T, while only some of nodes of {Depot} X R appear in T: those nodes in {Depot} X R define the active reload node set ACTIVE(T) of T;

(S5) -The root of T is the Depot node and the terminal nodes (leafs) of T must belong to the demand node set; (S6) -For any demand node k, its linearly ordered son set RELOAD(T, k) in T is made with active reload nodes and its father FATHER(T, k) is in ACTIVE(T); (S7) -For any reload node x, its linearly ordered son set DEMAND(T, x) in T is made with demand nodes and its father FATHER(T, x) is in K. (S8)

For such a bipartite ordered tree T, we may define a cost value Tree-Cost as follows:

- Theorem 2 : There is a one-to-one correspondence Tree between the strongly valid tours and the bipartite ordered tour which are consistent with X and K, which is such that, for any strongly valid tour , we have: Tree-Cost(Tree( )) = Cost( ).

Proof.

We first consider a strongly valid tour and perform the following construction, which make us get Tree( ) from :

-ACTIVE(T( )) is defined as the set of the nodes of {Depot} X R which appear in some labelled link of , and which are then said to be active reload nodes for ; -For any demand node k K, the son set Reload(T( ), k) is made with the reload nodes which appear in some labelled link of (k), ordered according to their appearance order in (k); -For any active reload node x in {Depot} X R , we denote by (x) = (x, y, 0) and by (x) = (z, x, 0) the two labelled links with label 0 which involve x in and which are such that: Rank( , (x)) < Rank( , (x)). Then we define the son set Demand(T( ), x) by setting that a demand k K is a son of x if the unique labelled link r(k) = (o(k), t, k) which appears in is such that: Rank( , (x)) < Rank( , r(k)) < Rank( , (x)) there exists no reload node y such that Rank( , (x)) < Rank( , (y)) < Rank( , r(k)) < Rank( , (y)) < Rank( , (x)).

Figure 3:

The bipartite tree Tree( ') which derives from the strongly valid tout ' of Figure 2.

Then it comes next that checking that the so defined Tree correspondence is as it is claimed in the statement of Theorem 2 is purely routine. END-THEOREM.

We deduce:

Corollary 1: Solving a APSCP instance (X, DIST, K) means finding a bipartite ordered tree T consistent with (X and K) such that Tree-Cost(T) is the smallest possible.

The interest of this last statement is clearly that it provides us with a bipartite tree formulation of the APSCP problem which is far less constrained that the original one.

III.3. A related Integer Linear Programming formulation of APSCP.

This ILP formulation is going to allow us to compare, in the case of small instances, the results obtained through the heuristic methods which will be described in Section IV with exact results. In order to get it, we first need to proceed to the following construction:

-An auxiliary network G = (X*, E). We first consider a copy X* R of the reload node set X R and a copy Depot* of the Depot node, and we set: X* = X X* R {Depot*}; For any node x in X R , we denote by x* its copy in X* R . Also, for any origin node x = o(k) in X O , we denote by x* the related node d(k). Then we define, on the node set X*, the arc set E as follows:

-

E = {(Depot, x), x X O } {(x, Depot*), x X D } {(o(k), d(k)), k K} {(d(k'), o(k)), k ≠ k' K} {(x, y), (y, x), x X O , y X R } {(x, y) (y, x), x X* R , y X D } {(x, y), x X* R , y X R }.
-every arc e in E is then provided with a length DIST*(e) which derives from the DIST distance matrix in a natural way.

Let us recall that a path of a such a network G is a node sequence such that, for any node x in , the pair (x, Succ( , x)) defines an arc of E. One easily checks that any strongly valid tour can be turned into a path * of the network G, in such a way that: -(S9): * starts from Depot and ends into Depot* and * is an elementary path, i.e, its visits any node at most once; -(S10): for every k in K, * visits o(k) and d(k), according to this order and for every x in X R , * visits x if and only if it visits x*, and, in case it does it, it does it according to this order; -(S11): for any pair x, y, x ≠ y, in X R X O the following implication is true: Rank( *, x) < Rank( *, y) and Rank( *, y) < Rank( *, x*) => Rank( *, y*) < Rank( *, x*). This condition is called the non overlapping condition.

-(S12): Cost( ) = x

x ≠ Depot* DIST*(x, Succ( *, x)) = Length of * for the DIST* length function. A path of the network G which satisfies (S9)…(S12) above will be said to be a strongly valid path.

Theorem 3.

For any strongly valid path , there exists a strongly valid tour such that * = .

Proof.

Let us first describe in an accurate the way is going to derives from . It will occur through the following RECONSTRUCT procedure: The instruction (I1) works here because of (S10) above, and because an arc of G which arrives on (x) must come from an origin node o(k) or from a node z in X* R . In this last case, a simple induction reasoning makes appear the fact that k is different from 0. We get our result while proceeding by induction on the length (the number of nodes) of . In case involves no node in X R , then the results come in a trivial way. Else, we consider x 0 X R which is the first node of X R which appears in . We notice that Pred( , x 0 ) must be some node o(k), k K. Thus the arc (Pred( , x 0 ), Succ( ,x 0 *)) belongs to the arc set E, and the removal of the subpath I( , x 0 , x 0 *) from provides us with an other path 1 of the graph G. Let us set: K 1 = { k K such that o(k) and d(k) are nodes of 1 }; K 2 = { k K such that o(k) and d(k) are nodes of 2 = I( , x 0 , x 0 *)}; K 1 and K 2 define a partition of K, and one sees that 2 can be viewed as a strongly valid path, if we restrict ourselves to K 2 as a demand set and if we consider that x 0 and x 0 * play the role of Depot and Depot*. Thus it comes from the induction hypothesis that it may be written, under this restriction, according to the form 2 = * 2 . By the same way, 1 is also a strongly valid path if we restrict the demand node set to K 1 , and it comes from the induction hypothesis that it may be written, under this restriction, according to the form 1 = * 1 . We only need to insert 2 between Pred( , x 0 ), x 0 , k) and (x 0 , Succ( , x 0 *), k) in 1 , in order to get such that = *. END-THEOREM.

RECONSTRUCT

Corollary 2: Solving a APSCP instance (X, DIST, K) means finding a strongly valid path with minimal length (for the DIST* length function) value in the network G.

This corollary allows us to set an Integer Linear Programming model as follows:

This model involves a {0, 1} vector flow z = (z e , e E), a Rank integral vector R = (R x , x X), as well as a positional {0, 1}vector t, which is indexed on the pairs (x, y), x ≠ y, x, y X R X O , with the following semantics:

-for any arc e in E, z e = 1 iff the arc e is in the strongly valid path ;

-for any node x in , R x will provide us with the rank of x in ; -for any pair (x, y), x ≠ y, x, y X R X O : t x,y = 1 iff Rank( , y) < Rank( , x*); t x,y = 0 iff Rank( , x*) < Rank( , y);

Then the translation of (S9)…(S12) into Integer Linear Programming constraints provides us with the following linear integer program:

APSCP Integer Linear Programming Formulation.

Unknown vectors.

z = (z e , e E), with values in {0,1}; R = (R x , x X) Integral and ≥0; t = (t x,y , x ≠ y, x, y X R X O ) with values in {0,1}.

Performance Criterion.

Minimize e E DIST*(e).z e (translation of (S12)) Constraints.

-z is a flow vector, which satisfies the usual Kirshoff law in any node but in Depot and Depot*; -the inflow induced by z in Depot (Depot*) is equal to 0 (1), while the related outflow is equal to 1 (0); (translation of (S9)) -in any node of X D X O , the inflow induced by z is equal to 1; (translation of (S10)) -in any node of X R X* R , the inflow induced by z is at most equal to 1, and the inflow value in x is equal to the inflow value in x*; (translation of (S10) and (S9)) -for any k in K, we have R o(k) ≤ R d(k) -1; (translation of (S10)) -for any x in X R , we have R x ≤ R x* -1; (translation of (S10)) -for any arc e = (x,y) in E, we have z e + (R x + 1 -R y )/Card(X*) ≤1; (translation of the implication z e = 1 -> (R x + 1 -R y ) ≤0) -for any pair (x, y), x ≠ y, x, y X R X O : (translation of the non overlapping condition) t x,y + (R y + 1 -R x* )/Card(X*) ≤1 ; t x,y + (R y -1 -R x* )/Card(X*) 0 ; t x,y + (R y* + 1 -R x* )/Card(X*) ≤1.

Corollary 2: Solving a APSCP instance (X, DIST, K) means solving the above Integer Linear Programming model.

IV. Tree Based Heuristics for the APSCP Problem.

The algorithms which we are going to describe here, and which will be tested in the next Section, derive in a straightforward way from the tree representation of the APSCP Problem which we got in Section III.2. These algorithms are very simple greedy insertion algorithms and descent algorithms, based upon the use of 2 classes of operators:

Insertion Operators: these operators act on some bipartite ordered tree T consistent with the node set X and with a subset K' of the demand set K, and insert some demand k K -K' in T. We use two operators:

INSERT-SIMPLE: its parameters are some active reload node x in {Depot} X R , and some cut (l 1 , l 2 ) of the sequence DEMANDE(T, x) = l 1 l 2 . It acts by inserting the segment {k} into this cut: DEMANDE(T, x) <-l 1 {k} l 2 .

INSERT-with-RELOAD: its parameters are some demand node k' in K', a cut c = (l 1 , l 2 ) of the sequence RELOAD(T, k'), and a non active reload node x. It acts by:

-inserting the segment {x} into the cut c: RELOAD(T, x) <-l 1 {x} l 2 ; -by making x be active and setting: DEMAND(T, x) <-{k}; RELOAD(T, k) <-Nil.

Local Transformation Operators: these operators act through side effect on some bipartite ordered tree T consistent with X and K, and they modify T. We use 6 operators:

-MOVE-RELOAD: its parameters are some active reload node x and some non active reload node y. It replaces x by y in T. -MOVE-RELOADS: its parameters are two different demand nodes k and k', a segment l of RELOAD(T, k) and a cut c = (l 1 , l 2 ) of RELOAD(T, k'). It removes l from RELOAD(T, k) and it inserts it into the cut c. Its precondition is that k does not dominate k' in the tree T, i.e, that k cannot be obtained from k' through a succession of applications of the FATHER operator. -MOVE-RELOADS1: its parameters are some demand node k, some segment l of RELOAD(T, k) which induces a decomposition RELOAD(T, k) = l 3 l l 4 , and a cut c = (l 1 , l 2 ) of l 3 l 4 . It first remove l and next insert it into the cut c: RELOAD(T, k') <-l 1 l l 2 . -MOVE-DEMANDS: its parameters are two different active reload nodes x and y, a segment l of DEMAND(T, x), and a cut c = (l 1 , l 2 ) of DEMAND(T, x'). It removes l from DEMAND(T, x) and it inserts it into the cut c. In case l 1 = l 2 = Nil, it remove the reload node x from T, which becomes non active. Its precondition is that x does not dominate y in the tree T. -MOVE-DEMAND1: its parameters are some reload node x, some segment l of DEMAND(T,

x) which induces a decomposition DEMAND(T, x) = l 3 l l 4 , and a cut c = (l 1 , l 2 ) of l 3 l 4 . It first remove l and next insert it and it into the cut c: DEMAND(T, x) <-l 1 l l 2 . -MOVE-DEMANDS-RELOAD: its parameters are an active reload node x, a non active reload node y, a demand node k, a segment l of DEMAND(T, x) and a cut c = (l 1 , l 2 ) of RELOAD(T, k). It first turns y into an active reload node, next removes l from DEMAND(T, x) and inserts it into DEMAND(T, y), and ends in inserting the segment {y}into the cut c. In case l = DEMAND(T, x), it turns x into a non active reload node. Its precondition is that k is dominated by not demand node k' in l.

Then we can propose a first insertion greedy algorithm for dealing with APSCP:

Algorithm APSCP-INSERTION:

Randomly define a linear ordering on the elements of K; T = {the tree reduced to the root node Depot}; For k K, K being scanned according to the linear order do Compute the insertion operator I (among INSERT-SIMPLE and INSERT-with-RELOAD) and the related parameter u (u = (x, (l 1 , l 2 )) in case I = INSERT-SIMPLE, u = (k', (l 1 , l 2 ), x) in case I = INSERT-with-RELOAD), such that the insertion of k through I(u) induces the smallest possible increase of Tree-Cost(T); Apply I(u) to T;

Filtering the search for the good value of the parameter u.

In case I = INSERT-SIMPLE, the related optimal value of u can be obtained in a very fast way through exhaustive scanning of the sequences DEMAND(T, x) for all the active reload nodes x. In case I = INSERT-with-RELOAD, one must deal with the search for the new reload node x, which may be time consuming in case X R is large. In order to avoid spending too much time while trying all the possible reload nodes x, we try to identify in a fast way those nodes which are likely to provide us with an efficient insertion by proceeding as follows:

-for every node x in X R , we keep in memory a set N(x) of neighbours of x, that means of nodes y which are such that DIST(x, y) ≤ R, where R is some threshold which is chosen in such a way that the induced neighbour graph be connected, and that the cardinality of any set N(x) remains small enough.

-by the same way, we keep in memory, for every pair of reload nodes (x, y), what we call the middle of x and y, that means some node z = MID(x, y) which is such that the difference between DIST(x, z) + DIST(z, y) -DIST(x, y) remains small, both quantities DIST(x, z) and DIST(y, z) being close to each other; -if we denote by y the last reload node in l 1 and by z the first reload node in l 2 , then we see that we should try to select x in such a way that DIST(y, Filtering the search for the good value of the parameter vectors. In the case of the (1) above instruction, the active nodes x are scanned in an exhaustive way, but the search for the y value is restricted to the neighbourhood set N(x).

In the case of the instruction (2) and (3), the threshold value H is involved as followed:

-the segment l is tried only if the decomposition RELOAD(T, k) = l 3 l l 4 is such that DIST(Last(l 3 ), First (l)) and DIST (Last(l), First(l 4 )) ≥ H; -the cut c = (l 1 , l 2 ) is tried only if DIST(Last(l 1 ), First (l 2 )) ≥ H.

In the case of the instruction (4) and ( 5), the threshold value H is involved as followed:

-the segment l is tried only if the decomposition DEMAND(T, x) = l 3 l l 4 is such that DIST(d(Last(l 3 )), o(First (l))) and DIST (d(Last(l)), o(First(l 4 ))) ≥ H; -the cut c = (l 1 , l 2 ) is tried only if DIST(d(Last(l 1 )), o(First (l 2 ))) ≥ H. In the case of the (6) instruction, the threshold value H is involved as in ( 4) and ( 5), and the search for y is performed, once x, k, l, c = (l 1 , l 2 ) have been determined inside the neighbour set N(t) of t = MID(MID(Last(l 1 ), o(First(l))), MID(d(last(l)), First(l 2 ))).

Remark: of course, it would be possible to improve the performance of the APSCP-DESCENT algorithm by casting it into a scheme like the Simulated Annealing scheme or the Tabu List scheme. But it is not really the purpose of the paper: as we shall see in Section V, our tree representation of the APSCP problem, together with the operators which we just described above, are sufficiently powerful to provide us, under small computing costs, with very good solutions for our APSCP problem.

V. Experiments.

We have been performing experiments, on PC IntelXeonwith 1.86 GHz, 3.25 Go Ram, while using a Visual Studio C++ compiler, and while focusing on several points:

-the ability of APSCP-INSERTION and APSCP-DESCENT to get solutions close to the optimal theoretical solutions; -the running time of those algorithms; -the characteristics of the solutions: number of reload nodes which appear in the solution, improvement of the Cost-Tree value induced by pre-emption; -the impact of the different local transformation operators on the behaviour of the algorithms.

In order to do this, we performed several tests, while using node sets X and distance matrices DIST proposed by the TSPLIB libraries, and by selecting origin/destination pairs (o(k), d(k), k K) in a random way inside the set X. We dealt with instances which involves from 20 to 300 nodes, and from 10 to 100 origin destination pairs, and, in case of small instances, got exact results through the use of the LIP formulation of Section III.3, augmented with cutting planes techniques (see [START_REF] Kerivin | Models for the single-vehicle preemptive pickup and delivery problem[END_REF][START_REF] Kerivin | The Eulerian closed walk with precedence path constraints problem[END_REF]). The results which we get may be summarized as follows:

Comments:

We see that the results which we get through application of a combination of a greedy scheme and a Monte-Carlo diversification scheme are most often rather goods, in the sense that its allows us to get in a fast way solutions which are not too far from the best theoretical ones. Still, we also notice that our greedy scheme is in trouble when it comes to pre-emption handling, that means when it comes to creating reload nodes. Thus, the results which derive from the application of APSCP-INSERTION get worse as soon as there is an increase of the gap between the optimal preemptive optimal value and the non pre-emptive one. The results which we get may be summarized as follows:

  for any consecutive pair of labelled links r, r'= Succ(G, r) in G, we have End(r) = Start(r'); -Start(First( )) = End(Last( )) = Depot; -Any node x in X O X D is involved in exactly two labelled links r and r' = Succ( , r): this means that V moves to o(k) (d(k)), k K, only when it comes to start (finish) dealing with demand k; -The Depot node is involved only in both labelled links r = First( ) and r' = Last( ); -For any demand k K, the deriving subsequence (k) related to k is such that: Start(First( (k)) = o(k); End(Last(d(k)) = d(k); For any consecutive labelled link pair r, r' = Succ( (k), r'), we have End(r) = Start(r'). The following figure 1 provides us with a visualization of a valid tour = {(Depot, o 1 , 0), (o 1 , x, 1), (x, o 2 , 0), (o 2 , y, 2), (y, o 3 , 0), (o 3 , x, 3), (x, d 1 , 1), (d 1 , y, 0), (y, d 2 , 2), (d 2 , 0, x), (x, d 3 , 3), (d 3 , Depot, 0)}.

Figure 1 :

 1 Figure 1: Visualizing a valid tour . All these definitions allow us to formally set our APSCP (Asymmetric Pre-emptive Stacker Crane Problem) Problem as follows:APSCP Problem:

- 1 =

 1 I( , First( ), r) ; -2 = I( , Succ( ,r), r') ; -3 = I( , Succ( ,r'), Pred( , ( , r))) ; -4 = I( , ( , r), Pred( , ( , r'))) ; -5 = I( , ( , r'), Last( )),a and we replace by the concatenation Aux = 1 4 3 2

- 1 =

 1 I( , First( ), r) ; -2 = I( , Succ( ,r), Pred( , ( , r))) ; -3 = I( , ( , r), r') ; -4 = I( ,Succ( ,r'), Pred( , ( , r'))) ; -5 = I( , ( , r'), Last( )), and we replace by the concatenation Aux =

Figure 2 :

 2 Figure 2: A derivation of the valid tour of figure 1 into a strongly valid tour '.

  for any demand node k K: we set: If k is not a terminal node then Cost-Demand(T, k) = DIST(o(k), First(Reload(T,k))) + DIST(Last(Reload(T, k)), d(k)) + x Reload(T, k), x ≠ Last(Reload(T,k)) DIST(x, Succ(Reload(T, k), x))) else Cost-Demand(T, k) = DIST(o(k), d(k)) -for any reload node x {Depot} X R , we set: Cost-Reload(T, k) = DIST(x, o(First(Demand(T,x)))) + DIST(d(Last(Demand(T, x))), x) + k Demand(T, x), k ≠ Last(Demand(T,x)) DIST(d(k), o(Succ(Reload(T, x), k))); -Tree-Cost(T) = Cost( ) = k K Cost-Demand(T, k) + x ACTIVE(T) Cost-Reload(T, x).

  Procedure x <-Depot; <-Nil; While x <> Depot* do y <-Succ( , x); If y may be written y = z*, with z {Depot} X R , then we set (y) = z, else we set (y) = y; If (multi-case branching instruction) 1. x = Depot then <-(x, y, 0)* ; 2. x = o(k), k K then <-(x, y, k)* ; 3. x X R then <-(x, (y), 0)* ; 4. x X D then <-(x, (y), 0)* ; 5. x X* R then <-(x, y, k)* , where k ≠ 0 is such that the labelled link (Pred( , (x)), (x), k) is already in ; (I1)

  x) + DIST(x, o(k)) + DIST(d(k), x) + DIST(x, z) be the smallest possible. Instead of trying all the possible nodes of X R we do it by selecting t = MID(MID(y, o(k)), MID(d(k), z)) and by trying all the nodes x in N(t). Of course, the algorithm APSCP-INSERTION may be used inside a Monte-Carlo Scheme as follows: Parameter : For i = 1 to run the APSCP-INSERTION Procedure; Keep the best result. This greedy insertion algorithm may now used in order to initialize the following APSCP-DESCENT descent algorithm: Algorithm APSCP-DESCENT: Initialize the tree T through APSCP-INSERTION; Initialize the filtering threshold value H; (1): Search (in a filtered way) parameter values x, y for the MOVE-RELOAD operator in such a way that applying MOVE-RELOAD(x, y) to T improves the Tree-Cost quantity; If Success then Go To (1); (2): Search (in a filtered way) parameter values k, k', l, c for the MOVE-RELOADS operator in such a way that applying MOVE-RELOADS(k, k', l, c) to T improves the Tree-Cost quantity; If Success then Go To (1); (3): Search (in a filtered way) parameter values k, l, c for the MOVE-RELOADS1 operator in such a way that applying MOVE-RELOADS1(k, l, c) to T improves the Tree-Cost quantity; If Success then Go To (1); (4): Search (in a filtered way) parameter values x, y, l, c for the MOVE-DEMANDS operator in such a way that applying MOVE-DEMANDS(x, y, l, c) to T improves the Tree-Cost quantity; If Success then Go To (1); (5): Search (in a filtered way) parameter values x, l, c for the MOVE-DEMANDS1 operator in such a way that applying MOVE-DEMANDS1(x, l, c) to T improves the Tree-Cost quantity; If Success then Go To (1); (6): Search (in a filtered way) parameter values x, y, k, l, c for the MOVE-DEMANDS-RELOAD operator in such a way that applying MOVE-DEMANDS-RELOAD(x, y, k, l, c) to T improves the Tree-Cost quantity; If Success then Go To (1); (7): If H is small enough then Stop else set H = H/2 and go to (1).

  Our first experiment is related to the procedure APSCP-INSERTION: we run the APSCP-INSERTION Monte-Carlo scheme with = 100, and we keep memory, for every test with name INST, of the following quantities: REF: Optimal theoretical Tree-cost value; MIN (MAX): Minimal (Maximal) Tree-cost value obtained through iterations of APSCP-INSERTION; MEAN: Maximal (worse) Tree-cost value obtained; EGI: Gap (in %) between REF and the global solution produced by the APSCP-INSERTION Monte Carlo scheme (MIN value). REL: Mean number of reload nodes involved in a solution produced by APSCP-INSERTION; DE/REL: Mean number of demands related to every reload node x (length of the list DEMAND(x)), for x different from Depot; CPU: CPU Mean Time (in milliseconds) for any iteration of APSCP-GREEDY-INSERTION.

  Our second experiment is related to APSCP-DESCENT. We run APSCP-DESCENT from a solution provided by only 1 application of APSCP-INSERTION, and we keep memory, for every test with name INST, of the following quantities: REF: Optimal theoretical value; VAL: The cost value obtained after application of APSCP-DESCENT; EGI: Gap (in %) between REF and the solution produced by APSCP-INSERTION; ED: Gap between REF and the solution produced by APSCP-DESCENT; ETDS: Part of the gap between EGI and ED which is induced by the operators MOVE-DEMANDS and MOVE-DEMANDS1 and MOVE-DEMAND-with-RELOAD; REL: Number of reloads involved in the solution produced by APSCP-DESCENT; TNB: Number of times a local transformation operator is effectively applied inside the APSCP-DESCENT process; CPU: Cpu Running Time (in milliseconds).

Table 1 :

 1 Tests

	INST	REF	MIN	MAX	MEAN	EGI	REL	DE/REL	CPU

performed on 10 instances which we got from the gr24 instance with 49 nodes of the TSPLIB library by randomly sorting 12 demands, and no reload nodes which is not the depot node, the copy of an origin node o(k) or the copy of a destination node d(k), k 1..12.

Table 3 :

 3 Tests performed on 10 instances which we got from the gr120 instance with 241 nodes of the TSPLIB library by randomly sorting 60 demands, and no reload nodes which is not the depot node,

	INST	REF	MIN	MAX	MEAN	EGI	REL	DE/REL	CPU
	Gr_120_v01	Unknown	327070	344139	335480	*	0.42	5.13	78
	Gr_120_v02	Unknown	332679	353240	343622	*	0.64	4.60	78
	Gr_120_v03	Unknown	326189	345859	335567	*	0.27	8.95	78
	Gr_120_v04	Unknown	366750	380409	374350	*	0.36	5.74	78
	Gr_120_v05	Unknown	310179	330002	320326	*	1.53	3.05	78
	Gr_120_v06	Unknown	319619	337921	327146	*	0.64	3.24	78
	Gr_120_v07	Unknown	285989	302602	293394	*	0.70	4.46	78
	Gr_120_v08	Unknown	345100	365379	355230	*	0.56	7.33	78
	Gr_120_v09	Unknown	333909	350171	342549	*	0.47	11.11	78
	Gr_120_v10	318099	337259	354809	345393	6.02	0.36	9.19	78

the copy of an origin node o(k) or the copy of a destination node d(k), k 1..60.

Table 4 :

 4 Tests performed on 10 instances RELn, , which we built in such a way that:

	-the related instance involves n nodes, p = (n/2 -1)/3 -1 reload nodes, (n -p -1)/2 demands;
	INST	REF	MIN	MAX	MEAN	EGI	REL	DE/REL	CPU

-its optimal value is the sum k K DIST(o(k), d(k)), and that the related optimal solution involves all the reload nodes.

Table 5 :

 5 Tests performed on 10 instances which we got from the gr24 instance with 49 nodes of the TSPLIB library by randomly sorting 12 demands, and no reload nodes which is not the depot node, the copy of an origin node o(k) or the copy of a destination node d(k), k 1..12.

	INST	REF	VAL	EGI	ED	ETDS	REL	TNB	CPU
	Gr24_v01	24654	24654	10.3	0	10.3	1	12	93
	Gr24_v02	21395	21914	5.0	2.4	2.6	1	3	47
	Gr24_v03	22834	22874	9.9	0.1	9.8	1	9	46
	Gr24_v04	23255	23435	3.9	0.7	3.2	2	4	31
	Gr24_v04	23993	23993	0	0	0	0	0	31
	Gr24_v06	23233	23763	8.3	2.2	6.1	0	3	31
	Gr24_v07	20224	20244	1.2	0.1	0.2	1	3	78
	Gr24_v08	20865	21084	4.9	1.0	3.9	1	6	31
	Gr24_v09	23054	23054	4.9	0	4.9	1	7	47
	Gr24_v10	26704	26754	4.9	0.2	3.7	1	7	63

Table 6 :

 6 Tests performed on 10 instances which we got from the hk48 instance with 97 nodes of the TSPLIB library by randomly sorting 24 demands, and no reload nodes which is not the depot node, the copy of an origin node o(k) or the copy of a destination node d(k), k 1..24.

	INST	REF	VAL	EGI	ED	ETDS	REL	TNB	CPU (s)

Table 7 :

 7 Tests performed on 10 instances which we got from the gr120 instance with 241 nodes of the TSPLIB library by randomly sorting 60 demands, and no reload nodes which is not the depot node, the copy of an origin node o(k) or the copy of a destination node d(k), k 1..60.

	INST	REF	VAL	EGI	ED	ETDS	REL	TNB	CPU (s)
	Gr_120_v01	Unknown	313315	*	*	5.2	6	89	130
	Gr_120_v02	Unknown	318870	*	*	8.9	1	116	189
	Gr_120_v03	Unknown	314493	*	*	6.2	4	96	125
	Gr_120_v04	Unknown	360902	*	*	3.5	3	107	145
	Gr_120_v05	Unknown	293613	*	*	9.8	4	99	125
	Gr_120_v06	Unknown	303441	*	*	7.7	2	90	80
	Gr_120_v07	Unknown	265713	*	*	6.4	4	84	86
	Gr_120_v08	Unknown	330823	*	*	6.9	4	77	94
	Gr_120_v09	Unknown	314073	*	*	8.5	4	105	70
	Gr_120_v10	318099	318913	9.5	0.25	9.1	4	96	103

Table 8 :

 8 Tests performed on 10 instances RELn, , which we built in such a way that:

	INST	REF	VAL	EGI	ED	ETDS	REL	TNB	CPU
	REL31	11843	11953	16.7	0.9	15.8	4	10	0.04
	REL55	17464	17583	27	0.7	36.3	8	28	0.34
	REL79	21053	21053	29.9	0	29.9	12	56	1.5
	REL91	22323	22605	20.8	1.2	19.6	14	46	2.9
	REL115	24142	24225	28.4	0.34	28.4	18	78	8.4
	REL163	26036	26279	19.3	0.93	18.3	26	79	25
	REL187	26494	26626	26.6	0.49	26.2	30	112	39
	REL211	26775	26851	33.1	0.28	30.1	34	134	78
	REL259	27042	27067	37	0.09	36.9	42	162	232
	REL283	27098	27211	33.1	0.4	32.4	46	202	328

-the related instance involves n nodes, p = (n/2 -1)/3 -1 reload nodes, (n -p -1)/2 demands; -its optimal value is the sum k K DIST(o(k), d(k)), and that the related optimal solution involves all the reload nodes.


Christophe Duhamel
email: christophe.duhamel@isima.fr

Andréa Cynthia Santos

Dario José Aloise

Multicommodity formulations for the prize collecting vehicle routing problem in the petrol industry

Keywords: Vehicle routing problem, prize-collecting, selective traveling salesman, multiobjective

The Mobile Oil Recovery (MOR) unit is a truck designed to pump marginal oil wells in a petrol field. The MOR optimization Problem (MORP) consists in optimizing both the oil extraction and the travel costs. In this article, we describe several formulations for the MORP using a single vehicle and we propose two formulations to the case where several vehicles are used. We strengthen the proposed formulations by taking advantage of the MORP characteristics, by improving the number of subtour elimination constraints and by using cuts. Computational results are presented for instances close to the reality and optimality is proved for instances with up to 200 nodes.

Résumé

Les unités mobiles de pompage sont des camions munis d'un système d'extraction de pétrole. Ils sont utilisés pour les puits ayant une production marginale. Le problème d'optimisation des unités mobiles de pompage consiste à maximiser la quantité totale de pétrole extrait en définissant des tournées quotidiennes tout en minimisant le temps total de parcours. Dans ce travail, nous décrivons plusieurs formulations pour ce problème avec un véhicule unique et nous proposons deux formulations pour le problème avec une flotte de véhicules. Les formulations sont renforcées par l'introduction d'inégalités valides, mais aussi par l'utilisation des caractéristiques du problème pour réduire la quantité de contraintes d'élimination de sous-cycles. Des résultats numériques sont présentés pour des instances proches de la réalité. L'optimalité est prouvée pour des instances avec 200 sommets.

Mots clés : Tournées de véhicules, prize-collecting, problème du voyageur de commerce sélectif, multiobjectif.

Introduction

Much effort has been made to increase the oil production in Brazil though the use of new technologies. As a consequence, the Brazilian oil production has met the country's need in 2006 and the country is globally self sufficient. The Rio Grande do Norte basin has been exploited for the last 30 years and about 98% of the oil wells are pumped using artificial lift systems. One such system is the Mobile Oil Recovery (MOR) unit.

The MOR unit is an artificial lift system which is used to exploit wells whose production is marginal. It consists of a truck (vehicle unit) equipped with an apparatus, see Figure 1. Considering a working day, the unit starts its tour at the depot, then it pumps several wells before returning to the depot. At each well, the driver spends some time to connect the unit to the well, to pump the oil and to disconnect the equipment. Whenever the unit's tank is full, an auxiliary vehicle is used to transfer the oil from the MOR unit to its own tank and to carry it back to the depot. Thus, the MOR unit does not need to stop its operations and its capacity can be considered unlimited. Typically, the level of oil within a well is not static: it raises over the time until its stabilization level is reached [START_REF] Assmann | Simulador de pistoneio para indução de surgência[END_REF]. This happens when the pressure between the wells and the rock formation is stabilized. We consider here the oil level is static for all wells. It is computed by performing a statistic analysis of its history (last time it was pumped versus estimated production). The volume extracted from the wells is usually composed of a mixture of longchain hydrocarbons (petroleo) and basic sediments and water (bsw). Thus, the "oil" production of a well referred in this work takes into account the oil volume and its quality, which depends on its bsw factor.

The MOR optimization Problem (MORP) is a multiobjective problem which consists in finding a set of wells to be pumped in a working day to maximize the oil extraction and to minimize the travel time. The two objectives are opposite, one pushing to increase profit and the other to reduce costs. With one MOR unit, the problem is close to the Selective Traveling Salesman Problem [START_REF] Laporte | The selective travelling salesman problem[END_REF] which is also called Orienteering Problem, see e.g. [START_REF] Feillet | Traveling salesman problems with profits[END_REF], or Maximum Collection Problem [START_REF] Kataoka | An algorithm for the single constraint maximum collection problem[END_REF]. With a fleet of vehicles, the problem becomes a Vehicle Routing Problem (VRP) close to the Prize-Collecting VRP [START_REF] Balas | The prize collecting traveling salesman problem[END_REF]. For further investigation on routing problems, readers are referred to the following works: the state of the art on exact and approximated methods for the VRP and its variants are found in [START_REF]The Vehicle routing problem[END_REF] and an overview covering about 500 papers on classical routing problems are found in [START_REF] Laporte | Routing problems: A bibliography[END_REF]. For multiobjective solutions strategies on routing problems, see [START_REF] Boffey | Multiobjective routing problems[END_REF][START_REF] Current | Multiobjective design of transportation networks: taxonomy and annotation[END_REF][START_REF] Keller | The multiobjective vending problem: a generalization of the travelling salesman problem[END_REF].

A mathematical formulation for the MORP is proposed in [START_REF] Santos | Um algoritmo GRASP reativo aplicado ao problema do emprego da unidade móvel de pistoneio[END_REF] for a single MOR unit. A two-step optimization is performed: the oil extraction is maximized in the first phase. The second phase aims at minimizing both the travel and the operation time such that the oil extraction is not smaller than the value obtained in the first phase. Heuristics applications of the MORP are presented in [START_REF] Aloise | Uma colnia de formigas para o problema de explotao de petrleo e otimizao de rotas de unidades mveis de pistoneio[END_REF][START_REF] Santos | Um algoritmo GRASP reativo aplicado ao problema do emprego da unidade móvel de pistoneio[END_REF].

In this work, we propose several formulations for the MORP with a single vehicle and for a fleet of vehicles. Our contributions for the MORP with a single vehicle are: (i) remove redundant constraints, (ii) simplify the flow conservation constraints, (iii) test different strategies to eliminate invalid subtours, and (iv) strengthen the subtour elimination constraints. When a fleet of vehicles are used two formulations are proposed. The proposed formulations are strengthened by taking advantage of the MORP characteristics, by reducing the number of subtour elimination constraints, by introducing generalized bounds, and by adding cuts. Instances with up to 200 nodes, simulating practical situations, are solved. Moreover, it is shown that a formulation with a weaker linear relaxation can be very useful in practice, since it proves optimality for real size instances faster than the multiflow formulation. The first ever results for the MORP with several vehicles are presented in this work.

The work is organized as follows: the problem is defined in Section 2. Formulations for one unit are presented in Section 3. Sections 4.1 and 4.2 are devoted to the MORP with several units. Computational results are shown in Section 5 and final remarks are made in Section 6.

The MORP problem

The geographical data (roads, wells and depot) are modeled as an undirected graph G = (N, E). G is preprocessed to build a complete digraph G ′ = (V, A) where V is the set of wells and the depot v 0 . Let d ij be the shortest distance from i to j, ∀(i, j) ∈ A, and let s be the MOR unit average speed. Thus, for every arc of G ′ , the travel time t ij is computed as

t ij = d ij /s.
The process of dynamically filled wells was studied in [START_REF] Assmann | Simulador de pistoneio para indução de surgência[END_REF]. Let H i be the static level for well i and let us assume that well i is filled with a factor τ i . Wells are filled until the static level is reached, according to Equation [START_REF] Aloise | Uma colnia de formigas para o problema de explotao de petrleo e otimizao de rotas de unidades mveis de pistoneio[END_REF], where h i (t) is the oil level within well i at instant t. The bigger τ i , the slower well i is filled. We consider that wells can be exploited any time even if the static level is not reached and that the MOR units can visit a well only one time a day as in the previously works [START_REF] Barros | Uma aplicação de grasp na otimização do emprego da unidade móvel de pistoneio[END_REF][START_REF] Neves | Uma aplicação de algoritmo genético na otimização do emprego da unidade móvel de pistoneio[END_REF][START_REF] Santos | Sgump: Um sistema computacional georeferenciado de apoio à decisão para otimização do gerenciamento de unidades móveis de pistoneio[END_REF]. In fact, the proposed formulations do not deal with dynamically filled wells. A production estimation based on the well history is calculated each day for each well and this value is used instead. Thus, the oil level is considered static.

h i (t) = H i • (1 -e -t/τ i ) ∀i ∈ V (1)
Let t ′ i be the total operation time at well i (time to connect the MOR unit, to pump, and to disconnect the unit). Denote by p i = υ i • (1 -bsw i ) the oil production of well i (its prize), where bsw i is the amount of basic sediments and water, and when υ i is the estimated volume to be extracted. Let P and T be respectively the total oil production and the total time of a MOR unit in a working day. Moreover, T is the maximal time a MOR unit can work in a day.

Given K, the total number of MOR units, the MORP consists in defining one circuit τ = {v 0 , v σ1 , v σ2 , ..., v σk , v 0 } for each MOR unit, where σ is the position of wells in the circuit to be exploited in a day, such that P is maximized and T is minimized. The time limit T ≤ T has to be satisfied.

As mentioned before, the MORP with a single vehicle is close to two classical problems, the selective traveling salesman problem (STSP) and the prize collecting traveling salesman problem (PCTSP). In the STSP, a prize is associated to each node visited. The problem consists in defining a circuit (the tour starts and ends at a fixed node, e.g., the depot) such that the total prize is maximized (this corresponds to the oil production for the MORP), and such that the tour does not exceed a given length (this corresponds to the working time in the MORP) [START_REF] Laporte | The selective travelling salesman problem[END_REF]. The main difference between MORP and STSP is that the MORP also requires to minimize. Furthermore, in the MORP there are also an operation time associated to each node visited. In PCTSP, a prize is given for each visited node and a penality is paid for each unvisited node. A minimum level prize must be reached in a tour. The objective is to define a tour such that the travel cost and the penalities are minimized [START_REF] Balas | The prize collecting traveling salesman problem[END_REF]. The differences between the PCTSP and the MORP are: instead of paying a penality for each unvisited node, a penality (operation time) is paid for each visited node in the MORP. The objective function for the MORP has opposite goals: the primary objective pushes to collect prizes (the oil production) and the secondary aims at reducing costs. The MORP with several vehicles is close to a generalization of the PCTSP with a fleet of vehicles, denoted here of the prize-collecting VRP (PCVRP).

Formulations using a single vehicle

As far as we know, only one formulation has been proposed in the literature to the MORP [START_REF] Aloise | Otimizao do emprego da unidade mvel de pistoneio na explotao de campos petrolferos de poos no surgentes[END_REF][START_REF] Aloise | Reactive grasp applied to the oil recovery mobile system problem[END_REF]. In this formulation, one vehicle is considered and the optimization is done in two phases: first, compute the maximal amount of oil (prize) that can be extracted in a working day and second, compute the shortest route to extract this amount. We propose a formulation for the MORP with a single vehicle based on the formulation proposed in [START_REF] Santos | Um algoritmo GRASP reativo aplicado ao problema do emprego da unidade móvel de pistoneio[END_REF]. The main differences rely on: we have removed the constraints ensuring the MOR unit returns to the depot since (constraints 5 and 6 guarantee the return to the depot). The flow conservation constraints were simplified. We have also tested several subtour elimination strategies. In [START_REF] Santos | Um algoritmo GRASP reativo aplicado ao problema do emprego da unidade móvel de pistoneio[END_REF], only the Gavish and Graves strategy to eliminate subtour was used. Moreover, an upper bound on the number of wells that can be exploited in a working day is defined. This bound is used to strengthen the proposed formulation.

Let f ij ∈ {0, 1} be the decision variable on the choice of arc (i, j) and let x i be the binary variables which specify if well i is exploited or not. The first optimization phase for the MORP is given as follows:

max P = i∈V \{v 0 } p i • x i s.t. ( 2
)
i∈V \{v 0 } t ′ i • x i + (i,j)∈A t ij • f ij ≤ T (3) j:(j,i)∈A f ji - j:(i,j)∈A f ij = 0 ∀i ∈ V \{v 0 } (4) j:(j,i)∈A f ji = x i ∀i ∈ V \{v 0 } (5)
j∈V

f 0j = 1 (6)
(subtour eliminations constraints) (7)

x i ∈ {0, 1} ∀i ∈ V (8)
f ij ∈ {0, 1} ∀(i, j) ∈ A (9)
The objective function (2) aims at minimizing the oil extraction. Inequality (3) limits the working time (travel and operation time) to T . The flow conservation constraints are (4) and [START_REF] Balas | The prize collecting traveling salesman problem[END_REF]. Restriction (6) guarantees the tour starts at the depot. Variables x i and f ij are respectively defined in Constraints ( 8) and [START_REF] Current | Multiobjective design of transportation networks: taxonomy and annotation[END_REF]. We discuss in Section 3.1 the use of several subtour elimination constraints: those of Miller and Tucker and Zemlin (MTZ) [START_REF] Desrochers | Improvements and extensions to the MillerTuckerZemlin subtour elimination constraints[END_REF][START_REF] Miller | Integer programming formulations and traveling salesman problems[END_REF], and those of Gavish and Graves using either aggregated (GGA) or disaggregated flow (GGD) [START_REF] Gavish | The travelling salesman problem and related problems[END_REF]. GGA constraints are used in [START_REF] Santos | Um algoritmo GRASP reativo aplicado ao problema do emprego da unidade móvel de pistoneio[END_REF].

The objective of the second optimization phase is to minimize the working time [START_REF] Desrochers | Improvements and extensions to the MillerTuckerZemlin subtour elimination constraints[END_REF] subject to Constraints (4)-( 9) and [START_REF] Feillet | Traveling salesman problems with profits[END_REF]. Constraint [START_REF] Feillet | Traveling salesman problems with profits[END_REF] restricts the total production to be equal to the total optimal prize P * obtained in the first phase.

min T = i∈V \{v 0 } t ′ i • x i + (i,j)∈A t ij • f ij s.t. ( 10
)
i∈V \{v 0 } p i • x i = P * (11)
Constraints ( 4)-( 9).

Subtour eliminations constraints

A subtour is defined by any ordered subset of vertices. For the MORP, only subtours starting and ending at the depot v 0 are valid. Subtour constraints have been evaluated in the literature for the TSP problem, see e.g. [START_REF] Wong | Integer programming formulations of the traveling salesman problem[END_REF]. MTZ, GGA and GGD subtour elimination constraints for the MORP, and some improvements are described below. An upper bound on the number M of wells that can be exploited in a working day can be computed. Considering the working time of the MOR unit, a simple procedure consists in computing M by sorting the wells in increasing order of operation time t ′ i [START_REF] Santos | Um algoritmo GRASP reativo aplicado ao problema do emprego da unidade móvel de pistoneio[END_REF]. Thus, M is such that:

M i=1 t ′ i ≤ T ≤ M +1 i=1 t ′ i . (12)
We propose to strengthen the value of M by taking also into account the travel time using the following argument: to arrive at a node, the vehicle has to use an arc. Thus, the minimum travel time to arrive at each node can be considered as follows:

M i=1 t ′ i + {min j∈V t ji } ≤ T ≤ M +1 i=1 t ′ i + {min j∈V t ji } (13)
Another argument can be used to strengthen even more M. Since the vehicle must return to the depot, the minimum travel time to arrive at the depot can also be considered. Then, M is such that:

M i=1 t ′ i + {min j∈V t ji } ≤ T -{min j∈V t j0 } ≤ M +1 i=1 t ′ i + {min j∈V t ji } (14)
3.1.1 Lifted Miller, Tucker and Zemlin Constraints.

The Miller, Tucker and Zemlin constraints define a topological order to eliminate invalid subtours. However, for the MORP, the depot appears twice (at the beginning and at the end). Thus, one can duplicate the depot and work on a support graph. We consider instead the depot only at the beginning of the topological design. This can be done since the flow structure defined by variables f ij and x i , and Constraints ( 4)-( 6) guarantees the return to the depot.

The corresponding MTZ constraints for the MORP is given in Equations ( 15)- [START_REF] Laporte | Routing problems: A bibliography[END_REF]. Variable u i states the order well i appears in the solution, thus defining a topological order.

u i -u j + |V | • f ij ≤ |V | -1 ∀(i, j) ∈ A, j = {v 0 } (15) 0 ≤ u i ≤ |V | -1 ∀i ∈ V (16)
There are O(|V | 2 ) of such MTZ constraints, which can be strengthened by using M as computed in 14, instead of |V |.

u i -u j + M • f ij ≤ M -1 ∀(i, j) ∈ A, j = {v 0 } (17) 0 ≤ u i ≤ M ∀i ∈ V (18)
These constraints can be lifted even more using the same ideas as Desrochers and Laporte [START_REF] Desrochers | Improvements and extensions to the MillerTuckerZemlin subtour elimination constraints[END_REF]. The idea of lifting consists in adding a valid non-negative term α ji f ji to the Inequalities (17)

u i -u j + M • f ij + α ji • f ji ≤ M -1 ∀(i, j) ∈ A, j = {v 0 } (19)
The larger the value of α ji , the stronger the reduction in the original solution space. If f ji = 0, then α ji may take any value. Suppose now f ji = 1. Then, the MOR unit exploits well j = v 0 before well i, u i = u j + 1. Thus, f ji = 1 implies f ij = 0, otherwise there is a subtour (i, j). By substitution, we obtain α ji ≤ M -2. The larger α ji , the stronger is the lift. Thus, α ji = M -2. A lifted version of Constraints ( 17) is given in Inequalities [START_REF] Ropke | Models and branch-and-cut algorithms for pickup and delivery problems with time windows[END_REF].

u i -u j + M • f ij + (M -2) • f ji ≤ M -1 ∀(i, j) ∈ A, j = v 0 (20)
3.1.2 Gavish and Graves Constraints.

The Gavish and Graves [START_REF] Gavish | The travelling salesman problem and related problems[END_REF] approach removes invalid subtours by building a network flow. A flow is sent to the nodes of the tour. Each node consumes one unit. In disaggregated flow, a specific flow is sent from the source to each node [START_REF] Corte-Real | Network flow models for the local access network expansion problem[END_REF][START_REF] Magnanti | Optimal trees in network models[END_REF]. Otherwise, if the flow is not specified, it is an aggregated flow. Let y ij be the flow variable on arc (i, j). Thus, GGA constraints for the MORP are given in Equations ( 21)- [START_REF]The Vehicle routing problem[END_REF]. Constraints [START_REF] Santos | Sgump: Um sistema computacional georeferenciado de apoio à decisão para otimização do gerenciamento de unidades móveis de pistoneio[END_REF] are the flow conservation constraints. Inequalities [START_REF] Santos | Um algoritmo GRASP reativo aplicado ao problema do emprego da unidade móvel de pistoneio[END_REF] state a flow uses the arc (i, j) if it is selected. These constraints are strengthened by using M. In this strategy, there are O(|V | 2 ) constraints and variables. j:(j,i)∈A y jij:(i,j)∈A

y ij = x i ∀i ∈ V \{v 0 } (21)
y ij ≤ M • f ij ∀(i, j) ∈ A, (22)
y ij ≥ 0 ∀(i, j) ∈ A (23)
The GGD version is given in Constraints ( 24)-( 27). Let y l ij be the variable specifying if flow for node l traverses arc (i, j) or not. Constraints [START_REF] Wong | Integer programming formulations of the traveling salesman problem[END_REF] are the flow conservation constraints. Equations (25) state that a flow unit is sent from the source to each node l. Restrictions (26) specify that flow for node l traverses arc (i, j) if and only if it is used.

j:(i,j)∈A y l ij - j:(j,i)∈A y l ji = 0 ∀l ∈ V \{v 0 }, ∀i ∈ V \{v 0 , l} (24)
j:(0,j)∈A

y l 0j = x l ∀l ∈ V \{v 0 } (25) y l ij ≤ f ij ∀l ∈ V \{v 0 }, ∀(i, j) ∈ A (26) y l ij ≥ 0 ∀l ∈ V \{v 0 }, ∀(i, j) ∈ A (27)
This strategy implies O(|V | 3 ) of such constraints and variables and it produces better linear relaxation than using the aggregated flow ( 21)- [START_REF]The Vehicle routing problem[END_REF]. However, we show in the computational experiments that under some conditions, the aggregated flow strategy can produce a better linear relaxation as Constraints ( 22) are stronger with a lifted value for M.

Formulations using several vehicles

In this section, we present two formulations for the MORP using several MOR units. The first one is a three-indexed formulation. The second formulation is more compact. Both formulations are designed for the two optimization phases as presented in Section 3.

A three-indexed formulation using several vehicles

In the three-indexed formulation, we specify explicitly which MOR unit exploits each well. Thus, new variables are defined as follows: let x k i be a decision variable that specifies if well i is exploited by the vehicle k or not. Variables f k ij ∈ {0, 1} state if vehicle k exploits well j after well i or not. P (K) is the total profit collected using the K MOR units. All other terms are defined in Section 3. The three-indexed formulation is as follows:

max P (K) = i∈V \{v 0 } p i • K k=1 x k i s.t. ( 28
) i∈V \{v 0 } t ′ i • x k i + (i,j)∈A t ij • f k ij ≤ T ∀k = 1, ..., K (29)
j:(j,i)∈A f k ji - j:(i,j)∈A f k ij = 0 ∀k = 1, ..., K, ∀i ∈ V \{v 0 , k} (30)
j:(0,j)∈A

f k 0j ≤ 1 ∀k = 1, ..., K (31)
j:(j,i)∈A f k ji = x k i ∀k = 1, ..., K, ∀i ∈ V \{v 0 } (32) K k=1 x k i ≤ 1 ∀i ∈ V \{v 0 } (33) j:(j,i)∈A y ji - j:(i,j)∈A y ij = K k=1 x k i ∀i ∈ V \{v 0 } ( 34
)
y ij ≤ M • K k=1 f k ij ∀(i, j) ∈ A, j = v 0 ( 35
)
y ij ≥ 0 ∀(i, j) ∈ A (36) x k i ∈ {0, 1} ∀k = 1, ..., K, ∀i ∈ V \{v 0 } (37) f k ij ∈ {0, 1} ∀k = 1, ..., K, ∀(i, j) ∈ A (38)
Restrictions (29) limit the units work in a day. The flow conservation constraints are defined in (30) and (31). Constraints (32) ensure that unit k pass though an arc (i, j) only if it exploits well i. Inequalities (33) specify that at most one unit visits well i. Constraints (34) and ( 35) are the GGA subtour elimination constraints. Constraints (36)-( 38) are the variables definition. This formulation contains O(|V 3 |) variables and constraints and can be seen as a generalization of formulation ( 2)-( 9) to several vehicles. The GGA constraints are chosen according to the computational results for one vehicle (Section 5). Obviously, other strategies could be used as well.

Valid inequalities for the three-indexed formulation

In the previous model, any vehicle can be assigned to any route in the solution. Thus, the number of feasible solutions are multiplied by up to K (ways to assign K vehicles to the K routes). This can dramatically slow down the effectiveness of the model. Figure 2 illustrates the symmetry problem. Figure 2 (a) and(b) have the same tours, but in Figure 2 (a) vehicle k 1 does the tour {v 0 , v 3 , v 6 , v 7 , v 1 , v 0 } and in the Figure 2 (b), it is done by the vehicle k 2 . The same applies to the second tour. Thus, the idea is to use valid inequalities to avoid these situations.

To remove the symmetry, so-called symmetry-breaking constraints can be defined. A first way is to specify that the first vehicle does the tour with the highest oil production, the second one does the tour with the second highest oil production, and so on. This is done through Inequalities (39).

v 0 (a) (b) v 6 v 7 v 1 v 2 v 5 v 4 v 3 k 1 v 8 v 9 k 2 v 0 v 6 v 7 v 1 v 2 v 5 v 4 v 3 k 2 v 8 v 9 k 1
i =v 0 p i x k i ≥ i =v 0 p i x k+1 i ∀k = 1, ..., K -1 (39)
Another way to remove symmetry consists in applying a lexicographical order on the routes: given two circuits τ 1 and τ 2 , τ 1 < τ 2 if the first node visited in τ 1 has a lower identifier than the first node visited in τ 2 . Constraints (40) define such a lexicographical order.

i =v 0 i • x k 0i ≤ i =v 0 i • x k+1 0i ∀k = 1, ..., K -1 (40)
The lexicographical order can be defined another way as in Constraints (41). The sum of the first n ′ variables on the outgoing arcs is greater in the first circuits.

n ′ i=1 x k 0i ≥ n ′ i=1 x k+1 0i ∀k = 1, ..., K -1, ∀n ′ < n (41)
Finally, the last way to remove symmetry is to use Constraints (42). Vehicle k goes to node i only if vehicle k -1 goes to at least one of the nodes 1, ..., i -1). Consequently, all variables x k i such that k > i can be set to zero.

x k i ≤ i-1 j=1 x k-1 j ∀i = 1, ..., n, ∀k = 2, ..., min{K, i} (42)

A two-indexed formulation using several vehicles

In the two-indexed formulation, we do not explicitly define which unit exploits well i since the MOR units have the same characteristics (homogeneous fleet). A similar idea was previously used in the literature, for example, in [START_REF] Ropke | Models and branch-and-cut algorithms for pickup and delivery problems with time windows[END_REF]. Instead of explicitly assigning a MOR unit k to a subset of wells, a restriction is set on the date (time) well i is exploited using a unit, whatever it is. We consider the variables f ij and x i as defined in the Section 3 for the single vehicle formulation. Additionally, variables d i specify the date (time) well i is visited by a MOR unit in a working day.

max P = i∈V \{v 0 } p i • x i s.t. ( 43
) j:(j,i)∈A f ji - j:(i,j)∈A f ij = 0 ∀i ∈ V \{v 0 } (44) j:(j,i)∈A
f ji = x i ∀i ∈ V \{v 0 } (45)
j∈V

f 0j = K (46)
d i -d j +(T +t ′ i +t ij )•f ij +(T -t ′ j -t ji )•f ji ≤ T ∀(i, j) ∈ A, i, j = v 0 (47) d i ≥ t 0i • f 0i + j =v 0 (t 0j + t ′ j + t ji ) • f ji ∀i ∈ V \{v 0 } ( 48
)
d i ≤ T -(t ′ i + t i0 ) • f i0 - j =v 0 (t ′ i + t ij + t ′ j + t j0 ) • f ij ∀i ∈ V \{v 0 } (49) x i ∈ {0, 1} ∀i ∈ V (50)
f ij ∈ {0, 1} ∀(i, j) ∈ A ( 51
)
d i ≥ 0 ∀i ∈ V \{v 0 } (52)
The flow conservation is given in Constraints (44). Restrictions (45) ensure arc (j, i) is used if well i is exploited. Equations (46) state all the K MOR units are used. Constraints (47) link the time node j is visited, to the time node i is visited, and to the selection of arc (i, j). This is an adaptation of the lifted MTZ constraints (see Section 3.1). Inequalities (48) and (49) define generalized lower and upper bounds on the time node i is visited. Inequalities (48) link the time node i is visited to variables f ji . At most one of the arcs entering node i is used. Thus, d i is at least equal to the minimal time required to arrive at node i, either by going from v 0 to i or by going from j to i. The same idea applies to the Inequalities (49). Variables

Computational results

The computational experiments were carried out on an Intel Core 2 Duo with 2.66 GHz clock and 4Gb of RAM memory, using CPLEX 11 under default parameters. Instances were generated using a geographical information system to simulate real situations. Comparison among the proposed formulations are measured in terms of time to prove optimality and of linear relaxation.

In the Tables 1 and2, each line corresponds to an instance. For each instance, the working day length (L) in minutes, the number of wells (|V |) and its optimal production (P * ) are given. For each formulation, the linear relaxation value (RL * ), the time (T ) spent by the unit in the optimal solution, the time (time(s)) required to prove optimality in seconds (rounded up), and the number of nodes (nodes) explored in the branch-and-bound tree are presented. The symbol (-) means the solver did not prove optimality because it ran out of memory. When the optimal solution is unknown, the best integer solution found so far is identified by "(≥ value)".

Table 1 summarizes the results for the formulations for one vehicle using the MTZ, GGA or GGD subtour elimination constraints. From the computational results, GGA proves optimality faster than MTZ and GGD for 9 instances. MTZ proves optimality faster than GGA and GGD for 7 instances. In spite of having the worst linear relaxation, MTZ is able to prove optimality for instances with up to 200 nodes. GGD consumes a lot of time even if it produces good linear relaxation. An interesting result on the linear relaxation is found for L = 480 and |V | = 20: the GGA linear relaxation is better than the linear relaxation of GGD. This happens here because the value of M is equal to the optimal amount of wells exploited in a day.

We performed an experiment to illustrate the evolution of the GGA lin-ear relaxation values when M varies. We compare it with the GGD linear relaxation values. Results are presented for the instance with |V | = 20 wells in Figure 3. The x and y axis represent respectively the values of M, and the oil production. For this instance, the MOR unit can exploit at most M = 4 wells. The results show that the smaller M the tighter the linear relaxation of GGA. For some instances, in terms of linear relaxation, a value of M close to the optimal amount of wells exploited in a day is sufficient to make GGA better than GGD can even help GGA having better linear relaxation than GGD. We have also run the second optimization phase for all instances presented in Table 1. The time T was only improved for the instance with L = 480 and |V | = 120 (T * = 479.5 instead of T = 480). Thus, results of the second optimization phase were not tested for several vehicles. Even so, it remains valuable since it takes place in the global decision process of the problem.

The results for the two-indexed and the three-indexed formulations are presented in Table 2. The number of vehicles used (K) and the sum of the total time spent by all the MOR units (T ′ ) are given. The three-indexed formulation produces a better linear relaxation than the two-indexed formulation. However, the two-indexed formulation performs better to compute the optimal solution. In addition to the number of wells, the problem becomes more difficult when the number of vehicles increases. Moreover, the working day limit also contributes to the difficulty of the problem. Results suggest it is suitable to use a small time window (480 minutes). The three-indexed formulation found sometimes a smaller value of T ′ as shown in bold.

Concluding remarks

Several formulations for the MORP are proposed in this work and the first ever results using several vehicles are presented. Additionally, we proposed to improve the subtour constraints by taking advantage of the time window. Thus, instances close to reality (up to 200 wells) are solved. Among the formulations for one vehicle, GGA performs globally better than MTZ and GGD to prove optimality. For several vehicles, the two-indexed formulation is faster to prove optimality in spite of weaker linear relaxations. Computational experiments show that the time window restriction plays a key role in computing an optimal solution: the smaller the time window, the easier the problem to solve. Optimal solutions can be computed for mediumsized instances with two MOR units. When using three vehicles, this does not hold as the CPU time increases dramatically for small instances.

The larger instances used here are larger than the problems considered by the company in the Rio Grande do Norte Basin. Consequently, the oil company is now able to compute the optimal solution for the MORP instead of using solutions given by heuristics.

For future work, there are several promising research to be developed for this problem. For example, investigating situations where the second optimization phase becomes really useful. Moreover, the experiments show that for large time windows the problem becomes more difficult. Thus, we could study an approach to split large time windows. In terms of algorithms, several strategies could be explored such as exact algorithms as a branch-andcut. It could be applied for difficult cases when a large number of vehicles are needed. Furthermore, pareto-based strategies to deal with the multiobjective function can be tested. Finally, the dynamic oil filling in the well can be integrated the models.

Figure 1 :

 1 Figure 1: An example of a MOR vehicle.

Figure 2 :

 2 Figure 2: An example of symmetry between two solutions.

  x k i , f k ij and d i are defined respectively by Constraints (50) to (52). The twoindexed formulation has O(|V | 2 ) variables and constraints. MTZ is used since the time constraints definition is straightforward and the formulation has still O(|V | 2 ) variables.

Figure 3 :

 3 Figure 3: Evolution of the linear relaxation for the formulations GGA and GGD.

Table 1 :

 1 The first optimization phase for the MORP using one vehicle.

	MTZ	GGA	GGD

Acknowledgement / Remerciements

We thank Petrobras staff for providing valuable informations about the MOR unit and its usage in a petrol field. We thanks also Haugland Dag for the suggestion of using the symmetry-breaking constraints to improve the three-indexed formulation.