Constructive heuristics for the multi-compartment vehicle routing
problem with stochastic demands
Jorge E. Mendozaa Bruno Castanierb Christelle Gu¶eretb Andr¶es L. Medagliac
Nubia Velascoc
a) ¶Equipe Optimisation des Systµemes de Production et Logistiques, LISA (EA CNRS 4094), Universit¶e Catholique de l'Ouest, 3 Place
Andr¶e Leroy 49008 Angers, France
b) ¶Equipe Systµemes Logistiques et de Production, IRCCyN (UMR CNRS 6597), ¶Ecole des Mines de Nantes, B.P. 20722 F-44307 Nantes
Cedex 3, France
c) Centro para la Optimizaci¶on y Probabilidad Aplicada (COPA), Departamento de Ingenier¶³a Industrial, Universidad de los Andes, Cr 1
Este No. 19A-10 Bogot¶a, Colombia
Abstract
The vehicle routing problem with stochastic demands (VRPSD) consists of designing transportation
routes of minimal expected cost to satisfy a set of customers with random demands of known probability
distribution. This paper tackles a generalization of the VRPSD known as the multi-compartment VRPSD
(MC-VRPSD), a problem in which each customer demands several products that, because of incompatibility
constraints, must be loaded in independent vehicle compartments. To solve the problem we propose three
simple and e®ective constructive heuristics based on a stochastic programming with recourse formulation.
One of the heuristics is an extension to the multi-compartment scenario of a savings-based algorithm for
the VRPSD; the other two are di®erent versions of a novel look-ahead heuristic that follows a route-¯rst,
cluster-second approach. In addition, to enhance the performance of the heuristics these are coupled with
a post optimization procedure based on the classical 2-Opt heuristic. The three algorithms were tested
on instances of up to 200 customers from the MC-VRPSD and VRPSD literature. The proposed heuristics
unveiled 26 and 12 new best known solutions for a set of 180 MC-VRPSD problems and a 40-instance testbed
for the VRPSD, respectively.
1 Introduction
Vehicle routing problems (VRPs) consist of designing a set of minimal-cost vehicle routes to serve the demands
for goods or services of a group of geographically spread customers, satisfying a group of operational constraints.
1

Toth and Vigo (2002), Cordeau et al. (2006), Golden et al. (2008), and Laporte (2009) present comprehensive
surveys on VRP variants and solution methods. Classical VRPs normally assume perfect knowledge of problem
parameters such as customer locations and demands. However, in practice, the nature of these parameters is
often stochastic, meaning that they are a®ected by di®erent sources of uncertainty that should be considered
during the route planning process. Consequently, in recent years a number of research e®orts have been ad-
dressed to solve VRP variants, known as stochastic VRPs (SVRPs), that incorporate uncertainty in the problem
parameters. Among SVRPs, probably the most studied variant is the VRP with stochastic demands (VRPSD),
a problem where customer demands are modeled as random variables. For a complete overview of SVRPs the
reader is referred to Gendreau et al. (1996a).
The main e®ect of the stochastic demands is that if their realizations (actual values) turn out to be larger
than expected, the vehicle capacity may be exceeded leading to what is commonly known as a route failure.
Di®erent solution frameworks can be applied to deal with uncertain demands in the VRPSD. A taxonomy
of these frameworks classi¯es them into dynamic or static (Secomandi and Margot, 2009). Routing decisions
in dynamic approaches are made in multiple stages and based on demand realizations (at former stages).
Depending on the remaining vehicle capacity and the set of unvisited customers, routes are re-optimized at each
stage. On the other hand, routing decisions in static approaches once made remain unchanged regardless of the
demand realizations. While dynamic strategies produce more accurate solutions, static strategies are preferred
from the computational point of view since the problem, known to be computationally intractable, is solved
only once (Laporte et al., 2002). Static strategies are also particularly useful when a stable solution is sought
(Bertsimas, 1992) or when re-optimization during route execution is impossible due to the lack of information.
Most of the research in the VRPSD has been based on static strategies.
Among the static approaches, the most widely explored in the VRPSD literature is the stochastic pro-
gramming with recourse (SPR) framework. SPR solves the problem in two phases: in the ¯rst phase a set
of routes is planned, while in the second phase the planned routes are executed. In case of a route failure, a
problem-dependent prede¯ned corrective action, known as recourse, is taken to recover the solution feasibility.
In general, the route failures generate an extra cost only known after the second phase. Thus, the objective of
the problem is to design during the ¯rst phase a set of routes that minimizes the sum of the cost of the planned
routes and the expected cost of the route failures.
Exact methods for the VRPSD based on SPR include that of Laporte et al. (2002) who proposed a branch-
and-bound method based on the integer L-shaped algorithm for a formulation where the recourse action is a
round trip to the depot that restores the vehicle capacity. Later, Christiansen and Lysgaard (2007) proposed a
branch-and-price approach that proved to be competitive in tight instances of the problem, that is, when the
ratio between the sum of the expected demands of all customers and the maximum number of available vehicles
2

is close to the vehicle capacity. More recently, Rei et al. (2007) tackled the single-route VRPSD (SVRPSD).
Their method consists of using local branching (Fischetti and Lodi, 2003) to generate optimality cuts on an
integer L-shaped algorithm. Although successful, these approaches are limited to solve instances of up to 100
customer nodes. Consequently, industrial-scale problems are out of reach and should be targeted using heuristics
(Gendreau et al., 1996a).
The pioneering SPR-based heuristic for the VRPSD is the savings algorithm of Dror and Trudeau (1986);
an extension of the Clarke and Wright heuristic (CW) for a formulation where the recourse action is to perform
upon failure an emergency trip from the depot to every customer remaining on the planned route. Later,
Savelsbergh and Goetschalckx (1995) developed a two-phase approach to solve an SPR formulation that allows
only one failure per route. In the ¯rst phase, their method extends the generalized assignment heuristic of Fisher
and Jaikumar (1981) to solve the deterministic counterpart of the problem using the expected demand for each
customer. This deterministic solution is improved in the second stage using post-optimization procedures
specially tailored to handle the stochastic nature of the demands. More recently, Yang et al. (2000) tackled
a variant of the problem with restocking. The principle behind the restocking policy is to decide upon the
remaining capacity after completing the service of each customer and perform preventive trips to the depot to
avoid future route failures. These authors introduced two constructive heuristics based on insertion procedures
and a dynamic programming algorithm for tour splitting. Lately, Liu and Lai (2009) introduced a sweep
algorithm for an SPR formulation where the recourse action considers cooperation between routes and the
demand uncertainty is modeled using fuzzy sets. Some authors have addressed the VRPSD and its variants
using metaheuristic approaches. Gendreau et al. (1996b) proposed tabustoch, a tabu search (TS) algorithm
designed to tackle a variant of the problem where in addition to the demands customers are also stochastic,
that is, they are present (or not) with a given probability. Ak and Erera (2007) tackled a stochastic routing
problem that allows cooperative strategies using TS. Bianchi et al. (2006) proposed a set of algorithms for
the SVRPSD with restocking comprising simulated annealing (SA), iterated local search (ILS), ant colony
optimization (ACO), evolutionary algorithms (EA), and TS.
This paper tackles a generalization of the VRPSD known as the Multi-Compartment VRPSD (MC-VRPSD),
a problem where each customer has a stochastic demand for di®erent products that due to incompatibility
constraints must be transported in independent vehicle compartments. The MC-VRPSD naturally arises in
several practical situations. For instance, dairies often use vehicles with multiple compartments to collect milk
of di®erent types (e.g., from cows and goats) and qualities (e.g., di®erent suckling dates); petroleum companies
deliver di®erent types of fuel to outlet retailers using multi-compartment tankers; public utilities use trucks
with compartments to perform selective waste collection; and food companies distribute in compartmentalized
vehicles groceries that require di®erent levels of refrigeration.
3

Despite the numerous practical applications, research in vehicle routing problems with multiple compart-
ments is rather scarce and much of the e®ort has concentrated on deterministic demands. Research along this
line has mainly focused on two problem categories based on two di®erent service policies commonly found in
practical settings. Under the ¯rst policy, the demand for each product can be served by di®erent routes as
long it is fully served in a single visit. Some of the solution approaches for problems that fall into this cate-
gory include the TS and memetic algorithm (MA) of El Fallahi et al. (2008), the guided local search (GLS)
algorithms of Muyldermans and Pang (2010a; 2010b), the heuristic components of Derigs et al. (2010) and the
multi-start local search procedures recently proposed by Caramia and Guerriero (2010). On the other hand,
under the second policy, the demand for every product must be served by the same route (i.e., service splitting
is not allowed). Early work on problems falling under the latter category include the constructive heuristic
of Van der Bruggen and Salomon (1995), followed by the more recent Lagrangean heuristic of Chajakis and
Guignard (2003), the implicit enumeration algorithm of Ruiz et al. (2004) and the TS and MA of El Fallahi
et al. (2008). Likewise, the few attempts that have been made to solve the MC-VRPSD also fall under the
latter service policy. Tatarakis and Minis (2009) tackled the problem of selecting the optimal restocking plan
in a single-route MC-VRPSD in which the sequence of customers (route) is ¯xed beforehand. The authors
propose a set of dynamic programming algorithms and solve to optimality problems of up to 15 customers.
Mendoza et al. (2010) proposed two MAs based on an SPR formulation of the problem. The ¯rst algorithm
(MA/SCS) solves the deterministic counterpart of the problem using the expected demands, while trying to
avoid route failures by preserving on each compartment some spare capacity in case of high demand realizations.
On the other hand, the second memetic algorithm (simply labeled MA) incorporates explicitly the stochastic
nature of the demands in the route planning process. Their computational experiments showed that MA out-
performs MA/SCS in terms of solution quality (accuracy), while MA/SCS proved to be more e±cient in terms
of computational performance (speed).
To expand the existing toolbox for solving the MC-VRPSD, under the single-vehicle service scenario, this
research introduces three new constructive heuristics that focus not only on accuracy and speed, but also on
simplicity and °exibility. The latter two design issues are commonly neglected in the academic literature,
but often hinder the VRP technology transfer to commercial software and real-world applications (Cordeau
et al., 2002). One of the proposed heuristics is an extension to the multi-compartment scenario of the savings-
based approach for the VRPSD of Dror and Trudeau (1986), while the other two algorithms are variants of a
novel look-ahead heuristic. The look-ahead heuristic integrates the pilot method (Vo¼ et al., 2005) with a route-
¯rst, cluster-second procedure. The routing phase uses simple constructive heuristics for the traveling salesman
problem (TSP) to produce a giant tour with a sequence of customer visits. Then, the clustering phase, splits
the giant tour into feasible routes for the MC-VRPSD using the dynamic programming principle. Finally, to
4

