Active guided evolution strategies for large-scale vehicle routing problems with time windows

Computers & Operations Research 32 (2005) 1593 – 1614

www.elsevier.com/locate/dsw

Active guided evolution strategies for large-scale vehicle routing problems with time windows David Mestera;∗ , Olli Br,aysyb a

Institute of Evolution, Mathematical and Population Genetics Laboratory, University of Haifa, 31905 Haifa, Israel b SINTEF Applied Mathematics, Department of Optimization, P.O. Box 124 Blindern, N-0314 Oslo, Norway

Abstract We present a new and e/ective metaheuristic algorithm, active guided evolution strategies, for the vehicle routing problem with time windows. The algorithm combines the strengths of the well-known guided local search and evolution strategies metaheuristics into an iterative two-stage procedure. More precisely, guided local search is used to regulate a composite local search in the 2rst stage and the neighborhood of the evolution strategies algorithm in the second stage. The vehicle routing problem with time windows is a classical problem in operations research, where the objective is to design least cost routes for a 4eet of identical capacitated vehicles to service geographically scattered customers within pre-speci2ed time windows. The presented algorithm is speci2cally designed for large-scale problems. The computational experiments were carried out on an extended set of 302 benchmark problems. The results demonstrate that the suggested method is highly competitive, providing the best-known solutions to 86% of all test instances within reasonable computing times. The power of the algorithm is con2rmed by the results obtained on 23 capacitated vehicle routing problems from the literature. ? 2003 Elsevier Ltd. All rights reserved. Keywords: Heuristics; Vehicle routing; Time windows; Evolution strategies; Guided local search

1. Introduction Vehicle routing problems are important and well-known combinatorial optimization problems occurring in many transport logistics and distribution systems of considerable economic signi2cance. The vehicle routing problem with time windows (VRPTW) has recently received a lot of attention in the literature. This is mainly because of the wide applicability of time window constraints in real-world cases. In the VRPTW, the objective is to 2nd a set of minimum-cost vehicle routes ∗

Corresponding author. E-mail address: [email protected] (D. Mester).

0305-0548/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2003.11.017

1594

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

which start at a central depot, service a set of customers with known demands, and return to the depot. Each customer must be serviced once by a vehicle, and the total demands of the customers serviced by the vehicle must not exceed the capacity of the vehicle. Moreover, each customer must be serviced within a speci2ed time window. If a vehicle arrives at a customer earlier than the lower bound of the customer’s time window, the vehicle must wait until the service is possible. In some applications, service after the upper bound of the time window is also allowed, but a penalty is set for the service (e.g. Koskosidis et al. [1] and Taillard et al. [2]). The depot has also a time window, and all the vehicles must return by the closing time of the depot. The objective is to minimize the number of tours or routes, and then for the same number of routes, to minimize the total traveled distance. VRPTW is an NP-hard problem, and although optimal solutions can be obtained using exact methods, the computation time required to solve the VRPTW to optimality is prohibitive (Desrochers et al. [3]). Since heuristic methods often produce near optimal solutions in a reasonable amount of computational time, most of the research has focused to the design of heuristics and metaheuristics. For surveys on exact, heuristic and metaheuristic methods, we refer to Desrosiers et al. [4], Cordeau et al. [5], and Br,aysy and Gendreau [6,7], respectively. The heuristic methods for the VRPTW can be divided into construction heuristics, improvement heuristics (local searches) and metaheuristics. Construction heuristics build a feasible solution by inserting unrouted customers iteratively into current partial routes according to some speci2c criteria, such as minimum additional distance or maximum savings, until the route’s scarce resources (e.g. capacity) are depleted. Sequential route construction heuristics build routes one at a time whereas parallel construction heuristics create several routes simultaneously. For examples of construction heuristics, see insertion heuristics of Solomon [8], Potvin and Rousseau [9], Bramel and Simchi-Levi [10], Mester [11], Ioannou et al. [12] and Dullaert and Br,aysy [13]. Improvement heuristics modify the current solution iteratively by performing local searches for better neighboring solutions. Generally, a neighborhood comprises the set of solutions that can be reached from the incumbent solution by replacing a set of k arcs or customer nodes. For the most successful applications of improvement heuristics for the VRPTW, see Thompson and Psaraftis [14], Potvin and Rousseau [15], Russell [16], Shaw [17], Caseau and Laburthe [18], Schrimpf et al. [19], Cordone and Wol4er Calvo [20] and Br,aysy et al. [21]. Metaheuristics are general solution procedures that explore the solution space to identify good solutions and they often embed some of the standard route construction and improvement heuristics. In contrast to classical approaches, metaheuristics allow deteriorating and even infeasible intermediate solutions in the course of the search process to escape local minima. So far tabu searches [22], genetic algorithms [23] and evolution strategies [24] have shown the best performance for vehicle routing problems. For example Rochat and Taillard [25], Taillard et al. [2], Chiang and Russell [26], Schulze and Fahle [27] and Cordeau et al. [28] have applied tabu search, while Berger et al. [29] and Jung and Moon [30] report good results with genetic algorithms. An agent-based methodology drawing on both genetic algorithms and tabu search is introduced in Le Bouthillier and Crainic [31]. Recently Homberger and Gehring [32–34] and Mester [11,35] have obtained excellent results with evolution strategies algorithms. For other recent eNcient metaheuristics for the VRPTW, see for example Kontoravdis and Bard [36] who describe a two-phase greedy randomized adaptive search procedure (GRASP) based on a penalty guided random insertion procedure and standard improvement heuristics. Liu and Shen [37] propose a two-stage metaheuristic based on a new neighborhood structure focusing on the

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

1595

relationship between routes and nodes. Guided local search [38] is applied in Kilby et al. [39], and Gambardella et al. [40] introduces an ant colony system [41] consisting of two separate ant colonies with di/erent objectives, guiding well-known construction and improvement heuristics. Br,aysy [42] and Rousseau et al. [43] experiment with variable neighborhood search [44], and simulated annealing [45] is applied in Chiang and Russell [46], Li and Lim [47] and Bent and Van Hentenryck [48]. Genetic algorithms, tabu search and simulated annealing are hybridized in Thangiah et al. [49], and Ibaraki et al. [50] describes di/erent multi-start local search strategies. For further details on the metaheuristics for the VRPTW, we refer to Br,aysy and Gendreau [7]. According to Br,aysy and Gendreau [7] the previously proposed heuristics for the VRPTW show signi2cant variability in performance, and fail to convincingly provide a single robust and successful solution technique. Several recent studies of real-life VRPTW applications indicate a signi2cant potential for improving existing optimization techniques (Hasle [51]). Given the practical importance of the VRPTW, further research on this topic is therefore justi2ed. The main contribution of this paper is development of a more robust and eNcient solution technique for large scale VRPTW. The novelty of the new metaheuristic concept consists of combining the strengths of guided local search (GLS) [38] and evolution strategies (ES) [24] into an iterative two-stage procedure. More precisely, GLS is used to regulate a composite local search in the 2rst stage, and the objective function and the neighborhood of the modi2ed ES local search algorithm of Mester [11,35] in the second stage. We call the suggested new metaheuristic active guided evolution strategies (AGES). To our best knowledge this kind of method has not been used before. Our experiments on 302 large-scale VRPTW benchmarks and on 23 capacitated vehicle routing problem (CVRP) instances from the literature demonstrate that the proposed algorithm is fast, cost-e/ective and highly competitive. It found the best-known solutions to 86% of the well-known large-scale VRPTW benchmark problems and to 21 out of the 23 tested CVRP instances. The remainder of this paper is outlined as follows. Sections 2 and 3 introduce the main concepts of evolution strategies and guided local search metaheuristics. Then, the di/erent components of the AGES procedure are described in Section 4, including general principles, initial solution procedure and improvement heuristics. Section 5 presents the results of computational experiments to assess the value of the proposed approach and reports a comparative performance analysis to alternate methods. Finally, in Section 6 conclusions are drawn. 2. Evolution strategies ES were founded in 1970s by Rechenberg [24] and Schwefel [52] to solve optimization problems with real-value variables. Together with genetic algorithm [23] and evolutionary programming [53], evolution strategies form the class of evolutionary algorithms [54]. Evolutionary algorithms (EAs) are heuristic search techniques that mimic biological evolution and natural selection. In general, an EA [55,56] maintains a population of individuals that undergoes a selection process for mating purposes, and is manipulated by genetic operators. Each individual represents a possible solution to the problem at hand. Here one must note that an individual is not a decision vector but rather encodes it to a chromosome based on an appropriate structure, e.g. a bit vector or a real-valued vector. Selection is a process in which solutions are selected for recombination (parents) based on their 2tness values. Here 2tness refers to a measure of pro2t, utility or goodness to be maximized while exploring the solution

