a network flow-based tabu search heuristic for the ... - Semantic Scholar

A NETWORK FLOW-BASED TABU SEARCH HEURISTIC FOR THE VEHICLE ROUTING PROBLEM Jiefeng Xu

Graduate School of Business, University of Colorado at Boulder, Boulder, CO 80309-0419 Email: [email protected]

and James P. Kelly

Graduate School of Business, University of Colorado at Boulder, Boulder, CO 80309-0419 Email: [email protected]

Last Revision: February, 1996.

Abstract. We develop a new local search approach based on a network

ow model that is used to simultaneously evaluate several customer ejection and insertion moves. We use this approach and a direct customer swap procedure to solve the well-known Vehicle Routing Problem. The capacity constraints are relaxed using penalty terms whose parameter values are adjusted according to time and search feedback. Tabu Search is incorporated into the procedure to overcome local optimality. More advanced issues such as intensi cation and diversi cation strategies are developed to provide eective enhancements to the basic tabu search algorithm. Computational experience on standard test problems is discussed and comparisons with best-known solutions are provided. Key words. Vehicle Routing Problem, Network Flow Model, Local Search Heuristic, Tabu Search.

2

Jiefeng Xu and James P. Kelly

Introduction The Vehicle Routing Problem (VRP) is an important management problem in the eld of physical distribution and logistics. A typical VRP can be described as the problem of designing least cost routes from one depot to a set of geographically scattered points, that can be cities, stores, warehouses, schools, customers, etc. The route must be designed in such a way that each point is visited once and only once by exactly one vehicle, all routes start and end at the depot, and the total demands of all points on one particular route cannot exceed the capacity of the vehicle that is assigned to execute that route. The VRP may contain other real-world constraints, such as maximum route length, time windows, maximum number of points visited in one route, etc. These side constraints complicate an already dicult problem. In practice, the VRP has been recognized as one of the great success stories of operations research. It is not dicult to nd VRP implementations resulting in millions of dollar savings (see Fisher et al 1982, Bell et al 1983, Brown et al 1987, etc.) . The economic importance of VRP applications cannot be overestimated. Given the fact that in 1991 approximately $506 billion was invested in the transportation and public utilities sector in the United States (Economic Report of the President, 1994), even a minor improvement in route design methodology may result in signi cant savings. The VRP has received extensive attention from operations researchers in the past three decades. Most of the research was reviewed in the survey by

A Network Flow-Based Tabu Search Heuristic for the VRP

3

Bodin, et al. (1983), the volume edited by Golden and Assad (1988), and an article by Laporte (1992). The VRP is a well-known NP-hard problem. Thus, the size of solvable instances of the VRP is limited. For example, Christo des (1985) reports solving a symmetric problem with 53 customers using a k-degree center tree algorithm, and Laporte, Mercure and Nobert (1986) have solved a randomly generated asymmetrical capacity constrained VRP with 260 customers using an assignment lower bound in a branch and bound procedure. More recently, Fisher (1994) has solved a 100 customer problem in Christo des, Mingozzi and Toth (1979) and another ve real problems containing 25 to 71 customers using a branch and bound algorithm for the associated mininum Ktree problems. Since most practical problems are large, heuristics that attempt to quickly nd good feasible solutions have become the focus of VRP research. Furthermore, heuristics can handle most practical constraints in a more exible way. Therefore, much of the current research eorts have been devoted to developing and re ning VRP heuristics. Most early heuristics for the VRP concentrated on constructing routes based on savings criteria. Among them are Clarke and Wright (1964), Gillett and Miller (1974), and Mole and Jameson (1976). These heuristics work fast, but the results could be easily improved by local search procedures such as 2-opt, 3-opt (Lin, 1965, Lin and Kernighan, 1973) and Or-opt (Or, 1976). Later heuristics, that employed mathematical programming techniques to solve relaxations of the VRP (For example, Fisher and Jaikumar 1981, Foster and

4


Ryan 1976, Cullen, Jarvis and Ratli 1981, and Stewart and Golden 1984, etc.), usually produced better solutions than the early heuristics. Recently, Harche and Raghavan (1993) employed a node swapping local search heuristic iterated over varying capacity constraints to produce good solutions for some problems. Noon, Mittenthal and Pillai (1994) used a Traveling Salesman subset tour relaxation to achieve good results. One of the milestones in heuristic design for the VRP was the introduction of Tabu Search (TS). TS is a metaheuristic for combinatorial optimization that was pioneered by Glover (1986). It explores the neighborhood of a solution and employs elegant recency-based memory structures to avoid being trapped in local optima. TS has received a growing awareness in the operations research community since its debut in 1986. There are numerous reports on successful applications of TS to VRP problems in the literature (i.e., Willard 1989, Pureza and Franca 1991, Osman 1993, Taillard 1993, Gendreau, Hertz, and Laporte 1994, Rochat and Taillard, 1995). Among them, several most recent papers, for example, Rochat and Taillard (1995), Taillard (1993) and Gendreau, Hertz, and Laporte (1994) use the advanced diversi cation strategy embedded within a Tabu Search structure to outperform the previous heuristics. In a recent working paper, Stewart, Kelly, and Laguna (1994) proposed a new algorithm that produced solutions that are better or as good as but more robust than that of Gendreau, Hertz, and Laporte using integer distances. They adapted Fisher and Jaikumar's formulation for VRP (Fisher and


5

Jaikumar 1981), modi ed their seed selection, and employed a more powerful heuristic that incorporates Tabu Search and ejection chains to quickly solve the generalized assignment problem to near-optimality. A signi cant contribution of their paper is the development of a network ow model to implement a diversi cation procedure. Their network ow model is further strengthened by TS through both short-term and long-term memory structures. This paper advocates the same network ow model used in Stewart, Kelly and Laguna (1994) as a general local search strategy to solve the VRP. Instead of employing the network ow model as a machanism to produce diversi ed moves and relying on solving the associated generalized assignment problems to obtain feasiblity as in Stewart, Kelly and Laguna (1994), our use of the network ow model is more straightforward. We relax the hard side constraints by introducing a dynamic penalty system, and eciently update and frequently solve the network ow model to nd the best customers to insert into new routes without the use of the generalized assignment algorithm. We present a new mechanism that couples the network ow model with the direct customer swap procedure. TS restrictions are applied to the operations of dropping a customer from its current route, inserting a customer into new route, and swapping two customers between dierent routes. The penalty parameters are changed such that the feasibility of the search is controlled. Computational experience is discussed and the results obtained from a set of well-known benchmark problems are presented.

6


This paper is organized as follows. The network ow model and corresponding local search heuristic are described in Section 1. Computational experience from our heuristic and the comparisons with other heuristics are provided in Section 2. We summarize our ndings in Section 3 and report the new bestknown solutions found by our algorithm in the appendix.