enhance the performance of the three constructive heuristics we propose a simple post-optimization procedure
based on the classical 2-Opt improvement algorithm.
The remainder of the paper is organized as follows. Section 2 formally states the problem and formulates it
as an SPR program. Section 3 introduces an expression to e±ciently approximate the objective function in the
SPR formulation. Sections 4 and 5 outline the proposed heuristics and Section 6 describes the post-optimization
procedure. Section 7 presents the computational experiments conducted on a set of benchmark instances and
compares the heuristics against a competing memetic algorithm proposed in the literature for the MC-VRPSD.
Also in Section 7 we compare the heuristics with previously published results for the VRPSD. Finally, Section
8 concludes the paper and outlines future research perspectives.
2 Problem formulation
Formally, the MC-VRPSD can be de¯ned on a complete and undirected graph G = (V; E) where V = f0; : : : ; ng
is the vertex set and E the edge set. Vertices v = 1; : : : ; n represent the customers and vertex v = 0 represents
the depot. A distance de is associated to edge e = (u; v) = (v; u) 2 E and it represents the travel cost between
vertices u and v. There exists a set P = f1; : : : ; p; : : : ;mg of products that must be transported in independent
compartments of ¯xed capacity Qp. All vehicles are identical and the °eet size is unlimited. For product p
customer v has an independent random demand »v;p following a known probability distribution with mean
¹v;p and standard deviation ¾v;p. The demand realizations ~» are nonnegative and less than the capacity of
the corresponding compartment Qp, yet only known upon the vehicle's arrival to the customer location. Each
customer must be visited only once by exactly one vehicle (route) and the total length of each route cannot
exceed a maximum distance L. Henceforth, without loss of generality the discussion is restricted to the case of
collection routes, nonetheless the case of delivery routes is equivalent.
The MC-VRPSD formulation as a two-stage stochastic programming model follows. In the ¯rst stage, a setR
of a priori or planned routes is designed. Each route r 2 R is a sequence of vertices r = (0; v1; : : : ; vi; : : : ; vnr ; 0),
where vi 2 V n f0g and nr is the number of customers serviced by the route. Depending on the context, we
refer to the route r 2 R also as an ordered set of arcs r = f(0; v1); : : : ; (vi¡1; vi); : : : ; (vnr ; 0)g. In the second
stage, each planned route is executed until a route failure occurs, that is, whenever the capacity of at least one
compartment is exceeded. Upon failure, the compartment is loaded up to its capacity and the recourse action
takes place. The recourse action is de¯ned as a return trip to the depot to unload the compartments, followed
by a trip back to the customer location to complete the service. After service completion, the route is resumed
from that point on as originally planned. It is assumed that there is no practical way to keep track of the
capacity available on each compartment, thus preventive trips to the depot (restocking) are not allowed. The
5

second stage solution is then the actual set of routes traveled by the vehicles. The problem is to determine in
the ¯rst stage the set of planned routes R that minimizes the expected cost E [C] of the second stage solution
given by:
E [C] =
X
r2R
E
h
lr +Gr
³~»
´i
(1)
=
X
r2R
lr +
X
r2R
E
h
Gr
³~»
´i
(2)
where lr denotes the planned length (planned cost) and E
h
Gr
³~»
´i
the expected length of the returning trips
to the depot, or cost of recourse, caused by route failures for each route r 2 R. Since the total travel distance of
each route is a random variable which value is only known when the vehicle returns to the depot after completing
the route, we model the distance constraint as
lr + E
h
Gr
³~»
´i
· L 8r 2 R (3)
3 Estimating the total expected cost of a MC-VRPSD solution
In contrast to deterministic VRPs, the overall cost expression in (2) is not straightforward to calculate. While
the planned cost of route r is easily computed by the sum of the distance covered by the arcs traversed in the
route (lr =
P
(u;v)2r d(u;v)), the computation of the cost of recourse E
h
Gr
³~»
´i
is far more complex. The cost
of recourse of route r is given by
E
h
Gr
³~»
´i
=
nrX
i=1
2£ dvi;0 £ Pr (vi) (4)
where Pr(vi) is the failure probability while servicing each customer vi 2 r.
In the multi-compartment scenario and under the selected recourse action, the failure probability while
servicing a given customer vi, depends on the location (customer) where the last failure occurred and the
demands collected for each product since that last failure. As demonstrated in Appendix A the value of Pr(vi)
is given by:
Pr(vi) =
i¡1X
j=0
2
4Y
p2P
Pr
0
@»0vj ;p +
i¡1X
k=j+1
»vk;p · Qp
1
A¡
Y
p2P
Pr
0
@»0vj ;p +
iX
k=j+1
»vk;p · Qp
1
A
3
5£ Pr(vj) (5)
where vj is the customer being serviced when the last failure occurred and »0vj ;p is customer vj 's unserviced
demand for product p due to the failure. Consequently, the ¯rst product denotes the probability that the
6

collected demand after the failure at customer vj and up to customer vi¡1 does not exceed the capacity of the
corresponding compartment Qp. Similarly, the second product denotes the probability that the capacities are
not exceeded after servicing customer vi.
The calculation of Pr(vi) in (5) requires the evaluation of the probability functions of all unserviced demands
(i.e., »0vj ;p; 8j < i; p 2 P), which may be impossible to perform analytically depending on the probability function
of the demands. Alternatively, one could use estimations calculated by simulation as suggested by Tan et al.
(2007); SÄorensen and Sevaux (2009); and Mendoza et al. (2009) for related routing problems with stochastic
demands. Although e®ective, the simulation approach may be overly expensive, from a computational point of
view, to be embedded within fast heuristics like the ones addressed in this work. To overcome this di±culty, we
use a fast approximation that can be analytically computed for a number of demand probability distributions.
The approximation is based on a take-all policy (TAP) which assumes that even if a route failure occurs at
customer vj the whole demand for all products is serviced before performing the return trip, that is, making
»0vj ;p = 0;8p 2 P in (5). Nonetheless, after unloading the compartments at the depot, the vehicle must return
to the customer being served when the failure occurred.
Extensive computational experiments show that if the capacity of the compartments is su±ciently large
compared to the largest expected demand for the corresponding product, which is the case in many real-world
applications, the true cost of recourse in (4) can be estimated within an average error of 3.00% by the TAP
approximation. This small gap suggests that the approximation is an acceptable proxy of the overall cost and
can be used to guide the search in the proposed heuristics. For a thorough discussion on the computational
experiments conducted to validate the TAP approximation the reader is referred to Mendoza et al. (2008).
4 Stochastic Clarke and Wright heuristic
The Stochastic Clarke and Wright (SCW) heuristic extends the savings-based algorithm introduced by Dror
and Trudeau (1986) for the VRPSD in several ways, namely, SCW: (i) works in a multi-compartment scenario,
(ii) considers the round trip to the depot recourse action, and (iii) deals with the distance constraint of the
routes. Similarly to the classical savings algorithm, the SCW heuristic starts from a trivial solution comprised
of n round trips from the depot to each customer. Then it tries at each iteration to merge two routes by
their extreme vertices (those connected to the depot) aiming to generate the largest possible savings in the
overall cost. In the SPR formulation introduced in x2, the savings generated by the merge of two routes having
customers v and u as extreme points is given by sv;u = E [Cv] +E [Cu]¡E [Cu;v], where E [Cv] and E [Cu] are
the expected costs of the routes servicing customers v and u, and E[Cu;v] is the expected cost of the merged
route that visits customer u before customer v. Because the failure probability at customer v depends on the
7

