The periodic vehicle routing problem with intermediate facilities

The periodic vehicle routing problem with intermediate facilities

European Journal of Operational Research 137 (2002) 233±247 www.elsevier.com/locate/dsw The periodic vehicle routing problem with intermediate facil...

153KB Sizes 0 Downloads 55 Views

European Journal of Operational Research 137 (2002) 233±247

www.elsevier.com/locate/dsw

The periodic vehicle routing problem with intermediate facilities Enrico Angelelli *, Maria Grazia Speranza Department of Quantitative Methods, University of Brescia, C.da S.Chiara 48/b, 25122 Brescia, Italy

Abstract In this paper, we study an extension of the PVRP where the vehicles can renew their capacity at some intermediate facilities. Each vehicle returns to the depot only when its work shift is over. For this problem we propose a tabu search (TS) algorithm and present computational results on a set of randomly generated instances and on a set of PVRP instances taken from the literature. Ó 2002 Elsevier Science B.V. All rights reserved. Keywords: Periodic vehicle routing; Intermediate facilities; Renewable capacity; Heuristics; Tabu search

1. Introduction In the periodic vehicle routing problem (PVRP) a set of customers have to be visited, one or several times, on a given time horizon. The set of dates in which a vehicle serves a customer is not ®xed a priori, but instead, a list of possible sets of dates (visiting schedules) is associated with each customer. A ¯eet of vehicles is available and each vehicle leaves the depot, serves a set of customers and, when its work shift or its capacity is over, returns to the depot. The objective of the problem is the minimization of the total length of the routes travelled by the vehicles on the time horizon. Solving the problem requires assigning a visiting schedule to each customer and, for each day of the time horizon, de®ning the routes of the vehicles in such a way that all the customers, whose assigned

*

Corresponding author.

schedule includes that day, are served. This model is a periodic version of the classical vehicle routing problem and the periodic characteristic of the PVRP is essential in various applications. For instance, in waste collection problems, each customer has to be typically served with a given periodicity, say once a week, and the working days include all the days of a week, except Sunday. In this case the list of visiting schedules is such that each schedule includes one of the working days of the week. Unfortunately, the routing problems in waste collection applications cannot typically be modelled with the PVRP, due to the presence of the waste treatment plants. The vehicles return to the depot only when the work shift is over. When the capacity limit is reached the vehicles renew their capacity by unloading the waste at one of the treatment plants. In the case of distribution problems, an intermediate facility may be a warehouse where the goods are stored and where a vehicle loads the goods any time it starts a new route.

0377-2217/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 2 2 1 7 ( 0 1 ) 0 0 2 0 6 - 5

234

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

In this paper, we present a generalization of the PVRP in which some intermediate facilities exist where the vehicles renew completely their capacity. Each vehicle can renew its capacity at possibly di€erent intermediate facilities and returns to the depot only when the work shift is over. The objective of the problem is the minimization of the total length of the routes travelled by the vehicles. We propose a tabu search (TS) heuristic for the solution of the PVRP with intermediate facilities (PVRP-IF). We test di€erent versions of the algorithm on PVRP instances taken from the literature and on randomly generated instances. The paper is organized as follows. In the rest of this section, we discuss the related literature. In Section 2, we de®ne a model for the problem. In Section 3, we present the solution algorithm and in Section 4 we report the results obtained when the algorithm is run on PVRP instances taken from the literature and on new random instances speci®c for the problem. 1.1. Literature Early formulations of the PVRP were developed by Beltrami and Bodin [2] and by Russell and Igo [15] who proposed heuristics applied to waste collection problems. The problem we discuss in this paper is a generalization of the PVRP as stated, among others, in Christo®des and Beasley [5] and Cordeau et al. [6]. A number of papers on the PVRP exists in the literature. Christo®des and Beasley [5] present a heuristic that assigns a priority level to each customer and, according to this priority, assigns a visiting schedule to each customer so that the increase of an estimation of the costs of the vehicle routing problems (VRPs) related to the schedules is minimized. Then, an improving phase follows where an attempt is made to decrease the estimation of the VRPs by changing the visiting schedules of the customers. Finally, a VRP is solved for every day. Tan and Beasley [17] use the idea of the generalized assignment method proposed by Fisher and Jaikumar [8] and assign a visiting schedule to each vertex. Eventually a heuristic for the VRP is applied to each day. Russell and Gribbin [14] de-

velop a heuristic organized in four phases. Chao et al. [4] combine a number of techniques to assign visiting schedules, build and improve routes and reinitialize the optimization process. In particular they relax the capacity constraints in order to give more freedom to the minimization of routes during the construction phase. Cordeau et al. [6] present an algorithm that is, to our knowledge, the best heuristic for the PVRP so far developed. The solution algorithm is a TS heuristic which, di€erently from the above heuristics, may allow infeasible solutions during the search process. Infeasibility is allowed only with respect to capacity of vehicles and duration of routes. The algorithm adopts a ``sweep'' technique to calculate the initial solution, moves to a solution which minimizes a penalized cost function and iterates. The penalized cost depends on possible violation of capacity and time constraints and on the frequency of repeated moves. The algorithm implements two kinds of moves: move a customer from a route to another route in the same day, and change the schedule of a customer. A heuristic for a variant of the PVRP, where the objective is to minimize the number of vehicles, is proposed by Gaudioso and Paletta [9]. 2. Model formulation In this section, we brie¯y recall the basic elements of the vehicle routing and the PVRP. Then we present the PVRP with intermediate facilities. 2.1. The vehicle routing and the periodic vehicle routing problems The basic version of the VRP is de®ned on a graph G ˆ …V ; A†, where V ˆ fv0 ; v1 ; . . . ; vn g is the set of vertices and A  V  V is the set of arcs. The vertex v0 corresponds to the depot. A set of m vehicles is available where each vehicle has a capacity Qi , i ˆ 1; . . . ; m (vehicles may be identical). Every other vertex vi corresponds to a customer whose demand qi is known. A cost, or distance, cij is associated with each arc …vi ; vj † 2 A. If there is no arc connecting a pair …vi ; vj † the corresponding cost is set to in®nity. The objective is to de®ne m routes on

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

