General solutions to the single vehicle routing problem with pickups and deliveries

European Journal of Operational Research 180 (2007) 568–584 www.elsevier.com/locate/ejor

Discrete Optimization

General solutions to the single vehicle routing problem with pickups and deliveries Irina Gribkovskaia

a,*

, Øyvind Halskau sr. a, Gilbert Laporte b, Martin Vlcˇek

a

a

b

Molde University College, Post Box 2110, N-6402 Molde, Norway Canada Research Chair in Distribution Management, HEC Montre´al, 3000 chemin de la Coˆte-Sainte-Catherine, Montre´al, Canada H3T 2A7 Received 30 October 2005; accepted 8 May 2006 Available online 7 July 2006

Abstract The single vehicle routing problem with pickups and deliveries (SVRPPD) is defined on a graph in which pickup and delivery demands are associated with the customer vertices. The problem consists of designing a least cost route for a vehicle of capacity Q. Each customer is allowed to be visited once for a combined pickup and delivery, or twice if these two operations are performed separately. This article proposes a mixed integer linear programming model for the SVRPPD. It introduces the concept of general solution which encompasses known solution shapes such as Hamiltonian, double-path and lasso. Classical construction and improvement heuristics, as well as a tabu search heuristic, are developed and tested over several instances. Computational results show that the best solutions generated by the heuristics are frequently nonHamiltonian and may contain up to two customers visited twice.  2006 Elsevier B.V. All rights reserved. Keywords: Combinatorial optimization; Heuristics; Metaheuristics; Tabu search

1. Introduction The single vehicle routing problem with pickups and deliveries (SVRPPD) is defined as follows. Let G = (V, A) be a graph where V = {0, 1, . . . , n} is the vertex set, A = {(i, j): i, j 2 V, i 5 j} is the arc set, and C = (cij)(n+1)·(n+1) is a cost matrix defined on G. Vertex 0 is a depot at which is based a vehicle of capacity Q, while the remaining vertices represent customers. With each vertex i 2 Vn{0} are associatedP a non-negative n pickup demand p and a non-negative delivery demand d , with p + d > 0. It is assumed that i i i i i¼1 p i 6 Q and Pn d 6 Q, for otherwise the problem is infeasible. In practice the products to be picked up are different from i i¼1 those delivered. The SVRPPD consists of designing a minimum cost vehicle route starting and ending at the depot, making all pickups and deliveries, and such that the vehicle load never exceeds Q along the route. We *

Corresponding author. Tel.: +47 712 14211; fax: +47 712 14100. E-mail addresses: [email protected] (I. Gribkovskaia), [email protected] (Ø. Halskau sr.), [email protected] (G. Laporte), [email protected] (M. Vlcˇek). 0377-2217/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2006.05.009

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

569

allow each customer to be visited once or twice. In the first case, both the pickup and the delivery are performed during the same visit. In the second case, one visit is used for the delivery and the other one for the pickup. However, neither the delivery demand nor the pickup demand can be split between two visits. The SVRPPD is NP-hard since it reduces to the traveling salesman problem (TSP) whenever pi P di for all i 2 Vn{0} or di P pi for all i 2 Vn{0}. The SVRPPD arises naturally in the planning of reverse logistics operations involving the delivery of full bottles and the collection of empty ones (Dethloff, 2001; Tang and Galva˜o, 2002, 2006). Real-life applications encountered in the beverage industry are described by Halskau and Løkketangen (1998) and Prive´ et al. (in press). Min (1989) has investigated an application related to the public library distribution system, while Mosheiov (1994) provides yet another application arising in the transportation of children in a community program. Applications encountered in mail transportation are described by Wasner and Za¨phel (2004). There exists an extensive literature on routing problems with pickups and deliveries, but more often than not these are ‘‘many-to-many’’ problems involving deliveries between customer locations. Such problems, which are beyond the scope of the present study, are reviewed in Cordeau et al. (in press). Our problem combines ‘‘one-to-many’’ and ‘‘many-to-one’’ features: all customer deliveries originate at a unique depot, and all collections are destined to the depot. The vehicle routing problem with pickups and deliveries (VRPPD) in which pickups and deliveries are made with a fleet of m vehicles based at the depot generalizes our problem. As far as we are aware, most published VRPPD models (see, e.g., Angelelli and Mansini, 2002; Nagy and Salhi, 2005) assume that each customer is visited only once, except for Nagy and Salhi (1998) who have developed a large scale model allowing two visits per customer for a more general many-to-many pickup and delivery locationrouting problem. Some heuristics (Gribkovskaia et al., 2001; Hoff and Løkketangen, in press; Nagy and Salhi, 2005) allow two visits to the same customer. As suggested by Halskau and Gribkovskaia (2002), it is preferable to design models and heuristics capable of constructing general solutions in the sense that some customers may be visited once while others may be visited twice, and no a priori shape is imposed on the solution. Most known algorithms for the SVRPPD and the VRPPD have been developed in contexts where the solution has a predefined shape. All are heuristics that include a construction phase followed by an improvement phase which moves customers between different routes or operates on a single route at a time. Typical movements include customer relocations in different positions, 2-opt and 3-opt moves (Min, 1989; Halse, 1992; Mosheiov, 1994; Salhi and Nagy, 1999; Gendreau et al., 1999; Dethloff, 2001, 2002; Tang and Galva˜o, 2002, 2006; Nagy and Salhi, 2005). Nagy and Salhi (2005) also consider the possibility of splitting customer visits into two, or of merging the two visits during the solution process, but this flexibility is not allowed by their problem formulation. Tabu search has been applied by Gendreau et al. (1999), and Hoff and Løkketangen (in press). Our aim is to provide new construction and improvement heuristics, including tabu search, yielding general solutions to the SVRPPD. The remainder of this paper is organized as follows. In Section 2, we present a mathematical model for the SVRPPD, followed by some properties of SVRPPD solutions in Section 3. Classical construction and improvement heuristics are described in Section 4. A tabu search heuristic is presented in Section 5, followed by computational results in Section 6 and by conclusions in Section 7. 2. Mathematical model To model the SVRPPD we associate two vertices i and i + n with each customer i, where pi+n = pi. The model allows two service possibilities for each customer i. The pickup and delivery operations may be performed simultaneously, in which case vertex i is visited and vertex i + n is not visited. Otherwise, customer i is visited twice: delivery is made at vertex i and pickup at vertex i + n. We define an extended cost matrix C ¼ ðcij Þð2nþ1Þð2nþ1Þ where 8 cij > > > < cin;j cij ¼ > ci;jn > > : cin;jn