demands collected at the preceding customers (see (5)), and as it is often the case that demands collected from
the depot to a given customer v di®er from those collected from customer v to the depot, the failure probability
at customer v may di®er upon the traveling direction, leading to di®erent costs of recourse. Consequently, when
calculating the savings generated by a merger both traveling directions should be considered as follows:
sv;u = E [Cv] + E [Cu]¡minfE [Cu;v] ; E [Cv;u]g (6)
Algorithm 1 outlines the logic of our SCW heuristic. Once the single-node routes are generated (steps 2{5),
the merging process begins. At each iteration, the savings of all feasible mergers leading to routes that satisfy
the distance constraint are computed and stored in the savings list S (steps 7{17). Then, the largest saving is
selected from the list, and the corresponding merger is committed (steps 18{29). The process is repeated until
the largest feasible saving has a negative value or no more feasible savings are available (steps 19{21). If the
largest saving has a zero value, the corresponding merger is undertaken to save one vehicle, although the value
of the overall cost does not improve. In our implementation, the savings in (6) are computed using the TAP
approximation in x3.
INSERT ALGORITHM 1
5 Look-ahead heuristic
The next method is a look-ahead heuristic that embeds a route-¯rst, cluster-second (RC) approach. The RC
approach for VRPs builds a single giant tour starting and ending at the depot after visiting every customer
(routing phase), then it splits the tour into routes that satisfy the particular constraints of the VRP on hand
(clustering phase). In the proposed heuristic, the routing phase uses di®erent construction algorithms for the
deterministic TSP (routing procedures), while the clustering phase applies a dynamic programming algorithm
specially designed to tackle the stochastic nature of demands in the multi-compartment scenario (clustering
procedure). One distinctive feature of our heuristic is that the routing phase is guided by a look-ahead mech-
anism that avoids traps often posed by the greedy choices of the TSP algorithms by evaluating the impact of
local moves in the ¯nal solution of the problem. The remainder of this section outlines the building blocks of
the RC approach, namely, the routing and clustering procedures, and explains how they are coupled with the
look-ahead mechanism to generate the proposed heuristic.
5.1 Routing procedures
Rei et al. (2010) point out that in stochastic vehicle routing problems the planned cost of a solution normally
dominates the cost of recourse. For instance, in the solutions by Mendoza et al. (2010) for a set of 40 MC-
8

VRPSD instances the planned cost represents 90% of the total cost of the solutions (on average). To build good
sequences of customers that translate into good planned routes during the clustering phase, our RC approach
uses simple but e®ective constructive heuristics for the geometric TSP, namely, the nearest neighbor (NN) and
the nearest insertion (NI) heuristics (Bentley, 1992). Starting by linking the depot to its closest customer, at
each iteration the NN heuristic grows a directed path connecting the last added node to the closest unvisited
customer, thus greedily minimizing the cost of adding one more customer to the path. Once every customer is
included in the path, the last added vertex is connected to the depot generating the TSP tour. On the other
hand, the NI heuristic starts directly from a partial TSP tour beginning and ending at the depot after visiting
its farthest customer. The tour grows at each iteration by adding the customer whose insertion generates the
smallest increment on the cost of the tour. The insertion procedure is repeated until all customers are visited.
We chose NN and NI over the large set of heuristics available for the TSP based on their tradeo® between
implementation simplicity, accuracy, and computational performance. For an extensive study on TSP heuristics
the reader is referred to Bentley (1992).
5.2 Clustering procedure
The TSP tours generated in the routing phase are not necessarily good nor even feasible solutions for the MC-
VRPSD. On the one hand, a single tour visiting every customer will face a large number of route failures, thus
leading to a large cost of recourse; on the other hand, a single tour will very likely violate the distance constraint.
In our RC approach, the giant tours generated in the routing phase are optimally partitioned into feasible
solutions for the MC-VRPSD using the stochastic split (s-split) procedure (Mendoza et al., 2010) as follows.
From the giant tour, an auxiliary graph G0 is built to ¯nd the optimal partition of the tour into feasible routes.
The directed and acyclic graph G0 = (V 0;A) is composed of the ordered vertex set V 0 = f0; v1; : : : ; vi; : : : ; vng
and the arc set A. Vertex 0 is an auxiliary vertex and vertices v1; : : : ; vn 2 V 0 nf0g in graph G0 are arranged so
that they follow the same order of the permutation (tour) in the routing phase. Arc (vi; vi+nr ) 2 A represents
a feasible route r with expected cost E[Cr] starting and ending at the depot and traversing the sequence of
customers from vi+1 to vi+nr . The s-split procedure ¯nds the set of arcs (i.e., routes) along the shortest path
connecting 0 and vn in G0. Figure 1 illustrates the s-split procedure where a single tour of customers is split into
two feasible routes (L = 70) of minimum expected cost. For simplicity, the example assumes that the problem
parameters are such that independent of the route sequence, the failure probability at a given customer increases
by 0.20 with respect to the failure probability of the preceding customer.
INSERT FIGURE 1 (¯g1.eps)
9

To implement the s-split procedure, our RC approach uses Algorithm 2. Based on Bellman's algorithm for
directed acyclic graphs, the algorithm splits the TSP tour without generating G0 explicitly. For each vertex
vi 2 V 0 n f0g the algorithm maintains two labels, namely, Zi and Bi. Label Zi holds the expected cost of the
shortest path from 0 to vi and label Bi holds the predecessor of vi in the shortest path. The loops in steps 5 and
9 generate all arcs (feasible routes) in A. The expected cost of each route is evaluated (steps 10 to 24), and its
value is used to update labels Zi and Bi (steps 25 to 31). The value of Zn is the minimal overall expected cost
of the MC-VRPSD solution and the corresponding routes can be built using the labels in Bi starting backwards
from Bn.
INSERT ALGORITHM 2
It is worth mentioning at this point that besides the original s-split procedure we tried two di®erent versions.
The ¯rst extension splits the single permutation in both directions, trying to capture the e®ect of the route
traveling directions described in x4. The second follows the intuition of Bertsimas' cyclic heuristic (Bertsimas,
1992) and generates from the single permutation n di®erent partitions into MC-VRPSD solutions by splitting
the cyclic TSP tour starting from every vertex. Based on some preliminary experimentation, we discarded both
variants simply because the gain on accuracy did not pay o® the loss of simplicity and speed.
5.3 Look-ahead mechanism
Since the s-split procedure only transforms the tours generated by NN or NI into feasible MC-VRPSD solutions
without a®ecting the customer visiting order, it is worth highlighting that the most important decisions are
made during the earlier routing phase. Moreover, as pointed out by Prins (2004) there exists at least one TSP
tour that if known, s-split would translate into the optimal MC-VRPSD solution. Thus the challenge of our RC
approach is to ¯nd such tour during the routing phase. This task, however, is far from trivial, specially using pure
constructive heuristics like NN and NI. Like most of the constructive heuristics for combinatorial optimization,
NN and NI build the ¯nal solution (i.e., a TSP tour) by expanding a partial solution at each iteration according
to a simple and intuitive decision rule. This greedy approach, however, has a major drawback: the decision
rules are based on the impact of decisions on the partial solutions and not on their impact on the ¯nal output
of the algorithm, often leading to myopic moves. To try to overcome this di±culty we coupled the proposed RC
approach with the pilot method (Vo¼ et al., 2005).
The preferred iterative look ahead technique or pilot method works on top of constructive heuristics as a
look-ahead mechanism that helps the subordinate procedures to avoid myopic moves. The basic idea of the
pilot method is to iteratively build the solution to a problem (master solution) by performing K di®erent local
moves (instead of only one) at each iteration and then evaluating their impact in the ¯nal output by using the
10

underlined heuristic to complete the new K partial solutions (pilot solutions) and their merit functions. Then,
the local move leading to the best pilot solution is performed in the master solution and a new iteration begins.
This process is repeated until the master solution is fully built. The pilot method has been successfully applied
to a number of combinatorial optimization problems including Steiner tree problems (Duin and Vo¼, 1999),
the general ring network design problem (Fink et al., 2000), and the °ow shop scheduling problem (Fink and
Vo¼, 2003); however, to the best of our knowledge this is the ¯rst time that a pilot-based approach is applied
to a VRP.
In the proposed look-ahead heuristic, the pilot method works on top of the RC approach as shown in Figure
2. At every iteration, the pilot method makes K copies of the master solution, that is, the partial TSP tour
being built by the routing procedure (NN or NI) and stores them as partial pilot solutions. Then, the pilot
method applies the K best local moves on the partial pilot solutions according to the decision rule of the
routing procedure, that is, the K nearest neighbors for NN and the K nearest insertions for NI. Using the
partial pilot solution as the starting point, the pilot method makes a callback to the routing procedure to
complete the generation of the TSP tour and then invokes s-split to perform the clustering phase. Next, the
partial pilot solution leading to the best MC-VRPSD solution becomes the new master solution and a new
iteration begins. Pilot iterations are carried out until every customer is included in the master solution, then
the s-split procedure performs the clustering phase and the resulting MC-VRPSD solution is the one reported
by the look-ahead heuristic. Since the use of NN and NI as routing procedures leads to two di®erent versions
of our heuristic, henceforth we will refer to each version using the notation Pilot(RP ), where RP 2 fNN;NIg
is the routing procedure used for the routing phase.
INSERT FIGURE 2 (¯g2.eps)
Even though NN and NI were the two only heuristics tried on the routing procedure, any of the constructive
heuristics for the TSP discussed by Bentley (1992) is also suitable to be embedded into our look-ahead heuristic.
Furthermore, s-split inherits from the original Prins' proposition its capability to be easily extended to handle
richer side constraints (Prins, 2004). Indeed, a number of split procedures for vehicle routing reported in the
literature (Villegas et al., 2008; Prins et al., 2009) can replace s-split in the clustering phase, leaving the other
building blocks of the look-ahead heuristic untouched. These facts highlight the °exibility and wide potential
of the proposed look-ahead heuristic.
6 Stochastic 2-Opt
Constructive heuristics for combinatorial optimization, particulary those for VRPs, build feasible solutions
trying to control the objective function without exploring deeply the vast solution space. Consequently, the
11

solely use of constructive algorithms is not advisable to generate competitive solutions for complex problems
like the MC-VRPSD and should be complemented with post-optimization or improvement methods. Therefore,
to enhance the performance of the three proposed heuristics we implemented a post-optimization procedure
based on the well known 2-Opt improvement heuristic for the TSP.
The 2-Opt heuristic evaluates all possible arc exchanges in a route r by deleting two non-adjacent arcs (u; v)
and (u0; v0) and reconnecting the remaining route segments inserting arcs (u; u0) and (v; v0). If an improving
move (arc exchange) is found, that is, the resulting route is feasible and has a cost lower than the original
route, the move is committed and the procedure restarted. The whole operation is repeated until no more
improving arc exchanges are found. For further details on the 2-Opt heuristic the reader is referred to Hansen
and Mladenovi¶c (2006). However, applying 2-Opt to MC-VRPSD solutions is not straightforward because each
arc exchange reverses the visiting order of a subsequence of customers, thus a®ecting the failure probability of
all the customers in the tour. Consequently, the impact of a move on the cost and the feasibility of a solution
can only be evaluated after performing the computationally expensive reestimation of the expected cost of the
route. To address this issue, we designed a special procedure that avoids expensive evaluations of unpromising
2-Opt moves while still explores a su±ciently large number of exchanges.
The proposed post-optimization procedure, henceforth called stochastic 2-Opt (S2-Opt), is described in
Algorithm 3. At every iteration routes r and r0 from the MC-VRPSD solution (R) are merged into a single tour
and all possible arc exchanges in the resulting route are explored (loop in 5). To avoid excessive computations,
every exchange is ¯rst evaluated solely on the basis of the deterministic part of the overall cost, which can
be done in constant time as in the classical 2-Opt procedure. If the exchange is promising, meaning that the
planned cost of the MC-VRPSD solution improves (step 7), the move is then evaluated in terms of the whole
overall cost using s-split (step 8). The tour partitioning process generates a temporary set of feasible routes T
serving the customers originally visited by r and r0. If the expected cost of the new set T is lower than the
sum of those of the two original routes r and r0 (step 9), the move improves the current MC-VRPSD solution
because the expected cost of the other routes (i.e., R n fr; r0g) remains unchanged. In such case, we deem the
move successful and the MC-VRPSD solution is updated (step 10). Finally, the procedure is restarted from the
top following a ¯rst improvement con¯guration. The S2-Opt procedure repeats until it cannot ¯nd any more
improvements.
INSERT ALGORITHM 3
12

7 Computational experiments
We implemented the proposed heuristics in Java, compiled the code using the Java 1.6 compiler embedded
in the Eclipse Java Development Tools (JDT) build I20080617-2000, and run all the experiments on the Java
Runtime Environment (JRE) build 1:6:0 11-b03. Our implementation uses Stochastic Simulation in Java (SSJ)
by L'Ecuyer et al. (2002) for the calculation of probabilities. On its current version, our code supports Normal
and Poisson distributions for the demands; nonetheless, it can be easily extended to handle other probability
distributions having the cumulative property. This property states that the sum of two or more independent
and ª distributed random variables is also ª distributed, as it is the case for the Normal, Poisson, and Gamma
distributions (Christiansen and Lysgaard, 2007; Law and Kelton, 2000). The remainder of this section details
the computational experiments carried out to assess the proposed heuristics on instances from the MC-VRPSD
and VRPSD literature.
7.1 MC-VRPSD instances
To evaluate the power of the proposed heuristics to solve the MC-VRPSD we used the two sets of randomly
generated instances proposed by Mendoza et al. (2010), namely, tuning set and test set. The former is used to
¯ne tune the value of K in the pilot-based heuristics, while the latter is used to compare the performance of
the three algorithms. For the sake of completeness the two sets of instances are brie°y described below.
The tuning set consists of 180 instances generated with the following characteristics. First, 50, 100, and
200 customers and the depot were randomly distributed over a 100 £ 100 Euclidean space. Each client v has
a demand for three di®erent products (p = 1; 2; 3) following a Normal distribution N(¹v;p; ¾v;p), where ¹v;p
takes the value of 10 or 30 and ¾v;p is set such that the coe±cient of variation CV = ¾v;p=¹v;p is either 0.1
or 0.3. The capacity of the compartments was ¯xed by de¯ning a tightness ratio (Pv2Vnf0g ¹v;p)=Qp equal to
10 (Bianchi et al., 2006). Finally, the maximum distance per route was set to L = ¯ £maxv2Vnf0g d0;v, where
¯ is uniformly distributed between 3 and 4. Thirty instances were generated for each combination of number
of customers and coe±cient of variation. Similarly, the test set comprises 180 di®erent instances generated
with the same parameters, except for the mean values of the demands ¹v;p, which were set using a Uniform
distribution U [10; 100].
7.2 Tuning parameter K in the look-ahead heuristics
While SCW is a parameterless heuristic, the pilot-based approaches use parameter K to decide upon the number
of pilot solutions to explore at each iteration. A common strategy in the pilot-method literature is to set K = n,
but in our case it translates into one pilot solution per each customer yet to be visited. Thus this strategy leads
13

to unacceptable execution times for the MC-VRPSD. Vo¼ et al. (2005) overcome this di±culty by limiting what
they call the evaluation depth, that is, to set K = n, but to stop the pilot method when the master solution
reaches a prede¯ned level of completeness. At this point they continue the construction of the solution using
only the subordinate heuristic. Alternatively, we decided not to limit the evaluation depth but to ¯ne tune the
value of K looking after a good tradeo® between accuracy and speed.
To ¯ne tune the value of K we conducted two experiments on the 180 instances in the tuning set. In the
¯rst experiment we conducted Pilot(NN) and Pilot(NI) runs setting K = 1; 2; : : : ; 10. The intuition behind
the selection of these low values is that the decision rules of NN and NI are good in principle, so if better
decisions can be made, they should lie in the vicinity of the ¯rst decision. However, to assess the bene¯ts of
signi¯cant increments in the number of pilot solutions, we tested in the second experiment larger values of K
(K = 15; 20; 25; 30).
Figures 3(a) and 3(b) summarize the results of the 2520 (= 14 £ 180) runs conducted for each heuristic
under the light of two performance metrics. The ¯rst metric is the average gap (over the 180 instances) with
respect to the solution provided by the SCW heuristic, given by
¢Pilot(RP )jSCW =
µE[C]Pilot(RP ) ¡ E[C]SCW
E[C]SCW
¶
£ 100% (7)
where E[C]Pilot(RP ) and E[C]SCW are the total overall cost on a given instance for the corresponding pilot-
based method and SCW, respectively. Henceforth we use Equation (7) to calculate the gap ¢Alg1jAlg2 between
two algorithms Alg1 and Alg2. The second performance metric is the average execution time on a 3GHz Intel
Xeon X5450 CPU with 8Gb of RAM running Windows Vista Ultimate 64bits.
INSERT FIGURE 3 (¯g3a.eps,¯g3b.eps)
The results plotted in Figure 3(a) show the signi¯cant contribution of the pilot method to the accuracy of
the heuristics. When the look-ahead mechanism is turned-o® (i.e., K = 1), the simple RC procedures perform
poorly compared to the SCW heuristic, with average gaps of 9:10% and 8:73% for Pilot(NN) and Pilot(NI),
respectively. However, just exploring one extra pilot solution at each iteration, that is, setting K = 2, reduces
the gap to 3:54% for Pilot(NN) and 4:10% for Pilot(NI) with a minor increment in the average CPU time of
6 and 9 seconds, respectively. As expected, the accuracy of the pilot-based heuristics increases with the value
of K, but from K = 7 the marginal gap °attens while the execution time continues to grow linearly. For
instance, raising K from 1 to 4 reduces the average gap by 6:61% (= 9:10%¡ 2:49%) for Pilot(NN) and 5.82%
(= 8:73% ¡ 2:91%) for Pilot(NI), whereas increasing the same 3 levels from 7 to 10 reduces the gap only by
0:04% and 0:26%, respectively. On the other hand, for each marginal unitary change in K over the range from
1 to 10 the average execution time increases around 5.7 seconds for Pilot(NN) and 8.3 seconds for Pilot(NI).
14

In conclusion, the ¯rst experiment suggests that the most interesting tradeo®s between accuracy and speed are
found with relatively small values of K (i.e., 4 · K · 7). Furthermore, the results of the second experiment
support this conclusion. In Figure 3(b) the gap curve for Pilot(NN) is virtually °at for values of K between
10 and 30. On the contrary, the accuracy of Pilot(NI) slightly increases with the value of K over the same
interval but at a high computational price: the average gap reduces by 1:21% (=1.97%-0.76%) at a threefold
increase in CPU time. In summary and based on these observations, we set the value of K to 6 for Pilot(NN)
and Pilot(NI).
7.3 Results
To assess the performance of the proposed heuristics, we tested SCW, Pilot(NN) and Pilot(NI) on the test
set and compared the results with those reported by the MA of Mendoza et al. (2010). For each of the 180
instances we ran the three heuristics with and without S2-Opt and the pilot-based heuristics with K = 6 (see
x7.2). Table 1 summarizes the results through ¯ve performance metrics for each instance group: the number
of best known solutions (BKSs) found by each approach; the mean gap with respect to the BKS (¢AlgjBKS);
the standard deviation of ¢AlgjBKS ; the third and fourth quartiles for ¢AlgjBKS ; and the execution time (in
seconds) on a 3GHz Intel Xeon X5450 CPU with 8Gb of RAM running Windows Vista Ultimate 64bits.
INSERT TABLE 1
The results in Table 1 show that the post-optimization procedure S2-Opt has a signi¯cant e®ect on the
performance of the three proposed heuristics. By applying the S2-Opt procedure, the average gap with respect
to the BKS was reduced by nearly 25% (from 4:45% to 3:26%) for SCW, and close to 60% for Pilot(NN) (from
6.75% to 2:73%) and Pilot(NI) (from 6:91% to 3:00%). A close look at these results, raises an interesting point
about the computational cost of achieving these improvements. In the worst case, that is, Pilot(NN) running
on the large 200-customer instances, the largest increment in the average execution time due to S2-Opt is just
7 seconds, a marginal computing cost of less than 7.3% (= ((103¡ 96)=96)£ 100%). In conclusion, the S2-Opt
procedure achieves its goal of enhancing the accuracy of the heuristics with a low impact on their simplicity
and speed. Consequently, in the remainder of the paper we will focus on the results obtained by the proposed
algorithms enhanced by S2-Opt, meaning that whenever there is a reference to one of the heuristics, it is implied
that S2-Opt was performed, unless otherwise stated.
The results indicate that the proposed heuristics are able to produce competitive solutions for the MC-
VRPSD with a reduced computational e®ort. A look at the number of BKSs achieved by each approach reveals
that the proposed heuristics set 26 new best known solutions for the 180-instance test set. Out of these 26
new BKSs, 3 (11.5%) were contributed by SCW, 10 (38.5%) by Pilot(NN), and the remaining 13 (50.0%) were
15

found by Pilot(NI). Moreover, the average gap with respect to the BKS in the testbed is only 3:26%, 2:73%,
and 3:00% for SCW, Pilot(NN), and Pilot(NI), respectively. These results are encouraging due to the moderate
computational e®ort required by the proposed heuristics. While the MA of Mendoza et al. (2010) invests between
1.5 minutes (50 customers, CV = 0:1) and 48 minutes (200 customers, CV = 0:3) to generate solutions with
average gaps between 0.00% and 0:44% (high quality), the proposed heuristics only require from 1 to 144 seconds
to obtain solutions with average gaps ranging from 1.48% to 6.18%. After analyzing the performance of the
heuristics by instance size, we found that while in the small 50-customer problems the gap between the proposed
approaches and the BKSs is on average 4:23% (= (6:18%+5:41%+3:65%+3:07%+3:52%+3:54%)=6), the average
gap in the 200-customer instances reduces to 2:07% (= (1:60%+ 1:48%+ 1:60%+ 2:58%+ 1:92%+ 3:22%)=6).
This observation tips the balance in favor of the proposed approaches to solve the MC-VRPSD in practice,
where both fast and accurate solutions are needed.
After comparing the results from the proposed heuristics we conclude that the three approaches are able
to produce solutions of similar quality. A close look at the third and forth quartiles for the gap reported in
Table 1 shows that only for the 50-customer instances there seems to be a di®erence in the accuracy of SCW,
Pilot(NN), and Pilot(NI). In those instances, the pilot-based approaches clearly dominate the SCW heuristic.
To illustrate, for the 50-customer instances with CV = 0:1, Pilot(NN) and Pilot(NI) achieve average gaps of
3.65% and 3.52%, that are clearly superior to the larger average gap of 6.18% by SCW. In addition, the worst
quality gap (fourth quartile) for Pilot(NN) and Pilot(NI), namely, 8:27% and 8:09%, respectively, fall below
the third quartile for SCW (8:29%), providing evidence of the superiority of the pilot-based approaches over
SCW. On the other hand, the SCW heuristic proved to be faster than the pilot-based approaches. Although for
the 50-customer instances the three algorithms report similar execution times (less than 1 second), the better
computing scalability of SCW becomes evident with the larger instance sizes.
7.4 VRPSD instances
Because the VRPSD is a special case of the MC-VRPSD where m = 1, we tested further the e®ectiveness of the
proposed approaches on the set of benchmark problems proposed by Christiansen and Lysgaard (2007). The set,
henceforth called VRPSD instances, comprises 40 problems ranging from 16 to 60 customers adapted from the
classical instances for the capacitated VRP by Augerat (1995) and Christo¯des and Eilon (1969). Christiansen
and Lysgaard (2007) assumed for each customer v a Poisson distributed demand with mean ¹v, equal to the
deterministic demand in the original instance, and kept the remaining data unchanged (i.e., vehicle capacity
and vertex locations). Based on the better performance of the pilot-based heuristics over SCW shown in x7.3
on similarly sized problems (50 customers), we decided to focus solely on pilot-based heuristics for the VRPSD.
To solve the VRPSD instances, we slightly modi¯ed our pilot-based heuristics. On one hand, in the single
16

product scenario the exact value of the expected cost of a route E[Cr] is analytically tractable (Laporte et al.,
2002). Therefore, we modi¯ed steps 10{24 in Algorithm 2 to allow the s-split procedure to evaluate the arcs
in A (feasible routes) using the exact value of the expected cost instead of the proxy provided by the TAP
approximation. On the other hand, previously published results for the VRPSD instances consider a formulation
with an additional constraint that forbids the total expected demand served by a route to exceed the vehicle's
capacity. To compare against these results, we took advantage of the °exibility of s-split and adapted Pilot(NN)
and Pilot(NI) to handle this new side constraint. As in Laporte et al. (2002) and Christiansen and Lysgaard
(2007) the maximum expected-demand constraint can be modeled as:
X
v2r
E[»v] · Q 8r 2 R (8)
Therefore, just by adding the condition Pv2r ¹v · Q in steps 9 and 25 of Algorithm 2, s-split and consequently
the pilot-based heuristics are easily adapted to handle the additional constraint.
We applied Pilot(NN) and Pilot(NI) with the standard con¯guration (i.e., K = 6) on the VRPSD instances
and compared the results with those reported by Christiansen and Lysgaard (2007). Table 2 summarizes the
results on a subset of 19 instances (out of the 40) for which Christiansen and Lysgaard (2007) reported optimal
solutions. Column 1 presents the instance following the template name S-n-V, where S 2 fA;P;Eg is the source
(A and P for Augerat sets A and P, and E for Christo¯des and Eilon), n is the number of nodes, and V the
minimum number of vehicles required; column 2 reports the total expected cost of the optimal solution (OPT);
column 3 shows the CPU time in seconds originally reported by the authors; column 4 reports the CPU time
scaled to our machine using the benchmarks of Dongarra (2009); and, columns 5 to 10 summarize the results for
the proposed methods Pilot(NN) and Pilot(NI) under the light of the total expected cost of the solution (E[C]),
the CPU time in seconds, and the gap with respect to the optimal solution (¢OPT ). Similarly, Table 3 presents
the results on the 21 instances without proven optimum. For these instances, we assess the quality of our
approaches by calculating the gaps with respect to the lower bound (LB) and the best integer solution (BCL)
reported by Christiansen and Lysgaard (2007). The latter is given either by: (i) the best solution returned by
their algorithm under a 1200-seconds time limit (396 seconds in scaled CPU) or (ii) the expected cost of the
optimal solution for the deterministic instance. To the best of our knowledge, the BCLs are the best known
solutions for these 21 VRPSD instances.
INSERT TABLE 2
INSERT TABLE 3
Similarly to the results obtained on the test set for the MC-VRPSD (see x7.3), our experiments on the
VRPSD show that our pilot-based heuristics are able to generate good quality solutions in short CPU times.
17

On the instances solved to optimality by Christiansen and Lysgaard (2007), Pilot(NN) and Pilot(NI) reported
average optimality gaps of 5.06% and 2.53% in short average execution times of 1.20 and 2.42 seconds, respec-
tively. Moreover, the data shows that in 58% (11/19) of the instances Pilot(NN) delivers a solution within 5%
of the optimal cost. Pilot(NI) reports even better results, achieving gaps lower than 5% on 95% (18/19) of the
instances. On the instances without proven optimum, the comparison with respect to the best integer solutions
reported by Christiansen and Lysgaard (2007) shows that our look-ahead approaches set 12 new best known
solutions (bold numbers in column E[C]) and reported small average gaps with respect to the previously best
known solutions of 0.36% and 0.31% for Pilot(NN) and Pilot(NI), respectively. In addition, the data con¯rms
that our pilot-based heuristics produce integer solutions for the instances without provably optimal solutions
within an average gap of just 4.5% with respect to the lower bounds reported by Christiansen and Lysgaard
(2007).
Even though the parameter tuning study (see x7.2) suggested that K = 6 o®ers an interesting tradeo®
between accuracy and speed for our pilot-based heuristics, it also showed that higher values of K may lead
to better solutions if one is willing to pay the extra computational burden. Based on the latter observation,
we selected Pilot(NI), our best performing approach, and conducted an additional experiment setting K =
minf20; ng. The results reported in Table 4 show that under the new con¯guration for K, Pilot(NI) can
improve further its results on the VRPSD instances. For the instances with provably optimal solutions, the
average optimality gap reduces by 0.57% (= 2:53% ¡ 1:96%) with respect to the results obtained under the
default con¯guration. Even better improvements are found on the instances without provably optimal solutions,
where the average gap with respect to the lower bounds reduces by 1% (= 4:54% ¡ 3:53%). In addition, for
the instances without provably optimal solutions, Pilot(NI) with K = minf20; ng improves 11 of the solutions
reported by Christiansen and Lysgaard (2007) {negative values in column ¢BCL{ and unveils 6 new best known
solutions for the subset.
In summary, the results obtained both on the MC-VRPSD and VRPSD instances are remarkable taking into
account the simplicity of the pilot-based heuristics, the good quality of the solutions, and the low computational
requirements.
INSERT TABLE 4
8 Concluding remarks and future research
This paper proposes simple and e®ective constructive heuristics for the MC-VRPSD. The problem consists of
designing collection routes to serve a set of customers who have independent stochastic demands for di®erent
products that must be transported in independent compartments due to product incompatibility constraints.
18

Despite its practical applicability the MC-VRPSD has not received the deserved attention in the literature.
We proposed a set of three constructive heuristics for the MC-VRPSD. The Stochastic Clarke and Wright
heuristic (SCW) extends the savings-based approach of Dror and Trudeau (1986) for the VRPSD to: (i) work
on the multi-compartment scenario, and (ii) handle a di®erent recourse action, that is, a round trip to the
depot to restore the capacity of the vehicle compartments whenever it is exceeded. The other two heuristics
are based on the principle of coupling route-¯rst, cluster-second (RC) procedures with a look-ahead mechanism
known as the pilot method (Vo¼ et al., 2005). To the best of our knowledge, this is the ¯rst proposal of a
pilot-based method in the ¯eld of vehicle routing. The RC procedures use for the routing phase two fast and
easy-to-implement construction algorithms for the deterministic TSP, namely, the nearest neighbor (NN) and
the nearest insertion (NI) heuristics. The clustering phase uses a dynamic programming algorithm (s-split)
specially tailored to handle the stochastic nature of demands in the complex multi-compartment scenario. The
pilot method iterates with the RC procedures to explore further the evaluation of each move in the NN and NI
heuristics by (i) considering K alternatives instead of only one, and (ii) measuring the impact of each of the K
choices in the ¯nal output of the solution rather than in the partial tour. The two resulting heuristics, namely,
Pilot(NN) and Pilot(NI), along with the SCW heuristic were complemented with S2-Opt, a new and e®ective
post-optimization procedure for the MC-VRPSD.
The three proposed heuristics were tested on a large set of instances for the MC-VRPSD and their results
compared with those obtained by the memetic algorithm (MA) of Mendoza et al. (2010). The experiments
indicate that the heuristics are competitive, both in terms of accuracy and speed. The proposed heuristics
followed by S2-Opt were able to produce in short execution times, ranging from 1 to 144 seconds (on average),
solutions with an average gap (with respect to the BKSs) of 3:26%, 2:73%, and 3:00% for SCW, Pilot(NN) and
Pilot(NI), respectively. Moreover, the three proposed algorithms unveiled new best known solutions for 26 out
of the 180 instances in the test set.
In addition to the experiments conducted on the MC-VRPSD, we tested the pilot-based heuristics on a
40-instance testbed for the VRPSD, a particular case of the MC-VRPSD in which each customer has a random
demand for a single product. The solutions found by the two heuristics were compared with those reported by
Christiansen and Lysgaard (2007). The results showed that Pilot(NN) and Pilot(NI) set 12 new best known
solutions for the testbed and reported gaps within 5% with respect to lower bounds reported in (Christiansen
and Lysgaard, 2007) on 24 (60%) and 38 (95%) instances, respectively. In conclusion, thanks to their accuracy,
speed and simplicity the proposed approaches can be used as standalone solution methods or to induce a warm
startup for metaheuristics and exact methods.
Research currently underway includes (i) the extension of the look-ahead heuristics to solve other VRP
variants, (ii) development of new pilot-based methods to work on top of other constructive heuristics for VRPs,
19

and (iii) the implementation of a version of s-split based on simulation that allows the constructive heuristics to
handle the MC-VRPSD with correlated demands at the customer site. Some problems that remain unsolved but
pose interesting challenges are: the derivation of a lower bound for the MC-VRPSD; an extension to the scenario
where preventive trips to the depot (restocking) are allowed; the inclusion of a constraint over the maximum
number of routes (limited °eet); the extension to handle di®erent recourse actions such as cooperation between
routes upon route failures; and the study of the MC-VRPSD under the scenario in which each product can be
delivered to a customer by a di®erent route.
Acknowledgments:
The authors would like to thank the referees and editors for their insightful comments and suggestions. This
research was partially funded by grant 471-2007 from COLCIENCIAS/SENA (Colombia) and project Lig¶eRO
funded by the Pays de la Loire Regional Council (France).
References
Ak, A. and Erera, A. (2007). A paired-vehicle recourse strategy for the vehicle-routing problem with stochastic
demands, Transportation Science 41(2): 222{237.
Augerat, P. (1995). Approche polyµedrale du problµeme de tourn¶ees de v¶ehicules, PhD thesis, Institut National
Polytechnique de Grenoble-INPG.
Bentley, J. B. (1992). Fast algorithms for geometric traveling salesman problem, INFORMS Journal on Com-
puting 4: 387{411.
Bertsimas, D. (1992). A vehicle routing problem with stochastic demand, Operations Research 40(3): 574{585.
Bianchi, L., Birattari, M., Chiarandini, M., Manfrin, M., Mastrolilli, M., Paquete, L., Rossi-Doria, O. and
Schiavinotto, T. (2006). Hybrid metaheuristics for the vehicle routing problem with stochastic demands,
Journal of Mathematical Modelling and Algorithms 5(1): 91{110.
Caramia, M. and Guerriero, F. (2010). A milk collection problem with incompatibility constraints, Interfaces
40(2): 130{143.
Chajakis, E. D. and Guignard, M. (2003). Scheduling deliveries in vehicles with multiple compartments, Journal
of Global Optimization 26(1): 43{78.
Christiansen, C. and Lysgaard, J. (2007). A branch-and-price algorithm for the capacitated vehicle routing
problem with stochastic demands, Operations Research Letters 35(6): 773{781.
20

Christo¯des, N. and Eilon, S. (1969). An algorithm for the vehicle dispatching problem, Operations Research
20(3): 309{318.
Cordeau, J. F., Gendreau, M., Laporte, G., Potvin, J. Y. and Semet, F. (2002). A guide to vehicle routing
heuristics, Journal of the Operational Research Society 53(5): 512{522.
Cordeau, J. F., Laporte, G., Savelsbergh, M. and Vigo, D. (2006). Vehicle routing, in C. Barnhart and G. Laporte
(eds), Handbooks in Operations Research and Management Science: Transportation, Vol. 14, Elsevier,
Amsterdam, pp. 367{428.
Derigs, U., Gottlieb, J., Kalko®, J., Piesche, M., Rothlauf, F. and Vogel, U. (2010). Vehicle routing with
compartments: Applications, modelling and heuristics, OR Spectrum. In press. DOI: 10.1007/s00291-
010-0194-3.
Dongarra, J. J. (2009). Performance of various computers using standard linear equations software, Technical
Report CS-89-85, Computer Science Department, University of Tennessee.
Dror, M. and Trudeau, P. (1986). Stochastic vehicle routing with modi¯ed savings algorithm, European Journal
of Operational Research 23(2): 228{235.
Duin, C. and Vo¼, S. (1999). The Pilot method: a strategy for heuristic repetition with application to the
Steiner problem in graphs, Networks 34: 181{191.
El Fallahi, A., Prins, C. and Wol°er-Calvo, R. (2008). A memetic algorithm and a tabu search for the multi-
compartment vehicle routing problem, Computers & Operations Research 35(5): 1725{1741.
Fink, A., Schneidereit, G. and Vo¼, S. (2000). Solving general ring network design problems by meta-heuristics,
inM. Laguna and J. L. Gonz¶alez-Velarde (eds), Computing tools for modeling, optimization, and simulation,
Kluwer Academic Publishers, Boston, pp. 91{113.
Fink, A. and Vo¼, S. (2003). Solving the continuous °ow-shop scheduling problem by metaheuristics, European
Journal of Operational Research 151(2): 400{414.
Fischetti, M. and Lodi, A. (2003). Local branching, Mathematical Programming 98(1{3): 23{47.
Fisher, M. L. and Jaikumar, R. (1981). A generalized assignment heuristic for the vehicle routing problem,
Networks 11(2): 109{124.
Gendreau, M., Laporte, G. and S¶eguin, R. (1996a). Stochastic vehicle routing, European Journal of Operational
Research 88(1): 3{12.
21

Gendreau, M., Laporte, G. and S¶eguin, R. (1996b). A tabu search heuristic for the vehicle routing problem
with stochastic demands and customers, Operations Research 44(3): 469{477.
Golden, B., Raghavan, S. and Wasil, E. (eds) (2008). The Vehicle Routing Problem: Latest Advances and New
Challenges, Springer US.
Hansen, P. and Mladenovi¶c, N. (2006). First vs. best improvement: An empirical study, Discrete Applied
Mathematics 154(5): 802{817.
Laporte, G. (2009). Fifty years of vehicle routing, Transportation Science 43(4): 408{416.
Laporte, G., Louveaux, F. and Van Hamme, L. (2002). An integer L-shaped algorithm for the capacitated
vehicle routing problem with stochastic demands, Operations Research 50(3): 415{423.
Law, A. M. and Kelton, W. D. (2000). Simulation Modeling and Analysis, McGraw-Hill, Singapore.
L'Ecuyer, P., Meliani, L. and Vaucher, J. (2002). SSJ: A framework for stochastic simulation in Java., in
E. YÄucesan, C. H. Chen, J. L. Snowdon and J. M. Charnes (eds), Proceedings of the 2002 Winter Simulation
Conference, IEEE Press, pp. 234{242.
Liu, C.-S. and Lai, M.-Y. (2009). The vehicle routing problem with uncertain demand at nodes, Transportation
Research Part E: Logistics and Transportation Review 45(4): 517{524.
Mendoza, J. E., Castanier, B., Gu¶eret, C., Medaglia, A. L. and Velasco, N. (2008). Approximating the expected
cost of recourse on a multi-compartment vehicle routing problem with stochastic demands, Technical Report
08/03/AUTO, Institut de Recherche en Communications et Cybern¶etique de Nantes (IRCCyN). ¶Ecole des
Mines de Nantes, France.
Mendoza, J. E., Castanier, B., Gu¶eret, C., Medaglia, A. L. and Velasco, N. (2009). A simulation-based MOEA
for the multi-compartment vehicle routing problem with stochastic demands, Proceedings of the VIII Meta-
heuristics International Conference (MIC). Hamburg, Germany.
Mendoza, J. E., Castanier, B., Gu¶eret, C., Medaglia, A. L. and Velasco, N. (2010). A memetic algorithm for the
multi-compartment vehicle routing problem with stochastic demands, Computers & Operations Research
37(11): 1886{1898.
Muyldermans, L. and Pang, G. (2010a). A guided local search procedure for the multi-compartment capacitated
arc routing problem, Computers & Operations Research 37(9): 1662{1673.
Muyldermans, L. and Pang, G. (2010b). On the bene¯ts of co-collection: Experiments with a multi-compartment
vehicle routing algorithm, European Journal of Operational Research 206(1): 93{103.
22

Prins, C. (2004). A simple and e®ective evolutionary algorithm for the vehicle routing problem, Computers &
Operations Research 31(12): 1985{2002.
Prins, C., Labadi, N. and Reghioui, M. (2009). Tour splitting algorithms for vehicle routing problems, Interna-
tional Journal of Production Research 47(2): 507{535.
Rei, W., Gendreau, M. and Soriano, P. (2007). Local branching cuts for the 0-1 integer L-shaped algorithm,
Technical Report CIRRELT-2007-23, Centre Interuniversitaire de Recherche sur les R¶eseaux d'Entreprise,
la Logistique et le Transport (CIRRELT).
Rei, W., Gendreau, M. and Soriano, P. (2010). A hybrid Monte Carlo local branching algorithm for the single
vehicle routing problem with stochastic demands, Transportation Science 44(1): 136{146.
Ruiz, R., Maroto, C. and Alcaraz, J. (2004). A decision support system for a real vehicle routing problem,
European Journal of Operational Research 153(3): 593{606.
Savelsbergh, M. and Goetschalckx, M. (1995). A comparison of the e±ciency of ¯xed versus variable vehicle
routes, Journal of Business Logistics 16: 163{187.
Secomandi, N. and Margot, F. (2009). Reoptimization approaches for the vehicle-routing problem with stochastic
demands, Operations Research 57(1): 214{230.
SÄorensen, K. and Sevaux, M. (2009). A practical approach for robust and °exible vehicle routing using meta-
heuristics and Monte Carlo sampling, Journal of Mathematical Modeling and Algorithms 8(4): 387{407.
Tan, K. C., Cheonga, C. Y. and Goh, C. K. (2007). Solving multiobjective vehicle routing problem with
stochastic demand via evolutionary computation, European Journal of Operational Research 177(2): 813{
839.
Tatarakis, A. and Minis, I. (2009). Stochastic single vehicle routing with a prede¯ned customer sequence and
multiple depot returns, European Journal of Operational Research 197(2): 557{571.
Toth, P. and Vigo, D. (eds) (2002). The Vehicle Routing Problem, Monographs on discrete mathematics and
applications, SIAM.
Van der Bruggen, L. Gruson, R. and Salomon, M. (1995). Reconsidering the distribution structure of gasoline
products for a large oil company, European Journal of Operational Research 81(3): 460{473.
Villegas, J., Medaglia, A. L., Mendoza, J. E., Prins, C., Prodhon, C. and Velasco, N. (2008). A split-based
framework for the vehicle routing problem, XIV Latin Ibero-American Congress on Operations Research
(CLAIO). Cartagena, Colombia.
23

Vo¼, S., Fink, A. and Duin, C. (2005). Looking ahead with the pilot method, Annals of Operations Research
136(1): 285{302.
Yang, W. H., Mathur, K. and Ballou, R. (2000). Stochastic vehicle routing with restocking, Transportation
Science 34(1): 99{112.
24

Appendix
A Calculation of the failure probability in the multi-compartment
scenario
Note that in the formulation given by (2) and (3) in x2 no constraints restrict the number of failures in route
r. However, since demand realizations ~» are nonnegative and less than the capacity of the corresponding
compartment, if a route failure occurs at customer vi it is in the worst case the (i ¡ 1)th failure on route r.
For notational simplicity, let us de¯ne z(vi) as a function that takes the value of 1 if a route failure occurs in
route r during the service of customer vi and the value of 0 otherwise. Since all compartments are emptied after
each visit to the depot, traveling a route is a process that is renewed after each failure. Therefore, the value of
Pr¡z(vi) = 1
¢ can be written as:
Pr
¡
z(vi) = 1
¢
= Pr
¡
z (vi) = 1; z (vk) = 0; 0 < k < i
¢
+
i¡1X
j=1
Pr (z (vi) = 1; z(vj) = 1; z (vk) = 0; j < k < i) (9)
where the ¯rst term represents the probability that the ¯rst route failure occurs at customer vi and the second
term the probability that at least one failure occurs between the start of the route at depot and customer vi.
Because the ¯rst renewal process occurs at the depot before visiting the ¯rst customer, Pr¡z(0) = 1¢ = 1, and
(9) can be rewritten as:
Pr¡z (vi) = 1
¢ =
i¡1X
j=0
Pr (z (vi) = 1; z(vj) = 1; z (vk) = 0; j < k < i) (10)
Considering that the last failure occurs at customer vj , an expression to calculate each term in (10) can be
derived as follows. Let A(p) be the event that capacity Qp for product p is exceeded while servicing customer vi,
and B the event that no route failures occurred between vertices vj and vi given that the last failure occurred
at customer vj (i.e., z(vj) = 1; z(vk) = 0; 8j < k < i). The probabilities of events A(p) and B are given by
equations (11) and (12) respectively.
Pr¡A (p) ¢ = Pr
0
@»0vj ;p +
iX
k=j+1
»vk;p > Qp
1
A (11)
Pr¡B¢ = Pr
0
@»0vj ;p +
i¡1X
k=j+1
»vk;p · Qp;8 p 2 P
1
A (12)
were »0vj ;p is the unattended demand of customer vj for product p after the last route failure. Figure 4 describes
the probability tree at arrival to customer vi.
25

INSERT FIGURE 4 (¯g4.eps)
From the probability tree it can be seen that
Pr¡z(vi) = 0jz(vj) = 1; z(vk) = 0; j < k < i
¢ = ¡1¡ Pr(A(1)jB)¢£ ¡1¡ Pr(A(2)jB)¢£ : : :
£ ¡1¡ Pr(A(m)jB)¢
(13)
and
Pr¡z(vi) = 1jz(vj) = 1; z(vk) = 0; j < k < i
¢ = 1¡ Pr¡z(vi) = 0jz(vj) = 1; z(vk) = 0; j < k < i
¢
= 1¡
h¡1¡ Pr(A(1)jB)¢£ ¡1¡ Pr(A(2)jB)¢£ : : :£
¡1¡ Pr(A(m)jB)¢
i
(14)
From Bayes' theorem, for product 1 we have
Pr¡A(1)jB¢ = Pr
¡A(1); B¢
Pr(B)
=
Pr
0
@»0vj ;1 +
i¡1X
k=j+1
»vk;1 · Q1 < »0vj ;1 +
iX
k=j+1
»vk;1 ; »0vj ;p +
i¡1X
k=j+1
»vk;p · Qp ; p 2 P n f1g
1
A
Pr
0
@»0vj ;p +
i¡1X
k=j+1
»vk;p · Qp ; p 2 P
1
A
=
Pr
0
BB@»0vj;1+
i¡1X
k=j+1
»vk;1·Q1<»0vj;1+
iX
k=j+1
»vk;1
1
CCA£Pr
0
BB@»0vj;2+
i¡1X
k=j+1
»vk;2·Q2
1
CCA£:::£Pr
0
BB@»0vj;m+
i¡1X
k=j+1
»vk;m·Qm
1
CCA
Pr
0
BB@»0vj;1+
i¡1X
k=j+1
»vk;1·Q1
1
CCA£Pr
0
BB@»0vj;2+
i¡1X
k=j+1
»vk;2·Q2
1
CCA£:::£Pr
0
BB@»0vj;m+
i¡1X
k=j+1
»vk;m·Qm
1
CCA
=
Pr
0
BB@»0vj;1+
i¡1X
k=j+1
»vk;1·Q1<»0vj;1+
iX
k=j+1
»vk;1
1
CCA
Pr
0
BB@»0vj;1+
i¡1X
k=j+1
»vk;1·Q1
1
CCA
(15)
26

where Pr
³
»0vj ;1 +
Pi¡1
k=j+1 »vk;1 · Q1
´
is the cumulative probability that the total demand for product 1 col-
lected between customers vj and vi¡1 after failing at vertex vj does not exceed the capacity of the compartment
Q1. The numerator in (15) can be rewritten as:
Pr
0
@»0vj ;1 +
i¡1X
k=j+1
»vk;1 · Q1 < »0vj ;1 +
iX
k=j+1
»vk;1
1
A = Pr
0
@»0vj ;1 +
i¡1X
k=j+1
»vk;1 · Q1
1
A¡Pr
0
@»0vj ;1 +
iX
k=j+1
»vk;1 · Q1
1
A (16)
Re-writing expression (15) in terms of (16) and simplifying we obtain
1¡ Pr¡A(1)jB¢ =
Pr
³
»0vj ;1 +
Pi
k=j+1 »vk;1 · Q1
´
Pr
³
»0vj ;1 +
Pi¡1
k=j+1 »vk;1 · Q1
´ (17)
The result on (17) extends to any p 2 P. Hence, (13) can be rewritten as
Pr¡z(vi) = 0jz(vj) = 1; z(vk) = 0; j < k < i
¢ =
mY
p=1
Pr
³
»0vj ;p +
Pi
k=j+1 »vk;p · Qp
´
Pr
³
»0vj ;p +
Pi¡1
k=j+1 »vk;p · Qp
´ (18)
and (14) as
Pr¡z(vi) = 1jz(vj) = 1; z(vk) = 0; j < k < i
¢ = 1¡
mY
p=1
Pr
³
»0vj ;p +
Pi
k=j+1 »vk;p · Qp
´
Pr
³
»0vj ;p +
Pi¡1
k=j+1 »vk;p · Qp
´ (19)
From Bayes' theorem
Pr¡z(vi) = 1jz(vj) = 1; z(vk) = 0; j < k < i
¢ = Pr
¡z(vi) = 1; z(vj) = 1; z(vk) = 0; j < k < i
¢
Pr¡z(vj) = 1; z(vk) = 0; j < k < i
¢ (20)
where
Pr¡z(vj) = 1; z(vk) = 0; j < k < i
¢ =
mY
p=1
Pr
0
@»0vj ;p +
i¡1X
k=j+1
»vk;p · Qp
1
A£ Pr¡z(vj) = 1
¢ (21)
replacing (19) and (21) in (20) and simplifying we obtain
Pr¡z(vi) = 1; z(vj) = 1; z(vk) = 0; j < k < i
¢ =
" mY
p=1
Pr
Ã
»0vj ;p +
i¡1X
k=j+1
»vk;p · Qp
!
¡
mY
p=1
Pr
Ã
»0vj ;p +
iX
k=j+1
»vk;p · Qp
!#
£ Pr¡z(vj) = 1
¢
(22)
27

¯nally, replacing (22) in (10) we have that the probability of having a failure in route r while servicing customer
vi 2 r is equal to:
Pr¡z(vi) = 1
¢ =
i¡1X
j=0
2
4
mY
p=1
Pr
0
@»0vj ;p +
i¡1X
k=j+1
»vk;p · Qp
1
A¡
mY
p=1
Pr
0
@»0vj ;p +
iX
k=j+1
»vk;p · Qp
1
A
3
5£Pr¡z(vj) = 1
¢
(23)
¥
28

A
rc
(r
ou
te
)
ev
al
ua
ti
on
Ar
c
Se
qu
en
ce
l r
E
h G
r
³ ~ »
´i
E
[C
r]
E
[C
r]
·
L
a,
a
0,
a,
0
10
+
10
=
20
2£
(0
£
10
)
=
0
20
tr
ue
a,
b
0,
a,
b,
0
10
+
15
+
15
=
40
2£
(0
£
10
+
0:
2£
15
)
=
6
46
tr
ue
a,
c
0,
a,
b,
c,0
10
+
15
+
20
+
20
=
65
2£
(0
£
10
+
0:
2£
15
+
0:
4£
20
)
=
22
87
fa
lse
a,
d
0,
a,
b,
c,d
,0
10
+
15
+
20
+
10
+
30
=
85
2£
(0
£
10
+
0:
2£
15
+
0:
4£
20
+
0:
6£
30
)
=
58
14
3
fa
lse
a,
e
0,
a,
b,
c,d
,e,
0
10
+
15
+
20
+
10
+
10
+
20
=
85
2£
(0
£
10
+
0:
2£
15
+
0:
4£
20
+
0:
6£
30
+
0:
8£
20
)
=
90
17
5
fa
lse
b,
b
0,
b,
0
15
+
15
=
30
2£
(0
£
15
)
=
0
30
tr
ue
b,
c
0,
b,
c,0
15
+
20
+
20
=
55
2£
(0
£
15
+
0:
2£
20
)
=
8
63
tr
ue
b,
d
0,
b,
c,d
,0
15
+
20
+
10
+
30
=
75
2£
(0
£
15
+
0:
2£
20
+
0:
4£
30
)
=
32
10
7
fa
lse
b,
e
0,
b,
c,d
,e,
0
15
+
20
+
10
+
10
+
20
=
75
2£
(0
£
15
+
0:
2£
20
+
0:
4£
30
+
0:
6£
20
)
=
56
13
1
fa
lse
c,c
0,
c,0
20
+
20
=
40
2£
(0
£
20
)
=
0
40
tr
ue
c,d
0,
c,d
,0
20
+
10
+
30
=
60
2£
(0
£
20
+
0:
2£
30
)
=
12
72
fa
lse
c,e
0,
c,d
,e,
0
20
+
10
+
10
+
20
=
60
2£
(0
£
20
+
0:
2£
30
+
0:
4£
20
)
=
28
88
fa
lse
d,
d
0,
d,
0
30
+
30
=
60
2£
(0
£
30
)
=
0
50
tr
ue
d,
e
0,
d,
e,0
30
+
10
+
20
=
60
2£
(0
£
30
+
0:
2£
20
)
=
8
68
tr
ue
e,e
0,
e,0
20
+
20
=
40
2£
(0
£
20
)
=
0
40
tr
ue
Fi
gu
re
1:
S-
sp
lit
pr
oc
ed
ur
ef
or
clu
ste
rin
g
29

Figure 2: Pilot method on top of the RC procedure
30

(a) Tradeo® between accuracy and speed for K = 1; :::; 10
(b) Tradeo® between accuracy and speed for K = 10; 15; 20; 25; 30
Figure 3: Tuning the value of K on the proposed pilot-based heuristics
31

Figure 4: Probability tree of events at arrival to customer vi given that the last failure occurred at customer vj
32

Algorithm 1 The stochastic Clarke and Wright heuristic
1: RÃ ;
2: for every customer v 2 V n 0 do
3: create route r = f0; v; 0g
4: RÃ R[ frg
5: end for
6: create the savings list S Ã ;
7: for every route r 2 R do
8: u á v1 2 r and v á vnr 2 r
9: for every route r0 2 R n frg do
10: u0 á v1 2 r0 and v0 á vnr0 2 r0
11: calculate su;u0 ; su;v0 ; sv;u0 ; sv;v0
12: if minfE[Cu;u0 ]; E[Cu0;u]g · L then S Ã S [ fsu;u0g
13: if minfE[Cu;v0 ]; E[Cv0;u]g · L then S Ã S [ fsu;v0g
14: if minfE[Cv;u0 ]; E[Cu0;v]g · L then S Ã S [ fsv;u0g
15: if minfE[Cv;v0 ]; E[Cv0;v]g · L then S Ã S [ fsv;v0g
16: end for
17: end for
18: smax = max(u;v)2Vnf0g£Vnf0g fsu;vjsu;v 2 Sg
19: if smax < 0 or S = ; then
20: stop
21: end if
22: if smax ¸ 0 then
23: (u¤; v¤) Ã argmax(u;v)2Vnf0g£Vnf0g fsu;vjsu;v 2 Sg
24: merge routes having u¤ and v¤ as extreme vertices r à merge (r(u¤); r(v¤))
25: RÃ R n fr(u¤)g
26: RÃ R n fr(v¤)g
27: RÃ R[ frg
28: go to 6
29: end if
33

Algorithm 2 S-split algorithm for tour partitioning
1: set Z0 = 0 and Zi = +1, 1 · i · n
2: set Bi = 0, 1 · i · n
3: set Pr(v0) = 1
4: initialize dummy variables x = 1 and y = 1
5: for i = 1 to n do
6: lr = 0
7: E
h
Gr
³¡!»
´i
= 0
8: j = i
9: while j · n and lr + E
h
Gr
³¡!»
´i
· L do
10: if j = i then
11: lr = d0;vj + dvj ;0
12: else
13: lr = lr ¡ dvj¡1;0 + dvj¡1;vj + dvj ;0
14: end if
15: Pr(vj) = 0
16: for k = 0 to j ¡ 2 do
17: x = 1, y = 1
18: for product p = 1 to m do
19: x = x£ Pr
³Pj¡1
l=k+1 »vl;p · Qp
´
20: y = y £ Pr
³Pj
l=k+1 »vl;p · Qp
´
21: end for
22: Pr(vj) = Pr(vj) + (x¡ y)£ Pr(vk)
23: end for
24: E
h
Gr
³~»
´i
= E
h
Gr
³~»
´i
+ 2£ dvj ;0 £ Pr(vj)
25: if lr + E
h
Gr
³~»
´i
· L then
26: if Zi¡1 + lr + E
h
Gr
³~»
´i
· Zj then
27: Zj = Zi¡1 + lr + E
h
Gr
³~»
´i
28: Bj = i¡ 1
29: end if
30: j = j + 1
31: end if
32: end while
33: end for
34

Algorithm 3 Stochastic 2-Opt procedure for post-optimization
1: for every route r 2 R do
2: for every route r0 2 R n r do
3: concatenate copies of r and r0 to generate r00
4: let W be the set of all possible 2-Opt exchanges in r00
5: for every exchange w 2 W do
6: make a copy of r00 and perform w to generate r000
7: if lr000 < lr00 then
8: call s-split and partition r000 to generate T
9: if E[C(T )] < E[Cr] + E[Cr0 ] then
10: RÃ T [R n fr; r0g
11: go to 1
12: end if
13: end if
14: end for
15: end for
16: end for
35

50
Cu
st
om
er
s
10
0
Cu
st
om
er
s
20
0
Cu
st
om
er
s
To
ta
l/
Av
g.
M
et
ric
Al
go
rit
hm
C.
V.
=
0.
1
C.
V.
=
0.
3
C.
V.
=
0.
1
C.
V.
=
0.
3
C.
V.
=
0.
1
C.
V.
=
0.
3
#
BK
S
SC
W
0
0
0
0
0
0
0
Pi
lo
t(
NN
)
0
0
0
0
0
0
0
Pi
lo
t(
NI
)
0
0
0
0
0
0
0
SC
W
+
S2
-O
pt
0
1
0
0
0
2
3
Pi
lo
t(
NN
)+
S2
-O
pt
0
0
1
1
6
2
10
Pi
lo
t(
NI
)+
S2
-O
pt
1
0
4
3
4
1
13
M
A
29
29
25
26
20
25
15
4
Av
g.
¢
A
lg
jB
K
S
SC
W
7.
60
7.
66
2.
97
3.
53
2.
54
2.
38
4.
45
Pi
lo
t(
NN
)
6.
50
5.
67
7.
33
5.
49
7.
40
8.
08
6.
75
Pi
lo
t(
NI
)
6.
39
5.
91
7.
34
5.
97
7.
74
8.
13
6.
91
SC
W
+
S2
-O
pt
6.
18
5.
41
2.
29
2.
57
1.
60
1.
48
3.
26
Pi
lo
t(
NN
)+
S2
-O
pt
3.
65
3.
07
2.
88
2.
63
1.
60
2.
58
2.
73
Pi
lo
t(
NI
)+
S2
-O
pt
3.
52
3.
54
2.
99
2.
80
1.
92
3.
22
3.
00
M
A
0.
05
0.
00
0.
15
0.
11
0.
44
0.
07
0.
14
St
d.
D
ev
.¢
A
lg
jB
K
S
SC
W
2.
94
2.
84
1.
31
1.
48
1.
25
1.
66
3.
04
Pi
lo
t(
NN
)
2.
05
1.
97
1.
99
1.
57
1.
46
1.
90
2.
05
Pi
lo
t(
NI
)
2.
42
1.
66
2.
10
1.
82
1.
81
1.
44
2.
07
SC
W
+
S2
-O
pt
2.
60
2.
63
1.
23
1.
40
1.
20
1.
48
2.
61
Pi
lo
t(
NN
)+
S2
-O
pt
2.
06
1.
41
1.
84
1.
64
1.
35
1.
81
1.
79
Pi
lo
t(
NI
)+
S2
-O
pt
2.
03
1.
43
1.
76
1.
77
1.
83
1.
82
1.
84
M
A
0.
29
0.
00
0.
45
0.
43
0.
78
0.
18
0.
45
¢
A
lg
jB
K
S
SC
W
[9
.3
4,
15
.4
4]
[9
.7
8,
12
.8
7]
[3
.8
6,
6.
43
]
[4
.5
5,
7.
02
]
[3
.3
2,
5.
88
]
[2
.6
8,
8.
14
]
-
[3
th
,4
th
]q
ua
rt
ile
s
Pi
lo
t(
NN
)
[8
.1
5,
10
.0
1]
[6
.4
7,
11
.9
0]
[8
.5
1,
11
.3
8]
[6
.2
8,
9.
73
]
[8
.3
7,
10
.9
2]
[8
.9
0,
14
.0
1]
-
Pi
lo
t(
NI
)
[7
.8
1,
11
.9
0]
[6
.3
8,
13
.0
8]
[8
.2
3,
14
.0
0]
[6
.9
4,
10
.5
0]
[8
.9
5,
11
.2
8]
[8
.9
8,
11
.5
2]
-
SC
W
+
S2
-O
pt
[8
.2
9,
11
.7
5]
[6
.6
5,
12
.1
9]
[2
.9
5,
4.
98
]
[3
.1
8,
6.
78
]
[2
.2
4,
5.
44
]
[1
.7
0,
6.
50
]
-
Pi
lo
t(
NN
)+
S2
-O
pt
[4
.7
7,
8.
27
]
[3
.9
7,
7.
10
]
[3
.3
0,
8.
49
]
[3
.4
6,
6.
08
]
[2
.6
0,
4.
78
]
[3
.8
7,
5.
69
]
-
Pi
lo
t(
NI
)+
S2
-O
pt
[4
.6
1,
8.
09
]
[4
.2
3,
6.
39
]
[4
.4
2,
6.
25
]
[4
.5
3,
5.
51
]
[3
.2
1,
6.
83
]
[4
.5
3,
6.
10
]
-
M
A
[0
.0
0,
1.
57
]
[0
.0
0,
0.
00
]
[0
.0
0,
2.
04
]
[0
.0
0,
2.
28
]
[0
.8
5,
2.
69
]
[0
.0
0,
0.
77
]
-
Av
g.
CP
U
SC
W
<
1
<
1
8
8
70
77
27
Pi
lo
t(
NN
)
<
1
<
1
8
8
96
96
35
Pi
lo
t(
NI
)
<
1
<
1
11
11
13
7
13
8
50
SC
W
+
S2
-O
pt
<
1
<
1
8
8
73
80
28
Pi
lo
t(
NN
)+
S2
-O
pt
<
1
<
1
9
9
10
3
10
3
37
Pi
lo
t(
NI
)+
S2
-O
pt
<
1
<
1
12
12
14
4
14
2
52
M
A*
11
5
11
9
50
3
50
4
25
54
28
82
11
13
Ta
ble
1:
Pe
rfo
rm
an
ce
m
et
ric
s
fo
r
th
e
te
st
se
t.
#
BK
S:
nu
mb
er
of
be
st
kn
ow
n
so
lut
ion
sf
ou
nd
.A
vg
.¢
A
lg
jB
K
S
(%
):
av
era
ge
ga
p
wi
th
res
pe
ct
to
th
eb
est
kn
ow
n
so
lut
ion
.S
td
.D
ev
.¢
A
lg
jB
K
S:
sta
nd
ar
d
de
via
tio
n
of
th
eg
ap
wi
th
res
pe
ct
to
th
eb
est
kn
ow
n
so
lut
ion
.
¢ A
lg
jB
K
S
[3t
h,
4t
h]
qu
ar
til
es:
th
ird
an
d
fou
rth
qu
ar
til
es
for
th
eg
ap
wi
th
res
pe
ct
to
th
eb
est
kn
ow
n
so
lut
ion
.
Av
g.
CP
U
(s)
:
av
era
ge
ex
ec
ut
ion
tim
e.
*O
n
aP
C
wi
th
a2
.4G
Hz
Co
re
2D
uo
pr
oc
ess
or
,4
GB
of
RA
M
,a
nd
ru
nn
ing
W
ind
ow
sX
P
Pr
ofe
ssi
on
al
36

Christiansen and Lysgaard (2007) Pilot(NN) Pilot(NI)
Instance OPT(=LB) CPU* CPU** E[C] CPU ¢OPT E[C] CPU ¢OPT
P-16-8 512.82 0.00 0.00 515.00 0.03 0.43 516.47 0.06 0.71
P-22-8 681.06 0.00 0.00 701.14 0.06 2.95 698.69 0.13 2.59
E-22-4 411.57 1.00 0.33 419.15 0.16 1.84 428.28 0.20 4.06
P-23-8 619.52 1.00 0.33 635.80 0.08 2.63 630.80 0.13 1.82
A-39-5 869.18 3.00 0.99 897.98 3.02 3.31 902.78 10.80 3.87
P-21-2 218.96 5.00 1.65 240.38 0.45 9.78 228.21 0.34 4.22
A-33-5 704.20 8.00 2.64 747.34 0.44 6.13 726.21 0.61 3.13
P-40-5 472.50 9.00 2.97 506.55 1.42 7.21 504.89 2.22 6.85
A-33-6 793.90 49.00 16.17 846.74 0.36 6.66 801.15 0.58 0.91
E-33-4 850.27 86.00 28.38 925.54 1.27 8.85 866.27 10.58 1.88
P-19-2 224.06 153.00 50.49 224.06 0.23 0.00 224.06 0.23 0.00
P-22-2 231.26 219.00 72.27 273.38 0.30 18.22 237.12 0.38 2.53
A-39-6 876.60 279.00 92.07 897.95 1.06 2.44 889.98 1.50 1.53
A-32-5 853.60 282.00 93.06 889.74 0.44 4.23 863.47 0.80 1.16
P-20-2 233.05 352.00 116.16 237.06 0.33 1.72 234.26 0.28 0.52
P-51-10 809.70 430.00 141.90 871.87 0.81 7.68 831.84 2.14 2.73
P-55-15 1068.05 792.00 261.36 1106.66 0.75 3.61 1097.81 1.67 2.79
A-45-7 1264.83 882.00 291.06 1341.87 10.38 6.09 1318.44 11.08 4.24
P-60-15 1085.49 1348.00 444.84 1111.86 1.30 2.43 1113.01 2.27 2.54
Avg.(%) 257.84 85.09 1.20 5.06 2.42 2.53
Min.(%) 0.00 0.00 0.03 0.00 0.06 0.00
Max.(%) 1348.00 444.84 10.38 18.22 11.08 6.85
< 3% 8/19 13/19
< 5% 11/19 18/19
Table 2: Performance metrics for the VRPSD instances with proven optimal solution. OPT(=LB): optimal
solution reported by Christiansen and Lysgaard (2007). CPU* (s): execution time on a PC with a 1.5GHz Pentium
Centrino processor and 480MB of RAM (operational system not reported). CPU** (s): scaled execution time. E[C]:
total expected cost of the solution. CPU (s): execution time on a PC with a 2.4GHz Core 2 Duo processor, 4GB of
RAM, and running Windows XP. ¢OPT (%): optimality gap
37

Christiansen and Pilot(NN) Pilot(NI)
Lysgaard (2007)
Instance LB BCL E[C] CPU ¢LB ¢BCL E[C] CPU ¢LB ¢BCL
A-34-5 825.26 827.87 826.87 0.52 0.19 -0.12 849.83 0.78 2.98 2.65
A-36-5 852.09 907.55 897.04 0.89 5.28 -1.16 892.54 1.11 4.75 -1.65
A-37-5 707.54 708.34 753.26 1.38 6.46 6.34 727.48 1.38 2.82 2.70
A-37-6 1030.44 1030.75 1051.80 0.80 2.07 2.04 1048.70 0.98 1.77 1.74
A-38-5 761.12 778.09 777.82 10.30 2.19 -0.04 782.87 10.69 2.86 0.61
A-44-6 1021.29 1025.48 1065.12 10.41 4.29 3.87 1065.58 11.05 4.34 3.91
A-45-6 1006.88 1096.19 1050.49 1.88 4.33 -4.17 1050.49 1.66 4.33 -4.17
A-46-7 999.87 1002.41 1008.75 3.02 0.89 0.63 1031.57 2.61 3.17 2.91
A-48-7 1180.22 1248.27 1279.56 1.42 8.42 2.51 1256.85 2.55 6.49 0.69
A-53-7 1109.34 1180.10 1190.52 10.69 7.32 0.88 1151.23 11.89 3.78 -2.45
A-54-7 1279.93 1342.87 1334.63 4.28 4.27 -0.61 1338.39 4.80 4.57 -0.33
A-55-9 1173.56 1264.18 1227.57 10.56 4.60 -2.90 1202.94 11.77 2.50 -4.84
A-60-9 1503.65 1608.40 1547.27 4.34 2.90 -3.80 1568.74 12.83 4.33 -2.47
E-51-5 544.63 553.26 588.63 3.63 8.08 6.39 594.58 12.16 9.17 7.47
P-45-5 527.77 546.05 569.82 2.25 7.97 4.35 562.53 11.34 6.59 3.02
P-50-10 757.44 792.20 783.93 1.03 3.50 -1.04 783.10 1.69 3.39 -1.15
P-50-7 581.19 606.41 613.14 10.56 5.50 1.11 626.63 2.75 7.82 3.33
P-50-8 666.90 724.69 690.98 0.98 3.61 -4.65 702.97 2.19 5.41 -3.00
P-55-10 741.25 797.21 767.90 1.77 3.60 -3.68 768.09 11.78 3.62 -3.65
P-55-7 585.47 616.44 630.39 5.00 7.67 2.26 620.60 12.33 6.00 0.68
P-60-10 799.58 831.24 826.42 3.13 3.36 -0.58 835.97 12.48 4.55 0.57
Avg.(%) 4.23 4.60 0.36 6.71 4.54 0.31
Min.(%) 0.52 0.19 -4.65 0.78 1.77 -4.84
Max.(%) 10.69 8.42 6.39 12.83 9.17 7.47
< 3% 5/21 17/21 5/21 17/21
< 5% 13/21 19/21 15/21 20/21
Table 3: Performance metrics for the VRPSD instances without proven optimal solution. LB: lower bound
reported by Christiansen and Lysgaard (2007). BCL: best integer solution reported by Christiansen and Lysgaard (2007).
E[C]: total expected cost of the solution. CPU (s): execution time on a PC with a 2.4GHz Core 2 Duo processor, 4GB
of RAM, and running Windows XP. ¢LB (%): Gap with respect to LB. ¢BCL (%): Gap with respect to BCL
38

Instances with proven optimal solution Instances without proven optimal solution
Instance E[C] CPU ¢OPT Instance E[C] CPU ¢LB ¢BCL
P-16-8 515.00 0.09 0.43 A-34-5 826.87 1.53 0.19 -0.12
P-22-8 700.08 0.23 2.79 A-36-5 890.13 2.27 4.46 -1.92
E-22-4 419.15 0.39 1.84 A-37-5 724.62 2.75 2.41 2.30
P-23-8 628.66 0.28 1.47 A-37-6 1056.76 1.83 2.55 2.52
A-39-5 895.24 2.78 3.00 A-38-5 782.17 17.75 2.77 0.52
P-21-2 228.21 0.63 4.22 A-44-6 1054.32 18.7 3.23 2.81
A-33-5 716.71 1.34 1.78 A-45-6 1050.49 3.25 4.33 -4.17
P-40-5 501.78 3.41 6.20 A-46-7 1006.49 3.77 0.66 0.41
A-33-6 801.15 1.30 0.91 A-48-7 1246.76 3.98 5.64 -0.12
E-33-4 862.71 1.92 1.46 A-53-7 1159.61 5.8 4.53 -1.74
P-19-2 224.06 0.44 0.00 A-54-7 1311.12 7.66 2.44 -2.36
P-22-2 237.12 0.70 2.53 A-55-9 1198.96 20.55 2.16 -5.16
A-39-6 879.52 2.41 0.33 A-60-9 1555.63 23.23 3.46 -3.28
A-32-5 859.62 1.45 0.71 E-51-5 564.15 21.77 3.58 1.97
P-20-2 234.26 0.52 0.52 P-45-5 553.04 19.55 4.79 1.28
P-51-10 817.90 4.27 1.01 P-50-10 781.05 3.64 3.12 -1.41
P-55-15 1076.31 3.99 0.77 P-50-7 608.11 7.39 4.63 0.28
A-45-7 1327.79 18.84 4.98 P-50-8 702.97 4.34 5.41 -3.00
P-60-15 1110.26 5.53 2.28 P-55-10 759.35 5.16 2.44 -4.75
P-55-7 622.97 22.06 6.40 1.06
P-60-10 839.57 22.38 5.00 1.00
Avg.(%) 2.66 1.96 10.45 3.53 -0.66
Min.(%) 0.09 0.00 1.53 0.19 -5.16
Max.(%) 18.84 6.20 23.23 6.40 2.81
< 3% 15/19 8/21 21/21
< 5% 18/19 17/21 21/21
Table 4: Performance of Pilot(NI) with K = minf20; ng on the VRPSD instances. E[C]: total expected cost
of the solution. CPU (s): execution time on a PC with a 2.4GHz Core 2 Duo processor, 4GB of RAM, and running
Windows XP. ¢OPT (%): Gap with respect to the optimal solution reported by Christiansen and Lysgaard (2007). ¢LB
(%): Gap with respect to the lower bound reported by Christiansen and Lysgaard (2007). ¢BCL (%): Gap with respect
to the best integer solution reported by Christiansen and Lysgaard (2007).
39