graph G, one per vehicle, so that the global cost, or length, of the routes is minimized. A feasible solution is such that every route starts and ends at the depot v0 , each customer vi …i 6ˆ 0† is assigned to one and only one route and the load of each vehicle does not exceed the vehicle capacity. Variants of the above formulation have been studied. If a time window is assigned to each customer, the time windows VRP is obtained. Another generalization of the VRP is the multi-depot VRP in which there is more than one depot and each vehicle starts and ends its route at the depot to which it is assigned. We refer to [3,7,11] for surveys on VRP and to [12,13] for bibliographies. All the variants of the classical VRP aim at planning the routes of a ¯eet of vehicles for a single time instant, say a day. The demand of the customers is de®ned for one day only. Situations where the day in which each customer has to be served is ®xed can be modelled with VRP as, for each day, a VRP has to be solved. A di€erent multi-day model is required when there is the possibility of choosing when to serve a customer. The PVRP is such a model. A T days horizon and a list of feasible visiting schedules for each customer are given (the customer's demand is speci®ed for each day of any possible visiting schedule). Each customer must be assigned one and only one feasible visiting schedule according to which it has to be visited and then, for each day of the planning horizon, a VRP has to be solved. Obviously, the choice of the visiting schedules and the de®nition of the routes are interrelated problems. The objective of the PVRP is the minimization of the global length of the routes. In the PVRP formulated, for instance, by Cordeau et al. [6], in addition to the classical characteristics of the VRP, a service time is given for each customer and a travelling time is given for each arc, together with a maximum duration of a route (work shift). 2.2. The PVRP with intermediate facilities The PVRP-IF is a generalization of the PVRP, where a set of intermediate facilities are introduced to serve as loading or unloading facilities.

235

In the case of a distribution problem, the intermediate facilities represent the warehouses where the goods are stored. A vehicle, which is empty when leaving the depot, goes directly to an intermediate facility for loading the goods, serves a set of customers, may return to a facility for the loading and returns to the depot only when the work shift is over. Let us note that no customer is allowed on the way from the depot to the ®rst facility because the vehicle is empty. In the case of a collection problem, the intermediate facilities represent the sites where the vehicles are unloaded. A vehicle leaves the depot and starts collecting the goods from the customers. When the vehicle is full, it reaches an intermediate facility for the unloading operation, then it may start for another collection tour and returns to the depot when the work shift is over. In case of a collection problem, di€erently from the distribution problem, customers are allowed on the path from the depot to the ®rst facility, but are not on the path from the last facility to the depot. In this paper, we describe a model for a collection problem. However, the model can be easily applied to a distribution problem as follows. First of all, the matrix of the distances and the matrix of the travelling times of the distribution problem have to be transposed. Then we can apply the model for the collection problem interpreting the demand of a customer as the quantity to be collected. Finally, when a solution for the collection problem is found, the orientation of all the arcs in the solution has to be reversed. The solution obtained in this way is the solution of the original distribution problem. Let us now formally de®ne the problem. The number of days in the planning horizon is denoted by T. A set of m identical vehicles is given and we denote by Q and D the vehicle capacity and the duration of the working shift, respectively. A graph G ˆ …V ; A† is de®ned as follows. A depot v0 , a set of n customers fv1 ; . . . ; vn g and a set of k facilities fvn‡1 ; . . . ; vn‡k g are identi®ed with the vertices of the graph, so that jV j ˆ n ‡ k ‡ 1. A service time ti is assigned to each vertex vi 2 V , which represents the time for the vehicle to leave the depot, to serve the customer or to use a facility. Each customer can be served at most once per day.

236

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

A set Hi ˆ fSi1 ; . . . ; Sihi g of visiting schedules is assigned to each customer, i.e. to each vertex vi (i ˆ 1; . . . ; n). A visiting schedule S is de®ned by a vector of T non-negative components S ˆ …q1 ; . . . ; qT †, where qd denotes the demand on day d (qd ˆ 0 means no visit on day d). We denote by cij the distance from vertex vi to vertex vj , or the cost to travel from vi to vj . Besides, a travelling time tij from vertex vi to vertex vj is associated with each arc …vi ; vj † 2 A. The problem consists in assigning to each customer a visiting schedule and, for each day of the planning horizon, ®nding for each vehicle a route in such a way that the total length of the routes travelled by the vehicles on the horizon is minimized. Obviously, a solution must satisfy all the problem constraints, that is the demand of the customers must be satis®ed, the vehicle capacity constraint and the working shift constraint must be satis®ed. We now give some de®nitions which will be useful in the algorithm description. A simple route ~r is de®ned as an oriented path on the graph where the ®rst and the last vertices are the depot and an intermediate facility, respectively, or an intermediate facility and the depot, or a pair of intermediate facilities. The other vertices, if any, are customers, but no customer is allowed on a simple route if the last vertex is the depot. A vehicle route is de®ned as a sequence (possibly empty, in the case a vehicle does not work in a speci®c day) of simple routes, such that the ending facility of a simple route is the starting facility of the next simple route. On a vehicle route no customer is visited more than once, the depot is visited exactly twice (the ®rst and the last vertex of the route) and each intermediate facility can be visited an arbitrary number of times. The vehicle route is the route travelled by a vehicle in a day. A daily solution is a collection of m vehicle routes for the same day. Finally, a solution is an ordered collection of T daily solutions such that each customer vi , i ˆ 1; . . . ; n, is visited according to a visiting schedule in Hi . The objective function to be minimized is the total length of the routes included in a solution which satis®es the capacity and time constraints, i.e., the total demand of the customers on any

simple route does not exceed the vehicle capacity Q and the total duration of any vehicle route does not exceed the time limit D. 3. The TS algorithm In this section, we describe the TS algorithm we designed and implemented to solve the PVRP-IF problem. As the PVRP-IF is a generalization of the PVRP, our starting point has been the TS algorithm presented by Cordeau et al. in [6] which is, to our knowledge, the best algorithm for the PVRP. The basic idea of a TS algorithm is to iteratively move from a current solution s to another solution in a neighborhood N …s† of the current solution. The neighborhood is implicitly de®ned by a set of rules that explain how to reach any solution in N …s† by means of ``small'' changes of s; such changes are called moves. The new solution is generally chosen so that a penalized cost function closely related to the objective function is minimized. Moreover, in order to avoid cycling on a local optimum, a list of forbidden solutions (tabu list) is updated at each iteration. The search space is usually larger than the set of feasible solutions, i.e., the search process might go through some infeasible, but ``promising'' solution. Here we report a general scheme of a TS algorithm (see [1] for references): · build the ®rst current solution s (initial solution); · if s is feasible, then s :ˆ s, else s ˆ ; (best feasible solution); · repeat the following steps until a stopping condition is satis®ed:  let N …s† be the neighborhood of s;  let s be the solution in N …s† which is not in a tabu list and minimizes a penalized cost function;  set the new current solution s :ˆ s and update the tabu list;  if s is feasible and its cost is less than s , then set s :ˆ s. In this section, we describe our implementation of the crucial points of such algorithm: the construction of the initial solution, the structure of the

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

solution neighborhood and the selection of a new current solution from the neighborhood which is based on the tabu list and a penalized cost function of candidate solutions. The stopping condition is satis®ed when a number g of iterations are performed. The parameter g is ®xed a priori. As already mentioned, during the searching process, the capacity and time constraints are relaxed and some infeasible solution may be selected. No relaxation is made on the structure of the solution. Each solution s in the search space is assigned a capacity and a time penalty such that both penalties are null P if and only if s is feasible. Let Pc …~ r† ˆ max…0; qi Q† bePthe capacity penalty of a simple route ~ r, where qi is the total demand of customers vi in ~ r (note that in a solution all customers are assigned a visiting schedule and their demand is well de®ned on every day). The capacity penalty Pc …s† of a solution s is de®ned as the sum of the capacity penalties Pc …~ r† of all the simple routes ~ r which compose the solution. Let Pt …r† ˆ max…0; s D† be the time penalty of a vehicle route r where s is the total duration of r (travelling time plus service times at customers and facilities on r). The time penalty of a solution s is de®ned as the sum of the time penalties Pt …r† of its vehicle routes r. The algorithm performs at any stage insertions and extractions of customers from simple routes. The removal of a customer consists in deleting the customer from the path and connecting its predecessor to its successor. The insertion of a customer i in a route is performed by removing an arc …j ; j‡ † from the route and adding the arcs …j ; i† and …i; j‡ †. The arc …j ; j‡ † is chosen so that the quantity cj i ‡ cij‡ cj j‡ is minimized. 3.1. Initial solution For the construction of the initial solution we designed and implemented procedures di€erent from the ``sweep'' technique, well known for the VRP and adopted in [6]. There are two main reasons for this choice. The ®rst is that, while in the VRP there is a single attraction point (the depot), in a problem with intermediate facilities there are more points that can be central for a route. The

237

second reason is that we looked for a more general algorithm which could be applied to non-euclidean instances like those arising in applications on urban networks. We ®rst describe a basic procedure for the construction of the initial solution I1 and then two variants I2 , I3 . Algorithm I1 assigns ®rst a visiting schedule, randomly chosen in Hi , to each customer vi . Then an iterative procedure is used to de®ne the vehicle routes on each day. Three nested loops allow the building of a complete solution. The ®rst level loop de®nes daily solutions from day 1 to day T. The second level loop de®nes vehicle routes from vehicle 1 to m and the third level loop de®nes concatenated simple routes. At any time the partial solution contains, for some days, the daily solutions and for a day the vehicle route of some vehicles and a sequence of simple routes for the vehicle for which the vehicle route is under construction. We describe how to build the next simple route for this vehicle. We ®rst see how the next facility and the ®rst customer of the route are selected, then how other customers are selected and inserted in the route. Let us de®ne C…f † the set of not-yet-visited customers for which f is the closest facility. Let f 0 be the facility (or the depot in the case the vehicle route has to be completely built) which is the starting point of the simple route and let f 00 be the facility closest to f 0 (possibly f 0 , but not the depot) such that the set of close customers C…f 00 † 6ˆ ;. The facility f 00 is the ending point of the simple route under construction. If cf 0 f 00 > 0, the ®rst customer i selected to be included in the simple route is the one which minimizes cif 0 over C…f 0 † [ C…f 00 †. If cf 0 f 00 ˆ 0, we proceed as follows. Let F be the set of all the facilities and the depot di€erent from f 0 and f 00 . Then, if F 6ˆ ;, for each customer i 2 C…f 00 † and for each facility f 2 F , the quantity Z…i; f † ˆ cif cif 00 is calculated. Z…i; f † P 0 represents a ``measure'' of how much i is attracted by f 00 with respect to f . The ®rst customer i selected to enter into the simple route is the one which maximizes minf 2F …Z…i; f †† over C…f 00 †, that is the customer which has maximum attraction to f 00 with respect to any other facility. If F ˆ ;, the

238

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

®rst selected customer i is the one which maximizes cif 0 over C…f 00 †. If the attempt to ®nd the customer i fails because no customer is available or the choice would imply some constraints violation, the procedure stops and the simple route …f 0 ; v0 ˆ depot†, which contains no customer, is appended to complete the vehicle route under construction. When the ®rst customer i of the simple route has been selected, the procedure tries to insert other customers, by identifying the customer that minimizes the insertion cost in the vehicle route under construction and does not generate any constraint violation. If no such customer is available, the simple route is complete and the procedure continues with the construction of the next simple route. Due to possible constraint violations, the previous procedure may be unable to ®nd a feasible solution. This happens when the procedure stops before all the customers have been assigned to a vehicle on each day of their chosen schedule. These customers are afterwards inserted in a route so that the insertion cost is minimized. The variant I2 follows the same scheme as I1 , but when the ®rst customer i is selected, the procedure computes the insertion cost Ci in ‰f 0 ; i; f 00 Š for each other customer i, and inserts customers in the simple route under construction according to the increasing values of Ci . If a customer generates some constraint violation, it is ignored and the construction of the simple route stops when no other customer is available. Variant I3 also follows the same scheme as I1 but, after the ®rst customer i has been selected, the customer closest to i is selected and inserted after it. The procedure keeps on appending the closest customer to the last inserted until no other customer is available. 3.2. Neighborhood structure (moves) The neighborhood of the current solution s is de®ned by a set of moves applied to s. All moves require that customers are removed from their simple routes and inserted in new routes. In the following lines we describe the four moves which de®ne the neighborhood N …s†.

Move a customer in the same day. A customer is removed from its simple route and inserted in a di€erent simple route at minimum insertion cost. Change the visiting schedule of a customer. The visiting schedule of a customer i is changed and a new visiting schedule, taken from Hi , is assigned to the customer. The new schedule may prescribe that the customer is served in some days which are in common with the current schedule. For each of the days in common, the daily solution remains unchanged. The customer is removed from all the daily solutions associated with the days of service of the current schedule which were not days of service in the new schedule. For each new day of service, the customer is inserted in the associated daily solutions at minimum insertion cost. These ®rst two moves, which we call Basic for simplicity, are taken from [6] and modi®ed with di€erent insertion and extraction heuristics. Redistribution of customers (redistribution). We consider any pair of simple routes r~1 and r~2 of any daily solution and remove all the customers belonging to the two simple routes. Then, the two simple routes are rebuilt with the same customers distributed on the two routes in a di€erent way. Consider the insertion cost of the customers in the empty routes and insert in route r~1 the customer i1 which maximizes the insertion cost in the route r~2 and vice versa insert in r~2 the customer i2 6ˆ i1 which maximizes the insertion cost in the route ~r1 . Then the following procedure is iterated until no customer has to be reinserted. Each customer i is assigned a temporary label …r; c†, where r is the route which minimizes the insertion cost of i and c is the insertion cost in the other route. Then the customer with the maximum insertion cost c is permanently inserted in the route r. The redistribution move has a complexity not greater than O…b3 †, where b is the total number of the customers in the two simple routes. Simpli®cation of intersections (intersection). With this move, an attempt is made to include in the neighborhood solutions which simplify the structure of a pair of simple routes. We consider two types of such move. The ®rst one is applied to pairs of simple routes that end at the same facility. An arc is deleted from each route. Then, the ®rst part of a route is connected to the second part of

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

the other route, while the ®rst part of the second route is connected to the second part of the ®rst route. The second type of move is applied to pairs of simple routes such that the ending facility of a route is the starting facility of the other route. Then, an arc is deleted from each route and the ®rst part of the ®rst route is connected to the ®rst part of the second route, while the second part of the ®rst route is connected to the second part of the second route. This second type of move implies an inversion of the visiting order of the customers in the second part of the ®rst route and in the ®rst part of the second route. Neither of the moves modi®es the sequence of the facilities visited by the vehicle routes of the current solution. The move has a O…b3 † complexity, where b is the total number of the customers in the two simple routes. 3.3. Selection of a solution from the neighborhood The penalized cost of a candidate solution s is given by the value of the objective function plus three penalization terms: capacity, time and diversi®cation terms. The idea of capacity and time penalization terms is taken from [6,10], while the diversi®cation term was ®rst proposed for VRP by Taillard [16]. The capacity penalization term is de®ned as aPc …s†, where a is a variable parameter. At each iteration the value of a is modi®ed. It is multiplied by a ®xed quantity …1 ‡ d† > 1, if the current solution is infeasible with respect to the capacity constraint, while it is divided by …1 ‡ d† if the current solution is feasible with respect to the capacity constraint. The time penalization term is de®ned as bPt …s† where b plays the same role as a in the capacity penalization term. It is multiplied or divided by a ®xed quantity …1 ‡ d†, if the current solution is infeasible or feasible with respect to the time constraint, respectively. Thus, if the solutions obtained during the algorithm execution remain infeasible for several iterations, the contribution of the penalization terms to the penalized cost of the solution becomes more and more relevant and an infeasible solution becomes less and less attractive.

239

The diversi®cation term aims at penalizing frequent insertions of customers in the same simple route. Each move which modi®es the current solution s into a new solution s implies that a number of customers (at least one) are moved from a simple route to a di€erent simple route. If the new potential solution has a larger value of the objective function with respect to the current solution, a diversi®cation term is added to the objective function. p P The cost of s is increased by c nmT fv~r , where c is a ®xed parameter and n, m, T are the number of customers, the number of vehicles and the number of days, respectively. The quantity fv~r is the frequency of insertion of customer v in the simple route r~, where v is a customer which is inserted in the route ~r by the move which transforms s into s. The frequency of insertion is measured over all the moves previously executed by the algorithm. When moving from the current solution s to a new solution s, a tabu list is updated. The tabu list keeps record of all the insertions of any customer i in any simple route r~ performed within the last h log…n† iterations where h is a ®xed parameter. The cost ri~r of the best feasible solution found up to the iteration at which i was inserted in ~r is stored with the pair …i; ~r†. A potential new current solution s is declared tabu ± and rejected from the selection process ± if the insertion of each customer i in a simple route r~ is in the tabu list. An exception is made if s is feasible and it costs less than ri~r for any insertion performed by the move. 4. Computational results Since no computational results are available in the literature for the PVRP-IF, we ®rst tested di€erent versions of the TS algorithm on PVRP instances taken from the literature. Then we made computational experiments on randomly generated instances speci®c for the PVRP-IF. In order to check the individual and combined in¯uence of the intersection and redistribution moves and the e€ectiveness of the diversi®cation strategy, we tested several di€erent algorithms. All the implemented algorithms use the basic moves and are named according to the following criterion. The

240

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

Table 1 Percentage deviations from Cordeau et al.'s results with 1000 iterations Instances p

P01 P02 P03 P04 P05 P06 P07 P08 P09 P10 P11 P12 P13 P14 P15 P16 P17 P18 P19 P20 P21 P22 P23 P24 P25 P26 P27 P28 P29 P30 P31 P32 PR01 PR02 PR03 PR04 PR05 PR06 PR07 PR08 PR09 PR10

Tested algorithms v

3 3 1 5 6 1 4 5 1 4 4 3 9 2 2 2 4 4 4 4 6 6 6 3 3 3 6 6 6 9 9 9 2 4 6 8 10 12 3 6 9 12

d

2 5 5 2 5 10 2 5 8 5 5 5 7 4 4 4 4 4 4 4 4 4 4 6 6 6 6 6 6 6 6 6 4 4 4 4 4 4 6 6 6 6

Average deviations

c

50 50 50 75 75 75 100 100 100 100 126 163 417 20 38 56 40 76 112 184 60 114 168 51 51 51 102 102 102 153 153 153 48 96 144 192 240 288 72 144 216 288

Computing times (seconds)

A1-D

A1RI-D

A2RI-D

A3RI-D

I1 + basic (%)

I1 + basic + redistribution + intersection (%)

I2 + basic + redistribution + intersection (%)

I3 + basic + redistribution + intersection (%)

3.40 7.61 1.60 5.89 8.12 10.74 5.73 7.67 4.33 8.39 4.41 4.52 8.70 0.00 0.00 7.69 5.75 2.95 7.06 17.23 5.75 7.00 5.00 2.70 7.21 7.39 3.70 5.83 )2.17 6.35 5.53 2.53 3.12 4.94 9.89 8.30 9.46 5.34 7.29 12.01 6.20 4.60

4.29 1.95 1.60 7.63 5.75 10.74 0.24 2.93 4.33 2.92 )1.55 5.41 9.81 0.00 0.00 1.23 1.12 0.73 )0.92 3.25 1.22 3.13 2.18 1.40 0.09 2.34 )0.23 2.62 )3.07 1.64 )0.71 2.77 0.97 1.29 1.30 1.07 5.63 2.74 2.54 2.72 1.49 2.17

