The multiple vehicle pickup and delivery problem with LIFO constraints

Discrete Optimization

Accepted Manuscript

The multiple vehicle pickup and delivery problem with LIFO constraints E. Benavent, M. Landete, E. Mota, G. Tirado PII: DOI: Reference:

S0377-2217(14)01047-9 10.1016/j.ejor.2014.12.029 EOR 12688

To appear in:

European Journal of Operational Research

Received date: Accepted date:

11 April 2014 15 December 2014

Please cite this article as: E. Benavent, M. Landete, E. Mota, G. Tirado, The multiple vehicle pickup and delivery problem with LIFO constraints, European Journal of Operational Research (2015), doi: 10.1016/j.ejor.2014.12.029

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • Approach a new pickup and delivery problem with LIFO and time constraints

CR IP T

• Propose two formulations for the problem, together with several set of valid inequalities to improve them • Design an efficient branch & cut solution method to solve the problem to optimality • Develop a fast tabu search heuristic to find feasible solutions within short running times

AC

CE

PT

ED

M

AN US

• Evaluate the performance of the proposed solution methods through a comprehensive computational study

1

ACCEPTED MANUSCRIPT

The multiple vehicle pickup and delivery problem with LIFO constraints

CR IP T

E. Benaventa , M. Landeteb , E. Motaa , G. Tiradoc,∗ a Departamento

de Estadística e Investigación Operativa, Universidad de Valencia de Estadística, Matemáticas e Informática, Universidad Miguel Hernández c Departamento de Estadística e Investigación Operativa, Universidad Complutense de Madrid b Departamento

Abstract

M

AN US

This paper approaches a pickup and delivery problem with multiple vehicles in which LIFO conditions are imposed when performing loading and unloading operations and the route durations cannot exceed a given limit. We propose two mixed integer formulations of this problem and a heuristic procedure that uses tabu search in a multi-start framework. The first formulation is a compact one, that is, the number of variables and constraints is polynomial in the number of requests, while the second one contains an exponential number of constraints and is used as the basis of a branchand-cut algorithm. The performances of the proposed solution methods are evaluated through an extensive computational study using instances of different types that were created by adapting existing benchmark instances. The proposed exact methods are able to optimally solve instances with up to 60 nodes.

PT

1. Introduction

ED

Keywords: vehicle routing, pickup and delivery, LIFO constraints, mixed-integer programming, branch and cut

AC

CE

Vehicle routing problems play a very important role in logistics and transportation and have been widely studied during the last decades. Many variants with different characteristics have been considered and a variety of mathematical and computational techniques have been proposed to solve them (see, for example, [26] and [13] for a survey on vehicle routing problems). The Capacitated Vehicle Routing Problem (CVRP) is the reference model for these problems and is defined as follows: given a set of customers, each one demanding a certain amount of a product, and a fleet of vehicles with limited capacity that are based at a given depot, the problem consists of finding a set of minimum cost routes for the vehicles so that all customers are ∗ Corresponding author: G. Tirado. Postal address: Departamento de Estadística e Investigación Operativa, Facultad de Matematicas, UCM, Plaza de Ciencias 3, 28040, Madrid, Spain. Tel.: +34 913944377, Fax: +34 913944606. Email addresses: [email protected] (E. Benavent), [email protected] (M. Landete), [email protected] (E. Mota), [email protected] (G. Tirado)

Preprint submitted to European Journal of Operational Research

January 7, 2015

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

visited exactly once, their demands are fulfilled and the capacity of the vehicles is not exceeded. The problem approached in this paper is a variant of the so called Vehicle Routing Problem with Pickups and Deliveries (VRPPD) in which the vehicles have to satisfy a set of customer requests where each request specifies the size of the load to be transported, the origin location and the destination location. According to the classification proposed by Berbeglia et al. [3] this problem belongs to the one-to-one class of pickup and delivery problems. Note that the CVRP can be considered as a particular case of the VRPPD where all the requests have the depot as the origin location. The VRPPD has been extensively studied in the last three decades (see [3] and [21, 22]), incorporating additional constraints that appear in practical situations. In this paper we consider a LIFO (Last-In-First-Out) rule of service, that means that only the last pickup customer request that has been loaded can be delivered. Hence, the vehicle servicing an origin location must continue to the corresponding destination or visit another origin location. This condition arises naturally when the vehicles used for transportation are rear-loaded and they have a single access point to their container. The rearrangement of the load when performing deliveries may be possible, but if it is too time consuming, it could be simply forbidden. This is the case, for example, of hazardous materials or very heavy or fragile items. Moreover, we also require that the total duration of each route does not exceed a given limit, where route duration includes traveling times among different locations and service times at each origin/destination for loading and unloading operations. The objective is to find a set of routes with minimum total duration satisfying all the constraints. We refer to this problem as the Multiple Vehicle Pickup and Delivery Problem with LIFO constraints and Maximum Time but it will be abbreviated simply by PDPLT. To our knowledge, the PDPLT has not been previously studied. However, algorithms have been proposed for related problems. If the capacity and maximum time constraints are relaxed the resulting problem is known as the Traveling Salesman Problem with Pickup, Delivery and LIFO constraints (TSPPDL). Carrabs et al. [4] develop a branch-and-bound algorithm for this problem and later Cordeau et al. [8] propose a branch-and-cut algorithm which is able to solve instances with up to 25 requests within a reasonable computing time. Battarra et al. [1] consider a version of the TSPPDL in which the LIFO rule can be violated, but in this case rehandling costs are incurred. They present models and exact algorithms based on branch-and-cut for this problem, while Erdogan et al. [10] have later developed a tabu search metaheuristic. Cheang et al. [5] propose several heuristics for a multi-vehicle case where the route length of each vehicle cannot exceed a predetermined limit and the vehicles have unlimited capacity (see also [11]). This problem is called Multiple Pickup and Delivery Traveling Salesman Problem with LIFO and Distance constraints (MTSPPD-LD). Their primary objective is minimizing the number of vehicles used, while the secondary objective is minimizing the total distance traveled. Note that the MTSPPDLD can be solved as a PDPLT, since minimizing the total distance is equivalent to minimizing the total time, and minimizing the number of vehicles as the primary objective can be achieved by simply adding a big number to the travel times of the 3

ACCEPTED MANUSCRIPT

ED

M

AN US

CR IP T

arcs leaving the depot. Another variant that is closely related to the PDPLT was introduced recently by Cherkesly et al. [6]. In that paper, in addition to imposing LIFO constraints, time windows are associated to every pickup and delivery location. However, they do not consider a maximum time duration for the routes. The authors develop three related branch-price-and-cut algorithms using a three-index formulation of the problem. Their exact solution procedures involve column generation techniques and constrained shortest path subproblems. The main purpose of this paper is to propose, for the first time, methods to solve the PDPLT. We have developed two mixed integer formulations of the PDPLT and a heuristic procedure that uses tabu search in a multi-start framework. The first formulation is a compact one, that is, the number of variables and constraints grows polynomially with the number of requests, while the second one contains an exponential number of constraints. Separation procedures for several families of the constraints of the non-compact formulation are implemented and embedded in a branch-and-cut algorithm. The compact formulation is easy to use because all the constraints can be included at once, but it is outperformed by the proposed branch-and-cut based on the non-compact formulation for larger instances. The proposed methods are able to solve medium size instances of the PDPLT and provide tight lower and upper bounds for larger instances. The remainder of the paper is organized as follows. Section 2 introduces the main notation that will be used through the paper. Section 3 presents the compact formulation and Section 4 is devoted to the non-compact formulation, separation algorithms for several families of valid inequalities and the associated branch-and-cut algorithm. Section 5 describes a fast tabu search heuristic for the PDPLT which is embedded into a multi-start framework. The performance of the proposed solution methods is evaluated and analyzed in Section 6, where an extensive computational study is presented. Finally, Section 7 concludes with some final remarks and ideas for future research.

PT

2. Problem notation

The PDPLT is defined on a directed graph G = (V, A).

CE

• V = P ∪ D ∪ {0} is the set of nodes, where P = {1, . . . , n} is the set of pickup locations, D = {n + 1, . . . , 2n} is the set of delivery locations and 0 is the depot. We assume that the delivery location of the load picked up at i ∈ P is i+n ∈ D.

AC

• A is the set of arcs (i, j) between every pair of nodes i, j ∈ V , i 6= j.

• The travel time to traverse each arc (i, j) is denoted by cij , for all (i, j) ∈ A. It is assumed that travel times satisfy the triangular inequality. Travel times will also be called costs throughout the paper. • The service time (pickup or delivery) at each node i ∈ P ∪ D is si . • The maximum duration of the route of any vehicle is T. 4

ACCEPTED MANUSCRIPT

• For each i ∈ P , di > 0 is the load of the item that must be picked up at i and delivered at i + n. We define, for each i + n ∈ D, di+n = −di < 0. • The number of available vehicles is denoted by m and their capacity by Q.

