Local search with annealing-like restarts to solve ... - Semantic Scholar
European Journal of Operational Research 150 (2003) 115–127 www.elsevier.com/locate/dsw
Discrete Optimization
Local search with annealing-like restarts to solve the VRPTW Haibing Li, Andrew Lim
*
Department of Computer Science, National University of Singapore, 3 Science Drive 2, Singapore 117543, Singapore Received 25 April 2001; accepted 3 May 2002
Abstract In this paper, we propose a metaheuristic based on annealing-like restarts to diversify and intensify local searches for solving the vehicle routing problem with time windows (VRPTW). Using the SolomonÕs benchmark instances for the problem, our method obtained seven new best results and equaled 19 other best results. Extensive comparisons indicate that our method is comparable to the best in published literature. This approach is flexible and can be extended to handle other variants of vehicle routing problems and other combinatorial optimization problems. Ó 2002 Elsevier Science B.V. All rights reserved. Keywords: Routing; Vehicle routing problem with time windows; Local search; Diversification; Intensification
1. Introduction Local searches are commonly used to solve hard combinatorial optimization problems. A local search algorithm starts with an initial solution and searches within neighborhoods for better solutions. Being trapped in local optima is a significant factor affecting the quality of solutions obtained by local searches. Advanced search algorithms, including simulated annealing (SA), tabu search (TS) and evolutionary algorithms (EA), were devised to help escaping from local optimal traps. These approaches are classified as metaheuristics. From the viewpoint of local search, one of the important issues in designing
* Corresponding author. Tel.: +65-874-6891; fax: +65-7794580. E-mail addresses: [email protected] (H. Li), alim@ comp.nus.edu.sg (A. Lim).
these algorithms lies in diversifying the search spaces using respective strategies, so that the whole solution space can be explored as widely as possible. Another important issue in designing local search algorithms is intensification. A local search algorithm is more efficient if it is smart enough to intensify the search on occasion, for checking if there are better solutions. With regards to the diversification and intensification issues, in this paper, we propose a tabu-embedded simulated annealing (TSA) metaheuristic which runs with a best-based K-restarts strategy during local searches within three constraint-based neighborhoods, for solving the vehicle routing problem with time windows (VRPTW). The sensitivity of K, which affects the solution quality and computational effort, is analyzed by experimentation. The rest of the paper is organized as follows. In Section 2, the notation is presented and then the VRPTW is briefly described with a review of
0377-2217/03/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. doi:10.1016/S0377-2217(02)00486-1
116
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
related work. Section 3 describes the important issues for designing our local search algorithms. In Section 4, we present the algorithm in detail. Computational experiments are conducted in Section 5, in addition, results on 56 SolomonÕs benchmark problem instances are reported and our approach is compared with the best methods for this problem. Finally, we present our conclusions in Section 6.
2. Notations and the VRPTW 2.1. Notations For the convenience of description, we define the following notation: GðV ; AÞ G is a digraph generated from problem data, V is the node set of G and A is the arc set of G n, m number of customers, number of vehicle used, respectively L the maximum number of consecutively non-improving iterations in the TSA procedure. It is proposed to be a small number K the maximum number of restarts of the TSA procedure on the current best solution CostðSÞ objective function of solution S of local search SACostðSÞ objective function of solution S in tabuembedded Metropolis procedure in the SA algorithm DistðSÞ, STðSÞ, WTðSÞ travel distance, schedule time and waiting time of solution S, respectively N) ðSÞ, N () ðSÞ, N,! ðSÞ constraint-based neighborhoods of solution S obtained by shift operator, exchange operator and rearrange operator respectively N ðSÞ N ðSÞ ¼ N) ðSÞ [ N () ðSÞ are the neighborhoods where local searches are conducted / the depth of K-ary tree due to K restarts T0 , T initial and global annealing temperature respectively
d cooling ratio of T a, b, c, k penalty factor for m, DistðSÞ, STðSÞ, WTðSÞ respectively the customer node at the position of pos vpos in a route Dpos1;pos2 the distance between customers at position pos1 and pos2 in a route 2.2. The VRPTW The VRPTW is a generalization of classical vehicle routing problems with the additional complexity of time window constraints. Let G ¼ ðV ; AÞ be a digraph. V ¼ VN [ fv0 g is the node set where VN ¼ fvi 2 V ji ¼ 1; 2; . . . ; ng represents the customers and node v0 denotes the depot where a fleet of vehicles are housed. Each node vi 2 V has an associated customer demand qi (q0 ¼ 0), a service time si (s0 ¼ 0) and a service-time window ½ei ; li . For each pair of nodes hvi ; vj i (i 6¼ j, i; j ¼ 0; 1; 2; . . . ; n), a non-negative distance dij and a non-negative travel time tij are known. Due to the time window constraints, arcs may not exist between some node pairs. Therefore, the arc set can be defined as A ¼ fhvi ; vj ijvi ; vj 2 V ; vi 6¼ vj ; t0i þ si þ tij 6 lj g. If a vehicle reaches a customer vi before ei , it needs to wait until ei in order to service the customer. The schedule time of a route is the sum of the waiting time, the service and travel time. The objective of the VRPTW is to service all customers without violating vehicle capacity constraints and time window constraints with a minimum number of vehicles and, for the same number of routes with the minimum travel distance, followed by the minimum schedule time and the minimum waiting time. Savelsbergh [23] proved that the VRPTW is NP-hard. Much work has been done to derive exact algorithms to solve the problem. Desrochers et al. [7] proposed a column generation approach that solved seven of SolomonÕs benchmark instances exactly proposed in [25]. Fisher et al. [10] proposed a K-tree relaxation approach that solved two of SolomonÕs benchmark instances. Other exact approaches can be found in [9,12,13]. Because of the exponential time complexity of these approaches, it is unlikely that these algorithms can produce optimal solutions for practical-sized
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
VRPTW with reasonable computational time. Hence, much work has focused on generating near optimal solutions using heuristic approaches. In general, heuristics for VRPTW can be classified as classical heuristics or metaheuristics. Classical heuristics include construction heuristics, improvement heuristics and composite heuristics. Route construction heuristics build feasible solutions by inserting an unrouted customer into a current partial route until all customers are routed. The sequential insertion algorithm proposed by Solomon [25] belongs to construction heuristics. Its parallel version was implemented in [17]. These methods often produce solutions very quickly, but the quality of these solutions may not be good. For this reason, construction heuristics are usually applied to obtain initial solutions for improvement heuristics or other two-phase heuristics. Route improvement heuristics modify a solution by performing local searches for better neighboring solutions within its neighborhoods generated by node/edge-swapping operators. Edge-exchange heuristics are referred to [18,23,24]. Efficient implementations of edge-exchange for speeding up the pruning of inferior solutions are reported in [23,24,26]. Russell [22] proposed a composite heuristic which mixes route construction and route improvement procedures. Metaheuristics are a new generation of heuristics developed in the recent decade. Near optimal solutions for the VRPTW were produced by TS [3,5,16,19,20,26], SA [2], EA [11,15,27] and hybrid algorithms of the above three algorithms [22,28]. Surveys on the VRPTW are referred to [4,6,8].
3. Issues for designing local search algorithms During designing the local search algorithms, we take several important issues into consideration. Firstly, the neighborhood generating mechanisms, which usually have a strong connection with the problems and the solution representation, provide spaces within which local searches are conducted. For example, for the vehicle routing problems, the solutions are naturally represented by the routes with sequences of customer nodes
117
which are visited by the vehicles to offer a service. So it is a natural idea to obtain neighboring solutions by swapping nodes/edges between routes. Secondly, it is also important to evaluate the quality of solutions during local searches. Different problems have different cost functions which faithfully reflect the problemsÕ objectives and act as solution evaluators. For problems with multiple objectives, penalty functions are commonly used for this purpose. Thirdly, the way to obtain and accept a neighboring solution, affects the exploring trail of local searches. This reflects an important issue in local searches, i.e. diversification. For example, a descent local search (DLS) systematically exploring the whole neighborhood of a seed solution, and discards infeasible/inferior solutions immediately after solution evaluation while it keeps better solutions as new seed solutions for repeating the search procedure. The procedure stops at a dead end when no more better solutions can be found within a seed solutionÕs neighborhood. This shortsighted strategy inevitably causes local searches to be trapped in local optima very quickly. To help jumping out of the dead ends of DLS, in addition to systematically explore the neighborhoods, TS intelligently goes a step further to label along its exploring trail and to accept unexploited local optima no matter whether they are better or worse than seed solutions. In contrast, SA obtains a neighboring solution randomly but accepts it with a probability reflecting annealing property. Therefore, both TS and SA diversify local searches to escape from local optima. However, TS does so in a deterministic way so that its exploring trail is fixed, while SA does so in a non-deterministic way and its exploring trails vary in different tries. In this work, we propose a TSA procedure to diversify local searches. Finally, another important local search issue is intensification. Both TS and SA do not deal specially with this issue. In past researches, postprocessing on the best solution found so far was used frequently for this purpose. However, it is a pity that this intensification was only used after local searches were finished, instead of being embedded in the local search procedure. In this work,
118
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
we propose a K-restart strategy which guides the TSA to search around the current best solution. If a new better solution is found within K restarts of TSA, the current best solution is updated and a new series of restarts are initiated. Otherwise, if K consecutive restarts are finished without finding a better solution, this procedure stops. Therefore, the intensification which is in favor of the current best solution always plays its role during local searches. Implementation details for these issues will be presented in Section 4.
4. The algorithm The algorithm is proposed as a two-phase framework: Phase 1: Obtain an initial solution for the VRPTW using a heuristic. Phase 2: Start local search from the initial solution, using TSA with K-restart strategy. In phase 1, we use the insertion heuristic proposed by Solomon [25] to produce initial solutions. The initial solutions are then used to feed the local search procedure in phase 2. Of course, the local search in phase 2 are of our main interest. We deal with them in details with regards to the issues discussed in Section 3.
change operator are used for generating neighborhoods for local searches in the TSA procedure, while the rearrange operator is only used for postprocessing on the local optima obtained by local searches. 4.1.1. The shift operator Using the shift operator, non-null route segments (with at least one node) can be moved from one route to another route under the condition of not violating the capacity constraints and time window constraints. For each pair of selected routes, say, Route i and Route j (we assume 1 6 i < j 6 m, i 6¼ j without loss of generality), the shift operator is used in two directions, so that a sub-neighborhood Ni)j ðSÞ can be obtained by shifting route segments from Route i to Route j and another sub-neighborhood Nj)i ðSÞ can be obtained by shifting route segments from Route j to Route i. We have denoted the neighborhood of solution S obtained by the shift operam1 tor as N) ðSÞ, thus, N) ðSÞ ¼ [i¼1 [mj¼iþ1 fNi)j ðSÞ[ Nj)i ðSÞg. Fig. 1 shows an example of shifting a route segment from Route 1 to Route 2. To accelerate the local searches, three simple conditions are considered to prune infeasible or inferior solutions in constant time, during local searches within N) ðSÞ. We use the example shown in Fig. 1 to explain this. Firstly, if the total-capacity of the route segment from pos11 to pos12 to be shifted in
4.1. Neighborhood generating mechanisms In our method, three edge-swapping operators are used to generate local search neighborhoods. These operators are, the shift operator, the exchange operator and the rearrange operator. Savelsbergh [23,24] and Taillard et al. [26] also dealt with similar edge-exchange operators for the VRPTW. But they restricted the length of edge segments to be swapped to certain maximum values. We do not impose such restrictions. Further, we use three simple conditions that can be calculated in constant time to prune infeasible/inferior solutions, so that the neighborhoodsÕ sizes are reduced to accelerate the local searches. Among these operators, the shift operator and the ex-
Fig. 1. The shift operator.
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
Route 1 exceeds the left-capacity of Route 2, the solution obtained by shifting is infeasible so that it can be discarded. We record the used-capacityso-far in each route node so that computing total-capacity and left-capacity can be done in constant time. Secondly, if either hvpos prev ; vpos11 i or hvpos12 ; vpos i is an infeasible arc in digraph G defined in Section 2, the shift operation also results in an infeasible solution that can be discarded. Since the arc set A is initialized when the graph data is generated, checking for this condition can be done in constant time. Finally, since we are only interested in searching for better solutions, if an edge-shift incurs an increment in the total travel distance of the two routes, then this shift operation results in an inferior solution that can be discarded too. To simplify the checking procedure such that if an increment in travel distance incurs, we can use the following condition, let: D ¼ ðDpos11 prev;pos12 next þ Dpos prev;pos11 þ Dpos12;pos Þ ðDpos11 prev;pos11 þ Dpos12;pos12 next þ Dpos prev;pos Þ ð1Þ Here, ðDpos11 prev;pos12 next þ Dpos prev;pos11 þ Dpos12;pos Þ is the total distance of the three newly connected edges, while ðDpos11 prev;pos11 þ Dpos12;pos12 next þ Dpos prev;pos Þ is the total distance of the three old edges removed due to shifting. D P 0 means that the travel distance will not be decreased, hence the solution obtained by the operation does not result in a better solution. However, as our first objective is to minimize the number of vehicles before excluding a shift operation under the condition of D P 0, we must make sure that this shift operation will move the whole route segment out of a route so that we can eliminate one vehicle. Therefore, we will only discard a solution if D P 0 and the shifting operation will not eliminate a vehicle. As determined by Eq. (1), this check can be performed in constant time. Using these three conditions, jN) ðSÞj can be reduced efficiently. We name the solution space N) ðSÞ constrained by the three conditions as a constraint-based neighborhood.
119
4.1.2. The exchange operator Using the exchange operator, non-null route segments can be exchanged between two routes under the condition of not violating capacity and time window constraints. We have denoted the neighborhood of solution S obtained by this operator as N () ðSÞ. Fig. 2 shows an example of exchanging two route segments between Route 1 to Route 2. Like the shift operator, infeasible or inferior solutions obtained by exchange operation can also be pruned during local searches using the same three conditions. Of course, since more arcs and route segments are involved in the exchange operator, the conditions are a little more complex. However, all three conditions can also be calculated in constant time. Actually, it is possible to obtain an exchange operation by two shift operations. However, as we discussed previously, each shift operation requires checking of three conditions to prune infeasible/ inferior solutions. Thus, two shift operations need to check six conditions. This results in additional computational time and implementation complexity in contrast to the exchange operator. 4.1.3. The rearrange operator Using the rearrange operator, route segments can be repositioned within the same route under the condition of not violating the time window constraints. We have denoted the neighborhood of solution S obtained by the rearrange operator as
Fig. 2. The exchange operator.
120
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
Fig. 3. The rearrange operator.
N,! ðSÞ. Fig. 3 shows an example of repositioning a route segment to another position within a route. Similarly, we can also get a constraint-based neighborhood N,! ðSÞ using the second condition and the third condition discussed above. The first condition is not needed because the capacity remains unchanged when the rearranging operation is applied within one route.
with four fields. Since the probability that two different solutions having same record values is very small, it is reasonable to regard two solutions as the same if their record values are not different. We store the records of the solutions in a tabu list which differs from some references that use edgemoves in their tabu structure. Firstly, the pseudocode of the TSA procedure is described as follows:
4.2. Objectives and cost function in local search The priority-order of the objectives for solving the VRPTW is: (1) to minimize the number of vehicles used; (2) to minimize the travel cost; (3) to minimize the schedule time; (4) to minimize the driversÕ total waiting time to begin service. To reflect this order of priority, we adopt the following cost function: CostðSÞ ¼ am þ bDistðSÞ þ cSTðSÞ þ kWTðSÞ ð2Þ The penalty weight factors are set to be: a b c k. Using this setting, even if the initial solutions obtained by the route construction heuristic contain large number of vehicles, the vehicle size will be reduced by local searches. In addition, we can change the setting of these factors to get different priority orders of the four objectives. 4.3. Tabu-embedded simulated annealing with Krestart strategy To introduce the TS into local searches, we record the four objective values in a data structure
ALGORITHM. TSA(Solution S) 1. Sb S; 2. NonImproving 0 3. WHILE NonImproving < L DO 3.1 Use TMPðSÞ to obtain a non-tabu neighbor S 0 2 N ðSÞ 3.2 Sb0 Perform TDLS within N ðS 0 Þ / local search 3.3 Sb0 Perform TDLS within N,! ðSb0 Þ / postprocessing 3.4 IF CostðSb0 Þ < CostðSb Þ THEN Sb Sb0 ; NonImproving 0 3.5 ELSE NonImproving NonImproving þ 1 3.6 S Sb0 4. RETURN Sb In the above procedure, TDLS is a tabuembedded DLS procedure, in which only nontabu neighbors are accepted. The procedure TMP is the tabu-embedded metropolis procedure reflecting the annealing property. It is described as follows:
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
ALGORITHM. TMP(Solution S) 1. WHILE NOT yet accept a neighbor S 0 2 N ðSÞ DO 1.1 Get a random non-tabu neighbor S 0 2 N ðSÞ 1.2 D SACostðS 0 Þ SACostðSÞ 1.3 IF D 6 0 THEN prob 1 1.4 ELSE prob eD=T 1.5 IF randomð0; 1Þ 6 prob THEN 1.5.1 Accept S 0 1.5.2 Update annealing temperature: T dT 1.5.3 Record S 0 into tabu list 2. RETURN S 0
121
Fig. 4. Diversification by TSA.
In the procedure, SACostðsÞ is a cost function to evaluate solutions on behalf of SA. It is defined to adjust the cooling procedure with setting of T0 , T and d. It can be defined as Eq. (3). SACostðSÞ ¼ DistðSÞ þ uWTðSÞ; where u ¼ 0:01 DistðSÞ
ð3Þ
Finally, the algorithm which combines TSA with the best-based K-restart strategy is as follows: INPUT: an initial solution S0 OUTPUT: a best solution Sg found so far ALGORITHM. K-TSA(Solution S0 ) 1. Configure parameters: K, L, d and T0 ; T T0 2. Set tabu list to be empty 3. Sg S0 ; gNonImproving 0 4. WHILE gNonImproving < K DO / K-restart 4.1 S Sg ; 4.2 Sg0 TSAðSÞ 4.3 IF CostðSg0 Þ < CostðSg Þ THEN / Update new current best solution Sg Sg0 ; gNonImproving 0 4.4 ELSE gNonImproving gNonImproving þ 1 5. OUTPUT Sb As shown in Fig. 4, TSA is in fact a depth-like search which starts from a neighbor in N ðSg Þ, and stretches the search away from N ðSg Þ by TMP and TDLS. Actually, only with this depth-like search way, can solutions outside of N ðSg Þ be explored in
Fig. 5. Intensification by K-restart strategy.
each restart of TSA on Sg . Thus, the local searches around Sg is diversified if L is not too small. On the other hand, the K-restart strategy guides TSA to start local searches within N ðSg Þ in a breadth-like way. It intensifies local searches in favor of N ðSg Þ. The search behavior of K-TSA is shown in Fig. 5. 4.4. Configuration of K, L, T0 and d First of all, K and L are the key factors that affect the diversification and intensification behavior. It would be preferable to have large value for both K and L, but the computational cost would be very expensive. The largest possible value of K is jN ðSg Þj. If jN ðSg Þj is not too large, it would be preferable to accept each neighbor of Sg systematically. Otherwise, it would be wise to have an appropriate
122
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
Table 1 Our solutions vs best published solutions Test cases
Best published
Our best results
NV/TD
Ref.
NV/TD
CT (s)
C101* C102* C103* C104* C105* C106* C108* C109*
10/828.94 10/828.94 10/828.06 10/824.78 10/828.94 10/828.94 10/828.94 10/828.94
RT RT RT RT RT RT RT RT
10/828.94 10/828.94 10/828.06 10/824.78 10/828.94 10/828.94 10/828.94 10/828.94
140 164 180 412 170 161 203 175
C201* C202* C203* C204* C205* C206* C207* C208*
3/591.56 3/591.56 3/591.17 3/590.60 3/588.88 3/588.49 3/588.29 3/588.32
RT RT RT RT RT RT RT RT
3/591.56 3/591.56 3/591.17 3/590.60 3/588.88 3/588.49 3/588.29 3/588.32
1356 1102 1334 1348 1084 1124 1159 1249
R101* R102 R103** R104** R105 R106 R107** R108 R109 R110 R111 R112
19/1650.80 17/1486.12 13/1292.85 9/1013.32 14/1377.11 12/1252.03 10/1113.69 9/964.38 11/1194.73 10/1124.40 10/1099.46 9/1003.73
RT RT HG HG RT RT CLM CLM HG RGP HG HG
19/1650.80 17/1486.41 13/1292.68 9/1007.31 14/1381.37 12/1269.72 10/1104.66 9/986.25 11/1208.96 10/1159.35 11/1066.32 10/967.88
835 1078 1926 2200 975 2296 908 1263 1356 1344 2016 1492
R201* R202 R203 R204** R205 R206 R207** R208 R209 R210 R211
4/1252.37 3/1191.70 3/942.70 2/849.62 3/994.42 3/912.97 2/914.39 2/730.771 3/909.86 3/939.373 2/910.09
HG RGP HG CLM RGP RT RT BFSKP RGP BFSKP HG
4/1252.37 4/1084.767 3/949.396 2/849.05 3/1032.55 3/931.62 2/905.13 2/732.8 3/930.59 3/1018.95 3/801.81
3325 3355 1510 4637 4470 3196 4856 6500 2805 3565 4481
RC101 RC102 RC103 RC104 RC105 RC106** RC107 RC108
14/1696.94 12/1554.75 11/1262.02 10/1135.48 13/1633.72 11/1427.13 11/1230.54 10/1139.82
TBGGP TBGGP RT CLM RGP CLM TBGGP TBGGP
15/1658.62 13/1513.60 11/1319.99 10/1141.09 13/1637.62 11/1424.73 11/1240.66 10/1147.42
1006 1631 505 683 748 1201 282 1278
RC201 RC202**
4/1046.94 3/1389.57
CLM HG
4/1425.21 3/1374.27
2039 3578
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
123
Table 1 (continued) Test cases RC203 RC204 RC205 RC206 RC207 RC208
Best published
Our best results
NV/TD
Ref.
NV/TD
CT (s)
3/1060.45 3/799.12 4/1302.42 3/1153.93 3/1062.05 3/829.69
HG HG HG RGP CLM RGP
3/1088.53 3/818.66 4/1304.64 3/1159.03 3/1107.16 3/862.34
1956 3185 1753 2475 3970 2398
value of K to stop the algorithm in reasonable computational time. It is not difficult to derive that a small K is not good for both diversification and intensification. The reason is that a small K allows few opportunities to explore within N ðSg Þ itself, not mentioning the solutions which are outside of N ðSg Þ but can be explored in TSA in each restart. Too small a value of L is not good for diversification either, since a small L results in that, only a small number of solutions which are outside of N ðSg Þ can be explored in TSA. Thus, the search is almost confined within N ðSg Þ. However, a large L results in heavy computational efforts. We propose L ¼ 4 for solving the VRPTW. d and T0 affect the cooling process. A commonly used cooling rule specifies that Tk ¼ dTk1 , k ¼ 1; 2; . . . ; here d is less than but close to 1. Since higher d and T0 values result in slower cooling processes. We propose d 0:95 and T0 ¼ 50 for solving the VRPTW.
5. Computational experiments 5.1. Benchmark problem instances We conducted computational experiments on SolomonÕs 56 benchmark problems [25]. These problems were generated in six classes: C1, C2, R1, R2, RC1 and RC2. All problems have 100 customers, a central depot, capacity constraints and time window constraints. The customers in C1 and C2 problems are clustered, and are randomly distributed in R1 and R2 problems. In RC1 and RC2 problems, customers are partially clustered and partially randomly distributed. The C1, R1 and
RC1 problems have short scheduling horizon, while C2, R2 and RC2 have longer scheduling horizon. 5.2. Experiments on Solomon’s problems The hardware used to run our experiments was the PC Pentium III 545 MHz. The algorithm was developed in the C++ programming language with the standard template library. The OS used was the linux kernel 2.4.0. The values of the parameters used were T0 ¼ 50, d ¼ 0:95, K ¼ 20dm=3e, L ¼ 4. Table 1 shows the best results obtained by our algorithms and the best published results obtained so far by heuristics proposed in [20] (RT), [3] (CR), [26] (TBGGP), [11] (HG), [5] (CLM), [1] (BFSKP) and [21] (RGP). In the table, NV means number of vehicles used, TD means travel distance, Ref. means the reference which obtained the corresponding best result; CT means computational time in seconds. As shown in the table by CT, our computational time is reasonable. On all of the 56 problem instances, our approach obtained seven new best solutions and matched 19 best-known solutions. The new best solutions obtained by our algorithms are for: R103, R104, R107, R204, R207, RC106 and RC202. Our other solutions are very close to the best-known solutions. In Table 1, a problem marked with (*) indicates that we matched the best known solution, and a problem marked with (**) indicates that our solution is better than the bestknown solution. Table 2 shows 7 new best solutions obtained by our approach. Here, TD and NV have the same
124
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
Table 2 New best solutions obtained by our approach Prob
NC
Details of routes
TD
NV
R103
5 12 7 8 6 6 10 7 6 7 8 8 10
0–40–53–12–68–80–0 0–92–98–14–44–38–86–16–61–85–91–100–37–0 0–36–64–49–63–90–32–70–0 0–71–65–78–34–35–81–77–28–0 0–50–33–76–79–10–31–0 0–94–96–95–97–87–13–0 0–7–19–11–8–46–47–48–82–18–89–0 0–52–62–88–84–17–93–59–0 0–60–45–83–5–99–6–0 0–2–22–75–56–4–25–54–0 0–26–39–23–67–55–24–29–3–0 0–27–69–30–9–66–20–51–1–0 0–42–43–15–57–41–74–72–73–21–58–0
1292.676
13
R104
12 11 11 13 10 10 11 10 12
92–98–14–44–38–86–16–61–85–91–100–37 42–43–15–87–57–41–22–74–73–21–26 52–88–62–11–64–63–90–32–10–31–27 6–94–96–59–93–99–84–17–45–83–5–60–89 18–82–48–46–8–47–36–49–19–7 70–30–20–66–9–35–65–71–51–50 28–1–69–76–53–40–2–13–95–97–58 72–75–56–23–67–39–55–4–25–54 68–80–24–29–79–78–34–81–33–3–77–12
1007.31
9
R107
10 11 11 8 11 8 12 9 10 10
0–26–21–39–23–67–55–24–29–3–77–0 0–27–69–1–76–79–78–34–68–80–12–28–0 0–52–7–62–11–64–49–19–48–82–18–89–0 0–40–53–6–96–59–95–13–58–0 0–60–83–45–5–99–87–57–2–74–72–73–0 0–20–65–9–66–71–35–81–50–0 0–94–92–98–44–14–38–86–16–91–100–37–97–0 0–47–36–46–8–84–17–61–85–93–0 0–42–43–15–41–22–75–56–4–25–54–0 0–33–51–30–88–31–10–63–90–32–70–0
1111.313
10
R204
50
0–27–52–18–7–82–48–19–11–62–88–31–69–76–12–67–39–53–28–79–81–9–51– 20–66–65–71–35–34–78–29–24–55–4–72–74–73–21–2–13–97–37–98–100–91–85– 93–59–96–6–89–0 0–94–95–92–42–57–22–75–56–23–41–15–43–14–44–38–86–16–61–99–87–60–83– 8–84–5–17–45–46–47–36–49–64–63–90–32–10–30–70–1–50–33–3–77–68–80–54– 25–26–40–58–0
849.05
2
0–1–50–33–81–65–34–29–24–39–67–23–15–43–14–44–38–86–16–84–5–99–6–18– 8–48–47–49–19–10–63–90–32–66–71–35–20–70–31–7–82–17–61–85–91–100–37– 98–93–96–0 0–42–92–59–60–83–45–46–36–64–11–62–88–52–27–69–30–51–9–78–79–3–76– 28–53–40–2–87–57–41–22–73–21–72–74–75–56–4–25–55–54–80–68–77–12–26– 58–13–97–95–94–89–0
905.13
2
1424.734
11
50
R207
49
51
RC106
9 8 9 10 9 11 8
0–69–98–88–53–12–10–9–13–17–0 0–14–11–87–59–75–97–58–74–0 0–95–62–63–85–76–51–84–56–66–0 0–42–44–39–40–36–38–41–43–37–35–0 0–72–71–67–30–32–34–50–93–80–0 0–2–45–5–8–7–6–46–4–3–1–100–0 0–92–61–81–90–94–96–54–68–0
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
125
Table 2 (continued) Prob
RC202
NC
Details of routes
10 10 7 9
0–82–52–99–86–57–22–49–20–24–91–0 0–65–83–64–19–23–21–18–48–25–77–0 0–33–31–29–27–28–26–89–0 0–15–16–47–78–73–79–60–55–70–0
32
0–65–82–12–14–47–15–11–69–64–19–23–48–18–76–51–22–86–87–9–57–52–10– 97–59–74–13–17–7–4–60–100–70–0 0–98–45–5–3–1–42–36–37–39–44–61–88–73–16–99–53–78–79–8–6–46–2–55–68– 43–35–54–96–93–94–80–0 0–92–95–85–63–33–28–26–27–29–31–30–62–67–71–72–38–40–41–81–90–91–84– 49–20–83–66–56–50–34–32–89–24–21–25–77–75–58–0
31 37
TD
NV
1374.267
3
Table 3 Comparisons of best averages on all data sets (VRPTW) Data sets
RT
CR
TBGGP
HG
CLM
BFSKP
Ours
C1 MNV C1 MTD
10.00 828.38
10.00 828.38
10.00 828.38
10.00 828.38
10.00 828.38
10.00 828.38
10.00 828.38
C2 MNV C2 MTD
3.00 589.86
3.00 591.42
3.00 589.86
3.00 589.86
3.00 589.86
3.00 589.86
3.00 589.86
R1 MNV R1 MTD
12.25 1208.50
12.17 1204.19
12.17 1209.35
11.92 1220.97
12.08 1210.14
12.41 1200.54
12.08 1215.14
R2 MNV R2 MTD
2.91 961.72
2.73 986.32
2.82 980.27
2.73 968.55
2.73 969.58
RC1 MNV RC1 MTD
11.88 1377.39
11.875 1397.44
11.50 1389.22
11.50 1388.24
11.50 1389.78
12 1383.21
11.75 1385.47
RC2 MNV RC2 MTD
3.38 1119.59
3.25 1229.54
3.38 1117.44
3.25 1140.43
3.25 1134.52
3.38 1116.51
3.25 1142.48
3 936.509
2.91 953.428
meanings as those in Table 1, NC means number of customers.
proach is competitive with those heuristics on solving the VRPTW.
5.3. Comparison with other heuristics
5.4. Sensitivity analysis of K by experimentation
Table 3 compares the mean number of vehicles (MNV) and the mean travel distance (MTD) obtained by our algorithms and other heuristics for the VRPTW. This table shows that, our results match the best results on data sets C1 and C2. We get best MNV on R1 with CLM. On R2 and RC1, our results are moderate. Finally, on RC2, we obtained best MNV with CR and HG. While the MTD on this data set is only slightly inferior to that of HG. The comparison shows that our ap-
To show the sensitivity of K which affects the solution quality and the computational time, we run the algorithm on C104 by setting K ¼ 1; 2; . . . Fig. 6 shows that a Kc exists. When K < Kc , the travel distances fluctuate very much, and the computational time increases quickly. When K P Kc , the travel distance tends to change between a constant, while the tendency of computational time changes not so abruptly as does for K < Kc , if we analyze them by regression models.
126
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
Fig. 6. Sensitivity analysis of K by experimentation.
The main reason may lie in that a K > Kc introduces (K Kc ) restarts for additional computational time, but has little influence on convergence behavior of solutions after Kc .
6. Conclusion In this paper, we proposed a TSA metaheuristic which runs with best-based K-restarts strategy, to diversify and intensify local searches for solving combinatorial optimization problems as a general framework. The diversification by TSA and the intensification by a K-restart strategy was discussed with details. The sensitivity of K was analyzed by experimentation. To show the effectiveness of the algorithm, we implemented it to solve the VRPTW. The computational results show that our algorithm is comparable to other approaches. The approach can be extended to solve other variants of VRP [14] and other combinatorial optimization problems.
References [1] B. Backer, V. Furnon, P. Shaw, P. Kilby, P. Prosser, Solving vehicle routing problem using constraint programming and metaheuristics, Journal of Heuristics 6 (2000) 501–523.
[2] W. Chiang, R. Russell, Simulated annealing metaheuristics for the vehicle routing problem with time windows, Annals of Operations Research 63 (1996) 3–27. [3] W. Chiang, R. Russell, A reactive tabu search metaheuristic for the vehicle routing problem with time windows, INFORMS Journal on Computing 9 (1997) 417–430. [4] J. Cordeau, G. Desaulniers, J. Desrosiers, M. Solomon, F. Soumis, The VRP with time windows, in: Paolo Toth, Daniele Vigo (Eds.), The Vehicle Routing Problem, SIAM Monographs on Discrete Mathematics and Applications, 2002 pp. 157–193 (Chapter 7). [5] J. Cordeau, G. Laporte, A. Mercier, A unified tabu search heuristic for vehicle routing problems with time windows, Technical report CRT-00–03, Center for Research on Transportation, Montreal, Canada, 2000. [6] J. Desrochers, J. Lenstra, M. Savelsbergh, F. Soumis, Vehicle routing with time windows: Optimization and approximation, in: Vehicle Routing: Methods and Studies, Elsevier Science Publishers, North-Holland, 1988. [7] M. Desrochers, J. Desrosiers, M. Solomon, A new optimization algorithm for the vehicle routing problem with time windows, Operations Research 40 (1992) 342–354. [8] J. Desrosiers, Y. Dumas, M. Solomon, F. Sormis, Time constrained routing and scheduling, in: Handbooks in Operations Research and Management Science: Network Routing, Elsevier Science Publishers, Amsterdam, 1995. [9] Y. Dumas, J. Desrosiers, E. Gelinas, M. Solomon, An optimal algorithm for the traveling salesman problem with time windows, Operations Research 43 (1995) 367–371. [10] M. Fisher, K. J€ ornsten, O. Madsen, Vehicle routing with time windows: Two optimization algorithms, Operations Research 45 (1997) 488–492. [11] J. Homberger, H. Gehring, Two evolutionary metaheuristics for the vehicle routing problem with time windows, INFOR 37 (1999) 297–318. [12] N. Kohl, O. Madsen, An optimization algorithm for the vehicle routing problem with time windows based on lagrangian relaxation, Operations Research 45 (1997) 395– 406. [13] A. Kolen, A. Rinnooy, H. Trienekens, Vehicle routing with time windows, Operations Research 35 (1987) 266–273. [14] H. Li, A. Lim, A metaheuristic for the pickup and delivery problem with time windows, in: ICTAI 2001, November 2001, pp. 151–158. [15] J. Potvin, S. Bengio, The vehicle routing problem with time windows-part 2: Genetic search, INFORMS Journal on Computing 8 (1996) 165–172. [16] J. Potvin, B. Garcia, L. Rousseau, The vehicle routing problem with time windows––Part 1: Tabu search, INFORMS Journal on Computing 8 (1996) 158–164. [17] J. Potvin, J. Rousseau, A parallel route building algorithms for the vehicle routing and scheduling problem with time windows, European Journal of Operational Research 66 (1993) 331–340. [18] J. Potvin, J. Rousseau, An exchange heuristic for routing problems with time windows, Journal of Operational Research Society 46 (1995) 1433–1446.
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127 [19] Y. Rochat, F. Semet, A tabu search approach for delivering pet food and flour in switzerland, Journal of Operational Research Society 45 (1994) 1233–1246. [20] Y. Rochat, E. Taillard, Probabilistic diversification and intensification in local search for vehicle routing, Journal of Heuristics 1 (1995) 147–167. [21] L. Rousseau, M. Gendreau, G. Pesant, Using constraintbased operators to solve the vehicle routing problem with time windows, Journal of Heuristics 8 (2002) 43–58. [22] R. Russell, Hybrid heuristics for the vehicle routing problem with time windows, Transportation Science 29 (1995) 156–166. [23] M. Savelsbergh, Local search for routing problems with time windows, Annals of Operations Research 4 (1985) 285–305. [24] M. Savelsbergh, The vehicle routing problem with time windows: Minimizing route duration, ORSA Journal on Computing 4 (1992) 146–154.
127
[25] M. Solomon, Algorithms for the vehicle routing and scheduling problems with time window constrains, Operations Research 35 (1987) 254–264. [26] E. Taillard, P. Badeau, M. Gendreau, F. Geurtin, J. Potvin, A tabu search heuristic for the vehicle routing problem with time windows, Transportation Science 31 (1997) 170–186. [27] S. Thangiah, Vehicle routing with time windows using genetic algorithms, Technical report SRU-CPSC-TR-93– 23, Computer Science Department, Slippery Rock University, Slippery Rock, PA, 1993. [28] S. Thangiah, I. Osman, T. Sun, Hybrid genetic algorithm, simulated annealing and tabu search methods for vehicle routing problems with time windows, Technical report UKC/OR94/4, Institute of Mathematics and Statistics, University of Kent, Canterbury, UK, 1994.
Discrete Optimization
Local search with annealing-like restarts to solve the VRPTW Haibing Li, Andrew Lim
*
Department of Computer Science, National University of Singapore, 3 Science Drive 2, Singapore 117543, Singapore Received 25 April 2001; accepted 3 May 2002
Abstract In this paper, we propose a metaheuristic based on annealing-like restarts to diversify and intensify local searches for solving the vehicle routing problem with time windows (VRPTW). Using the SolomonÕs benchmark instances for the problem, our method obtained seven new best results and equaled 19 other best results. Extensive comparisons indicate that our method is comparable to the best in published literature. This approach is flexible and can be extended to handle other variants of vehicle routing problems and other combinatorial optimization problems. Ó 2002 Elsevier Science B.V. All rights reserved. Keywords: Routing; Vehicle routing problem with time windows; Local search; Diversification; Intensification
1. Introduction Local searches are commonly used to solve hard combinatorial optimization problems. A local search algorithm starts with an initial solution and searches within neighborhoods for better solutions. Being trapped in local optima is a significant factor affecting the quality of solutions obtained by local searches. Advanced search algorithms, including simulated annealing (SA), tabu search (TS) and evolutionary algorithms (EA), were devised to help escaping from local optimal traps. These approaches are classified as metaheuristics. From the viewpoint of local search, one of the important issues in designing
* Corresponding author. Tel.: +65-874-6891; fax: +65-7794580. E-mail addresses: [email protected] (H. Li), alim@ comp.nus.edu.sg (A. Lim).
these algorithms lies in diversifying the search spaces using respective strategies, so that the whole solution space can be explored as widely as possible. Another important issue in designing local search algorithms is intensification. A local search algorithm is more efficient if it is smart enough to intensify the search on occasion, for checking if there are better solutions. With regards to the diversification and intensification issues, in this paper, we propose a tabu-embedded simulated annealing (TSA) metaheuristic which runs with a best-based K-restarts strategy during local searches within three constraint-based neighborhoods, for solving the vehicle routing problem with time windows (VRPTW). The sensitivity of K, which affects the solution quality and computational effort, is analyzed by experimentation. The rest of the paper is organized as follows. In Section 2, the notation is presented and then the VRPTW is briefly described with a review of
0377-2217/03/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. doi:10.1016/S0377-2217(02)00486-1
116
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
related work. Section 3 describes the important issues for designing our local search algorithms. In Section 4, we present the algorithm in detail. Computational experiments are conducted in Section 5, in addition, results on 56 SolomonÕs benchmark problem instances are reported and our approach is compared with the best methods for this problem. Finally, we present our conclusions in Section 6.
2. Notations and the VRPTW 2.1. Notations For the convenience of description, we define the following notation: GðV ; AÞ G is a digraph generated from problem data, V is the node set of G and A is the arc set of G n, m number of customers, number of vehicle used, respectively L the maximum number of consecutively non-improving iterations in the TSA procedure. It is proposed to be a small number K the maximum number of restarts of the TSA procedure on the current best solution CostðSÞ objective function of solution S of local search SACostðSÞ objective function of solution S in tabuembedded Metropolis procedure in the SA algorithm DistðSÞ, STðSÞ, WTðSÞ travel distance, schedule time and waiting time of solution S, respectively N) ðSÞ, N () ðSÞ, N,! ðSÞ constraint-based neighborhoods of solution S obtained by shift operator, exchange operator and rearrange operator respectively N ðSÞ N ðSÞ ¼ N) ðSÞ [ N () ðSÞ are the neighborhoods where local searches are conducted / the depth of K-ary tree due to K restarts T0 , T initial and global annealing temperature respectively
d cooling ratio of T a, b, c, k penalty factor for m, DistðSÞ, STðSÞ, WTðSÞ respectively the customer node at the position of pos vpos in a route Dpos1;pos2 the distance between customers at position pos1 and pos2 in a route 2.2. The VRPTW The VRPTW is a generalization of classical vehicle routing problems with the additional complexity of time window constraints. Let G ¼ ðV ; AÞ be a digraph. V ¼ VN [ fv0 g is the node set where VN ¼ fvi 2 V ji ¼ 1; 2; . . . ; ng represents the customers and node v0 denotes the depot where a fleet of vehicles are housed. Each node vi 2 V has an associated customer demand qi (q0 ¼ 0), a service time si (s0 ¼ 0) and a service-time window ½ei ; li . For each pair of nodes hvi ; vj i (i 6¼ j, i; j ¼ 0; 1; 2; . . . ; n), a non-negative distance dij and a non-negative travel time tij are known. Due to the time window constraints, arcs may not exist between some node pairs. Therefore, the arc set can be defined as A ¼ fhvi ; vj ijvi ; vj 2 V ; vi 6¼ vj ; t0i þ si þ tij 6 lj g. If a vehicle reaches a customer vi before ei , it needs to wait until ei in order to service the customer. The schedule time of a route is the sum of the waiting time, the service and travel time. The objective of the VRPTW is to service all customers without violating vehicle capacity constraints and time window constraints with a minimum number of vehicles and, for the same number of routes with the minimum travel distance, followed by the minimum schedule time and the minimum waiting time. Savelsbergh [23] proved that the VRPTW is NP-hard. Much work has been done to derive exact algorithms to solve the problem. Desrochers et al. [7] proposed a column generation approach that solved seven of SolomonÕs benchmark instances exactly proposed in [25]. Fisher et al. [10] proposed a K-tree relaxation approach that solved two of SolomonÕs benchmark instances. Other exact approaches can be found in [9,12,13]. Because of the exponential time complexity of these approaches, it is unlikely that these algorithms can produce optimal solutions for practical-sized
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
VRPTW with reasonable computational time. Hence, much work has focused on generating near optimal solutions using heuristic approaches. In general, heuristics for VRPTW can be classified as classical heuristics or metaheuristics. Classical heuristics include construction heuristics, improvement heuristics and composite heuristics. Route construction heuristics build feasible solutions by inserting an unrouted customer into a current partial route until all customers are routed. The sequential insertion algorithm proposed by Solomon [25] belongs to construction heuristics. Its parallel version was implemented in [17]. These methods often produce solutions very quickly, but the quality of these solutions may not be good. For this reason, construction heuristics are usually applied to obtain initial solutions for improvement heuristics or other two-phase heuristics. Route improvement heuristics modify a solution by performing local searches for better neighboring solutions within its neighborhoods generated by node/edge-swapping operators. Edge-exchange heuristics are referred to [18,23,24]. Efficient implementations of edge-exchange for speeding up the pruning of inferior solutions are reported in [23,24,26]. Russell [22] proposed a composite heuristic which mixes route construction and route improvement procedures. Metaheuristics are a new generation of heuristics developed in the recent decade. Near optimal solutions for the VRPTW were produced by TS [3,5,16,19,20,26], SA [2], EA [11,15,27] and hybrid algorithms of the above three algorithms [22,28]. Surveys on the VRPTW are referred to [4,6,8].
3. Issues for designing local search algorithms During designing the local search algorithms, we take several important issues into consideration. Firstly, the neighborhood generating mechanisms, which usually have a strong connection with the problems and the solution representation, provide spaces within which local searches are conducted. For example, for the vehicle routing problems, the solutions are naturally represented by the routes with sequences of customer nodes
117
which are visited by the vehicles to offer a service. So it is a natural idea to obtain neighboring solutions by swapping nodes/edges between routes. Secondly, it is also important to evaluate the quality of solutions during local searches. Different problems have different cost functions which faithfully reflect the problemsÕ objectives and act as solution evaluators. For problems with multiple objectives, penalty functions are commonly used for this purpose. Thirdly, the way to obtain and accept a neighboring solution, affects the exploring trail of local searches. This reflects an important issue in local searches, i.e. diversification. For example, a descent local search (DLS) systematically exploring the whole neighborhood of a seed solution, and discards infeasible/inferior solutions immediately after solution evaluation while it keeps better solutions as new seed solutions for repeating the search procedure. The procedure stops at a dead end when no more better solutions can be found within a seed solutionÕs neighborhood. This shortsighted strategy inevitably causes local searches to be trapped in local optima very quickly. To help jumping out of the dead ends of DLS, in addition to systematically explore the neighborhoods, TS intelligently goes a step further to label along its exploring trail and to accept unexploited local optima no matter whether they are better or worse than seed solutions. In contrast, SA obtains a neighboring solution randomly but accepts it with a probability reflecting annealing property. Therefore, both TS and SA diversify local searches to escape from local optima. However, TS does so in a deterministic way so that its exploring trail is fixed, while SA does so in a non-deterministic way and its exploring trails vary in different tries. In this work, we propose a TSA procedure to diversify local searches. Finally, another important local search issue is intensification. Both TS and SA do not deal specially with this issue. In past researches, postprocessing on the best solution found so far was used frequently for this purpose. However, it is a pity that this intensification was only used after local searches were finished, instead of being embedded in the local search procedure. In this work,
118
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
we propose a K-restart strategy which guides the TSA to search around the current best solution. If a new better solution is found within K restarts of TSA, the current best solution is updated and a new series of restarts are initiated. Otherwise, if K consecutive restarts are finished without finding a better solution, this procedure stops. Therefore, the intensification which is in favor of the current best solution always plays its role during local searches. Implementation details for these issues will be presented in Section 4.
4. The algorithm The algorithm is proposed as a two-phase framework: Phase 1: Obtain an initial solution for the VRPTW using a heuristic. Phase 2: Start local search from the initial solution, using TSA with K-restart strategy. In phase 1, we use the insertion heuristic proposed by Solomon [25] to produce initial solutions. The initial solutions are then used to feed the local search procedure in phase 2. Of course, the local search in phase 2 are of our main interest. We deal with them in details with regards to the issues discussed in Section 3.
change operator are used for generating neighborhoods for local searches in the TSA procedure, while the rearrange operator is only used for postprocessing on the local optima obtained by local searches. 4.1.1. The shift operator Using the shift operator, non-null route segments (with at least one node) can be moved from one route to another route under the condition of not violating the capacity constraints and time window constraints. For each pair of selected routes, say, Route i and Route j (we assume 1 6 i < j 6 m, i 6¼ j without loss of generality), the shift operator is used in two directions, so that a sub-neighborhood Ni)j ðSÞ can be obtained by shifting route segments from Route i to Route j and another sub-neighborhood Nj)i ðSÞ can be obtained by shifting route segments from Route j to Route i. We have denoted the neighborhood of solution S obtained by the shift operam1 tor as N) ðSÞ, thus, N) ðSÞ ¼ [i¼1 [mj¼iþ1 fNi)j ðSÞ[ Nj)i ðSÞg. Fig. 1 shows an example of shifting a route segment from Route 1 to Route 2. To accelerate the local searches, three simple conditions are considered to prune infeasible or inferior solutions in constant time, during local searches within N) ðSÞ. We use the example shown in Fig. 1 to explain this. Firstly, if the total-capacity of the route segment from pos11 to pos12 to be shifted in
4.1. Neighborhood generating mechanisms In our method, three edge-swapping operators are used to generate local search neighborhoods. These operators are, the shift operator, the exchange operator and the rearrange operator. Savelsbergh [23,24] and Taillard et al. [26] also dealt with similar edge-exchange operators for the VRPTW. But they restricted the length of edge segments to be swapped to certain maximum values. We do not impose such restrictions. Further, we use three simple conditions that can be calculated in constant time to prune infeasible/inferior solutions, so that the neighborhoodsÕ sizes are reduced to accelerate the local searches. Among these operators, the shift operator and the ex-
Fig. 1. The shift operator.
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
Route 1 exceeds the left-capacity of Route 2, the solution obtained by shifting is infeasible so that it can be discarded. We record the used-capacityso-far in each route node so that computing total-capacity and left-capacity can be done in constant time. Secondly, if either hvpos prev ; vpos11 i or hvpos12 ; vpos i is an infeasible arc in digraph G defined in Section 2, the shift operation also results in an infeasible solution that can be discarded. Since the arc set A is initialized when the graph data is generated, checking for this condition can be done in constant time. Finally, since we are only interested in searching for better solutions, if an edge-shift incurs an increment in the total travel distance of the two routes, then this shift operation results in an inferior solution that can be discarded too. To simplify the checking procedure such that if an increment in travel distance incurs, we can use the following condition, let: D ¼ ðDpos11 prev;pos12 next þ Dpos prev;pos11 þ Dpos12;pos Þ ðDpos11 prev;pos11 þ Dpos12;pos12 next þ Dpos prev;pos Þ ð1Þ Here, ðDpos11 prev;pos12 next þ Dpos prev;pos11 þ Dpos12;pos Þ is the total distance of the three newly connected edges, while ðDpos11 prev;pos11 þ Dpos12;pos12 next þ Dpos prev;pos Þ is the total distance of the three old edges removed due to shifting. D P 0 means that the travel distance will not be decreased, hence the solution obtained by the operation does not result in a better solution. However, as our first objective is to minimize the number of vehicles before excluding a shift operation under the condition of D P 0, we must make sure that this shift operation will move the whole route segment out of a route so that we can eliminate one vehicle. Therefore, we will only discard a solution if D P 0 and the shifting operation will not eliminate a vehicle. As determined by Eq. (1), this check can be performed in constant time. Using these three conditions, jN) ðSÞj can be reduced efficiently. We name the solution space N) ðSÞ constrained by the three conditions as a constraint-based neighborhood.
119
4.1.2. The exchange operator Using the exchange operator, non-null route segments can be exchanged between two routes under the condition of not violating capacity and time window constraints. We have denoted the neighborhood of solution S obtained by this operator as N () ðSÞ. Fig. 2 shows an example of exchanging two route segments between Route 1 to Route 2. Like the shift operator, infeasible or inferior solutions obtained by exchange operation can also be pruned during local searches using the same three conditions. Of course, since more arcs and route segments are involved in the exchange operator, the conditions are a little more complex. However, all three conditions can also be calculated in constant time. Actually, it is possible to obtain an exchange operation by two shift operations. However, as we discussed previously, each shift operation requires checking of three conditions to prune infeasible/ inferior solutions. Thus, two shift operations need to check six conditions. This results in additional computational time and implementation complexity in contrast to the exchange operator. 4.1.3. The rearrange operator Using the rearrange operator, route segments can be repositioned within the same route under the condition of not violating the time window constraints. We have denoted the neighborhood of solution S obtained by the rearrange operator as
Fig. 2. The exchange operator.
120
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
Fig. 3. The rearrange operator.
N,! ðSÞ. Fig. 3 shows an example of repositioning a route segment to another position within a route. Similarly, we can also get a constraint-based neighborhood N,! ðSÞ using the second condition and the third condition discussed above. The first condition is not needed because the capacity remains unchanged when the rearranging operation is applied within one route.
with four fields. Since the probability that two different solutions having same record values is very small, it is reasonable to regard two solutions as the same if their record values are not different. We store the records of the solutions in a tabu list which differs from some references that use edgemoves in their tabu structure. Firstly, the pseudocode of the TSA procedure is described as follows:
4.2. Objectives and cost function in local search The priority-order of the objectives for solving the VRPTW is: (1) to minimize the number of vehicles used; (2) to minimize the travel cost; (3) to minimize the schedule time; (4) to minimize the driversÕ total waiting time to begin service. To reflect this order of priority, we adopt the following cost function: CostðSÞ ¼ am þ bDistðSÞ þ cSTðSÞ þ kWTðSÞ ð2Þ The penalty weight factors are set to be: a b c k. Using this setting, even if the initial solutions obtained by the route construction heuristic contain large number of vehicles, the vehicle size will be reduced by local searches. In addition, we can change the setting of these factors to get different priority orders of the four objectives. 4.3. Tabu-embedded simulated annealing with Krestart strategy To introduce the TS into local searches, we record the four objective values in a data structure
ALGORITHM. TSA(Solution S) 1. Sb S; 2. NonImproving 0 3. WHILE NonImproving < L DO 3.1 Use TMPðSÞ to obtain a non-tabu neighbor S 0 2 N ðSÞ 3.2 Sb0 Perform TDLS within N ðS 0 Þ / local search 3.3 Sb0 Perform TDLS within N,! ðSb0 Þ / postprocessing 3.4 IF CostðSb0 Þ < CostðSb Þ THEN Sb Sb0 ; NonImproving 0 3.5 ELSE NonImproving NonImproving þ 1 3.6 S Sb0 4. RETURN Sb In the above procedure, TDLS is a tabuembedded DLS procedure, in which only nontabu neighbors are accepted. The procedure TMP is the tabu-embedded metropolis procedure reflecting the annealing property. It is described as follows:
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
ALGORITHM. TMP(Solution S) 1. WHILE NOT yet accept a neighbor S 0 2 N ðSÞ DO 1.1 Get a random non-tabu neighbor S 0 2 N ðSÞ 1.2 D SACostðS 0 Þ SACostðSÞ 1.3 IF D 6 0 THEN prob 1 1.4 ELSE prob eD=T 1.5 IF randomð0; 1Þ 6 prob THEN 1.5.1 Accept S 0 1.5.2 Update annealing temperature: T dT 1.5.3 Record S 0 into tabu list 2. RETURN S 0
121
Fig. 4. Diversification by TSA.
In the procedure, SACostðsÞ is a cost function to evaluate solutions on behalf of SA. It is defined to adjust the cooling procedure with setting of T0 , T and d. It can be defined as Eq. (3). SACostðSÞ ¼ DistðSÞ þ uWTðSÞ; where u ¼ 0:01 DistðSÞ
ð3Þ
Finally, the algorithm which combines TSA with the best-based K-restart strategy is as follows: INPUT: an initial solution S0 OUTPUT: a best solution Sg found so far ALGORITHM. K-TSA(Solution S0 ) 1. Configure parameters: K, L, d and T0 ; T T0 2. Set tabu list to be empty 3. Sg S0 ; gNonImproving 0 4. WHILE gNonImproving < K DO / K-restart 4.1 S Sg ; 4.2 Sg0 TSAðSÞ 4.3 IF CostðSg0 Þ < CostðSg Þ THEN / Update new current best solution Sg Sg0 ; gNonImproving 0 4.4 ELSE gNonImproving gNonImproving þ 1 5. OUTPUT Sb As shown in Fig. 4, TSA is in fact a depth-like search which starts from a neighbor in N ðSg Þ, and stretches the search away from N ðSg Þ by TMP and TDLS. Actually, only with this depth-like search way, can solutions outside of N ðSg Þ be explored in
Fig. 5. Intensification by K-restart strategy.
each restart of TSA on Sg . Thus, the local searches around Sg is diversified if L is not too small. On the other hand, the K-restart strategy guides TSA to start local searches within N ðSg Þ in a breadth-like way. It intensifies local searches in favor of N ðSg Þ. The search behavior of K-TSA is shown in Fig. 5. 4.4. Configuration of K, L, T0 and d First of all, K and L are the key factors that affect the diversification and intensification behavior. It would be preferable to have large value for both K and L, but the computational cost would be very expensive. The largest possible value of K is jN ðSg Þj. If jN ðSg Þj is not too large, it would be preferable to accept each neighbor of Sg systematically. Otherwise, it would be wise to have an appropriate
122
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
Table 1 Our solutions vs best published solutions Test cases
Best published
Our best results
NV/TD
Ref.
NV/TD
CT (s)
C101* C102* C103* C104* C105* C106* C108* C109*
10/828.94 10/828.94 10/828.06 10/824.78 10/828.94 10/828.94 10/828.94 10/828.94
RT RT RT RT RT RT RT RT
10/828.94 10/828.94 10/828.06 10/824.78 10/828.94 10/828.94 10/828.94 10/828.94
140 164 180 412 170 161 203 175
C201* C202* C203* C204* C205* C206* C207* C208*
3/591.56 3/591.56 3/591.17 3/590.60 3/588.88 3/588.49 3/588.29 3/588.32
RT RT RT RT RT RT RT RT
3/591.56 3/591.56 3/591.17 3/590.60 3/588.88 3/588.49 3/588.29 3/588.32
1356 1102 1334 1348 1084 1124 1159 1249
R101* R102 R103** R104** R105 R106 R107** R108 R109 R110 R111 R112
19/1650.80 17/1486.12 13/1292.85 9/1013.32 14/1377.11 12/1252.03 10/1113.69 9/964.38 11/1194.73 10/1124.40 10/1099.46 9/1003.73
RT RT HG HG RT RT CLM CLM HG RGP HG HG
19/1650.80 17/1486.41 13/1292.68 9/1007.31 14/1381.37 12/1269.72 10/1104.66 9/986.25 11/1208.96 10/1159.35 11/1066.32 10/967.88
835 1078 1926 2200 975 2296 908 1263 1356 1344 2016 1492
R201* R202 R203 R204** R205 R206 R207** R208 R209 R210 R211
4/1252.37 3/1191.70 3/942.70 2/849.62 3/994.42 3/912.97 2/914.39 2/730.771 3/909.86 3/939.373 2/910.09
HG RGP HG CLM RGP RT RT BFSKP RGP BFSKP HG
4/1252.37 4/1084.767 3/949.396 2/849.05 3/1032.55 3/931.62 2/905.13 2/732.8 3/930.59 3/1018.95 3/801.81
3325 3355 1510 4637 4470 3196 4856 6500 2805 3565 4481
RC101 RC102 RC103 RC104 RC105 RC106** RC107 RC108
14/1696.94 12/1554.75 11/1262.02 10/1135.48 13/1633.72 11/1427.13 11/1230.54 10/1139.82
TBGGP TBGGP RT CLM RGP CLM TBGGP TBGGP
15/1658.62 13/1513.60 11/1319.99 10/1141.09 13/1637.62 11/1424.73 11/1240.66 10/1147.42
1006 1631 505 683 748 1201 282 1278
RC201 RC202**
4/1046.94 3/1389.57
CLM HG
4/1425.21 3/1374.27
2039 3578
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
123
Table 1 (continued) Test cases RC203 RC204 RC205 RC206 RC207 RC208
Best published
Our best results
NV/TD
Ref.
NV/TD
CT (s)
3/1060.45 3/799.12 4/1302.42 3/1153.93 3/1062.05 3/829.69
HG HG HG RGP CLM RGP
3/1088.53 3/818.66 4/1304.64 3/1159.03 3/1107.16 3/862.34
1956 3185 1753 2475 3970 2398
value of K to stop the algorithm in reasonable computational time. It is not difficult to derive that a small K is not good for both diversification and intensification. The reason is that a small K allows few opportunities to explore within N ðSg Þ itself, not mentioning the solutions which are outside of N ðSg Þ but can be explored in TSA in each restart. Too small a value of L is not good for diversification either, since a small L results in that, only a small number of solutions which are outside of N ðSg Þ can be explored in TSA. Thus, the search is almost confined within N ðSg Þ. However, a large L results in heavy computational efforts. We propose L ¼ 4 for solving the VRPTW. d and T0 affect the cooling process. A commonly used cooling rule specifies that Tk ¼ dTk1 , k ¼ 1; 2; . . . ; here d is less than but close to 1. Since higher d and T0 values result in slower cooling processes. We propose d 0:95 and T0 ¼ 50 for solving the VRPTW.
5. Computational experiments 5.1. Benchmark problem instances We conducted computational experiments on SolomonÕs 56 benchmark problems [25]. These problems were generated in six classes: C1, C2, R1, R2, RC1 and RC2. All problems have 100 customers, a central depot, capacity constraints and time window constraints. The customers in C1 and C2 problems are clustered, and are randomly distributed in R1 and R2 problems. In RC1 and RC2 problems, customers are partially clustered and partially randomly distributed. The C1, R1 and
RC1 problems have short scheduling horizon, while C2, R2 and RC2 have longer scheduling horizon. 5.2. Experiments on Solomon’s problems The hardware used to run our experiments was the PC Pentium III 545 MHz. The algorithm was developed in the C++ programming language with the standard template library. The OS used was the linux kernel 2.4.0. The values of the parameters used were T0 ¼ 50, d ¼ 0:95, K ¼ 20dm=3e, L ¼ 4. Table 1 shows the best results obtained by our algorithms and the best published results obtained so far by heuristics proposed in [20] (RT), [3] (CR), [26] (TBGGP), [11] (HG), [5] (CLM), [1] (BFSKP) and [21] (RGP). In the table, NV means number of vehicles used, TD means travel distance, Ref. means the reference which obtained the corresponding best result; CT means computational time in seconds. As shown in the table by CT, our computational time is reasonable. On all of the 56 problem instances, our approach obtained seven new best solutions and matched 19 best-known solutions. The new best solutions obtained by our algorithms are for: R103, R104, R107, R204, R207, RC106 and RC202. Our other solutions are very close to the best-known solutions. In Table 1, a problem marked with (*) indicates that we matched the best known solution, and a problem marked with (**) indicates that our solution is better than the bestknown solution. Table 2 shows 7 new best solutions obtained by our approach. Here, TD and NV have the same
124
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
Table 2 New best solutions obtained by our approach Prob
NC
Details of routes
TD
NV
R103
5 12 7 8 6 6 10 7 6 7 8 8 10
0–40–53–12–68–80–0 0–92–98–14–44–38–86–16–61–85–91–100–37–0 0–36–64–49–63–90–32–70–0 0–71–65–78–34–35–81–77–28–0 0–50–33–76–79–10–31–0 0–94–96–95–97–87–13–0 0–7–19–11–8–46–47–48–82–18–89–0 0–52–62–88–84–17–93–59–0 0–60–45–83–5–99–6–0 0–2–22–75–56–4–25–54–0 0–26–39–23–67–55–24–29–3–0 0–27–69–30–9–66–20–51–1–0 0–42–43–15–57–41–74–72–73–21–58–0
1292.676
13
R104
12 11 11 13 10 10 11 10 12
92–98–14–44–38–86–16–61–85–91–100–37 42–43–15–87–57–41–22–74–73–21–26 52–88–62–11–64–63–90–32–10–31–27 6–94–96–59–93–99–84–17–45–83–5–60–89 18–82–48–46–8–47–36–49–19–7 70–30–20–66–9–35–65–71–51–50 28–1–69–76–53–40–2–13–95–97–58 72–75–56–23–67–39–55–4–25–54 68–80–24–29–79–78–34–81–33–3–77–12
1007.31
9
R107
10 11 11 8 11 8 12 9 10 10
0–26–21–39–23–67–55–24–29–3–77–0 0–27–69–1–76–79–78–34–68–80–12–28–0 0–52–7–62–11–64–49–19–48–82–18–89–0 0–40–53–6–96–59–95–13–58–0 0–60–83–45–5–99–87–57–2–74–72–73–0 0–20–65–9–66–71–35–81–50–0 0–94–92–98–44–14–38–86–16–91–100–37–97–0 0–47–36–46–8–84–17–61–85–93–0 0–42–43–15–41–22–75–56–4–25–54–0 0–33–51–30–88–31–10–63–90–32–70–0
1111.313
10
R204
50
0–27–52–18–7–82–48–19–11–62–88–31–69–76–12–67–39–53–28–79–81–9–51– 20–66–65–71–35–34–78–29–24–55–4–72–74–73–21–2–13–97–37–98–100–91–85– 93–59–96–6–89–0 0–94–95–92–42–57–22–75–56–23–41–15–43–14–44–38–86–16–61–99–87–60–83– 8–84–5–17–45–46–47–36–49–64–63–90–32–10–30–70–1–50–33–3–77–68–80–54– 25–26–40–58–0
849.05
2
0–1–50–33–81–65–34–29–24–39–67–23–15–43–14–44–38–86–16–84–5–99–6–18– 8–48–47–49–19–10–63–90–32–66–71–35–20–70–31–7–82–17–61–85–91–100–37– 98–93–96–0 0–42–92–59–60–83–45–46–36–64–11–62–88–52–27–69–30–51–9–78–79–3–76– 28–53–40–2–87–57–41–22–73–21–72–74–75–56–4–25–55–54–80–68–77–12–26– 58–13–97–95–94–89–0
905.13
2
1424.734
11
50
R207
49
51
RC106
9 8 9 10 9 11 8
0–69–98–88–53–12–10–9–13–17–0 0–14–11–87–59–75–97–58–74–0 0–95–62–63–85–76–51–84–56–66–0 0–42–44–39–40–36–38–41–43–37–35–0 0–72–71–67–30–32–34–50–93–80–0 0–2–45–5–8–7–6–46–4–3–1–100–0 0–92–61–81–90–94–96–54–68–0
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
125
Table 2 (continued) Prob
RC202
NC
Details of routes
10 10 7 9
0–82–52–99–86–57–22–49–20–24–91–0 0–65–83–64–19–23–21–18–48–25–77–0 0–33–31–29–27–28–26–89–0 0–15–16–47–78–73–79–60–55–70–0
32
0–65–82–12–14–47–15–11–69–64–19–23–48–18–76–51–22–86–87–9–57–52–10– 97–59–74–13–17–7–4–60–100–70–0 0–98–45–5–3–1–42–36–37–39–44–61–88–73–16–99–53–78–79–8–6–46–2–55–68– 43–35–54–96–93–94–80–0 0–92–95–85–63–33–28–26–27–29–31–30–62–67–71–72–38–40–41–81–90–91–84– 49–20–83–66–56–50–34–32–89–24–21–25–77–75–58–0
31 37
TD
NV
1374.267
3
Table 3 Comparisons of best averages on all data sets (VRPTW) Data sets
RT
CR
TBGGP
HG
CLM
BFSKP
Ours
C1 MNV C1 MTD
10.00 828.38
10.00 828.38
10.00 828.38
10.00 828.38
10.00 828.38
10.00 828.38
10.00 828.38
C2 MNV C2 MTD
3.00 589.86
3.00 591.42
3.00 589.86
3.00 589.86
3.00 589.86
3.00 589.86
3.00 589.86
R1 MNV R1 MTD
12.25 1208.50
12.17 1204.19
12.17 1209.35
11.92 1220.97
12.08 1210.14
12.41 1200.54
12.08 1215.14
R2 MNV R2 MTD
2.91 961.72
2.73 986.32
2.82 980.27
2.73 968.55
2.73 969.58
RC1 MNV RC1 MTD
11.88 1377.39
11.875 1397.44
11.50 1389.22
11.50 1388.24
11.50 1389.78
12 1383.21
11.75 1385.47
RC2 MNV RC2 MTD
3.38 1119.59
3.25 1229.54
3.38 1117.44
3.25 1140.43
3.25 1134.52
3.38 1116.51
3.25 1142.48
3 936.509
2.91 953.428
meanings as those in Table 1, NC means number of customers.
proach is competitive with those heuristics on solving the VRPTW.
5.3. Comparison with other heuristics
5.4. Sensitivity analysis of K by experimentation
Table 3 compares the mean number of vehicles (MNV) and the mean travel distance (MTD) obtained by our algorithms and other heuristics for the VRPTW. This table shows that, our results match the best results on data sets C1 and C2. We get best MNV on R1 with CLM. On R2 and RC1, our results are moderate. Finally, on RC2, we obtained best MNV with CR and HG. While the MTD on this data set is only slightly inferior to that of HG. The comparison shows that our ap-
To show the sensitivity of K which affects the solution quality and the computational time, we run the algorithm on C104 by setting K ¼ 1; 2; . . . Fig. 6 shows that a Kc exists. When K < Kc , the travel distances fluctuate very much, and the computational time increases quickly. When K P Kc , the travel distance tends to change between a constant, while the tendency of computational time changes not so abruptly as does for K < Kc , if we analyze them by regression models.
126
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127
Fig. 6. Sensitivity analysis of K by experimentation.
The main reason may lie in that a K > Kc introduces (K Kc ) restarts for additional computational time, but has little influence on convergence behavior of solutions after Kc .
6. Conclusion In this paper, we proposed a TSA metaheuristic which runs with best-based K-restarts strategy, to diversify and intensify local searches for solving combinatorial optimization problems as a general framework. The diversification by TSA and the intensification by a K-restart strategy was discussed with details. The sensitivity of K was analyzed by experimentation. To show the effectiveness of the algorithm, we implemented it to solve the VRPTW. The computational results show that our algorithm is comparable to other approaches. The approach can be extended to solve other variants of VRP [14] and other combinatorial optimization problems.
References [1] B. Backer, V. Furnon, P. Shaw, P. Kilby, P. Prosser, Solving vehicle routing problem using constraint programming and metaheuristics, Journal of Heuristics 6 (2000) 501–523.
[2] W. Chiang, R. Russell, Simulated annealing metaheuristics for the vehicle routing problem with time windows, Annals of Operations Research 63 (1996) 3–27. [3] W. Chiang, R. Russell, A reactive tabu search metaheuristic for the vehicle routing problem with time windows, INFORMS Journal on Computing 9 (1997) 417–430. [4] J. Cordeau, G. Desaulniers, J. Desrosiers, M. Solomon, F. Soumis, The VRP with time windows, in: Paolo Toth, Daniele Vigo (Eds.), The Vehicle Routing Problem, SIAM Monographs on Discrete Mathematics and Applications, 2002 pp. 157–193 (Chapter 7). [5] J. Cordeau, G. Laporte, A. Mercier, A unified tabu search heuristic for vehicle routing problems with time windows, Technical report CRT-00–03, Center for Research on Transportation, Montreal, Canada, 2000. [6] J. Desrochers, J. Lenstra, M. Savelsbergh, F. Soumis, Vehicle routing with time windows: Optimization and approximation, in: Vehicle Routing: Methods and Studies, Elsevier Science Publishers, North-Holland, 1988. [7] M. Desrochers, J. Desrosiers, M. Solomon, A new optimization algorithm for the vehicle routing problem with time windows, Operations Research 40 (1992) 342–354. [8] J. Desrosiers, Y. Dumas, M. Solomon, F. Sormis, Time constrained routing and scheduling, in: Handbooks in Operations Research and Management Science: Network Routing, Elsevier Science Publishers, Amsterdam, 1995. [9] Y. Dumas, J. Desrosiers, E. Gelinas, M. Solomon, An optimal algorithm for the traveling salesman problem with time windows, Operations Research 43 (1995) 367–371. [10] M. Fisher, K. J€ ornsten, O. Madsen, Vehicle routing with time windows: Two optimization algorithms, Operations Research 45 (1997) 488–492. [11] J. Homberger, H. Gehring, Two evolutionary metaheuristics for the vehicle routing problem with time windows, INFOR 37 (1999) 297–318. [12] N. Kohl, O. Madsen, An optimization algorithm for the vehicle routing problem with time windows based on lagrangian relaxation, Operations Research 45 (1997) 395– 406. [13] A. Kolen, A. Rinnooy, H. Trienekens, Vehicle routing with time windows, Operations Research 35 (1987) 266–273. [14] H. Li, A. Lim, A metaheuristic for the pickup and delivery problem with time windows, in: ICTAI 2001, November 2001, pp. 151–158. [15] J. Potvin, S. Bengio, The vehicle routing problem with time windows-part 2: Genetic search, INFORMS Journal on Computing 8 (1996) 165–172. [16] J. Potvin, B. Garcia, L. Rousseau, The vehicle routing problem with time windows––Part 1: Tabu search, INFORMS Journal on Computing 8 (1996) 158–164. [17] J. Potvin, J. Rousseau, A parallel route building algorithms for the vehicle routing and scheduling problem with time windows, European Journal of Operational Research 66 (1993) 331–340. [18] J. Potvin, J. Rousseau, An exchange heuristic for routing problems with time windows, Journal of Operational Research Society 46 (1995) 1433–1446.
H. Li, A. Lim / European Journal of Operational Research 150 (2003) 115–127 [19] Y. Rochat, F. Semet, A tabu search approach for delivering pet food and flour in switzerland, Journal of Operational Research Society 45 (1994) 1233–1246. [20] Y. Rochat, E. Taillard, Probabilistic diversification and intensification in local search for vehicle routing, Journal of Heuristics 1 (1995) 147–167. [21] L. Rousseau, M. Gendreau, G. Pesant, Using constraintbased operators to solve the vehicle routing problem with time windows, Journal of Heuristics 8 (2002) 43–58. [22] R. Russell, Hybrid heuristics for the vehicle routing problem with time windows, Transportation Science 29 (1995) 156–166. [23] M. Savelsbergh, Local search for routing problems with time windows, Annals of Operations Research 4 (1985) 285–305. [24] M. Savelsbergh, The vehicle routing problem with time windows: Minimizing route duration, ORSA Journal on Computing 4 (1992) 146–154.
127
[25] M. Solomon, Algorithms for the vehicle routing and scheduling problems with time window constrains, Operations Research 35 (1987) 254–264. [26] E. Taillard, P. Badeau, M. Gendreau, F. Geurtin, J. Potvin, A tabu search heuristic for the vehicle routing problem with time windows, Transportation Science 31 (1997) 170–186. [27] S. Thangiah, Vehicle routing with time windows using genetic algorithms, Technical report SRU-CPSC-TR-93– 23, Computer Science Department, Slippery Rock University, Slippery Rock, PA, 1993. [28] S. Thangiah, I. Osman, T. Sun, Hybrid genetic algorithm, simulated annealing and tabu search methods for vehicle routing problems with time windows, Technical report UKC/OR94/4, Institute of Mathematics and Statistics, University of Kent, Canterbury, UK, 1994.