5.66 1.56 1.60 10.59 3.62 7.48 2.56 3.04 4.05 3.24 )0.31 1.36 7.70 0.00 0.00 6.44 1.12 0.61 1.51 5.14 2.08 0.89 1.24 0.95 0.00 0.91 )1.73 3.02 )1.95 1.54 2.45 2.00 1.11 1.10 1.73 0.56 3.30 3.77 1.28 2.29 0.12 1.90

3.89 0.49 3.60 7.00 4.17 6.55 1.95 2.63 4.13 2.80 0.39 5.12 7.43 0.00 3.55 0.74 1.48 1.66 )0.03 4.75 3.99 4.17 2.46 2.91 4.74 5.97 )0.81 3.85 )0.94 3.42 0.96 1.30 )0.08 0.74 1.64 3.16 4.04 2.55 3.42 3.00 1.05 1.53

5.95

2.30

2.27

2.75

A1-D

15 41 22 31 84 66 39 104 71 96 146 129 1089 9 16 25 27 56 91 181 54 116 191 43 45 45 136 136 136 276 277 288 24 77 156 271 418 582 73 242 506 888

A2RI-D

32 69 24 78 173 65 128 252 72 192 268 258 1706 12 29 57 51 149 300 836 117 354 728 57 60 61 228 234 234 516 531 538 49 249 592 1104 1735 2599 195 896 2179 3922

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