1. Network Flow-Based TS Algorithm In this section, we describe the network ow model, the direct swap procedure, TSP components, TS restrictions, diversi cation strategies, and the dynamic search procedure. We illustrate this methodology using the classical capacityconstrained VRP problem, though it can be readily extended to solve other variants of the VRP. We now describe the capacity-constrained VRP. Let v1 ; v2; : : : ; vn be n customers each with demand q1; q2 ; : : : ; qn respectively, and v0 be the depot where m identical vehicles with capacity Q are based. Denote R1 ; R2; : : : ; Rm as the vehicle routes for vehicle 1; 2; : : : ; m, each contains a speci c permutation of a disjoint subset of customers and the depot. The distance dij is de ned over the arc (vi; vj ) for any i; j 2 f0; 1; : : : ; ng. The object is to minimize the total distances on all routes such that the sum of customers' demands on each route does not exceed Q.

7


1.1. Network Flow Model. This model is an extension of the ejection chains model for Traveling Salesman Problem (TSP) that appeared in Glover (1992). Stewart, Kelly and Laguna (1994) also use it as a part of a diversi cation and improvement mechanism. Its purpose here is to provide a mechanism for exchanging customers between routes such that distance and feasibility are simultaneously considered. The network model is shown in Figure 1. Customer Nodes

Vehicle Nodes

1 Vehicle Nodes

1

1 2 2

2 C

3

. . .

S

T

C

. . . n-1

m m Level 1

n

Level 2 Current Assignments All arcs are bounded by [0,1]

Level 4

Level 3 Potential Assignments All arcs are bounded by [0,1]

Figure 1: Network Flow Model As shown in Figure 1, this model is composed of four distinct levels of arcs. On the leftmost level (Level 1), a source node S with the input ow of c units, is connected with m nodes, that represent the m vehicle routes.

8


The ows on these arcs express the number of customers that are removed from current routes, while the input ow c speci es the total number of such moves. On the next level (Level 2), all m vehicle nodes are connected with n nodes that represent the n customers. The ows on these arcs imply that the corresponding customer is removed from its current route. Similarly, the ows along arcs on Level 3 indicate customers that are inserted into a particular route other than the current one. On Level 4, the ows represent the number of customers adding to the routes. Finally, the model terminates at sink T with the output ow of c units. The ows on Levels 2 and 3 are bounded by [0,1] and any nonzero ow indicates a move. The foregoing network ow model may be further simpli ed by connecting the source node S directly to each of the customer nodes. However, the current model is more generic and exible. For example, we can impose a bound on arcs of Level 1 to specify the maximum or minimum number of customers which are allowed to be removed from that vehicles route. We adhere to the complete model though the above \advanced" features are not implemented in this paper. Let fijl be the ow along the arc from node i to node j on level l, Cijl be the cost on that arc, Llij and Uijl be the lower bound and upper bound on that arc respectively, then the mathematical representation of this model is given


below: minimize subject to:

nX m X j =1i=1

Cij2 fij2

+

m X i=1

fsi1 ?

n X j =1

n X m X

9

Cjk3 fjk3

j =1k=1

fsi1 = c

fij2 = 0

m m X X fij2 ? fjk3 = 0 i=1 n k=1 X 3 4 fjk ? fkt = 0 j =1

(1) (2)

i = 1; : : : ; m;

(3)

j = 1; : : : ; n

(4)

k = 1; : : : ; m

(5)

m X fit4 = c i=1

(6)

L1si fsi1 Usi1 i = 1; : : : ; m

(7)

L2ij fij2 Uij2 i = 1; : : : ; m j = 1; : : : ; n

(8)

L3jk fjk3 Ujk3 k = 1; : : : ; m j = 1; : : : ; n

(9)

L4kt fkt4 Ukt4 k = 1; : : : ; m There are several interesting features of this model. Given that c and all the lower and upper bounds are integer, then the optimal ows must also be integer. The ows fij2 and fji3 are bounded by [0,1], and represent decision variables on the respective deletion and insertion operations. From a ow balance point of view, if a customer is to be removed from its current route, then it must also be assigned to another route. If c is xed to one, the model actually

(10)

10


represents the best ejection/insertion move. For c larger than one, the model activates a composite move, that is a combination of multiple ejection/insertion moves. Non-permissible ejections/insertions can be modeled by either assigning a prohibitively large cost to the corresponding arcs or by changing the appropriate upper bounds to zero. For example, if customer j is currently in route i, then all arcs from route nodes other than i to customer node j in Level 2 are \infeasible", as well as the arc from customer node j to route node

i in Level 3. Thus, this model can represent any customer-route con guration by appropriately changing the arc costs and/or upper bounds. Additionally, the model can be updated eciently to accommodate the local search moves. Since subsequent changes of customer-route con gurations only alter the cost coecients and upper bounds, the standard primal simplex and dual simplex methods can be applied directly to obtain the new optimal ows. One improvement in the eciency of the heuristic could be obtained through the use of a network solver in place of the Linear Programming (LP) solver. The authors did not have access to a good network solver and chose to use the well-known CPLEX LP solver instead of a generic network solver. The cost terms Csi1 and Ckt4 are set to zero. All costs on \infeasible" moves are set to BIGM , a prohibitively large number and associated upper bounds are set to zero. The other cost terms Cij2 and Cjk3 for i; k = 1; : : : ; m and

j = 1; : : : ; n incorporate the net distance change by the respective move and a penalty term used to motivate satisfaction of the capacity constraints. More


11

speci cally, let d2ij denote the net distance change by removing customer j from its route i, and d3jk be the net distance change by inserting customer j to route k, then

Cij2 = d2ij + '1(

X

ql ) (ejection cost): l2R X Cjk3 = d3jk + '2( ql ) (insertion cost): i

l2Rk

where '1 (x) and '2(x) are piecewise functions such that 8 >< p1(H1 ? x) if x > H1 '1 (x) = > p1(L1 ? x) if x < L1 :0 otherwise 8 >< p2(x ? Q ? q) if x > Q + q if Q < x Q + q '2 (x) = > p3(x ? Q) :0 otherwise and p1 , p2 , p3 are positive penalty parameters, while H1 and L1 (H1 > L1 ) are high and low threshold parameters for the total demand on each route respectively, and q is the average demand for all customer (q = Pnj=1 qj =n).

d2ij can be calculated by removing customer j from route i. For d3jk , we need to nd the least cost insertion position in route k for customer j and the associated net distance change. If only a single customer is removed or inserted, then the associated net distance changes are exact. If several deletions/insertions occur in a single route simultaneously, then these calculations may only approximate the actual move costs. The above cost calculations have an eective encouragement or discouragement mechanism realized by the penalty terms. For example, in the calculation of Cij2 , if Pl2R ql > H1, then the capacity of route i is \over-saturated" i

12


or \nearly-saturated" and will encourage the deletion of customer j from its

current route i since the penalty term '1(Pl2R ql ) is negative. If Pl2R ql < L1 , i

i

the vehicle is under-loaded and the calculation will discourage the deletion of customer j from route i by adding a positive penalty. Similarly, for Cjk3 , when customer j is inserted into route k when the vehicle is highly or moderately

over-loaded, then '2(Pl2R ql ) is positive, and discourages the insertion through k

the cost term. The system provides a mechanism that considers the tradeo between distance minimization and feasibility. The feasibility can be controlled if the values of the penalty parameters are carefully selected and altered.