CR IP T

The goal of the PDPLT is to find a set of routes with minimum total duration that services all the requests, that is, each requested item must be picked up at its origin and delivered at its destination. The duration of the routes include the traversing times among locations and service times and cannot exceed T . The routes must start and end at the depot, the capacity cannot be exceeded, and the LIFO rule must be respected. 3. A compact formulation

CE

PT

ED

M

AN US

In this section we propose a compact formulation for our problem using three families of variables. Binary variables xij indicate if arc (i, j) is traversed by any vehicle, or not. Continuous variables αij are typical flow-like variables that in our case indicate the arrival time at node j when the vehicle travels from i to j. Binary variables uik , defined for any pair of nodes i ∈ V and k ∈ P , are equal to 1 if the item corresponding to node k is in the vehicle when it leaves node i. These variables are used to prevent subtours, to express the capacity constraints and to assure that the pickup and the delivery nodes associated to a given request belong to the same route. This last condition is not easy to impose because variables xij do not allow to identify which vehicle is traveling from i to j. Variables uik and their corresponding constraints are inspired by the Miller Tucker Zemlin formulation for the Travelling Salesman Problem (see [20]) that were improved later and used for the Vehicle Routing Problem by Kara et al. [15]. Let us define coefficients βik as follows: βkk = 1 and βk+n,k = −1 for all k ∈ P , and βik = 0 for any i ∈ V, i 6= k, n + k. Note that coefficients βik indicate that item k can only be picked up in node k (βkk = 1) and delivered in node k + n (βk+n,k = −1). For instance, if n = 3, it is   1 0 0  0 1 0     0 0 1   β=  −1 0 0     0 −1 0  0 0 −1

AC

Proposition 1. Constraints (1) are valid for the PDPLT. uik − ujk + xij + (1 − βik − βjk )xji ≤ 1 − βjk

∀(i, j) ∈ A, k ∈ P

(1)

Proof : Let us first prove the validity of constraints (2). uik − ujk + xij ≤ 1 − βjk

5

∀(i, j) ∈ A, k ∈ P

(2)

ACCEPTED MANUSCRIPT

j6=i

AN US

Proposition 2. Constraints (3) are valid for the PDPLT. X uik − βjk xji ≥ βik ∀i ∈ V, k ∈ P

CR IP T

Note that the definition of variables uik implies that ukk = 1 and un+k,k = 0 for all k ∈ P . If xij = 0, then the possible values of uik − ujk are −1, 0, 1 and, since the RHS is always greater than or equal to zero, we only have to analyze the case where ujk = 0 and uik = 1, but then βjk 6= 1 and the constraint holds. If xij = 1, then uik − ujk ≤ −βjk . If uik = 1, then either ujk = 1 and βjk = 0 or ujk = 0 and βjk = −1 (j = k + n). The case where xij = 1 and uik = 0 is proved similarly. These constraints can be improved by introducing an extra term, in a similar way as in [15] : uik − ujk + xij + γji xji ≤ 1 − βjk . It can be observed that the largest value that can be given to γji is (1 − βik − βjk ) ≥ 0.  Following the same lines as in [15], Proposition 2 presents another set of lifted constraints.

(3)

Proof : Note P that in a feasible solution there is only one arc entering each node and therefore j6=i βjk xji is equal to βi0 k where i0 is the node from which the vehicle enters node i. If uik = 1 then: either i = k and βi0 k cannot be 1 or i 6= k and therefore item k is already on the vehicle when it leaves node i0 (in the case that i0 = k the constraint still holds). The case when uik = 0 is proved similarly. 

M

The following notation will be used. For any subsets S, S 0 ⊆ V , we denote (S : S ) = {(i, j) ∈ A : i ∈ S, j ∈ S 0 }, δ + (S) = (S : V \S), δ − (S) = (V \S : S), δ(S) = δ + (S) ∪ δ − (S), and γ(S) = (S : S). For simplicity, P we will write i instead of {i} when using this notation. Finally, we denote x(F ) = (i,j)∈F xij , for any F ⊆ A. Similar notation will be used with α-variables.

AC

CE

PT

ED

0

6

ACCEPTED MANUSCRIPT

Then, the basic formulation of our problem is: P min (i,j)∈A (cij + si )xij

uik − ujk + xij + (1 − βik − βjk )xji ≤ 1 − βjk ∀(i, j) ∈ A, k ∈ P P uik − j6=i βjk xji ≥ βik ∀i ∈ V, k ∈ P ukk = 1,

u0k = 0

un+k,k = 0 uik = ui+n,k P k∈P dk uik ≤ Q

(5) (6)

CR IP T

s.t.

(4)

∀k ∈ P ∀k ∈ P

∀i ∈ P, k ∈ P, i 6= k

(7) (8)

(9)

∀i ∈ V

(10)

∀i ∈ P ∪ D

(12)

∀j ∈ P

(11)

∀(i, j) ∈ A

(13)

∀i ∈ P ∪ D

(15)

αij ≥ 0

∀(i, j) ∈ A

(18)

xij ∈ {0, 1}

∀(i, j) ∈ A

(20)

α0j = c0j x0j P

α(δ + (i)) − α(δ − (i)) =

j∈V :j6=i(cij +si )xij

−

AN US

αij ≤ T xij

x(δ (j)) = 1 x(δ + (i)) = 1 −

x(δ (0)) ≤ m x(δ + (0)) ≤ m

∀i ∈ V, k ∈ P

M

uik ∈ {0, 1}

∀j ∈ P ∪ D

(14)

(16) (17) (19)

AC

CE

PT

ED

It can be observed that constraints (5)-(8) avoid subtours disconnected from the depot and force any pair of nodes i ∈ P, i + n ∈ D to belong to the same route. Constraints (9) guarantee that the items loaded when leaving node i are the same that when leaving node i + n for all k ∈ P, i 6= k because, as imposed by the LIFO requirement, any item loaded when traveling form i to n + i must be delivered before the item corresponding to i is delivered in node n + i. Constraints (10) assure that the total load on the vehicle when leaving a node cannot exceed the vehicle capacity Q. The values of the α variables (arrival time at each node) when departing from the depot are stated by constraints (11). Flow conservation constraints (12) define the values of the other α variables and constraints (13) introduce the route duration limit. Constraints (14) and (15) impose that any node different from the depot is visited exactly once, while constraints (16) and (17) set an upper bound on the number of vehicles used. Finally, (18)-(20) are the non-negativity and integrality constraints respectively. It can be observed that (4)-(20) define a valid formulation of the PDPLT. 3.1. Improving the compact formulation We now consider other sets of valid inequalities, that in some cases are improved versions of the above constraints. Capacity constraints are valid for any node i, but

7

ACCEPTED MANUSCRIPT

if we consider a pick-up node i ∈ P we can write: X dk uik + (Q − di )x0i ≤ Q ∀i ∈ P

(21)

k∈P

AN US

CR IP T

P because if x0i = 1 then k∈P dk uik = di and the inequality holds. In a similar way we can improve the trivial constraint uik ≤ 1 ∀i ∈ V, k ∈ P obtaining: uik + (1 − βik )x0i ≤ 1 ∀i, k ∈ P (22) P Given a node i, k∈P dk uik corresponds to the load in the vehicle leaving it and P d x is a lower bound of the load of the vehicle entering it. Therefore, the j6=i j ji difference between these two quantities is an upper bound on the load picked up at node i, if i ∈ P , or minus the load delivered at i, if i ∈ D, so the following inequality is valid: X X dk uik − dj xji ≥ di ∀i ∈ V (23) k∈P

j6=i

M

The load of a vehicle when leaving node i cannot be greater than Q; on the other hand the load when leaving a node j is alwaysPat least equal to dj and therefore, the difference cannot be greater than Q − dj , i.e.: k∈P dk (uik − ujk ) ≤ Q − dj ∀(i, j) ∈ A. These constraints can be improved and it can be observed the following constraints are also valid: X dk (uik − ujk ) + Qxij + (Q − di − dj )xji ≤ Q − dj ∀(i, j) ∈ A (24) k∈P

PT

ED

Both constraints (23) and (24) are obtained following the same lines as in [16]. The following constraints ensure that node n + i is visited at least ci,n+i + si units of time after visiting node i, which also implies that n + i is always visited after i. Note that these constraints are valid since we assume that travel times satisfy the triangle inequality. X X αj,n+i − αji ≥ ci,n+i + si ∀i ∈ P (25) j∈V :j6=n+i

j∈V :j6=i

AC

CE

On the other hand, the following x-variables can be fixed.

xi0 = 0 ∀i ∈ P

(26)

xij = 0 ∀i ∈ P, j ∈ D : j 6= i + n

(28)

xij = 0 ∀i, j ∈ D : di + dj < −Q

