Vehicle routing with time windows: An overview of exact, heuristic and metaheuristic methods

2.1 Mathematical model of the VRPTW

The VRPTW is given by a fleet of homogeneous vehicles denoted by V, a set of customers C and a directed graph G = (V, C). The graph consists of |C| + 2 vertices, where the customers are denoted 1, 2, …, n and the depot is represented by the vertex 0 (the driving-out depot) and n + 1 (the returning depot).

The VRPTW has multiple objectives in that the goal is to minimize not only the number of vehicles required, but also the total travel time, waiting time, and total travel distance incurred by the fleet of vehicles. The set of arcs denoted by A represents connections between the depot and the customers and among the customers. No arc terminates in vertex 0, and no arc originates from vertex n + 1. With each arc (i, j), where i ≠ j, we associate a cost c_ij and a time t_ij, which may include service time at customer i.

Each vehicle has a capacity q and each customer i a demand d_i. Each customer i has a time window [a_i, b_i]. A vehicle must arrive at the customer before b_i. It can arrive before a_i but the customer will not be serviced before. The depot also has a time window [a₀, b₀]. Vehicles may not leave the depot before a₀ and must be back before or at time b_n+1.

It is assumed that q, a_i, b_i, d_i, c_ij are non-negative integers, while the t_ij are assumed to be positive integers. It is assumed that the triangular inequality is satisfied for both the c_ij and the t_ij. The model contains two sets of decision variables x_ijk and s_ik. For each arc (i, j), where $i\ne j$ , $i\ne n+1$ , $j\ne 0$ , and each vehicle k we define x_ijk = 1, if and only if the optimal solution, arc (i, j) is traversed by vehicle k and equal 0, otherwise.

The decision variable s_ik is defined for each vertex i and each vehicle k and denotes the time vehicle k starts to service customer i. In case the given vehicle k does not service customer i, s_ik does not mean anything. We assume a₀ = 0 and therefore s_0k = 0, for all k. We want to design a set of minimal cost routes, one for each vehicle, such that each customer is visited exactly once, are every route originates at vertex 0 and ends at vertex n + 1, and the time windows and capacity constraints observed. We can state the VRPTW mathematically as follows:

The constraints (2.2) states that each customer is visited exactly ones, and (2.3) means that no vehicle is loaded with more than its capacity allows it to. The next three equations (2.4), (2.5), and (2.6) ensures that each vehicle leaves the depot 0, after arriving at a customer the vehicle leaves again, and finally arrives at the depot n + 1. The inequalities (2.7) states that a vehicle k cannot arrive at j before s_ik + t_ij if it is travelling from i to j. Here K is a large scalar. Finally constraints (2.8) ensure that time windows are observed, and (2.9) are the integrality constraints. Note that an unused vehicle is modeled by driving the empty route (0, n + 1).

2.2 Generalizations of the VRPTW

A number of additional constraints or properties of more complex routing problems can be modelled using the framework just developed. In this section we will briefly discuss how to allow non-identical vehicles; work with more than one depot, multi-compartment vehicles, the use of multiple time windows or soft time windows, and pick-up and delivery.

3.1 Lagrange relaxation-based methods