letter ``A'' is a common pre®x, a ®gure ``1'', ``2'' or ``3'' indicates that the initial solution is found by means of algorithm I1 , I2 or I3 . The letters ``R'' and ``I'' mean that redistribution and intersection moves are implemented. The sux ``-D'' is used if the algorithm adopts the diversi®cation strategy. Thus, for example, algorithms A1RI-D A2RI-D and A3RI-D implement the basic moves, use the algorithms I1 , I2 and I3 , respectively, for the construction of the initial solution, include redistribution and intersection moves and adopt the diversi®cation strategy. The algorithms were implemented in C++ and run on a personal computer with a 500 MHz Pentium III processor. 4.1. PVRP instances The PVRP instances, which have been taken from [6], are named (p01±p32; pr01±pr10). Instances p01±p32 were previously available in the literature, while instances pr01±pr10 were randomly generated by Cordeau et al. We ran the three algorithms A1RI-D, A2RI-D and A3RI-D a few times and decided to ®x the following values for the parameters: d ˆ 0:5, h ˆ 7:5, c ˆ 0:5 and set the number of iterations g ˆ 1000. In order to evaluate the performance of the algorithms A1RID, A2RI-D, A3RI-D we compared our results with those obtained in [6] with their standard values for the parameters. The algorithm by Cordeau et al. was run for 15,000 iterations. They also ran their algorithm several times, with di€erent values of the parameters, on the same instance improving the results obtained with the standard values of the parameters. A comparison with these latter results is behind the scope of this paper. The results obtained are shown in Table 1. Columns 1±4 of the table show the basic characteristics of the instances: the problem name (p), the number of vehicles in the ¯eet (v), the number of days in the planning horizon (d) and the number of customers (c). In columns 5±8 the percentage deviations of A1-D, A1RI-D, A2RI-D, A3RI-D with respect to Cordeau et al.'s algorithm are given. The best result found is reported in bold for each instance. At the bottom of the table the average percentage deviation of each algorithm is reported.