if i 6 n and j 6 n; if i > n and j 6 n; if i 6 n and j > n; if i > n and j > n:

570

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

We also define the following variables: 8 > < 1 if the vehicle travels directly from i to j xij ¼ ði; j ¼ 0; . . . ; 2n; i 6¼ j; j 6¼ i þ n if 1 6 i 6 n; j 6¼ i  n if i > nÞ; > : 0 otherwise;  1 if pickup and delivery are performed simultaneously at customer i ði ¼ 1; . . . ; nÞ; yi ¼ 0 otherwise; wi ¼ an upper bound on the vehicle load upon leaving vertex i ði ¼ 0; . . . ; 2nÞ: The SVRPPD model is to minimize

2n X 2n X i¼0

ð1Þ

cij xij

j¼0

subject to 2n X

xij ¼ 1

ði ¼ 0; . . . ; nÞ;

ð2Þ

xij ¼ 1

ðj ¼ 0; . . . ; nÞ;

ð3Þ

j¼0 2n X i¼0 2n X

xij ¼ 1  y in

ði ¼ n þ 1; . . . ; 2nÞ;

ð4Þ

xij ¼ 1  y jn

ðj ¼ n þ 1; . . . ; 2nÞ;

ð5Þ

j¼0 2n X i¼0

w0 ¼

n X

ð6Þ

d i;

i¼1

0 6 wi 6 Q ði ¼ 1; . . . ; 2nÞ; wj P wi  d j þ pj y j  ð1  xij ÞQ

ði ¼ 0; . . . ; 2n; j ¼ 1; . . . ; nÞ;

wj P wi þ pj ð1  y jn Þ  ð1  xij ÞQ ði ¼ 0; . . . ; 2n; j ¼ n þ 1; . . . ; 2nÞ; X xij 6 jS j  1 ðS  f0; . . . ; 2ng; jS j P 2Þ;

ð7Þ ð8Þ ð9Þ ð10Þ

i;j2S

xij 2 f0; 1g ði; j ¼ 0; . . . ; 2n; i 6¼ j; j 6¼ i þ n if 1 6 i 6 n; j 6¼ i  n if i > nÞ; y i 2 f0; 1g ði ¼ 1; . . . ; nÞ:

ð11Þ ð12Þ

In this formulation, constraints (2) and (3) mean that the first vertex associated with each customer is visited once, either for a single delivery or for a simultaneous pickup and delivery. Constraints (4) and (5) express the fact that the second vertex associated with a customer is visited only if a combined pickup and delivery does not occur at the first vertex. Constraint (6) defines the vehicle load upon leaving the depot while constraints (7) guarantee that the vehicle load will never exceed the vehicle capacity. Constraints (8) and (9) define wj in terms of wi whenever j is visited immediately after i. In constraints (8), j is the first vertex associated with customer j: wj P wi  dj + pj if pickup and delivery are performed simultaneously at j, and wj P wi  dj otherwise. In constraints (9), j is the second vertex associated with customer j: wj P wi if pickup and delivery are performed simultaneously at j, and wj P wi + pj otherwise. Note that the first of the last two inequalities is redundant because it is not suboptimal to set xij = 0 whenever yjn = 1 for j P n + 1. Constraints (10) are the standard subtour elimination constraints (Dantzig et al., 1954). Contrary to what happens in several routing problems (Desrochers and Laporte, 1991), constraints (8) and (9) are not sufficient to eliminate subtours. This is because vehicle load does not always along the route. For example, subtours could P vary monotonically P be created over subsets S of V for which i2S d i ¼ i2S pi . Constraints (11) and (12) define binary conditions on the variables.

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

571

It is interesting to observe that in the case of non-simultaneous pickup and delivery at a customer, the model does not specify an order for these two operations. In practice, however, it is never suboptimal to perform the delivery first. The value of the model just described is mostly descriptive as its linear relaxation yields a rather poor lower bound and only small instances can be solved to optimality. However, the model can be used to analyze properties and shapes of optimal solutions. Constraints (8) and (9) can be strengthened by applying a lifting process. Proposition 1. The constraints wj P wi  d j þ pj y j  ð1  xij ÞQ þ ðQ þ d i þ d j  pi  pj Þxji wj P wi  d j þ pj y j  ð1  xij ÞQ þ ðQ þ d j  pi  pj Þxji

ði ¼ 0; . . . ; n; j ¼ 1; . . . ; nÞ;

ði ¼ n þ 1; . . . ; 2n; j ¼ 1; . . . ; nÞ

are valid for the SVRPPD. Proof. Constraints (8) can be rewritten as wj P wi  d j þ pj y j  ð1  xij ÞQ þ aji xji ;

ð13Þ

where currently aji = 0. We seek the largest value of aji so that (8) remains valid. If at the optimum xji = 0, then aji can take any value. If xji = 1, then xij = 0 and (13) becomes wj P wi  d j þ pj y j  Q þ aji :

ð14Þ

Two cases must be considered. (1) If i 6 n, because xji = 1, it follows from (8) that: wi P wj  d i þ pi y i :

ð15Þ