(30)

x0j = 0 ∀j ∈ D

(27)

xij = 0 ∀i, j ∈ P : di + dj > Q

(29)

xij = 0 ∀i, j ∈ V \ {0} : c0i + cij + cj0 > T

(31)

Equations (26) and (27) forbid, respectively, to travel from a pickup node to the

8

ACCEPTED MANUSCRIPT

αij ≥ cij xij

∀(i, j) ∈ A

CR IP T

depot and from the depot to a delivery node. Equations (28) are valid as a consequence of the LIFO strategy: if an item is picked up in a node, the next node visited cannot be the delivery node of a different item. Equations (29)-(30) state that two pickup or delivery nodes with a combined demand exceeding the capacity of the vehicle cannot be visited consecutively. Equations (31) indicate that if the total time required to go from the depot to a node i, from i to j and from j back to the depot exceeds the maximum driving time T , nodes i and j cannot be visited consecutively. Furthermore, (32)

AN US

are also valid. Finally, given that the objective function expresses the total time needed to visit all the pickup and delivery nodes, if LB is a valid lower bound on the optimal cost, dLB/T e is a lower bound of the number of vehicles that have to be used, so the following constraint is valid: x(δ + (0)) ≥ dLB/T e

(33)

M

The compact formulation that was finally implemented includes constraints (5)(20), where constraints (10) were replaced by (21) for the pick-up nodes, and constraints (22)-(28). Furthermore, relationships (7) to (9) were used to reduce the number of variables ujk (see the Appendix for further details). 4. A non-compact formulation

AC

CE

PT

ED

In this section we present a formulation that uses only the integer variables xij to define the routes and the continuous variables αij to control the duration of the routes. These variables have been already defined in Section 3. The main drawback of this new formulation is that it contains an exponential number of constraints, so it has to be used in a cutting plane scheme that will be embedded in a branch-and-cut procedure in order to find the optimal solution. This formulation includes the same objective function (4) as the compact formulation, constraints (14)-(17) that control the node degrees with respect to x-variables, additional constraints (26)-(33), and constraints (11)-(13) that concern α-variables. Finally, constraints (18) and (20) are also included in this formulation. In what follows we will introduce several families of constraints that are necessary to complete the formulation. Additional constraints are also included to strengthen the formulation so that the linear relaxation provides a better lower bound on the optimal solution cost. Given a vector of binary x-variables x, that represents a feasible solution of the PDPLT, we denote by Gx the graph whose node set is V and contains the arcs (i, j) such that xij = 1, for all (i, j) ∈ A.

9

ACCEPTED MANUSCRIPT

4.1. Capacity constraints The vehicle capacity constraints are expressed in a similar way as it is done for other capacitated routing problems: (34)

CR IP T

P

x(δ + (S)) ≥ max{1, dD(S)/Qe} ∀S ∈ V \{0}

where D(S) = | i∈S di | is the net demand of the nodes in S. These constraints can also be written as: x(δ(S)) ≥ 2 max{1, dD(S)/Qe} or, equivalently x(γ(S)) ≤ |S| − max{1, dD(S)/Qe}. Note that if S = {i, j} and D({i, j}) ≤ Q the capacity constraint can be written as xij + xji ≤ 1. 4.2. Precedence constraints

AN US

The condition that the pickup is an activity preceding the delivery can be imposed as follows: x(γ(S)) + x(S : i) ≤ |S| − 1 ∀S ⊂ V \{0}, i ∈ / S, n + i ∈ S, i ∈ P

(35)

Note that a solution containing a path that visits node n + i before node i would violate this constraint for a set S containing all the nodes in the path between n + i and the node preceding i. The above precedence constraints can be strengthened because a path visiting every node in S cannot visit afterwards any pickup node whose corresponding delivery node is in S.

M

Proposition 3. Constraints (36) are valid for the PDPLT.

x(γ(S)) + x(S : S 0 ) ≤ |S| − 1 ∀S ⊂ V \{0} S 0 = {i ∈ P \S : n + i ∈ S} = 6 ∅ (36)

CE

PT

ED

Proof : LetP us consider a feasible solution x. Since x(δ + (i)) = 1 for all i ∈ V \{0}, we have that i∈S x(δ + (i)) = x(γ(S)) + x(S : V \S) = |S|. Note that the subgraph of Gx induced by S consists of a number of node-disjoint paths. Let t ≥ 1 be the number of paths, so that x(γ(S)) = |S| − t. Therefore x(S : S 0 ) ≤ x(S : V \S) = t and all we need to show is that x(S : S 0 ) < t. Note that for each arc (i, j) ∈ (S : S 0 ) used by the solution, there should exist a path following that arc that returns to S without visiting the depot because the delivery node corresponding to node j is in S. Therefore, if x(S : S 0 ) = t, the solution would contain at least one cycle not visiting the depot, which is not allowed. Thus, the constraint is satisfied by x. 

4.3. Same-route constraints

AC

These constraints impose that the same vehicle services nodes i and n + i:

x(0 : S) + x(γ(S)) ≤ |S| − 1 ∀S ⊂ V \{0} such that ∃n + i ∈ S ∩ D with i ∈ / S (37)

If a route contains node n + i but not node i, then define S as the set of nodes visited by the route from the depot, without including it, to node n + i, and the above constraint is violated.

10

ACCEPTED MANUSCRIPT

Proposition 4. Constraints (37) are valid for the PDPLT.

CR IP T

Proof : Let i∗ ∈ P \ S such that n + i∗ ∈ S, and let us consider a feasible solution x. As in the proof of (36), note that x(γ(S)) = |S| − t, t ≥ 1 implies that there are t node-disjoint paths in S; given that there is no arc from the depot to the path in S containing node n + i∗ , x(0 : S) ≤ t − 1 holds and the proof is complete. 

There is a symmetric family to these inequalities that avoid routes visiting a pickup node and not visiting the corresponding delivery node. These constraints are not necessary for the integer formulation if (37) are present, but they may be useful to strengthen the linear relaxation and to cut fractional solutions. x(γ(S)) + x(S : 0) ≤ |S| − 1 ∀S ⊂ V \{0} such that ∃i ∈ S ∩ P with n + i ∈ / S (38)

AN US

If there is a route servicing node i ∈ S ∩ P but not node n + i, then considering S as the set of nodes in the route from i to the depot, without including it, the above constraint is violated. The proof of validity is similar to the one above.

M

4.4. LIFO constraints Cordeau et al. [8] introduced the following constraints that assure the LIFO policy in the TSPPDL: x(i : S) + x(γ(S)) + x(S : n + i) ≤ |S|

∀S ∈ Ω, ∀i, n + i ∈ / S, i ∈ P

(39)

ED

where Ω is the family of all subsets S ⊂ P ∪ D such that there is j ∈ S ∩ P and n+j ∈ / S or there is n + j ∈ S ∩ D and j ∈ / S. Propositions 5 and 6 introduce new strengthened LIFO constraints that dominate and extend constraints (39).

PT

Proposition 5. The following constraints are valid for the PDPLT, x(γ(S)) + x(S : S 0 ) ≤ |S| − 1

(40)

CE

where subset S ⊂ V \{0} is such that there exists k ∈ S ∩ P for which n + k ∈ / S, and S 0 = {n + j : n + j ∈ D\S and j ∈ / S} ∪ {j : j ∈ P \S and n + j ∈ S}. Note that n+k ∈ / S0.

AC

Proof: Let us consider a feasible solution x. As in the proof of Proposition 3, we know that x(γ(S)) + x(S : V \S) = |S|. Then, all we need to show is that there is at least one arc in Gx from a node in S to a node not in S 0 . Note that the subgraph of Gx induced by S consists of a number of node-disjoint paths, let P1 be the path of the subgraph of Gx induced by S that contains node k. Let (l, y), l ∈ S, y ∈ V \S, be the arc leaving S that continues path P1 . If y = n + k, we are done, since n + k ∈ / S0. Otherwise, given that the corresponding vehicle has not yet visited node n + k, y cannot be a delivery node in S 0 . If y = j ∈ S 0 ∩ P , since n + j ∈ S, the vehicle must 11

ACCEPTED MANUSCRIPT

CR IP T

visit S again (without visiting n + k), in order to perform the delivery at node n + j, and leave again set S through a new arc, say (l0 , y 0 ). If y 0 = j 0 ∈ S 0 ∩ P , the same argument applies, and so on. Given that the number of arcs is finite, some arc leaving S will enter a node not in S 0 . This is illustrated in Figure 1. 

Figure 1: Illustration of the proof of Proposition 5

AN US

Note that constraints (40) dominate constraints (39) for the case where set S ∈ Ω is such that there exists j ∈ S ∩ P such that n + j ∈ / S, because x(i : S) ≤ 1, and n + i ∈ S 0 if i, n + i ∈ / S. Proposition 6. The following constraints are valid for the PDPLT, x(i : S) + x(γ(S)) + x(S : S 0 ) ≤ |S|

