Available online at www.sciencedirect.

com

European Journal of Operational Research 195 (2009) 791–802

www.elsevier.com/locate/ejor

A variable neighborhood search heuristic for periodic routing

problems
Vera C. Hemmelmayr, Karl F. Doerner *, Richard F. Hartl
Department of Business Administration, University of Vienna, Bruenner Strasse 72, 1210 Vienna, Austria
Received 16 January 2007; accepted 8 August 2007
Available online 13 November 2007

Abstract

The aim of this paper is to propose a new heuristic for the Periodic Vehicle Routing Problem (PVRP) without time windows. The
PVRP extends the classical Vehicle Routing Problem (VRP) to a planning horizon of several days. Each customer requires a certain num-
ber of visits within this time horizon while there is some flexibility on the exact days of the visits. Hence, one has to choose the visit days
for each customer and to solve a VRP for each day. Our method is based on Variable Neighborhood Search (VNS). Computational
results are presented, that show that our approach is competitive and even outperforms existing solution procedures proposed in the
literature. Also considered is the special case of a single vehicle, i.e. the Periodic Traveling Salesman Problem (PTSP). It is shown that
slight changes of the proposed VNS procedure is also competitive for the PTSP.
© 2007 Elsevier B.V. All rights reserved.

Keywords: Periodic vehicle routing problem; Periodic traveling salesman problem; Metaheuristics; Variable neighborhood search

1. Introduction tics are developed by Christofides and Beasley (1984), Tan

and Beasley (1984) and Russel and Gribbin (1991). Gaud-
Vehicle Routing Problems (VRPs) have received consid- ioso and Paletta (1992) present a heuristic for the PVRP
erable attention both in theoretical research and in real with the objective to minimize the number of vehicles.
world applications. Forming the basis of many routing Chao et al. (1995a) provide a two phase heuristic. To
models, they have been extended in various directions. In obtain an initial solution they solve an integer linear pro-
this paper we focus on the Periodic Vehicle Routing Prob- gram to assign visit day combinations to the customers.
lem (PVRP) in which a planning period of several days is In a second phase, they use several improvement operators
considered and customers must be visited more than once. while they relax the capacity of the vehicles. When getting
Different customers usually require different numbers of stuck, re-initializations are performed.
visits in a certain time horizon. Customers with larger Cordeau et al. (1997) propose a tabu search heuristic for
demands or smaller storage capacities require more visits the PVRP that can also be used to solve the Multi-Depot
than customers with smaller demands or larger storage Vehicle Routing Problem and the Periodic Traveling Sales-
capacities. This type of problem occurs e.g. in grocery dis- man Problem (PTSP). The neighborhood consists of mov-
tribution (Carter et al., 1996), soft drink industry (Golden ing a customer from one route to another route of the same
and Wasil, 1987), waste collection (Beltrami and Bodin, day or assigning a new visit combination to a customer.
1974; Russel and Igo, 1979) and others. Insertions and removals of customers are performed using
Early heuristics for the PVRP are proposed by Beltrami the GENI operator (Gendreau et al., 1992). The tabu
and Bodin (1974) and Russel and Igo (1979). Other heuris- search algorithm allows for infeasible solutions during
the search process using an adaptive penalty function.
*
Recently a scatter search procedure was developed by Ale-
Corresponding author. Tel.: +43 1 427738113; fax: +43 1 427738094. E- gre et al. (2007) for solving a problem of periodic pick-up
mail addresses: [email protected] (V.C. Hemmel- mayr), of raw materials for a manufacturer of auto parts. They
[email protected] (K.F. Doerner), Richard.Hartl@
univie.ac.at (R.F. Hartl). use a two-phase approach, that first assigns orders to days

0377-2217/$ - see front matter © 2007 Elsevier B.V. All rights

reserved. doi:10.1016/j.ejor.2007.08.048
79 V.C. Hemmelmayr et al. / European Journal of Operational Research 195 (2009) 791–802
2
and then constructs the routes of each day. Finally approach and compares it to state of the art metaheuristics.
Drummond et al. (2001) proposed parallel genetic Finally, Section 5, concludes the paper.
algorithms.
Furthermore special implementations for real-world
problems were provided by Hadjiconstantinou and Baldac- 2. Solution procedure for the PVRP
ci (1998) who propose a heuristic for a Multi-Depot Period
Vehicle Routing Problem that arises in the utility sector. We first describe our VNS algorithm for the PVRP. In Section
Angelelli and Speranza (2002) suggest a tabu search for a 3 we will propose minor adaptations of this method in
special application, a periodic vehicle routing problem with order to be used for solving the PTSP. The basic idea of
intermediate facilities, where vehicles can replenish their VNS is a systematic change of neighborhoods within a
capacity at intermediate facilities. Blakeley et al. (2003) local search procedure. Starting from any initial solution,
developed an optimization system relying on several algo- a so-called shaking step is performed by randomly selecting
rithms for planning the maintenance of escalators and a solution from the first neighborhood. This is followed by
elevators. applying an iterative improvement algorithm. This proce-
The PTSP is a special case of the PVRP where only one dure is repeated as long as a new incumbent solution is
vehicle is available every day and tour length or duration found. If not, one switches to the next neighborhood
(which is typically larger) and performs a shaking step fol-
constraints are not considered. A mathematical formula-
lowed by the iterative improvement. If a new incumbent
tion of the PVRP and the PTSP can be found in Cordeau solution is found one starts with the first neighborhood;
et al. (1997). Heuristics for the PTSP are provided by otherwise one proceeds with the next neighborhood, etc.
Christofides and Beasley (1984) and Paletta (1992). The steps of the basic VNS are shown in Fig. 1, where
These earlier solution techniques are outperformed by Nj (j = 1, ... , jmax) is the set of neighborhoods. The stop-
more recent metaheuristic approaches. Chao et al. ping condition can be a limit on CPU time, a limit on the
(1995b) start from an initial solution and exchange visit number of iterations, or a limit on the number of iterations
day combinations by using the record-to-record approach between two improvements. See Mladenovic´ and Hansen
of Dueck (1993). New solutions are accepted if their cost (1997) and Hansen and Mladenovic´ (2000, 2001) for a
is below a specified threshold, being the cost of the best more thorough description of VNS (see Fig. 2).
solution found plus a certain deviation. Local Search and re- In the next subsection, we describe the different compo-
initializations are performed afterwards. Cordeau et al. nents of the VNS implemented for the PTSP and the
(1997) develop a tabu search method based on the GENI
PVRP.
operator (Gendreau et al., 1992). They apply their method
also to the PVRP and to the multi-depot vehicle routing
problem. Paletta (2002) proposes an improvement proce- 2.1. Initial solution
dure within a tour construction procedure and a new ver-
sion of the improvement procedure is proposed by For obtaining an initial solution each customer is
Bertazzi et al. (2004). assigned a visit day combination randomly. Routes are
We develop another metaheuristic solution approach for constructed by solving a vehicle routing problem for each
the PVRP and the PTSP that is based on Variable Neigh- day using the Clarke and Wright savings algorithm (Clarke
borhood Search (VNS). VNS is a local search based meta- and Wright, 1964). The Clarke and Wright savings algo-
heuristic first proposed by Mladenovic´ (1995), Mladenovic´ rithm terminates when no two routes can feasibly be
and Hansen (1997) and Hansen and Mladenovic´ (1997). merged, i.e., no two routes can be merged without violating
The VNS approach has already been successfully applied the route duration or capacity constraints. As a result, the
to other variants of the VRPs (see e.g. Braysy, 2003; Pola- number of routes may exceed the number of available vehi-
cek et al., 2004, 2007). However, to the best of our knowl- cles. In that case, a route with the fewest customers is iden-
edge VNS has not been applied to periodic routing tified and the customers in this route are moved to other
problems so far. routes (minimizing the increase in costs). Note that this
The paper has two main contributions. First, from a may result in routes that no longer satisfy the duration
technical point of view, it presents the first application of or capacity constraints. This step is repeated until the num-
a VNS to periodic routing problems. Second, from a prob- ber of routes is equal to the number of vehicles. Since the
lem oriented point of view the computational results show initial solution may not be feasible the VNS needs to incor-
that the approach is competitive with the existing tech- porate techniques that drive the search to a feasible solu-
niques. The developed algorithm is simple, flexible and tion. See Fig. 2 for a more detailed description of the
accurate and yields several new best solutions. procedure for constructing the initial solution.
The remainder of the paper is organized as follows. In
Section 2 we describe the proposed algorithm for the PVRP
and in Section 3 we discuss the appropriate adaptations for 2.2. Shaking
applying it to the PTSP. In Section 4, we present a compu-
tational study that analyzes the proposed solution The set of neighborhoods used for shaking is at the heart
of the VNS. Each neighborhood should strike a proper
 Fig. 1. Steps of the VNS (cf. Hansen and Mladenovic´, 2001).