Combining (14) and (15) yields aji 6 Q + di + dj  piyi  pjyj. The right-hand side is minimized for yi = yj = 1. Therefore setting aji = Q + di + dj  pi  pj is valid. (2) If i P n + 1, it follows from (9) that: wi P wj þ pi ð1  y in Þ:

ð16Þ

Combining (14) and (16) yields aji 6 Q + dj  pi  pjyj + piyin. The right-hand side is minimized for yj = 1 and yin = 0. Therefore setting aji = Q + dj  pi  pj is valid. h Proposition 2. The constraints wj P wi þ pj ð1  y jn Þ  ð1  xij ÞQ þ ðQ þ d i  pi  pj Þxji wj P wi þ pj ð1  y jn Þ  ð1  xij ÞQ þ ðQ  pi  pj Þxji

ði ¼ 0; . . . ; n; j ¼ n þ 1; . . . ; 2nÞ;

ði ¼ n þ 1; . . . ; 2n; j ¼ n þ 1; . . . ; 2nÞ

are valid for the SVRPPD. Proof. The proof is similar to that of Proposition 1. h 3. Solution properties The model just described allows for the generation of feasible solutions in which any customer may be visited once or twice. Some special Pn cases are of particular interest. For example, forcing all customers to be visited only once (by P setting i¼1 y i ¼ n) yields a Hamiltonian solution. Similarly, a double-path solution is n obtained by setting i¼1 y i ¼ 1. This means that all customers will be visited by means of two Hamiltonian paths: one from the depot to some customer in, and a second one from in to the depot. In this type of solution only customer in has a simultaneous pickup and delivery, while all other customers are visited twice. These two types of solutions are in fact extreme cases of a lasso solution, a term first used by Halskau and Løkketangen (1998), and Gribkovskaia et al. (2001). In a lasso solution, the vehicle first performs deliveries to a subset S of customers, then visits all remaining customers to perform a combined delivery and pickup, and finally performs collections at the vertices of S (see Fig. 1). Setting jSj = 0 yields a Hamiltonian solution while setting jSj = jVj  2 yields a double-path solution. A solution having no predetermined shape and no constraint

572

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

S 0

Fig. 1. A lasso solution.

L

H

D G

Fig. 2. Feasible SVRPPD solutions.

1

1

1

0

1

2 (1,1)

(1,1)

Q = 26

1

3

5 (6,3) 1/2

3 /2

(6,3)

4 (6,15) 1/2

1

6 (6,3)

Fig. 3. Euclidean instance for which a non-lasso solution is optimal. The pickup and delivery demands are indicated by (pi, di).

on the number of customers visited twice is called general. Denote by G, L, H and D the sets of feasible general, lasso, Hamiltonian and double-path solutions, respectively. Fig. 2 depicts lasso, Hamiltonian and doublepath solution types as subsets of all general solutions to our problem. An optimal SVRPPD solution may belong to GnL, such as the solution to the Euclidean instance depicted in Fig. 3.pIn ffiffiffi this example the best lasso solution, which is also Hamiltonian, is (0, 1, 2, 4, 5, 6, 3, 0) and has a cost of 8:5 þ 3=2 ¼ 9:36; the best non-Hamiltonian lasso solution is (0, 1, 2, 3, 5, 4, 6, 3, 2, 1, 0) and has a cost of 11; the best double-path solution is (0, 1, 2, 3, 5, 4, 6, 4, 5, 3, 2, 1, 0) and has a cost of 12; the optimal solution, which is non-lasso, is (0, 1, 2, 3, 5, 4, 6, 3, 0) and has a cost of 9. To compare the solution shapes just described, denote by zP the optimal SVRPPD solution value for shape P, where P can be D (double-path), H (Hamiltonian), L (lasso), or G (general). Then obviously zG 6 zL 6 zD and zG 6 zL 6 zH . Proposition 3. If C satisfies the triangle inequality, then zG 6 zL 6 zH 6 zD . Proof. Let (0, i1, . . . , in1, in, in+1, . . . , i2n1, 0) be an optimal double-path solution of cost zD . This solution consists of visiting all vertices of V by first making deliveries on a Hamiltonian path from 0 to in, and then making pickups at the same vertices on a return Hamiltonian path from in to 0. If the cost matrix is symmetric these two paths will be the reverse of each other, i.e., in+k = ink for k = 1, . . . , n  1. A feasible Hamiltonian solution of cost zH with combined pickup and delivery at each customer is obtained by first following the subpath (0, i1, . . . , in1, in) but skipping all customers with di < pi, and then following the subpath (in, in+1, . . . , i2n1, 0) but skipping all customers with di P pi. The triangle inequality implies that zH 6 zD . The thesis follows since zH 6 zH . h Proposition 4. If C satisfies the triangle inequality, then zG 6 zL 6 zH 6 zD 6 2zG .

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

573