(41)

M

where sets S and S 0 and node i ∈ P satisfy one of the following sets of conditions: C1. S ⊂ V \{0} is such that there exists k ∈ P \S for which n + k ∈ S ∩ D, and i ∈ P \S. S 0 = {i, k} if n + i ∈ S, and S 0 = {n + i} if n + i ∈ / S.

PT

Proof:

ED

C2. S ∈ / Ω (this means that for all j ∈ P either j, n + j ∈ S or j, n + j ∈ / S), i ∈ P \S, n + i ∈ D\S, and S 0 = {i} ∪ {n + j : n + j ∈ / S, j ∈ / S, j 6= i}

AC

CE

C1. Let us consider a feasible solution x. Again, we know that x(γ(S)) + x(S : V \S) = |S|. Then, if x(i : S) = 0 the constraint is obviously satisfied. If x(i : S) = 1, all we need to show is that there is at least one arc in Gx leaving set S and not entering S 0 . Note that the subgraph of Gx induced by S consists of a number of node-disjoint paths, let P1 be the path of the subgraph of Gx induced by S that follows the arc entering S from node i. We consider two cases: (a) n + i ∈ / S. Then n + k cannot be in P1 because before delivering the item at n + k, the item corresponding to i must be delivered at n + i. Therefore, x(γ(S)) ≤ |S| − 2, so x(S : V \S) ≥ 2, and then given that S 0 = {n + i} in this case, the constraint is satisfied. (b) n + i ∈ S. Let (l, y) be the arc leaving S that continues path P1 . If n + k is in the path P1 , then obviously y ∈ / S 0 = {i, k}. On the other hand, if n + k is not in the path P1 , it is clear that y 6= i, and if y 6= k, we are done. 12

ACCEPTED MANUSCRIPT

CR IP T

If y = k, the vehicle that follows path P1 must enter S again in order to deliver the item at node n + k, therefore it will have to leave S again, this time entering to a node different from nodes i and k, so the constraint is satisfied.  C2. Let us consider a feasible solution x. If x(i : S) = 0 the constraint is obviously satisfied. Assume x(i : S) = 1, and let P1 be the path of the subgraph of Gx induced by S that follows the arc entering S from node i. It is then clear that the arc that continues path P1 to leave S cannot enter any node in S 0 , so the constraint is also satisfied in this case. 

AN US

Note that constraints (39), for the case where set S ∈ Ω is such that there exists n + j ∈ S ∩ D so that j ∈ / S, are exactly the same as constraints (41) with conditions C1, if n + i ∈ / S. On the other hand, conditions C2 introduce new cases that were not considered by Cordeau et al. [8]. Any integer solution that satisfies constraints (11)-(18), (20), capacity constraints (34), precedence constraints (36), same-route constraints (37, 38), and LIFO constraints (40,41) defines a set of routes that are a feasible PDPL solution. Additional constraints (26)-(33) are also used in this formulation.

CE

PT

ED

M

4.5. Branch-and-cut The non-compact formulation is the base of a branch-and-cut algorithm that is described next. The initial linear programming (LP) relaxation at the root node of the branch-and-cut tree contains the objective function (4), constraints (11)-(18), (26)-(32), and xij ≥ 0, ∀(i, j) ∈ A. At any node of the branch-and-cut tree a cutting plane algorithm is used that works as follows. At each iteration, the current LP is solved and several separation procedures are used to search for inequalities that are violated by the optimal LP solution. These inequalities belong to the families described in the previous section that contain an exponential number of constraints, namely (34), (36)-(38), and (40,41). After the addition of violated constraints, the LP is reoptimized using the dual simplex algorithm. The cycle of LP optimization and addition of violated constraints iterates until no more violated inequalities are found, in which case, the branch-and-cut node is branched. The strong branching strategy of CPLEX [14] is used to explore the branch-and-cut tree. Inequality (33) is added to the LP after the first LP is solved and a valid lower bound on the total duration of the routes is known. The right-hand-side of (33) is checked at each iteration and updated whenever it is possible.

AC

Separation algorithms In what follows we describe the separation procedures that have been used for each family of constraints. Let (x, α) be the solution of the linear relaxation of the current LP. We define the directed weighted graph G[x] as the graph induced by arcs (i, j) with xij > 0, which will be called the support graph. In G[x] the weight of an arc (i, j) ∈ A is given by xij . Additionally, we consider the undirected graph Gu [x] where each pair of arcs (i, j) and (j, i) are replaced by an edge e = {i, j} with weight xe = xij + xji . 13

ACCEPTED MANUSCRIPT

CE

PT

ED

M

AN US

CR IP T

Most of the separation procedures work by generating "promising" sets of nodes that are used to check the violation of several types of constraints. Given a subset of nodes S ⊂ V \{0}, the capacity constraints (34) for subsets S and V \{0}\S are checked for possible violation. In the case that D(S) ≤ Q the precedence constraint is checked instead, because it is stronger than the capacity constraint. Concerning the LIFO constraints (40,41), each time a subset S ⊆ V \{0} is generated, we determine the type of LIFO constraint associated to it, and then, for each i ∈ P \S such that x(i : S) > 0, the appropriate subset S 0 is generated and the corresponding LIFO constraint is checked. Precedence constraints (36) for a given S ⊂ V \{0} are checked if the corresponding subset S 0 is not empty. Finally, same route constraints (37) and (38) are easily checked for a given S ⊂ V \{0} by first computing x(0 : S) or x(S : 0), respectively. Roughly speaking the "promising" sets S are those for which x(γ(S)) is as large as possible. They are generated by using the following procedures. First, we generate the node sets of the connected components of Gu [x] and Gu [x]\{0}, respectively, and capacity constraints are checked. The shrinking procedure is then used. This procedure iteratively shrinks any edge of Gu [x] with weight 1, where shrinking an edge {i, j} means replacing nodes i and j by a new super-node, say vij , that is linked to any node u that was adjacent to i or j by an edge {vij , u} with weight xiu + xju . Every time a new super-node is created, the corresponding subset of original nodes is used to check capacity, precedence, LIFO, and same route constraints. If no violated constraint is generated by the above procedures, two general checking algorithms that work on the directed graph G[x] are called. The first checking algorithm works by iteratively enlarging a path that goes out of the depot, say (i0 = 0, i1 , . . . , ir ) such that xit it+1 = 1 for all t = 0, . . . , r − 1. Each time a new arc, say (is−1 , is ), is added to the path the set {i1 , . . . , is } is used to check the capacity, LIFO and same route constraints. The second checking algorithm is similar but it uses paths entering the depot. Finally, if no violated constraint has been found, the growing sets procedure is called. This procedure starts with an initial subset S = {i}, i ∈ V \{0}, and iteratively adds to it a new node j ∈ / S ∪{0} such that x(S : j) > 0 is as large as possible, until no node can be added. Every generated subset is used to check the capacity, precedence, LIFO and same route constraints. About 90% of the nodes are randomly chosen to be used as the initial subset in this procedure. 5. Heuristic approach

AC

In this section the heuristic approach used to find feasible solutions for the problem, and thus upper bounds, within a short computational time, will be described. It is based on a tabu search metaheuristic, following [12], that is embedded into a multistart method. Several initial solutions are first built by different simple constructive methods, and then the tabu search improving heuristic is applied to each of them. 5.1. Constructive algorithms The following three constructive algorithms are proposed: 14

ACCEPTED MANUSCRIPT

• Clarke and Wright:

CR IP T

This method is an adaptation of the algorithm proposed in [7] for the Capacitated Vehicle Routing Problem. Initially, one vehicle is assigned to each request i, following a route of the form (0, i, n + i, 0). At each iteration, two routes are merged into one single route in order to improve the overall cost of the solution, meeting LIFO and maximum capacity and duration constraints at all times. The savings achieved by each possible merge are calculated to decide which pair of routes are merged in the next iteration. This process is iterated until no more routes can be merged. This algorithm generates one single solution with a reasonably good quality.

• Random solution with consecutive pickups and deliveries:

AN US

Initially, all requests are unrouted and one route (the current one) is initialized with the depot node (0). At each iteration, one unrouted request i is selected at random. If feasible, the sequence (i, n + i) is added at the end of the current route; otherwise, the current route is closed, returning to the depot, and a new route is initialized with the sequence (0, i, n + i).

This process is iterated until all requests are routed. If it is run several times, this algorithm may generate multiple solutions, where each solution contains routes of the form (0, i, n + i, j, n + j, · · · , k, n + k, 0). • Random solution with pickups before deliveries:

PT

ED

M

