An iterated local search meta-heuristic
for a new class of pickup and delivery
problem
An iterated local search meta-heuristic for a new
class of pickup and delivery problem
Zedi Liu s1796917
August 13, 2015
Abstract
Contents
1 Introduction 2
2 Literature review 4
3 Problem description 6
4 An introduction to iterated local search 10
5 The proposed heuristic 14
5.1 Initial solution construction . . . 15
5.2 Local search operators . . . 18
5.3 Perturbation mechanism . . . 19
5.4 Feasibility check . . . 21
5.5 Integration of components . . . 22
5.6 Streamlining the ILS . . . 24
6 Computational experimentation and results 26 6.1 Test data . . . 26
6.2 Parameter determination . . . 28
6.3 Experiment settings . . . 31
6.4 The influence of side constraints . . . 33
6.5 Results . . . 33
1
Introduction
The vehicle routing problem (VRP), first introduced half a century ago by Dantzig and Ramser [6], has been extensively studied by academics in the fields of oper-ations research and transportation over the past few decades. During this time, a number of variants of the basic VRP, along with numerous sophisticated exact algorithms and heuristics have been put forward. Although various modifications can be made to the basic VRP in order to better model different real life scenar-ios, at its core, the same fundamental question remains: what is the optimal set of routes for a fleet of vehicles to traverse in order to service a given set of cus-tomers? Clearly, the VRP has many potential applications, especially in today’s globalized economy where large quantities of goods often need to be transported among geographically dispersed locations via vast road networks. Optimizing the routes can lead to substantial reductions in transportation costs.
The VRP and all its variants are generalizations of the well-known traveling sales-man problem, which is NP-hard, therefore all VRP variants too are NP-hard. Even with the help of powerful modern computers, existing exact algorithms are only capable of solving small instances of the VRP optimally. In practice, heuristics and meta-heuristics are often the only options available for obtaining decent solutions to the VRP and its variants of any realistic size [1].
more sophisticated heuristic could potentially yield better results.
In this paper, we consider an aspect the GPDP that was not part of the main focus in the study of Buijs et al. Specifically, we ignore the collaboration between different depots, and direct our effort to improving the solution quality around one isolated depot, while keeping most of the assumptions and problem specifications of the GPDP. Thus, the objectives of this thesis are:
In the context of the single-depot GPDP, we attempt to develop an effective meta-heursitic algorithm to find high quality solutions. Using the algorithm, we analyze the characteristics of the single-depot GPDP. The idea of applying heuristics to solve complicated problems is not a new one. In the literature, articles treating heuristic solution methods abound. Meta-heuristics are generic solution schemes designed to be adaptable to a variety of related prob-lems. Among popular meta-heuristics for VRP related problems, local search based methods such as variable neighborhood search, tabu search, and simulated anneal-ing have a well-established track record. After scouranneal-ing the literature and noticanneal-ing some similarities between the single-depot GPDP and the vehicle routing problem with simultaneous pickup and delivery (VRPSPD), we decided to construct and apply a heuristic based on the iterated local search (ILS), which is a well-known meta-heuristic that has been shown to work very well on the VRPSPD [18]. The embedded local search metohd in the ILS is adapted from the variable neighbor-hood descent (VND) method.
2
Literature review
In this section, we present a very brief summary of some relevant works in the lit-erature. Regarding the general classifications of the VRP and its variants, we refer the reader to Toth and Vigo (2002)[19]. For the PDP and its variants, we suggest the survey papers by Parrah et al. (2008) [15], Berbeglia [1] and Savelsbergh and Sol (1995)[17]. Savelsbergh and Sol (1995) divide the General Pickup and Delivery Problem into four categories, which are Static Single -Vehicle PDP, Static Multi-Vehicle PDP, Dynamic Single -Multi-Vehicle PDP and Dynamic Multi-Multi-Vehicle PDP. For an overview of various proposed meta-heuristic solution approaches for the VRP variants, a categorized bibliography by Gendreau et al. (2007)[7] is a good source to start with.
The most similar problem to the problem we study in this thesis is the VRPSPD defined as follows (see [14], [18] and[20]): a set of customers are located on a trans-portation network; each customer requires either a delivery or a pick-up operation (or both) of a certain amount and each must be visited once for both operations. The service is provided by a set of vehicles of limited capacity; each vehicle leaves the depot carrying an amount of goods equal to the total amount it must deliver and returns to the depot carrying an amount equal to the total amount it has picked up. In each point along its tour each vehicle cannot carry a total load greater than its capacity. The goal is to minimize the overall length of the tours, or the number of vehicles used, or both.
Lau and Liang [12] studied the pickup and delivery problem with time windows (PDPWT) in which each request specifies a quantity to be picked up at an origin and delivered to a destination, constituting a pickup-delivery pair. Each route starts and end at a single depot and serves these pickup-delivery pairs, a vehicle leaves the depot empty and returns empty.
3
Problem description
The single-depot GPDP based on the GPDP in Buijs et al. (2014)[4] has the following elements:
• a fleet of homogeneous vehicles, i.e. all vehicles are identical with regard to capacity and maximum range. The number of vehicles in the fleet is suf-ficiently large as to preclude situations that some requests are left unfulfilled;
• a single depot serves as both the origin and the destination of all vehicles in the fleet. In other words, all feasible routes start from and end at the depot. It also serves as a warehouse where goods can be stored overnight;
• a fixed number of clients located at distinct locations around the depot.
To model these elements, we define a complete directed graph G = (V, E), where V = {0, 1, . . . , n} is the set of n clients and 1 depot, with 0 designating the depot, and E = {(i, j)|i, j ∈ V, i 6= j} is the set of edges.
Each client i ∈ V, i 6= 0 is associated with 3 types of requests:
• type A: a fixed quantity ai ≥ 0, to be delivered to i from the depot;
• type B: a fixed quantity bi ≥ 0, to be picked up at i and delivered to the
depot;
• type C: one or more fixed quantities ci,j ≥ 0, j 6= 0, j 6= i, to be picked up at
i and delivered to client j without transshipment.
i is identified by a unique pair of coordinates (oi, pi).
The vehicles all have a maximum capacity of Q and a maximum range of dmax.
Each edge (i, j) ∈ E represents the direct path between clients i and j, its length di,j is the Euclidean distance between i and j:
di,j =
q
(oi− oj)2+ (pi− pj)2
The use of Euclidean distances implies that the triangular inequality holds: di,j < di,k+ dk,j ∀k ∈ V and k 6= i, j ,
Also
di,j = dj,i ∀i, j ∈ V
di,i= 0 ∀i ∈ V
These distances are collected in a non-negative symmetric matrix D = (di,j).
If all non-zero requests are type A requests or if all non-zero requests are type B requests, the problem is reduced to the classical VRP. If all clients have either non-zero type A requests or type B requests, the problem is VRP with backhauls, if some clients have both non-zero type A and type B requests, the problem can be viewed as the VRPSPD. In the single-depot GPDP, a client can simultaneously have non-zero requests of all three types. Furthermore, clients are served in arbi-trary orders.
Since each tour begins and ends at the depot and we assume no transshipment, all type C requests will need to be fulfilled by one vehicle in one single tour, thus a non-zero ci,j induces a precedence relation between clients i and j: i and j need to
in later section Test data under the chapter computational experimentation and reuslts.
The objective is to minimize the total length of routes in a solution, we do not explicitly consider the minimization of the number of vehicles in our objective. We define the following decision variables:
xijk =
1, if edge (i, j) belong to route k 0, otherwise
m : the number of routes/vehicles
yij : the total amount of requests picked up from clients up to and including i,
to be transported via edge (i, j).
zij : the total amount of requests to be delivered to clients after i,
to be transported via edge (i, j).
Minimize m X k=1 X i∈V X j∈V dijxijk (1) subject to X i∈V m X k=1 xijk= 1 ∀j = 1, . . . , n (2) X i∈V xijk− X i∈V xjik = 0 ∀j = 1, . . . , n ∀k = 1, . . . , m (3) X j∈V xk0j ≤ 1 ∀k = 1, . . . , m (4) X i∈V X j∈V dijxijk≤ dmax ∀k = 1, . . . , m (5) X i∈V yji− X i∈V yij = bj ∀j 6= 0 (6) X i∈V zij − X i∈V zji = aj ∀j 6= 0 (7) yij + zij ≤ Q m X k=1 xijk ∀i, j ∈ V (8) xijk ∈ {0, 1} ∀i, j ∈ V ∀k = 1, . . . , m (9) yij ≥ 0 ∀i, j ∈ V (10) zij ≥ 0 ∀i, j ∈ V (11)
Constraints (2) ensure that each client is visited exactly once by exactly one ve-hicle. Constraints (3) guarantee that the same vehicle arrives and departs from every client it serves. Constraints (4) define that at most m vehicles are used. Constraints (5) are the vehicle range constraints. (6) and (7) are flow constraints that ensure all type A and type B requests are satisfied. Constraints (8) establish that requests will only be transported via edges included in the solution, they further impose that the total amount of requests on a vehicle does not exceed the vehicle’s capacity on any given section of the route. (Remaining task: define precedence orders constraints)
A solution is represented as a collection of routes s = {r1, r2, . . . , rm}, where each
routes in a solution is necessarily equal to the number of vehicles dispatched. We do not assume a fixed number of available vehicles nor do we impose any condition on the number of vehicles used in a feasible solution, in principle 1 ≤ m ≤ n − 1. The problem instances considered throughout this thesis are static and determin-istic, all the pickup and delivery requests are known before the initialization of our solution method.
Solving the single-depot GPDP involves two decisions: assignment decision that determines which clients are assigned to which vehicle, routing decision that deter-mines the best order in which the clients assigned to each vehicle should be served.
4
An introduction to iterated local search
For many classic combinatorial optimization problems including the VRP and its variants, the number of possible solutions grows exponentially with the problem size. Exhaustively checking every one of these possible solutions to find an opti-mal solution usually makes the problem intractable. Often the best one can hope for is trying to find a solution that is ”good enough” for practical purposes with-out exceeding the capabilities of computing tools. A meta-heuristic is a general algorithmic framework for addressing this issue.Typically, meta-heuristics are ap-proximation algorithms that cannot always produce provable optimal solutions, rather they are used for obtaining good solutions in short amounts of time. Meta-heuristics can be classified into two broad categories: trajectory based and population based. Trajectory based meta-heuristics create and maintain one sin-gle solution at each iteration, while population based meta-heuristics maintain a population of several solutions simultaneously [2].
[10]).
In the remainder of this chapter, we give a general introduction to the iterated local search method. The information presented here is covered in much greater detail with a few extensions in [11] and [10].
Let f be the cost function of a generic minimizing combinatorial optimization problem. Let s denote feasible solutions and S the set of all s, we call S the solution space of the problem. A solution so ∈ S is optimal if
f (so) ≤ f (s) ∀s ∈ S A locally optimal minimum s∗ is such that
f (s∗) ≤ f (s) ∀s ∈ N (s∗) ∩ S
where N (s∗) denotes a neighborhood of s∗. A neighborhood structure is defined for all s ∈ S, in the context of discrete optimization, it usually consists of all solutions that can be obtained from s by a specific type of simple modification. If there are many local minima, the range of values they span may be large. And the global minimum so may differ substantially from the average value of local
minimums obtained by some simple optimization routine/heuristic.
Let S∗ be the set containing all local minimums s∗. Obviously S∗ ⊂ S. An opti-mization routine can be seen as a mapping from S to the much smaller S∗. The basin of attraction of a local minimum B(s∗) is defined as the set of s that lead an optimization routine to s∗. The ILS helps its embedded optimization routines escape basins of attraction of local minimums so that other prospective areas of S∗ can be explored until a satisfactory solution is found.
Algorithm 1 Basic iterated local search 1: s0 ← GenerateInitialSolution 2: s∗ ← LocalSearch(s0) 3: repeat 4: s0 ← Perturbation(s∗, history) 5: s∗0 ← LocalSearch(s0) 6: s∗ ← AcceptanceCriterion(s∗,s∗0 , history)
7: until Termination condition met
GenerateInitialSolution is responsible for generating an initial solution. It is usu-ally a greedy construction heuristic, popular candidates for VRP variant problems include the nearest neighbor heuristic, insertion heuristics, savings heuristics, etc. LocalSearch is made up of one or more optimization routines. It systematically explores neighborhoods of a given solution in the hope of improving upon it. The ultimate goal of the LocalSearch is to find a locally optimal solution among the neighborhoods of a given solution.
solu-tion is a global minimum. This is the main drawback of the ILS and any other meta-heuristic. However, the ILS will stochastically converge to a global minimum.
Figure 1: An illustration of ILS
Figure 1 (taken from [11]) gives a visual representation of one iteration of the ILS. A local minimum s∗ has been found. It is then perturbed to obtain s0, leaving the basin of attraction of s∗. Applying local search to s0 will eventually give another local minimum s∗0.
5
The proposed heuristic
5.1
Initial solution construction
For the initial solution generator, we have chosen to use a modified nearest neighbor heuristic. There are two other types of heuristics for generating initial solutions, the cheapest insertion heuristic [9], and the sweeping heuristic [8]. The fact we need to respect precedences of nodes presents difficulties in implementing these types of heuristics, we therefore resort to the nearest neighbor heuristic.
Our nearest neighbor heuristic constructs new routes sequentially, starting a new route only when the previous route is completed. When constructing a route, all clients that are not in a route are considered, the decision of which client to add to the route is based on both a candidate client’s distance to the last previously added client and the feasibility of adding the candidate client.
Because the vehicle’s load fluctuates in an unpredictable manner as new clients are put in its route, simply adding the nearest neighbor of the last client in a par-tial route might violate the vehicle load constraint, therefore we sort all unrouted clients {i|i ∈ V, i /∈ rk, f ork = 1, . . . , m} by their distances to the last client in the
Algorithm 2 Initial Solution Construction
1: procedure GenerateInitialSolution
2: s ← ∅
3: CL ← set of all clients
4: rp ← {0} . rp is the working partial route
5: while CL 6= ∅ do
6: sort CL
7: for i ∈ CL do
8: if r0 = rp + {i} feasible then
9: if ci,j > 0 then 10: if r00 = r0 + {j} feasible then 11: rp ← r00 12: CL ← CL − {i} − {j} 13: else 14: rp ← rp 15: CL ← CL 16: end if 17: else 18: rp ← r0 19: CL ← CL − {i} 20: end if 21: else 22: rp ← rp 23: CL ← CL 24: end if 25: end for
26: rc← rp + {0} . complete the new route
27: s ← s + {rc} . add new route to solution
28: r ← {0}
29: end while
5.2
Local search operators
The local search routine in our ILS algorithm has an intricate structure of its own consisting of 6 operators. An operator is a specialized procedure that applies small changes to a given solution. The 6 operators are organized into the local search routine to form a variable neighborhood descent (VND) algorithm, which system-atically explores and modifies the neighborhood structures in a deterministic way. We have chosen 7 operators found in to building the variable neighborhood descent local search routine. Of the 7 chosen local search operators, 6 are inter-route op-erators, meaning they perform movements of clients between a pair of routes; the 6th operator is an intra-route operator, meaning it moves clients within a given route. We emphasize that all chosen operators are of the so-called edge exchange operators. Each operator defines a neighborhood structure. Their descriptions are given below. The same operators have been used in [20] and [18]:
(1,0)Shift: A client i on route r1 to a route r2. First the edges joining i, (i − 1, i),
(i, i + 1) on route r1 and (j, j + 1) on route r2 are deleted, then clients i − 1 and
i + 1 are joined by adding a new edge (i − 1, i + 1), i is inserted onto r2 by two
new edges (j − 1, i) and (i, j).
Crossover: The edge between two adjacent clients i and i + 1, (i, i + 1) on route r1, and the edge between two adjacent clients j and j + 1, (j, j + 1) on route r2
are deleted. The severed parts are then reconnected to form two new routes by linking i and j + 1 by (i, j + 1), j and i + 1 by (j, i + 1).
(1,1)Swap: An exchange of a client i from route r1 and j from r2. The edges on
route r1 (i − 1, i) and (i, i + 1) are replaced by (i − 1, j) and (j, i + 1) respectively.
Likewise on route r2, the edges (j − 1, j) and (j, j + 1) are replaced by (j − 1, i)
and (i, j + 1) respectively.
(2,0)Shift:Two adjacent clients i and j from the same route r1 are transferred to
(2,1)Swap: Two adjacent clients i and j from route r1 are exchanged with an
isolated client l from route r2. The edge modifications are similar to those in (1,1)
swap.
(2,2)Swap: Two adjacent clients i and j from route r1 are exchanged with two
adjacent clients l and k from r2. The edge modifications are similar to those in
(1,1) swap.
If an improvement of the current incumbent solution is found, the routes in the better solution are subject to further optimization using the 2-Opt operator: 2-Opt: Two nonadjacent arcs in a route are removed and another two are added to form a new route.
Although minimizing the number of vehicles is not part of our objective function, our edge-exchange operators never increase the number of vehicles produced by Algorithm 2, because none of the moves performed by these operators involves cre-ating new routes. In fact, crossover, (1,0)shift, and (2,0)shift can sometimes eliminate a route if transferring all its clients to other routes decreases the total distance traveled.
The pseudo-code in Algorithm 3 shows how these operators are integrated into our VND local search routine.
5.3
Perturbation mechanism
The perturbation mechanism is made up of its own operators that change the structure of a given locally optimal solution (local minimum) so that nearby neigh-borhoods can be explored.
Algorithm 3 Variable Neighborhoods Descent
1: procedure VND
2: s ← InputSolution
3: k ← 1
4: while k ≤ r do . r is the number of inter-route neighborhoods
5: find the best neighor s0of s in neighborhood k
6: if f (s0) < f (s) then
7: s ← s0
8: f (s) ← f (s0)
9: k ← 1
10: intensification in the modified routes
basin of attraction of the incumbent solution, on the other hand, if the step size is too large, perturbation jumps become like a random walk and many promising neighborhood will likely be skipped.
The essential difference between a perturbation operator and a local search oper-ator is that a perturbation operoper-ator does not seek improvements, a new solution is accepted as long as it is feasible.
We make use of two perturbation operators in our perturbation mechanism: Ejection chain Used in [16] for the classic version of the VRP, the ejection chain works as follows. First a client is randomly removed from each route in a solution. Next each of the removed clients is randomly reinserted into a random route dif-ferent than the route to which it originally belonged. Only one client is inserted to each route and clients with predecessors or successors are left untouched. If no feasible solution is found within 50 attempts, we abort the ejection chain and move the next perturbation operator.
Double swap Two (1,1) swaps are performed consecutively at random. The max-imum number of attempts is also 50.
Each time the perturbation mechanism is called, one of these operators is used with a probability of 0.5.
5.4
Feasibility check
The last piece of building block is the Feasibility Check, whose job it is to ensure no constraint is violated at any point on any given route. It is called during the construction of an initial solution, as well as each time a change is made by an operator.
simplest constraints are examined first, only when no violation is found does the feasibility check proceed to more complex constraints.
First, given a route rk, it makes sure the load on each vehicle when leaving the
depot and when arriving at the depot are no greater than the vehicle’s max-imum capacity. Define the initial load on a vehicle which serves route rk as
Linit =
P
i∈rk,i6=0ai, and its final load as Lf in =
P
i∈rk,i6=0bi, this condition is
expressed by:
Linit ≤ Q and Lf in ≤ Q
Next, the total length of rk needs to be less than or equal to the vehicle’s range:
X
i∈V
X
j∈V
dijxijk ≤ D
It should be pointed out that the algorithm does not use xijk as a decision
vari-able, the edges in rkare directly read from its representation. For instance, a route
rk = {0, 1, 2, 3, 0} includes edges (0, 1),(1, 2),(2, 3),and (3, 0).
After that, the vehicle’s load at any given point on rk also cannot exceed the
vehicle capacity: Linit− X j∈rk,j≤i aj + X j∈rk,j≤i bi+ X j∈V cij − X j∈V cji ≤ Q
Finally, the precedence orders need to be maintained:
for all i ∈ rk, i 6= 0 ,if cij > 0 for some j ⇒ j ∈ rk
conversely,
for all j ∈ rk, j 6= 0 ,if cij > 0 for some i ⇒ i ∈ rk
A route is flagged as infeasible if it fails any one of these checks. The feasibility check is the most frequently called procedure in the algorithm. It is crucial that everything be done to ensure its fast performance.
5.5
Integration of components
Algorithm 4 Iterated Local Search
1: procedure ILS
2: s ← GenerateInitialSolution
3: s0 ← s
4: IterILS ← 0
5: while IterILS < MaxIterILS do
6: s ← VND(s) 7: if f (s) < f (s0) then 8: s0 ← s 9: f (s0) ← f (s) 10: IterILS ← 0 11: end if 12: s ← Perturb(s0) 13: IterILS ← IterILS + 1 14: end while 15: s∗ ← s0 16: return s∗ 17: end procedure
The acceptance criterion is to accept any solution better than the incumbent solu-tion s0. Once accepted, the new solution replaces s0as the new incumbent solution. IterILS is a counter that records the number of iterations. When the VND ex-hausts the neighborhoods of a given solution, the perturbation mechanism is called. This is counted as one iteration.
5.6
Streamlining the ILS
We conclude this chapter with a discussion of techniques for speeding up the algo-rithm. A major challenge we face when using the local search operators to explore solution neighborhoods is that a large number of infeasible moves are inevitably attempted, evaluating all these infeasible moves wastes a lot of computing time. However, we can use to alleviate this problem to some extent with the help of a few simple techniques . First, observe the following propositions:
• Proposition 1 : the removal from a feasible route of any client i for which cij = 0 for all j and cki = 0 for all k results in a new feasible route. This
operation doesn’t violate any precedence orders. The vehicle load decreases by ai on the section before i and by bi on the section after i, so the vehicle
load constraint is not violated. By the triangular inequality, the route length is not increased. The same holds for the removal of two consecutive clients under similar circumstances.
• Proposition 2 : A vehicle would have to carry all type A requests before leaving the depot, and be loaded with all type B requests of all clients on its route before arriving at the depot. A new insertion i would render the route infeasible if the capacity constraint is violated by adding ai to the starting
load Linit or by adding bi to the finishing load Lf in, irrespective of i’s position
on the route.
Proposition 1 suggests that we only need to check whether a client i has any prece-dence relations before removing it from a feasible route. Proposition 2 tells us we can perform a prilimiary “feasibility check” without calling the feasibility check procedure by simply verifying:
Linit+ ai ≤ Q and Lf in+ bi ≤ Q
need only consider the clients and edges affected in each step of a local search. Let us consider the following example:
Example 1. Suppose we would like to evaluate a (1,1)exchange between i from route r1 and j from route r2, assuming such a move is feasible. Denote the
im-mediate predecessor and successor of i on route r1 by i − 1 and i + 1 respectively,
j − 1 and j + 1 are similarly defined for j. Client i is to be inserted between j − 1 and j + 1 on route r2 and j between i − 1 and i + 1 on route r1, given that these
insertions are feasible, we calculate the change in the total length of r1 and r2 by:
∆1 = (di−1,j + dj,i+1) − (di−1,i+ di,i+1)
∆2 = (dj−1,i+ di,j+1) − (dj−1,j+ dj,j+1)
∆ = ∆1+ ∆2
In the equations above, ∆1measures the change in the length of r1, ∆2 measures the
change in the length of r2 caused by the exchange of i and j between the two routes.
If their sum ∆ < 0, the combined length of r1 and r2 is shortened, the move is said
to produce an improvement. The change is then accepted and stored as a candidate.
This is less costly than evaluating the entire lengths of the new routes after a change is made. In our implementation of the algorithm, all route altering changes are evaluated in a similar fashion.
Algorithm 5 Accepting first improvement
1: Given a solution s, and a neighborhood structure N
2: repeat
3: Find the next solution s0 ∈ N (s)
4: if f (s0) < f (s) then
5: s ← s0
6: end if
7: until All eligible s0 are examined Algorithm 6 Accepting best improvement
1: Given a solution s, and a neighborhood structure N
2: repeat
3: Find s0 = argmins∈N (s)f (s)
4: if f (s0) < f (s) then
5: s ← s0
6: end if
7: until The neighborhood N is thoroughly searched
6
Computational experimentation and results
6.1
Test data
The clients are located on a 200 × 200 grid, on which the range of coordinates is (−100, 100) both horizontally and vertically; the depot occupies the center of the grid with coordinates (0, 0). A client i (i ∈ V i 6= 0) has integer valued coordinates (oi, pi), the coordinates are sampled from the uniform distribution U (−100, 100),
no two clients are allowed to share the same coordinates.
at a time. Given the way our initial solution is constructed, all type C requests appear in the initial solution in consecutive pairs. So this restriction allows pickup-delivery pairs to be reassigned to a different route than the on they are on in the initial solution.
We give an illustrative example of a precedence list in Table 1:
Table 1: Example of type C requests origin (i) destination (j) quantity (ci,j)
1 2 2
4 5 3
6 7 1
8 9 3
10 11 2
In Table 1, the first column is the set of origins for type C requests, the second the destinations and the third the quantities to be delivered between each pair of origin of destination.
In our randomly generated instances, the type C request pairs are constructed by running the nearest neighbor heuristic once and choosing (not necessarily consec-utive) pairs from the resulting routes. The chosen pairs are then assigned type C quantities.
All type A and type B requests quantities are sampled from the uniform distribu-tion U (0, 15). The vehicle capacity Q are varied in our experiments.
Because the locations of the clients are randomly generated, we find that the ve-hicle’s range dmax needs to be at least twice the diagonal of the grid, otherwise
We choose the number of clients n = 50 and n = 100 to represent medium sized and large problem instances respectively.
6.2
Parameter determination
The only parameter in our algorithm is MaxIter, the maximum number of it-erations without improvement to the incumbent solution before the algorithm is terminated. In this section, we present our analysis for determining this parameter. On a randomly generated instance with n = 50 clients without type C reqeusts, we choose 4 values of MaxIter : 10, 30, 50 and 100. For each value of MaxIter, the algorithm is run 30 times. The distributions of solutions for different values of MaxIter are given in Figure 2. Even at a low MaxIter, the algorithm occa-sionally produces very good results, but our goal is for the algorithm to output results consistently. At MaxIter =100, the average solution is on par with the best solutions found when using lower MaxIter values, although the range of solution distribution is not better than when MaxIter =50, roughly half the solutions are below the average.
On the n = 50 instance (see Figure 2), as MaxIter increases, the average solution decreases slightly, however the gap between the worst and the best solutions nar-rows. For MaxIter =50, the worst solution is 1.8% worse than the average solution and the best solution is 1.92% better than the average solution. For MaxIter =100, the worst is 2.45% worse than the average solution and the gap between the best and the average solutions is 1.35%. We observe the gaps between the extreme solutions and the average solution stabilize at point, even though the
Figure 2: Boxplot of solution quality for different values of MaxIter on the instance with n=100 clients.
par with the best solutions obtained with lower MaxIter values, and the distribu-tion of soludistribu-tions is more or less centered around the average. We fix MaxIter =50 in the rest of our analysis. In Table 2 and Table 3, we summarize the rest of the results not presented in the boxplots, the average runtime, the average solution, the best and the worst solutions and ∆1 = ((worst − average)/average) × 100 is the percentage gap between the worst solution and the average solution, ∆2 = ((average − best)/average) × 100 is the percentage gap between the best solution and the average solution.
Figure 3: Boxplot of solution quality for different values of MaxIter on the instance with n=100 clients.
higher, it is beyond any doubt our algorithm presents very significant improve-ments over the nearest neighbor solutions. Therefore in the subsequent analysis, instead of comparing the final solutions to the initial solutions, we compare the final solutions to the local minimums produced by the algorithm before the first iteration.
A possible explanation for why the algorithm exhibits a wider range of solutions on the smaller instance is that the number of routes and the number of clients in a small instance are both lower, giving the perturbation mechanism a smaller set of feasible moves.
Table 2: Results with different MaxIter for n=50
MaxIter average t average best worst ∆1 ∆2
10 9.9 1514.8 1470.8 1555.4 2.70 2.90
30 27.9 1504.8 1460.8 1534.8 2.00 2.93
50 49.0 1501.1 1472.8 1528.1 1.80 1.92
100 80.5 1490.9 1470.8 1527.6 2.45 1.35
Table 3: Results with different MaxIter for n=100
MaxIter average t average best worst ∆1 ∆2
10 63.7 2414.4 2357.5 2436.7 0.93 2.36
30 242.9 2405.6 2356.7 2430.7 1.02 2.03
50 439.2 2354.0 2330.0 2382.0 1.19 1.00
6.3
Experiment settings
In the literature, when testing approximate solution approaches for VRP and PDP variants, many authors are focused on finding the best solutions from multiple runs. We feel that this practice somewhat defeats the purpose of meta-heuristics, since in most practical situations, we are not afforded the luxury of repeatedly run-ning a meta-heuristic in order to find the best solutions among many candidates. Therefore, it is useful to examine not only the best-case performance, but also the average performance as well as the worst-case performance. A good approxima-tion algorithm should be able to output consistent soluapproxima-tions, too much variaapproxima-tion is generally not a desirable feature.
In the absence of readily available optimal or best solutions in the literature, we resign ourselves to testing the influences of side constraints (vehicle capacity, travel range) and the percentage of type C requests on the quality of solutions.
Nevertheless, in order to gain some limited insight on the performance of our al-gorithm, we construct pickup-delivery pairs in the following way: first set cij = 0
∀i, j ∈ V , reducing the problem to a vehicle routing problem with simultaneous pickup and delivery (VRPSPD). We run the algorithm 30 times and select the best solution that we name snull, then we examine the routes in snull and pick pairs of
locations to form pickup-delivery pairs. To limit the impact of introducing addi-tional requests, if clients i and j are chosen to form a pickup-delivery pair, we set cij ≤ min(ai, bj) and subtract cij from both i’s type A request ai and j’s type B
request bj. This way, we can guarantee snull is at least a feasible, and possibly good
solution to the problem with the additional type C requests. The instances with the highest percentage of type C requests are constructed first, then type C re-quests are removed randomly when needed to generate lower percentage instances, so that the instances in each subgroup share some common type C requests. This method is partly inspired by a similar approach found in [12].
We test our algorithm on 2 groups of 24 instances. One group contains instances with n = 50 clients, the other is made up of instances with n = 100 clients. In each group we differentiate 3 subgroups with regard to the percentage of clients with type C requests: 10%, 20% and 30% (rounded to the nearest even integer). Within each group, the instances share the same coordinates and the same type A and type B requests, that is, only the percentage and distribution of type C requests differ from one subgroup to another.
phase and requires restarts.
6.4
The influence of side constraints
Of the two side constraints: the vehicle load constraint and the vehicle range con-straint, most routes in a feasible solution appear to be load constrained, although we set very generous vehicle range dmax to encourage feasible initial solutions and
facilitate the perturbation mechanism. Decreasing dmax often results in the
algo-rithm freezing after a few iterations, even when the number of clients are low. That said, in all our solutions, the route lengths are always far below dmax, especially
in large instances, the
It is also worth noting that when the distance constraints are restrictive, the routes in good quality incumbent solutions generally have lengths close to the maximum travel distance, even changing nodes between two adjacent routes tend to render the solution infeasible, thus making the use of perturbation mechanism difficult.
6.5
Results
We present the results in Tables 4, 5, 6, 7, 8, 9, 10 and 11. In each test, the solutions are compared against the snull associated with that subgroup of problem
instances, δ1, δ2, δ3 record the percentage increase over snull of the average
solu-tion, the best solution and the worst solution respectively.
Each table contains the results obtained on a subgroup of 3 instances, The in-stances in each consecutive pair of tables (4 and 5),(6 and 7),(8 and 9) (10 and 11) are identical expect for type C requests.
mechanism, whereas when there’s no type C requests, all clients are eligible to be moved. Secondly, as we have discussed in the Heuristic section it also stands to reason that the inclusion of pickup-delivery pairs also hinders the operations of the local search operators to some extend. This conjecture is supported by some evidence when we examine the recorded objective history of the extreme case in the last row in Table 5.
Of 30 runs, the solution with the objective value 1470.8 was found 7 times, 1471.3 was found 14 times, 1473.9 was found 4 times, 1510.9, 3 times, and 1507.2, 1507.6 once each. Following our argument in subsection Experiment settings, we know the solutions 1471.3 and 1470.8 are not the best solutions for this instance. There-fore, one can deduce the algorithm is likely trapped in a basin of attraction. The usually short average runtime t=40.3, which is a sign of premature termination of the algorithm, lends some further credence to this conjecture.
Also, comparing the results on instances with different number of clients n, this variation in solution qualities seems to be more pronounced on smaller instances with n = 50 clients.
Unfortunately, we don’t believe this drawback can be fixed in any simple way. Af-ter a few iAf-terations, the pickup-delivery pairs are usually no longer consecutive to each other, so even if the perturbation mechanism targets consecutive pairs rather than single clients, it is very unlikely many if any pickup-delivery pairs will be selected for a move.
Figure 4: The null solution for Q=200, n=50 instances
Table 4: Experiment 1, Q=100, snull= 1459.1.
n type C % average δ1 best δ2 worst δ3 average t
50 10 1513.6 3.73 1464.5 0.37 1566.2 7.34 67.2
50 20 1568.8 7.51 1546.2 5.97 1638.8 12.3 55.7
50 30 1548.8 6.14 1539.4 5.50 1561.6 7.02 53.3
In Figure 5b, we present a visual comparison between the best solution from the last row in Table 4 where 30% of clients have type C requests. Its associated null solution snull (with an objective value of 1459.1) is shown in Figure 5b.
The algorithm was written in M atLab, the experiments were performed on an Intel i7 4790K CPU with 16 GB of RAM.
7
Conclusion and suggestions for future research
Table 5: Experiment 2, Q=100, snull= 1459.1.
n type C % average δ1 best δ2 worst δ3 average t
50 10 1500.7 2.85 1471.3 0.84 1509.8 3.47 75.6
50 20 1488.6 2.02 1470.8 0.80 1524.8 4.50 72.7
50 30 1477.9 1.29 1470.8 0.80 1510.9 3.55 40.3
Table 6: Experiment 3, Q=200, snull= 1240.3.
n type C % average δ1 best δ2 worst δ3 average t
50 10 1242.5 0.18 1240.3 0 1244.0 0.29 58.6
50 20 1279.1 3.12 1266.2 2.08 1335.1 7.64 54.2
50 30 1269.9 2.38 1264.7 1.97 1287.0 3.77 62.1
Table 7: Experiment 4, Q=200, snull= 1240.3.
n type C % average δ1 best δ2 worst δ3 average t
50 10 1243.0 0.21 1240.3 0 1244.0 0.29 61.2
50 20 1259.1 1.51 1255.7 1.24 1303.1 5.05 56.0
50 30 1257.1 1.35 1252.2 0.96 1280.8 3.27 65.3
Table 8: Experiment 5, Q=100, snull= 2330.0.
n type C % average δ1 best δ2 worst δ3 average t
100 10 2461.4 5.64 2427.7 4.19 2496.6 7.15 537.5
100 20 2487.1 6.74 2443.4 4.87 2580.2 10.72 481.4
100 30 2450.7 5.18 2433.4 4.43 2507.7 7.62 500.2
optimal solutions to the single-depot GPDP. However, from the results of our lim-ited experiments, its performance is highly problem dependent, on average, the solutions it produces are only a few percent worse than the best solutions we used to generate our test instances.
Table 9: Experiment 6, Q=100, snull= 2330.0.
n type C % average δ1 best δ2 worst δ3 average t
100 10 2433.3 4.43 2419.7 3.85 2461.3 5.63 501.4
100 20 2506.6 7.58 2485.4 6.67 2553.2 9.57 533.2
100 30 2500.7 7.33 2490.4 6.88 2511.4 7.78 614.7
Table 10: Experiment 7, Q=200, snull = 1885.8.
n type C % average δ1 best δ2 worst δ3 average t
100 10 1905.0 1.01 1902.0 0.86 1914.5 1.52 561.2
100 20 1963.4 4.11 1955.4 3.70 1988.3 5.43 557.6
100 30 1991.5 5.60 1987.2 5.37 2002.7 6.20 562.9
Table 11: Experiment 8, Q=200, snull = 1885.8.
n type C % average δ1 best δ2 worst δ3 average t
100 10 1905.7 1.05 1899.6 0.73 1935.0 2.60 571.2
100 20 1974.8 4.71 1971.5 4.54 2000.3 6.07 665.7
100 30 1982.3 5.11 1978.5 4.91 2011.1 6.64 652.3
(a) The null solution without
(b) With 30% clients having type C requests
Figure 5: Comparing sull and the modified problem with 30% type C requests.
References
2007.
[2] Birattari, M., Paquete, L., St¨utzle, T., Varrentrapp, K. “Classification of metaheuristics and design of experiments for the analysis of components,” Technical report at AIDA-01-05. 2001.
[3] Bianchessi, N., Righini, G. “ Heuristic algorithms for the vehicle routing prob-lem with simultaneous pick-up and delivery,” Computers and Operations Re-search, 34, pp. 578–594. 2007.
[4] Buijs, P., Roodbergen,K.J., Veenstra,M.,Lopez Alvarez, J.A. “ Improved col-laborative transport planning at Dutch logistics service provider Fritom,” Uni-versity of Groningen (not yet published) 2014.
[5] Br¨aysy, O., Gendreau, M. “Vehicle routing problem with time windows. part I: route construction and local search algorithms,” Transportation Science, 39, pp. 104–118. 2005.
[6] Dantzig, G.B., Ramser, J.H. “The truck dispatching problem,” Management Science, 6 , pp. 80–91. 1959.
[7] Gendreau, M.,Potvin, J.,Br¨asy, O.,Hasle, G.,Lokketangen, A. “Metaheuristics for the vehicle routing problem and its extensions: a categorized bibliogra-phy,” 2007.
[8] Gribkovskaia, I.,O. Halskau, G. Laporte, and M. Vlcek, “General solutions to the single vehicle routing problem with pickups and deliveries,” European Journal of Operational Research, 180, pp. 568–584. 2007.
[9] Groer, C., Golden, B., Wasil, E. “A library of local search heuristics for the vehicle routing problem,” Mathematical Programming Computation, 2, pp. 97–101. 2010.
[10] Hansen, P., Mladenovi´c, N. “Variable neighborhood search, a chapter of Hand-book of Metaheuristics,” Search Methodologies, pp. 211–238. 2005.
[12] Lau, H.C and Liang, Z. “Pickup and delivery problem with time windows, algorithms and test case generation,” International Journal on Artificial In-telligence Tools, 11, pp. 455–462. 2002.
[13] Montan´e, F., Galv˜ao, R. “A tabu search algorithm for the vehicle routing problem with simultaneous pick-up and delivery service,” Computers & Op-erations Research, 33, pp.595–619. 2006.
[14] Min, H. “The multiple vehicle routing problem with simultaneous delivery and pickup points,” Transportation Research, Part A: General, 23, pp.377– 386. 1989.
[15] Parragh, S., Doerner, K., Hartl, R. “ A survey on pickup and delivery prob-lems Part I: Transportation between customers and depot,” Journal f¨ur Be-triebswirtschaft, 58, pp. 21–51. 2008.
[16] Rego, C., Roucaiirol, C. “A parallel tabu search algorithm using ejection chains for the vehicle routing problem,” Metaheuristics Theory and Applica-tions, pp. 253–295. 1996.
[17] Savelsbergh, M., Sol, M. “The general pickup and delivery problem,” Trans-port Science, 29, pp. 17–29. 1995.
[18] Subramanian, A., Ochi, L.S., Cabral, L.A.F. “ An efficient ILS heuristic for the vehicle routing problem with simultaneous pickup and delivery,” Comput-ers and Operations Research, 37, pp. 1899-1911. 2010.
[19] Toth, P., Vigo, D. “An overview of vehicle routing problems,” The vehicle routing problem. SIAM monographs on discrete mathematics and applications 2002.