Proof. Consider an optimal general solution of cost zG given by the path (0, . . . , ia, . . . , ib, . . . , 0), where a is the largest index for which the subpath (ia, . . . , 0) contains all vertices of V, and b (Pa) is the smallest index for which the subpath (0, . . . , ib) contains all vertices of V. Now construct the path (0, . . . , ia, . . . , ib, . . . , 0, . . . , ia, . . . , ib, . . . , 0) of cost 2zG obtained by following the general solution twice. Then (1) follow the subpath (0, . . . , ia, . . . , ib) by skipping any vertex already visited, so that upon reaching ib every vertex has been visited exactly once; (2) proceed immediately to ia in the second part of the path by skipping all vertices between ib and ia; (3) follow the subpath (ia, . . . , ib, . . . , 0) by skipping ib and any vertex already visited twice. This yields a double-path solution (0, . . . , ia, . . . , ib, ia, . . . , 0) of cost zD 6 2zG because of the triangle inequality. The thesis follows since zD 6 zD . h 4. Classical heuristics We first describe a new construction heuristic with several variants. A simple improvement heuristic is then described. 4.1. Construction heuristic The construction heuristic consists of three phases: (1) initial Hamiltonian cycle; (2) initial general solution; (3) customer merging. It always yields a feasible solution. 4.1.1. Initial Hamiltonian cycle We construct an initial Hamiltonian cycle by means of one of the following three procedures. The last two only apply to planar instances. N: Nearest neighbour procedure (Rosencrantz et al., 1977). Starting with the depot, iteratively insert the customer closest to the last customer inserted. S: Sweep procedure (Gillett and Miller, 1974). Starting with the depot, construct a Hamiltonian cycle by iteratively including the customers in increasing order of the polar angle they make with an arbitrary radius rooted at the depot. M: Modified sweep procedure. Proceed as in the sweep procedure but initialize the tour not with the depot, but with the first customer having the smallest angle with the depot and an arbitrary radius rooted at the depot. Then insert the depot in the cycle in order to minimize the insertion cost. Note that any of these three initial Hamiltonian cycles may be infeasible for the SVRPPD. 4.1.2. Initial general solution Construct two Hamiltonian circuits (0, i1, . . . , in, 0) and (0, in, . . . , i1, 0) by following the initial Hamiltonian cycle (0, i1, . . . , in, 0) in both directions. Then generate 2n new solutions by successively removing n arcs from each circuit. 1. Consider the circuit (0, i1, . . . , in, 0). (a) Remove arc (0, i1) and construct the double-path solution (0, in, . . . , i1, i2, . . . , in, 0). (b) Remove arc (it, it+1) where 1 6 t 6 n  1, and construct the solution (0, i1, . . . , it, it1, . . .. , i1, in, in1, . . . , it+1, it+2, . . . , in,0). 2. Consider the circuit (0, in, . . . , i1, 0). (a) Remove arc (in, 0) and construct the double-path solution (0, i1, . . . , in, in1, . . . , i1, 0). (b) Remove arc (it+1, it) where 1 6 t 6 n  1, and construct the solution (0, in, in1, . . . , it+1, it+2, . . . , in, i1, . . . , it, it1, . . . , i1, 0). Examples of these solutions are depicted in Fig. 4. All solutions are general, according to the terminology introduced in Section 3. Note that only the two double-path solutions are guaranteed to be feasible. All infeasible solutions generated in this phase are discarded.

574

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

Initial Hamiltonian circuits 1

i1

2

i2

0

i4

i2

i1

i4

i1

i2

0

i3

b) Remove arc (it , it +1 ) 1≤t≤n-1

i4

a) Remove arc (0, in )

0

i5

i3

i5

a) Remove arc (0, i1 ) i1

i2

0

i3

i5

i1

i3

i5

i4

b) Remove arc (it +1 , it ) 1≤t≤n-1

,

i2

i1

,

i2

t=1 0

0

i5 i1

i4 i2

i3

i5

i4

i1

i2

t=2 0

i3

i5

i4

i1

i2

0

i3

i5

i4

i1

i2

t=n-1 0

i3

i5

i4

0

i3

i5

i4

Fig. 4. Construction of 2n initial general solutions.

4.1.3. Customer merging In the feasible solutions just generated some customers are visited once while the others are visited twice. We attempt to reduce the number of double visits by means of a forward or a backward merging procedure. As illustrated in Fig. 4, all initial general solutions consist either of a double-path, or of a double-path incident to a three-vertex circuit including the depot, or of two double-paths incident to such a circuit.

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

575

Initial general solution

F: Forward merging

B: Backward merging

Fig. 5. Forward and backward merging procedures.

F: Forward merging procedure. Starting with the first vertex of a double-path, combine the two visits made at that vertex into a single visit if this is feasible. Proceed to the next vertex of the double-path, moving away from the depot, and continue until all vertices of the double-path have been scanned. Then go to the first vertex of the other double-path if there are two double-paths, and continue until all vertices have been considered (see Fig. 5). B: Backward merging procedure. Proceed as in the forward case starting this time with the penultimate vertex of each double-path, and scanning vertices moving this time towards the depot (see Fig. 5). Each merging procedure is applied in turn to all feasible initial general solutions and returns the overall best final solution. 4.2. Improvement heuristic The improvement heuristic iteratively relocates a vertex i in a feasible route s as long as a cost reduction can be achieved. Convergence is controlled by keeping track of the last vertex k whose relocation has resulted in a cost reduction. The steps of this heuristic are: Step 1: Initialization Let s be a feasible route. Set i :¼ 1 and k :¼ 1. Step 2: Vertex removal If i is visited twice or if removing i from s would change the number of visits of another vertex (see, e.g., Fig. 6), go to Step 4. Otherwise, remove i from s and connect its predecessor and successor in s to obtain an infeasible route s.

576

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

i

j

j

Fig. 6. Removing i changes the number of visits to j from 2 to 1.

Step 3: Vertex insertion Insert i in s, possibly in its original position, using a cheapest insertion criterion, while preserving feasibility. If the resulting route s 0 yields a cost reduction with respect to s, set s :¼ s 0 and k :¼ i. Step 4: Next vertex Let j be the vertex following i in s. Stop if j = k. Otherwise, set i :¼ j and go to Step 2.