241

In columns 9±10 the running times, in seconds, of algorithms A1-D and A2RI-D are indicated. Algorithms A1RI-D and A3RI-D have a running time similar to A2RI-D. The computational experiments show that the algorithms A1RI-D, A2RI-D and A3RI-D produce similar average results, especially A1RI-D and A2RI-D which are the best ones. The use of a di€erent initial solution substantially in¯uences the search space and the ®nal solution of the algorithm. The moderate average positive di€erence between our algorithms and Cordeau et al.'s algorithm may be explained by the fact that we ran a substantially smaller number of iterations, partially due to the use of a less powerful computer, and we used di€erent initial solutions. We will later show that an increase in the number of the iterations signi®cantly improves the quality of the solutions. We used the best perAlgorithm A1-D, which implements the basic moves only and uses the same initial solution as A1RI-D, gives de®nitely worse results with respect to the algorithms which include the redistribution and intersection moves. The improvement of the results obtained by A1RID, with respect to A1-D, measures the e€ectiveness of the latter moves. Table 2 Percentage deviations from Cordeau et al.'s results with 1000 and 10,000 iterations Number of iterations

Results obtained by algorithm A2RI-D instances

1000 (%)

10,000 (%)

p01 p02 p03 p04 p06 p09 p14 p15 p16 p17 p24 p25 p26 pr01

5.66 1.56 1.60 10.59 7.48 4.05 0.00 0.00 6.44 1.12 0.95 0.00 0.91 1.11

0.58 )0.44 0.25 2.79 3.33 0.99 0.00 0.00 5.20 0.00 0.56 )0.01 0.24 )0.84

2.96

0.90

Average deviations

242

Table 3 Values obtained on PVRPI-IF random instances Instances p

2 2 2 2 2 2 4 4 4 4 3 3 5 5 5 5 5 5 2 2 2 2 2 2 4 4 4 4 3 3 5 5 5 5 5 5 2

d

2 2 4 4 6 6 2 2 4 4 6 6 2 2 4 4 6 6 2 2 4 4 6 6 2 2 4 4 6 6 2 2 4 4 6 6 2

c

50 50 50 50 50 50 100 100 100 100 100 100 150 150 150 150 150 150 50 50 50 50 50 50 100 100 100 100 100 100 150 150 150 150 150 150 50

f

1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1

A2-D

A2R-D

A2I-D

A2RI-D

A2I

A2RI

I2 + basic

I2 + basic + redistribution

I2 + basic + intersection

I2 + basic + redistribution + intersection

A2I-D without diversi®cation

A2RI-D without diversi®cation