1596

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

space. Recombination (or crossover) and mutation are genetic operators aiming at generating new solutions within the search space from existing ones. The crossover operators use a certain number of parents to create a certain number of children (o/spring) by recombining the parents. By contrast, the mutation operator modi2es individuals by changing (typically) small parts in the associated vectors according to a given mutation rate. Based on the above concepts, natural evolution is simulated by an iterative computational process. In the beginning a population of candidate solutions to a problem at hand is initialized. This is often accomplished by randomly sampling from the space of possible solutions. Then a loop consisting of evaluation (2tness assignment), selection, recombination and/or mutation is executed a certain number of times. Each loop iteration is called a generation, and the search is stopped if the convergence criteria or conditions are met. Such criteria might, for instance, refer to a maximum number of generations or the convergence to a homogeneous population composed of similar individuals. Two important features separate ES from the other evolutionary algorithms. First, ES use a real coding of the decision vector, and model the organic evolution at the level of individual’s phenotypes. Second, ES depend on deterministic selection of the individuals to the next generation, and the search is mainly driven by a high mutation rate. An individual in ES is represented as a pair of real vectors, v = (x; ). The 2rst vector, x, represents a point in the search space (solution) and consists of a number of real-valued variables. The second vector, , represents so-called strategy parameters. In the real-value case, the strategy parameters represent standard deviations of normally distributed random variables. Mutation is performed by replacing x by xt+1 = xt + N (0; ) where N (0; ) is a random Gaussian number with a mean of zero and standard deviation . For discrete combinatorial optimization problems such as the VRPTW, the above-described ES for real-valued optimization problems has to be modi2ed in three ways. First, the representation is modi2ed to consider a set of discrete elements. The strategy parameters  are used to adjust the mutation step size, but now contain encoded information on the transition from a given solution S to a solution S’ in the neighborhood N (S) of S. Third, the neighborhood N (S) of a solution S that is searched through mutation consists here of a set of feasible solutions that can be generated from S by changing one or several attributes of S in one step (iteration). The evolution strategies can be categorized according to the population size ( ) and the number of o/spring ( ) created in each generation (often ¿ ¿ 1). As in Mester [11,35] we use a single individual ( = 1) to create a single o/spring ( = 1), a so-called (1+1)-evolution strategy. The (1+1)-evolution strategy does not apply recombination, but generates the o/spring exclusively through mutation without considering the vector of strategy parameters . The size of the mutation step is random, i.e., in the mutation an attempt is made to modify the incumbent (parent) solution for a given number of times, and the number and type of modi2cations done each time are random. The modi2cations are done one at a time, and after each modi2cation, the modi2ed solution (o/spring) is evaluated according to the 2tness value. If the 2tness of the o/spring is better than the 2tness of the parent, the o/spring replaces the parent solution. For more information we refer to B,ack et al. [56], Homberger and Gehring [32] and Mester [11,35]. 3. Guided Local Search Guided local search (GLS), developed by Voudouris [38] and Voudouris and Tsang [57], is a conceptually simple memory-based metaheuristic that includes some properties of the well-known

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

1597

tabu search and simulated annealing algorithms. It has been shown to be e/ective for numerous combinatorial optimization problems, such as partial constraint satisfaction [58] and routing and scheduling problems [59,39,60]. GLS operates by penalizing particular solution features (usually one) it considers should not occur in a near-optimal solution (e.g. long arcs in a routing problem) weighted by the number of times the feature has already been penalized. The more often a feature appears and is penalized, the less likely it is to be penalized further. GLS associates a cost and penalty to each feature. The cost refers often to the values of the original objective function, and the penalties are initially set to 0 and are updated during the search each time a local minimum is reached. Given an objective function g that maps every candidate solution s to a numerical value, GLS replaces g with a new objective function h
h(s) = g(s) +

(pi Ii (s)); (1) where s is a candidate solution, is a parameter controlling the impact of penalties, pi is the penalty counter for feature i, which holds the number of times i has been penalized, and Ii is a indicator function, which tests whether s exhibits feature i: Ii (s) = 1 if s exhibits feature i; and 0 otherwise:

(2)

When a search is trapped in a local optimum, the goal is to systematically fade out features that contributed most to the current state (the local optimum) through penalizing them. However, to avoid constantly penalizing the most signi2cant features, one must consider the number of times each present feature was penalized in the past. For this purpose we use the utility function (2) such that in the local minimum S ∗ the feature i with maximum Ui is selected to be penalized. Ii (s∗ )ci ; (3) Ui (s∗ ) = (1 + p1 ) where ci is the cost and pi is the current penalty value of feature i. In other words, if a feature is not exhibited in the local optimum the utility of penalizing it is zero, and the higher the cost of the feature ci , the greater the utility of penalizing it. In a local optimum the penalty of feature k with the greatest Uk value is modi2ed as follows: pk = pk + 1:

(4)

By taking the cost and the current penalty into consideration in selecting the feature to penalize, we are distributing the search e/ort in the search space. Candidate solutions that exhibit “good features”, i.e., features involving lower cost, are given more e/ort in the search, but penalties help to prevent all e/orts to be directed to the best features. Here we selected arcs as the feature to penalize. 4. Active guided evolution strategies This section 2rst describes the algorithm and the general principles of the AGES metaheuristic. Then the initial solution procedure as well as the two improvement stages of AGES are presented in detail.

1598

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

4.1. The search strategy The suggested metaheuristic for the VRPTW consists of two phases. First an initial solution is created with the hybrid insertion heuristic of Mester [11], followed by an attempt to minimize the number of routes with a modi2ed version of the 2lling procedure of Bent and Van Hentenryck [48]. The procedures are explained in more detail in Section 4.2. In the second phase an attempt is made to improve the created solution with the proposed new AGES metaheuristic. AGES combines the ES and GLS metaheuristics described in the previous sections into an iterative two-stage procedure. The algorithm is: Phase 1: Step 1.1. Step 1.2. Step 1.3. Step 1.4. Phase 2: Step 2.1. Step 2.2. Step 2.3. Step 2.4. Step 2.5. Step 2.6. Step 2.7. Step 2.8. Step 2.9. Step 2.10. Step 2.11. Step 2.12. Step 2.13.

Create an initial solution S 0 with hybrid insertion heuristic of Mester [3] (see Section 4.2) Initialize penalty parameters and set c1 = 0, c2 = 0, i = 0 While the number of routes in S 0 is greater than kS or until no more routes can be eliminated Apply the modi2ed 2lling procedure (see Section 4.2) to S 0 Set the global best solution S ∗∗ = S 0 2 Re-initialize penalty parameters, set i = 0, c2 = 0, S i = S ∗∗ and de2ne cmax randomly between 5000 and 50000 De2ne the currently penalized arc in S i and recalculate the penalties De2ne the current PVN (see Section 4.3.1) within S i Create a new solution S ∗ using composite local search (see Section 4.3.1) within PVN of S i If S ∗ is better than S ∗∗ Set S ∗∗ = S ∗ If S ∗ is better than S i according to GLS-based objective function Set i = i + 1, c1 = 0, c2 = 0, S i = S ∗ and goto step 2.2. 1 If c1 ¡ cmax or if size of PVN ¡ 25 1 Set c = c1 + 1, c2 = c2 + 1 and goto step 2.2. Create a new solution S ∗ by applying the ES local search (see Section 4.3.2) within the PVN of S i If S ∗ is better than S ∗∗ Set S ∗∗ = S ∗ If S ∗ is better than S i according to GLS-based objective function Set i = i + 1; c2 = 0; S i = S ∗ and goto step 2.8. If the stopping criterion is met Stop the search and return the global best solution S ∗∗ Set c1 = 0; c2 + c2 + 1 2 If c2 ¡ cmax Goto step 2.2. Else Goto step 2.1

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

1599