5. Tabu search heuristic We have also developed a tabu search (TS) heuristic for the SVRPPD. This metaheuristic was first proposed by Glover (1986) and has since been successfully applied to a variety of problems. Tabu search is essentially an improvement method. Starting with an initial solution, it explores the solution space by moving at each iteration from a solution s to the best solution s 0 in a subset of the neighbourhood N(s) of s. The process stops after a termination criterion has been satisfied. If c(s) denotes the cost of s, then c(s 0 ) may be larger than c(s) and it is therefore necessary to implement an anti-cycling mechanism. More specifically, if B(s) denotes the set of all solutions possessing a given attribute of solution s, then any solution s 0 having an attribute in B(s) is forbidden, or tabu, for h iterations, unless c(s 0 ) is less than the cost of the best known solution encountered in B(s), in which case it is clear that cycling cannot occur. This last rule is referred to as an aspiration criterion. A number of other features are normally necessary in order to obtain a successful TS algorithm. In our implementation, we use a continuous diversification mechanism in order to provide a broad exploration of the search space, we perform a local reoptimization of each new best known solution encountered, and we allow intermediate infeasible solutions during the search. In what follows we describe the main features of our TS heuristic, which is an adaptation of the Unified Tabu Search Algorithm proposed by Cordeau et al. (2001). 5.1. Main features of the tabu search heuristic 5.1.1. Initial solution We have considered three different construction heuristics to generate an initial solution. HC: If the customers are represented by Euclidean coordinates they are first sorted according with the angle they make with an arbitrary radius rooted at the depot. Otherwise the customers are considered in an arbitrary order. A Hamiltonian circuit is constructed by means of a cheapest insertion procedure. This initial solution may be infeasible. PD2T: Construct a feasible Hamiltonian circuit by means of Mosheiov’s PD2T algorithm (1994) which can be summarized as follows: Step 1: Using a cheapest insertion heuristic, determine a Hamiltonian circuit (i1, . . . , in, i1) on a graph G 0 from which vertex 0 has been removed, disregarding pickup and delivery demands.

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

577

Pr Step 2: Let ir be such that t¼1 ðpit  d it Þ is maximized. Step 3: Then the Hamiltonian circuit (0, ir+1, . . . , in, i1, . . . , ir, 0) is feasible. Con*: Consider the best feasible solution generated by means of the construction heuristics described in Section 4.1.

5.1.2. Penalized objective function In order to allow a mix of feasible and infeasible solutions, the algorithm works with a penalized function f(s) = c(s) + aq(s), where c(s) is the cost of solution s, q(s) is the total vehicle capacity violation of s, and a is a self-adjusting parameter. At each iteration, a is divided by 1 + d > 1 if the current solution is feasible, and multiplied by 1 + d otherwise, where d is a user-controlled parameter. This mechanism is identical to that of Cordeau et al. (1997). 5.1.3. Neighbourhood structure and attribute set At each iteration, the neighbourhood N(s) of a solution s is defined as the set of all solutions reachable from s by changing the number v of visits of a customer. With s is associated an attribute set B(s) = {(i, v): i 2 Vn{0} and v = 1 or 2}. A transition from s to a neighbour s 0 is called a move, which can be defined as the removal of an attribute (i, v) from B(s) and the inclusion of (i, v 0 ) in B(s 0 ), where v 0 5 v. There are two possible moves. 1. If v = 1, a second visit of i is inserted in the current solution so as to yield the smallest increase of f(s). Visiting i twice will typically increase c(s) and decrease q(s). 2. If v = 2, the second occurrence of vertex i in the solution is deleted and its predecessor and successor are linked together. As a result pickup and delivery are now performed simultaneously at vertex i.

5.1.4. Tabu status of an attribute If (i, v) is replaced with (i, v 0 ), then the reinclusion of (i, v) in B(s) is forbidden for h iterations, where h is a user-controlled parameter. 5.1.5. Aspiration criterion The aspiration level riv of attribute (i, v) is initially defined as the cost of the initial solution s if (i, v) 2 B(s) and s is feasible; otherwise, riv = 1. Every time a feasible solution is identified, the aspiration level of attribute (i, v) is updated to min{riv, c(s)}. When considering a solution s 0 obtained by the inclusion of a tabu attribute (i, v) in B(s 0 ), the tabu status of (i, v) is revoked if q(s 0 ) = 0 and c(s 0 ) < riv. Let M(s)  N(s) be the subset of neighbour solutions satisfying the aspiration criterion, i.e., solutions for which (i, v 0 ) 2 B(s 0 )nB(s) and s 0 is non-tabu, or s 0 is feasible and c(s 0 ) < riv. In other words, the search only moves to a tabu solution if this yields a new best solution among those possessing attribute (i, v). 5.1.6. Diversification In order to diversity the search, any solution s 0 2 M(s) such that f(s 0 ) P f(s) is penalized by a term p(s 0 ) = fqivc(s), where f is a positive parameter, and qiv is the relative frequency of all iterations for which attribute (i, v) has been in B(s 0 ); if f(s 0 ) < f(s), then p(s 0 ) = 0. The best s 0 2 M(s) is the solution for which g(s 0 ) = f(s 0 ) + p(s 0 ) is minimized. This type of diversification scheme was introduced by Glover (1989) and later fine tuned by Taillard (1993). 5.1.7. Local reoptimization Route reoptimization is performed every u iterations, where u is a user-controlled parameter, or whenever a new best feasible solution is encountered. This is done by means of the improvement heuristic of 4.2 in which the penalized function f(s) is used instead of c(s). 5.1.8. Termination criterion The algorithm is applied for a fixed number g of iterations, where g is a user-controlled parameter.

578

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

5.2. Step-by-step description of the algorithm We now provide a step-by-step description of the TS algorithm. It is applied for g iterations. The other user-controlled parameters, already defined, are d, h, f and u. If the algorithm is applied from a feasible initial solution, it returns the best feasible solution s* identified. Note that the PD2T and Con* construction procedures always generate a feasible initial solution. We use the following notation: (i, v) B(s) c(s) q(s) N(s) M(s) s; es s0 s* a d f g h k qiv riv siv u

attribute: number of visits v at customer i; attribute set of solution s; routing cost of solution s; total load violation of the vehicle during the route for given solution s; neighbourhood of solution s; a subset of N(s); solutions; initial solution; best solution identified; penalty factor for overload; parameter used to update a; factor used to adjust the intensity of the diversification; total number of iteration to be performed; tabu duration; iteration counter; number of times attribute (i, v) has been added to the solution; aspiration level of attribute (i, v); last iteration for which attribute (i, v) is declared tabu; parameter used for route reoptimization.

