A new hybrid ant colony optimization algorithm for the vehicle routing problem

Pattern Recognition Letters 30 (2009) 848–855

Contents lists available at ScienceDirect

Pattern Recognition Letters journal homepage: www.elsevier.com/locate/patrec

A new hybrid ant colony optimization algorithm for the vehicle routing problem Xiaoxia Zhang, Lixin Tang * The Logistics Institute, Northeastern University, Shenyang, China

a r t i c l e

i n f o

Article history: Available online 12 June 2008 Keywords: Scatter search Ant colony optimization Cyclic transfer Hybrid ant colony optimization algorithm

a b s t r a c t This paper presents a novel hybrid ant colony optimization approach called SS_ACO algorithm to solve the vehicle routing problem. The main feature of the hybrid algorithm is to hybridize the solution construction mechanism of the ant colony optimization (ACO) with scatter search (SS). In our hybrid algorithm, we use ACO algorithm and greedy heuristic to generate the initial solutions which are then formed the reference set. Within the scatter search framework, after two-solution combination method for the reference set has been applied, we employ ACO method to generate new solutions through updating the common arc pheromone mechanism. Moreover, during implementing the hybrid algorithm, cyclic transfers, a new class of neighborhood search algorithm can also be embedded into the scatter search framework as neighborhood search to improve solutions. Despite the size of the cyclic transfer neighborhood is very large, a restricted subset of the cyclic transfer neighborhood is adopted to reduce the computational requirements to reasonable levels. Finally, the experimental results have shown that the proposed hybrid method is competitive to solve the vehicle routing problem compared with the best existing methods in terms of solution quality. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction The vehicle routing problem (VRP) has been an important problem in the field of distribution and logistics. The VRP can be defined as a complete graph G = (V, A) where V = {0, . . ., n} is a vertex set, and A = {(i, j)ji, j 2 V} is an edge set. Vertex 0 is a central depot, while each other vertex represents a customer with a known nonnegative demand to be delivered. The distance dij is associated with each edge (i, j) 2 A and represents the distance from customer i to customer j. Customers geographically scattered are visited by a homogeneous fleet of vehicles with limited capacity Q and initially located at the central depot, and routes are assumed to start and end at the depot. The objective is to minimize total traveled distance and the number of vehicles, such that each customer is visited only once by a single vehicle, and vehicle capacity and the total distance must not be violated. The VRP is clearly NP-hard combinatorial optimization problem and difficult to solve. There have been important advances in the development of exact and approximate algorithms. Exact solution methods can only be used for very small instances, so for real-world problems, researchers have to rely on and resort to approximate or heuristic methods in solving the problem. Most published research for the VRP has concentrated on the development of heuristics. The majority of the research has emphasized particularly on designing efficient and effective heuristics. * Corresponding author. Tel./fax: +86 24 83680169. E-mail address: [email protected] (L. Tang). 0167-8655/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.patrec.2008.06.001

Earliest constructive heuristics for the VRP build routes successively adding an unvisited customer at each step according to cost-saving rule, as represented by Clarke and Wright (1964). These heuristics can generate feasible solutions very quickly, and these solutions can be easily improved by 2-opt, 3-opt procedures. Most recent research has focused on employing the advanced meta-heuristics such as simulated annealing and tabu search. Taillard (1993) and Rochat and Taillard (1995) have obtained the best known results to benchmark VRPs using the tabu search implementations while Osman (1993) have reported similar results obtained using simulated annealing. However, to obtain a high quality of solutions, these meta-heuristic algorithms usually consume substantial computing times and several parameter settings (Rego, 1998). Ant colony optimization (ACO) is a relatively new meta-heuristic and a very effective local search algorithm for a large number of combinatorial optimization problems. It simulates the behavior of ant colonies in nature as they forage for food and find the most efficient routes from their nests to food sources. The first ACO algorithm was developed by Clolrni et al. (1991) and successfully applied to the traveling salesman problem (TSP) based on the path-finding abilities of real ants (Dorigo and Gambardella, 1997). Bulleneimer et al. (1997) have designed the first ant system for the vehicle routing problem and then proposed an improved ant system algorithm (Bulleneimer et al., 1999). Because the vehicle routing problem is very complicated, the basic ACO algorithms cannot directly apply to the problem with satisfaction performance need. Many researchers have proposed new methods to improve the original ACO and applied them successfully to a whole range

X. Zhang, L. Tang / Pattern Recognition Letters 30 (2009) 848–855