The 2rst stage (steps 2.2–2.7) consists of a composite local search procedure, where the three well-known improvement heuristics relocate [61], 1-interchange [62] and 2-opt* [15] are applied cyclically (2-opt* is applied only to the CVRP problem instances, see Section 5.2). The composite local search is guided by the above-described GLS by penalizing long arcs appearing often in local 1 minima. If no more improvements have been found for a user-de2ned number of iterations (cmax =25), the second stage (steps 2.8–2.10) is started. The local search procedure of the second stage is based on the modi2ed ES local search procedure (mutation) of Mester [11]. The local search procedure is similar to the large neighborhood search (LNS) of Shaw [17], and here the neighborhood considered in the local search process is determined by the GLS. Moreover, the GLS-guided objective function of the 2rst stage is used also in the second stage. The second stage is continued until no more improvements can be found and then the search goes back to the 2rst stage. The two stages are repeated iteratively, and the search is stopped by a user if no more improvements can be found (step 2.11.). As deteriorating is allowed, the algorithm maintains in memory the best solution found during the entire search, S ∗∗ , (the global best solution) and returns that solution in the end. The best solution found is used as a restart point (steps 2.1. and 2.13.) if the best solution has not been improved for a randomly de2ned number of trials. After the restart the penalty counters and penalties are set to their original values. In the experimental tests the multi-restart mechanism was found to have a signi2cant role in 2nding better solutions and speeding up the algorithm. The improvement procedures are described in more detail in Section 4.3. 4.2. Creation of the initial solution The hybrid insertion heuristic of Mester [11] is used to create the initial solution. It draws upon the savings heuristic of Clarke and Wright [63] and the composite heuristic of Russell [16]. A solution is represented as a vector of customer indices grouped into separate routes. As in Clarke and Wright [63] the heuristic begins with a solution in which each customer is supplied individually by a separate route. The cheapest reinsertions of single customers are then repeated as long as one can reduce the number of routes. The construction heuristic is of order O(n) and for a customer k currently serviced between customers i and j in the solution vector it considers only positions between the adjacent customers l and m, k ¡ l ¡ m 6 n using a modi2ed version of the insertion criterion of Osman [62] (dik + dkj − dlm ) − (dlk + dkm − dij );

(5)

where n refers to the number of points in the solution vector, dij is the distance between the customers at positions i and j, and  is a parameter whose values depend on the problem size. For the CVRP and 200-customer VRPTW instances all values within the range 0.6 –1.4 are tried in increments of 0.2 units. Similarly, for the 400 – 600-customer problems the range is 0.8–1.2 and for 800 –1000-customer problems only value 1 is tried. Therefore, 5, 3 and 1 initial solution(s) are created in total for problem sizes of 200, 400 – 600 and 800 –1000 customers, respectively. Each time 20% of the customers has been inserted, the relocate [61] and 1-interchange [62] improvement procedures are applied cyclically until no improvement can be found. The basic idea of the relocate operator is to reinsert a single customer in another position. Here the procedure does both intra- and inter-route reinsertions, i.e., it considers insertion positions both in the selected customer’s current route and in an alternate route. 1-interchange swaps simultaneously the position of two customers

1600

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

in two separate routes. The best solution obtained with di/erent values of  is selected in the end as the initial solution. Then an attempt is made to reduce the number of routes in the created initial solution. For reducing the number of routes, a reinsertion procedure using the 2lling strategy of Bent and Van Hentenryck [48] is used. The procedure is applied until no more improvements can be found or S Here kS is set equal to the the number of routes in the current solution is equal to a given value k. number of routes in the best-known solutions reported in earlier papers. The use of the best-known number of routes can be justi2ed by the fact that in practice the number of routes is often known, i.e., since the company has a limited number of vehicles and because we assume that each route is served by a separate vehicle. However, when the number of routes/vehicles is not known in advance, there are also several ways to handle the problem algorithmically. For example one can use another (meta)heuristic to estimate the number of routes or use lower bounds e.g. based on solving Bin Packing Problems [37]. The basic idea of the 2lling procedure is to relocate customers such that the sum of the squares of the route sizes |r| is maximized. The procedure reinserts some points from routes with smaller |r| value to the routes with greater |r| value. We modi2ed the original 2lling procedure of Bent and Van Hentenryck [48] by using also a random 1-interchange procedure [62] with 0.1 probability. This often worsens the solutions according to the distance criterion, but computational testing revealed that it sometimes creates a better basis for reducing the number of routes. 4.3. Improvement phase The improvement phase is divided in two stages that are applied iteratively. In the 2rst stage a simple GLS is used to guide a composite local search consisting of the well-known relocate [61], 1-interchange [62] and 2-opt* [15] improvement heuristics and determine a restricted neighborhood for the improvement procedures. The second stage is based on a modi2ed version of the ES local search of Mester [35]. The two stages are described in the next two subsections. 4.3.1. GLS composite local search The GLS composite local search procedure is based on the standard improvement heuristics relocate [61], 1-interchange [62] and 2-opt* [15]. These operators are applied cyclically with the best-accept strategy such that each operator is applied only once at a time, and then the neighborhood is changed. Here one must note that 2-opt* is applied only for the CVRP problem instances. The composite local search is guided by the above-described simple GLS by penalizing an arc with maximum utility according to Eq. (3) every time the algorithm gets stuck in a local minimum. Here the cost in Eq. (3) refers to the length of the arc, and the penalty parameter is set to = 0:05 Moreover, the neighborhood considered by the composite local search is restricted to the set of geographically closest routes to the currently penalized arc. This restricted neighborhood is called a penalty variable neighborhood (PVN). Its size is adjusted dynamically during the search such that the minimum number of the routes in the PVN (Nmin ) equals three. The maximum size of the PVN (Nmax ) is limited by the minimum of the number of customers in the selected close routes and 60% of the problem size. The current size N k of the PVN is determined randomly in the beginning of every iteration N k = Nmin + (Nmax − Nmin )a2 ;

(6)

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

1601

where a is a random value between 0 and 1. This de2nition favors small sizes for the PVN. To reduce the number of the routes, we periodically (with probability 10=n) apply a PVN that consists of all customers. If no improvements have been found for a user-de2ned number of iterations (25), the 2rst stage is stopped and the second stage (see Section 4.3.2.) is started. The PVN is used also in the second stage to restrict the neighborhood considered, and to control launching the second stage: the second stage is started only if the current size of the PVN is bigger than 25. Thus, the ES local search is not used in every iteration. In case of large PVN (more than 100 customers or more than 5 routes), the PVN is decomposed in the second stage into sets of two di/erent routes. After the decomposition each pair of routes make up their own PVN and there are k(k − 1)=2 PVNs, where k is the number of routes in a solution. These two-route PVNs are then considered in a loop that is restarted from the 2rst pair every time an improvement is found. 4.3.2. ES local search The ES local search of Mester [11,35] draws upon the large neighborhood search of Shaw [17] and the remove–insert scheme of Mester [11]. The basic idea is to remove a selected set of customers from the current solution and then reinsert the removed customers at optimal cost. Three strategies are used to select the customers to be removed. The customers to be removed are selected from the current PVN. In the 2rst two strategies the number of removed customers are also limited to = (0:2 + 0:5a)n;

(7)

where n is the number of customers in the current PVN and a is random value (uniformly distributed) between 0 and 1. The 2rst strategy randomly selects the customers from the PVN. The second strategy randomly ejects customers from the PVN within a ring created by two circles with random radiuses, centered at the depot. The third strategy removes all customers from the current PVN. After the removal, the selected customers are reinserted using the hybrid insertion heuristic of Mester [11] described in Section 4.2 with the modi2ed objective function determined by the GLS in the 2rst stage, varying the value of parameter  from 0.2 to 1.4 in increments of 0.2 units. Once all removed customers have been inserted, the solution is improved with the GLS composite local search procedure. The ES local search is repeated for 14 iterations at a time. During the 2rst 7 iterations only the random and radial removal strategies (strategies 1 and 2) are applied randomly with 50% probability. In the last 7 iterations, the algorithm uses only the third removal strategy to remove all customers from the current PVN. During both sets of 7 iterations, the value of parameter  is varied from 0.2 to 1.4 in increments of 0.2 units. In each iteration, a comparison is made between the o/spring and the parent solution using the GLS-based objective function, and if the o/spring is better than the parent, it replaces the parent solution. The second stage is continued until no more improvements can be found according to the GLS-based biased objective function (2tness value), and then the search goes back to the 2rst stage. If the algorithm cannot improve the current solution at all in the second stage, the currently penalized arc is determined randomly instead of heuristic rule based on Eq. (3) in the 2rst stage to protect the algorithm from cycling. Compared to the 2rst stage, the ES local search performs signi2cantly more thorough search, but it requires higher executing time

1602

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