1.2. Direct Customer Swap Procedure. The direct customer swap is a commonly used local improvement approach. A \swap" simply exchanges two customers between two routes. The best swap is de ned as the move that most reduces the total distance. Extending this concept, the cost associated with a swap may include the net distance change and a penalty term. Let Cij be the cost of swapping customer i and j , and dij be the net distance change caused by the swap, k and l be the corresponding indices of the routes in which customers i and j currently reside, thereby

X X Cij = dij + p4(max(0; qr ? qi + qj ? Q) + max(0; qr ? qj + qi ? Q)): r2Rk

r2Rl

where the p4 is a positive penalty parameter that activates the encouragement/discouragement mechanism.


13

In many cases, dij can be derived by taking advantage of the pre-computed d2ik , d2jl , d3il , d3jk . More speci cally, the net distance change associated with swap of customer i and j , is the sum of the net distance changes resulted from removing customer i and j from route k and l respectively, and adding customer i to route l and customer j to route k. That is, dij = d2ik + d2jl + d3il + d3jk : However, if the least cost insertion position of the inserting customer is next to the other customer that is to be dropped from this route, then the above calculation is no longer accurate. In this case, we can calculate the net distance change directly by removing the two customers from their current routes and then inserting each into their respective routes. Of course, all of these calculations are approximate since the sequence of other customers in the route remains unchanged. Note that it is possible to swap two customers within the same route. We disable such moves since these moves are considered in the TSP component that follows in Section 1.3. As with the move costs described above, the swap cost Cij has a term that penalizes \over-loaded" routes. Compared with the ejection moves produced by the network ow model, it has a more \balanced" performance. If the penalty parameter p4 is large enough, and the current solution is feasible or \nearly" feasible, then the swap with the minimum cost is likely to produce a solution that is still feasible or \nearly feasible". For those solutions that are far from

14


feasibility, the swap will drive the search towards feasibility. In contrast, if the penalty term is small, then the swap may lead to a more distance-based customer clustering. One of the ndings of this work is that oscillating between swap and ejection/insertion moves is highly eective.

1.3. TSP Component. Most VRP heuristics use TSP components to nd the best sequence of customers that are clustered into a route. In our case, we already consider optimizing the move costs (although they are approximate for some instances), so the TSP is merely employed as an improvement method. Two types of TSP algorithms, namely the 3-opt and a stand-alone TS heuristic for TSP, are applied in our algorithm. The 3-opt is a well-known local search improvement method for the TSP. Its moves delete three edges from a current existing route and reconnect them in dierent ways, attempting to nd a better route. Our 3-opt is preceded with a farthest insertion heuristic and a 2-opt procedure. Empirical experience shows that such a component works well for small TSP, however for larger TSP, the gap between the heuristic solution and the optimal route may be considerable. We developed a simple TS component for the TSP (TSTSP) to reduce this gap. In the TSTSP, two dierent moves are considered: ejection and swap. The former move will eject one customer to another position while the later will exchange two customers' positions. The cost of these moves are de ned by the net distance change. The TSTSP proceeds with the initial input sequence.


15

In each iteration, it evaluates all eligible candidate moves and selects the best move (not necessarily an improving one). Once a move is made, an opposite move that may lead the search back to the previous solution is tabu for a certain number of iterations (tabu tenure). However, the tabu status of a move may be overridden if such a move would lead to a new best solution (aspiration criterion).

1.4. Tabu Search Memory. The network ow-based method and swap procedure are local search approaches that may be easily trapped in local optima if no TS restrictions are imposed. TS restrictions with randomly generated tabu tenures are applied to three dierent neighborhood moves: dropping a customer from its current route, inserting a customer into a dierent route, and swapping two customers between routes. For a swap, in addition to tabu restrictions on future swaps, the associated ejections and insertions are also subject to tabu restrictions. Additionally, when a customer is moved to a new route, a tabu restriction that prevents its removal from that route is only activated when there are only a few customers (less than a pre-determined number) in the route. When a move is marked as tabu, the associated costs are changed to BIGM so that such a move will be forbidden. When the tabu restriction of a move is released, the recalculation of its associated cost (changed from BIGM to its approximated move cost) is required. A recency-based memory structure is

16


used to trace the tabu status for the tabu tenure, that is, when a move is to be tabu, we set

tabu(move) = iter + tabu tenure: where iter is the current iteration number. The tabu status can be determined by if iter tabu(move); then the move is tabu; else the move is permissible: In this implementation, the tabu status is not overridden (aspiration criterion) until the search is restarted as a part of the diversi cation strategy. The local search will terminate at a pre-de ned maximum number of iterations (iter max).

1.5. Intensi cation and Diversi cation Strategies. Intensi cation and diversi cation may often improve the power of a metaheuristic search method. In this research, we develop and implement a diversi cation strategy based on long term frequency memory and a intensi cation strategy based on advanced restart/recovery procedure. (1) Frequency memory implementation: A frequency memory-based penalty is added to the move costs to break ties. This memory records the frequency of the customer being assigned to a speci c route. The frequency is normalized by dividing it by the maximum frequency and then is added to the move cost. It


17

breaks ties by giving preference to those customers that appear less frequently in a speci c route. This strategy is easily implemented. However, since routes cannot be uniquely identi ed over the entire search, it is only used as a minor diversi cation strategy. (2) Advanced restart/recovery strategy: The set of best feasible solutions produced by the search are de ned as elite solutions. A repository of elite solutions is maintained. Advanced restart is executed periodically during the late stages of the search. When restarting, the current solution is obtained from the repository and all tabu restrictions are released and all penalty parameters are initialized. This strategy is based on the assumption that there may exist short relinking paths in the search process from the restart points to new local or global optima. However, these paths may not be detected during the prior search due to the tabu restrictions and inappropriate penalty parameter values (these parameters are dynamically changed as described later). The advanced restart/recovery strategy may nd these paths and thereby lead the search to new local or global optima. It is crucial to select the right restart points. We keep the repository as a priority queue. The queue has the key value de ned by the objective value of the elite solution. Once a new solution is found better than the worst solution in the repository, the worst one is replaced by this solution. Furthermore, this repository is kept sorted so the worst solution can be immediately identi ed. For the purposes of better diversifying and sampling,

18


the restart points in the pool should be somewhat distinct from each other. If a solution has the same key value as one in the repository, we compare their customer-route con guration to accept only the suciently distinct solution (i.e., with more than two dierent routes). Thus, the restart points are more likely to be independent and the search process is more eective.

1.6. Dynamic Search Procedure. The search procedure starts from any feasible or infeasible customer-route con guration. During the dynamic search, penalty parameters and search parameters are changed based on various scenarios which in turn depend on the current solution. Generally, if the capacity violation is high, the penalty parameters should be set to drive the search back to feasibility. If a solution is found close (may be infeasible with low capacity violation) to the current best solution, swaps move may be useful for producing feasible solutions with high quality. Our search oscillates between network ow moves and direct swap moves. While network ow moves dominate the search, swap moves are executed once every s1 iterations, or are triggered by some speci c conditions described later. Let cur sol and best sol be the objective values of the current solution and best feasible solution in the search respectively, and be the sum of the capacity constraint violations for the current solution, that is =