of different problems. It is a trend to combine ACO with other algorithms to solve very large-scale of the vehicle routing problem. The development of modern meta-heuristics has led to considerable progress, but each meta-heuristic has its own strength and weakness. Therefore, much research has tried to develop the quest for the performance of hybrid algorithms expecting to achieve the effectiveness and efficiency. The distinctive features of this paper differs from all the papers described above on following points: first, the main idea of this paper is to hybridize the solution construction mechanism of the ACO with scatter search (SS), an evolutionary method, which results in a novel approach that we call SS_ACO. Scatter search (SS), first introduced by Glover (1997), is a population-based meta-heuristic that combines solutions from a reference set to create improved solutions. Within the scatter search framework, after two-solution combination method for the reference set has been applied, we employ ACO method to generate new solutions through updating the common arc pheromone mechanism. Second, cyclic transfers, a new class of neighborhood search algorithm (Thompson and Orlin, 1989) should be embedded into the scatter search framework as neighborhood search to improve solutions. The cyclic transfer is based on the multi-exchange exchange neighborhood structure developed by Thompson and Orlin (1989). Third, despite the size of the cyclic transfer neighborhood is very large (Thompson and Psaraftis, 1993), a restricted subset of the cyclic transfer neighborhood is adopted to reduce the computational requirements to reasonable levels. Moreover, since some of the parameters in our approach can be adjusted dynamically, the hybrid algorithm will explore different parts of the solution space, and the method search is not trapped at the local optimum. Finally, the experimental results have shown that the SS_ACO algorithm is to be very efficient and competitive in terms of solution quality. The paper is organized as follows: in Section 2, we review the basics of ant colony optimization and scatter search algorithms. In Section 3, we describe SS_ACO algorithm to tackle the vehicle routing problem. The speed-up search strategies are described in Section 4. The computational results are reported in Section 5. Finally, the conclusion is given in Section 6. 2. Ant colony optimization and scatter search Ant colony optimization (ACO) simulates the behavior of ant colonies in nature as they forage for food and find the most efficient routes from their nests to food sources. As some ants travel, they deposit a constant amount of pheromone trail that other ants are attracted to follow them. The increase in pheromone enhances the probability of the next ants selecting the path. Over time, as more ants are able to complete the shorter route, pheromone accumulates faster on shorter paths and longer paths are less reinforced. The ants are capable of not only finding the shortest path from a food source to the nest, but also adapting to changes in the environment once the old one is no longer feasible due to a new obstacle. This natural behavior of ants can be used to explain reason that they can find the shortest path. Bell and McMullen (2004) applies ant colony optimization to an established set of vehicle routing problems. The ACO method includes the two basic steps of construction of solutions and pheromone updating. Scatter search (SS) is a novel evolutionary method. It is a population-based meta-heuristic that combines solutions from a reference set to create improved solutions. The reference set is a set of feasible solutions and is usually updated by combining these existing solutions to obtain new ones. In the reference set, scatter search maintains some high quality solutions during the search process and also ensures diversity to explore other regions when the process is trapped in a local minimum. The main purpose behind scat-

849

ter search is that the newly combined solutions will explore various regions of the solution space where the better solution may possibly be found. Scatter search has regulations and strategies that are still not emulated by other evolutionary methods, and this method has been applied successfully to many combinatorial optimization problems. Marti et al. (2006) and Glover (1998) provide detailed descriptions and an implementation framework of the scatter search approach. Russell and Chiang (2006), and Greistorfer (2003) applied scatter search to the standard vehicle routing problem and combined solutions from the reference set through a common arc mechanism. The scatter search method is not restricted to a single uniform design, and it is very flexible and effective, since each of its elements can be implemented in variety of ways. 3. SS_ACO algorithm In this section, we describe our solution approaches to the vehicle routing problem. Because the vehicle routing problem is extremely complex, this makes the solution methodology more difficult. Quality solutions require an elaborate design of the algorithm. We proposed a SS_ACO algorithm, which is to hybridize the solution construction mechanism of ACO and scatter search (SS) to solve the vehicle routing problem. Meanwhile, cyclic transfer procedure is embedded into scatter search framework as neighborhood search to improve solutions by transferring customers among routes in a cyclic manner. The method has both the advantages of ant colony optimization, the ability to find the higher performance solutions, and that of scatter search, the ability to explore different parts of the solution space and to find better solutions. Furthermore, to search for an improvement neighbor, we propose an approximately dynamic programming to identify the most negative cost cycle on the improvement graph. The SS_ACO procedure steps should be summarized as follows: Step 1: Generate initial solutions. Initialize data and generate different initial solutions with the ACO algorithm and the greedy heuristic. Step 2: Build the reference set (RefSet). Select b1 best solutions using ACO algorithm and choose b2 best solutions using the greedy heuristic, and then put jRefSetj = b1 + b2 solutions in the reference set. Step 3: Find the best solution BestSol in the reference set. Step 4: Combine RefSet solutions. The common arc pheromone values are updated and new solutions are generated by two-solution combination. Step 5: Improvement solutions. After Nnum feasible solutions have been generated, the 2-opt heuristic and the cyclic transfer procedure are applied to improve solutions. Compute the objective function values of Nnum solutions and find high quality solutions. Step 6: Update the reference set and the best solution Bestsol. Search for the best solution in new reference set. The best solution Bestsol is updated if the improvement generates a better solution than BestSol. Step 7: Repeat Steps 4–6 until a terminating criterion is met. The terminating criterion is the maximum number of iterations. The procedure is repeated until this number is achieved. The framework of our proposed hybrid algorithm consists of the following sections. 3.1. Solution construction In the VRP ant colony optimization, the artificial ants simulate vehicles and successively select customers to visit, until all

850

X. Zhang, L. Tang / Pattern Recognition Letters 30 (2009) 848–855

customers have been visited. Initially, each ant starts at the central depot and the initial routes are empty. To select the next customer j, the ant k at the current position of customer i uses the following probabilistic formula:

j¼

8 < argmaxf½sil ½gil  eil b g; if q 6 q0 ; :

lRM k

S;

ð1Þ

otherwise;

where sil is equal to the amount of pheromone trail, b is a parameter which determines the relative importance of pheromone versus distance. The value q is a random uniform variable in [0, 1], and the value q0 (0 6 q0 6 1) is a parameter which determines the relative importance of exploitation versus exploration. The visibility gil is defined as the inverse of the distance dil and the savings eil of combining any two customers on one tour as opposed to serving them on two different tours are computed as eil = di0 + d0l  dil. Customers already visited by an ant k are stored in the ants working memory Mk and are not allowed to choose. If q 6 q0 then the best edge is chosen (exploitation) depending on Eq. (1); otherwise an edge is selected according to S (biased exploration). S represents a random variable selected according to the probabilistic distribution as the following equation:

P kij ¼

8 b < P sij ½gij eij  s

:

½g e  lRMk il il il

0;

b

; j R Mk ;

ð2Þ

otherwise;

3.2. Pheromone trail updating The pheromone trail is updated both locally and globally. Local updating is performed during solutions construction while global updating is performed at the end of the constructive phase. Local updating is intended to avoid a very strong edge being selected by all the ants, and it is accomplished according to the following rule:

ð3Þ

where 0 < c < 1 is a local pheromone decay parameter, L0 is the length of the initial solution produced by one of the construction heuristics and n is the number of customers. On the other side, global updating is used to intensify the search in the neighborhood of the best solution computed. In original ACO, only the best solution is used to globally modify the pheromone trail. Global trail updating is performed according to the following relationship:

sij ¼ ð1  qÞ  sij þ q=Lbest ;

ð4Þ

where 0 < q < 1 is a global pheromone decay parameter and Lbest is the length of the shortest path generated by ants since the beginning of the computation. Moreover, in addition to the local updating and global updating, when the two-solution combination method is used to generate new solutions, we also use a temporary rule to update common arc pheromone values according to the following relationship:

sij ¼ ð1  a0 Þ  sij þ a0 s0

3.3. A reference set generation method A feasible solution consists of r least-cost routes such that each customer is visited only once by a single vehicle, and vehicle capacity and the total distance must not be violated to minimize total traveled distance. The reference set is a set of feasible solutions, so we use two different methods to generate the initial solution set. The first one is ACO algorithm, and the second one is a greedy heuristic. The greedy heuristic is similar to the nearest neighbor heuristic for determining the sequence of a route is near at hand. The greedy heuristic adopts a parallel strategy here to simultaneously schedule multiple routes from a global optimal view. These procedures continue until the initial solution set is generated. The reference set consists of a total of b1 best solutions using ACO algorithm and b2 best solutions chosen using the greedy heuristic. Scatter search does not allow duplications in the reference set. The initialization reference set phase guarantees that various regions of the solution space will be explored. 3.4. A solution combination method

where pkij gives the probability with which ant k chooses to move from customeri to customer j. The solution construction process in ACO algorithm is equivalent to a greedy heuristic, except that the choice of next customers at each step is done probabilistically instead of deterministically. If we only consider the visibility gil, the inverse of the distance between customers, the objective function value may be not optimal. Our experimental results have shown that the inclusion of savings can lead to better results.

sij ¼ ð1  cÞ  sij þ c=ðn  L0 Þ

where a0 (0 < a0 6 1) and s0 is parameters. Common arc pheromone values should be stored using temporary variables and only take effect in each two-solution combination process so that the new solution cannot influence the other solutions. This updating strategy has been proved to be more efficient.

ð5Þ

New solutions are generated by a solution combination method from in the reference set. We adopt two-solution combination method to generated new solutions. To combine solutions from the reference set, common arcs are generated (Kelly and Xu, 1999), which should increase a probability of being chosen (Corberan et al., 2002). Let the reference set RefSet = {s1, s2, s3, . . . , sjRefSetj}. Two-solution combination subsets consist of {s1, s2}, {s1, s3}, . . . , {s2, s3}, . . . , {s—RefSet—1, s—RefSet—}. Each subset can generate a feasible solution based on these common arc pheromones at most, so the total number of solutions all subsets generate is equal to   jRefSetj . The combination method starts by matching the routes 2 from one solution to the routes of the other solution. The core behind route matching is to find the common arcs of the two solutions. We start from the principle that, if a solution includes common arcs that belong to a better solution, then the objective function value is probably better than that of a solution that does not contain such arcs. Since the shorter paths can deposit higher pheromone, the probability of ants following these shorter paths would be higher than that of those following the longer ones. The combination procedure is applied to all pairs of solutions in the current reference set Refset. For example, we would like to give the process combining the following two solutions (A and B) with 10 customers and three routes of combing solutions procedure (see Fig. 1). We suppose that Ai is the ith route for solution A and Bi is the ith route for solution B:

A1 ¼ fd01 ; d13 ; d30 g; A2 ¼ fd02 ; d25 ; d54 ; d40 g; A3 ¼ fd07 ; d76 ; d68 ; d89 ; d90 g: B1 ¼ fd01 ; d13 ; d32 ; d20 g;

B2 ¼ fd05 ; d54 ; d47 ; d70 g;

B3 ¼ fd06 ; d68 ; d89 ; d90 g: After the solutions are matched, the common arcs (i.e., \ Bj}, H = {(i,j)jAi \ Bj – /, i = 1, 2, 3, j = 1, 2, 3}) can be generated between routes. Pheromone trails accumulate faster on shorter paths and longer paths are less reinforced, so there would be the more probability of the shorter paths if ants follow these common arcs with combination method. After two solutions being S