(factor 15 –20). Therefore ES is applied only in iterations where the randomly determined size of the PVN is large (more than 25 customers). 5. Experimental results A computational experiment has been conducted to analyze the performance of the proposed algorithm. In this section we describe the test data sets and present the results of the suggested AGES metaheuristic along with a comparative analysis with the state-of-the-art methods from the literature. 5.1. Problem data and the experimental setting As the focus is on large-scale problems, the algorithm has been tested on the 300 extended VRPTW benchmarks of Gehring and Homberger [33] and two real-life problems of Russell [16]. To demonstrate the robustness of our algorithm, we experimented also with 23 well-known CVRP benchmark problems introduced by Christo2des et al. [64], Fisher [65] and Taillard [66,25]. The problem description of the CVRP is the same as for the VRPTW except for the fact that there are no time window constraints. We focused on problems with only capacity constraints and used the same parameter setting as described earlier for the VRPTW. An earlier version of our algorithm has been applied successfully also to the well-known traveling salesman problem and some genetic mapping problems. For details, we refer to Mester et al. [67]. The Gehring and Homberger test instances consist of 200, 400, 600, 800 and 1000 customers. For each problem size, the 60 problems are divided into R1, C1, RC1, R2, C2 and RC2 groups. Travel times separating two customers corresponds to their relative Euclidean distance. The problems vary in 4eet size, vehicle capacity, spatial and temporal distribution of customers, time window density and width and customer service times. In problems belonging to the C group, the customers are clustered geographically or based on their time windows. In the R group the customers are uniformly distributed, and the RC group combines clustered and randomly distributed customers. Problem sets R1, C1 and RC1 have a narrow scheduling horizon. Hence, only a few customers can be serviced by the same vehicle. Conversely, problem sets R2, C2 and RC2 have a larger scheduling horizon, and more customers can be served by the same vehicle. The two benchmark instances D417 and E417 of Russell [16] are extracted from a real-life fast-food routing application in the Southeastern United States. There are 417 customers in both problems and the problems di/er only in that the problem E417 has a higher percentage of tight time windows. The CVRP instances vary also in problem and 4eet size, vehicle capacity and spatial distribution of customers. The number of customers in each problem is described in the second column of Table 1. In the table the problems of Christo2des et al. [64] are listed 2rst with notations C1–C5 and C11–C12. In the 2rst 2ve problems the customers are randomly distributed, whereas problems C11 and C12 are structured in the sense that customers appear in clusters. For more details on these problems we refer to Christo2des et al. [64] and Reimann et al. [68]. The three benchmarks of Fisher [65] are taken from real vehicle routing applications. More precisely, F44 and F134 represent a day of grocery deliveries in Ontario and F71 is concerned with the delivery of tires, batteries and accessories to gasoline stations. The problem T385 described in Taillard [66] is also based on

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

1603

Table 1 The results for the capacitated vehicle routing benchmarks Problem

Problem size

Best published

C1 C2 C3 C4 C5 C11 C12

50 75 100 150 199 120 100

524.61 835.26 826.14 1028.42 1291.45 1042.11 819.56

F44 F71 F134

44 71 134

T75a T75b T75c T75d T100a T100b T100c T100d T150a T150b T150c T150d T385

AGES

% above best-known

CPU/s

524.61 835.26 826.14 1028.42 1291.29* 1042.11 819.56

0 0 0 0 −0.01 0 0

0.4 13 40 400 39 575 12 1.5

723.54 241.97 1162.96

723.54 241.97 1162.96

0 0 0

75 75 75 75 100 100 100 100 150 150 150 150

1618.36 1344.64 1291.01 1365.42 2041.34 1940.61 1406.20 1581.25 3055.23 2727.77 2341.84 2645.40

1618.36 1344.64 1291.01 1365.42 2041.34 1939.90* 1406.20 1581.25 3055.23 2727.67* 2343.11 2645.40

385

24431.44

24855.32

0 0 0 0 0 −0.04 0 0 0 −0.004 0.05 0 1.73

0.08 2 24 2 1 5 3 130 16 33 81 60 420 2 13 900 840

real-life data. The problem was generated by selecting the most important towns or villages in a canton in Switzerland and by basing the demands on the number of inhabitants in each commune. The T75, T100 and T150 problem sets are introduced in Rochat and Taillard [25]. The customers in these problems are spread in several clusters with variable size and compactness and the quantities ordered by the customers are exponentially distributed. We used a hierarchic objective function in our experiments, where the primary objective is to minimize the number of routes. For the same number of routes, the secondary objective is the minimization of the total traveled distance. The proposed algorithm has been implemented in Visual Basic 6.0 and it includes a graphical user interface that both have a considerable negative impact on the computation time compared to an eNcient C++ implementation. The experimental tests consist only of one simulation run for each problem with a Pentium IV 2 GHz (512 Mb RAM) computer. Because of limited computational resources, parameter values described in the previous sections were determined experimentally over a few intuitively selected combinations, choosing the one that yielded the best average output. This is justi2ed by the fact that the sensitivity of the results with respect to changes in the parameter values was found to be generally quite small.

1604

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