m X i=1

maxf0;

X

j 2Ri

qj ?Qg.


19

Then the current solution is feasible if and only if = 0. To illustrate the oscillation between the network ow model and direct swaps, we de ne iter swap such that the search will perform swaps if and only if iter swap iter. Otherwise, the network ow model is executed. In addition, de ne iter net as the most recent iteration in which the network ow model was used, and best iter as the iteration counter at which the current best sol is obtained. Finally, let

1 , 2 and be pre-de ned tolerances of the solution values , and 1 through 3 be pre-de ned measures of the capacity violation, n1 through n9 represent pre-determined iteration counters, and 4 ; 4; _4 and 4 be the possible values for p4 (penalty parameter associated with a swap move). The search procedure is de ned as follows:

Pseudo-Code Listing 1 Search Procedure p4 = 4 : iter = 0: iter swap = 0: While (iter < max iter) do If (iter mod s1 = 0) or (iter swap > iter) do

Else

Evaluate the swap moves, update the elite solution list if necessary, and select the best move as cur sol;

Evaluate the network ow model, update the elite solution list if necessary, and select the best move as cur sol; Improve the current solution using 3-opt or TSTSP if necessary (see TSP components).

20


If (cur sol ? best sol 1) and ( 1 ) do If If If If

iter swap = iter + n1 ; p4 = _4 : (1 < 2) do p4 = 4 : ( > 2 ) do p4 = 4 : ( = 0) and (cur sol ? best sol ) do p4 = 4 ; iter swap = iter swap + n2 ; (iter > iter net + n3 ) do iter swap = iter ? 1:

Change parameters if necessary (see dynamic parameter adjustment pseudo-code). Apply the TS restriction by updating the move costs and upper bounds. If (cur sol < best sol) do best sol = cur sol; best iter = iter:

Recover the current solution by a selected elite solution if necessary.

iter = iter + 1: The search starts at iter = 0 and p4 is initialized to a small number 4 . At each iteration, a best move is determined by the network ow model or the swap procedure, and a new solution is found. If the search nds a solution that is within 1 of the best, and the capacity violation is limited ( < 1 ), running the direct swap procedure for n1 iterations is then mandated, and its penalty parameter p4 is set to an average value _4. If 1 < 2 (indicates the capacity violation is sucient but not extremely high), p4 is again reduced to a small number 4. The search focuses on the distance minimization in this case. If the capacity violation is extremely high ( > 2 ), p4 is set to a large number 4 that will drive the search towards feasibility. If a feasible solution


21

is found within of the best (surpass or very close to the best), then n2 more swap executions are added to enhance the search, and p4 is set to an higher value 4. Finally, after n3 consecutive direct swap executions, the search will switch to the network ow model. We now describe the mechanism that dynamically changes the penalty parameters for the network ow model. This leads the search to distance minimization with controlled feasibility. Let iter penalty represent the iteration counter such that the penalty parameters p1 through p3 are changed if and only if iter > iter penalty. De ne i as initial setting for pi, and r1 and r2 as multipliers of i where i = 1; 2; 3 . The pseudo-code of the dynamic parameters adjustment is de ned as follows.

Pseudo-Code Listing 2 Dynamic Parameter Adjustment Initialize pi = i for i = 1; 2; 3: Initialize iter penalty = n4: At each iteration, do If ( > 2 ) do iter penalty = iter ? 1: If (iter > iter penalty) do If ( 3 ) do

iter penalty = iter + n5 ; pi = 0; If iter > best iter + n6 do iter penalty = iter penalty + n7 : If (3 < 2 ) do iter penalty = iter + n8 ; pi = r1i for i = 1; 2; 3: If ( 2 ) iter penalty = iter + n9 ; pi = r2i for i = 1; 2; 3:

22


In the dynamic parameters adjustment procedure, the parameters are set to their initial value settings for n4 iterations. In each iteration, their values are adjusted based on various scenarios of the capacity violation. If the capacity violation is extremely high ( > 2 ), then an immediate adjustment of parameters is mandated. If the capacity violation is not high ( 3 ), then all parameters are set to zero. Since this setting represents the most diversi ed case, it can be further reinforced by n7 iterations when the search cannot improve the solution for the most recent n6 iterations. If is moderate (3 < 2), then the parameters (pi) are assigned to the average values r1 i for n8 iterations. If is extremely high ( 2 ), the parameters are set to r2 i for n5 iterations. In addition to the penalty parameters, the total number of network moves (value of c) is changed dynamically. Its value is reset periodically and can be determined by moving the pointer in a pre-de ned circular list. The larger the value of c, the more diverse the move is. However, extremely large values for c deteriorate the solution quality since cost estimations become less accurate as c is made larger. We now brie y describe the TSP components of the search. First, the 3opt are executed every n9 iterations or when a feasible solution that is close to the best solution available is found (cur sol ? best sol 2 ). After 3opt is completed, the TSTSP is executed only when cur sol ? best sol . Repetitive executions of 3-opt and TSTSP can be avoided by using a checkthen-run strategy. To perform this, we label the routes that have been changed since the last execution of 3-opt and TSTSP. TSP heuristics are only performed on labeled routes.


23

2. Computational Experience The computational tests were conducted on a set of benchmark problems with capacity constraints. These problems are taken from the literature and their characteristics are listed in Table 1. Brie y, problems 1 through 7 are from Christo des, Mingozzi and Toth (1979). In problems 1 through 5, customers are randomly generated in the plane, while in problems 6 and 7, they are in clusters. Problems 2a through 2d and 3b are problems 2 and 3 using dierent vehicle capacities and numbers of vehicles as suggested by Russell (1977). The last three problems (f44, f71 and f134) having 44, 71 and 134 customers respectively, are taken from Fisher (1994). All fourteen problems have very

tight capacity constraints, as the ratios of demand to capacity (Pni=1 qi=(mQ)) in Table 1 are close to 1.0. The last column of Table 1 displays the standard deviation of customer demands for each problem, indicating that the customer demands in Fisher's problems vary much more than in the other benchmark problems.

24


Test Customer Vehicle Vehicle Ratio of Demand Standard Deviation Problem Number Number Capacity to Capacity of Customer Demand 1 50 5 160 0.97 8.06 2 75 10 140 0.97 7.96 3 100 8 200 0.91 8.87 4 150 12 200 0.93 8.59 5 199 17 200 0.94 8.51 6 100 10 200 0.91 10.41 7 120 7 200 0.98 5.38 2b 75 14 100 0.97 7.96 2c 75 8 180 0.95 7.96 2d 75 7 220 0.89 7.96 3b 100 14 112 0.93 8.87 f44 44 4 2010 0.90 258.55 f71 71 4 30000 0.96 3009.05 f134 134 7 2210 0.95 187.32