1229.95 1127.96 2076.02 2468.42 3432.08 3259.72 2107.24 2157.40 NA NA 4804.25 5109.76 2750.27 2655.53 4588.62 4625.66 6924.70 NA NA 1167.76 2148.51 2274.13 3094.62 3752.98 NA 2165.23 NA 3393.53 4934.98 4952.10 2665.75 2626.17 4639.26 4484.75 6998.73 6385.05 1132.99

1229.95 1151.99 2063.60 2459.31 3432.08 3259.72 2038.57 2155.33 NA NA 4760.77 5109.76 2750.27 2655.53 4588.62 4625.66 6947.34 NA NA 1167.76 2148.51 2274.13 3094.62 3733.46 NA 2151.42 NA 3393.53 4915.26 4974.10 2665.75 2606.58 4639.26 4484.75 7054.60 6444.63 1126.28

1195.14 1113.58 1962.96 2475.51 3270.20 3259.72 2025.62 2102.44 3116.61 3390.67 4603.25 5071.87 2711.32 2676.44 4546.82 4489.01 6767.18 6574.27 1217.62 1161.85 2072.33 2359.87 3018.16 3720.23 2052.65 2097.08 3241.42 3317.92 4711.71 4857.04 2627.53 2623.65 4600.73 4393.60 6756.53 6362.04 1071.50

1198.16 1156.99 1974.79 2475.51 3246.95 3259.72 2021.32 2106.09 3120.96 3390.67 4637.61 5071.87 2711.32 2691.14 4599.17 4489.01 6677.44 6731.23 1217.62 1161.85 2072.33 2359.87 3027.95 3720.23 2028.49 2097.08 3223.23 3317.92 4710.71 4835.77 2625.08 2603.23 4566.93 4331.42 6793.61 6311.63 1070.37

1180.76 1150.08 1960.80 2439.12 3239.67 3259.72 1994.69 2105.77 3095.80 3315.37 4585.38 5007.06 2693.51 2688.92 4549.15 4554.83 6646.26 6569.07 1197.68 1186.87 2092.91 2300.17 3015.91 3726.37 2008.59 2102.24 3245.66 3272.18 4686.76 4821.00 2598.68 2561.89 4557.57 4345.28 6703.09 6318.09 1086.65

1196.17 1150.08 1961.46 2439.12 3239.67 3259.72 1978.59 2095.45 3117.91 3313.49 4585.38 5007.06 2689.24 2688.92 4516.27 4554.83 6646.26 6567.71 1197.68 1186.87 2092.91 2300.17 3015.91 3726.37 2032.13 2102.24 3169.16 3259.21 4662.12 4812.73 2598.68 2627.04 4548.12 4345.28 6736.25 6301.85 1086.65

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

R01 R02 R03 R04 R05 R06 R07 R08 R09 R10 R11 R12 R13 R14 R15 R16 R17 R18 R19 R20 R21 R22 R23 R24 R25 R26 R27 R28 R29 R30 R31 R32 R33 R34 R35 R36 R37

Tested algorithms v

R38 R39 R40 R41 R42 R43 R44 R45 R46 R47 R48 R49 R50 R51 R52 R53 R54

2 2 2 2 2 4 4 4 4 3 3 5 5 5 5 5 5

2 4 4 6 6 2 2 4 4 6 6 2 2 4 4 6 6

50 50 50 50 50 100 100 100 100 100 100 150 150 150 150 150 150

4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4

1134.92 2196.62 1991.07 3199.64 3413.19 2103.74 2168.26 NA 3204.99 5145.96 5103.87 NA 2627.10 4796.71 4505.28 7174.81 6430.11

1134.92 2196.62 1991.07 3199.64 3424.69 2100.15 2168.26 NA 3204.99 5100.79 5103.87 NA 2620.85 4796.71 4505.28 7174.81 6424.09

1134.69 2178.63 1965.07 3101.75 3516.77 2050.45 2152.97 3186.06 3123.79 4945.41 4907.36 2673.38 2614.01 4599.79 4364.31 7019.60 6517.04

1134.69 2190.60 1965.07 3090.39 3365.77 2039.84 2152.97 3179.76 3123.79 4943.25 5029.13 2673.38 2604.76 4667.76 4366.58 7007.75 6507.58

1077.98 2187.96 1974.05 3098.77 3460.37 2067.26 2086.11 3172.46 3161.77 4947.74 4908.01 2679.50 2634.18 4610.39 4267.50 6988.23 6265.79

1077.98 2163.04 1974.05 3098.77 3460.37 2051.20 2086.11 3174.99 3161.77 4947.74 4834.97 2662.23 2579.42 4594.99 4274.67 6947.17 6265.79

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

243

We used the best performing algorithm A2RID to investigate the impact of the number of iterations on the quality of the solutions. We ran A2RI-D for 10,000 iterations on the set of instances which took less than 100 seconds of computational time. The results are listed in Table 2. Columns 2 and 3 give the results obtained with 1000 and 10,000 iterations, respectively and clearly show that the increase in the number of iterations improves the quality of the ®nal solution. We also ran algorithm A1-D for 10,000 iterations on each instance. In this case the computational time has been observed to be larger than the time required to run A2RI-D for 1000 iterations, but no relevant improvement has been achieved in terms of solution quality. The average deviation of A1-D, when run for 10,000 iterations, is 3.58%, while A2RI-D run for 1000 iterations achieves 2.27%. This con®rms the e€ectiveness of the redistribution and intersection moves. 4.2. Random PVRP-IF instances A set of 54 new instances was generated for the PVRP-IF. Table 3 summarizes the features of these instances. The columns 1±5 give the problem number (p), the number of vehicles in the ¯eet (v), the number of days in the planning horizon (d), the number of customers (c) and the number of available facilities (f). The instances were generated by the combination of 2, 3, 4 and 5 vehicles; 2, 4 and 6 days; 50, 100 and 150 customers and 1 or 4 facilities. The capacity Q and the duration D were ®xed at 55 and 400, respectively, on all instances. The customers, the depot and the facilities are spread over a square ‰0; 100Š  ‰0; 100Š. The vehicles travel through the square at speed 1 and spend 1 unit time for service at any facility. Each customer requires the same amount of product at each visit. Such demand is generated by an integer random variable uniformly distributed in ‰1; 9Š. Similarly, the service time at each customer is an integer uniformly distributed in ‰1; 9Š. In the case that four facilities are available, the square is partitioned into four identical subsquares and one facility is

244

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