Initially, all requests are unrouted and one route (the current one) is initialized with the depot node (0). At each iteration, one unrouted request i is selected at random. If feasible, the sequence (i, n + i) is inserted on the current route after the last pickup node and before its corresponding delivery node; if this is not possible (due to capacity or maximum duration constraints), the sequence (i, n + i) is added at the end of the route, if feasible; otherwise, the current route is closed, returning to the depot, and a new route is initialized with the sequence (0, i, n + i).

CE

This process is iterated until all requests are routed. This algorithm may generate multiple solutions whose routes are of the form (0, i, j, k, · · · , n + k, n + j, n + i, · · · , 0).

5.2. Improving heuristic: tabu search

AC

In order to improve the quality of the solutions built by the constructive algorithms proposed in the previous section, an improving heuristic based on a tabu search scheme will be applied to each initial solution generated. Tabu search has been previously applied successfully to similar problems. See, for example, [9], where it is used to solve a vehicle routing problem with time windows, or [17], where it is applied to a ship routing and scheduling problem. The tabu search framework is used to introduce diversification in the search, so that it is able to escape from local optima. Still, it includes a local search process that depends on the definition of a neighborhood for the solutions. In our case, the 15

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

neighborhood of a solution S is formed by the solutions that can be obtained by removing the pickup and delivery nodes of one request from S and reinserting them at different positions, on the same or on a different route. The search is started from an initial solution built by using any of the procedures introduced in the previous section, and at each iteration a new solution belonging to the neighborhood of the current solution is selected. The move to be performed is chosen mostly based on the cost of the new solution obtained, but other indicators aimed to diversify the search and allow the escaping from local optima are also considered. Following this idea, when one request c is moved from one route r1 to another route r2 , this move, together with the reverse move, are considered tabu during a certain number of iterations. As a result, none of these two moves can be implemented during the following iterations, unless one of them leads to a solution that improves the best solution found so far, in which case the tabu status is revoked. Furthermore, diversification is encouraged by adding a penalty to the evaluation function that is proportional to the number of times the considered move was performed in past iterations, thus penalizing frequently used moves. The tabu search process described above is based on removal and insertion moves, which usually lead to solutions in which one route visits two more nodes (the pickup and delivery nodes of the request that is moved). In situations where the maximum duration constraint is tight, there may exist very few moves (if any) verifying this constraint, so the search process is very restricted and it is very likely to get stuck in a local optimum. This could be solved by considering additional moves based on interchanges of two requests belonging to two different routes; however, the number of possible moves of this type is much higher, thus increasing too much the size of the neighborhood, and complicating the evaluation of the moves. As a result, we decided not to use interchange moves and expand the search space by temporarily relaxing the maximum duration constraints and penalizing their violation when evaluating the candidate moves. In this way we introduce diversification by allowing intermediate infeasible solutions, but at the same time the search is guided towards feasible solutions. Candidate moves leading to a solution S are thus evaluated according to a function of the form C(S) + λV (S) + µD(S). C(S) is the cost of S. V (S) is a measure related to the violation of the maximum duration constraint, calculated as the sum of the amount of time exceeding the maximum duration limit over all routes. D(S) is the diversification penalty applied to the candidate move and is calculated as the number of times one insertion has been performed at the same position of the same route as in the candidate move during the execution of the algorithm divided by the number of iterations performed up to that moment. It represents the proportion of accepted moves that are similar to the candidate move, aiming to diversify the search by penalizing moves performed repeatedly. The diversification intensity µ is kept constant during the execution of the algorithm, but λ is dynamically updated during the search: if the current solution is unfeasible, λ is increased in order to give more importance to the verification of the maximum duration constraint; otherwise, it is decreased to increase diversification. This mechanism is aimed to guide the search towards feasible solutions but at the same time increase its flexibility by allowing certain constraint violations. Note also 16

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

that the calculation of the cost of the new solution obtained after performing one move does not require the re-evaluation of the whole solution, because the cost saving achieved by the move can be calculated efficiently in constant time by performing a few additions and subtractions. As stated above, tabu search includes several mechanisms to diversify the search in order to escape from local optima. However, even using all these mechanisms, the search might get trapped inside a small neighborhood in which different solutions can actually be attained, but they are all very similar, not allowing other promising regions of the solution space to be visited. To avoid this, when the relative variation of the cost of the solutions generated during a certain number of iterations is less than a certain value ε, a perturbation operator is applied to move the search to a different region of the solution space. This perturbation operator consists on performing several random consecutive moves on the solution, probably increasing its cost but also diversifying the search. Furthermore, to avoid reversing them too early after the perturbation phase, the performed random moves are added to the tabu list, and their tabu status is maintained for twice as many iterations as for the moves performed during the standard tabu search. After some preliminary experiments, and also following [25] and [17], the values of the main parameters used for the tabu search are the following: the number of iterations that a move is considered tabu is d7.5 log10 (n)e, where n is the number of customer requests; the diversification intensity is set to µ = 0.025; the weight λ related to the maximum duration violation is given initially value 1 and increased or decreased dynamically by a factor 0.5; finally, if for 300 consecutive iterations the cost of the solutions found has not varied at least ε = 0.01%, 10 perturbation moves are performed to further diversify the search.

ED

5.3. Multi-start framework

PT

The heuristic proposed follows a multi-start framework in which several initial solutions generated using the methods introduced in Section 5.1 are improved by the Tabu Search procedure described in Section 5.2. A schematic pseudocode of this heuristic procedure is shown in Algorithm 1, where Tabu(S, t) represents a call to the Tabu Search method that performs t iterations using S as the initial solution.

CE

Algorithm 1. Heuristic procedure

AC

Input • nsols: Number of random initial solutions to be generated. • itCW : Number of tabu search iterations from the C-W solution. • itR: Number of tabu search iterations from each initial random solution. • itB: Number of tabu search iterations from the best solution. Output • S: Feasible solution produced. Pseudocode 1. Initialize: Set i = 0. 2. Clarke and Wright: Generate an initial solution S using the adaptation of the Clark and Wright method.

17

ACCEPTED MANUSCRIPT

CR IP T

3. Improve: S¯ = Tabu(S, itCW ). 4. Initial solution: Set i = i + 1. If i is odd, generate a new initial solution S with consecutive pickups and deliveries. Otherwise, generate a new initial solution S with pickups before deliveries. 5. Improve: S˘ = Tabu(S, itR). ˘ ¯ update best S¯ = S. ˘ 6. Update best: If cost(S)
Some preliminary experiments were done to tune the parameters of the multistart framework with the aim of achieving a trade-off between computational effort and solution quality. After those experiments, the following parameter values were selected: nsols = 100, itCW = itR = 3000 and itB = 20000.

AN US

6. Computational results

In this section the computational study aimed to analyze the performance of the solution methods proposed in previous sections is presented. The algorithms have been implemented in C++ and executed on a 2.66GHz Core2Quad processor with 3Gb RAM running Windows XP. 6.1. Benchmark

AC

CE

PT

ED

M

To the best of the authors knowledge, it does not exist any benchmark available in the literature for the PDPLT. Then, we created a new set of test instances for this problem based on the Li & Lim [18] benchmark for the PDPTW that consists of 56 instances with about 50 customer requests each and 60 instances with about 200 customer requests each. Li & Lim instances were generated from the 56 Solomons’ benchmark instances for the Vehicle Routing Problem with Time Windows ([24]). The instances are divided into six classes lc1, lc2, lr1, lr2, lrc1 and lrc2, where the nodes in lc instances are in clusters, in lr instances the nodes are randomly distributed, and in lrc instances the nodes are partially clustered and partially randomly distributed. The lc1, lr1 and lrc1 problems have a narrow time window associated to the depot, while lc2, lr2 and lrc2 have a wider time window. We modified the Li & Lim instances by setting the maximum route duration equal to the width of the time window associated to the depot and dropping all other time windows. Furthermore, since the smallest Li & Lim instances had about 50 customer requests, we have generated 5 additional sets of smaller instances by selecting randomly 10, 15, 20, 25 and 30 customer requests, respectively, from each Li & Lim instance with about 50 requests. The resulting instances have 20, 30, 40, 50 and 60 nodes (excluding the depot) respectively. Hence, for our experiments we have used a total of 396 instances. The original PDPTW instances are available at [19] and the new PDPLT instances generated by us are available at [2].

18

ACCEPTED MANUSCRIPT

6.2. Main results

AC

CE

PT

ED

M

AN US

CR IP T