Table 1: Benchmark Problems Computational results that have appeared in the literature are often varied since rounded, truncated, or real Euclidean distances are used. Generally, most earlier studies relied on integer distance while more recent authors incline to report their nding using real distances (normally in single precision). Nevertheless, there are still some new integer solution values reported recently (for example, Stewart, Kelly and Laguna, 1994). There are no widely accepted rules on distance measures. On the other hand, we nd that the performance dierence due to the dierent distance calculations are quite signi cant, especially for those heuristics attempting to produce extremely high quality solutions, i.e., new best solutions. Therefore, to make direct comparisons with the most recent results, we report solutions using both real and rounded distances. Our algorithm can start with any feasible or infeasible customer-route con guration. We use a sweep method to produce the initial con guration. Starting with any customer, it sweeps customers into a route by their polar angles


25

until the capacity is saturated, and then continues with a new route. For problems in which the number of available vehicles is given, when the maximum number of vehicles are used, the remaining unclustered customers are assigned to neighboring routes. After the clusters are formed, we compute the center of gravity for each cluster (representing a route), and then assign each customer to its nearest route (in terms of the distance between the customer and the center of gravity of that route) in which the sum of demands (not including that of the inserting customer) on that route does not exceed the capacity. This simple heuristic often produces inferior infeasible solutions. However, our algorithm is insensitive to the initial feasibility and solution quality. The network ow model is solved using CPLEX 3.0. Once the initial LP model is solved, we only need to change the cost coecients and upper bounds for the subsequent moves. Then, the revised model is solved by the dual simplex method. As mentioned previously, replacing the LP solver with a specialized network solver may increase the eciency of the algorithm.

2.1. Parameter Settings. We use a dynamic tabu tenure that is generated from a uniform distribution from an interval [a; b]. We use the interval [3,9] for dierent moves. After inserting a customer into a route, the customer is not

tabu unless there are no more than two customers in the route. The search terminates when iter reaches iter max = 100; 000. The advanced restart/recovery strategy begins at iteration 20,000, and for every 750 iterations thereafter. Each recovered solution is selected by moving a pointer in the repository that stores

26


37 sorted elite solutions. The pointer is moved circularly in the repository, and starts from the worst solution. The penalty parameters p1 through p4 are initialized to 1 through 4 , and may be reset to dierent settings on speci c conditions. Their values are listed in Table 2. Parameter Value Parameter Value 1 0.2 4 1 2 0.2 4 8 3 1 _ 4 + d iter= 5000e 4 r1 5 4 6 + diter=5000e r2 10 Table 2: Penalty Parameters Settings Time related parameters and the control parameters regarding tolerances of solution, gap, capacity violations, and capacity thresholds are listed in Table 3. Parameter Value Parameter Value n1 35 s1 4 n2 5 1 15 n3 100 2 minf2Q + 20diter=500e; 400g n4 20 3 100+15diter=5000e n5 5+2diter=5000e 0.1 n6 1010 1 20 n7 5 2 10 n8 150+10diter=5000e L1 0:8Q n9 15 H1 0:95Q Table 3: Additional Parameters Settings


27

Finally, the number of network moves c is changed every 210 iterations. Its value is read from a circular list f1, 2, 3, 1, 2, 3, 1, 2, 3, 4g. Clearly, the algorithm employs dynamic penalty functions that vary with dierent scenarios based on the feasibility and solution quality. Therefore many parameters are required. However, most of them are relatively insensitive to their values and were set via common sense. A few \sensitive" parameters are chosen through a handful of experiments on several problem instances (problem 1 through 7 using real distances). Using systematical ne-tuning procedure based on statistical tests (for example, Xu, Chiu and Glover, 1996), one may nd more appropriate parameter value settings and further improve the performance of the procedure.

2.2. Comparisons. We conducted all testing on a DEC ALPHA workstation (DEC OSF/1 v3.0) and reported computation times in minutes. We rst present the performance of our algorithm on seven problems (problem 1 through 7) using real distances in Table 4. We list times required by our algorithm, times required to obtain our nal solution, the best solutions published, times required to obtain the solutions within 5% and 1% of the best. The best solutions published are obtained from Taillard (1993) and Rochat and Taillard (1995). Out of the seven best known solutions, the solution values for problem 1 (524.61) and problem 6 (819.56) were proven to be optimal by Hadjiconstantinou, Christo des and Mingozzi (1995) and Fisher (1994) respectively.

28

Jiefeng Xu and James P. Kelly Test Problem 1 2 3 4 5 6 7

Time to Best Time to Time to Our Total Obtain Our Solution Obtain Obtain Solution Time Solution Published 5% of Best 1% of Best 524.61 29.92 4.89 524.61 0.65 0.77 835.26 48.8 37.94 835.26 0.05 6.72 826.14 71.93 47.45 826.14 0.15 2.38 1029.56 149.9 100.81 1028.42 4.09 38.41 1298.58 272.52 207.80 1291.45 6.22 201.18 819.56 56.61 5.61 819.56 0.09 2.39 1042.11 91.23 2.80 1042.11 0.52 1.18

Table 4: Computational Results on Seven problems with Real Distances The results in Table 4 show that our algorithm can nd high quality solutions in reasonable times. Out of the seven problems, we nd the best known solution for ve problems with up to 120 customers. For the two large problems with 150 and 199 customers (problem 4 and 5), the gaps between our solutions and the best known solutions are small (0.11% for problem 4 and 0.55% for problem 5). In addition, our algorithm can produce solutions within 1% of the best in relatively short times except in problem 5. We compare our algorithm with several recent results on real distance from literature, namely, Osman (1993), Taillard (1993) and Gendreau, Hertz and Laporte (1994), as shown in Table 5. Problem 1 2 3 4 5 6 7

Osman 524.61 844.00 835.00 1044.35 1334.55 819.59 1042.11

Taillard Gendreau, Hertz and Laporte Our Solution 524.61 524.61 524.61 835.26 835.32 835.26 826.14 826.14 826.14 1028.42 1031.07 1029.56 1298.79 1311.35 1298.58 819.56 819.56 819.56 1042.11 1042.11 1042.11

Table 5: Comparison with Other Methods on Real Distances


29

Table 5 clearly shows that our algorithm is very competitive with the other recent TS heuristics. Compared with that of Osman (1993), our algorithm signi cantly improves in four cases and ties in the remaining three problems. For solutions produced by Taillard (1993), we obtain one slightly better solution and one slightly worse solution, and tie in the remaining ve cases. In contrast to Gendreau, Hertz and Laporte (1994), we improve three cases and tie in four. However, we need to point out that all these solutions by other heuristics were obtained by running with multiple parameter settings, therefore, our results yielded by a single pass appear to be more robust. We now extend the tests to eleven problems with integer distances using the same algorithm and parameter settings as for the real distances. The results are reported in Table 6. Similarly, we also list the time required by our algorithm, time required to obtain our nal solution, the best published solutions, and the times for obtaining solutions within 5% and 1% of the best. Again, the best solution for problem 1 (521) is known to be optimal by the analysis of Cornuejols and Harche (1993). Moreover, in two cases (problems 4 and 6), the best-known solutions were reported without indicating the speci c distance metric. Therefore, their true solution values for the rounded distance metric cannot be veri ed (see Gendreau, Hertz and Laporte, 1994).

30