The steps of the algorithm are: Step 1: Initial solution a. Construct an initial solution s0. b. Set s :¼ s0 and a :¼ 1. If s is feasible, set s* :¼ s. Step 2: Attribute set and functions initialization For every (i, v), do a. Set qiv :¼ 0, siv :¼ 0. b. If (i, v) 2 B(s) and s is feasible, set riv :¼ c(s); else, set riv :¼ 1. Step 3: Iterative search For k = 1, . . . , g, do a. Update d, f, h: • d :¼ random value from the interval (0, 1). • f :¼ random value from the interval (0, 1). • h :¼ random integer from the interval [1, [10 log10 n]]. b. Set N(s) :¼ ;. c. For each attribute (i, v 0 ) 62 B(s) do • Create a solution s 0 applying a move definition. • Replace the corresponding attribute (i, v) 2 B(s) by the attribute (i, v 0 ), i.e., move from (i, v) 2 B(s) to (i, v 0 ) 2 B(s 0 ). • Set N(s) :¼ N(s) [ {s 0 }. d. Set M(s) :¼ ;. e. For each s 0 2 N(s) do • For (i, v) 2 B(s 0 )nB(s) such that siv < k or such that c(s) + f(s 0 )  f(s) < riv, set M(s) :¼ M(s) [ {s 0 }.

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

579

f. For each s 0 2 M(s), do • If f(s 0 ) P f(s), set g(s 0 ) :¼ f(s 0 ) + fc(s)qiv/k; else, set g(s 0 ) :¼ f(s 0 ). g. Identify a solution s 0 2 M(s) minimizing g(s 0 ). h. For attribute (i, v) 2 B(s 0 )nB(s) do • Set qik :¼ qik + 1 and siv :¼ k + h. • If s 0 is feasible (q(s 0 ) = 0), do – If c(s 0 ) < c(s*), set s* :¼ s 0 . – For each (i, v) 2 B(s 0 ), set riv :¼ min{riv, c(s 0 )}. – Set a :¼ a/(1 + d). • If s 0 is infeasible, set a :¼ a(1 + d). j. Set s :¼ s 0 . k. If q(s) = 0 and c(s) = c(s*), or k = ku, k = 1, 2, . . . , do route reoptimization: • For each (i, 1) 2 B(s) do – Remove customer i from the route in solution s. – Reinsert customer i in the route in solution es minimizing f ðes Þ ¼ cðes Þ þ aqðes Þ such that BðsÞ ¼ Bðes Þ. – Set s :¼ es . • If both q(s) = 0 and c(s) < c(s*), set s* :¼ s. l. Set k :¼ k + 1. 6. Computational results All heuristics were coded in C using Pelles for Windows, version 2.90.8. All tests were run on an AMD Athlon XP 1.8 GHz (Intel 1.5 GHz), with a 256 MB memory, and using the operating system Microsoft Windows XP Home Edition. 6.1. Test instances We have created 34 test instances, 17 of type A and 17 of type B. These instances are derived from the capacitated vehicle routing problem (CVRP) instances of VRPLIB: http://www.or.deis.unibo.it/research_ pages/ORinstances/VRPLIB/VRPLIB.html. We have selected the following 17 CVRP instances containing between 16 and 101 vertices: E-016-03, E-021-04, E-022-04, E-023-03, E-026-08, E-030-03, E-031-09, E-03304, E-036-11, E-041-14, E-045-04, E-048-04, E-051-05, E-072-04, E-076-07, E-101-08, E-101-10. The letter E means that the instances are generated in the plane with Euclidean distances, the next three digits represent the number of vertices (including the depot), while the last two digits correspond to the number of vehicles. This number is irrelevant in the SVRPDP. In the A instances delivery and pickup demands were generated as follows. Let qi be the demand of customer i in the original CVRP instance. We then set di = qi and pi = [(1  b)qi] if i is odd or pi = [(1 + b)qi] if i is even, where 0 6 b < 1. We selected b = 0.20. The B instances have the same coordinates as the CVRP instances, and di = qi for i = 1, . . . , n. To generate the pickup demands, the customers were first ranked in increasing distance from the depot. For the first third of customers we set pi = 1.2qi. For the second third pi P was generated as for the P  n n . A instances. For the last third pi = 0.8qi. In all A and B instances, Q ¼ max p ; d i i¼1 i i¼1 Because optimal solutions and even good lower bounds cannot be generated for these instances, we can only compare heuristics between themselves and report deviations from the best known solution values. 6.2. Overview of heuristics used in the computational comparisons We have first compared the six construction heuristics described in Section 4.1: NB: nearest neighbour procedure with backward merges; NF: nearest neighbour procedure with forward merges;

580

SB: SF: MB: MF:

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

sweep procedure with backward merges; sweep procedure with forward merges; modified sweep procedure with backward merges; modified sweep procedure with forward merges.

We have also tested the following combinations: Con* + improvement: solve the SVRPPD with Con* (see Section 5.1.1) and apply the improvement heuristic of Section 4.2; PD2T + improvement: solve the SVRPPD with Mosheiov’s PD2T heuristic and apply to it the improvement heuristic of Section 4.2; HC + TS: apply tabu search starting with a Hamiltonian circuit obtained by means of a cheapest insertion procedure; PD2T + TS: apply tabu search starting with a Hamiltonian circuit generated by Mosheiov’s PD2T heuristic; Con* + TS: apply tabu search starting with the best solution generated by the six construction heuristics. In Tables 1–4 we report average gap values, where the gap is defined as the percentage deviation with respect to the best known solution identified in all our experiments. The gap is equal to 100[c(s*)  c(s**)]/ c(s**), where s* is the best solution found by the algorithm for a given instance, and s** is the overall best known solution for the same instance. Computational times are only reported when they are significant. 6.3. Assessment of the classical algorithms We first compare our six construction heuristics in Table 1. Computational times are omitted because they are on average less than one second. The overall best construction algorithm is MB. The results show that Table 1 Average gap values for six construction heuristics Instances

NB

NF

SB

SF

MB

MF

A B Average

54.83 55.06 54.95

207.60 205.88 206.74

54.55 55.00 54.78

207.01 206.24 206.63

50.89 48.84 49.87

51.91 49.53 50.72

Table 2 Average gap values for the Con* and PD2T heuristics Instances

A B Average

Without improvement

With improvement

Con*

PD2T

Con*

PD2T

47.14 45.81 46.48

24.11 23.41 23.76

5.67 4.61 5.14

9.88 7.16 8.52

Table 3 Average gap values and computational time over A and B instances obtained by the Con* + TS with different u values and g = 105 u

Gap

Seconds

2 3 5 10 20 50

1.27 1.28 1.26 1.22 1.23 1.50

1060 826 638 495 420 378

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

581

Table 4 Average gap values obtained by the TS algorithm with three initialization procedures, u = 10 and g = 104, 105 HC + TS

PD2T + TS

Con* + TS

g = 10 A B Average

2.52 2.50 2.51

3.95 3.78 3.87

2.24 1.35 1.80

g = 105 A B Average

2.04 1.98 2.01

3.33 2.91 3.12

1.44 0.99 1.22

4

irrespective of the first part of the construction heuristic, the backward merging procedure is better than the forward procedure because it tends to merge customers located far away from the depot, thus yielding larger cost savings. In the same spirit, the MF combination performs very well compared with NF and SF because after having inserted the depot in its best possible position in the route, the unmerged customers tend to be close to the depot. Similarly MB outperforms NB and SB. Overall there is no significant difference between A and B instances. In Table 2, we compare Con* and PD2T with or without the improvement phase. Again computation times are insignificant. When no improvement phase is applied, Con* is significantly outperformed by PD2T. However, the result is the opposite when an improvement phase is applied. In fact, Con* + improvement is quite good, with an average gap of only 5.14% with respect to the best known solutions. In all cases slightly better results are obtained on the B instances. 6.4. Assessment of the tabu search algorithm As indicated in Section 5.2 the TS algorithm works with five user-controlled parameters: d h f u g

1 + d is the multiplier or the divisor of a in the objective function f(s); tabu tenure; multiplier used in the penalized objective function p(s); frequency of route reoptimization; number of iterations.

The first three of these parameters have already been the subject of intensive testing. Cordeau et al. (2001) report that their best values are relatively stable, irrespective of the routing application under consideration. We therefore use their recommended values: d and f are randomly selected in the interval (0, 1), and h is an integer randomly selected in [1, [10 log10 n]]. The parameter u is particular to the SVRPPD and has been the subject of some investigation. Finally, g really depends on how long the user is willing to run the algorithm. Typically solution quality is correlated with g. We have performed some experiments on u. Recall that this parameter controls the frequency of route reoptimization by means of the improvement heuristic of Section 4.2. Table 3 indicates that u = 10 yields the smallest average gap observed for all values of this parameter used in the comparison. As expected, computational time is inversely related to the size of u. Finally, we present in Table 4 average gap values obtained by applying the TS algorithm starting with three different construction heuristics, using u = 10, and g = 104 and 105. These three versions of TS algorithm yield smaller average gaps on B instances than on A instances. It can be seen that applying TS starting with Con* yields a much better final solution than if HC or PD2T are used to generate the initial solution. Using g = 105 instead of 104 reduces the average gap by about 32% when the Con* + TS heuristic is used. Note that running Con* + TS for 104 iterations produces a smaller gap than running HC + TS or PD2T + TS for 105 iterations. Computational times are similar for the three algorithms. As expected, computational times are about 10 times longer with g = 105 than with g = 104 (about 500 seconds compared with 50 seconds).

582

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

6.5. Solution shapes Finally, we depict in Fig. 7 the shapes of all solutions encountered during our experiments. The first digit in the codes represents the number of vertices visited twice. This number is always at most 3. The shapes of all best solutions produced by our algorithms are of type 0, 1a, 1b, 1c, 2c or 2d. Overall, the percentage of Hamiltonian solutions is 62% over all instances, 32% of the instances contain one customer visited twice, and 6% contain two customers visited twice. It is also interesting to compare the frequency of the shapes of the best solution identified for A and B instances. These results are summarized in Table 5. Hamiltonian solutions (shape 0) are less frequently encountered in B instances than in A instances. This can be explained by the fact that in B instances, customers located close to the depot have a large pickup demand compared to their deliv-

0

1a 1b

1c

2a

2h

2b 2i 2c 2j

2d

2e 2k 2f

2l

2g

3a

3b

3f

3c

3d

3e Fig. 7. Solution shapes generated during the experiments, where arcs are shown by full lines and multi-arc paths by dotted lines.

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

583

Table 5 Shapes of the overall best known solutions Instance

A

B

E-016-03 E-021-04 E-022-04 E-023-03 E-026-08 E-030-03 E-031-09 E-033-04 E-036-11 E-041-14 E-045-04 E-048-04 E-051-05 E-072-04 E-076-07 E-101-08 E-101-10

0 0 0 0 0 0 0 0 0 0 0/1a 0 0 0 1a 1b 1b

0 0 0 0 1b 0 1a 0 1a 1a 1a/2c 1a 1a 0 0 1c 2d

ery demand. It is therefore often advantageous to first make a delivery at the customers and perform the pickup in a second visit, on the way back to the depot. Overall our results clearly show that it is beneficial to construct models and algorithms capable of producing general solutions for the SVRPPD, instead of imposing a predefined solution shape such as Hamiltonian or lasso. 7. Conclusion We have proposed a mixed integer programming model and heuristics for the single vehicle routing problem with pickups and deliveries in which each customer may be visited once or twice, giving rise to general solutions that encompass known solution shapes such as Hamiltonian, double-path and lasso. We have also designed classical construction and improvement heuristics, as well as a tabu search heuristic for the problem. Our results show that the best known solutions generated by these heuristics are frequently non-Hamiltonian and may contain up to two customers visited twice. This indicates that it is advantageous to develop models and algorithms capable of producing such general solutions. Acknowledgements This work was partially supported by the Norwegian Research Council and by the Canadian Natural Sciences and Engineering Council under grant OGP0039682. This support is gratefully acknowledged. Thanks are due to the referees for their valuable comments. References Angelelli, E., Mansini, R., 2002. The vehicle routing problem with time windows and simultaneous pick-up and delivery. In: Klose, A., Speranza, M.G., Van Wassenhove, L.N. (Eds.), Quantitative Approaches to Distribution Logistics and Supply Chain Management. Springer, Berlin, pp. 249–267. Cordeau, J.-F., Gendreau, M., Laporte, G., 1997. A tabu search algorithm for periodic and multi-depot vehicle routing problems. Networks 30, 105–119. Cordeau, J.-F., Laporte, G., Mercier, A., 2001. A unified tabu search algorithm for vehicle routing problems with time windows. Journal of the Operational Research Society 52, 928–936. Cordeau, J.-F., Laporte, G., Potvin, J.-Y., Savelsbergh, M.W.P., in press. Transportation on demand. In: Barnhart, C., Laporte, G. (Eds.), Transportation, Elsevier, Amsterdam. Dantzig, G.B., Fulkerson, D.R., Johnson, S.M., 1954. Solutions of a large-scale traveling-salesman problem. Operations Research 2, 393– 410.

584

I. Gribkovskaia et al. / European Journal of Operational Research 180 (2007) 568–584

Desrochers, M., Laporte, G., 1991. Improvements and extensions to the Miller–Tucker–Zemlin subtour elimination constraints. Operations Research Letters 10, 27–36. Dethloff, J., 2001. Vehicle routing and reverse logistics: The vehicle routing problem with simultaneous delivery and pick-up. Operations Research Spectrum 23, 79–96. Dethloff, J., 2002. Relation between vehicle routing problems: An insertion heuristic for the vehicle routing problem with simultaneous delivery and pick-up applied to the vehicle routing problem with backhauls. Journal of the Operational Research Society 53, 115–118. Gendreau, M., Laporte, G., Vigo, D., 1999. Heuristics for the traveling salesman problem with pickup and delivery. Computers & Operations Research 23, 501–508. Gillett, B.E., Miller, L.R., 1974. A heuristic algorithm for the vehicle dispatch problem. Operations Research 22, 340–349. Glover, F., 1986. Future paths for integer programming and links to artificial intelligence. Computers & Operations Research 13, 533–549. Glover, F., 1989. Tabu search – Part I. ORSA Journal on Computing 1, 190–206. Gribkovskaia, I., Halskau, Ø., Myklebost, K.N.B., 2001. Models for pick-up and deliveries from depots with lasso solutions. In: Stefansson, G., Tilanus, B. (Eds.), Proceedings of the 13th Annual Conference on Logistics Research NOFOMA 2001. Chalmers University of Technology, Go¨teborg, pp. 279–293. Halse, K., 1992. Modeling and solving complex vehicle routing problems, Ph.D. dissertation, IMSOR, Technical University of Denmark. Halskau, Ø., Gribkovskaia I., 2002. Pick-up and deliveries from depots with different sub-graphs in the solution, presented at the IFORS Conference, Edinburgh, 2002. Halskau, Ø., Løkketangen, A., 1998. Analyse av distribusjonsopplegget ved Sylte Mineralvannfabrikk A/S, Report M9812 Møreforskning, Molde. Hoff, A., Løkketangen, A., in press. Creating lasso-solutions for the TSPPD-problem by tabu search, Central European Journal of Operations Research. Min, H., 1989. The multiple vehicle routing problem with simultaneous delivery and pick-up points. Transportation Research A 23, 377– 386. Mosheiov, G., 1994. The travelling salesman problem with pick-up and delivery. European Journal of Operational Research 79, 299–310. Nagy, G., Salhi, S., 1998. The many-to-many location-routing problem. TOP 6, 261–275. Nagy, G., Salhi, S., 2005. Heuristic algorithms for single and multiple depot vehicle routing problems with pickups and deliveries. European Journal of Operational Research 162, 126–141. Prive´, J., Renaud, J., Boctor, F.F., Laporte, G., in press. Solving a vehicle routing problem arising in soft drink distribution, Journal of the Operational Research Society. Rosencrantz, D.J., Stearns, R.E., Lewis II, P.M., 1977. An analysis of several heuristics for the traveling salesman problem. SIAM Journal on Computing 6, 563–581. Salhi, S., Nagy, G., 1999. A cluster insertion heuristic for single and multiple depot vehicle routing problems with backhauling. Journal of the Operational Research Society 50, 1034–1042. Taillard, E´.D., 1993. Parallel iterative search for vehicle routing problems. Networks 23, 661–673. Tang, F.A., Galva˜o, R.D., 2002. Vehicle routing problems with simultaneous pick-up and delivery service. Opsearch 39, 19–33. Tang, F.A., Galva˜o, R.D., 2006. A tabu search algorithm for the vehicle routing problem with simultaneous pick-up and delivery service. Computers & Operations Research 33, 595–619. Wasner, M., Za¨phel, G., 2004. An integrated multi-depot hub-location vehicle routing model for network planning of parcel service. International Journal of Production Economics 90, 403–419.