(i,j)2H{Ai

X. Zhang, L. Tang / Pattern Recognition Letters 30 (2009) 848–855

5

3

2

5

4

3.5. The improvement method

3

2

4 1

0 7 9

6

1

0

7

9

6

8

8

a. Solution 1

5

2

3

b. Solution 2

4 1

0

7

9

6

851

8 c. Common arcs Fig. 1. Common arcs generated with combination method.

combined, one pheromone update strategy should be adopted according to modifying all the common arc pheromone mechanism (see Fig. 1c). We update the common arc pheromone values according to Eq. (5). Let mat(si, sj) be the number of the common edges between si and sj solutions in the reference set. The bigger the value of the mat(si, sj) is, the more similar these two solutions are. After the common arc pheromone values are updated via the pheromone update strategy, the ACO algorithm can find better solutions based on these common arc pheromones. In a combination subset, the ACO algorithm can generate a feasible solution based on these common arc pheromones at most. In special cases, if there are a few common arcs, i.e., with minimum mat(si, sj) in a combination subset, the combination subset may produce a diverse solution. In order that a new solution generated by a combination subset cannot influence the other solutions, common arc pheromone values should be stored using temporary variable and only take effect in the combination subset. Based on these common arc pheromone values, as far as the ACO algorithm process goes, at the end, Nnum solutions are generated from two-solution combination subsets.

After Nnum solutions generated, we attempt to improve the solution quality by applying improvement strategies. The 2-opt heuristic is applied at all times to act as an improvement heuristic within the route, while the cyclic transfer procedure is used for the between-route improvement. The central concept behind cyclic transfers, as applied to the VRP, is the attempt to find improved solutions by simultaneously transferring small customers among routes in a cyclic manner. Let I be a feasible solution to the vehicle routing problem, Ik be the set of customers assigned to vehicle k. A cyclic transfer with respect to I is a customer sequence W = {j1, j2, . . . , jm, j1} in which each customer belongs to distinct routes. The cyclic transfer W denotes the following simultaneous transfer of customers: customer j1 moves from I1 to I2, customer j2 moves from I2 to I3, and so on until customer jm moves from Im to I1. If I0 represents the obtained solution, then cost of the transfer is F(I0 )  F(I), and W is defined as an improvement cyclic transfer if F(I0 )  F(I) < 0. Note that r denotes the total number of vehicles and 2 6 m 6 r. Thompson and Orlin (1989) have found improvement moves in the cyclic transfer neighborhood by finding a negative cost cyclic transfer search, and proved that the search for negative cost cyclic transfers transforms into a search for negative cost cycles on the improvement graph. In order to construct a cyclic transfer, the improvement graph related to I is the auxiliary graph Gk = (Vk, Ek, Ck) defined as

V k ¼ fsets of customers in the same route k; k ¼ 1; 2; . . . ; mg; Ek ¼ fði; jÞ :; i 2 V k ; j 2 V kþ1 ; k ¼ 1; 2; . . . ; m  1g; C k ¼ fcki;j : ði; jÞ 2 EK ; k ¼ 1; 2; . . . ; mg are edge cost set: Let Ik1(i) and Ik(j) be the different routes, in which customers i 2 Ik1(i) and j 2 Ik(j). Let cki;j be the cost of each transfer (i, j). The cost cki;j of is equal to the corresponding changes in the cost of route Ik(j), owing to simultaneously adding customer i to Ik(j) and removing customer j from Ik(j), i.e., cki;j ¼ FðIk ðjÞ þ i  jÞ  FðIk ðjÞÞ. We propose an approximately dynamic programming to search for an improvement neighbor in the restricted cyclic transfer neighborhood. The iteration relationship of the best transfers is

Procedure cyclic transfer () { Initialize routes I(i) ,I(j) , the cycle length, k, i, j while (terminating criterion is not met) { repeat { Choose i from I(i), j node from I(j) If (the current route I(j) satisfies capacity constraints, or capacity constraints and the route length restrictions after node i replaces every node j) then Compute cost cik, j , fk (i, j) } end if until every node in route I(j) has been computed. Compute and record the corresponding nodes. I(i)=I(j),i=j,j=j+1, k=k+1. Update I(i) and I(j).} } if (the minimum negative cost is found in the last route) then If (there exists the negative cost cycle) then Backtrack to all the nodes constituting a cycle and record it else Increase the cycle length end if else Delete the cost cik, j end if end while Output the improved solution Fig. 2. Cyclic transfer procedure.

852

X. Zhang, L. Tang / Pattern Recognition Letters 30 (2009) 848–855

f1 ði; jÞ ¼ c1ij ; fk ði; jÞ ¼ min ffk1 ðp; iÞg þ ðp;iÞ2Ek1

ð6Þ ckij ;

k ¼ 2; . . . ; m:

ð7Þ

It is the key to identify the most negative cost cycle using approximately dynamic programming. If there exists the negative cost cycle, exchange customers forming the improved solution, backtrack to all the customers constituting a cycle and record it. Otherwise, we attempt to simultaneously reduce the cycle’s cost by modifying its length. If the improved solutions are better than some Nnum solutions, replace these solutions with the improved solutions. This procedure will be iteratively performed for a given time to improve the solution. Fig. 2 shows the cyclic transfer procedure. 3.6. A reference set update method The reference set is constructed by considering the different solutions obtained using two-solution combination method, the current reference set and the general ACO method. The best b solutions should be chosen to update the reference set. The new reference set consists of the best [b/2]  1 solution from the set of Nnum combined–improved solutions, [b/2]  1 solutions from the current reference set, a solution with minimum mat(si, sj) from Nnum combined–improved solutions and a solution generated with the general ACO to extend the diversity of the solutions. The solutions in the reference set are sorted by increasing their objective function values, where the best solution is the first one in the list. 4. Speed up the search strategy In general, the two-exchange neighborhood has O(n2) neighbors, whereas, for fixed value of m, there are O(nm) cyclic neighbors (Ahuja et al., 2000). For a VRP with r vehicles, the neighborhood size of the cyclic transfer involving m distinct vehicles is r  nm . The first term is the number of cyclic permutaðm  1Þ! m tion of m vehicles; the second one is the combination number of m vehicles chosen from a set of r vehicles; and the third one is the number of combinations of one customer from m vehicles. It is obviously shown that the cyclic transfer neighborhood is generally very large, and it is unlikely that one can search the neighborhood completely in polynomial time. Nevertheless, we can identify carefully selected special customers and search this restricted neighborhood to enhance the algorithm without worsening the time needed to find negative cost cycles. To speed up the algorithm, we should use a scheme to approximate the neighborhood search problem, and this scheme can search only a restricted subset of the cyclic transfer neighborhood. We do this in three ways. The first one is we limit the initial search to small negative costs cycle. We experimentally show there is no probability of finding negative cost cycles, if the value of cycle length is very higher. So we mainly adopt 2-cyclic 1-transfer or 3-cyclic 1-transfer in the cyclic transfer procedure, and then we use a variable depth approach to attempt to increase the cycle length. The second one is the sequence of routes is fixed, that is, we only consider transferring adjacent moves of customers that are belonged to adjacently routes in the cyclic transfer. For example, assume the vehicle sequence is fixed as {I1, I2, . . . , Im} (2 6 m 6 M) and customers can only be permitted to move from Ik1 to Ik (2 6 k 6 m) or from Im to I1. If the sequence of routes is fixed and m is small,   the size of the cyclic transfer r  nm . For a  small m, such as 3neighborhood is reduced to m r cyclic 1-transfer, this quantity reduces to  n3 . The third one 3 is the use of a candidate list for selecting the next customer in each vehicle route. Most of the good solutions have relatively shorter paths, and it is very rare for good solutions to connect very long

paths. If we can omit calculations for unimportant customers about long paths, we can pay more attention to important customers for finding better solutions. Each customer has a candidate list based on the distance to all other customers. Only the closest customers are assigned the candidate list for the current customer and are made available for the next customer to choose. The size of the candidate list can be given by restricting its size to a fraction of the number of customers. It is believed that this restriction can reduce calculations for wasting attempt considering moves to customers that are a very long distance from the current customer. Moreover, the cyclic transfer procedure is embedded into scatter search process and is executed to improve these solutions by transferring customers among routes in a cyclic manner. To save computation time, it is not necessary to call the cyclic transfer procedure at every solution or iteration. We restrict the certain frequencies of using the cyclic transfer procedure such as at every 10 iterations, or the cyclic transfer procedure is executed only if high quality solutions are obtained. The combination of these methods reduces the computational requirements to reasonable levels. For illustration, an example will be described to explain process of the cyclic transfer. Consider a 13-customer problem with customer positions (xi, yi) and demand qi shown in Table 1. Let {I1, I2, I3, I4} be a feasible solution, with I1 = {1, 2, 3}, I2 = {4, 5, 6}, I3 = {7, 8, 9, 10}, I4 = {11, 12, 13} and a vehicle capacity Q = 200. The feasible objective value is 183. Assume the vehicle sequence for dynamic programming is (I1, I2, I3, I4). The whole process of the cyclic transfer is described as follows: Step 1: Calculate the set of exchange customers for each vehicle. After exchange customers, vehicle capacity and the total distance must not be violated. On the first stage (k = 1), the set of exchange customers is E1 = {(1, 5), (1, 6), (2, 4), (2, 5), (2, 6)} from I1 to I2. Similarly, we have

E2 ¼ fð5; 7Þ; ð5; 9Þ; ð6; 7Þ; ð6; 8Þ; ð6; 9Þ; ð6; 10Þg; E3 ¼ fð7; 11Þ; ð7; 12Þ; ð8; 11Þ; ð8; 12Þ; ð9; 11Þð9; 12Þg: On the fourth stage (k = 4), the set of exchange customers is from I4 and I1.

E4 ¼ fð11; 1Þ; ð11; 2Þ; ð11; 3Þ; ð12; 1Þ; ð12; 2Þð12; 3Þg: Step 2: Compute cki;j of all in Ek for each vehicle k independently. By cki;j ¼ FðIk ðjÞ þ i  jÞ  FðIk ðjÞÞ, compute C1.

c11;5 ¼ 7;

c11;6 ¼ 4;

c12;6 ¼ 11;

c12;4 ¼ 9;

c12;5 ¼ 0;

C 1 ¼ f7; 4; 9; 0; 11g:

ð8Þ

Similarly, we have

C 2 ¼ f5; 6; 6; 5; 2; 3g;

C 3 ¼ f5; 0; 21; 1; 23; 1g;

C 4 ¼ f27; 29; 24; 15; 18; 10g: Step 3: By Eq. (6), compute f1(i, j).

f1 ð1; 5Þ ¼ 7;

f 1 ð1; 6Þ ¼ 4;

f 1 ð2; 5Þ ¼ 0;

f 1 ð2; 6Þ ¼ 11:

f 1 ð2; 4Þ ¼ 9;

Step 4: By Eq. (7), compute f2(i, j).

f2 ð5; 7Þ ¼ c25;7 þ minff1 ð1; 5Þ; f1 ð2; 5Þg ¼ 5; f2 ð5; 9Þ ¼ 6;

f 2 ð6; 7Þ ¼ 5;

f2 ð6; 8Þ ¼ 6;

f 2 ð6; 9Þ ¼ 9;

Step 5: f3 ð7; 11Þ ¼ 5,

c37;11

f 2 ð6; 10Þ ¼ 8:

þ minff2 ð5; 7Þ; f2 ð6; 7Þg ¼ 10; f 3 ð7; 12Þ ¼

f3 ð8; 11Þ ¼ 27; f 3 ð9; 12Þ ¼ 10:

f 3 ð8; 12Þ ¼ 7; f3 ð9; 11Þ ¼ 32;

853

X. Zhang, L. Tang / Pattern Recognition Letters 30 (2009) 848–855 Table 1 Data for the example problem Customer

1

2

3

4

5

6

7

8

9

10

11

12

13

xi yi qi

0 12 48

6 5 36

7 15 43

9 12 92

15 3 57

20 0 16

17 4 56

7 4 30

1 6 57

15 6 47

20 7 91

7 9 55

2 15 38

Step 5: f4 ð11; 1Þ ¼ c411;1 þ minff3 ð7; 11Þ; f3 ð8; 11Þ; f3 ð9; 11Þg ¼ 5,

f4 ð11; 2Þ ¼ 3; f4 ð11; 3Þ ¼ 8; f4 ð12; 1Þ ¼ 5; f4 ð12; 2Þ ¼ 8; f4 ð12; 3Þ ¼ 0: Step 6: Backtrack to all the customer sets constituting negative cost cycle. There exist negative cost value, f4(11,1) = 5, f4(11,2) = 3, f4(11,3) = 8. The customer sets constituting negative cost are {3,11,9,6,2}, {1,11,9,6,2}, {2,11,9,6,2}. Transfers (2,11), (11,9), (9,6) and (6,2) can form a negative cost cycle {2-11-9-6-2}. The improved solution is {(1, 11, 3), (4, 5, 2), (7, 8, 6, 10), (9, 12, 13)}, and the distance value after performing the cyclic transfer is 180. From this, we can see that some transfers such as (3,11), (11,9), (9,6), (6,2) and (1,11), (11,9), (9,6), (6,2) can cause the most negative transfer cost, but they cannot form a cycle. Fig. 3 shows a 4cyclic transfer for the VRP. As can be seen from the steps of the SS_ACO algorithm described above in the beginning of Section 3, the complexity of the SS_ACO algorithm is analyzed as follows. First, (Step 1), SS_ACO algorithm starts by calling the ACO algorithm and the greedy heuristic to generate initial solutions. The greedy heuristic is similar to the nearest neighbor heuristic and its complexity is of O(n2  r). The complexity of the ACO algorithm is O(n2  r). Moreover, (Step 5), the complexity of the 2-opt heuristic is O(n2). For a VRP with r vehicles,

I1

1

I1

3

2

1

I2

I1

3

11

I2

I1

4

9

12

5

12

5

13

6

13

2

11

I3

7

8

9

I3

10

4

7

8

6

10

Fig. 3. The effect of a 4-cyclic 1-transfer.

each serving n customers, of cyclic transfer proce complexity  the  r 2  m  n . So, the total operations tadure is evaluated to O m     r  m  n2 . ken by SS_ACO are bounded by O n2  r þ r  n2 þ m If m is a fixed value and it is a very small number,  theoverall com r  m  n2 . plexity of SS_ACO algorithm is equivalent to O m 5. Computational results In this section, we present computational results of our proposed algorithm, which was coded in the visual C++ and executed on an IBM computer with 512MB RAM and 1600Mhz CPU Speed. To evaluate validity of our proposed algorithm for the vehicle routing problem, the performance of our algorithms was tested on a set of 14 benchmark instances designed by Christofides et al. and can be downloaded from the OR Library at the website with URL:http:// mscmga.ms.ic.ac.uk/jeb/orlib/vrpinfo.html. The instances have 50– 199 customers in addition to the depot. The first 10 instances have customers that are randomly distributed around the depot. In the last four instances, the customers appear in clusters and the depot is not centered. All test instances are subjected to capacity constraints. Additionally, instances 6-10,13 and 14 have the route length restrictions as well. Many parameters exist in the SS_ACO algorithm, and the values of these affect directly or indirectly the final solution quality. We tested several values for each parameter while all the others were held constant. For these 14 test instances, the best known solution for C1 have been proven optimal by Fisher (1994). Therefore, most of the parameter values have been determined on the C1 instance by the numerical experiments. To better understand how b parameter affects the efficiency of the solution process. Different b values were tested in simulation. Fig. 4 shows that SS_ACO with parameter b = 3 converges faster than with b = 1. That is, a better solution can be obtained in the same number of or fewer iterations. Fig. 5 shows the number of optimal solution after 500 iterations with different b values. Simulation results reveal that the parameter b 2 {3,4,5,6} can yield more optimal solutions.

Best length Best length

620

620

600

600

580

580

560

560

540

540

520 Iteration 500

50

125

175

225

275

a. β =1

325

375

425

475

520 500

Iteration 50

125

175

225

275

b. β = 3

Fig. 4. The performances of different b parameters for the C1 problem.

325

375

425

475

854

X. Zhang, L. Tang / Pattern Recognition Letters 30 (2009) 848–855 Table 2 The comparison of SS_ACO with ACO algorithm

Number of optimal solutions

284

Prob.

n

m

Q

D

d

ACO

SS_ACO

% Gap

282

C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 C13 C14

50 75 100 150 199 50 75 100 150 199 120 100 120 100

5 10 8 12 17 6 11 9 14 18 7 10 11 11

160 140 200 200 200 160 140 200 200 200 200 200 200 200

1 1 1 1 1 200 160 230 200 200 1 1 720 1040

0 0 0 0 0 10 10 10 10 10 0 0 50 90

533.15 860.26 840.51 1073.00 1350.24 566.00 927.00 889.45 1198.12 1451.00 1073.00 846.00 1613.00 890.00

524.61 835.26 830.14 1038.20 1307.18 559.12 912.68 869.34 1179.4 1410.26 1044.12 824.31 1556.52 870.26

1.60 2.91 1.23 3.24 3.18 1.21 1.54 2.26 1.56 2.81 2.69 2.56 3.50 2.22

280 278 276 274 272

β

270 1

2

3

4

5

6

7

8

Fig. 5. The number of optimal solution after 500 iterations with different b values.

To determine q0 values, a test is run 20 times for C1 instance, and the test is terminated until an optimum is found. Fig. 6 illustrates the average time to find the optimum solutions for each q0. The figure indicate that the setting of q0 does not affect more to the final results, but q0 2 [0.65, 0.90] takes least average time in finding the optimums. So the range of q0 is set q0 2 [0.65, 0.90] for the following experiments. In the same way, the q and c values were set in the [0.04, 0.06] range. As described above, experiments have been done with the following parameter settings: q0 2 [0.65, 0.90], b 2 {3,4,5,6} and c,q 2 [0.04, 0.06]. The default value of the parameters was q0 = 0.90, b = 6, q = 0.04, c=0.04. With some parameter settings we observed that, after several iterations, all the ants followed the same routes because of a much higher trail level on the edges, which made new possibilities that ants explored very low. Thus, some of the parameters in our approach are adjusted dynamically as follows: when SS_ACO algorithm produces the same routes during successive three iterations, it may be fall into a local optimum, and hence, the parameter q0 was decreased from 0.90 to 0.65 with a step by 0.05 and b from 6 to 3 with a step by 1 in order to explore different regions and escape from the local minimum. Parameter c and q were increased from 0.04 to 0.06 with a step by 0.005. Other settings were: jRefSetj = 10, s0 = 0.01. The experimental results have shown that a lowers0 value applied in SS_ACO is better than a high value. All the tests have been carried out for maximum iterations Nmax = 2000. Table 2 reports the comparison between SS_ACO and the general ACO algorithm without combination with other algorithms. It has shown that particular construction mechanism of ACO can find the better solutions, but SS_ACO performs consistently better than the general ACO. It is demonstrated that the hybridization technique we made to the ACO procedure could improve the effectiveness of the algorithm. The average improvement of the algorithm is about 2.32% in the distance. Since scatter search

Note: n, number of customers; m, number of vehicles; Q, vehicle capacity; D, maximum tour length; d, service time.

(SS) is a population-based meta-heuristic which takes much longer time to create improved solutions by combining solutions from a reference set, the computation speed of SS_ACO was slower than that of ACO algorithm, the results are presented in Table 4. We compare SS_ACO with a number of the better methods available for the VRP, and the results of some problems are described in Table 3, where SA refers to Simulated Annealing by Osman (1993), TS to Tabu Search by Osman, GA to genetic algorithm by Baker and Ayechew (2003), ACO to pure ACO, and SS_ACO is the algorithm we proposed. The SS_ACO algorithm has shown to be competitive with the best existing methods in terms of solution quality. Moreover, during this experiment, some solutions are as close as the best solution published so far. Table 3 also shows the gap between the best known solution and our hybrid algorithm. The gap is defined as the percentage increase in the objective function value relative to the best known. A zero gap indicates that best known solutions should be obtained. Overall, the hybrid ACO performs better, with distances averaging only 0.62% above the best known values. In comparisons to the best known solutions, SS_ACO can produce very good solutions for the seven instances with capacity constrains, and an average deviation is less than 0.5%. However, its performance for the seven route length-constrained instances is less impressive with an average deviation of 0.76%. We should remark that there exist some tabu search algorithms, such as Tailard’s algorithm (Taillard, 1993), which have found most of the best known solutions. However, these best solutions can be gotten on an unspecified number of runs with different settings, and no computation times has been mentioned by the authors.

Table 3 Comparison of heuristics for the vehicle routing problem

Average time (seconds) 34 33 32 31 30 29 28 27 0. 1

0. 2

0. 3

0. 4

0. 5

0. 6

0. 7

0. 8

0. 9

1

Fig. 6. The average time to find the optimum solutions for each q0.

q0

Prob.

Best published

SA

TS

GA

SS_ACO

% Gap

C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 C13 C14

524.61 835.26 826.14 1028.42 1291.45 555.43 909.68 865.94 1162.55 1395.85 1042.11 819.56 1541.14 866.37

528 838 829 1058 1376 555 909 866 1164 1418 1176 826 1545 890

524 844 835 1052 1354 555 913 866 1188 1422 1042 819 1547 866

524.81 849.77 840.72 1055.85 1378.73 560.29 914.13 872.82 1193.05 1483.06 1060.24 877.8 1562.25 872.34

524.61 835.26 830.14 1038.20 1307.18 559.12 912.68 869.34 1179.4 1410.26 1044.12 824.31 1556.52 870.26

0 0 0.48 0.95 1.22 0.66 0.33 0.39 1.45 1.03 0.19 0.58 1.00 0.45

X. Zhang, L. Tang / Pattern Recognition Letters 30 (2009) 848–855 Table 4 The comparison of CPU times (in seconds) with published results Prob.

SA

TS

GA

ACO

SS_ACO

C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 C13 C14

167.4 6434.3 9334 5012.3 2318.1 3410.2 626.5 957.2 84301.2 5708 315.8 632 7622.5 305.2

114 178.7 1543 3560 3246 173 1056.7 2998 4755.8 4561 1445.4 892.2 2834 1175.9

213 765 1148 2475 3999 217 786 1134 2258 3687 1633 1160 1694 1197

40.20 50.37 102.56 145.7 523.16 53.37 60.31 150.06 386.80 745.24 127.5 113.73 306.62 280.36

55.23 70.43 120.25 250.76 707.80 65.17 90.42 210.46 520.52 1012.23 232.46 156.47 467.24 368.72

Table 4 shows the computation times that have been reported by various algorithms. Among them, Osman (SA and TS) reported the CPU times to obtain the best solutions on a VAX 8600. The algorithm of GA was run on a Pentium 266 MHz. Despite the cyclic transfer neighborhood is a very large-scale, we adopt the speedup search strategy to search only a restricted subset of the cyclic transfer neighborhood. The experimental results have shown that the method is effective to reduce the computational requirements to reasonable levels in practice. 6. Conclusions In this paper, we propose a SS_ACO algorithm, a hybrid ACO, which hybridizes the solution construction mechanism of ACO with the scatter search. In this procedure, in order to obtain higher solution, 2-opt heuristic is used to improve within-route, while the cyclic transfer is selected for the between-route improvement. Although the size of the cyclic transfer neighborhood is very large, a restricted subset of the cyclic transfer neighborhood is adopted to reduce the computational requirements to reasonable levels. The experimental results have shown that the SS_ACO algorithm produces uniformly higher performance solutions relative to the other competing heuristics on the vehicle routing problem. Acknowledgements This research is partly supported by National Natural Science Foundation for Distinguished Young Scholars of China (Grant

855

No. 70425003), National 863 High-Tech Research and Development Program of China through approved No. 2006AA04Z174 and National Natural Science Foundation of China (Grant No. 60674084). References Ahuja, R.K., Orlin, J.B., Sharma, D., 2000. Very large-scale neighborhood search. Internat. Trans. Oper. Res. 7, 301–317. Baker, B.M., Ayechew, M.A., 2003. A genetic algorithm for the vehicle routing problem. Comput. Oper. Res. 30, 787–800. Bell, J.E., McMullen, P.R., 2004. Ant colony optimization techniques for the vehicle routing problem. Adv. Eng. Inform. 18, 41–48. Bulleneimer, B., Hartl, R.F., Strauss, C., 1997. Strauss, applying the ant system to the vehicle routing problem. In: Proc. of the Second Internat. Conf. on Metaheuristics, pp. 1–12. Bulleneimer, B., Hartl, R.F., Strauss, C., 1999. An improved ant system algorithm for the vehicle routing problem. Ann. Oper. Res. 89, 319–328. Clarke, G., Wright, J.W., 1964. Scheduling of vehicles from a central depot to a number of delivery points. Oper. Res. 12, 568–581. Clolrni, A., Dorigo, M., Maniezzo, V., 1991. Distributed optimization by ant colonies [A]. In: Proc. of the First European Conf. of Artificial Life (ECAL’91). Elsevier, pp. 134–142. Corberan, A., Fernandez, E., Laguna, M., Marti, R., 2002. Heuristic solutions to the problem of routing school buses with multiple objectives. J. Oper. Res. Soc. 53 (4), 427–435. Dorigo, M., Gambardella, L.M., 1997. Ant colonies for the traveling salesman problem. BioSystems 43, 73–81. Fisher, M.A., 1994. Optimal solution of vehicle routing problems using minimum ktrees. Oper. Res. 42, 626–642. Glover, F., 1997. Heuristics for integer programming using surrogate constraints. Decision Sci. 8, 156–166. Glover, F., 1998. A template for scatter search and path relinking. In: Hao, J.-K., Lutton, E., Ronald, E., Schoenauer, M., Snyers, D. (Eds.), Artificial Evolution. Lecture Notes in Computer Science, vol. 1363. Springer-Verlag, Berlin, Heidelberg, New York, pp. 13–54. Greistorfer, P., 2003. A tabu scatter search metaheuristic for the arc routing problem. Comput. Ind. Eng. 44 (2), 249–266. Kelly, J.P., Xu, J., 1999. A set-partitioning-based heuristic for the vehicle routing problem. INFORMS 11 (2), 161–172. Marti, M., Laguna, M., Glover, F., 2006. Principles of scatter search. Eur. J. Oper. Res. 169, 359–372. Osman, I.H., 1993. Metastrategy simulated annealing and Tabu search algorithms for the vehicle routing problem. Oper. Res. 41, 421–451. Rego, C., 1998. A subpath ejection method for the vehicle routing problem. Manage Sci. 44 (10), 1447–1458. Rochat, Y., Taillard, E.D., 1995. Probabilistic diversification and intensification in local search for vehicle routing. J. Heuristics 1, 147–167. Russell, R.A., Chiang, W.C., 2006. Scatter search for the vehicle routing problem with time windows. Eur. J. Oper. Res. 169, 606–622. Taillard, R.E., 1993. Parallel iterative search methods for vehicle routing problems. Networks 23, 661–673. Thompson, P.M., Orlin, J.B., 1989. The theory cyclic transfers. Working Paper, Operations Research Center, MIT. Thompson, P.M., Psaraftis, H.N., 1993. Cyclic transfers algorithms for multivehicle routing and scheduling problem. Oper. Res. 41 (5), 935–946.