Fig. 2. Build initial solution.

balance between perturbing the incumbent solution and route one. In our algorithm we relocate up to three custom-
retaining the good parts of the incumbent solution. ers. The cross exchange operator exchanges two segments
Two popular and effective neighborhoods for vehicle of different routes. Fig. 4 shows that the segment from cus-
routing are based on the move and the cross-exchange tomer x01 to y1 of route one is exchanged with the segment
operators. A classification of move and cross-exchange x20 to y2 of route two. We consider a segment length of up to
can be found in Van Breedam (1994) and Kindervater six customers. The orientation of the segment(s) and of the
and Savelsbergh (1997). The operator move inserts a seg- route(s) is preserved by the move and cross-exchange oper-
ment of one route into a different route. For example in ators. The move and cross-exchange operators are used to
Fig. 3 customers x02 and y2 are moved from route two to define a set of neighborhoods that allow the exploration of

Fig. 3. The Move operator. Fig. 4. The cross exchange operator.

Table 1
2.4. Acceptance decision
Set of neighborhood structures with jmax = 15
j Operator Min. customers Max. customers
After the shaking and the local search procedures have
1 Change combination 1 1 been performed, the solution thus obtained has to be com-
2 Change combination 1 2 pared to the incumbent solution to be able to decide
3 Change combination 1 3 whether or not to accept it. The acceptance criterion in
4 Change combination 1 4
5 Change combination 1 5
the basic VNS is to accept only improvements. However
6 Change combination 1 6 that way the search can easily get stuck in a local optimum.
Thus in most cases it is essential to also have a strategy of
Min. segment Max. segment length accepting non-improving solutions under certain condi-
7 Move length 1 min(1, n) tions. We implement a scheme that is inspired by Simulated
8 Move 1 min(2, n) Annealing (SA) (Kirkpatrick et al., 1983). Hence, our
9 Move 1 min(3, n) method could be considered a hybrid of VNS and SA.
1 Cross 1 min(1, n) However since the SA part is rather small we prefer to
10 Cross 1 min(2, n) regard it as a VNS. More specifically, improving solutions
1 Cross 1 min(3, n)
are always accepted and inferior solutions are accepted
12 Cross 1 min(4, n) —ðf ðx00Þ—f ðxÞÞ
31 Cross 1 min(5, n) with a probability exp where f(x) is the cost
14 Cross 1 min(6, n) of solution x possibly including penalty costs if it is infeasible.
T
5 The acceptance of inferior solutions depends on a given
temperature T and the difference between the costs of the
increasingly distant solutions from the incumbent to over- new solution and the incumbent solution. The temperature
come local optimality and strive for global optimality. The T is linearly decreased in g/k stages during the search pro-
metric to measure the increasing size of a neighborhood is cess, where g represents the total number of iterations exe-
given by the maximum number of customers in the route cuted. Thus, every k iterations T is decreased by an amount
T mk
segments used within the operators. The cross-exchange . Different cooling schedules like exponential cooling and
operator is shown in Fig. 4. a constant temperature have also been considered but the
For the PVRP it is essential to also have a neighborhood linear
g annealing scheme provided the best results.
that changes the visit combinations for customers. We use An alternative strategy would be the so called Skewed
neighborhoods in which the visit combinations of a limited VNS, an extension of the basic VNS proposed by Hansen
number of customers are changed. For each of these cus- and Mladenovic´ (2000). In this approach a solution is
tomers a visit combination is chosen randomly. Here, the not only evaluated by its objective value but also by its dis-
metric to measure the increasing size of a neighborhood tance to the incumbent solution, favoring more distant
is given by the maximum number of customers for which solutions. Let the function q(x, x00) measure the distance
the visit day combination is changed. Table 1 shows the between the incumbent solution x and the new solution
maximum segment length considered for each neighbor- x00. A new solution is accepted if f(x00) — aq(x, x00)< f(x).
hood j. Another approach to accept non-improving solutions is
In each neighborhood all the possible segment lengths based on Threshold Accepting (TA). A solution yielding an
and numbers of customers are equally likely to be chosen. improvement is always accepted. Moreover ascending
Hence our choice of neighborhoods is biased toward smal- moves are accepted after a minimum number of iterations
ler changes to focus the search rather close to the incum- counted from the last accepted move, but only if the cost
bent solution. increase is below a certain threshold. We implemented
these three possibilities (SA, skewed, TA). In the Skewed
2.3. Local search VNS (SVNS) approach we measure the distance q(x, x00)
by using the number of customers that are exchanged in
A solution obtained through shaking is submitted to a the move and cross operator and the number of routes that
local search procedure to come up with a locally optimal are changed in the change combination operator. Our
solution. We apply one of the most popular iterative implementation of the TA approach is based on the one
improvement procedures, namely 3-opt, which was intro- described in Polacek et al. (2004). Computational experi-
duced by Lin (1965). This heuristic tries all shifts of some ments show that SA delivers 2.71% better solutions than
subsequence to different positions in the same route. More SVNS and 3.61% than TA. In what follows, we will only
precisely, three edges are deleted and replaced by three report results based on the SA acceptance criterion.
other edges. The tour is 3-optimal, when it cannot be As mentioned at the start of this section, the VNS has to
improved by such a change. In our algorithm 3-opt with- be able to handle infeasible solutions. Infeasibility occurs if
out sequence inversion, also often denoted as 3-opt*, is the total capacity or tour duration exceed the specified lim-
used. Only individual routes are improved, so that only its. We use a weighted, linear penalty function for viola-
the routes that have changed during shaking have to be re- tions of this constraint. This penalty function is added to
optimized. The local search restarts immediately after the objective function before the solution is evaluated for
an improving move was found.
acceptance. The weights are adjusted dynamically. If the restricted. The 10 instances of the new data set were pro-
total capacity or tour duration of any tour is exceeded vided by Cordeau et al. (1997). The old data were solved
the respective weight is increased, if it is feasible the weight by Christofides and Beasley (1984) (CB), Tan and Beasley (1984)
is decreased. However the weights can only be adjusted (TB), Russel and Gribbin (1991) (RG), Chao et al. (1995a)
within predefined upper and lower bounds. Hence if a tour (CGW), Cordeau et al. (1997) (CGL) and Alegre et
is infeasible, but the weight would exceed the upper bound al. (2007) (ALP) and we compare our results with their results.
it will not be increased any more and vice versa. The weight Results from Drummond et al. (2001) were not taken into
is initialized with its upper bound in order to lead the account because according to Alegre et al. (2007) they
search toward feasible solutions in the beginning. are not comparable. The algorithms CGL and ALP
delivered the best solution values so far and CGW delivers
3. Solution procedure for the PTSP some ties. Results for the ‘‘new data” set were only
given by Cordeau et al. (1997) (CGL). A descrip- tion of the
The solution procedure for the PTSP is based on the one different instances of the old data set can be found in Table
for the PVRP. But in order to solve the PTSP more effi- 2 and for the new data set can be found in Table 3.
ciently some adaptations are appropriate. The algorithm
for the PVRP uses Clarke and Wright savings algorithm
(Clarke and Wright, 1964) for building an initial solution. 4.1.2. Parameter settings
For the PTSP the savings measure is irrelevant and best Compared to other metaheuristics only few parameters
insertion is used to build a starting solution. have to be determined and tuned. In our implementation
In the PVRP algorithm the shaking phase is composed these are the initial temperature for accepting deteriorating
of the operators move, cross and change combination,
where move and cross are used inter route only, i.e. cus-
tomers are exchanged between routes. The shaking phase Table 2
Instance description of the ‘‘old data” set for the PVRP, where n is the
in the PTSP is similar, but the operators move and cross number of customers, m is the number of vehicles that can be used, t is the
are now used intra route because there is only one route number of days in the planning horizon, D is the maximum duration of a
for every day. route, Q is the maximum capacity of the vehicles, and fi is the number of
For the PVRP the local search phase consists of 3-opt. customers that must be visited i times
But typically, in the standard benchmark instances each Instance n m t D Q Service frequencies
PTSP tour consists of considerably more customers than
f1 f f f f f
a PVRP tour. This is why an efficient implementation for 2 3 4 5 6
the PTSP should employ a faster local search. In order to v-p01 5 3 2 1 5
v-p02 50 3 5 16 10 2 7
save computation time the 2-opt operator is applied for 05 61 75 6
v-p03 1 5
the PTSP instead of 3-opt as used for the PVRP. 2-opt v-p04 70 2 5 16 70
was introduced by Croes, 1958. This operator deletes two v-p05 57 6 5 41 53 3 1
edges of a tour and reconnects those paths in the other pos- v-p06 75 1 10 14 70 4 1
sible way. The local search restarts immediately after an v-p07 105 4 2 42 15
improving move was found. v-p08 0
10 5 5 20 04 4 1
Chao et al. (1995b) imposed the constraint that the trav- v-p09 010 1 8 02 10 6 4
eling salesman has to visit at least one city each day. Our v-p10 0
10 4 5 20 04 4 1
v-p11 013 4 5 02 10 62 41 1 1
algorithm handles this by penalizing an empty day with a
v-p12 9
16 3 5 13 10 28 27
constant penalty. The penalty is added to the objective 341
v-p13 9 7 204 43 4
function before the solution is evaluated. 72 00 2 78 08
v-p14 2 4 4
v-p15 03 2 4 03 1 1 6
4. Computational experiments v-p16 58 2 4 40 26 26 8
v-p17 64 4 4 02 41 41 8
4.1. Results on PVRP instances v-p18 70 4 4 30 36 36 1
v-p19 116 4 4 04 24 24 21
v-p20 2
18 4 4 60 8 8 26
4.1.1. Test instances v-p21 46 6 4 02 02 02 41
We tested our algorithm on instances taken from the lit- v-p22 110 6 4 30 4 4 12
erature. There are two data sets available. The so-called v-p23 416 6 4 04 87 87 82
‘‘old data” set contains 32 instances. Instances p1–p10 v-p24 85 3 6 20 32 29 4 6
were proposed by Eilon et al. (1971) for the VRP and adapted v-p25 15 3 6 02 63 9 6
to the PVRP by Christofides and Beasley (1984). Russel v-p26 51 3 6 20 36 9 6
and Igo (1979) proposed instance p11. Instance p12 v-p27 101 6 6 02 67 1 1
v-p28 2
10 6 6 20 72 18 12
and p13 are taken from Russel and Gribbin (1991). 2
v-p29 10 6 6 20 72 18 12
Instances p14–p32 were introduced by Chao et al. (1995). In v-p30 2
15 9 6 20 12 28 12
the lat- ter data set D = 1, which means that tour duration is v-p31 3
15 9 6 20 10 27 18
not v-p32 3
15 9 6 20 10 27 18
3 0 0 7 8
Table 3 Table 4
Instance description of the ‘‘new data” set for the PVRP, where n is the Results for PVRP old data, compared to Tan and Beasley (1984) (TB),
number of customers, m is the number of vehicles that can be used, t is the Christofides and Beasley (1984) (CB), Russel and Gribbin (1991) (RG),
number of days in the planning horizon, D is the maximum duration of a Chao et al. (1995) (CGW), Cordeau et al. (1997) (CGL) and Alegre et al.
route, Q is the maximum capacity of the vehicles, and fi is the number of (2007) (ALP)
customers that must be visited i times
Instance TB CB RG CGW CGL ALP VNS
Instance n m t D Q Service frequencies v-p1 – 54 537 524 524. 531. 524.
f f f f f6 v-p2 148 7.4
144 .3
1355 .6
1337 61
1330. 1324.
02 61
1332.
v-p3 1.3
– 3.1
54 .4
– .2524 09524. 74537. 01528.
v-pr01 4 2 4 5 2 1 24 2 12 3 4 12
v-p4 – 6.7
84 867 .6
860 61
837. 37
845. 97
847.
v-pr02 98 4 4 40 10 48 24 24
v-p5 219 3.9
218 .8
2141 .9
2089 93
2061. 97
2043. 48
2059.
v-pr03 146 6 4 48 19 72 36 36
4 v-p6 –2.5 7.3
938.2 –.3 881.1 36
840.3 75
840.1 74
884.69
v-pr04 19 8 4 46 19 9 4 4
2 v-p7 – 839.2 833.6 832 829.37 829.65 829.92
v-pr05 24 1 4 4 18 16 68 68
0 v-p8 2281.8 2151.3 2108.3 2075.1 2054.9 2052.21 2058.36
v-pr06 28 10 4 42 18 12 70 70
87 23 v-p9 – 875 – 829.9 829.45 829.65 834.92
v-pr07 6 50 27 41 12 1 2 1
v-p10 1833.7 1674 1638.5 1633.2 1629.96 1621.21 1629.76
v-pr08 142 6 6 40 10 38 38 38 38
4 v-p11 878.5 847.3 820.3 791.3 817.56 782.17 791.18
v-pr09 21 9 6 47 19 56 56 56 56
6 v-p12 – – 1312 1237.4 1239.58 1230.95 1258.46
v-pr10 28 1 6 45 18 74 74 74 74
8 2 2 7 2 2 2 2 v-p13 – – 3638.1 3629.8 3602.76 – 3835.9
v-p14 – – – 954.8 954.81 954.81 954.81
v-p15 – – – 1862.6 1862.63 1862.63 1862.63
solutions and the initial values for the penalty terms in the v-p16 – – – 2875.2 2875.24 2875.24 2875.24
objective function. Several experiments were made to find a good v-p17 – – – 1614.4 1597.75 1597.75 1601.75
initial temperature T for SA. It turned out that the initial v-p18 – – – 3217.7 3159.22 3157 3147.91
temperature should be higher the larger the average distance v-p19 – – – 4846.5 4902.64 4846.49 4851.41
between customers is. In order not to have differ- ent v-p20 – – – 8367.4 8367.4 8412.02 8367.4
initial temperatures for all instances we group the instances v-p21 – – – 2216.1 2184.04 2173.58 2180.33
into those with large average distances between customers v-p22 – – – 4436.4 4307.19 4330.59 4218.46
v-p23 – – – 6769 6620.5 6813.45 6644.93
(p27–p32) and small average distances (all other instances). The
v-p24 – – – 3773 3704.11 3702.02 3704.6
initial temperature is set to 125 in the for- mer case and v-p25 – – – 3826 3781.38 3781.38 3781.38
to 7 in the latter case. The temperature is decreased every v-p26 – – – 3834 3795.32 3795.33 3795.32
1000 iterations, in a way that it becomes 0 in the last v-p27 – – – 23401.6 23017.45 22561.33 22153.31
1000 iterations. The weights for penalizing tour length and v-p28 – – – 23105.1 22569.4 22562.44 22418.52
duration violations are adjusted dynamically. They are set to v-p29 – – – 24248.2 24012.92 23752.15 22864.23
a value of 1000 in the beginning and multi- plied by a factor v-p30 – – – 80982.1 77179.33 76793.99 75579.23
1.001 if the solution is infeasible or divided by this factor v-p31 – – – 80279.1 79382.35 77944.79 77459.14
if the solution is feasible. But these adapta- tions are only v-p32 – – – 83838.7 80908.95 81055.52 79487.97
applied if the weights remain in the interval between 10 and
1000.
have frequency one. Since all customers need to be visited
4.1.3. Numerical results only once and this visit can take place on any day of the
The algorithm was coded in ANSI C and compiled with planning period, these instances reduce to simple VRPs
the GNU C compiler version 3.4.4. Our experiments were (or more precisely multiple TSPs with given vehicle num-
performed on a PC with 3.2 GHz. Preliminary tests were ber). Moreover v-p3, v-p6 and v-p9 have only one vehicle
performed in order to verify that all 15 neighborhoods con- available. It can be seen that our algorithm delivers com-
tribute to the performance of the algorithm. Clearly neigh- petitive results also on these degenerate instances but out-
borhoods with lower number found new incumbent performs the other algorithms on instances that have
solutions more often since they are used more frequently. higher visit frequencies, i.e. instances v-p14–v-p32. Espe-
However all neighborhoods found a significant number cially for larger instances our algorithm provides very good
of new incumbent solutions (see Table 15 below). results.
Table 4 reports results for the old data set. It shows the As the algorithms CGL and ALP provided the best
objective function values for CB, TB, RG, CGW, CGL, results on the old data set so far, we report the develop-
7
ALP and our algorithm with 10 iterations (in Table 7 ment of our solution approach in comparison to these
7
we will show that the run times for 10 iterations are com- two algorithms. Table 5 shows average results over 10 runs
parable with those of the best benchmark algorithms). The 6 7 8 9
with 10 , 10 , 10 and 10 iterations and compares them to
best results are always marked in bold. the results delivered by CGL and ALP. The last 2 rows
Our algorithm was developed and tuned for solving show the average per cent differences between our results
truly PVRPs. In the old data set, however there are some and CGL and ALP, respectively. More precisely, we com-
degenerate instance that are not really periodic; e.g. in puted the per cent difference of the average of 10 runs to
instances v-p1, v-p3, v-p4, v-p6, v-p7 and v-p9 all orders the benchmark approaches for each instance and reported
the average aver all instances. It is seen that after a very
 Table 5