located at the center of each of them, namely at points …25; 25†, …25; 75†, …75; 75†, …75; 25†. In the other case, that is when a single facility is available, the facility is located at the center of the main square that is at point …50; 50†. The depot is always located at the left-bottom corner of the main square (point …0; 0†). The customers are equally distributed among the facility subsquares. The position of a customer with respect to the assigned facility is determined as follows. A pair of coordinates are generated according to the uniform distribution on the facility's subsquare and such point is accepted with probability e dl0:05 , where d is the distance from the facility and l is the number of facilities in the instance. The set Hi of visiting schedules to be assigned to customer vi is determined as follows. Let q be the demand of customer vi each time it is visited by a vehicle. If d ˆ 2 (two days in the planning horizon), then customer vi is assigned either Hi ˆ f…q; 0†; …0; q†g (this schedule is assigned to 80% of the customers) or Hi ˆ f…q; q†g (20% of customers). If d ˆ 4, then customer vi is assigned either Hi ˆ f…q; 0; 0; 0†; …0; q; 0; 0†; …0; 0; q; 0†; …0; 0; 0; q†g (30%) or Hi ˆ f…q; 0; q; 0†; …0; q; 0; q†g (60%) or Hi ˆ f…q; q; q; q†g (10%). If d ˆ 6, then customer v is assigned either Hi ˆ f…q; 0; 0; 0; 0; 0†; …0; q; 0; 0; 0; 0†; …0;0;q; 0;0;0†; …0;0;0; q;0;0†;…0; 0;0;0;q;0†; …0;0; 0;0;0; q†g (20%) or Hi ˆ f…q; 0;0;q;0; 0†;…0;q;0; 0; q;0†; …0;q; 0;0; q;0†g (20%) or Hi ˆ f…q; 0;q;0; q;0†; …0; q;0;q; 0;q†g (50%) or Hi ˆ f…q; q;q;q;q; q†g (10%). A number of computational experiments were performed on this set of instances. We investigate the individual and combined e€ectiveness of the intersection and redistribution moves and the impact of the diversi®cation strategy on the quality of the solution. We implemented ®ve di€erent versions of algorithm A2RI-D. Algorithm A2-D implements only the basic moves, algorithm A2ID implements the basic and intersection moves and algorithm A2R-D implements the basic and redistribution moves. Finally, algorithms A2RI and A2I do not include any diversi®cation strategy, that is the parameter c is set to zero. In Figs. 1 and 2 the computational times of A2I-D and A2RI-D are plotted with respect to the product of the number of vehicles, days and customers. The costs of the solutions are reported in Table 3. The

Fig. 1.

Fig. 2.

symbol NA means that no feasible solution was found, while the best solution for each instance is reported in bold. We note that A2R-D and A2-D obtain poor results with respect to A2RI-D and A2I-D, since rarely ®nd the best solution and in eight cases out of 54 they are not able to ®nd any feasible solution. By comparing the versions with and without the redistribution move, the redistribution move seems to give a little contribution to the search of a good solution, while it is time consuming. Computational times of A2RI-D and A2I-D are compared in Fig. 1 and one can note that A2I-D is up to two times faster than A2RI-D. Computational times of A2I-D are shown in Fig. 2. Instances with one facility have the tendency to require more computational time than similar instances with four facilities, probably due to more opportunities of applying the intersection move. Computational time required by A2RI and A2I are similar to A2RI-D and A2I-D, respectively, but the absence of diversi®cation allows the algorithms to ®nd better solutions in more than 60% of the instances.

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

245

Table 4 Percentage deviations from the best solution found Instances Tested algorithms p

R01 R02 R03 R04 R05 R06 R07 R08 R09 R10 R11 R12 R13 R14 R15 R16 R17 R18 R19 R20 R21 R22 R23 R24 R25 R26 R27 R28 R29 R30 R31 R32 R33 R34 R35 R36 R37 R38 R39 R40 R41 R42 R43 R44 R45 R46 R47 R48 R49

A2-D

A2R-D

A2I-D

A2RI-D

A2I

A2RI

I2 + basic (%)

I2 + basic + redistribution (%)

I2 + basic + intersection (%)

I2 + basic + redistribution + intersection (%)

A2I-D without diversi®cation (%)

A2RI-D without diversi®cation (%)

4.17 1.29 5.88 1.20 5.94 0.00 6.50 2.96 NA NA 4.77 2.05 2.27 0.00 1.60 3.04 4.19 NA NA 0.51 3.68 0.00 2.61 0.88 NA 3.25 NA 4.12 5.85 2.90 2.58 2.51 2.00 3.54 4.41 1.32 5.85 5.28 1.55 1.32 3.54 1.41 3.13 3.94 NA 2.60 4.10 5.56 NA

4.17 3.45 5.24 0.83 5.94 0.00 3.03 2.86 NA NA 3.82 2.05 2.27 0.00 1.60 3.04 4.53 NA NA 0.51 3.68 0.00 2.61 0.36 NA 2.59 NA 4.12 5.43 3.35 2.58 1.74 2.00 3.54 5.24 2.27 5.22 5.28 1.55 1.32 3.54 1.75 2.96 3.94 NA 2.60 3.19 5.56 NA

1.22 0.00 0.11 1.49 0.94 0.00 2.38 0.33 0.67 2.33 0.39 1.29 0.82 0.79 0.68 0.00 1.82 0.10 1.66 0.00 0.00 3.77 0.07 0.00 2.19 0.00 2.28 1.80 1.06 0.92 1.11 2.41 1.16 1.44 0.80 0.96 0.11 5.26 0.72 0.00 0.37 4.49 0.52 3.21 0.43 0.00 0.04 1.50 0.42

1.47 3.90 0.71 1.49 0.22 0.00 2.16 0.51 0.81 2.33 1.14 1.29 0.82 1.34 1.84 0.00 0.47 2.49 1.66 0.00 0.00 3.77 0.40 0.00 0.99 0.00 1.71 1.80 1.04 0.48 1.02 1.61 0.41 0.00 1.35 0.16 0.00 5.26 1.27 0.00 0.00 0.00 0.00 3.21 0.23 0.00 0.00 4.02 0.42

0.00 3.28 0.00 0.00 0.00 0.00 0.81 0.49 0.00 0.06 0.00 0.00 0.16 1.26 0.73 1.47 0.00 0.02 0.00 2.15 0.99 1.15 0.00 0.17 0.00 0.25 2.41 0.40 0.53 0.17 0.00 0.00 0.21 0.32 0.00 0.26 1.52 0.00 1.15 0.46 0.27 2.81 1.34 0.00 0.00 1.22 0.09 1.51 0.65