This section is devoted to the analysis and comparison of the performance of the compact formulation, the branch-and-cut method based on the non-compact formulation and the heuristic procedure. Table 1 summarizes the results obtained by these three approaches on the instances with 20 to 60 nodes (for a total of 280 instances). To simplify the presentation of the results and improve the readability of the table, only average results are presented for each group of instances according to their size (20, 30, 40, 50 and 60 nodes) and their type (lc1, lc2, lr1, lr2, lrc1 and lrc2), leading to a total of 30 groups. The first two columns show the name and the total number of instances of each group, respectively, and the rest of the table is divided in three parts, one associated to each solution method. The heuristic part contains five columns: the first one shows the average number of vehicles used in the heuristic solution (#veh), the average and maximum gaps are given in the following two columns (av. gap and max. gap), the number of instances for which the optimal solution was found (#opt) is given in the fourth column and the last column shows the running time in seconds. The gaps for the heuristic method are calculated as 100(UBH-UBB)/UBB , where UBH is the cost of the solution provided by the heuristic and UBB the cost of the best known feasible solution. The other two parts of the table (columns 8 to 15) show, for each exact method and for each group of instances, the average and maximum gap (columns av. gap and max. gap), the number of instances solved to optimality (#opt) and the average running time in seconds (time). Gaps for the exact methods are calculated as 100(UBBLB)/UBB, where LB denotes the best lower bound provided by the respective exact method. The exact methods were run with a time limit of 3 hours. If this limit is reached before the branch-and-cut tree was completed when solving an instance, then the method provides a lower bound that is equal to the minimum of the lower bounds among all the open nodes and the best feasible solution found so far. The average running time was computed considering only those instances for which the exact method was completed. An ** indicates that no instance in the group could be solved to optimality. Since the compact formulation was not able to solve any of the lc instances with 50 nodes, nor any lr or lcr instances with 40 nodes, it was not applied to larger instances, which is indicated by — in the table. As it can be observed in Table 1, the heuristic method is able to obtain good quality solutions with short running times, quickly providing good upper bounds to be used by the exact methods. For the instances that were optimally solved by any of the exact methods, the heuristic procedure has reached a solution that is optimal or very close to the optimal. Concerning the exact methods, the branch-and-cut algorithm based on the non-compact formulation shows a better overall performance than the compact formulation, providing more exact solutions and smaller gaps with less computational effort. It is also noticeable how lr and lcr instances seem harder to solve than lc instances. All but one of the lc instances with up to 40 nodes were solved to optimality by the branch-and-cut procedure, as well as 7 instances with 50 nodes and 5 instances with 60 nodes. However, none of the lr or lcr instances with 60 nodes could be solved to optimality. In general, the average gaps for lr and lcr instances are

19

name lc1_20 lc2_20 lc1_30 lc2_30 lc1_40 lc2_40 lc1_50 lc2_50 lc1_60 lc2_60 lr1_20 lr2_20 lr1_30 lr2_30 lr1_40 lr2_40 lr1_50 lr2_50 lr1_60 lr2_60 lrc1_20 lrc2_20 lrc1_30 lrc2_30 lrc1_40 lrc2_40 lrc1_50 lrc2_50 lrc1_60 lrc2_60

#inst 9 8 9 8 9 8 9 8 9 8 12 11 12 11 12 11 12 11 12 11 8 8 8 8 8 8 8 8 8 8

AC #veh 2.00 1.00 3.00 1.00 3.78 2.00 4.56 2.00 5.00 2.25 2.83 1.00 4.00 1.00 4.67 1.00 5.50 2.00 6.25 2.00 3.00 1.00 3.88 1.00 4.88 1.75 5.75 2.00 6.50 2.00

CE ED

Heuristic max. gap % #opt 0.00 9 0.00 8 0.00 9 0.05 7 0.38 6 0.00 8 0.00 0 0.00 7 0.00 0 0.00 5 0.00 12 0.00 11 1.55 10 0.00 10 0.00 2 0.14 5 0.00 0 0.00 0 0.00 0 0.00 0 0.00 8 0.00 8 0.00 6 0.00 8 0.00 2 0.00 4 0.00 2 0.00 1 0.00 0 0.00 0

PT

av. gap % 0.00 0.00 0.00 0.01 0.05 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.13 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

M

time 0.54 0.22 1.23 0.86 870.48 3.64 ** 299.30 ** 2970.37 5.69 4.39 565.14 251.37 2014.24 3480.93 ** ** ** ** 1.61 1.25 1202.64 93.51 1016.27 1678.01 3120.05 107.70 ** **

Compact Formulation av. max. gap % gap % #opt time 0.00 0.00 9 29.00 0.00 0.00 8 3.00 0.00 0.00 9 910.00 0.00 0.00 8 108.00 0.21 1.13 0 ** 0.00 0.00 8 3338.00 0.56 1.13 0 ** 0.30 0.50 0 ** — — — — — — — — 0.00 0.00 12 60.00 0.00 0.00 11 62.00 0.00 0.00 8 1391.00 0.04 0.13 9 1647.00 0.77 1.65 0 ** 4.18 13.41 0 ** — — — — — — — — — — — — — — — — 0.00 0.00 8 48.00 0.00 0.00 8 23.00 0.54 1.83 3 1562.00 0.00 0.00 8 2441.00 3.26 17.60 0 ** 1.82 3.24 0 ** — — — — — — — — — — — — — — — —

CR IP T

Branch-and-cut max. gap % #opt 0.00 9 0.00 8 0.00 9 0.00 8 0.26 8 0.00 8 0.65 0 0.46 7 0.54 0 0.44 5 0.00 12 0.00 11 2.09 11 1.68 10 6.00 2 2.86 7 5.59 0 8.47 0 5.71 0 9.72 0 0.00 8 0.00 8 6.95 6 0.00 8 4.27 2 9.39 4 10.51 2 7.36 1 7.82 0 12.38 0

AN US

time 12.98 24.11 17.72 40.41 21.17 50.16 26.58 62.94 31.05 72.31 8.62 22.71 14.21 42.12 18.94 60.47 24.82 79.52 28.63 104.03 6.86 21.51 12.06 39.13 16.43 50.81 0.00 67.20 26.71 84.67

av. gap % 0.00 0.00 0.00 0.00 0.03 0.00 0.46 0.06 0.25 0.10 0.00 0.00 0.17 0.15 1.96 0.80 3.36 4.90 3.67 6.26 0.00 0.00 1.14 0.00 2.27 2.80 4.42 3.62 5.11 7.13

Table 1: Results for instances with 20, 30, 40, 50 and 60 nodes

ACCEPTED MANUSCRIPT

20

ACCEPTED MANUSCRIPT

CE

PT

ED

M

AN US

CR IP T

significantly higher than the ones for lc instances, which shows that random instances seem to be harder than clustered ones. Recall that instances of groups lc1, lr1, lcr1 have tighter maximum route duration than instances of groups lc2, lr2, lcr2, which implies that a larger number of vehicles is used. The results shown in Table 1 seem to indicate that instances with more vehicles are harder. Since solving to optimality instances with more than 60 nodes was intractable, we decided to apply only the heuristic to the larger instances generated by Li & Lim and use the branch-and-cut method just to provide a valid lower bound. These instances are grouped in two sets: the ones belonging to the first set have between 100 and 110 nodes and the ones in the second set between 402 and 422 nodes (in both cases excluding the depot). The results for each instance are presented in Tables 2 and 3, respectively. These tables show, for each instance, its name, the cost of the heuristic solution (ub), the number of vehicles (#veh), the gap with respect to the lower bound (gap), and the running time in seconds required by the heuristic (time). The average gap and time are given in the last row of the table. In Table 2 it can be observed that lc instances, with an average gap of 0.89%, seem to be easier than lr and lrc ones, with average gaps of 10.51% and 13.41%, respectively. The results of this table show that instances with 100-110 nodes have been managed properly by the heuristic procedure within short computational times, and thanks to the branch-and-cut we have a guarantee of the quality of the solutions, that in the case of the lc instances are almost optimal. Note that for this set, instances with a larger number of vehicles show smaller running times and smaller gaps. The results presented in Table 3 for instances with more than 400 nodes exhibit a similar pattern, although the running times and gaps are obviously larger than those for the 100 nodes instances. Note that lc instances show gaps that are about 4%, which are quite small taking into account the large size of these instances. The closest variant to the PDPLT in the literature on pickup and delivery problems is the recent paper by Cherkesly et al. [6] where time windows are associated to every pickup and delivery location, while no maximum duration time for the routes is explicitly considered. They succeeded in solving to optimality instances involving up to 75 requests (150 nodes) with a branch-price-and-cut algorithm based on solving constrained shortest path subproblems. The instances used in [6] are very different to the ones we have used because the existence of time windows drastically reduces the number of feasible routes, and thus a direct comparison of the performance of both methods is not possible.

AC

7. Conclusions and future reseach This paper deals with a particular version of the pickup and delivery problem with multiple vehicles that has not been considered up to now. In our problem, loading and unloading operations are subject to LIFO conditions and maximum duration constraints are applied to the routes. In this work we propose two mixed integer formulations for the problem: a compact formulation, in which the number of variables and constraints grows polynomially with the number of requests, and a non-compact

21

CR IP T

ACCEPTED MANUSCRIPT

Table 2: Results for instances with 100-110 nodes

PT

2289.98 2336.21 2213.32 2191.28 2335.36 2238.28 2218.77 2201.42

time 65.91 59.19 56.47 62.22 54.31 60.28 62.25 60.89 52.73 59.36 61.36 62.24 49.53 47.72 51.80 58.52 52.13 45.59 60.38 54.05 60.55 53.89 54.81 50.88 48.58 61.28 57.22 60.52 57.77 50.31 51.73 54.79

inst lc201 lc202 lc203 lc204 lc205 lc206 lc207 lc208

ub 9825.64 9860.47 9873.71 9832.31 9794.84 9882.15 9802.90 9803.00

#veh 3 3 3 3 3 3 3 3

gap % 0.83 1.24 1.34 0.89 0.49 1.17 0.58 0.58

time 125.00 115.66 107.53 121.44 103.66 116.55 98.61 99.86

lr201 lr202 lr203 lr204 lr205 lr206 lr207 lr208 lr209 lr210 lr211

1967.24 2130.35 2102.09 2184.12 2079.14 2111.64 2146.86 2090.03 2066.29 2073.73 2027.38

3.00 3 3 3 3 3 3 3 3 3 3 3

0.89 8.46 13.91 14.13 16.06 12.24 13.70 15.67 12.64 11.92 12.83 11.42

111.04 188.69 188.00 208.83 189.06 202.53 168.78 188.17 170.94 211.49 193.95 216.20

lrc201 lrc202 lrc203 lrc204 lrc205 lrc206 lrc207 lrc208

2203.00 2250.50 2192.66 2018.41 2261.11 2281.59 2426.00 2176.96

3.00 3 3 3 3 3 3 4 3 3.13

13.00 16.11 16.32 15.45 8.30 17.81 16.93 22.41 14.59 15.99

193.33 161.02 158.69 157.36 158.92 171.17 145.84 169.05 154.55 159.57

AN US

gap % 0.94 0.86 0.88 0.83 0.83 0.97 0.94 0.90 0.79 0.88 9.17 9.21 8.55 6.89 8.09 8.92 9.14 8.22 9.59 5.97 8.02 6.93 8.23 14.75 11.97 10.47 8.80 13.10 10.13 9.62 7.82 10.83

M

2140.15 2125.61 2113.11 2072.80 2114.97 2147.80 2188.71 2150.33 2165.48 2041.50 2114.81 2102.96

#veh 10 9 9 10 9 9 9 9 9 9.22 11 10 10 10 10 10 11 10 11 10 10 10 10.25 10 10 11 10 11 10 10 10 10.25

ED

ub 9867.61 9862.70 9867.84 9857.95 9848.78 9866.18 9874.48 9862.42 9855.31

AC

CE

inst lc101 lc102 lc103 lc104 lc105 lc106 lc107 lc108 lc109 avg. lr101 lr102 lr103 lr104 lr105 lr106 lr107 lr108 lr109 lr110 lr111 lr112 avg. lrc101 lrc102 lrc103 lrc104 lrc105 lrc106 lrc107 lrc108 avg.

22

CR IP T

ACCEPTED MANUSCRIPT

Table 3: Results for instances with 402-422 nodes

10502.64 10624.45 10515.92 9910.92 10190.46 10414.40 10321.30 9290.86 10447.86 10041.44 9700.01 9651.76 9667.97 9015.70 9555.94 9363.88 9645.08 9528.92 9634.73 9639.88

time 405.08 382.77 402.63 417.38 395.00 437.97 391.61 400.83 377.06 406.77 401.71 840.63 827.59 805.61 850.89 762.41 829.00 851.63 856.17 827.13 848.27 829.93 818.39 838.09 759.05 772.53 791.08 814.67 870.91 829.53 837.19 838.55 817.00

AC

CE

PT

LRC1_4_1 LRC1_4_2 LRC1_4_3 LRC1_4_4 LRC1_4_5 LRC1_4_6 LRC1_4_7 LRC1_4_8 LRC1_4_9 LRC1_4_10

gap % 2.72 2.77 3.21 2.89 2.68 2.62 2.55 2.75 3.36 3.07 2.86 23.80 23.72 23.56 21.26 21.90 23.20 23.79 16.57 24.16 20.42 22.24 22.24 21.76 21.08 16.39 20.69 19.43 21.65 21.45 22.72 21.66 20.91

inst LC2_4_1 LC2_4_2 LC2_4_3 LC2_4_4 LC2_4_5 LC2_4_6 LC2_4_7 LC2_4_8 LC2_4_9 LC2_4_10

ub 40714.04 41140.66 41316.85 41035.84 40996.57 40727.48 41001.45 40925.04 40810.30 41060.15

AN US

LR1_4_1 LR1_4_2 LR1_4_3 LR1_4_4 LR1_4_5 LR1_4_6 LR1_4_7 LR1_4_8 LR1_4_9 LR1_4_10

#veh 30 30 30 30 30 29 30 30 30 30 29.90 18 18 18 17 16 16 20 15 18 17 17.30 17 16 16 14 16 17 18 16 17 16 16.30

M

ub 42336.75 42366.74 42599.50 42482.09 42312.18 42284.78 42257.38 42409.92 42657.31 42539.79

ED

inst LC1_4_1 LC1_4_2 LC1_4_3 LC1_4_4 LC1_4_5 LC1_4_6 LC1_4_7 LC1_4_8 LC1_4_9 LC1_4_10

23

LR2_4_1 LR2_4_2 LR2_4_3 LR2_4_4 LR2_4_5 LR2_4_6 LR2_4_7 LR2_4_8 LR2_4_9 LR2_4_10

12095.03 12087.73 11890.15 11445.21 11753.71 11549.06 11454.11 11684.98 11695.28 12228.19

LRC2_4_1 LRC2_4_2 LRC2_4_3 LRC2_4_4 LRC2_4_5 LRC2_4_6 LRC2_4_7 LRC2_4_8 LRC2_4_9 LRC2_4_10

10733.52 10444.30 10704.17 10857.97 10288.07 10596.35 11160.72 10882.58 11054.00 10855.83

#veh 12 14 13 13 13 12 13 13 13 14 13.00 7 10 7 10 6 6 8 6 9 8 7.70 6 5 7 12 7 10 7 11 7 7 7.90

gap % 3.20 4.07 4.50 3.85 3.95 3.33 3.99 3.83 3.63 3.97 3.83 33.73 34.23 33.49 31.47 32.41 31.80 31.37 33.14 32.89 34.87 32.94 32.85 31.01 31.45 32.57 29.31 31.92 35.39 32.97 33.85 34.00 32.53

time 732.05 737.39 741.09 764.05 772.70 744.81 752.27 737.17 788.56 735.99 750.61 1598.11 1528.31 1578.91 1612.69 1622.72 1540.38 1658.30 1472.44 1574.20 1589.20 1577.53 1617.91 1670.25 1577.23 1489.13 1612.36 1638.67 1551.36 1478.17 1578.86 1590.75 1580.47

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

formulation with an exponential number of constraints. Several sets of additional valid inequalities that improve the original formulations are proposed. New separation procedures were designed and implemented for those families of constraints with exponential size. These separation algorithms were used to implement an efficient branch-and-cut algorithm based on the non-compact formulation. Besides, a metaheuristic procedure based on a tabu search process embedded in a multi-start framework is also proposed. The performance of these solution methods was evaluated through extensive computational experiments performed on randomly generated instances with different characteristics and sizes. Overall, the branch-and-cut method based on the non-compact formulation exhibits a better performance than the method based on the compact one, and it was able to solve instances with up to 60 nodes (30 requests). The heuristic procedure reached the optimal solution in most of the instances for which it is known. When applied to larger instances (with more than 200 requests) the heuristic procedure requires relatively short running times and obtains good solutions, as it has been guaranteed by the lower bound given by the branch-and-cut method. Nevertheless, further research is needed to consistently solve medium size instances (with about 60 nodes). In order to achieve this, new valid cuts are needed to strengthen the non-compact formulation, as well as better branching strategies. The development of alternative exact methods, for example based on column generation techniques, could be also interesting ideas for future work in this problem, as well as the design of hybrid heuristics based on other metaheuristics that also worked successfully for other vehicle routing problems, such as variable neighborhood search or simulated annealing.

PT

ED

Acknowledgements.. This work was partially supported by projects by the Generalitat Valenciana and the Government of Spain: projects GVPROMETEO2013049, MTM2009-14039-C06-02 and MTM2012-36163-C06-02, (E. Benavent and E. Mota), projects MTM2009-14039-C06-03 and TIN2012-32482 (G. Tirado) and project MTM2012-36163-C06-04, RDEF funds and Fundación Séneca 08716/PI/08 (M. Landete). Project MTM2012-36163-C06-02 was co-funded by FEDER. References

CE

[1] Battarra M, Erdogan G, Laporte G, Vigo D (2010) The traveling salesman problem with pickups, deliveries, and handling costs, Transportation Science 44: 383399.

AC

[2] Benavent, Mota, Landete and Tirado PDPLT benchmark (2014). Last accessed: 10-11-2014. http://www.mat.ucm.es/∼gregoriotd/PDPLT.htm [3] Berbeglia G, Cordeau JF, Gribkovskaia I, Laporte G (2007) Static pickup and delivery problems: A classifcation scheme and survey, TOP 15(1): 1-31. [4] Carrabs F, Cerulli R, Cordeau JF (2007) An additive branch-and-bound algorithm for the pickup and delivery traveling salesman problem with LIFO or FIFO loading, INFOR 45: 223-238. 24

ACCEPTED MANUSCRIPT

[5] Cheang B, Gao X, Lim A, Qin H, Zhu W (2012) Multiple pickup and delivery traveling salesman problem with last-in-first-out loading and distance constraints, European Journal of Operational Research 223: 60-75.

CR IP T

[6] Cherkesly M, Desaulniers G, Laporte G (2014) Branch-Price-and-Cut Algorithms for the Pickup and Delivery Problem with Time Windows and Last-in-First-Out Loading, Transportation Science, to appear.

[7] Clark G, Wright JW (1964) Scheduling of Vehicles from a Central Depot to a Number of Delivery Points, Operations Research 12(4): 568-581. [8] Cordeau JF, Iori M, Laporte G, Salazar JJ (2010) A Branch-and-Cut Algorithm for the Pickup and Delivery Traveling Salesman Problem with LIFO Loading, Networks 55: 46-59.

AN US

[9] Cordeau JF, G Laporte, A Mercier (2001) A unified tabu search heuristic for vehicle routing problems with time windows, Journal of the Operational Research Society 52: 928-936. [10] Erdogan G, Battarra M, Laporte G, Vigo D (2012) Metaheuristics for the traveling salesman problem with pickups deliveries and handling costs, Computers and Operations Research 39: 1074-1086.

M

[11] Gao X, Lim A, Qin H, Zhu W (2011) Multiple pickup and delivery TSP with LIFO and distance constraints: a VNS approach, Lecture Notes in Computer Science 6704: 193-202.

ED

[12] Glover F (1989) Tabu Search - Part I, INFORMS Journal on Computing 1(3): 190-206. [13] Golden BL, Raghavan S, Wasil EA (2008) The Vehicle Routing Problem: Latest Advances and New Challenges. Springer.

PT

[14] IBM-ILOG (2011) CPLEX 12.4. User’s Manual,

CE

[15] Kara I, G Laporte, T Bektas (2004) A note on the lifted Miller-Tucker-Zemlin subtour elimination constrints for the capacitated vehicle routing problem, European Journal of Operational Research 158: 793-795. [16] Kara I (2008) Two indexed polynomial size formulation for vehicle routing problems, Working paper.

AC

[17] Korsvik JE, K Fagerholt, G Laporte (2010) A tabu search heuristic for ship routing and scheduling, Journal of the Operational Research Society 61: 594-603. [18] Li H, Lim A (2001) A MetaHeuristic for the Pickup and Delivery Problem with Time Windows, In Proceedings of the 13th International Conference on Tools with Artificial Intelligence, Dallas, TX, USA. [19] Li & Lim benchmark (2013). Last accessed: 10-4-2014. http://www.sintef.no/Projectweb/TOP/PDPTW/Li--Lim-benchmark/ 25

ACCEPTED MANUSCRIPT

[20] Miller CE, AW Tucker, RA Zemlin (1960) Integer programming formulation of travelling salesman problems, Journal of the Association for Computer Machinery 7: 326-329.

CR IP T

[21] Parragh S, Doerner KF, Hartl RF (2008a) A survey on pickup and delivery problems, Part I: Transportation between customers and depot, Journal für Betriebswirtschaft 58: 21-51. [22] Parragh S, Doerner KF, Hartl RF (2008b). A survey on pickup and delivery problems, Part II: Transportation between pickup and delivery locations, Journal für Betriebswirtschaft 58: 81-117.

[23] Ropke S, Cordeau JF (2009) Branch-and-cut-and-price for the pickup and delivery problem with time windows, Transportation Science 43: 267-286.

AN US

[24] Solomon MM, (1987) Algorithms for the vehicle routing and scheduling problems with time window constrains, Operations Research 35: 254-264.

[25] Tirado G, LM Hvattum, K Fagerholt, JF Cordeau (2013) Heuristics for dynamic and stochastic routing in industrial shipping, Computers and Operations Research 40: 253-263.

AC

CE

PT

ED

M

[26] Toth P, Vigo D (2002) The Vehicle Routing Problem. SIAM.

26

ACCEPTED MANUSCRIPT

Appendix

−ujk + xkj − βjk xjk ≤ −βjk uik + xik − βik xki ≤ 1

−uj−n,k + xk+n,j + (2 − βjk )xj,k+n ≤ 1 − βjk ui−n,k + xi,k+n + (2 − βik )xk+n,i ≤ 2

(43)

∀i ∈ P, k ∈ P : k 6= i

∀i, j ∈ P, i 6= j, k ∈ P : k 6= i, j

∀j ∈ D, k ∈ P : k 6= j − n

uik + xi,k+n + (2 − βik )xk+n,i ≤ 2

uik − uj−nk + xij + (1 − βik − βjk )xji ≤ 1 − βjk −ujk + xk+n,j + (2 − βjk )xj,k+n ≤ 1 − βjk

ui−n,k − ujk + xij + (1 − βik − βjk )xji ≤ 1 − βjk

∀j ∈ D, k ∈ P : k 6= j − n

(48)

∀i ∈ P, j ∈ D, k ∈ P : k 6= i, j − n

(50)

∀i ∈ D, k ∈ P : k 6= i − n

(52)

∀i, j ∈ D, i 6= j, k ∈ P : k 6= i − n, j − n (47) ∀i ∈ P, k ∈ P : k 6= i

(49)

∀j ∈ P, k ∈ P : k 6= j

(51)

∀i ∈ D, j ∈ P, k ∈ P : k 6= j, i − n

(53)

(54)

M

ED P

k∈P :k6=i

dk uik ≤ Q − di

PT

CE

AC

27

(55)

(56)

∀i ∈ P

Finally, (57). . . (60) may replace (24), and (61) and (62) may replace (23), P k∈P :k6=i,j dk (uik − ujk ) + Qxij + (Q − di − dj )xji ≤ Q − di P dk (ui−n,k − uj−n,k ) + Qxij + (Q − di − dj )xji ≤ Q − dj Pk∈P :k6=i−n,j−n k∈P :k6=i,j−n dk (uik − uj−n,k ) + Qxij + (Q − di − dj )xji ≤ Q − dj − di P k∈P :k6=i−n,j dk (ui−n,k − ujk ) + Qxij + (Q − di − dj )xji ≤ Q P P k∈K:k6=i dk uik − j6=i dj xji ≥ 0 P P k∈K:k6=i−n dk ui−n,k − j6=i dj xji ≥ di Thus, the proposed model reads as follows:

(45) (46)

Analogously, constraints (54) and (55) may replace (6) P uik − j6=i βjk xji ≥ βik ∀i, k ∈ P : i 6= k P ui−n,k − j6=i βjk xji ≥ βik ∀i ∈ D, k ∈ P : k 6= i − n and (56) may replace (10),

(44)

∀i ∈ D, k ∈ P : k 6= i − n

AN US

ui−n,k − uj−n,k + xij + (1 − βik − βjk )xji ≤ 1 − βjk

ui−n,k + xik − βik xki ≤ 1

(42)

∀j ∈ P, k ∈ P : k 6= j

uik − ujk + xij + (1 − βik − βjk )xji ≤ 1 − βjk

−uj−n,k + xkj − βjk xjk ≤ −βjk

CR IP T

Constraints (7), (8) and (9) imply that we can replace ui+n,k by ui,k when k 6= i and ukk , un+k,k by their corresponding values, leading to a formulation with fewer variables. In particular, constraint (5) may be replaced by (42). . . (53).

∀i, j ∈ P, i 6= j

(57)

∀i ∈ P, j ∈ D

(59)

∀i ∈ P

(61)

∀i, j ∈ D, i 6= j (58) ∀i ∈ D, j ∈ P

(60)

∀i ∈ D

(62)

ACCEPTED MANUSCRIPT

min

(i,j)∈A (cij

+ si )xij

(11) − (17), (21), (22), (25) − (28), (42) − (62) αij ≥ 0 ∀(i, j) ∈ A

uik ∈ {0, 1} ∀i, k ∈ P : i 6= k

AC

CE

PT

ED

M

AN US

xij ∈ {0, 1} ∀(i, j) ∈ A

CR IP T

s.t.

P

28