Results for PVRP old data, compared to CGL (Cordeau et al., 1997) and ALP (Alegre et al., 2007
instancev-pr10 with 288 customers improved solutions of
_4.83%were possible. After 108 iterations already for all of the10
instances improved results were found on average. Theresults
confirm the findings for the old data. For truly periodicproblems
and in particular for larger instances our
algorithm performs very well.

Table 7 reports computation times on the instances of

the old data set for the so far best approaches CGL (Cordeauet al.,
1997) and ALP (Alegre et al., 2007) and run
times for 107 iterations of VNS. All run times are given
in seconds. If algorithms are run on different machines, a
direct comparison of computation times is always difficult.
To make the results comparable at least to the CGL algorithm we
ran the code of CGL on the same machine as
VNS which yielded considerably shorter times compared
to those originally reported by CGL. Although the results
by ALP were published in 2007, they use a rather old
machine which is – according to Dongarra (2005) – 7.17
times slower than our machine. The total computation time
44914.6 must therefore be corrected to 6264.24 seconds,
which is still larger than for CGL and our VNS.

Table 8 reports run times for the instances of the new

data set of CGL and 107 iterations of VNS. Again both
algorithms were run on the same machine. T refers to the
total time used to execute the algorithm and T* is the time
needed to obtain the best solution during search process. It
shows that the VNS algorithm is competitive to CGL and is faster
in large instances.

short running time, i.e. 106 iterations, the VNS is still about
1.5% worse. When comparing VNS after 107 iterations to
CGL and ALP the solution quality is already _0.17%
and _0.21% better, respectively. The results still improve
with increasing number of iterations. It should be noted
that these average values also include the results for the

degenerate (non-periodic) instances. If we would omit these or

adapt the VNS to these special situations the average results
would be better.
We applied our algorithm also to the new instances.
Table 6 shows the results for the new data after 106, 107,
108 and 109 iterations compared to CGL. We present average
results over 10 runs and we compare it to the results of the CGL
(1 run). Note that these instances were only solved by CGL so
far. Applying our algorithm 107 iterations we improved the
results by _0.92%. When applyingthe algorithm 109 iterations
the improvement is almost_2%. On some small instances there
are only moderateimprovements e.g. instance v-pr07 with 72
customers animprovement of _0.09% is possible whereas for
Table 7
seen in Tables 9 and 10. For the old data set 27 out of 32
PVRP old data – Reported run times in seconds executed on a Pentium III
600 MHz (ALG), and on a 3.2 GHz PC (CGL and VNS)
instances are improved or equal results were found. For
the new data set 9 out of 10 instances are improved and
Instance T T*
there is one tie. The best known results of the new data
ALP CGL VNS CGL VNS
have been improved by —1.50%.
v-p01 26 29.4 98.3 6.6 4
v-p02 49 8 35.4 81.6 11.4 9.
4
v-p03 445 32.4 100.5 13.2 67. 4.2. Results on PTSP instances
v-p04 1 46.8 67.2 2 2.
4
v-p05 14 52.8 6 432.4 59. 4.2.1. Test instances
v-p06 21 57 87 51.6 76
7 6
Our algorithm was tested with the standard benchmark
v-p07 19 76.8 183.2 1 147.
v-p08 39 123.6 142.9 5
119.4 812 instances proposed in the literature. Instances t-p1 to t-p10
v-p09 5
97 95.4 193.1 34.2 6.5
16 were given by Eilon et al. (1971) for the VRP and adapted
v-p10 90 123.6 170 111.6 2.1
15 to the PTSP by Christofides and Beasley (1984). Instances t-
v-p11 64 206.4 253.7 182.4 4.3
25 p11 to t-p23 were introduced by Chao et al. (1995b) and instances
v-p12 451 236.4 354.7 184.2 1.3
35 t-p24 to t-p34 are taken from Cordeau et al. (1997). Results
v-p13 5 1491.6 127.6 97.2 3.6
12 were given by Christofides and Beasley (1984), Paletta
v-p14 5 9.6 37.4 0 6.90 (1992), Chao et al. (1995b), Cordeau et al. (1997), Paletta
v-p15 1 23.4 93.9 1.2 0 (2002) and Bertazzi et al. (2004). A detailed description of the
v-p16 2 41.4 217.7 9.6 0 instances indicating the number of cities and the planning
v-p17 96 22.2 56.7 12.6 2.
horizon is given in Table 11.
v-p18 40 64.2 142.5 46.8 101.
v-p19 120 135.6 258.3 79.2 1.4
11 As in case of the PVRP also some of the PTSP instances
v-p20 60 232.2 889.1 91.8 9.6
36 are degenerated. More precisely, in instances t-p1, t-p3, t-
v-p21 373 33 72.5 21.6 49.7
2.3 p4, t-p6, t-p7 and t-p9, all customers have a visit frequency
v-p22 528 132 169.6 91.8 158.4
of one. Due to the constraint that at least one customer has
v-p23 42.09 342 341.4 306 310.3
v-p24 114.31 31.8 52.2 18.6 1
v-p25 6 31.2 46.9 5.4 18.
v-p26 9 7.53 31.2 45.2 7.8 38 Table 9
v-p27 2 82.8 6 82.8 6. Best known results for PVRP old data, found with different parameter
v-p28 41 80.4 664.6 3 64. settings
v-p29 3
19 76.2 59.3 29 54. Instance CB TB RG CGW CGL ALP VNS
v-p30 19.712 171 7 7
161.4 7.
76.9
v-p31 7650 160.8 877.1 144.6 v-p1 547.4 537.3 524.6 524.61 531.02
75.1
v-p32 8316 145.8 70.4 6 524.61
v-p2 1443.1 1481.3 1355.4 1322.9 1322.87 1324.74
67.8
9 1322.87 v-p3 54 524. 524. 537. 524.
Total 44914.6 4454.4 4755.6 2099.4 3225.8 v-p4 6.7
84 867 6
840. 61
835. 37
845. 61
835.
v-p5 3.9
218 219 .8
2141 2
2046. 43
2027. 97
2043. 26
2028.
v-p6 7.3
93 2.5 .3 20847. 99836. 75840. 02835.
Table 8 v-p7 8.2
83 833 2
831. 37
826. 1
829. 45
827.
v-p8 9.2
215 228 .6
2108 1
2042. 14
2034. 65
2052. 39
2034.
PVRP new data – Run times in seconds, both algorithms executed on a PC
v-p9 1.3
87 1.8 .3 0 828. 15826. 21829. 15827.
with 3.2 GHz
v- 5
167 183 1638 3
1611. 14
1595. 65
1621. 39
1593.
Instance T T* p10
v- 4 84 3.7
87 .5820 9 785. 84779. 21782. 45779.
CGL VNS CGL VNS p11
v- 7.3 8.5 .3
1312 7
1219. 29
1195. 17
1230. 06
1201.
v-pr01 34.2 180.3 22.8 5 p12
v- 3638 6
3538 88
3511. 95 79
3513.
v-pr02 97.2 283.3 84.6 239. p13
v- .1 954. 62954. 954. 69954.
v-pr03 210.6 278.2 163.2 8.2
25 p14
v- 8
1862. 81
1862. 81
1862. 81
1862.
v-pr04 331.2 314.9 284.4 5.9
29 p15
v- 6
2875. 63
2875. 63
2875. 63
2875.
v-pr05 576.6 264.9 5 5.3
24 p16
v- 2
1614. 24
1597. 24
1597. 24
1597.
v-pr06 1090.2 324.4 1
1031.4 8.9
31 p17
v- 4
3217. 75
3147. 75
3157 75
3136.
v-pr07 94.8 246.7 82.2 8.5
16 p18
v- 7
4846. 24
4834. 4846. 69
4834.
v-pr08 260.4 338.6 2 4.7
29 p19
v- 5
8367. 34
8367. 49
8412. 34
8367.
v-pr09 802.2 423.9 3733.2 7.4
41 p20
v- 4
2216. 4
2184. 02
2173. 4
2170.
v-pr10 2023.8 376.1 2013.6 7.3
37 p21
v- 1
4436. 04
4271. 58
4330. 61
4193.
5.1 p22
v- 4
6769 11
6602. 59
6813. 95
6420.
Total 5521.2 3031.3 5168.4 267 p23 59 45 71
v- 3773 3687. 3702. 3687.
0.7 p24
v- 3826 46
3777. 02
3781. 46
3777.
p25
v- 3834 15
3795. 38
3795. 15
3795.
p26 33 33 32
As also done by CGL, we collected all results obtained v-p27 23401.6 21956.46 22561.33 21956
for different run times and different parameter settings dur- v-p28 23105 22934. 22562. 22305.
ing the fine tuning phase in order to keep track of the best v-p29 .1
24248 71
22909. 44
23752. 34
22639.
solutions found. Our algorithm was able to improve almost v-p30 .2
80982 36
75016. 15
76793. 85
74464.
v-p31 .1
80279 58
78179. 99
77944. 26
76552.
all of the best known results reported in literature as can be
v-p32 .1
83838.7 89
80479.2 79
81055.52 25
78072.88
 Table 10
efficiently. Our code is designed to solve truly periodic
Best known results for PVRP new data, found with different parameter settings
problems.
Instance CGL VNS %
v-pr01 0.
4.2.2. Parameter settings
2209.02 2209.02
v-pr02 3799.28 3774.09 —0 In order to provide an (almost) generic solution
v-pr03 5218.13 5175.15 0.6
— approach, we kept all parameters for the PTSP same as
5914.93 0.8
—
v-pr04 6012.79 for the PVRP. So we again set the initial temperature T
v-pr05 6769.8 6618.95 1.6
—
2.2 for SA to 7 for all instances. The temperature is decreased
v-pr06 8422.64 8258.08 —
v-pr07 4997.41 4996.14 1.9
— linearly every 1000 iterations, in a way that it becomes 0 in
v-pr08 7094.52 6989.81 0.0
— the last iterations. The penalty for an empty day is set to
v-pr09 10370.45 10075.4 1.4
— 1000 and the stopping condition is a fixed number of
2.8
v-pr10 13370.04 12924.66 ——
3.3
iterations.
1.50

4.2.3. Numerical results

Table 12 reports the results of our VNS with indepen-
to be visited every day, the best solution in these cases is to 6 7 8
dent runs of 10 , 10 and 10 iterations compared to the
form a TSP on one day with all the customers except for
best results obtained by the different algorithms of Chao
T — 1 customers, that are close to the depot. Then on each et al. (1995b) (CGW), Cordeau et al. (1997) (CGL), Paletta
of the remaining T — 1 days one of these close customers is (2002) (P) and Bertazzi et al. (2004) (BPS). The average
visited. We report our results also for these degenerated in 6
solution quality is about the same after 10 iterations
instances, but no fine tuning was made to solve these more
Table 12
6 7 8
Results for the PTSP with 10 , 10 and 10 iterations compared to the best
Table 11 solution values of the algorithms of Chao et al. (1995b), Cordeau et al.
Description of the PTSP test instances (1997), Paletta (2002) and Bertazzi et al. (2004) with a recommended value
Instance n t Service frequencies of parameters
6 7 8
 f1 f f f f f6 Instance Min VNS 10 VNS 10 VNS 10
t-p1 50 5 2
2 3 4 5 t-p1 436. 433. — 432. — 432. —
t-p2 50 10 5 2 7 t-p2 50
1106. 18
1105. 0.76
— 1
1106. 1.01
0. 1
1105. 1.0
—
t-p3 70469. 81470. 0.08 0. 84467. —0 81466. 0.0
—
t-p3 50 57 5 6
t-p4 16
554. 77
556. 30. 42
552. 0.37
— 71
549. 0.5
—
t-p4 75 70 2
t-p5 20
1384. 68
1388. 4
0. 39
1384. 0.33
— 05
1384. 0.9
—
t-p5 75 35 5 3 1
t-p6 75643. 35666. 23. 58652. 0.01
1. 05645. 0.00.
t-p6 75
1 70 4 1
t-p7 59
646. 16
661. 52. 65
649. 40. 65
644. —3
t-p7 02
100 105
t-p8 65
1633. 82
1618. —3 17
1615. —3 53
1614. 0.3
—
t-p8 100 5 04 4 1
t-p9 92733. 91747. 0.92
1. 51729. 1.13
— 39723. 1.2
—
t-p9 100 8 100 6 4
t-p10 13
1240. 74
1247. 90. 33
1237. 0.52
— 08
1235. 1.3
—
t-p10 100 5 04 4 1
t-p11 01490. 37490. 50. 72490. 0.18
0. 01490. 0.40.
t-p11 65 4 40 16 4 5
t-p12 97
664. 97
664. 00. 97
664. 00. 97
664. 00.
t-p12 87 4 68 12 7
t-p13 10
830. 1830. 00. 1830. 00. 1830. 00.
t-p13 109 4 84 26 9
t-p14 80
994. 8994. 00. 8994. 00. 8994. 00.
t-p14 131 4 90 20 1
t-p15 60
1157. 6
1157. 00. 6
1157. 00. 6
1157. 00.
t-p15 153 4 116 24 1
t-p16 07660. 09660. 00. 07660. 00. 07660. 00.
t-p16 48 4 23 18 3
t-p17 12
776. 52
778. 00. 12
776. 00. 12
776. 00.
t-p17 66 4 42 26
t-p18 43
873. 82
879. 30. 71
875. 00. 43
874. 00.
t-p18 84 4 54 2
t-p19 70
958. 8964. 70. 82
965. 20. 42
960. 00.
t-p19 102 4 6 38
t-p20 51
1033. 1045.61 61. 54
1035. 70. 69
1035. 0. 2
t-p20 120 4 8 4
t-p21 58
1375. 45
1375. 10. 51
1375. 10. 27
1375. 10.
t-p21 77 4 50 10 7
t-p22 07
4319. 08
4312. —0 07
4312. —0 07
4312. —0
t-p22 154 4 116 24 1
t-p23 72
8390. 33
8384. 0.17
— 31
8349. 0.17
— 31
8315. 0.1
—
t-p23 231 4 2
16 48 24
t-p24 53
2064. 88
2065. 0.07 0. 26
2064. 0.49 0. 45
2064. 0.80.
t-p24 48 4 82 12 1
t-p25 84
3231. 03
3211. —0 84
3208. —0 84
3207. —0
t-p25 96 4 4 2 2
t-p26 50
4084. 93
4041. 0.61
— 49
4045. 0.71
— 88
4033. 0.7
—
t-p26 144 4 78 34 34
t-p27 75
4621. 66
4554. 1.05
— 73
4547. 0.96
— 36
4545. 1.2
—
t-p27 192 4 92 46 46
t-p28 36
4682. 13
4635. 1.45
— 77
4628. 1.59
— 29
4622. 1.6
—
t-p28 240 4 126 68 68
t-p29 54
5595. 26
5542. 1.01
— 24
5529. 1.16
— 16
5527. 1.2
—
t-p29 288 4 0
14 70 70
t-p30 45
4453. 09
4443. 0.95
— 68
4436. 1.18
— 39
4435. 1.2
—
t-p30 72 6 41 12 1 2 1
t-p31 15
5405. 23
5375. 0.22
— 31
5370. 0.38
— 39
5368. 0.4
—
t-p31 144 6 38 38 38 38
t-p32 40
7346. 4
7259. 0.56
— 59
7244. 0.64
— 45
7238. 0.6
—
t-p32 216 6 56 56 56 56
t-p33 32
8394. 14
8242. 1.19
— 02
8216. 1.39
— 45
8208. 1.4
—
t-p33 288 6 74 74 74 74
52 78 1.81 48 2.12 38 2.2
2 2 2 2 Average 0.05 —0.34 —0.52
The notation is the same as for the PVRP instances.
 V.C. Hemmelmayr et al. / European Journal of Operational Research 195 (2009) 791–802 801

Table 13 We refrain from reporting detailed run times here since

Best known solutions for PTSP instances, found with different parameter most algorithms did not solve all instances. But the run
settings 6
times for 10 iterations of the VNS are comparable or
Instance BKS PTSP BKS VNS %
shorter than those of the other metaheuristic approaches
t-p1 43 432.10 0. also after adjustment of run times w.r.t. the machine used.
t-p2 2.
1105.81 1105.81 0
0.
t-p03 466 466.71 0
0.
t-p4 .71
549 549.05 0
0.
4.3. Analysis of neighborhood structure
t-p5 .05
1382.33 1382.33 0
0.
t-p6 64 643.50 0
0. We also investigated the contribution of all shaking
t-p7 3.
64 643.80 0
0. operators to the performance of the algorithm and the best
t-p8 3.
1613.42 1611.96 —0 ordering of the operators. Table 14 reports the number of
t-p9 721 720.72 0.0
— times that a neighborhood structure lead to a new incum-
t-p10 .24
1237.77 1233.53 0.0
— 6
0.30. bent solution after local search within 10 iterations. Three
t-p11 490 490.97
.97 00. different orderings of the neighborhood structures are ana-
t-p12 66 664.10
t-p13 4.
83 830.80 00. lyzed: move-cross-change combination (cc), cross-cc-move
t-p14 0.
99 994.60 0
0. and cc-move-cross. From the average deviation to the best
t-p15 4.
1157.07 1157.07 0
0. solution reported in literature it can be seen that the cc-
t-p16 660 660.12 0
0. move-cross ordering delivers the best results. This can be
t-p17 .12
776 776.43 0
0. explained as follows. As VNS starts from the first neigh-
t-p18 .43
873 873.73 0
0. borhood again when it finds a new incumbent solution
t-p19 .73
958 958.51 00. (see Fig. 1), earlier neighborhood structures are used more
t-p20 .51
1033.58 1033.58 0
0. often, especially in the beginning of the search. The order-
t-p21 1375.07 1375.07 0
0.
0
ing cc-move-cross leads to the best results, because it is
t-p22 4312.31 4312.31 0.
t-p23 8308.51 8308.48 00. important that in the beginning the most suitable visit
t-pr01 2064.84 2064.84 0
0. day combination is selected. However, according to exten-
t-pr02 3207.44 3205.94 —0 sive tests the later neighborhoods, move and cross, also
t-pr03 4030.54 4027.71 0.0
— play an essential part for obtaining good results. These
t-pr04 4558.94 4538.19 0.0
— neighborhood structures become more important in the
t-pr05 4628.89 4613.58 0.4
— later phases of the optimization runs – when already good
t-pr06 5534.94 5521.24 0.3
— visit day combinations are assigned to the customers. Table
t-pr07 4435.39 4435.39 0.20.
14 shows that the number of new incumbent solutions
t-pr08 5376.11 5366.53 —0
t-pr09 7282.39 7234.35 0.1
—
found is generally decreasing with the index of the neigh-
t-pr10 8280.07 8199.55 0.6
— borhood if the ordering of the neighborhoods is correct.
0.9 If however the cc neighborhoods are scheduled after move
7
and/or cross they can show more improvements than ear-
compared to the best results found by the other authors.
7 lier neighborhoods. This is an indication that the ordering
Applying the VNS for 10 iterations the solution quality is not ideal.
can be improved by 0.34%. Note that this improvement
is over the best solution obtained by any of the approaches
mentioned. An individual comparison of our VNS with Table 14
any of the other approaches would yield even higher Use of different neighborhoods with different ordering of neighborhood
improvements. For longer run times further improvements structures
8 Nk move-cross-cc cross-cc-move cc-move-
are obtained. After 10 iterations the solution quality is
cross
0.52% better compared to the best results reported by the 1 m 1551 cr 1688 cc 1845
other authors. As mentioned in Section 4.2.1 there are 2 m 942 cr 936 cc 1249
some instances, that are no real PTSP instances but rather 3 m 622 cr 532 cc 820
standard TSP instances. If these degenerated instances, 4 cr 421 cr 311 cc 539
namely t-p1, t-p3, t-p4, t-p6, t-p7, and t-p9 are disregarded, 5 cr 276 cr 195 cc 363
6 cr 174 cr 121 cc 260
and the comparison is only made for the remaining ones,
7 cr 113 cc 703 m 155
then we already have an average improvement by our 8 cr 77 cc 506 m 107
6
VNS of —0.24% after 10 iterations and improvements of 9 cr 54 cc 345 m 71
7 8
—0.40% and —0.49% after 10 and 10 iterations, 10 cc 523 cc 254 cr 54
respectively. 11 cc 396 cc182 cr 36
Table 13 shows the best known results produced by our 12 cc 276 cc 132 cr 25
13 cc 198 m 20 cr 17
algorithm. Best known solutions means all solutions that
14 cc 138 m 17 cr 15
were found with different parameter settings. The results 15 cc 115 m 11 cr 9
are compared to the best known results given in the litera- Avg. dev. to 3.64 0.98 0.05
ture. New best results were found for 11 instances and for best solution
all the other instances tie were achieved.
 Table 15
strength of our algorithm is that on average it provides bet-
Use of different neighborhoods with cc-move-cross for PVRP and PTSP
ter results than the existing techniques especially when the problem
Nk PVRP PTSP
size increases. Another important aspect with
10
6
10
7
10
8
10
6
10
7
10
8
respect to runtime is that the algorithm scales quite well.
1 3 338 367 184 173 1915 The increase in runtime is much lower when the problem
2 10 9
216 3
244 5
124 4
122 1281 size increases compared to the other algorithms in the liter-
3 19 7
162 2
174 9 82 1 80 833 ature. For the PTSP our VNS finds 11 new best results on
41 3119 6133 053 152 528
4
368 the existing instances in the literature. Moreover also the
5 083 594 4
103 9
36 0
36 average results resemble or even outperform the results
6 5
66 777 9 84 326 125 259
7 5
39 0
45 7
51 0
15 7
17 177 published in the literature within a comparable runtime.
8 5
29 633 238 510 012 123 For future research we will extend our algorithm to be
9 6
25 7
29 9
33 77 27 81 applicable to Inventory Routing Problems. This algorithm
9 636 645 15 95 60
10 31
41 will then be applied to the periodic delivery of blood prod-
11 4
23 5
26 1
32 34 46 ucts to hospitals (see Hemmelmayr et al., 2006).
12 5
18 021 026 62 2 28
13 9
16 9
17 1
22 15 18 18
14 7
14 816 720 71 61 15 Acknowledgement
15 7
14 8
16 4
18 59 13 6
5 0 8 0
Financial support from the Oesterreichische National-
bank (OeNB) by grant #11187 and from the Austrian Sci-
While we report detailed tests results for different order- ence Fund (FWF) by grant #L286-N04 is gratefully
ings of the neighborhoods only for the PTSP (similar
acknowledged. We are grateful to the authors of CGL
results are obtained for the PVRP, where some limited test
were performed), we investigated for PVRP and PTSP, for providing the code to us.
whether all 15 neighborhoods are in fact reasonable, i.e.,
whether all contribute to the performance of the algorithm. References
Table 15 shows (for the best ordering of the neighborhood
structures, i.e., cc-move-cross), the usage of the neighbor- Alegre, J., Laguna, M., Pacheco, J., 2007. Optimizing the periodic pick-up
hood structures for a different number of iterations. The of raw materials for a manufacturer of auto parts. European Journal
reported results are the average over 10 runs, summed up of Operational Research 179, 736–746.
Angelelli, E., Speranza, M.G., 2002. The periodic vehicle routing problem
over all instances. It can be seen that most incumbent solu-
with intermediate facilities. European Journal of Operational Research
tions are found in an early phase of the optimization runs. 137 (2), 233–247.
Move and cross neighborhoods tend to find more improve- Beltrami, E.J., Bodin, L.D., 1974. Networks and vehicle routing for
ments also in the later stages of the optimization. municipal waste collection. Networks 4, 290–306.
Bertazzi, L., Paletta, G., Speranza, M.G., 2004. An improved heuristic for
the period traveling salesman problem. Computers & Operations
5. Conclusion Research 31, 1215–1222.
Blakeley, F., Bozkaya, B., Cao, B., Hall, W., Knolmajer, J., 2003.
We presented an (almost) generic Variable Neighbor- Optimizing periodic maintenance operations for Schindler elevator
hood Search heuristic for the Periodic Vehicle Routing corporation. Interfaces 33 (1), 67–79.
Braysy, O., 2003. A reactive variable neighborhood search for the vehicle-
Problem and the Periodic Traveling Salesman Problem that
routing problem with time windows. INFORMS Journal on Comput-
is competitive or even outperforms the existing methods. ing 15, 347–368.
The main features of this algorithm are a simple and flexible local Carter, M.W., Farvolden, J.M., Laporte, G., Xu, J., 1996. Solving an
search as well as an acceptance criterion for neighbor- ing integrated logistics problem arising in grocery distribution. INFOR 34, 290–
solutions inspired by Simulated Annealing. We show the 306.
robustness of our approach by applying the identical basic Chao, I.M., Golden, B.L., Wasil, E.A., 1995a. An improved heuristic for
algorithm to both problem classes. The only difference between the period vehicle routing problem. Networks 26, 25–44.
the implementations for PVRP and PTSP is the choice of Chao, I.M., Golden, B.L., Wasil, E.A., 1995b. A new heuristic for the
period traveling salesman problem. Computers & Operations Research
the local search. We also made no special efforts to 22, 553–565.
adapt the algorithm to degenerated problem instances (i.e. Christofides, N., Beasley, J.E., 1984. The period routing problem.
those without any periodic aspects). Nevertheless, the results Networks 14, 237–256.
obtained through an extensive numerical analysis showed that Clarke, G., Wright, J.W., 1964. Scheduling of vehicles from a central
the algorithm is competitive to other state of the art depot to a number of delivery points. Operations Research 12, 568–
approaches applied to these problem classes. Con- sidering the 581.
Cordeau, J.F., Gendreau, M., Laporte, G., 1997. A tabu search heuristic
best solutions found our algorithm for the PVRP for periodic and multi-depot vehicle routing problems. Networks 30, 105–
outperforms the existing techniques by finding 24 new best 119.
solutions and 13 ties. In detail, we improved 9 test instances of Croes, G.A., 1958. A method for solving traveling salesman problems.
the 10 instances of the new data and 15 test instances out Operations Research 6, 791–812.
of the 32 instances of the old data. The Dongarra JJ. Performance of Various Computers Using Standard Linear
Equations Software, Technical Report of the Computer Science
Department, University of Tennessee, Knoxville, 2005.
Dueck, G., 1993. New optimization heuristics: The great deluge algorithm Local Search in Combinatorial Optimization. Wiley, Chichester, pp. 337–
and the record-to-record travel. Journal of Computational Physics 360.
104, 86–92. Kirkpatrick, S., Gelatt, C.D., Vecchi, M.P., 1983. Optimization by
Drummond, L.M.A., Ochi, L.S., Vianna, D.S., 2001. An asynchronous Simulated Annealing. Science 220 (4598), 671–680.
parallel metaheuristic for the period vehicle routing problem. Future Lin, S., 1965. Computer Solutions of the Traveling Salesman Problem.
Generation Computer Systems 17, 379–386.
Bell Systems Technical Journal 44, 2245–2269.
Eilon S, Gandy W, Christofides N. Distribution Management: Mathematical
Mladenovic´, N. A variable neighbourhood algorithm – a new metaheu-
Modeling and Practical Analysis. Griffin: London; 1971.
ristic for combinatorial optimization. Abstract of Papers presented at
Gaudioso, M., Paletta, G., 1992. A heuristic for the periodic vehicle
Optimization Days 1995, 112.
routing problem. Transportation Science 26 (2), 86–92. Mladenovic´, N., Hansen, P., 1997. Variable Neighborhood Search.
Gendreau, M., Hertz, A., Laporte, G., 1992. New insertion and
Computers & Operations Research 24, 1097–1100.
postoptimization procedures for the traveling salesman problem.
Paletta, G., 1992. A multiperiod traveling salesman problem: heuristic
Operations Research 40, 1086–1094.
algorithms. Computers & Operations Research 19, 789–795.
Golden, B.L., Wasil, E.A., 1987. Computerized vehicle routing in the soft
Paletta, G., 2002. The period traveling salesman problem: a new heuristic
drink industry. Operations Research 35, 6–17.
algorithms. Computers & Operations Research 29, 1343–1352.
Hadjiconstantinou, E., Baldacci, R., 1998. A multi-depot period vehicle
Polacek, M., Hartl, R.F., Doerner, K.F., Reimann, M., 2004. A Variable
routing problem arising in the utilities sector. Journal of the Opera-
tional Research Society 49 (12), 1239–1248. Neighborhood Search for the Multi Depot Vehicle Routing Problem
Hansen, P., Mladenovic´, N., 1997. Variable Neighborhood Search for the p- with Time Windows. Journal of Heuristics 10, 613–627.
median. Location Science 5, 207–226. Polacek, M., Doerner, K.F., Hartl, R.F., Kiechle, G., Reimann, M., 2007.
Hansen, P., Mladenovic´, N., 2000. Variable Neighborhood Search. In: Scheduling periodic customer visits for a Travelling Salesperson.
Pardalos, P.M., Resende, M.G.C. (Eds.), Handbook of Applied European Journal of Operational Research 179, 823–837.
Optimization. Oxford University Press, New York, pp. 221–234. Russel, R.A., Gribbin, D., 1991. A multiphase approach to the period
Hansen, P., Mladenovic´, N., 2001. Variable Neighborhood Search: routing problem. Networks 21, 747–765.
Principles and Applications. European Journal of Operational Russel, R.A., Igo, W., 1979. An assignment routing problem. Networks 9, 1–
Research 130, 449–467. 17.
Hemmelmayr V., Doerner K.F., Hartl R.F., Savelsbergh M. Delivery Tan, C.C.R., Beasley, J.E., 1984. A heuristic algorithm for the periodic
Strategies for Periodic Blood Supplies. Submitted, 2006. vehicle routing problem. Omega 12 (5), 497–504.
Kindervater, G.A.P., Savelsbergh, M.W.P., 1997. Vehicle routing: Han- Van Breedam A. An Analysis of the behaviour of heuristics for the vehicle
dling edges exchanges windows. In: Aarts, E., Lenstra, J.K. (Eds.), routing problem for a selection of problems with vehicle-related,
customer-related and time-related constraints. Ph.D. dissertation,
University of Antwerp, 1994.