1.31 3.28 0.03 0.00 0.00 0.00 0.00 0.00 0.71 0.00 0.00 0.00 0.00 1.26 0.00 1.47 0.00 0.00 0.00 2.15 0.99 1.15 0.00 0.17 1.17 0.25 0.00 0.00 0.00 0.00 0.00 2.54 0.00 0.32 0.49 0.00 1.52 0.00 0.00 0.46 0.27 2.81 0.56 0.00 0.08 1.22 0.09 0.00 0.00

246

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

Table 4 (Continued) Instances Tested algorithms p

R50 R51 R52 R53 R54

A2-D

A2R-D

A2I-D

A2RI-D

A2I

A2RI

I2 + basic (%)

I2 + basic + redistribution (%)

I2 + basic + intersection (%)

I2 + basic + redistribution + intersection (%)

A2I-D without diversi®cation (%)

A2RI-D without diversi®cation (%)

1.85 4.39 5.57 3.28 2.62

1.61 4.39 5.57 3.28 2.53

1.34 0.10 2.27 1.04 4.01

0.98 1.58 2.32 0.87 3.86

2.12 0.34 0.00 0.59 0.00

0.00 0.00 0.17 0.00 0.00

1.16

1.17

0.58

0.45

Average deviations

In Table 4, the percentage deviations with respect to the best solution found for every instance are reported for each algorithm. At the bottom of the table the average deviations are given for algorithms A2I-D, A2RI-D, A2I and A2RI. Due to the large number of failures in ®nding a feasible solution, the average deviations of A2-D and A2R-D are not summarized. Nevertheless, it is clear that A2-D and A2R-D achieve poor results even on those instances where feasibility is reached. Table 5 provides information on the algorithms performance by directly comparing pairs of versions. It gives for each pair the number of instances on which the ®rst version has given a better than, equal to or worse result than the second version. In particular, if we compare the algorithms with and without the redistribution move, we see that such move improves the solution on a large number of instances, but in spite of this the average quality of the solutions does not improve signi®cantly (see Table 4). Moreover, in Table 5 we note that algorithms which di€er only in the implementation of the redistribution move, often reach the same cost of the solution. On the other hand if two algorithms di€er in the use of the diTable 5 Direct comparison of algorithms on 54 instances a () b

a
aˆb

a>b

A2RI-D () A2I-D A2RI-D () A2RI A2I-D () A2I A2RI () A2I

21 15 20 19

18 1 1 27

15 38 33 8

versi®cation strategy, they almost never get to the same solution. Summarizing the conclusions which can be drawn from the computational results, we can say that the intersection move is essential to obtain feasible and good solutions, while the redistribution move may be of help in some cases only. The absence of a diversi®cation strategy substantially improves the quality of the solutions and does not increase the computational times.

5. Conclusions In this paper, we de®ned and studied a PVRP with intermediate facilities which has applications in collection problems, such as those arising in waste collection, in reverse logistics and in distribution problems, where the starting and ending point of the tour of each vehicle is a depot but the goods are stored in one or several warehouses distributed on a geographic area (this is the case for instance of a transportation company). We designed and implemented a TS algorithm which extends the set of types of moves with respect to other TS algorithms designed for the VRP and the PVRP. The computational results show that the new moves are e€ective, as they allow an improvement on the quality of the solutions obtained when compared with the algorithm which only implements the basic moves on the same number of iterations. Moreover, we show that the algorithm with the only basic moves cannot reach the performance of the algorithms with the new

E. Angelelli, M. Grazia Speranza / European Journal of Operational Research 137 (2002) 233±247

moves, even when run for a larger number of iterations. References [1] E.H.L. Aarts, J.K. Lenstra (Eds.), Local Search in Combinatorial Optimization, Wiley, New York, 1997. [2] E.J. Beltrami, L.D. Bodin, Networks and vehicle routing for municipal waste collection, Networks 4 (1974) 65±94. [3] L.D. Bodin, B.L. Golden, A.A. Assad, M.O. Ball, Routing and scheduling of vehicles and crews, the state of the art, Computers & Operations Research 10 (2) (1983) 63±211. [4] I-M. Chao, B.L. Golden, E.A. Wasil, An improved heuristic for the period vehicle routing problem, Networks 26 (1995) 25±44. [5] N. Christo®des, J.E. Beasleym, The period routing problem, Networks 14 (1984) 237±256. [6] J.-F. Cordeau, M. Gendreau, G. Laporte, A tabu search heuristic for periodic and multi-depot vehicle routing problems, Networks 30 (1997) 105±119. [7] M. Fisher, Vehicle routing, in: M.O. Ball, et al. (Eds.), Network Routing: Handbooks in Operations Research and Management Science, vol. 8, North-Holland, Amsterdam, 1995.

247

[8] M.L. Fisher, R. Jaikumar, A generalized assignment heuristic for vehicle routing, Networks 11 (1981) 109±124. [9] M. Gaudioso, G. Paletta, A heuristic for the periodic vehicle routing problem, Transportation Science 26 (2) (1992) 86±92. [10] M. Gendreau, A. Hertz, G. Laporte, A tabu search heuristic for the vehicle routing problem, Management Science 40 (1994) 1276±1290. [11] G. Laporte, The vehicle routing problem: An overview of exact and approximate algorithms, European Journal of Operational Research 59 (1992) 345±358. [12] G. Laporte, I.H. Osman, Routing problems: A bibliography, Annals of Operations Research 61 (1995) 227±262. [13] Laporte, G., 1997. In: Dell'Amico, M., Maoli, F., Martello, S. (Eds.), Vehicle Routing, Annotated Bibliographies in Combinatorial Optimization, Wiley, New York, pp. 223±240. [14] R.A. Russell, D. Gribbin, A multiphase approach to the period routing problem, Networks 21 (1991) 747±765. [15] R.A. Russell, W. Igo, An assignment routing problem, Networks 9 (1979) 1±17. [16] E.D. Taillard, Parallel iterative search methods for vehicle routing problems, Networks 23 (1993) 661±673. [17] C.C.R. Tan, J.E. Beasley, A heuristic algorithm for the period vehicle routing problem, OMEGA International Journal of Management Science 12 (5) (1984) 497±504.