Jiefeng Xu and James P. Kelly Time to Best Time to Time to Test Our Total Obtain Our Published Obtain Obtain Problem Solution Time Solution Solution 5% of Best 1% of Best 1 521 32.67 3.50 521 0.22 2.56 2 830 49.97 34.35 831 0.30 1.66 3 822 75.06 35.84 815 0.33 9.52 4 1024 136.63 49.87 1014 0.78 136 63 5 1296 246.65 228.77 1307 2.58 18.08 6 820 57.13 2.35 816 0.40 1.37 7 1034 137.33 104.71 1035 0.89 92.12 2b 1023 51.12 4.15 1022 0.37 3.81 2c 736 41.92 18.59 739 0.14 2.89 2d 685 42.63 1.40 690 0.12 0.74 3b 1080 75.41 34.91 1082 0.37 9.04 >

:

Table 6: Computational Results on 11 Problems with Integer Distances We observe from Table 6 that our algorithm yields good solutions for the real distance problems and can also produce high-quality solutions for the integer distance cases. For the nine problems whose best-known solutions are reliable, our algorithm improves the best published solutions for six problems, ties in one problem (optimal solution in this case) and produces inferior solutions in two problems. We list the new best solutions in the appendix. However, the performance gap due to these two dierent distance calculations is significant. The one which yields the best-known solution with one distance metric for a particular problem does not necessarily produce the best-known solutions with the other distance metric for the same problem. For example, although our algorithm produces the best-known solution for problem 3 with real distance, its integer solution is relatively poor. With this nding, we feel that fair comparisons include analysis using both real and integer distances. To further evaluate our solution quality with integer distances, we compare them with the most recent integer results from the literature (Gendreau, Hertz


31

and Laporte (1994), Harche and Raghavan (1991), Noon, Mittenthal and Pillai (1994) and Stewart, Kelly and Laguna (1994)). The integer solutions of Gendreau, Hertz and Laporte were listed in the earlier version working paper, but were deleted from the same published paper. The comparisons are presented in Table 7. Test Gendreau,Hertz Harche and Noon,Mittenthal Stewart,Kelly Our Problem and Laporte Raghavan and Pillai and Laguna Solution 1 521 521 521 521 521 2 832 847 845 831 830 3 815 825 831 815 822 4 1024 1070 1057 1025 1024 5 1316 1397 |{ 1307 1296 6 824 |{ 824 820 820 7 1035 1045 1037 1035 1034 2b |{ 1042 1042 1022 1023 2c |{ 751 739 739 736 2d |{ 695 690 690 685 3b |{ 1113 |1082 1080

Table 7: Comparisons with the Most Recent Integer Results As shown in Table 7, we nd that our algorithm produces far better solutions than those by Harche and Raghavan (1991) and Noon, Mittenthal and Pillai (1994). In contrast to Gendreau, Hertz and Laporte (1994), we produce four better solutions, tie in two case, and fail to improve on the remaining one problem. Again, we need to point out that the solutions reported by Gendreau, Hertz and Laporte were obtained by running multiple parameter settings whereas our solutions are obtained from a single set of parameters. Among eleven problems where comparisons are made with Stewart, Kelly and

32


Laguna (1994), our method improves on seven problems, ties in two cases and fails to improve on the remaining two. The results indicate that the network

ow-based TS approach for the VRP is extremely eective and ecient in determining high-quality solutions for a host of benchmark problems. Although our algorithm is empirically ne-tuned based on the real distance results, it appears to perform better on integer cases. This may be due to the fact that some of the best heuristics such as Taillard (1993) and Rochat and Taillard (1995) did not report integer results. (For example, the integer solutions for problem 4 and 5 are 1016 and 1275 respectively by calculating the best real distance solutions using the rounded distance metric.) In addition, for integer distance cases, we can improve several solutions using dierent parameter settings. The improved solutions are 1017 for problem 4, 1021 for problem 2b and 1076 for problem 3b. This indicates the potential for further improvement of our algorithm. We also list the new best-known solutions for problem 2b and 3b in the appendix for future veri cation. While the computational times represent an important aspect of the algorithm, direct comparisons may be misleading. Dierent platforms, operating systems, computer languages, and precisions impact computation times. Furthermore, clever data structures and coding skills may also change computation times. For example, Gendreau, Hertz and Laporte (1994) reported computational times from 6 to 91 minutes on a Silicon Graphics workstation, and Stewart, Kelly and Laguna (1993) ran on a DEC 5000 machine from 40


33

minutes to more than a hundred hours. Most recently, Rochat and Taillard (1995) reported average computational times on a Silicon Graphics Indigo (100 Mhz). Our method requires moderate computational times that range from three minutes to four hours on a DEC ALPHA workstation, and far less times to obtain the solutions within the 1% of the best for most problems. We compare the times of obtaining solutions within 5% and 1% of the best with those reported by Rochat and Taillard (1995) in Table 8. Though these two algorithms were not run on a same machine and cannot be directly compared, the correlation of the computing times between them is interesting. Test Times by Our Algorithm Times by Rochat and Taillard Problem 5% of Best 1 % of Best 5% of Best 1% of Best 1 0.65 0.77 0.05 0.18 2 0.05 6.72 0.03 1.13 3 0.15 2.38 0.25 15.00 4 4.09 38.41 0.42 30.00 5 6.22 201.18 1.48 | 6 0.09 2.39 0.85 5.83 7 0.52 1.18 11.50 45.00 Table 8: Computational Times Compared with Those of Rochat and Taillard (1995) We further test the algorithm on the three planar problems that appear in Fisher (1994). These problems are generated from real applications and have 44, 71 and 134 customers respectively. All problems use real Euclidean distances. Fisher (1994) found the optimal solutions for problems f44 and f71 and Rochat and Taillard (1995) located a new best-known solution for f134. We compare the results in Table 9.

34


Problem Our Solution Time Fisher's Solution Optimal or Best-known Solution f44 723.54 1.43 723.54 723.54 f71 244.54 86.65 241.97 241.97 f134 1176.72 191.46 1163.60 1162.96

Table 9: Tests on Fisher's Problems From Table 9, we see that our algorithm could not match the optimal or best-known solution for the two larger problems. However, the gaps are relatively small (within 1.2% of the best). For the problem with 134 customers, we can improve the solution to 1168.39 if we prolong the search to 20,000 iterations. Two factors may be responsible for not producing the optimal or best-known solutions for these two problems. First, based on Table 1, the standard deviations of the customer demands for Fisher's problems are much larger than the uniform demands that we tested above. Thus, our parameter settings for changing penalty functions, which are based on the uniform problems, are not appropriate here. In addition, Fisher's problems use fewer vehicles. Again due to the huge variations of the demands, some vehicles have far more customers (more than 40) than they have in the other uniform problems. In this case, our TSP heuristics may not be nding optimal or near-optimal tours.

2.3. Performance Analysis. Our algorithm produces high-quality solutions for a set of benchmark problems. We attribute this success to several factors. First, neighborhood network-based moves provide a powerful search mechanism. Second, by combining the network ow-based moves and the simple customer swap moves, the algorithm produces much better solutions than those


35

obtained by either method alone. Third, the tabu memories, especially the short term memory in this application, are eective in preventing the search from being trapped at local optima. Fourth, the dynamic penalty functions can better handle the tradeos between feasibility and solution quality. Finally, our unique restart/recovery strategy provides an eective intensi cation strategy to locate better solutions in the late stages of the search. We conducted some tests to identify the relative contributions of the various algorithm components. We limit our analysis to the performance on seven problems (1 through 7) with real distances. First, we report the experiments in which we (one at a time) disable the network ow moves, swap moves, tabu short term memory, restart/recovery strategy and TSTSP respectively. We compare our solutions (in Table 4) with the results of these experiments in Table 10. Test No Network No Swap No Short No Restart Problem Flow Moves Moves Term Memory & Recovery 1 552.29 527.98 524.61 524.61 2 876.18 851.53 835.89 837.52 3 866.80 831.56 833.00 831.68 4 1137.86 1050.41 1051.17 1037.32 5 1447.21 1337.04 1318.91 1320.45 6 956.39 819.56 819.56 819.56 7 1091.20 1043.37 1042.11 1042.11

No TSTSP 524.61 835.26 831.66 1036.39 1304.45 819.56 1071.95

Table 10: Comparison with Dierent Strategies

Our Solution 524.61 835.26 826.14 1029.56 1298.58 819.56 1042.11

36


The outcomes in Table 10 clearly demonstrate the dominance of our comprehensive algorithm, validating that the network ow moves, swap moves, TS Memory, restart/recovery strategy and TSTSP are necessary and eective components of our algorithm. In particular, the TS memory and restart/recovery strategy eectively help to locate extremely good solutions, and TSTSP provides an eective enhancement over 3-opt to solve TSPs. Our algorithm combines both network ow moves and swap moves with a certain frequency. Currently, we execute the swap moves once for every four iterations and several consecutive iterations under some speci c conditions. Table 11 includes the results for dierent frequencies of executing the swap moves. Test Every Two Every Three Every Five Every Six Problem Iterations Iterations Iterations Iterations 1 524.61 524.61 524.61 524.93 2 835.26 835.26 835.26 835.26 3 835.31 830.82 831.88 828.26 4 1045.99 1041.14 1041.21 1043.83 5 1314.30 1312.00 1304.96 1310.86 6 819.56 819.56 819.56 819.56 7 1042.97 1044.16 1044.41 1042.51

Our Solution 524.61 835.26 826.14 1029.56 1298.58 819.56 1042.11

Table 11: Comparison with Dierent Frequencies of Swap Moves


37

We nd from Table 11 that the performance of our algorithm is sensitive to the frequency which is used to combine the two dierent moves, especially for the two large problems (problems 4 and 5). Our algorithm performs best using a medium frequency.

3. Summary In this paper, we develop a new local search approach based on a network

ow model. This network ow model can determine several ejection/insertion moves simultaneously. We use this approach and a swap procedure to solve the well-known Vehicle Routing Problem. The capacity constraints are relaxed using penalty terms. Tabu Search is incorporated in the search procedure to overcome local optimality. Diversi cation strategies and restart strategies are discussed. Moreover, we design a dynamic search procedure that changes penalty parameters based on time and search feedback. We present our computational experience on a set of benchmark test problems and compare them with the best-known solutions in the literature. The algorithm determine several new best-known solutions with one parameter setting on test problems based on integer distance. The algorithm is competitive or outperforms most recent heuristics on both real distance and integer distance. One of the interesting ndings of this paper is the eectiveness of the network ow model that was used to implement the generalized and more composite local search moves. Complex constraints can be relaxed by incorporating the penalties into the cost terms.

38


Acknowledgements The rst author is grateful for the support from a 1994 Hart Summer Fellowship, Graduate School of Business, University of Colorado at Boulder. We also thank Gilbert Laporte and the other two anonymous referees for their helpful comments that signi cantly improved this paper.


39

References W. Bell, L. Dalberto, M. L. Fisher, A. Green eld, R. Jaikumar, R. Mack and P. Prutzman, \Improving the Distribution of Industrial Gases with an On-Line Computerized Routing and Scheduling System," Interfaces 13, 4-23 (1983). L. D. Bodin, B.L. Golden, A.A. Assad and M.O. Ball, \Routing and Scheduling of Vehicles and Crews: The State of the Art," Computers and Operations Research 10, 69-211 (1983). G. Brown, C. Ellis, G. Graves and D. Ronen, \Real-Time Wide Area Dispatch of Mobil Tank Trucks," Interfaces 17, 107-120 (1987). N. Christo des, \Vehicle Routing," in: E.L. Lawler, J.K. Lenstra, A.H.G. Rinnoy Kan and D.B. Shmoys, (eds.), The Traveling Salesman Problem. A Guided Tour of Combinatorial Optimization, Wiley, Chichester, 431-448, 1985. E. Hadjiconstantinou, N. Christo des and A. Mingozzi, \A New Exact Algorithm for the Vehicle Routing Problem Based on q-Path and k-Shortest Path Relaxations, Annals of Operations Research 61, 21-44 (1995). G. Clarke and J.W. Wright, \Scheduling of Vehicles from a Central Depot to a Number of Delivery Points," Operations Research 12, 568-581 (1964). G. Cornuejols and F. Harche, \Polyhedral Study of the Capacitated Vehicle Routing Problem," Mathematical Programming 60, 21-52 (1993). F. Cullen, J.Jarvis and H.D. Ralti, \Set Partitioning Based Heuristics for Interactive Routing", Networks 11, 125-144 (1981). Economic Report of the President, United States Government Printing Oce, Washington, 282 (1994). M.L. Fisher, \Optimal Solution of Vehicle Routing Problems Using Minimum K-Trees," Operations Research 42, 626-642 (1994). M.L. Fisher, A. Green eld, R. Jaikumar and J. Lester, \ A Computerized Ve-

40


hicle Routing Application," Interfaces 12, 42-52 (1982). M.L. Fisher and R. Jaikumar, \A Generalized Assignment Heuristic for Vehicle Routing," Networks 11, 109-124 (1981). B.A. Foster and D.M. Ryan, \An Integer Programming Approach to the Vehicle Scheduling Problem," Operational Res. Quart. 27, 367-384 (1976). M. Gendreau, A. Hertz and G. Laporte, \A Tabu Search Heuristic for the Vehicle Routing Problem," Management Science 40, 1276-1290 (1994). B. Gillett and L. Miller, \A Heuristic Algorithm for the Vehicle Dispatch Problem," Operations Research, 22,340-349 (1974). F. Glover, \Future Paths for Integer Programming and Links to Arti cial Intelligence," Computers and Operations Research 5, 533-549 (1986). F. Glover, \Ejection Chains, Reference Structures and Alterating Path Methods for Travelling Salesman Problems," Working Paper, School of Business, University of Colorado at Boulder, 1992. B.L. Golden, B.L. and A.A. Assad, Vehicle Routing: Methods and Studies, North-Holland, Amsterdam, 1988. Harche, F. and P. Raghavan, \A Generalized Exchanged Heuristic for the Capacitated Vehicle Problem," Working Paper, Stern School of Business, New York University, 1991 G. Laporte, \The Vehicle Routing Problem: An Overview of Exact and Approximate Algorithms," European Journal of Operational Research 59, 345-358 (1992). G. Laporte, G., H. Mercure, and Y. Nobert, \An Exact Algorithm for the Asymmetrical Capacitated Vehicle Routing Problem," Networks 16, 33-46 (1986). S. Lin, \Computer Solutions of the Traveling Salesman Problem", Bell System


41

Computer Journal 44, 2245-2269 (1965). S. Lin and B.W. Kernighan, \An Eective Heuristic Algorithm for the Traveling Salesman Problem," Operations Research 21, 498-516 (1973). R.H. Mole and S.R. Jameson, \A Sequential Routing-Building Algorithm Employing a Generalised Savings Criterion", Operations Research Quarterly 27, 503-511 (1976). C.E. Noon, J. Mittenthal and R.Pillai, \A TSSP+1 Decomposition Approach for the capacity-constrained Vehicle Routing Problem," European Journal of Operational Research 79, 524-536 (1994). I. Or, \Traveling Salesman-Type Combinatorial Optimization Problems and their Relation to the Logistics of Regional Blood Banking," Ph.D. Dissertation, Northwestern University, 1976. I.H. Osman, \Metastrategy Simulated Annealing and Tabu Search Algorithms for the Vehicle Routing Problem," Annals of Operations Research 41, 421-451 (1993). V.M. Pureza and P.M. Franca, \Vehicle Routing Problems via Tabu Search Metaheuristic," Publication CRT-747, Centre de recherche sur les transports, Universite de Montreal, 1991. Y. Rochat and E .D. Taillard, \Probabilistic Diversi cation and Intensi cation in Local Search for Vehicle Routing," Journal of Heuristics 1, 147-167 (1995). R.A. Russell, \An Eective Heuristic for the M-Tour Travelling Salesman Problem with Some Side Constraints," Operations Research 25, 517-524 (1977). W.R. Jr. Stewart and B.L. Golden, \A Lagrangian Relaxation Heuristic for Vehicle Routing," European Journal of Operational Research 15, 84-88 (1984). W.R. Jr. Stewart, J.P. Kelly and M. Laguna, \Solving Vehicle Routing Problems Using Generalized Assignments and Tabu Search," Working Paper, College of Business Administration, University of Colorado at Boulder, 1994.

42


E .D. Taillard, \Parallel Iterative Search Methods for Vehicle Routing Problems," Networks 23, 661-673 (1993). J.A.G. Willard, \Vehicle Routing Using r-Optimal Tabu Search," M.Sc. Dissertation, The Management School, Imperial College, London, 1989. J. Xu, S. Y. Chiu and F. Glover, Fine-Tuning a Tabu Search Algorithm with Statistical Tests, Working Paper, School of Business, University of Colorado at Boulder, 1996.


43

Appendix The solutions listed below show the sequences of customers on each route. The route length, route demand and total tour length are also listed. The distance is reported in rounded distance form. The depot is omitted from the solutions. The solutions for problems 2b and 3b are obtained via dierent parameter settings. Problem 2 (Q=140) Vehicle No. 1 2 3 4 5 6 7 8 9 10

30 74 21 61 28 2 68 75 4 27 52 34 67 46 8 54 13 57 15 5 45 11 66 65 38 58 10 31 25 55 18 50 40 35 19 59 14 53 7 26 6 33 63 23 56 24 49 16 29 37 20 70 60 71 69 36 47 48 62 22 64 42 41 43 1 73 51 17 3 44 32 9 39 72 12 Total length = 830

Route Length 74 36 83 77 116 86 92 102 99 65

Route Demand 140 135 135 128 140 135 140 135 138 138

Problem 7 (Q=200) Vehicle No. 1 2 3 4 5 6 7

82 111 86 85 89 93 96 94 97 115 110 116 103 107 104 99 101 102 106 105 109 21 20 23 26 28 31 35 32 29 36 34 33 30 27 24 22 25 19 16 17 95 37 38 39 42 41 44 46 47 49 50 51 48 45 43 40 52 54 57 59 65 61 62 64 66 63 60 56 58 55 53 100 98 68 73 76 77 79 80 78 75 72 74 71 67 70 69 87 92 91 90 114 18 118 108 83 113 117 84 112 81 119 120 88 2 1 3 4 5 6 7 9 10 11 15 14 13 12 8 Total length = 1034

Route Length 69 207 199 213 145 69 132

Route Demand 195 197 200 199 200 185 199

44


Problem 2b (Q=100) Vehicle No. 1 2 3 4 5 6 7 8 9 10 11 12 13 14

61 69 36 71 60 70 20 37 35 14 59 53 30 48 47 21 74 12 25 55 18 50 32 16 49 24 3 44 17 7 11 66 65 19 54 13 57 15 5 29 45 27 52 8 46 28 22 64 42 62 26 39 9 40 6 51 33 2 68 63 23 56 41 43 1 73 75 4 34 67 58 38 10 31 72 Total length = 1021

Route Length 115 79 61 93 73 76 86 47 95 52 49 84 25 86

Route Demand 98 87 99 100 99 98 99 100 98 96 94 100 99 97

Problem 2c (Q=180) Vehicle No. 1 2 3 4 5 6 7 8

45 15 57 13 54 19 59 14 35 8 26 12 72 38 65 66 11 53 7 6 62 28 61 21 74 2 68 75 30 48 47 36 69 71 60 70 20 37 5 29 67 34 46 52 27 4 51 3 44 32 9 39 40 17 16 49 24 18 50 25 55 31 10 58 33 63 23 56 41 43 42 64 22 1 73 Total length = 736

Route Length 117 86 79 102 38 65 132 117

Route Demand 167 179 175 178 142 166 179 178

Problem 2d (Q=220) Vehicle No. 1 2 3 4 5 6 7

17 40 3 44 32 9 39 72 12 8 19 54 13 57 15 37 20 70 60 71 69 36 5 29 45 67 46 34 52 27 4 51 16 49 24 18 50 25 55 31 10 58 26 6 33 63 23 56 41 43 42 64 22 62 1 73 38 65 11 66 59 14 53 35 7 30 48 47 21 61 28 74 2 68 75 Total length = 685

Route Length 70 135 38 133 123 105 81

Route Demand 171 219 142 209 215 209 199


Problem 3b (Q=112) Vehicle No. 1 2 3 4 5 6 7 8 9 10 11 12 13 14

61 16 86 38 44 91 98 1 51 9 81 33 79 50 58 2 57 42 14 43 15 41 22 74 73 31 10 63 90 32 30 70 69 80 24 29 78 34 35 71 65 66 20 12 68 3 77 76 28 52 7 19 47 48 82 18 8 46 36 49 64 11 62 88 59 93 85 100 92 54 55 25 39 67 23 6 96 99 37 5 84 17 45 83 60 89 40 21 72 75 56 4 26 94 95 97 87 13 53 27 Total length = 1076

Route Length 90 75 102 78 134 50 74 115 52 102 82 64 40 18

Route Demand 112 109 112 111 112 111 110 106 110 111 111 105 108 30

45