5.2. Results for the capacitated vehicle routing benchmark problems The results for the 23 CVRP instances are presented in Table 1, in which the 2rst column to the left describes the problem names. The second column gives the number of customers (or cities) in each problem, as depicted earlier. The best published solutions are taken from Reimann et al. [68] and from the web page of Taillard (http://ina.eivd.ch/collaborateurs/etd/problemes.dir/vrp.dir/vrp.html). The values in Table 1 refer to the total distance over all routes in a solution that is often considered as the only objective in CVRP. Our results over a single test run are shown in the AGES column. Next to the absolute distance values the table indicates the di/erence of our results with respect to the best-known solutions in percentages. Finally, the rightmost column presents the CPU time in seconds. According to Table 1, AGES obtains the best-known solution to 21 out of the 23 test instances. For C5, T100b and T150b we found a new best-known solution. Correspondingly, to T150c and to T385 AGES did not 2nd the best-known solution. The main reason for this is that we used the same algorithm to solve both CVRP and VRPTW where the primary objective is to minimize the number of routes. In these two cases AGES found a solution that has a smaller number of routes than in the best-known solution, but a slightly higher total distance value. Here one must note that AGES was able to 2nd a solution with smaller number of routes than in the best-known solution also to problem T75b. More precisely, in the best-known solution to T75b there are 10 routes. AGES found a solution that requires 9 routes and whose total distance is 1355.15. On the average our results are within 0.075% from the best-known solutions and the computing times are reasonable. As one can see from Table 1, the computing times vary a lot over the di/erent problems examined. The reason for this is the stopping criterion where the search in continued as long as the algorithm can 2nd improvements. In some problems AGES seems to be able to continue improving the solutions for a very long time. When comparing our results with the results of the best-known methods presented in Reimann et al. [68], one can conclude that our method is competitive. 5.3. Results for large-scale VRPTW benchmark problems In this section we describe our experimental results for the large-scale VRPTW instances of Gehring and Homberger [33] and Russell [16] and present a comparative analysis with the state-ofthe-art methods from the literature. 5.3.1. Results for the benchmark problems of Gehring and Homberger The 2nal results of our experimental investigation on the 300 VRPTW benchmarks of Gehring and Homberger [33] are shown in Tables 2–6 in a uniform way. Each table presents results for a given problem size. The results obtained with AGES are compared with the state-of-the-art results from the literature. The authors described in the 2rst row of Tables 2–6 are: Gehring and Homberger [33] (GH99), Gehring and Homberger [34] (GH01), Bent and Van Hentenryck [48] (BVH), Li and Lim [46] (LL), Le Bouthillier and Crainic [31] (LC) and Br,aysy et al. [21] (BHD). For each dataset (R1, R2, C1, C2, RC1 and RC2), we report the average number of vehicles and the average total distance with respect to all problems in that data set. The notations CNV and CTD

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

1605

Table 2 The results for 200-customer benchmark problems

R1 R2 C1 C2 RC1 RC2 CNV CTD Computer CPU (min) Simulations

GH99

GH01

BVH

LL

LC

BHD

AGES

18.20 3705 4.00 3055 18.90 2782 6.00 1846 18.00 3555 4.30 2675 694 176 180 Pentium 200 4 × 10 1

18.20 3855.03 4.00 3032.49 18.90 2842.08 6.00 1856.99 18.10 3674.91 4.40 2671.34 696 179 328 Pentium 400 4 × 2:1 3

18.20 3677.96 4.10 3023.62 18.90 2726.63 6.00 1860.17 18.00 3279.99 4.50 2603.08 697 171 715 Sun Ultra 10 n/a n/a

18.30 3736.20 4.10 3023.00 19.10 2728.60 6.00 1854.90 18.30 3385.80 4.90 2518.70 707 172 472 Pentium 545 182.1 3

18.20 3676.95 4.00 2986.01 18.90 2743.66 6.00 1836.10 18.00 3449.71 4.30 2613.75 694 173 061 Pentium 933 5 × 10 1

18.20 3718.30 4.00 3014.28 18.90 2749.83 6.00 1842.65 18.00 3329.62 4.40 2585.89 695 172 406 AMD 700 2.4 3

18.20 3618.68 4.00 2942.92 18.80 2717.21 6.00 1833.57 18.00 3221.34 4.40 2519.79 694 168 573 Pentium 2 GHz 8 1

Table 3 The results for 400-customer benchmark problems

R1 R2 C1 C2 RC1 RC2 CNV CTD Computer CPU (min) Simulations

GH99

GH01

BVH

LL

LC

BHD

AGES

36.40 8925 8.00 6502 38.00 7584 12.00 3935 36.10 8763 8.60 5518 1390 412 270 Pentium 200 4 × 20 1

36.40 9478.22 8.00 6650.28 38.00 7855.82 12.00 3940.19 36.10 9294.99 8.80 5629.43 1392 428 489 Pentium 400 4 × 7:1 3

36.40 8713.37 8.00 6959.75 38.00 7220.96 12.00 4154.40 36.10 8330.98 8.90 5631.70 1393 410 112 Sun Ultra 10 n/a n/a

36.60 8912.40 8.00 6610.60 38.70 7181.40 12.10 4017.10 36.50 8377.90 9.50 5466.20 1414 405 656 Pentium 545 359.8 3

36.50 8839.28 8.00 6437.68 37.90 7447.09 12.00 3940.87 36.00 8652.01 8.60 5511.22 1390 408 281 Pentium 933 5 × 20 1

36.40 8692.17 8.00 6382.63 37.90 7230.48 12.00 3894.48 36.00 8305.55 8.90 5407.87 1391 399 132 AMD 700 7.9 3

36.30 8530.03 8.00 6209.94 37.90 7148.27 12.00 3840.85 36.00 8066.44 8.80 5243.06 1389 390 386 Pentium 2 GHz 17 1

correspond to cumulative number of vehicles and distance, respectively, over all test instances of a given problem size. The lowest row describes the computer, average CPU time, and number of independent simulations for each method, as reported by the authors.

1606

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

Table 4 The results for 600-customer benchmark problems

R1 R2 C1 C2 RC1 RC2 CNV CTD Computer CPU (min) Simulations

GH99

GH01

BVH

LL

LC

BHD

AGES

54.50 20 854 11.00 13 335 57.90 14 792 17.90 7 787 55.10 18 411 11.80 11 522 2082 867 010 Pentium 200 4 × 30 1

54.50 21 864.47 11.00 13 656.15 57.70 14 817.25 17.80 7 889.96 55.00 19 114.02 11.90 11 670.29 2079 890 121 Pentium 400 4 × 12:9 3

55.00 19 308.62 11.00 14 855.43 57.80 14 357.11 17.80 8 259.04 55.10 17 035.91 12.40 11 987.89 2091 858 040 Sun Ultra 10 n/a n/a

55.20 19 744.80 11.10 13 592.40 58.20 14 267.30 18.20 8 202.60 55.50 17 320.00 13.00 11 204.90 2112 843 320 Pentium 545 399.8 3

54.80 19 869.82 11.20 13 093.97 57.90 14 205.58 17.90 7 743.92 55.20 17 678.13 11.80 11 034.71 2088 836 261 Pentium 933 5 × 30 1

54.50 19081.18 11.00 13 054.83 57.80 14 165.90 18.00 7 528.73 55.00 16 994.22 12.10 11 212.36 2084 820 372 AMD 700 16.2 3

54.50 18 358.68 11.00 12 703.52 57.80 14 003.09 17.80 7 455.83 55.00 16 418.63 12.10 10 677.46 2082 796 172 Pentium 2 GHz 40 1

Table 5 The results for 800-customer benchmark problems

R1 R2 C1 C2 RC1 RC2 CNV CTD Computer CPU (min) Simulations

GH99

GH01

BVH

LL

LC

BHD

AGES

72.80 34 586 15.00 21 697 76.70 26 528 24.00 12 451 72.40 38 509 16.10 17 741 2770 1 515 120 Pentium 200 4 × 40 1

72.80 34 653.88 15.00 21 672.85 76.10 26 936.68 23.70 11 847.92 72.30 40 532.35 16.10 17 941.23 2760 1 535 849 Pentium 400 4 × 23:2 3

72.70 33 337.91 15.00 24 554.63 76.10 25 391.67 24.40 14 253.83 73.00 30 500.15 16.60 18 940.84 2778 1 469 790 Sun Ultra 10 n/a n/a

73.00 33 806.34 15.10 21 709.39 77.40 25 337.02 24.40 11 956.60 73.20 31 282.54 17.10 17 561.22 2802 1 416 531 Pentium 545 512.9 3

73.10 33 552.40 15.00 21 157.56 76.30 25 668.82 24.10 11 985.11 72.30 37 722.62 15.80 17441.60 2766 1 475 281 Pentium 933 5 × 40 1

72.80 32 748.06 15.00 21 170.15 76.30 25 170.88 24.20 11 648.92 73.00 30 005.95 16.30 17 686.65 2776 1 384 306 AMD 700 26.2 3

72.80 31 918.47 15.00 20 295.28 76.20 25 132.27 23.70 11 352.29 73.00 30 731.07 15.80 16 729.18 2765 1 361 586 Pentium 2 GHz 145 1

According to Table 2, AGES dominates other approaches for all problem groups except RC2 to which GH99 and LC report a smaller number of vehicles. Regarding the total distance traveled, AGES appears to be clearly the best in all subclasses, and the CTD of AGES is 1.9 – 6.4% better

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

1607

Table 6 The results for 1000-customer benchmark problems

R1 R2 C1 C2 RC1 RC2 CNV CTD Computer CPU (min) Simulations

GH99

GH01

BVH

LL

LC

BHD

AGES

91.90 57 186 19.00 31 930 96.00 43 273 30.20 17 570 90.00 50 668 19.00 27 012 3461 2 276 390 Pentium 200 4 × 50 1

91.90 58 069.61 19.00 31 873.62 95.40 43 392.59 29.70 17 574.72 90.10 50 950.14 18.50 27 175.98 3446 2 290 367 Pentium 400 4 × 30:1 3

92.80 51 193.47 19.00 36 736.97 95.10 42 505.35 30.30 18 546.13 90.20 48 634.15 19.40 29 079.78 3468 2 266 959 Sun Ultra 10 n/a n/a

92.70 50 990.80 19.00 31 990.90 96.30 42 428.50 30.80 17 294.90 90.40 48 892.40 19.80 26 042.30 3490 2 176 398 Pentium 545 606.3 3

92.20 55 176.95 19.20 30 919.77 95.30 43 283.92 29.90 17 443.50 90.00 49 711.36 18.50 26 001.11 3451 2 225 366 Pentium 933 5 × 50 1

92.10 50 025.64 19.00 31 458.23 95.80 42 086.77 30.60 17 035.88 90.00 46 736.92 19.00 25 994.12 3465 2 133 376 AMD 700 39.6 3

92.10 49 281.48 19.00 29 860.32 95.10 41 569.67 29.70 16 639.54 90.00 45 396.41 18.70 25 063.51 3446 2 078 110 Pentium 2 GHz 600 1

than for the other methods. AGES appears to be also very robust, considering that it is a stochastic procedure and the results in Tables 2–6 are based on only one test run. The computational e/orts required by AGES appear to be reasonable and comparable with the other approaches. Here one must note that AGES is implemented with Visual Basic and it uses a graphical user interface that clearly reduce its speed compared to the other methods. On the other hand, AGES found a new best-known solution to 79% of the test problems that partially justi2es high CPU times. AGES found solutions equal to the previous best-known solutions on the average about three times faster than reported in Tables 2–6. As for the other methods, GH99, GH01 and BHD seem to be most eNcient. As Table 3 shows for the 400-customer benchmarks, AGES dominates the other approaches for all other subclasses except RC2. The dominance is only because of lower total distance, except for problem class R1, where the number of vehicles is also lower. For RC2, GH99 and LC are superior in number of vehicles. According to Table 3, AGES obtained the lowest CNV (1389), and the CTD is 2.2–9.8% better than for the other methods. For problem group RC2, AGES is even 3.1% better than the second best method, BHD, in terms of distance. In the 600-customer benchmark cases AGES outperforms the other approaches in subclasses R1, R2, C2 and RC1, as can be seen from Table 4. For C1 and RC2, AGES is dominated regarding the number of vehicles. For C1, GH01 performs best, while the smallest number of vehicles for RC2 is reported in GH99 and LC. Nevertheless, the CNV of AGES is only 0.14% worse than the best value 2079 reported in GH01. On the other hand, AGES dominates all other methods due to lower total distance in all subclasses and the CTD of AGES is 8.9% and 11.8% better than the CTD of GH99 and GH01, respectively. According to Table 5, AGES dominates the other approaches in subclasses R2, C2 and RC2 of 800-customer instances. For R1 and C1, BVH reports the best results, and for RC1 LC seems to

1608

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614 Influence of the Filling procedure (C2-10-1)

Number of routes

50 45 40 35 30 25 0

500

1000

1500

2000

Time (seconds) Filling-ON - - - - - - Filling-Off

Fig. 1. In4uence of the 2lling procedure to the number of routes for problem C2-10-1.

perform best. AGES is competitive in terms of number of vehicles, and the di/erence to the smallest CNV reported in GH01 is only 0.18%. Moreover, AGES is superior to other methods in distance optimization. Table 6 shows that for 1000-customer instances AGES dominates all other methods in subclasses R2, C1, C2 and RC1. For R1 and RC2, GH99 and LC report the best output, respectively. AGES is competitive regarding the number of vehicles, and reports the lowest CNV. AGES also outperforms other approaches due to lower distance values in all problem groups and the CTD of AGES is 2.7–10.2% better than the CTD of the other algorithms. The reason for exceptionally high time consumption of AGES for 1000-customer problems is the same as for some problems in Table 1. In some cases AGES continues to 2nd improvements for a very long time. 5.3.2. Concluding analysis of the results for the Gehring and Homberger benchmark problems The results in Tables 2–6 show that AGES produces highly competitive results regarding both number of vehicles and total distance. The main reason for good performance in minimizing the number of routes is the 2lling strategy inspired by Bent and Van Hentenryck [48]. Its typical behavior and eNciency are illustrated in Fig. 1 on test problem C2-10-1 of Gehring and Homberger [33]. According to Fig. 1, the 2lling strategy helps the algorithm to 2nd a small number of routes considerably faster, and it also a/ects the 2nal output slightly. In other problems the 2lling strategy speeded up the minimization of the number of routes by a factor of 2–300. AGES performs particularly well for distance minimization and outperforms previous approaches in most cases. AGES seems also to scale well. Our approach obtained 259 best-known solutions to the 300 instances of Gehring and Homberger [33]. This equals to about 86% of the extended VRPTW benchmarks. About 7% of these best-known solutions are matches to previously obtained solutions. The results of AGES are presented in detail in Table 7 where the best-known solutions are presented in boldface. For the updated list of best-known solutions to the VRPTW benchmarks, we refer to www.top.sintef.no.

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

1609

Table 7 The detailed results of AGES for the 200 –1000-customer test problems of Gehring and Homberger [33] R1/200

R2/200

C1/200

C2/200

RC1/200

RC2/200

01/20/4761.61 02/18/4054.44 03/18/3382.65 04/18/3067.93 05/18/4112.88 06/18/3599.84 07/18/3151.42 08/18/2963.90 09/18/3784.33 10/18/3307.78

01/4/4501.80 02/4/3645.38 03/4/2932.44 04/4/1981.29 05/4/3381.85 06/4/2914.56 07/4/2453.62 08/4/1849.87 09/4/3111.41 10/4/2657.00

01/20/2704.58 02/18/2949.46 03/18/2708.08 04/18/2644.61 05/20/2702.05 06/20/2701.04 07/20/2701.04 08/18/2769.19 09/18/2642.82 10/18/2649.26

01/6/1931.44 02/6/1863.16 03/6/1777.57 04/6/1720.09 05/6/1879.19 06/6/1857.35 07/6/1849.46 08/6/1820.59 09/6/1830.24 10/6/1806.60

01/18/3691.99 02/18/3298.68 03/18/3025.90 04/18/2879.40 05/18/3419.81 06/18/3393.09 07/18/3266.48 08/18/3115.82 09/18/3083.41 10/18/3038.85

01/6/3103.48 02/5/2827.77 03/4/2617.90 04/4/2055.97 05/4/2912.57 06/5/2620.61 07/4/2550.58 08/4/2317.80 09/4/2175.61 10/4/2015.60

R1/400

R2/400

C1/400

C2/400

RC1/400

RC2/400

01/39/10553.50 02/36/9161.26 03/36/7941.53 04/36/7332.93 05/36/9512.25 06/36/8534.05 07/36/7710.41 08/36/7398.68 09/36/8878.19 10/36/8227.49

01/8/9257.92 02/8/7674.90 03/8/5988.02 04/8/4331.07 05/8/7143.55 06/8/6163.81 07/8/5082.10 08/8/4068.97 09/8/6493.13 10/8/5895.93

01/40/7152.06 02/37/7357.45 03/36/7151.17 04/36/6822.18 05/40/7152.06 06/40/7153.45 07/40/7149.39 08/38/7113.40 09/36/7524.32 10/36/6907.26

01/12/4116.05 02/12/3930.29 03/12/3768.63 04/12/3535.99 05/12/3939.42 06/12/3875.94 07/12/3894.13 08/12/3787.08 09/12/3876.10 10/12/3684.89

01/36/8960.82 02/36/8174.27 03/36/7737.99 04/36/7411.02 05/36/8499.15 06/36/8304.99 07/36/8051.71 08/36/7917.68 09/36/7890.45 10/36/7716.32

01/12/6527.76 02/10/5924.84 03/8/5114.76 04/8/3648.64 05/9/6063.46 06/9/5833.12 07/8/5519.25 08/8/4854.16 09/8/4628.26 10/8/4316.36

R1/600

R2/600

C1/600

C2/600

RC1/600

RC2/600

01/59/21131.09 02/54/19603.70 03/54/17400.60 04/54/15993.80 05/54/20395.00 06/54/18620.26 07/54/17107.91 08/54/15725.86 09/54/19372.96 10/54/18235.57

01/11/18325.60 02/11/15346.42 03/11/11663.06 04/11/8386.64 05/11/15640.60 06/11/12937.47 07/11/10536.84 08/11/8023.64 09/11/13567.84 10/11/12607.09

01/60/14095.64 02/56/14325.96 03/56/13898.99 04/56/13610.66 05/60/14085.70 06/60/14089.70 07/60/14085.70 08/58/14346.81 09/56/13733.56 10/56/13758.19

01/18/7774.10 02/18/7486.88 03/18/7238.70 04/17/7216.45 05/18/7576.35 06/18/7478.63 07/18/7560.53 08/18/7352.42 09/18/7350.94 10/17/7523.34

01/55/17454.39 02/55/16208.24 03/55/15524.33 04/55/15180.72 05/55/17468.57 06/55/17248.87 07/55/16454.79 08/55/16462.49 09/55/16153.00 10/55/16030.86

01/16/13407.75 02/14/11092.35 03/11/9978.25 04/11/7349.88 05/13/11919.72 06/12/11411.08 07/11/11687.04 08/11/10474.95 09/11/10113.82 10/11/9339.41

R1/800

R2/800

C1/800

C2/800

RC1/800

RC2/800

01/80/38324.53 02/72/33548.54 03/72/30151.90 04/72/26838.04 05/72/34741.53 06/72/31737.47 07/72/29538.40 08/72/28342.64 09/72/34231.38

01/15/28440.28 02/15/23335.67 03/15/17992.25 04/15/13625.25 05/15/24611.39 06/15/20697.06 07/15/17058.30 08/15/13053.31 09/15/22588.02

01/80/25030.36 02/76/24678.90 03/73/23978.67 04/72/24040.47 05/80/25166.30 06/80/25160.90 07/79/25518.85 08/76/25379.85 09/73/24713.38

01/24/11654.81 02/24/11422.34 03/23/11554.18 04/23/10963.49 05/24/11432.92 06/24/11357.86 07/24/11397.54 08/24/11206.32 09/24/11249.00

01/73/31590.23 02/73/31286.47 03/73/29772.25 04/73/28414.31 05/73/30454.15 06/73/29674.68 07/73/31778.09 08/73/31227.53 09/73/31746.43

01/20/19989.12 02/17/18099.68 03/15/15116.26 04/15/11392.25 05/16/19105.75 06/15/18882.30 07/15/17461.44 08/15/16529.24 09/15/15823.50

1610

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

Table 7 (continued) R1/800

R2/800

C1/800

C2/800

RC1/800

RC2/800

10/72/31730.45

10/15/21551.26

10/73/27655.00

10/23/11284.46

10/73/31366.51

10/15/14892.29

R1/1000

R2/1000

C1/1000

C2/1000

RC1/1000

RC2/1000

01/100/54145.31 02/92/54102.69 03/91/46621.19 04/91/43461.84 05/92/56068.20 06/91/49059.80 07/91/45847.84 08/91/42767.77 09/91/51391.80 10/91/49348.36

01/19/42977.52 02/19/34928.12 03/19/25693.70 04/19/18886.15 05/19/37276.22 06/19/30735.94 07/19/24399.92 08/19/18191.78 09/19/33787.84 10/19/31725.99

01/100/42478.95 02/93/40647.42 03/90/40934.87 04/90/40410.58 05/100/42469.20 06/100/42471.30 07/100/40998.21 08/96/42170.31 09/92/42221.51 10/90/40894.38

01/30/16879.24 02/29/17228.82 03/29/16367.59 04/29/17153.19 05/30/16587.24 06/30/16371.65 07/31/16578.42 08/30/16399.03 09/30/16651.96 10/29/16178.26

01/90/47143.90 02/90/44906.58 03/90/43782.57 04/90/41917.14 05/90/47632.31 06/90/46391.60 07/90/46157.71 08/90/45585.08 09/90/45405.54 10/90/45041.64

01/22/30701.20 02/20/26899.21 03/18/20843.54 04/18/16815.48 05/19/29399.72 06/18/27003.30 07/18/26202.46 08/18/25604.78 09/18/23582.89 10/18/23582.55

Table 8 Results for the two problem instances by Russell [16] Reference

Thangiah et al. [49] Kontoravdis et al. [36] Rochat et al. [25] Russell [16] Chiang et al. [45] Taillard et al. [2] Chiang et al. [26] Liu et al. [37] Homberger et al. [32] Gehring et al. [34] Br,aysy [42] Br,aysy et al. [21] AGES

Problem D417

Problem E417

Vehicles

Distance

Vehicles

Distance

54 55 54 55 55 55 55 54 54 54 54 54 54

4866 4273.4 6264.80 4964 4232.39 3439.8 3455.28 3747.52 4703 3512.01 3506.21 3387.93 3317.19

55 55 54 55 55 55 55 54 55 54 54 54 54

4149 4985.7 7211.83 6092 4397.49 3707.1 3796.61 4691.14 4732 3832.24 3801.64 3672.73 3542.73

CPU

NeXT 25 MHz, −, 26 min Sun Sparc 10, 5 runs, 11 min Silicon Graphics 100 MHz, −, − PC 486 66 MHz, 3 runs, 7 min PC 486 66 MHz, −, 25 min Sun Sparc 10, −,− Pentium 166 MHz, −, 37 min HP 9000/720, 3 runs, 45 min Pentium 200 MHz, 5 runs, 30 min 4 × Pentium 400 MHz, 5 runs, 142 min Pentium 200 MHz, 1 run, 378 min AMD 700 MHz, 3 runs, 47 min Pentium IV 2 GHz, 1 run, 167 min

5.3.3. Results for the benchmark problems of Russell The results for the two real-life problems D417 and E417 of Russell [16] are shown in Table 8. The 2rst column presents the authors, the next two columns show the number of vehicles and total distance for the two problems, and the last column describes the computer, number of computational runs and the CPU time in minutes. According to Table 8, AGES obtained the best-known solutions to both problems within reasonable computation time, and using only one simulation.

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

1611

6. Conclusion The vehicle routing problem with time windows is one of the classical research areas in operations research with considerable economic signi2cance. It has many real-life applications, especially in the transportation industry. For industrial problems, methods that are able to produce high-quality results in limited time, even for several hundreds of customers, are particularly important. In this paper we have focused on large-scale problems, and developed a novel two-phase metaheuristic, called active guided evolution strategies. In the 2rst phase an initial solution is generated with the hybrid cheapest insertion heuristic of Mester [11]. In the improvement phase the two powerful metaheuristics guided local search and evolution strategies are combined into a new iterative two-stage procedure. The guided local search is used to guide three standard improvement heuristics that are applied cyclically in the 2rst stage, and the local search operations of the modi2ed evolution strategies of Mester [11,35] in the second stage. The new method has been tested on 302 large-scale test problem instances from the literature. To illustrate the power and robustness of the suggested method, it was also tested on 23 capacitated vehicle routing problem benchmarks. The results show that the proposed metaheuristic is eNcient and very competitive in comparison to the previous heuristic solution methods, providing best-known solutions to 86% and 91% of the tested time window-constrained and classical vehicle routing problems, respectively. Acknowledgements This work was partially supported by the Ministry of Absorption of Israel, the Liikesivistysrahasto and Volvo Foundations, and the TOP project funded by the Research Council of Norway. This support is gratefully acknowledged. We thank also Wout Dullaert for his valuable comments. References [1] Koskosidis YA, Powell WB, Solomon MM. An optimization-based heuristic for vehicle routing and scheduling with soft time window constraints. Transportation Science 1992;26:69–85. [2] Taillard E, Badeau P, Gendreau M, Guertin F, Potvin J-Y. A tabu search heuristic for the vehicle routing problem with soft time windows. Transportation Science 1997;31:170–86. [3] Descrochers M, Desrosiers J, Solomon MM. A new optimization algorithm for the vehicle routing problem with time windows. Operations Research 1992;40:342–54. [4] Desrosiers J, Dumas Y, Solomon MM, Soumis F. Time constrained routing and scheduling. In: Ball MO, Magnanti TL, Monma CL, Nemhauser GL, editors. Handbooks in operations research and management science 8: network routing. Amsterdam: Elsevier; 1995. p. 35–139. [5] Cordeau J-F, Desaulniers G, Desrosiers J, Solomon MM, Soumis F. The VRP with time windows. In: Toth P, Vigo D, editors. The vehicle routing problem, SIAM monographs on discrete mathematics and applications. Philadelphia: SIAM; 2001. p. 157–94. [6] Br,aysy O, Gendreau M. Vehicle routing problem with time windows, part I: route construction and local search algorithms. Transportation Science 2002, in press. [7] Br,aysy O, Gendreau M. Vehicle routing problem with time windows, part II: metaheuristics. Transportation Science 2002, in press. [8] Solomon MM. Algorithms for the vehicle routing and scheduling problems with time window constraints. Operations Research 1987;35:254–65.

1612

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

[9] Potvin J-Y, Rousseau J-M. A parallel route building algorithm for the vehicle routing and scheduling problem with time windows. European Journal of Operational Research 1993;66:331–40. [10] Bramel J, Simchi-Levi D. Probabilistic analyses and practical algorithms for the vehicle routing problem with time windows. Operations Research 1996;44:501–9. [11] Mester D. A parallel dichotomy algorithm for vehicle routing problem with time windows. Working paper, Minerva Optimization Center, Technion, Israel, 1999. [12] Ioannou G, Kritikos M, Prastacos G. A greedy look-ahead heuristic for the vehicle routing problem with time windows. Journal of the Operational Research Society 2001;52:523–37. [13] Dullaert W, Br,aysy O. Routing with relatively few customers per route. Top 2002, in press. [14] Thompson PM, Psaraftis HN. Cyclic transfer algorithms for multivehicle routing and scheduling problems. Operations Research 1993;41:935–46. [15] Potvin J-Y, Rousseau J-M. An exchange heuristic for routeing problems with time windows. Journal of the Operational Research Society 1995;46:1433–46. [16] Russell RA. Hybrid heuristics for the vehicle routing problem with time windows. Transportation Science 1995;29: 156–66. [17] Shaw P. Using constraint programming and local search methods to solve vehicle routing problems. In: Maher M, Puget J-F, editors. Principles and practice of constraint programming—CP98, Lecture notes in computer science. New York: Springer; 1998. p. 417–31. [18] Caseau Y, Laburthe F. Heuristics for large constrained vehicle routing problems. Journal of Heuristics 1999;5: 281–303. [19] Schrimpf G, Schneider J, Stamm-Wilbrandt H, Dueck G. Record breaking optimization results using the ruin and recreate principle. Journal of Computational Physics 2000;159:139–71. [20] Cordone R, Wol4er Calvo R. A heuristic for the vehicle routing problem with time windows. Journal of Heuristics 2001;7:107–29. [21] Br,aysy O, Hasle G, Dullaert W. A multi-start local search algorithm for the vehicle routing problem with time windows. European Journal of Operational Research 2002, in press. [22] Glover F. Future paths for integer programming and links to arti2cial intelligence. Computers and Operations Research 1986;13:533–49. [23] Holland JH. Adaptation in natural and arti2cial systems. Ann Arbor: University of Michigan Press; 1975. [24] Rechenberg I. Evolutionsstrategie. Stuttgart: Fromman-Holzboog; 1973. [25] Rochat Y, Taillard E. Probabilistic diversi2cation and intensi2cation in local search for vehicle routing. Journal of Heuristics 1995;1:147–67. [26] Chiang WC, Russell RA. A reactive tabu search metaheuristic for the vehicle routing problem with time windows. INFORMS Journal on Computing 1997;9:417–30. [27] Schulze J, Fahle T. A parallel algorithm for the vehicle routing problem with time window constraints. Annals of Operations Research 1999;86:585–607. [28] Cordeau J-F, Laporte G, Mercier A. A uni2ed tabu search heuristic for vehicle routing problems with time windows. Journal of the Operational Research Society 2001;52:928–36. [29] Berger J, Barkaoui M, Br,aysy O. A route-directed hybrid genetic approach for the vehicle routing problem with time windows. INFOR 2003;41:179–94. [30] Jung S, Moon B-R. A hybrid genetic algorithm for the vehicle routing problem with time windows. In: Proceedings of Genetic and Evolutionary Computation Conference. San Francisco: Morgan Kaufmann, 2002. p. 1309–16. [31] Le Bouthillier A, Crainic TG. A cooperative parallel meta-heuristic for the vehicle routing problem with time windows. Working Paper, Centre for Research on Transportation, University of Montreal, Canada, 2003. [32] Homberger J, Gehring H. Two evolutionary meta-heuristics for the vehicle routing problem with time windows. INFOR 1999;37:297–318. [33] Gehring H, Homberger J. A parallel hybrid evolutionary metaheuristic for the vehicle routing problem with time windows. In: Miettinen K, M,akel,a M, Toivanen J. editors. Proceedings of EUROGEN99. Jyv,askyl,a; University of Jyv,askyl,a, 1999. p. 57– 64. [34] Gehring H, Homberger J. Parallelization of a two-phase metaheuristic for routing problems with time windows. Asia-Paci2c Journal of Operational Research 2001;18:35–47.

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

1613

[35] Mester D. An evolutionary strategies algorithm for large scale vehicle routing problem with capacitate and time windows restrictions. Working Paper, Institute of Evolution, University of Haifa, Israel, 2002. [36] Kontoravdis GA, Bard JF. A GRASP for the vehicle routing problem with time windows. INFORMS Journal on Computing 1995;7:10–23. [37] Liu F-H, Shen S-Y. A route-neighborhood-based metaheuristic for vehicle routing problem with time windows. European Journal of Operational Research 1999;118:485–504. [38] Voudouris C. Guided local search for combinatorial problems. Dissertation: University of Essex, UK; 1997. [39] Kilby P, Prosser P, Shaw P. Guided local search for the vehicle routing problem with time windows. In: Voss S, Martello S, Osman IH, Roucairol C, editors. Meta-heuristics advances and trends in local search paradigms for optimization. Boston: Kluwer; 1999. p. 473–86. [40] Gambardella LM, Taillard E, Agazzi G. MACS-VRPTW: a multiple ant colony system for vehicle routing problems with time windows. In: Corne D, Dorigo M, Glover F, editors. New ideas in optimization. London: McGraw-Hill; 1999. p. 63–76. [41] Dorigo M, Di Caro G, Gambardella LM. Ant algorithms for discrete optimization. Arti2cial Life 1999;5:137–72. [42] Br,aysy O. A reactive variable neighborhood search for the vehicle routing problem with time windows. INFORMS Journal on Computing 2003;15:347–68. [43] Rousseau L-M, Gendreau M, Pesant G. Using constraint-based operators to solve the vehicle routing problem with time windows. Journal of Heuristics 2002;8:43–58. [44] Mladenovic N, Hansen P. Variable neighborhood search. Computers & Operations Research 1997;24:1097–100. [45] Kirkpatrick S, Gelatt CD, Vecchi PM. Optimization by simulated annealing. Science 1983;220:671–80. [46] Chiang WC, Russell RA. Simulated annealing metaheuristics for the vehicle routing problem with time windows. Annals of Operations Research 1996;63:3–27. [47] Li H, Lim A. Large scale time-constrained vehicle routing problems: a general metaheuristic framework with extensive experimental results. AI Review 2001, in press. [48] Bent R, Van Hentenryck P. A two-stage hybrid local search for the vehicle routing problem with time windows. Transportation Science 2002, in press. [49] Thangiah SR, Osman I, Sun T. Hybrid genetic algorithm, simulated annealing and tabu search methods for vehicle routing problems with time windows. Working Paper, Institute of Mathematics and Statistics, University of Kent, UK, 1994. [50] Ibaraki T, Imahori S, Kubo M, Masuda T, Uno T, Yagiura M. E/ective local search algorithms for routing and scheduling problems with general time window constraints. Transportation Science 2002, in press. [51] Hasle G. Designing VRP software. Invited Talk at the Spring School on Logistics and Distribution Management, University of Montreal, Canada, 2002. [52] Schwefel H-P. Numerische optimierung von computer-modellen mittels der evolutionsstrategie. Basel: Birkh,auser; 1977. [53] Fogel L, Owens A, Walsh M. Arti2cial intelligence through simulated evolution. Chichester: Wiley; 1966. [54] Nissen V. Evolutionare Algorithmen. Wiesbaden: Deutscher Universitats-Verlag; 1994. [55] Fogel DB. Evolutionary computation: toward a new philosophy of machine intelligence. New York: IEEE Press; 1995. [56] B,ack T, Fogel DB, Michalewicz Z. Handbook of evolutionary computation. Oxford: Institute of Physics Publishing and Oxford University Press; 1997. [57] Voudouris C, Tsang E. Guided local search. European Journal of Operations Research 1998;113:80–119. [58] Voudouris C, Tsang E. Partial constraint satisfaction problems and guided local search. In: Proceedings of practical application of constraint technology. London: Practical Application Company; 1996. p. 337–56. [59] Voudouris C, Tsang E. Guided local search and its application to the traveling salesman problem. European Journal of Operational Research 1999;113:469–99. [60] Tsang E, Voudouris C. Fast local search and guided local search and their application to British telecom’s workforce scheduling problem. Operations Research Letters 1997;20:119–27. [61] Savelsbergh MWP. The vehicle routing problem with time windows: minimizing route duration. INFORMS Journal on Computing 1992;4:146–54. [62] Osman IH. Metastrategy simulated annealing and tabu search algorithms for the vehicle routing problems. Annals of Operations Research 1993;41:421–52.

1614

D. Mester, O. Br3aysy / Computers & Operations Research 32 (2005) 1593 – 1614

[63] Clarke G, Wright JW. Scheduling of vehicles from a central depot to a number of delivery points. Operations Research 1964;12:568–81. [64] Christo2des N, Mingozzi A, Toth P. The vehicle routing problem. In: Christo2des N, Mingozzi A, Toth P, Sandi C, editors. Combinatorial optimization. Chichester: Wiley; 1979. p. 315–38. [65] Fisher ML. Optimal solution of vehicle routing problems using minimum K-trees. Operations Research 1994;42: 626–42. [66] Taillard E. Parallel iterative search methods for vehicle routing problems. Networks 1993;23:661–73. [67] Mester D, Korol A, Nevo E. Guided evolution strategy algorithms for optimization of the several large scale genetic problems. Working Paper, Institute of Evolution, University of Haifa, Israel, 2003. [68] Reimann M, Doerner K, Hartl RF. D-ants: savings based ants divide and conquer the vehicle routing problem. Computers & Operations Research 2003, in press.