There are a number of papers using slightly different approaches. There is variable splitting followed by lagrange relaxation, K-tree approach followed by lagrange relaxation (Fisher et al., 1997; Holland, 1975, and in Kohl and Madsen (1997) presented shortest path with side constraints approach followed by lagrange relaxation. The relaxes of the constraints ensuring that every customer is served exactly once, that is: $$\sum _{k\in V}\sum _{j\in N}\mathit{xijk}=1\text{,}\forall i\in C$$ is relaxed and the objective function with the added penalty term then becomes $$\min {\sum }_{k\in V}{\sum }_{i\in N}{\sum }_{j\in N}(\mathit{cij}-\lambda j)\mathit{xijk}+{\sum }_{j\in V}\lambda j$$

Here λ_j is the lagrange multiplier associated with the constraint that ensures that customer j is serviced.

In Fisher et al. (1997) are presents an algorithm for solving the VRPTW optimally where the problem is formulated as a K-tree problem with degree 2K on the depot. A K-tree for a graph containing n + 1 vertices is a set of n + K edges spanning the graph. Informally, the VRPTW could be described as finding a K-tree with degree 2K on the depot, degree 2 on the customers and subject to time and capacity constraints. A K-tree with degree 2K on the depot therefore becomes equal to K routes.

The problem is defined and all constraints except the constraints ensuring that at most one arc is joining customers i and j are then lagrangian relaxed. The problem is then solved with a minimum degree-constrained K-tree problem as subproblem and the lagrange multipliers are set using the sub-gradient approach.

3.2 Column generation

In Agarwal et al. (1989) is used for the first column generation to solve the VRP problem, while (Desrosiers et al., 1984) used the column generation approach to solve the m-TSP with time windows. In Desrosiers et al. (1984), the column generation approach is used for the first time for solving the VRPTW, and more effective version of the same model with addition of valid inequalities solves more instances to optimality in Desrosiers et al. (1992).

If a linear program contains too many variables to be solved explicitly, then we can initialize the linear program with a small subset of the variables and compute a solution of this reduced linear program. Afterwards, we check if the addition of one or more variables, currently not in the linear program, might improve the linear program solution. This check can be done by the computation of the reduced costs of the variables. In this case, a variable of negative reduced cost can improve the solution.

3.3 Dynamic programming

The dynamic programming approach for VRPTW is presented for the first time in Kolen et al. (1987) and Christofides and Beasley (1984) are uses the dynamic programming paradigm to solve the VRP.

The algorithm of Kohl and Madsen (1997) use branch-and-bound to achieve optimality. Each node α in the branch-and-bound tree corresponds to three sets: F(α) which the set is of fixed feasible routes starting and finishing at the depot, P(α) which is a partially build route starting at the depot and C(α) denotes the set of customers forbidden to be next on P(α).

Branching is done by selecting a customer i that is not forbidden, that is iC(α), and that does not appear on any route, that is iF(α) ∪ P(α). Branching decisions are taken on route-customer allocations. Then two branches are generated: one in which the partially build route P(α) is extended by i and one where i is forbidden as the next customer on the route that is added to C(α). Customer i is chosen as the customer the partial route P(α) was extended with in the calculation that lead to the lower bound of node α. At each branch-and-bound node dynamic programming is used to calculate a lower bound on all feasible solution defined by F(α), P(α) and C(α).

First we discuss the case of the root node (F(α) =  $\varphi$ , C(α) =  $\varphi$ and P(α) = depot). Here we construct a directed graph with vertices v(i, q, k) for i = 0, 1, …, n; q = 0, 1, …, Q and k = 0, 1, …, m, where n is the number of customers, m is the number of vehicles and Q is the sum of all customer demands q_i. Hence, associated with each branch-and-bound node is a set of routes.

A directed path from v(0, 0, 0) to v(i, q, k) in the graph corresponds to a set of k routes with a total load of q and with different last visited customers (each one in{1, 2, …, i}). The arc lengths in the directed graph will be defined as the total length of the corresponding routes. The lower bound is then given by the minimum over k = 1, 2, …, m of the shortest paths lengths from v(0, 0, 0) to v(n, Q, k). Note that there are no constraints enforcing customers to be visited by any of the routes generated. Therefor the resulting minimum is a lower bound. Dynamically we try to extend a set of k routes with load q and last customers {1, 2, …, i} to last customers {1, 2, …, i + 1}. Here there are two possibilities:

• Customer i + 1 is not included as endpoint of any routes. This results in an arc from v(i, q, k) to v(i + 1, q, k) with length 0. Note that a customer i + 1 that is not an endpoint might still be member of one of the other routes generated by the function F(i, q) described below.

• Insert customer i + 1 as the last customer on a route of load q′. This generates an arc from v(i, q, k) to v(i + 1, q + q′, k + 1) of length F(i + 1, q′) for each possible value of q′. Generally, F(i, q) is defined as the minimum length of a feasible route with total load q and last customer i. The routes associated with a given vertex v(i, q, k) are the routes given from the computation of F(i, q) and the extension made using the arcs.

If we are at an arbitrary node in the branch-and-bound tree we distinguish between two cases: P(α) =  $\varphi$ and P(α) ≠  $\varphi$ . If P(α) =  $\varphi$ we just adjust the problem for the number of already generated routes kˆ, their load qˆ and the set of customers already used in these routes Iˆ.

In the case of P(α) ≠  $\varphi$ exactly one of the routes in the lower bound is an extension of P(α). Now, let F′(i, q) be the minimum length of such an extension (F′(i, q) is calculated in the same way as F(i, q)). As before, the problem can be reduced according to kˆ, qˆ and Iˆ. The directed graph is now extended to contain vertices v(i, q, k) and v′(i, q, k) and the following arcs:

• Arcs of length 0 from v(i, q, k) to v(i + 1, q, k) and from v′(i, q, k) to v′(i + 1, q, k). These correspond to not using the customer i + 1 in the routes.

• Arcs of length F(i + 1, q′) from v(i, q, k) to v(i + 1, q + q′, k + 1) and from v′(i, q, k) to v′(i + 1, q + q′, k + 1) for each possible value of q′.

In Kolen et al. (1987) problems up to 15 customers are solved by this method.

4.1 Route-building heuristics

The first paper on route-building heuristics for the VRPTW is (Baker and Schaffer, 1989). Their algorithm is an extension of the saving heuristic of the VRP (Clarke and Wright, 1964). The algorithm begins with all possible single-customer route (depot-i-depot). In every iteration we calculate which two routes can be combined with the maximum saving, where the saving between customers i and j are given by:

Also Landeghem (1988) has presented a heuristic based on the saving heuristic. His bi-criteria heuristic uses the time windows in (4.1) in order to get a measurement of how good a link between customers is in terms of timing.

Another heuristic described in Solomon (1987) is a time-oriented, nearest-neighbourhood heuristic. Every route in this heuristic is started by finding the unrouted customer that is closest to the depot. The closeness relation tries to take both geographical and temporal closeness of the customers into account. At every subsequent iteration the customer closest to the last customer added to the route is considered for insertion to the end of the route presently generated. When the search fails a new route is started.

Assigning customers to clusters is done by using a technique introduced in Gillett and Miller (1974) for the VRP. Here a centre of gravity is computed and the clusters are partitioned according to their polar angle. Scheduling the customers, one of the previously developed route-building heuristics is used to build a 1-route solution. Due to the time windows and/or capacity constraints some customers may now be unscheduled. In order to schedule these scheduled customers are removed and the process is repeated. The sweep heuristic typically performs better than other heuristics in cases where many customers can be assigned to each route.

In general the heuristics of Solomon (1987) and Landeghem (1988) return a solution fast. Their solution does however generally lack in quality. Most of the time the solutions are more than 10% from optimum.

A problem of building one route at a time is usually that the routes generated in the latter part of the process are of poor quality as the last unrouted customers tends to be scattered over the geographic area. In Potvin and Rousseau (1993) tries to overcome this problem of the insertion heuristics by building several routes simultaneously. The initialization of the routes is done by using the insertion heuristic of Solomon (1987). On each route the customer farthest away from the depot is selected as a seed customer. Then we proceed by computing the best feasible insertion place for each unserviced customer and insert the one with the largest difference between the best and the second best insertion place. This method is better than the Solomon heuristics but still the solutions are quite far away from optimum. In Russell (1995) elaborates further on the insertion approach.

In Antes and Derigs (1995) another approach builds upon the classical insertion idea. Here every unrouted customer requests and receives from every route in the schedule a prize for insertion, defined in a similar way as in the Solomon heuristics. Then the unrouted customers send a proposal to the route with the best offer, and each route accepts the best proposal among those customers with the fewest number of alternatives. Note that more customers can be inserted in each iteration. If a certain threshold of routes is violated a certain number of customers are removed and the process is initiated again. The results of Antes and Derigs (1995) are comparable to those presented in Potvin and Rousseau (1993). Generally building several routes in parallel results in better solutions than building the routes one by one.

4.2 Route-improving heuristics

The basis of almost every route-improving heuristic is the notion of a neighbourhood. The neighbourhood of a solution S is a set N(S) of solutions that can be generated with a single modification of S.

Checking some or all of the solutions in a neighbourhood might reveal solutions that are better with respect to objective function. This idea can be repeated from the better solution. At some point no better solution can be found an optimum has been reached. It is definitely a local optimum but it might even be global. This algorithm is called local search. The neighbourhood structures used in various VRPTW heuristics will be introduced, and then the algorithms in which they are used are described.

4.3 The neighbourhoods for the VRPTW

One of the most used improvement heuristics in routing and scheduling is the r-Opt heuristic. Here r arcs removed and replaced by r other arcs. A solution obtained using an r-Opt neighbourhood that cannot be improved further is called r-Optimal. Usually r is at most 3. Using the 3-Opt on the routes of a solution to the VRPTW problem is not without problems. For all possible 2-Opt interchanges and some of the exchanges in the 3-Opt neighbourhood, parts of the route is reversed. This may very likely lead to a violation of the time windows. In Potvin and Rousseau (1995), are presents two variants 2-Opt∗ and Or-Opt that maintain the direction of the route.

It is quite easy to see that Or-Opt’s are a subset of 3-Opt’s as we exchange three special arcs with three others. The size of the neighbourhood is although reduced from O(n³) to O(n²). Generally the size of a r-Opt neighbourhood is O(n^r). The 2-Opt∗ is exchanging one segment of one route with a segment of another route. These neighbourhood-operators are sometimes denoted crossover or simply cross.

The exchange operator swaps customers in different routes, thereby interchanging two customers simultaneously into the other routes. The k-node interchange by Christofides and Beasley (1984) is modified by several authors to take time windows into account. Sequentially each customer i is considered and the sets M₁ and M₂ are identified. M₁ is defined as the customer i and its successor j. Then the elements of M₂ are found as the two customers closest to i and j but not on the same route as i and j (found by minimizing insertion cost calculated by the Euclidean distance). A neighbourhood is then defined by removing the elements of M₁ and M₂ and inserting them in any other possible way. As this neighbourhood is quite large, only the k most promising candidates are checked.

Another neighbourhood is the λ-interchange. In Osman (1993) originally developed for the VRP. It is a generalization of the relocate operator. Here a subset of customers of size ⩽λ in one route is exchanged with a subset of size ⩽λ from another route. Typically there is also given an ordering in which the different set sizes are tested. For example using a 2-interchange scheme we would first try to move one element from one route to another, and one the other way. Then we would try the reverse situation; and then try to exchange one element from one route with one from the other. This would be written as (1, 0), (0, 1), (1, 1), (0, 2), (2, 0), (2, 1), (1, 2) and (2, 2).

Finally, a neighbourhood denoted shift-sequence is used in Schulze and Fahle (1999). A customer is moved from one route to another then checking all possible insertion positions. If an insertion is made feasible by removing another customer j, it is removed and inserted in another route. This procedure is repeated until feasibility is restored.

6.1 Simulated annealing

The name “simulated annealing” is due to the fact that conceptually it is similar to a physical process, known as annealing, where a material is heated into a liquid stat then cooled back into a recrystallized solid state. The simulated annealing was one of the first metaheuristics developed. When using simulated annealing one does not search for the best solution in the neighbourhood of the current solution. Instead one simply draws at random a solution from the neighbourhood. If the solution is better it is always accepted as a new current solution, but if the solution is worse than the present current solution it is only accepted with a certain probability. The acceptance probability is determined by a temperature which is gradually decreased. By reducing the temperature the selection becomes more and more selective in accepting new solution. The idea of simulated annealing comes from thermodynamics and metallurgy: when a metal in fusion is cooled slowly enough it tends to solidify in a structure of minimal energy.

In Chiang and Russell (1996) develop the three different simulated annealing methods: one using a modified version of the k-node interchange mechanism and the second one using the λ-interchange mechanism as proposed by Osman (1993), with λ = 1. The third algorithm borrows the concept of a tabu list from the tabu search metaheuristic. Using simulated annealing with the λ-interchange mechanism the tabu list contains moves that will not be allowed for the time being. The last two methods have faster convergence than the first, although the first yields slightly better results. The travel distances obtained by the three simulated annealing were between 7% and 11.5% from optimum.

In Thamgiah et al. (1995) a non-monotone probability function is used and are using the λ-interchange with λ = 2 scheme to define the neighbourhood. The temperature is decreased after every iteration. In case the entire neighbourhood has been explored without finding any accepting moves the temperature is increased. This is called a reset. The temperature is increased to the maximum of the temperature at which the best solution was found and half of the temperature at the last reset. After R resets without improving the best solution the algorithm terminates. The quality of the solutions obtain in Thamgiah et al. (1995) is about the same as those obtained by Chiang and Russell (1996).

6.2 Tabu search

The tabu search is one of the old metaheuristics. It was introduced by Glover (1989) and Fisher et al. (1997). At each iteration the neighbourhood of the current solution is explored and the best solution in the neighbourhood is selected as the new current solution. In order to allow the algorithm to escape from a local optimum the current solution is set to the best solution in the neighbourhood even if this solution is worse than the current solution. To prevent cycling visiting recently selected solutions is forbidden. This is implemented using a tabu list. Often, the tabu list does not contain illegal solutions, but forbidden moves. It makes sense to allow the tabu list to be overruled if this leads to an improvement of the current over-all best solution. Criteria such as this for overruling the tabu list are called aspiration criteria. The most used criteria for stopping a tabu search are a constant number of iterations without any improvement of the over-all best solution or a constant number of iteration in all.

In Solomon (1987) is used to generate the initial solution. The algorithm shifts between the two strategies. When one has not made improvement for a certain number of iterations the other strategy is used instead and vice versa. In order to minimize the number of routes the algorithm tries to move customers from routes with few customers to other routes. This is done using Or-Opt.

In Badeau et al. (1995), was first generate a series of solutions and then new solutions are composed by randomly selecting from the already generated routes. The selection is done based to the good routes. When one route is selected the remaining routes servicing customers from this route is excluded. This process is continued until all customers are serviced by a route, or the algorithm runs out of routes.

The solution is now decomposed into groups of routes. For each group a tabu search is performed using the exchange operator on segments. Furthermore segments of customers are moved around within each route. In order to force the algorithm to make a thorough exploration of the search space, frequently performed crossovers are penalized.

Another tabu search algorithm with similarities to Garcia et al. (1994) is developed in Schulze and Fahle (1999). It is based upon the local search methods discussed in Potvin and Rousseau (1995).

In Garcia et al. (1994) and Schulze and Fahle (1999) have been written about parallelization of the tabu search heuristic. Their efforts can be divided into three groups, partitioning of the neighbourhood, run parallel threads of tabu searches and decompose the problem into subproblems each solved by parallel tabu searches.

In Badeau et al. (1995) the tabu search heuristic for the VRPTW is essentially a parallelization of the tabu search. Here a combination of strategy 2 and 3 is used. After the generation of a solution at the master processor, the solution is decomposed into groups of routes. At each slave processor an independent tabu search tries to improve the over-all solution by improving the solution for the routes it received from the master processor.

Among the parallel implementations the algorithm developed in Schulze and Fahle (1999) performs slightly worse that the other parallel algorithms when working in problems with tight windows. As time windows get larger all algorithms perform equal with respect to the quality of the solution. As different parallel computers or network of computers are used it is difficult to compare the algorithms with respect to running time.

The reactive tabu search metaheuristic is applied to the parallel construction approach of Russell (1995). The tabu search route improvement procedure is invoked each time another 10% of the customers have been added to the emerging routes using the λ-interchange of Osman (1993) as the neighbourhood. If the same solution defined with the number of vehicle, the total accumulated distance, and the total accumulated travel time. Occurs to often the tabu list is increased by a constant factor. If waiting time is eliminated, the customer is fixed at its position for a number of iterations, and as in Badeau et al. (1995) too frequent switches of customers is penalized. On the other hand if no feasible solution is found by the tabu search, the size of the tabu list is decreased by constant factor. Generally it is among the best heuristics for the VRPTW. One of the conclusions in Badeau et al. (1995) is that diversification/intensification is just as important in obtaining good solutions as variable length tabu list.

6.3 The genetic algorithms

The genetic algorithms are probabilistic procedures of search, which took as a starting point the genetic evolution of a species. The principal idea is to reproduce the natural evolution of organisms, generation after generation, by respecting the phenomena’s heredity and law of survival stated by Darwin.

The first use of the genetic algorithms goes bake to the year 1950 when biologists simulated the evolution of the organisms. Later Holland (1975), Goldbeerg (1989) and Fisher (1994) adapted the genetic algorithms to the resolution of combinatorial optimization problems.

Contrary to simulated annealing and tabu search, the genetic algorithms operate in a population of solution and not only one solution (Pirlot, 1996). A population of structures, each one corresponding to a possible solution, represents the space of the solution of the problem. The new structures are generated by application of genetic operators (selection, crossing and mutation) on the potential parents chosen inside the population. The genetic algorithms are based on the principle that best parents produce best children; so the strongest members of the population have strong probability to be selected as the parents. The best parents are crossed to give the place to new children who replace the parents or the week elements of the population. The procedure is repeated until a population is obtained where all the individuals are very strong, corresponding to optimal or almost optimal solution of the problem.