A deterministic tabu search algorithm for the fleet size and mix vehicle routing problem

Available online at www.sciencedirect.com

European Journal of Operational Research 195 (2009) 716–728 www.elsevier.com/locate/ejor

A deterministic tabu search algorithm for the fleet size and mix vehicle routing problem Jose´ Branda˜o Departamento de Gesta˜o, Escola de Economia e Gesta˜o, Universidade do Minho, Largo do Pacßo, 4704-553 Braga, Portugal CEMAPRE, ISEG, Technical University of Lisbon, Portugal Received 12 January 2007; accepted 31 May 2007 Available online 22 November 2007

Abstract The fleet size and mix vehicle routing problem consists of defining the type, the number of vehicles of each type, as well as the order in which to serve the customers with each vehicle when a company has to distribute goods to a set of customers geographically spread, with the objective of minimizing the total costs. In this paper, a heuristic algorithm based on tabu search is proposed and tested on several benchmark instances. The computational results show that the proposed algorithm produces high quality results within a reasonable computing time. Some new best solutions are reported for a set of test problems used in the literature.  2007 Elsevier B.V. All rights reserved. Keywords: Tabu search; Heuristics; Vehicle routing; Heterogeneous fleet; Fleet mix; Logistics

1. Introduction The fleet size and mix vehicle routing problem (FSMVRP) can be described mathematically as follows. Let G = (V, E) be an undirected connected graph with a vertex set V = {0, 1, . . . , n} and an edge set E = {(i, j): i, j 2 V}. Vertex 0 represents the depot and each other vertex i 2 Vn{0} is a customer with a non-negative demand qi. A distance dij (dii = 0, "i 2 V) is associated to each edge (i, j) 2 E. There is a fleet of T different types of vehicles located at the depot, and the number of vehicles of each type is considered unlimited. A capacity Qk, a fixed cost Fk, which is incurred by simply using a vehicle, and a variable cost vk are associated to each type of vehicle k (k = 1 ,. . . , T). We assume that Q1 < Q2,   , < QT and F1 < F2,   , < FT. The travelling cost of each edge (i,j) 2 E by a vehicle of type k is cij = vkdij. The FSMVRP consists of defining a set of routes and the vehicles assigned to them so that the following constraints are taken into account: (i) satisfy customers’ demand; (ii) visit each customer exactly once; (iii) a vehicle route starts and finishes at the depot, and the capacity of E-mail address: [email protected] 0377-2217/$ - see front matter  2007 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2007.05.059

the vehicle is not exceeded. The objective of the FSMVRP is to minimize the sum of fixed and variable costs of all the routes subject to the previous constraints. A particular case of FSMVRP happens when T = 1, i.e., when the fleet is homogeneous, which is the classical vehicle routing problem (VRP), should the fixed costs not be taken into consideration. Therefore, we can conclude that, like the VRP (this is proven, for example, in Lenstra and Rinnoy Kan (1981)), the FSMVRP is a NP-hard combinatorial problem. The FSMVRP has been studied under several other designations like vehicle fleet composition (Etezadi and Beasley, 1983; Salhi and Rand, 1993), mix fleet vehicle routing problem (Wassan and Osman, 2002) and heterogeneous fleet vehicle routing problem (Gendreau et al., 1999; Yaman, 2006; Choi and Tcha, 2007). On the other hand, Taillard (1999), Tarantilis et al. (2004) and Li et al. (2007) studied a variant of the FSMVRP with a fixed number of vehicles for each type. Another variant of this problem, the FSMVRP with time windows has been studied by researchers such as Liu and Shen (1999) and Dullaert et al. (2002). In practice, many companies face distribution problems similar to the FSMVRP or, more commonly, with some

J. Branda˜o / European Journal of Operational Research 195 (2009) 716–728

additional constraints. The need for different types of vehicles is determined by the characteristics of the customers. Usually, larger vehicles are more appropriate for serving customers who require large orders and the smaller vehicles are more adequate to deliver small quantities or serve customers that have access restrictions. The FSMVRP may be present either when making strategic decisions or in the operational ones. The strategic (or tactical) decision is when the fleet is chosen and the company must decide which type of vehicles to buy and how many of each type, i.e., the company wishes to define the best composition of the fleet. In this case, the most relevant costs to be included in the objective function are the fixed costs of the vehicles. The operational decisions appear in the daily planning, once the fleet of vehicles has been defined, since the total costs will depend on the number and type of vehicles effectively used for the deliveries, as well as on the routes in which they are used. Besides, when the owned vehicles are not enough, a decision has to be made regarding which types of vehicles should be hired and how many of each. For the daily planning the most important costs to be considered in the objective function are the variable costs, although the fixed costs may also be included, especially in relation to the hired vehicles, depending on the kind of contract. The VRP has received a lot of attention and there is a vast amount of papers and books about it. As a consequence of this, there has been an enormous progress in its resolution. A survey of the existing resolution methods may be found, for example, in Laporte (1992), Golden et al. (1998) and Laporte et al. (2000), and an extensive bibliography is presented in Laporte and Osman (1995). Contrarily, the FSMVRP, in spite of its practical importance, has been much less studied. Among these, deserve mention the works of Semet and Taillard (1993), Rochat and Semet (1994) and Branda˜o and Mercer (1997), which are some examples of real-life problems that include the FSMVRP with several additional constraints. In the case of Branda˜o and Mercer (1997), the company owned a fleet of two types of vehicles and, whenever necessary, they could complement their fleet with hired vehicles. Some of the additional constraints were the following: some customers could only be served by one vehicle type, time windows for serving the customers, each vehicle could perform more than one route per day (a more detailed description can also be found in Branda˜o (1994)). Below, we will briefly review some of the literature on FSMVRP resolution methods. As far as we know, no exact methods have been developed for the FSMVRP. However, three papers dealing with the determination of lower bounds deserve mention: Golden et al. (1984), Yaman (2006) and, more recently and with better results, Choi and Tcha (2007). In the latter, the authors were able to find lower bounds, with a gap, in relation to the best known solution, ranging from 0.17% to 4.32% for a set of 36 problems from the literature, containing between 20 and 100 customers. The average gap was only 1.71%.

717

Considering the similarities between the VRP and the FSMVRP, no wonder that most of the methods developed for solving the latter are adaptations of methods already widely used for the former. Golden et al. (1984) have created several heuristics based on the savings method of Clarke and Wright (1964), as well as on the giant tour algorithm of Gillett and Miller (1974). Gheysens et al. (1986) adapted the generalised assignment heuristic of Fisher and Jaikumar (1981). Desrochers and Verhoog (1991) have adapted a matching based savings heuristic initially proposed by themselves for the VRP. The heuristic of Salhi and Rand (1993) is also based on a previous one developed by themselves for the VRP. This heuristic was later improved by Osman and Salhi (1996). All these methods, more or less sophisticated, can be classified as classical heuristics. The most recent heuristic of this type is due to Renaud and Boctor (2002) and is, among them, the one that produces better quality solutions, but also the one that requires more computing time. This heuristic starts by generating a large set of routes using different procedures and afterwards chooses those that satisfy the constraints of the problem at the lowest cost using an exact, but polynomial, set partitioning algorithm. Among those procedures of producing the initial set of routes is the petal method used by several authors, including Renaud et al. (1996), for the VRP. More recently, meta-heuristics have been proposed for solving the FSMVRP: Osman and Salhi (1996), Gendreau et al. (1999) and Wassan and Osman (2002), which will be designated from now on as OS, GLMT and WO, respectively. All these algorithms are based on tabu search. The OS algorithm takes an initial solution produced by the heuristic of Salhi and Rand (1993), which is then improved by a tabu search method with short term memory and with moves defined by a 1-interchange mechanism. The GLMT algorithm is rather complex and requires the use of GENIUS, developed by Gendreau et al. (1992) for the travelling salesman problem (TSP), as well as an adaptive memory procedure developed by Rochat and Taillard (1995). The WO algorithm comprises several variants obtained from the selection of different neighbourhood mechanisms, tabu restrictions and tabu tenure schemes. Very recently, Choi and Tcha (2007) developed a mathematically based resolution method for the FSMVRP. In this method, the FSMVRP is formulated as a set covering problem as it is done for VRP with time windows. Then, its linear programming relaxation is solved by the column generation technique. This method could not solve exactly any of the examples tested, but it generated good lower bounds and upper bounds and, therefore, can be used as a good heuristic method. The structure of the remainder of this paper is as follows. Section 2 describes the methods used to obtain initial feasible solutions in our algorithm. Section 3 describes the tabu search algorithm (TSA) for solving the FSMVRP, presenting its main features and a flow chart that shows how the different components interact. In Section 4, the proposed parametric setting is discussed, and are given

718

J. Branda˜o / European Journal of Operational Research 195 (2009) 716–728

the results obtained on a set of test problems. Final conclusions are presented in Section 5. 2. Initial solutions In order to be applied, the TSA requires an initial solution. Three different methods are used to obtain initial solutions. Each method produces a feasible solution that can be used as an initial solution for the TSA. These methods were designed to be fast to compute and to provide a variety of starting solutions so that the TSA may improve. The diversity provided by the different starting solutions was found to be useful in ultimately finding high quality solutions. In order to clarify the TSA description that will follow, we should take into consideration that, in any part of the algorithm, the ties in terms of distance, cost, frequency of move, etc., when they exist, are always arbitrarily broken. 2.1. Method 1 – nearest neighbour First, the vehicle type is selected, beginning with type 1. If none of the unrouted customers is admissible in terms of the load carrying capacity of the vehicle, then that type of vehicle is increased by a unit, and this is repeated until at least one of the unrouted customers can be served by this type of vehicle. Afterwards, the route that will be assigned to the chosen vehicle starts with the unrouted customer farthest from the depot. The next customer to be inserted in the route, will be the one who has not been served yet, and who is nearest to the customers of the route. That customer also has to be admissible in terms of the vehicle’s load carrying capacity. The selected customer is inserted in the route before or after its nearest neighbour according to the place that increases the travelling distance by the least amount. This process is repeated until no customer is admissible in the current route. When this happens a new vehicle is selected and the whole process is repeated until all the customers are routed. The GENIUS algorithm of Gendreau et al. (1992) is composed of two procedures: GENI and US. The GENI is used to construct a TSP tour and the US to improve this tour. In this method, after constructing all the routes, the US is applied to each of the routes in order to try to reduce the travel distance. 2.2. Method 2 – nearest neighbour plus insert This method operates exactly the same way as method 1, except in the two points that follow. After inserting the first customer in a route (the vehicle assigned to the route was already defined as in method 1), the admissibility of any other customer in this route requires as an additional condition that it belongs to the d-neighbourhood. This neighbourhood is constituted by the d-nearest neighbours of each of the customers already

in the route. The value of d is defined as the average estimated number of customers per route, which is obtained assuming that only vehicles of type 1 are used. Firstly, the number of vehicles required is estimated by dividing the total demand by the capacity of the vehicle. Then, we get the value of d by dividing the total number of customers by the number of vehicles. When no more customers are admissible in the current route, the insertion phase is applied. In this phase, the unrouted customer that increases least the cost of the current route is inserted each time, and that is admissible in terms of capacity. When no more customers can be inserted, the whole process is repeated with a new vehicle (route) until all the customers are routed. 2.3. Method 3 – giant tour Initially, a TSP tour including all the customers is constructed. The tour starts with the depot and each time a new customer is inserted, following their order number, using the GENI procedure. The partition of the tour into feasible routes is performed as follows: the partition starts with the customer adjacent to the depot (the tour is represented as a vector whose first element is the depot and, therefore, the customer selected is its second element); then, the vehicle type with the minimum capacity equal or greater than the customer’s demand is selected, and, therefore, the current route is started; in this route, the customers that are admissible in terms of capacity are included, one by one, following the sequence defined by the tour. When no more customers are admissible, the process is repeated with the shorter tour that results from eliminating the customers already routed. This process ends when there are no more customers in the tour. 3. The tabu search algorithm 3.1. Neighbourhood structure The tabu search algorithm is based on three types of neighbourhood moves: the first two are insertions (single and double) and the third one is a swap. In each iteration all the customers are candidates for being moved, except in the frequency phase that will be defined later on. In a single insertion move, a candidate customer xi is removed from its current route (origin) and a trial insertion is made in any other route (destination), containing at least one of the d-nearest neighbours of xi, between any two customers or near the depot. The value of d varies during the search and will only be defined later on. This d-neighbourhood restriction not only reduces the computing time but also contributes to find better solutions. We have used this concept, although with different settings, for the first time in Branda˜o and Mercer (1997). Gendreau et al. (1994), in their algorithm for solving the VRP, also imposed this neighbourhood condition to allow a customer to be moved to another route.

J. Branda˜o / European Journal of Operational Research 195 (2009) 716–728

In a double insertion move, the operation is similar to the single insertion one, except for the following: the candidates are any pair of customers (xi, xj) belonging to the same route, the destination route should have at least one of the d-nearest neighbours of xi or xj and the insertion place of xj takes into account also the previous (virtual) insertion of xi. Both the single and double insertion can originate a new route, subject, however, to the following two conditions: (i) the origin route should have at least 3 or 4 customers for the single and double insertion, respectively; (ii) the type of vehicle of the destination route should be the lowest capable of containing the demand of the moving customers or, in the case of the double insertion, the largest vehicle if the sum of the demand of the two customers exceeds the capacity of any type of vehicle. A swap move consists of exchanging two customers belonging to two different routes. A trial swap of a customer xi 2 ri (route ri) with a customer xj 2 rj exists if and only if: (i) at least one customer of rj belongs to the d-neighbourhood of xi; (ii) at least one of the routes ri or rj contains more than one customer. In this implementation, the frequencies of the different types of move may change according to the cycle and phase of the algorithm. Parameters FD and FS denote the frequencies of the double insertion and swap moves, respectively. For example, the value of FS implies that the swap move is only tested every FS iterations. The trial move chosen depends on the effect of the move on the objective function defined in the next section. In the designing of the TSA we took advantage of the knowledge that we acquired from our previous tabu search algorithms for other types of routing problems, namely Branda˜o (2004), Branda˜o (2006) and Branda˜o and Eglese (2006). Therefore, the TSA has several characteristics in common with these algorithms and, in particular, with the last one, because both are deterministic and both have the same type of neighbourhood moves (although being executed with different frequencies). However, the TSA contains two important features that do not appear in Branda˜o and Eglese’s (2006) algorithm: the d-neighbourhood restriction and the frequency phase. Furthermore, this algorithm is for the arc routing problem which implies that the route improvement procedures must be completely different. In the TSA, we also tried other types of moves, namely the exchange of more than one customer between each pair of routes. The results were worse in terms of both solution quality and computing time and, consequently, these types of moves were abandoned. 3.2. Objective function

f ðSÞ ¼

r X

ðF i þ ci þ Pli Þ:

719

ð1Þ

i¼1

where r is the total number of routes in S, Fi is the fixed cost of the vehicle assigned to route i, ci is the variable cost of route i, P is a penalty term and li is the load excess in route i. The strategic oscillation created by allowing the search to cross the boundary between the feasible and the infeasible space plays an important role in the quality of the final solutions. We tried with the same algorithm, but without allowing infeasible solutions to occur, and the results were rather worse. We should note that for any S the value of li depends on the type of vehicle assigned to route i. Therefore, in each trial move, the type of the vehicle can be changed to one with more or less load carrying capacity, and the vehicle type that implies the lowest value of f(S) will be the one chosen. All this takes place in the insertion (single or double), as follows. Let us suppose that ri is the origin route after removing the customer (or customers), and rj is the destination route after inserting the customer (or customers). If ri is infeasible nothing is done; otherwise, the vehicle type is reduced while ri is still feasible. If rj is feasible nothing is done; otherwise, the vehicle type is increased by one unit at a time until rj becomes feasible or type of vehicle with the highest load carrying capacity is reached. At the same time, for each vehicle type, including the initial one, the cost of the route is calculated by (1) and, in the end, the vehicle type to which the lowest cost corresponds is chosen. Regarding the trial swap move, after performing the move, the same procedure is applied to both routes rk and rm: if rk (rm) is feasible or infeasible, it is handled as ri or rj in the insertion, respectively. 3.3. Admissibility of moves A conventional tabu list is created to prevent a customer that has been moved to return to its original route for the next h iterations, where h is the tabu tenure. The tabu tenure changes systematically during the search according to the evolution of the results. The tabu restriction may be overridden if the move produces a solution that is better than the ones found in the past. This is referred to as the aspiration criterion. In the TSA, a trial move to solution S0 is regarded as admissible: (i) if it is not currently on the tabu list; (ii) or if it is tabu and infeasible, but the value of f(S0 ) is lower than the value of the best infeasible solution yet found; (iii) or if it is tabu and feasible, but the value of f(S0 ) is lower than the value of the best feasible solution yet found. 3.4. Route improvement procedure

For any candidate solution, S, the objective function is denoted by f(S). The objective function to be minimized by the TSA is defined by the next equation:

At the end of each iteration, the customers that had been moved are introduced into their new routes by the GENI

720

J. Branda˜o / European Journal of Operational Research 195 (2009) 716–728

procedure. This may result in a better solution because the route (or routes in the swap move) is reorganised, while the cost determined in the trial moves was based solely on a simple insertion. After this, if a new best feasible solution is found in this iteration, the US procedure is applied to each individual route, in order to try to reduce its cost further.

2

Perform the trial moves and choose the best admissible

3.5. Diversification and intensification

3

Improve changed routes Im using GENI

4

Update: tabu list, best feasible and infeasible solutions, iteration counters

The main sources of diversification in the TSA are: the use of several different initial solutions and phases; the strategic oscillation, as well as the change of the vehicle type driven by the objective function; the change of the tabu tenure, the neighbourhood and the frequency of the types of moves, during the search. The intensification is mainly achieved using the GENIUS algorithm and the restart, after a given number of iterations, with the best found solution. The fact that several parameter values, such as the duration of each cycle (defined in the next section), depend on the evolution of the search is also important.

Start phase 1

Generate an initial solution

1

S

Yes

Reached the stopping criterion? No Is in the middle of the cycle?

Yes

Change θ, δ, FD and FS

5

Empty tabu list

6

No No

Is in the end of the cycle?

Yes

Stop

3.6. Overview of the TSA

Fig. 1. Flow chart of the TSA.

In Fig. 1, the flow chart of the TSA is presented, which shows how the different elements, described in the previous sections, work together. A more detailed description is given in Appendix B. Each TSA procedure that is identified in Fig. 1 by a number has the following meaning:

Procedure 6 – empty the tabu list, update the iteration counters and set S equal to the best solution already found.

Procedure 1 – the initial solution is given by one of the algorithms described in Section 2. Procedure 2 – the execution of the trial moves follows the following sequence: single insertion, double insertion and swap, performing all the potential moves of each type and then proceeding to the next. The trial single insertion moves are perform in every iteration, while the other two types only when the current iteration is a multiple of the frequency set for that type. If in a given iteration a trial move generates a solution better than the best already discovered no further trial move is performed in that iteration. Procedure 3 – the customers that have been moved are introduced into their respective routes using the GENI procedure. Procedure 4 – the tabu list is updated as well as the best known feasible and infeasible solutions and several iteration counters: total number of iterations, number of iterations without improving the best feasible or infeasible solutions, number of consecutive feasible or infeasible solutions. Procedure 5 – the values of the parameters h, d, FD and FS are changed.

The set of procedures 2–6 that is executed in each iteration constitutes the basic structure of the TSA – the cycle. A cycle ends after a given number of iterations (KBL) without improving the best feasible or infeasible solution (SB). Therefore, in spite of KBL being fixed in each cycle, the duration of the cycle varies according to the search, because each time SB is improved the counter KB restarts from zero. Besides, the value of KBL is increased from one cycle to the next. The middle of the cycle is reached when KB = KBL/2 (note that in each cycle the middle can be reached more than once). A phase is a set of cycles. Any cycle and, therefore, the phase can be interrupted if the stopping criterion is reached. The TSA contains four phases. Phase 1 is just the one described in Fig. 1. The other three phases start with the best feasible solution found in the previous phase labelled in Fig. 1 as S. Apart from this, they only differ from phase 1 in the initial values of h and d and in the limit for the total number of iterations. Phase 4 is a frequency phase because it takes into account the frequency of moves of each customer during all the previous iterations of all the phases. So, phase 4 has two additional differences: in each iteration only a fraction of the customers is candidate for being moved and the neighbourhood restriction is removed. Two versions of the TSA have been applied to the experimental set of problems: Version 1 and Version 2.

J. Branda˜o / European Journal of Operational Research 195 (2009) 716–728

In Version 1 the initial solution is taken from the giant tour method as described in Section 2.3, because, on average, this method gave the best final results. In Version 2 the three initial solutions described in Section 2 are used. Therefore, for each initial solution, the four phases just described are carried out. Additionally, Version 1 is also included. In the end, another phase identical to phase 2 is executed, but taking as initial solution the best solution found before and with different values of h and d. 4. Computational experiments The definition of the value of the TSA parameters was based on our experience with other types of tabu search algorithms and on a set of experiments and principles described below. First of all, the most important feature that distinguishes different instances of the FSMVRP is the number of customers, n. Therefore, the value of most of the parameters is a function of n. This rule ensures to some extent a similar performance of the TSA with problems of different sizes. The two most influential parameters of the TSA are the tabu tenure (h) and the neighbourhood restriction (d). The initial value of h is given by n/c (a number is always rounded up to the nearest integer, if the parameter must be integer), where c is a constant. We tried with the following values of c: 1, 1.5, 2, 3 and 4. The best average results were obtained with c = 2. During the search h changes dynamically in the following way: whenever the middle of the cycle is reached, it is set as h = max(h/2,5) and at the end of the cycle it is defined as h = 2h + 3. In order to have influence in the search, d must be small. Its maximum value, dL, takes into account the number of customers and the number of vehicles required to serve them. So it was defined as dL = min(n  1, 2cr), where cr is the number of customers of the route of the initial solution with more customers. As it happens with the tabu tenure, d also changes dynamically during the search as follows: in the middle of the cycle, d = min(d + 1, dL) and at the end of the cycle, d = min(d + 2, dL). A set of experiments indicated that the initial values should be very low. In the case of Version 1 the best average results were obtained using d = 1 in phase 1. The initial value of d increases a bit in the next phases (for example, it is 2 in phase 2) in order to allow the exploration of different and larger regions of the solution space. This is also why in the last phase (frequency phase) no neighbourhood restriction is imposed; it allows the search of moves that could not be performed in the previous phases due to this restriction. The GENIUS algorithm depends on a parameter, P, that defines the neighbourhood around each customer, and we used P = 5 because according to the authors (Gendreau et al., 1992) it gives a good compromise between solution quality and computing time.

721

The value of P was defined simply according to our previous experience with other types of tabu search algorithms. Initially, P = 1. If for 10 consecutive iterations, all the solutions are infeasible, P is multiplied by 2. Conversely, if all the solutions are feasible for the last 10 consecutive iterations, P is divided by 2. The restart associated with the cycle plays an important role in the performance of the TSA because it promotes intensification around the best solution found and some diversification is also guaranteed by changing parameters’ values. If the cycles are too long there is not enough time for restarting and if they are too small the search is limited to a narrow region of the solution space. According to our experiments, good results are obtained when the cycle limit (KBL) is around 15n. This limit is increased by 2n at every restart to allow some more time to find better solutions because this becomes harder as the search progresses. The execution of each phase of the TSA is stopped when the limits of the total number of iterations, KL, and number of iterations without improving the best feasible solution, KBFL, are both reached. In the definition of these limits we tried to obtain a good compromise between solution quality and computing time. Furthermore, KBFL is the double of KBL in order to allow the search to continue during a reasonable number pffiffiffi of iterations after the restart. In phase pffiffiffi 1, K L ¼ 3000 n and in the other phases K L ¼ 2000 n. Usually, these phases contribute much less to improve the solution and, therefore, in general, would be a waste of computing time if the value of KL was identical to that of phase 1. The other component of the stopping criterion (KBFL) is adapted to the evolution of the search, allowing this to go on whenever new best feasible solutions are discovered, because the corresponding iteration counter, KBF, restarts from zero when this happens. The trial single insertion moves are executed in every iteration. In relation to the double insertion and swap moves, the following frequencies were defined, respectively: FD = 10 and FS = 5, in the beginning; FD = 5 and FS = 10, in the middle of the cycle; FD = 10 and FS = 7, in the end of the cycle. As a rule, these two types of trial moves are not performed in the same iteration (see Appendix B, step 2) so that they do not compete with each other. This variation of the frequencies during the cycle induces a diversification in the search. We observed that both these two types of moves give an important contribution for discovering good solutions. However, if their frequency is very high (for example, every iteration) not only the computing time increases a lot but also the final results are worse. This behaviour may be explained as follows: these moves drive the search towards a region of the solution space which is rather distant from the current one and, if this happens frequently, a high perturbation is created acting against the proximate optimality principle – POP – (Glover and Laguna, 1997, pp. 138–141). Conversely, the single insertion move allows the exploration of the same ‘‘level”, according to the concept defined by the POP. We should note that the perturbation created

722

J. Branda˜o / European Journal of Operational Research 195 (2009) 716–728

Table 1 Variable costs of the vehicles Vehicle type

A B C D E F

Problem number

vA vB vC vD vE vF

23

24

25

26

27

28

29

30

1.0 1.1 1.2 1.7 2.5 3.2

1.0 1.1 1.4

1.0 1.6 2.0

1.0 1.6 2.1

1.0 1.2 1.5 1.8

1.0 1.3 1.9 2.4 2.9 3.2

1.0 1.4 1.7

1.0 1.7 2.0

by those types of moves is boosted by the existence of different types of vehicles. The TSA has been programmed in C language and executed on a Compaq Presario X1000 Intel Pentium M at 1.4 GHz, 256 MB of RAM. The performance of our algorithm was tested using two sets of benchmark problems from the literature. The first set was taken from Golden et al. (1984) and contains 20 problems, numbered from 1 to 20. The second contains eight problems, numbered from 23 to 30, and their data are just the same as those for the corresponding problems of the first set from 13 to 20, except that in this second set, the vehicles have associated variable costs, instead of fixed costs. These variable costs, presented in Table 1, have been proposed by Taillard (1999). We compare the results produced by our algorithm with those given by the algorithms of Taillard (1999), Gendreau et al. (1999), Wassan and Osman (2002) and Renaud and Boctor (2002), because the first three seem to be the best published algorithms using meta-heuristics and the last one seems to be the best among the classical heuristics. The mathematically based algorithm of Choi and Tcha (2007) is also used in our comparisons. The first three algorithms were executed on Sun Sparc workstations at 50 MHz, while the fourth was executed on a Pentium II at 233 MHz, and Choi and Tcha (2007) used a Pentium IV at 2.6 GHz, 523 MB of RAM. In order to get a rough idea of the relative speeds of the different computers, we use Dongarra (2006) tables and, additionally, our own experience with similar computers, because there are many models that do not

appear in these tables. Therefore, being aware that this is just a rough estimate, we will assume that the speed of the computers referred to before, measured in millions of floatingpoint operations per second (Mflop/s), is the following: workstations – 10 Mflop/s, Pentium II – 40 Mflop/s, Pentium IV – 1190 Mflop/s and Pentium M – 250 Mflop/s. Taillard (1999) and Gendreau et al. (1999) did not solve the problems 1–2 and 7–12 because their algorithms were not designed to handle these problems, since the coordinates of the customers are not given in the data. On the other hand, Renaud and Boctor (2002) did not solve the problems 23–30, with variable costs associated to the vehicles. So, we decided to present the results in three separate tables. Both Wassan and Osman (2002) and Renaud and Boctor (2002) developed more than one version of their algorithms and, therefore, we will compare our results with the versions that yield better results: the RTS-F and the standard sweep, respectively. Moreover, these authors used different ways of presenting their results: Taillard (1999) executed the algorithm five times, giving the average solutions and computing times; Gendreau et al. (1999) also provide the averages from ten runs, and the computing time presented is for the best run, but we take this to be the average time because it should be similar; Wassan and Osman (2002) run the algorithm five times for the first set of problems and twice for the second, presenting the best solutions found, but the computing time is only the time of the best run, for each instance; Choi and Tcha (2007) executed the algorithm five times with different parameters’ values in each run and, for each problem, present the average solution cost and the average computing time. The tables give the following information: the problem number (No.); the number of customers; the best known solution costs from the literature (although, for seven instances, the costs have been recalculated, using double precision (64 bits), from the solutions published, because the values reported in the main text are not according to these solutions), with the indication of the source in the case the solution was found by only one algorithm; ‘‘our best” which is the best result obtained by our algorithm

Table 2 Results for eight problems without variable costs No.

n

1 12 2 12 7 30 8 30 9 30 10 30 11 30 12 30 Average Average deviation (%) NB

Best known

Our best

602 722 7273a 2347 2209 2355a 4755 4087a 3043.75

602 722 7273 2347 2209 2355 4755 4080 3042.88 0.03 8

Wassan and Osman

Renaud and Boctor

Version 1

Cost

CPU

Cost

CPU

Cost

CPU

Cost

CPU

602 722 7273 2347 2209 2355 4755 4087 3043.75 0 8

17 16 113 288 358 229 342 458 228 – –

602 722 7346 2347 2211 2362 4765 4092 3055.88 0.40 3

0 0 28 44 32 10 34 19 21 – –

602 722 7278 2347 2209 2355 4756 4092 3045.13 0.05 5

1 2 15 13 20 11 18 11 11 – –

602 722 7273 2347 2209 2355 4755 4080 3042.88 0.03 8

7 6 59 52 73 60 62 56 47 – –

Values in boldface – identical to the best known or better (in this last case they are also in italic). a Solution cost taken from Wassan and Osman (2002).

Version 2

21 22 20 25 145 220 110 111 322 267 438 601 192 – – 470 570 334 349 2072 2744 12528 2117 2648d – –

Values in boldface – identical to the best known or better (in this last case they are also in italic). a Solution cost taken from Choi and Tcha (2007). b Solution cost taken from Taillard (1999). c These problems were solved using a non-standard version of the algorithm. The CPU for solving these problems was between 20 and 40 seconds. d For the eight problems solved with the standard version of the algorithm.

Cost

961.03 6437.33 1007.05 6516.47 2406.36 9119.03 2586.84 2728.14 1736.09 2376.89 8667.26 4048.09 4049.21 0.04 7 4 5 4 6 34 43 24 24 73 47 93 132 41 – –

CPU Cost

961.03 6437.33 1007.05 6516.47 2406.36 9120.11 2586.84 2739.65 1744.41 2378.39 8668.24 4085.70 4054.30 0.16 6 4 6 5 9 50 160 45 28 652 1037 1110 307 284 – –

CPU Cost

965.6 6467.01 1009.14 6588.93 2412.51 9128.98 2618.03 2767.29 1757.72 2419.92 8687.31 4102.82 4077.11 0.73 0 0 1 1 0 10 51 10 11 207 70 1179 264 150 – –

CPU Cost

961.03 6437.33 1007.05 6516.84 2409.77 9119.13 2588.92 2731.08 1746.26 2375.78 8665.08 4057.33 4051.30 0.09 3 88 80 52 88 2084 1660 2349 689 1874 2261 8570 2692 1874 – –

CPU Cost

961.03 6437.33 1007.05 6516.47 2422.10 9119.86 2586.37 2730.10 1755.10 2385.52 8665.75 4061.64 4054.03 0.16 5 164 253 164 309 724 1033 901 815 1022 691 1687 1421 765 – –

CPU Cost

961.03c 6437.33c 1008.59c 6516.47c 2436.78 9123.60 2593.61 2744.25 1753.74 2382.80 8665.40 4063.18 4220.42d 0.34d 3

961.03 6441.01 1008.72 6517.98 2424.88 9121.98 2590.68 2743.96 1752.29 2392.57 8682.50 4100.20 4061.48 – 1

CPU Cost

961.03 6437.33 1007.05 6516.47 2406.36 9119.03 2586.37 2728.14 1734.53 2369.65 8661.81 4042.59 4047.53 0.003 10

Cost Cost

961.03 6437.33 1007.05 6516.47 2406.36a 9119.03 2586.37 2720.43a 1744.83a 2371.49a 8661.81b 4039.49a 4047.64

Version 1 Renaud and Boctor Choi and Tcha Wassan and Osman Gendreau et al. Taillard Our best Best known n No.

Table 3 Comparison of the results for 12 problems without variable costs

3 20 4 20 5 20 6 20 13 50 14 50 15 50 16 50 17 75 18 75 19 100 20 100 Average Average deviation (%) NB

Version 2

CPU

J. Branda˜o / European Journal of Operational Research 195 (2009) 716–728

723

during all experiments with different values of the parameters; the cost of the solution produced by the different algorithms and respective computing time in seconds (CPU) on the computers previously indicated; the average for each set of problems; the percentage difference of the average in relation to the average of the best known solution costs (Average deviation); the number of solutions found by each algorithm which are identical or better than the best published ones (NB). The results in Table 2 show that the TSA was able to find all the best known solutions and a new best one. These problems are rather easy to solve because they have a small number of customers but, in spite of this, a new best solution was discovered. In comparison with Wassan and Osman’s (2002) algorithm, Version 1 is slightly worse while Version 2 is slightly better, finding all the best known solutions and a new best one. In terms of computing time, considering the fact that Wassan and Osman’s (2002) results were found in five runs and the relative speeds of the computers, we can conclude that Version 1 is about four times faster and Version 2 is almost as fast. Comparing with Renaud and Boctor (2002), we can see that both versions of the TSA yield substantially better solutions, but the computing time is approximately three and fourteen times higher for Version 1 and Version 2, respectively. The results in Table 3 show again that the TSA performs better than any of the other algorithms. Comparing to Taillard (1999) (considering the average for the 8 problems solved with the standard version of the algorithm), Version 1 produces slightly better results and is a bit faster. However, Taillard’s (1999) algorithm is substantially faster than Version 1 on problems with 50 customers or less. The results of Version 2 are also better than those of Taillard (1999), though more than twice as slow. Gendreau et al.’s (1999) average costs are 0.18% and 0.30% higher than those of Version 1 and Version 2, respectively, while the computing time of Version 1 is a bit higher and the computing time of Version 2 is about six times higher. Comparing with the solutions of Wassan and Osman (2002), we can see that Version 1 produces almost identical average costs but the computing time is about nine times lower for Version 1, while Version 2 yields an average cost 0.12% lower and it is almost twice as fast. Choi and Tcha’s (2007) average costs are 0.07% lower and 0.05% higher than those of Version 1 and Version 2, respectively. However, Choi and Tcha’s (2007) computing times are substantially higher: about seventeen times and almost four times higher than those of Version 1 and Version 2, respectively. As far as the algorithm of Renaud and Boctor (2002) is concerned, both versions of the TSA found substantially better solutions. In terms of computing time, Renaud and Boctor’s (2002) algorithm is a bit slower than Version 1 and about four times faster than Version 2. The largest differences, in terms of computing time, are observed in the problems containing 75 customers or more. This allows us to conclude that for medium or large problems Version 1 not only yields better solutions but it is also faster than a good

724

J. Branda˜o / European Journal of Operational Research 195 (2009) 716–728

Table 4 Comparison of the results for eight problems with variable costs No.

23 24 25 26 27 28 29 30 Average Average NB

n

50 50 50 50 75 75 100 100

Best known

Our best

Taillard

Cost

Cost

Cost

1491.86 603.21 999.82 1131.00 1038.60 1800.80 1105.439 1530.43 1212.64 0.007 8

1494.58 603.21 1007.35 1144.39 1044.93 1831.24 1110.96 1550.36 1223.38 0.88 1

1491.86 603.21 999.82 1131.00 1038.60 1801.40 1105.443a 1530.43b 1212.72 deviation (%)

Gendreau et al.

Wassan and Osman

Choi and Tcha

Version 1

Version 2

CPU

Cost

CPU

Cost

CPU

Cost

CPU

Cost

CPU

Cost

CPU

473 575 335 350 2245 2876 5833 3402 2011 – –

1494.21 603.33 1001.80 1137.01 1046.36 1812.00 1117.09 1553.72 1220.69 0.66 0

626 669 736 852 1453 1487 1681 1706 1151 – –

1499.69 608.57 999.82 1131.00 1047.74 1814.11 1108.97 1536.24 1218.27 0.46 2

276 695 893 668 2222 2847 6009 3174 2098 – –

1493.36 603.21 1001.33 1135.62 1038.60 1801.40 1106.62 1536.1 1214.53 0.15 3

4 37 6 6 103 81 299 112 81 – –

1491.86 603.21 999.82 1131.00 1038.60 1811.79 1105.439 1539.62 1215.17 0.20 6

28 29 31 19 96 62 196 109 71 – –

1491.86 603.21 999.82 1131.00 1038.60 1801.40 1105.439 1531.83 1212.90 0.01 7

101 135 137 95 312 269 839 469 295 – –

Values in boldface – identical to the best known or better (in this last case they are also in italic). a Solution cost taken from Gendreau et al. (1999). b Solution cost taken from Choi and Tcha (2007).

classical heuristic, which is rather uncommon when comparing heuristics with metaheuristics. The results in Table 4 show a trend similar to the previous ones. Both versions of the TSA give better results than the first three algorithms. On average, Version 1 is a bit faster than Taillard’s (1999) algorithm on the problems with 75 customers or more and a bit slower on the problems with 50 customers. Version 1 is a bit slower than Gendreau et al.’s (1999) algorithm and more than twice faster than Wassan and Osman’s (2002) algorithm. On the other hand, Version 2 is more than three times slower than Taillard’s (1999) algorithm, about six times slower than Gendreau et al.’s (1999) algorithm and almost twice slower than Wassan and Osman’s (2002) algorithm. In relation to Choi and Tcha’s (2007) algorithm, average costs are 0.05% lower and 0.13% higher than those of Version 1 and Version 2, respectively. On the other hand, Choi and Tcha’s (2007) computing times are more than five times higher and a bit higher than those of Version 1 and Version 2, respectively. Overall, Version 1 found 17 solutions equal or better than the best known, while the algorithm of Wassan and Osman (2002) found 15, but Version 1 found two new best solutions while the latter found none. Version 2 found 22 solutions equal or better than the best known, but three are new best solutions. By using different combinations of parameters during our experiments, the TSA was able to find all the best known solutions but two, and five new best ones for the 28 problems, i.e., 18% are new best solutions. 5. Conclusions The paper has demonstrated that the TSA is able to provide high quality solutions to the fleet size and mix vehicle routing problem in a reasonable computing time. Five new best solutions are provided for the set of test problems that have been studied by other researchers. The results presented demonstrate the good performance of the TSA when compared to the algorithms of

Taillard (1999), Gendreau et al. (1999), Wassan and Osman (2002), Choi and Tcha’s (2007) and Renaud and Boctor (2002). In addition, the TSA is a deterministic algorithm, so all the results are fully reproducible. All the first three algorithms include several random elements, so different runs of these algorithms may produce different results. Moreover, Taillard’s (1999) algorithm is not appropriate for problems of large dimensions, as the author recognises, because it includes a subroutine that requires the exact resolution of a set partitioning problem, whose computing time increases exponentially with the dimension of the problem. On the other hand, Choi and Tcha’s (2007) algorithm is mathematical programming based. This means that, depending on the stopping criterion, the running time may become exponential in n or, if the computing time is polynomial in n due to the stopping criterion, the quality of the solutions may be reduced. Choi and Tcha (2007) do not explain which stopping criterion they have used. Renaud and Boctor’s (2002) algorithm is also deterministic, but its solutions are significantly worse. Among all the factors that constitute the TSA, we should emphasize the influence of the d-neighbourhood both in quality of the solutions and in the computing time. Acknowledgements This research was partially supported by Fundacßa˜o para a Cieˆncia e Tecnologia under the Program FEDER/POCI 2010. We also wish to thank Professor E´ric Taillard and the four anonymous referees for their valuable comments. We also thank the University of Minho Editing Program for revising the text. Appendix A. New best solutions Here are presented the new best solutions found by the TSA (note that for some of these problems other authors have found better solution costs, as indicated in Section 4,

J. Branda˜o / European Journal of Operational Research 195 (2009) 716–728

but they did not publish the corresponding solutions). All the calculations have been performed with a precision of 64 bits and the total cost of the solution is presented with three decimal places in order to avoid or, at least reduce, the possibility of different solutions being taken as identical as it happens with the solution for problem 29 given here and in Gendreau et al. (1999). In the tables T is the vehicle type. Problem 8 Route

T

Cost

0-22-3-4-28-24-0 0-18-19-10-8-25-20-2-0 0-1-12-26-27-14-9-15-7-13-11-16-5-6-29-23-0 0-30-21-17-0

C D D C

390 777 840 340

Solution cost = 2347. Problem 12 Route

T

Cost

0-17-21-30-8-0 0-25-10-26-27-18-0 0-22-4-28-11-16-7-23-0 0-3-29-6-5-15-9-14-12-0 0-24-13-19-0 0-20-1-2-0

E F E E E C

672 1007 714 715 710 262

Solution cost = 4080. Problem 13 Route

T

Cost

0-2-28-22-1-43-42-41-23-16-33-6-0 0-44-3-49-24-18-50-25-9-32-40-0 0-12-39-31-10-38-11-14-35-7-0 0-17-0 0-8-19-13-15-20-37-36-47-21-48-5-29-0 0-45-34-0 0-4-0 0-26-0 0-30-27-0 0-46-0

F F F A F C B A C B

510.271 509.680 506.598 36.125 525.981 84.142 49.142 32.166 94.896 57.361

Solution cost = 2406.361. Problem 16 Route

T

Cost

0-25-14-6-27-0 0-12-5-2-0 0-49-10-9-38-11-32-0 0-13-41-40-19-42-0 0-17-37-44-15-45-33-39-30-34-50-0 0-16-21-29-35-36-20-0 0-23-24-43-7-26-48-0 0-1-22-3-28-31-8-0 0-46-47-4-18-0

B B B B C B B B B

256.841 259.329 271.682 301.167 514.687 305.871 296.012 280.224 242.326

Solution cost = 2728.139.

725

Problem 17 Route

T Cost

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

C C C B C C

291.204 251.894 233.265 148.166 222.891 273.236

C 273.420 A 40.456

Solution cost = 1734.531. Problem 18 Route

T

Cost

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

A D C C D C C B B C B D D B B A

16.000 277.748 153.782 165.733 295.133 177.289 181.410 65.700 52.125 159.956 58.416 315.813 287.466 66.324 70.626 26.125

Solution cost = 2369.646. Problem 20 Route

T

Cost

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

B A A B A A B A A A A A A B A B A A

390.071 173.042 126.834 399.670 155.783 165.932 406.754 140.036 159.949 168.443 145.566 151.535 188.404 427.863 148.609 417.559 141.126 135.409

Solution cost = 4042.586.

726

J. Branda˜o / European Journal of Operational Research 195 (2009) 716–728

Problem 28 Route

T Cost

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

E 303.092 C 117.898 C 58.629 C 153.878 F 440.220 F 595.253 B 131.832

Solution cost = 1800.802.

Problem 29 Route

T Cost

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

A 88.834 A 71.377 B 134.386 A 68.613 B 167.841 B 141.547 A 131.192 B 81.842 B 131.690 A 88.117

Solution cost = 1105.439. Problem 30 Route

T

Cost

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

C B C A C A B A A A C A A

174.396 99.782 109.680 103.969 186.886 75.015 183.317 120.284 105.635 101.689 179.571 71.581 18.627

Solution cost = 1530.43. Appendix B. Description of the TSA This section describes the operational details of the TSA, with the same parameter values that produced the final results reported in the previous tables.

Phase number – F Initial solution – SI Current solution – S Best solution (feasible or infeasible) – SB Best feasible solution – SBF Iteration counter – K Cycle counter – KC Number of restarts with SBF – KS Number of iterations since the beginning of the cycle without improving SB – KB Number of iterations since SBF was found – KBF Limit of the total number of iterations – KL Limit of KB in the current cycle – KBL Limit of KBF – KBFL Number of consecutive iterations whose solution is feasible – KF Number of consecutive iterations whose solution is infeasible – KI Limit of d – dL Step 0. Beginning of phase 1: F=1 SI is the solution calculated by one of the methods described in Section 2. h = n/2 (any rounding of the value of the parameters is done to up to the nearest integer) d=3 pffiffiffi K L ¼ 3000 n Step 1. Initialise (beginning of the cycle): Set S = SI, SB = S and SBF = S, if S is feasible; otherwise set f(SBF) = 1 Set K = KC = KS = KB = KBF = KF = KI = 0 Set KBL = 15n, KBFL = 30n Set dL = min(n  1, 2cr), where cr is the number of customers of the route of SI with more customers Set tabu list to be empty Set frequency parameters, FD = 10, FS = 5 Set penalty parameter, P = 1. Step 2. Find neighbourhood move (S0 ) from S, using all the customers as a candidate list: Set f(S0 ) = 1. Step 2.1. Perform trial single insertion moves taking into account the neighbourhood condition as defined in Section 3.1 For each admissible move, s, do: If f(s) < f(S0 ), then: (1) S0 = s and f(S0 ) = f(s); (2) If s is feasible and f(s) < f(SBF), or f(s) < f(SB) go to Step 3. Repeat until all the potential moves have been explored. Step 2.2. If K is a multiple of FD perform trial double insertion moves taking into account the neighbourhood condition as defined in Section 3.1 For each admissible move, s, do: If f(s) < f(S0 ), then:

J. Branda˜o / European Journal of Operational Research 195 (2009) 716–728

(1) S0 = s and f(S0 ) = f(s); (2) If s is feasible and f(s) < f(SBF), or f(s) < f(SB) go to Step 3. Repeat until all the potential moves have been explored. Step 2.3. If K+1 is a multiple of FS perform trial swap moves taking into account the neighbourhood condition as defined in Section 3.1 For each admissible move, s, do: If f(s) < f(S0 ), then: (1) S = s and f(S0 ) = f(s); (2) If s is feasible and f(s) < f(SBF), or f(s) < f(SB) go to Step 3. Repeat until all the potential moves have been explored. Step 3. Improve changed routes: Introduce the customers that have been moved into their respective routes using the GENI procedure, resulting in the solution S00 : Step 4. Update: Update tabu list. If S00 is feasible and f(S00 ) < f(SBF) then: (1) apply the US algorithm (P = 5 if KC = 0, else P = 6) to each route individually to get a new solution S000 ; (2) set SBF = S000 , S00 = S000 and set KS = KB = KBF = 0. If f(S00 ) < f(SB) then set SB = S00 and KB = 0. If S00 is infeasible set KI = KI + 1; else set KF = KF + 1. Set K = K + 1, KB = KB + 1, KBF = KBF + 1. If K is a multiple of 10 do: (1) if KF = 10 set P = P/2; else if KI = 10 set P = 2P; (2) if KF = 10 or KI = 10 then recalculate f(SB) with the new value of P; (3) set KF = KI = 0. Step 5. Change of the parameter values: If KB = KBL/2, then set FD = 5, FS = 10, h = max (h/2, 5). Set d = min(d + 1, dL). If KB < KBL and KC = 0 and KL - K < KBL, then set KL = KL + 15KBL, KB = KBL. If KB < KBL and KC = 1 and KL - K < 1, then set KL = KL + 15KBL, KB = KBL. If KB < KBL and KS = 0 and KL - K < 1, then set KL = KL + 15KBL, KB = KBL. Step 6. Restart: If KB = KBL, then: (1) if f(SBF) < 1 then set S = SBF, else set S = SB; (2) set KC = KC + 1, KB = KF = KI = 0, KBL = KBL + 2n, KS = 1, P = 1, FD = 10, FS = 7, h = 2h + 3, d = min(d + 2, dL); (3) recalculate f(SB) with the new value of P and empty the tabu list.

727

(Note that the changes to the parameter values, as well as emptying the tabu list at this point mean that the sequence of solutions following S is different from the sequence previously generated.) Step 7. Stopping criterion:If K P KL and KBF P KBFL then if F = 1 go to Step 8, else if F = 2 go to Step 9, else if F = 3 go to Step 10, else stop; otherwise go to Step 2. Step 8. Beginning of phase 2: F=2 SI = SBF (best feasible solution from the previous phase). pffiffiffi Set h = n/3, d = 2, K L ¼ 2000 n. Go to Step 2. Step 9. Beginning of phase 3: F=3 SI = SBF (best feasible solution from the previous phase). pffiffiffi Set h = n/2, d = 5, K L ¼ 2000 n. Go to Step 2. Step 10. Beginning of phase 4 (frequency phase): This phase has a few differences from the previous ones: in each iteration only 0.6n customers are candidates (in the even cycles, the less frequently moved, and in the odd cycles, the most frequently moved) and the neighbourhood condition is not imposed. F=4 SI = SBF (best feasible solution from the previous phase). pffiffiffi Set h = n/2, K L ¼ 2000 n. Go to Step 2. Note: At the end of Version 2 there is an additional phase, identical to phase 2, that takes as initial solution the best solution found before and with h = 3n/2 and d = 7. References Branda˜o, J.C.S., 1994. A decision support system and algorithms for the vehicle routing and scheduling problem. Doctor of Philosophy Thesis, Department of Management Science, Lancaster University. Branda˜o, J., Mercer, A., 1997. A tabu search algorithm for the multi-trip vehicle routing and scheduling problem. European Journal of Operational Research 100, 180–191. Branda˜o, J., 2004. A tabu search algorithm for the open vehicle routing problem. European Journal of Operational Research 157, 552–564. Branda˜o, J., 2006. A new tabu search algorithm for the vehicle routing problem with backhauls. European Journal of Operational Research 173, 540–555. Branda˜o, J., Eglese, R., 2006. A deterministic tabu search algorithm for the capacitated arc routing problem. Computers and Operations Research. doi:10.1016/j.cor.2006.07.007. Choi, E., Tcha, D.-W., 2007. A column generation approach to the heterogeneous fleet vehicle routing problem. Computers and Operations Research 34, 2080–2095. Clarke, G., Wright, J.W., 1964. Scheduling of vehicles from a central depot to a number of delivery points. Operations Research 12, 568–581. Dullaert, W., Janssens, G.K., So¨rensen, K., Vernimmen, B., 2002. New heuristics for the fleet size and mix vehicle routing problem with

728

J. Branda˜o / European Journal of Operational Research 195 (2009) 716–728

time windows. Journal of the Operational Research Society 53, 1232–1238. Desrochers, M., Verhoog, T.W., 1991. A new heuristic for the fleet size and mix vehicle routing problem. Computers and Operations Research 18, 263–274. Dongarra, J., 2006. Performance of various computers using standard linear equations software. Report CS-89-85, University of Tennessee. Etezadi, T., Beasley, J.E., 1983. Vehicle fleet composition. Journal of the Operational Research Society 34, 87–91. Fisher, M., Jaikumar, M., 1981. A generalized assignment heuristic for vehicle routing. Networks 11, 109–124. Gendreau, M., Hertz, A., Laporte, G., 1992. New insertion and postoptimization procedures for the traveling salesman problem. Operations Research 40 (6), 1086–1094. Gendreau, M., Hertz, A., Laporte, G., 1994. A tabu search heuristic for the vehicle routing. Management Science 40, 1276–1290. Gendreau, M., Laporte, G., Musaraganyi, C., Taillard, E., 1999. A tabu search heuristic for the heterogeneous fleet vehicle routing problem. Computers and Operations Research 26, 1153–1173. Gheysens, F., Golden, B., Assad, A., 1986. A new heuristic for determining fleet size and composition. Mathematical Programming Study 26, 233–236. Gillett, B., Miller, L., 1974. A heuristic for the vehicle dispatching problem. Operations Research 22, 340–349. Glover, F., Laguna, M., 1997. Tabu Search. Kluwer Academic Publishers. Golden, B., Assad, A., Levy, L., Gheysens, F., 1984. The fleet size and mix vehicle routing problem. Computers and Operations Research 11, 49–66. Golden, B., Wasil, E., Kelly, J., Ghao, I-M., 1998. The impact of metaheuristics on solving the vehicle routing problem: Algorithms, problem sets, and computational results. In: Grainic, T., Laporte, G. (Eds.), Fleet Management and Logistics. Kluwer Academic Publishers, Norwell, MA, pp. 33–56. Laporte, G., 1992. The vehicle routing problem: An overview of exact and approximate algorithms. European Journal of Operational Research 59 (3), 345–358. Laporte, G., Osman, I.H., 1995. Routing problems: A bibliography. Annals of Operations Research 61, 227–262. Laporte, G., Gendreau, M., Potvin, J.-Y., Semet, F., 2000. Classical and modern heuristics for the vehicle routing problem. International Transactions on Operational Research 7, 285–300.

Li, F., Golden, B., Wasil, E., 2007. A record-to-record travel algorithm for solving the heterogeneous fleet vehicle routing problem. Computers and Operations Research 34, 2734–2742. Liu, F.-H., Shen, S.-Y., 1999. Fleet size and mix vehicle routing problem with time windows. Journal of the Operational Research Society 50, 721–732. Lenstra, J.K., Rinnoy Kan, A.H.G., 1981. Complexity of vehicle routing and scheduling problems. Networks 11, 221–227. Osman, I., Salhi, S., 1996. Local search strategies for the vehicle fleet mix problem. In: Rayward-Smith, V.J., Osman, I.H., Reeves, C.R., Smith, G.D. (Eds.), Modern Heuristic Search Methods. Wiley, New York, pp. 131–153. Renaud, J., Boctor, F.F., 2002. A sweep-based algorithm for the fleet size and mix vehicle routing problem. European Journal of Operational Research 140, 618–628. Renaud, J., Boctor, F.F., Laporte, G., 1996. An improved petal heuristic for the vehicle routing problem. Journal of the Operational Research Society 47, 329–336. Rochat, Y., Semet, F., 1994. A tabu search approach for delivering pet food and flour in Switzerland. Journal of the Operational Research Society 45, 1233–1246. Rochat, Y., Taillard, E., 1995. Probabilistic diversification and intensification in local search for vehicle routing. Journal of Heuristics 1, 147–176. Semet, F., Taillard, E., 1993. Solving real-life vehicle routing problems efficiently using tabu search. In: Glover, F., Laguna, M., Taillard, E., de Werra, J.C.D. (Eds.), . In: Annals of Operations Research, vol. 41. Baltzer AG, Science Publishers, Basel – Switzerland, pp. 469–488. Salhi, S., Rand, G.K., 1993. Incorporating vehicle routing into the vehicle fleet composition problem. European Journal of Operational Research 66, 313–330. Taillard, E.D., 1999. A heuristic column generation method for the heterogeneous fleet VRP. RAIRO 33, 1–34. Tarantilis, C.D., Kiranoudis, C.T., Vassiliadis, V.S., 2004. A threshold accepting metaheuristic for the heterogeneous fixed fleet vehicle routing problem. European Journal of Operational Research 152, 148–158. Wassan, N.A., Osman, I.H., 2002. Tabu search variants for the mix fleet vehicle routing problem. Journal of the Operational Research Society 53, 768–782. Yaman, H., 2006. Formulations and valid inequalities for the heterogeneous vehicle routing problem. Mathematical Programming 106, 365–390.