Annals of Operation Research manuscript No.

(will be inserted by the editor)

A Randomized Granular Tabu Search Heuristic for the Split Delivery Vehicle Routing Problem Leonardo Berbotto · Sergio Garc´ıa · Francisco J. Nogales

Received: date / Accepted: date

Abstract The Split Delivery Vehicle Routing Problem (SDVRP) is a variant of

the classical Capacitated Vehicle Routing Problem where multiple visits to each customer are allowed. It is an NP-hard problem where exact solutions are difficult to obtain in a reasonable time. This paper shows a tabu search heuristic with granular neighborhood called Randomized Granular Tabu Search that uses a tabu search technique in a bounded neighborhood (granular) defined by the most promising arcs and introduces some new local operators in the local granular tabu search. The algorithm also uses a random selection of the move to be introduced at the current solution. In addition, the local search procedures can explore infeasible neighborhoods in terms of vehicle capacity. These two ideas help to escape from local optima. After the local search process, the algorithm solves one traveling salesman problem per route to improve the solution. Finally, a computational study shows that the proposed method improves many of the best-known solutions for the benchmark instances of the SDVRP literature. Keywords Vehicle routing problem · Split deliveries · Tabu search · Granular

neighborhood Leonardo Berbotto Universidad Carlos III de Madrid Av. Madrid 126, 28903 Getafe, Spain E-mail: [email protected] Sergio Garc´ıa Universidad Carlos III de Madrid Av. de la Universidad 30, 28911 Legan´ es, Spain E-mail: [email protected] Francisco J. Nogales Universidad Carlos III de Madrid Av. de la Universidad 30, 28911 Legan´ es, Spain E-mail: [email protected]

2

Berbotto, Garc´ıa, Nogales

1 Introduction

Routing problems are among the most studied problems in Combinatorial Optimization. The so-called Split Delivery Vehicle Routing Problem (SDVRP) belongs to this family: a set of routes with minimum total cost must be determined for a fleet of equally capacitated vehicles based at one depot serving the demand of geographically dispersed customers. It differs from the classical Capacitated Vehicle Routing Problem (CVRP) in that a customer can be visited by more than one vehicle. The interest of this variant is that splitting the deliveries makes the model more realistic and can lead to savings with respect to the transportation cost of the CVRP (see Dror and Trudeau (1989) and Archetti et al. (2006a)). Moreover, individual demands larger than the vehicle capacity are allowed and there is always a feasible solution using the minimum number of vehicles, which does not happen in the CVRP. The SDVRP was introduced in Dror and Trudeau (1989) and Dror and Trudeau (1990) and it has been applied successfully to several real situations such as the management of a fleet of trucks in a feed distribution problem (Mullaseril et al. (1997)), the problem of determining the flight schedule for helicopters to off-shore platforms for exchanging crew people employed on these platforms (Sierksma and Tijssen (1998)), and a waste collection problem with vehicles with small capacities and customers with demands greater than the vehicle capacity (Archetti and Speranza (2004)), among others. A complete survey of the SDVRP can be found in Archetti and Speranza (2011). As the classical CVRP, the SDVRP is NP-hard (Dror and Trudeau (1990)). There are very few papers studying polyhedral properties of the SDVRP: some valid inequalities are derived in Dror et al. (1994) and a new family that defines facets of the polyhedron is introduced in Belenguer et al. (2000). Therefore, exact solution methods are few and cannot solve problems too large in size. In Jin et al. (2008), it is presented a column generation approach similar to the one in Sierksma and Tijssen (1998). They improve the error gap with respect to the results of Belenguer et al. (2000). In Jin et al. (2007), a two-stage algorithm is presented: the first stage creates clusters that cover all the demands and establishes a lower bound for the optimal value, and the second stage solves a traveling salesman problem at each cluster and establishes an upper bound. An exact algorithm for the SDVRP with Time Windows (where customers must be visited within a known time range) based on a set covering formulation and a column generation technique is proposed in Gendreau et al. (2006). In Desaulniers (2010), a branch-and-price-and-cut method for the variant of the SDVRP with time windows is proposed. Recently, it was developed in Archetti et al. (2011) the first branch-and-price-and-cut method for the simple SDVRP. They propose a decomposition of the problem similar to the one of Desaulniers (2010) but with two main differences: a) visited customers are always served, and b) the algorithm used to get the solution of the subproblem is different. The authors also provide a heuristic solution using the routes generated by the solution of the subproblem. Because of the already mentioned complexity of the problem, heuristic procedures are more often found in the literature. In Dror and Trudeau (1989) and Dror and Trudeau (1990), a first local search approach is presented. Tabu search heuristics are used in Ho and Haugland (2004), Archetti et al. (2006b) and Aleman and Hill (2010). In Ho and Haugland (2004) the authors compare the saving

A Randomized Granular Tabu Search Heuristic for the SDVRP

3

between solutions with and without split delivery in a model with time windows. A tabu search with only two procedures named Order Routes and Best Neighbor is implemented in Archetti et al. (2006b). They also add an improving phase using the GENIUS algorithm developed in Gendreau et al. (1992) and a k-split cycles elimination procedure. They compare their results with Dror and Trudeau (1989) and improve almost all the considered instances. Finally, in Aleman and Hill (2010) a tabu search and a learning procedure are developed. The method is based on a set of initial solutions that are used to build a new solution with higher quality. The generated solutions are improved with a variable neighborhood descendent procedure previously introduced in Aleman et al. (2008). Other metaheuristics based techniques can be found in Campos et al. (2008), where the first algorithm based on a scatter search methodology for this problem is proposed, and in Boudia et al. (2007), where a memetic algorithm with population management is implemented by combining a genetic algorithm with local search for intensification and diversification. In Archetti et al. (2008), the tabu search procedure of Archetti et al. (2006b) is used to identify which part of the solution space has a higher probability of containing a good solution. After the identification, an integer program is run to obtain improved feasible solutions. Another hybrid technique can be found in Chen et al. (2007), where a mixed integer program is combined with a record-to-record travel algorithm to produce high quality solutions. Finally, a review on different techniques for solving the SDVRP can be found in Aleman et al. (2009). This last paper also shows a new diversification methodology. In this paper, we present a tabu search heuristic that we call Randomized Granular Tabu Search (RGTS). The local search is done in a granular neighborhood based on the idea introduced in Toth and Vigo (2003), which seems that it has not been used too much later in local search techniques for vehicle routing problems. We use the granular neighborhood to reduce the computational time necessary to explore the neighborhoods of a solution. As far as the authors know, this paper is the first that uses this idea to solve the Split Delivery Vehicle Routing Problem. We define a granular threshold for each route by introducing information of the remaining vehicle capacity of the route. Our heuristic uses several local operators to improve a solution, two of them new in the literature. Some of these local operators focus on getting a good location of a node in a route and others focus on making “good” splits of the demands, whenever splitting is necessary. Also, local improving operators are used based on known properties of the problem and on the triangle inequality. An important idea is that this algorithm uses a random selection of the moves to be introduced at the current solution based on hierarchical probabilities. This concept is based on the idea that we can reach a better solution passing through non-best current solutions for some iterations of the granular tabu search. Moreover, unlike other algorithms, RGTS allows to search in neighborhoods of infeasible solutions in terms of the vehicle capacity. The infeasibility is fixed in a later corrective phase. Once the local search is over, a TSP is solved for each route to reorder the customers. Finally, we compare the performance of our method with the best heuristics found in the literature. The results show that RGTS is very competitive with the existing methods and improves many best-known solutions of the tested instances. This paper also presents a complete summary of the best known solutions for classic sets of instances. We put together and sort all these results to obtain the fairest possible comparison with our algorithm.

4

Berbotto, Garc´ıa, Nogales

The rest of the paper is structured as follows: Section 2 introduces the mathematical formulation for the SDVRP, Section 3 describes the RGTS and computational results are exposed in Section 4. Finally, some conclusions and further research directions are given in Section 5.

2 Mathematical Formulation

In this section, the notation and formulation of the SDVRP are introduced. Let C = {1, 2, . . . , n} be the set of customers, each customer i with a positive integer demand di , and let V = {0} ∪ C , where 0 is the depot (with no demand). Each arc (i, j ) ∀i, j ∈ V has a non-negative travel cost cij (usually distances). These costs are assumed to satisfy the triangle inequality. The SDVRP is defined on a directed graph G = (V, A), where A is the set of arcs {(i, j ) / i, j ∈ V, i 6= j}. A fleet of K identical vehicles, each with capacity Q ∈ Z+ , is available at the depot. We assume that K is equal to Kmin , the minimum number of vehicles needed to serve all the customers, which is determined as Kmin = dd(C )/Qe, where d(C ) is the sum of all the demands. Defining variables xijk as binary variables that take value one if vehicle k travels directly from i to j (and zero otherwise) and yik ≥ 0 as the demand of customer i ∈ C delivered by vehicle k = 1, 2, . . . , K , the vehicle flow formulation for the SDVRP (Archetti and Speranza (2008)) is: Min.

K X X

cij xijk

(i,j )∈A k=1

s.t.

K XX i∈I k=1 K XX

xijk ≥ 1, x0jk = K,

j∈I k=1 X

xipk −

i∈I X X

j ∈ C, (1)

X

(2)

xpjk = 0, p ∈ I, k = 1, 2, . . . , K, (3)

j∈I

xijk ≤ |S| − 1,

S ⊆ C, k = 1, 2, . . . , K, (4)

i∈S j∈S

yik ≤ di

X

xijk ,

i ∈ C, k = 1, 2, . . . , K, (5)

j∈I

X

yik ≤ Q,

k = 1, 2, . . . , K, (6)

yik = di ,

i ∈ C, (7)

i∈C K X k=1

xijk ∈ {0, 1}, yik ≥ 0,

(i, j ) ∈ A, k = 1, 2, . . . , K, i ∈ C, k = 1, 2, . . . , K.

The objective function is the total transportation cost of the K routes. Constraints (1) state that every customer must be visited by at least one vehicle while constraints (2) impose that exactly K routes are performed. Constraints (3)

A Randomized Granular Tabu Search Heuristic for the SDVRP

5

indicate that, if vehicle k visits node p, then it must leave it (flow conservation constraints) and (4) are the classical subtour elimination constraints. Constraints (5) guarantee that the customer can only be served if the vehicle visits it. Constraints (6) limit to Q the maximum load of each vehicle while constraints (7) ensure that the entire demand of each node is satisfied. Note that it is showed in Archetti et al. (2006a) that if the demands di and the capacity Q are integer values, then there exists always an optimal solution to the previous formulation where variables yik have integer values.

3 Solution method

In this section we present a tabu search heuristic with granular neighborhood that we call Randomized Granular Tabu Search (RGTS). First, some concepts about tabu search and local granular search are introduced. The neighborhood of a solution is a set of other solutions that can be obtained from a simple modification of the current solution. This modification (i.e., the transition from a solution to a new one) is called a move. In a simple local search, at each iteration, the best solution of the neighborhood that improves the current solution is considered. A weakness of this methodology is that, once the solution is close to a local optimum, it may be difficult to escape from it. The tabu search technique deals with this problem. A basic tabu search is based on operators designed to cross boundaries of local optimality or feasibility by guiding a local search procedure to explore the solution space beyond local optimality (Glover and Laguna (2001)). The main component of the tabu search is its use of adaptive memory: to escape from a local optimum and cycles, a tabu list or short term memory is used. This tabu list introduces the main attributes of a move and forbids new moves that imply at least one attribute from the list. Intensification and diversification components are also important in a tabu search. During the intensification stage, the algorithm searches on a neighborhood of the best solution by returning to attractive regions to search more thoroughly. On the diversification stage, the search is made on unvisited regions. After a given number of non-improving moves, the tabu search is stopped. Although this is an effective tool for avoiding local optima, the time required to explore all possible neighborhoods may be very large, especially in large instances. For this reason, a restricted neighborhood called granular neighborhood (Toth and Vigo (2003)) can be obtained from the standard one by removing the arcs that cannot belong to high quality feasible solutions, defining so an elite neighborhood. Because a granular neighborhood can be examined in much less time than the whole original neighborhood, the time required to reach a high quality solution is often much smaller. Considering the CVRP, Toth and Vigo (2003) observe that long arcs (i.e., expensive) have a low probability of being part of high quality solutions. In order to define what a “short” arc is, a granular threshold value g is defined as: g=β

z0

n+K

,

where β is a positive parameter and z 0 is the value of an initial CVRP solution obtained by a fast algorithm. The denominator represents the total number of arcs used in a solution with K vehicles. Note that the fractional term represents the

6

Berbotto, Garc´ıa, Nogales

average length of an arc in the initial solution. The granular neighborhood of a solution is defined on the arcs (i, j ) ∈ A such that cij ≤ g.

3.1 The Randomized Granular Tabu Search algorithm Our Randomized Granular Tabu Search heuristic has four phases: I) II) III) IV)

Initial solution. Granular tabu search. Capacity correction. Solution improvement.

3.1.1 Phase I: Initial solution

The initial solution is obtained with the well known Clarke and Wright savings algorithm (Clarke and Wright (1964)) based on the savings from merging two simple routes (0 − i − 0) and (0 − j − 0) into one route (0 − i − j − 0). The saving for nodes i, j ∈ C is defined as sij = c0i + c0j − cij . A decreasing ordered list of these savings is built by storing at each row the information of the nodes i and j and the saving sij . We use a parallel version of the Clarke and Wright algorithm for split delivery demands described in Aleman and Hill (2010), but slightly modified because we fix the number of vehicles to K = Kmin . We also evaluated the performance of RGTS using two different versions of the Clarke and Wright savings algorithm: a sequential version where the routes are built one by one, and a modified parallel version where splits at interior points are allowed (contrary to the parallel algorithm of Aleman and Hill (2010), which only allows splits at the extremes of the routes). We observed that the parallel version presented by Aleman and Hill (2010) gave better initial solutions for almost all the tested instances. Besides, in those cases where either the sequential algorithm or the modified parallel version improved the solution provided by Aleman and Hill (2010), the final solution of RGTS was not too much better. Therefore, for the sake of simplicity, we only use the parallel Clarke and Wright to obtain the initial solution. 3.1.2 Phase II: Granular tabu search

In this phase, seven different local search operators are considered and a tabu list is built. Each operator searches for an improving move in a granular neighborhood of the current solution and the best move (maximum saving) from among some promising moves randomly selected is used to obtain a new current solution. The move introduced at the current solution is stored in the tabu list and any other move that involves at least one attribute of the move in the list is tabu for the next t iterations. Thus, t is the maximum size of the tabu list. For the local search procedures we use a modified version of the granular threshold g that considers particular information regarding the remaining capacities of the vehicles. This granular threshold gk is defined for each vehicle k as: z C C gk = β 1 + k = gk0 + gk0 k , n+K +s C C

A Randomized Granular Tabu Search Heuristic for the SDVRP

7

where z is the value of the current solution, s is the number of splits of the current solution (which must be added to n + K to obtain the total number of arcs used in this solution), Ck is the remaining capacity of vehicle k, C represents the average of the remaining capacities of the routes that still have some unused capacity and gk0 = β z/(n + K + s). This new term in the definition of gk allows to create routes in a neighborhood of solutions that are not in the normal granular space but that are attractive because they have a relevant remaining capacity. Note that the granular threshold is defined for each route because it depends on the remaining capacity of each route k and it is calculated at each iteration. In the explanations of the local moves we use the following notation. Denoting with rk the route performed by vehicle k, we say that a node i ∈ rk stays in the 0 granular neighborhood of j ∈ rv if cij ≤ gv and we denote it as i ∈ gv (j ). Finally, yik 1 and yik are the demands of customer i served by vehicle k before introducing and after having introduced a move, respectively. The saving cost for each move is weighted by a factor that represents the remaining capacity of the route(s) where the node(s) will be inserted. Let D be the set of arcs to be deleted, F the set of arcs to be added, R the set of routes that receive new demand and H the set of nodes involved in the move. The weighted saving cost for move p is: ! X X X X ca − ca yj , sp = Ck − a∈D

a∈F

k∈R

j∈H

where Ck represents the remaining capacity of route k and yj represents the demand of node j to be relocated. As in the granular definition, the saving costs also favor those moves implying routes with a larger remaining capacity. The insertion cost of adding some nodes to a route is defined as the difference between the added arcs and the erased arcs. Similarly, the erasing cost of deleting some nodes from a route is the difference between the erased arcs and the added arcs. The operators used in our tabu search are the following: – Relocate Node (Relocate) : for customer nodes i ∈ rk and j ∈ rv , if either i ∈ gv (j ) or i ∈ gv (j + 1), and if the remaining capacity of rv allows the move (i.e., after the move, rv would deliver at most the allowed capacity limit), then move node i from rk to rv and locate it after node j on rv . Demands are set as 1 1 0 following: yik = 0 and yiv = yik (see Figure 1). – Exchange Node (Exchange) : for customer nodes i, t ∈ rk and j, p ∈ rv , if either i ∈ gv (p) or i ∈ gv (p +1), and either j ∈ gk (t) or j ∈ gk (t +1), and the remaining capacities of rk and rv allow the move, then move node i from rk to rv , locating it at after node p, and move j from rv to rk , locating it after node t. Demands 1 1 0 1 1 0 are set as following: yik = 0, yiv = yik and yjv = 0, yjk = yjv (see Figure 2). – 2-Opt (2Opt) : for customer nodes i ∈ rk and j ∈ rv , if i ∈ gv (j + 1) and j ∈ gk (i+1), switch the initial part of rk (rv ) until node i (node j ) and reconnect it with the part that starts at node j + 1 on rv (i + 1 on rk ). Demands are set 1 0 1 0 as following: yzv = yzk (yzk = yzv ) for every node z of rk (rv ) up to node i (j ). This move is done if the remaining capacities of vehicles k and v allow it (see

Figure 3). – Exchange Split (ExchSplit) : for a customer node i visited by routes rk and rv (i.e., i ∈ rk and i ∈ rv ) and customer nodes z, z + 1 ∈ rv and j ∈ rk , if either

8

Berbotto, Garc´ıa, Nogales 0 0 j ∈ gv (z ) or j ∈ gv (z + 1), and if yjk − yiv > 0, then erase i from rv and generate a new split at node j by inserting it after node z + 1 on rv . Demands are set as 1 1 0 0 1 0 0 1 0 following: yiv = 0, yik = yik + yiv , yjk = yjk − yiv and yjv = yiv (see Figure 4). – Delete Split and New Split (DelSplitNew) : for a customer node i visited by routes rk and rv and customer nodes j ∈ rk and z, z + 1 ∈ rt , if either j ∈ gt (z ) 0 0 or j ∈ gt (z + 1), yjk − yiv > 0 and the remaining capacity of rt allows the move, then erase i from rv and generate a new split on node j between route k and route t, locating it after node z on rt . Demands are set as following: 1 1 0 0 1 0 0 1 0 yiv = 0, yik = yik + yiv , yjk = yjk − yiv and yjt = yiv (see Figure 5). – Delete Split (DelSplit) : given a customer node i that is visited by more than one vehicle, let B be the set of routes that contain i in their paths and let B 0 be the set of routes without i in their paths. Let r∗ be the route in B 0 with the minimum insertion cost of node i. If for some node p ∈ r∗ , either i ∈ gr∗ (p) or i ∈ gr∗ (p + 1), and the remaining vehicle capacity of r∗ allows the move, then delete i from all routes in B and introduce the node in r∗ after node p. 1 1 Demands are set as following: yik = 0 ∀k ∈ B and yir ∗ = di (see Figure 6). – Delete Worst Split (DelWorstSplit) : let B be the set of routes that contain node i in their paths and let rk ∈ B be the route with the greatest saving cost of erasing i. Sort B by decreasing saving cost, delete i from the first route in B , 0 that is rk , and split yik among the other routes in B starting by the last one of

the list (the one with the lowest saving cost). When the vehicle is full, follow 0 with the next one and so on until having assigned the full amount yik (see Figure 7).

(a) Routes before Relocate.

(b) Routes after Relocate.

Fig. 1: Relocate operator.

The operators Relocate, Exchange and 2Opt are the classical ones developed for the traveling salesman problem and, jointly with with ExchSplit, are applied in Ho and Haugland (2004). DelWorstSplit is presented in Archetti et al. (2006b), but they use the insertion criterion in the inverse sense: starting by the route with the second greater saving cost. We use a different criterion because we consider that a route with a lower saving cost from erasing i has a smaller chance to delete this node in next moves and, thus, its demand must not be relocated. For example, think of a route k that visits the node i located on the segment that connect nodes j and z (that means, cji + ciz = cjz ) and assume that vehicle k serves a large demand yik . In this case, the saving cost from erasing i is equal to zero,

A Randomized Granular Tabu Search Heuristic for the SDVRP

(a) Routes before Exchange.

9

(b) Routes after Exchange.

Fig. 2: Exchange operator.

(a) Routes before 2Opt.

(b) Routes after 2Opt.

Fig. 3: 2Opt operator.

(a) Routes before ExchSplit.

(b) Routes after ExchSplit.

Fig. 4: ExchSplit operator.

the lowest possible because there is no extra cost for visiting i when the vehicle goes from j to z . This means that i is served by vehicle k with a cost of zero. On the other hand, if i presents a large saving cost from being erased from route k, it means that it is a promising node to be moved to other route in some future iteration and relocating yik might be complicated. The operators DelNewSplit and DelSplit are new and they focus on the relocation and deletion of splits. We introduce these operators because in a SDVRP with limited vehicle fleet size it is more important to optimize the split of the demands than when the fleet is unlimited. At every iteration, each operator searches for a move to improve the current solution based on the largest saving cost. For each operator, if there is more than one move with positive saving cost, they are sorted in a list in descendent saving

10

Berbotto, Garc´ıa, Nogales

(a) Routes before DelSplitNew.

(b) Routes after DelSplitNew.

Fig. 5: DelSplitNew operator.

(a) Routes before DelSplit.

(b) Routes after DelSplit.

Fig. 6: DelSplit operator, B={rk , rv }, B’={rt }.

(a) Routes before DelWorstSplit.

(b) Routes after DelWorstSplit.

Fig. 7: DelWorstSplit operator.

cost. Next, the algorithm uses randomness in the selection of the candidate moves for each operator. The idea is that, in order to obtain a better final solution, it is not necessary to pass through the best current solutions at each iteration. Thus, for each operator the algorithm selects randomly one move among the m best candidate moves by assigning higher probabilities to the moves with larger saving costs. The random selection of the candidate moves at each operator is defined as follows: 1. Fix the maximum candidate list size to m. 2. Obtain q ≤ m candidates.

A Randomized Granular Tabu Search Heuristic for the SDVRP

11

3. Assign the following probabilities to each candidate: – p1 = q/qˆ for the candidate with the largest saving cost, – p2 = (q − 1)/qˆ for the candidate with the second largest saving cost, – ... – pq = 1/qˆ for the candidate with the q -th largest saving cost, P where qˆ = qi=1 i. For example, if m = q = 5 for all the operators, the moves for each operator have the following probabilities: p1 = 5/15, p2 = 4/15, p3 = 3/15, p4 = 2/15, and p5 = 1/15. With these probabilities, one move is randomly selected for each operator. Once a candidate move has been selected for each operator, the one with the largest saving cost is implemented in the current solution. This move is stored in the tabu list and the reverse move is forbidden for t iterations. Table 1 shows an example of a possible iteration. Each operator presents a maximum of m = q = 5 candidate moves with positive saving costs and, with the probabilities explained above, the candidates in bold letter are selected. The best move among these selected candidates, ExchSplit, is implemented. Table 1: Tabu Phase with m=5 Moves Relocate Exchange 2opt ExchSplit DelSplitNew DelSplit DelWorstSplit ∗ By

saving cost;

10 7∗∗ 9 10 6 9 8 ∗∗ by

Candidates∗ 8 6∗∗ 4 6 3 2 8 7∗∗ 6 9∗∗ 3 2 5 2∗∗ 1 8∗∗ 7 4 7 3 2∗∗

I 2 1 3 1 0 2 0

⇒

Selected moves∗ 6 7 7 9 2 8 2

II

We choose

⇒

ExchSplit

random selection.

It must be noted that an alternative strategy consisting in selecting the best move for each process and then implementing randomly one of them was also tested but it was discarded because it performed worse. The intensification and diversification criterion are introduced as follows: if an improving move is found at an iteration, then the value of β in gk is reduced by a certain value φ% to intensify the search in this neighborhood. Otherwise, if there is no candidate move with positive saving cost at an iteration, the algorithm allows to introduce a non-positive move and a diversification strategy is applied: the value of β is increased by φ% for the next iteration to explore a new neighborhood. We do not fix a maximum value for β . Thus, at each iteration, either a diversification or an intensification strategy is applied. Finally, if a move is tabu but it leads to an improvement in the current best solution (incumbent solution), it is allowed to be a candidate move. Another important idea is that we allow the algorithm to explore neighborhoods of infeasible solutions in terms of vehicle capacity. The infeasibility in capacity is controlled by a parameter π such that the total demand served by a route cannot exceed the threshold (1 + π )Q, π ≥ 0. Thus, the tabu phase relaxes constraints (6) of the vehicle flow formulation to X yik ≤ (1 + π )Q, k = 1, . . . , K. i∈C

12

Berbotto, Garc´ıa, Nogales

The granular tabu phase stops when there are T consecutive iterations with no positive moves. Besides, this current solution may be infeasible in terms of the capacities of the vehicles, in which case this infeasibility is fixed in the next phase of the algorithm. 3.1.3 Phase III: Capacity correction

This phase is carried out only if there is some vehicle whose total served demand exceeds the capacity limit Q. It presents a set of operators combining up to three routes and up to three nodes that eliminate the excess of delivered product by combining different simple ways of making new splits or/and relocations. Given a route rk with excess of capacity and routes t and v with enough remaining capacity, the following operators are used to move product from the route with excess: 1. Two routes and one node: (a) For node i ∈ rk relocate i on rt . (b) For node i ∈ rk , rt relocate some demand of i from rk to rt . (c) For node i ∈ rk , introduce i in rt and split its demand. 2. Three routes and two nodes: (a) For nodes i, j ∈ rk , relocate them on rt and rv , respectively. (b) For nodes i, j ∈ rk , relocate i on rt , introduce j in rv and split yjk between rk and rv . If there is no operator among the previous ones that eliminates the whole excess of demand of the route, then we consider more complex operators that combine two or three of the ones in the previous list (considering a maximum of up to three nodes). If this new operator allows to remove the whole excess of demand, then it is implemented. The first of these operators (simple or complex) that eliminates totally the excess of delivered product is applied. If the excess of capacity of a route cannot be totally eliminated either with a simplex or a complex operator because there is no route with enough remaining capacity, then partial eliminations are allowed: the excess is reduced a ρ% (20% in our algorithm) and the procedure is restarted for the remaining excess. This guarantees that we will end with a solution where original capacity limits are not violated. 3.1.4 Phase IV: Solution improvement

After the capacity correction, the algorithm improves the obtained solution with three operators: 2-Split Cycles, Node-Position and Subtours-Elimination. The 2-Splits Cycle operator eliminates solutions where two routes share two split nodes. The motivation for this move is that it is shown in Dror and Trudeau (1989) that, if the costs satisfy the triangle inequality, then there is an optimal solution where any two routes share at most one customer. Subtour-Elimination is a new operator that deletes subtours in a given solution (based on constraints (4)) by erasing multiple visits of a vehicle to the same customer. In our method it is possible to obtain a solution with subtours in the previous phase because the local search is made in a bounded neighborhood. Thus, a relocation of node i in a route that already visits this node is something that may happen. Node-Position is another new operator where a node is moved to a different position if this leads to a reduction in the

A Randomized Granular Tabu Search Heuristic for the SDVRP

13

route cost. It tries to move a node only up to three positions backward or forward in the route. Next, the algorithm applies a granular local search with the same procedures of the tabu phase but setting π to 1 and selecting the best move at each step (i.e., the candidate list size is m = 1). Finally, a TSP is solved for every route looking for possible savings by reordering the nodes in the route. We use the standard TSP formulation with the MTZ subtour elimination constraints (Miller et al. (1960)). Only routes with at least four nodes are considered because for the routes with three nodes this is equivalent to the application of the procedure Node-position. We impose a time limit of 20 seconds per node for solving these problems. Because of the presence of randomness in the selection of the move at each local operator, the granular tabu search and capacity correction phases are executed M times. The tabu list is rebooted at each of these executions. To sum up, the structure of the algorithm is the following: 1. Phase I: Obtain an initial solution x0 with the Clarke and Wright algorithm adapted from Aleman and Hill (2010). Define BestSolution=IncumbentSolution:=x0 . 2. Repeat M times: (a) Define InitialSolution:=BestSolution. (b) Run Phase II with InitialSolution as starting solution. (c) Run Phase III and let BestSolution the final solution of this phase. (d) If value(BestSolution)

(will be inserted by the editor)

A Randomized Granular Tabu Search Heuristic for the Split Delivery Vehicle Routing Problem Leonardo Berbotto · Sergio Garc´ıa · Francisco J. Nogales

Received: date / Accepted: date

Abstract The Split Delivery Vehicle Routing Problem (SDVRP) is a variant of

the classical Capacitated Vehicle Routing Problem where multiple visits to each customer are allowed. It is an NP-hard problem where exact solutions are difficult to obtain in a reasonable time. This paper shows a tabu search heuristic with granular neighborhood called Randomized Granular Tabu Search that uses a tabu search technique in a bounded neighborhood (granular) defined by the most promising arcs and introduces some new local operators in the local granular tabu search. The algorithm also uses a random selection of the move to be introduced at the current solution. In addition, the local search procedures can explore infeasible neighborhoods in terms of vehicle capacity. These two ideas help to escape from local optima. After the local search process, the algorithm solves one traveling salesman problem per route to improve the solution. Finally, a computational study shows that the proposed method improves many of the best-known solutions for the benchmark instances of the SDVRP literature. Keywords Vehicle routing problem · Split deliveries · Tabu search · Granular

neighborhood Leonardo Berbotto Universidad Carlos III de Madrid Av. Madrid 126, 28903 Getafe, Spain E-mail: [email protected] Sergio Garc´ıa Universidad Carlos III de Madrid Av. de la Universidad 30, 28911 Legan´ es, Spain E-mail: [email protected] Francisco J. Nogales Universidad Carlos III de Madrid Av. de la Universidad 30, 28911 Legan´ es, Spain E-mail: [email protected]

2

Berbotto, Garc´ıa, Nogales

1 Introduction

Routing problems are among the most studied problems in Combinatorial Optimization. The so-called Split Delivery Vehicle Routing Problem (SDVRP) belongs to this family: a set of routes with minimum total cost must be determined for a fleet of equally capacitated vehicles based at one depot serving the demand of geographically dispersed customers. It differs from the classical Capacitated Vehicle Routing Problem (CVRP) in that a customer can be visited by more than one vehicle. The interest of this variant is that splitting the deliveries makes the model more realistic and can lead to savings with respect to the transportation cost of the CVRP (see Dror and Trudeau (1989) and Archetti et al. (2006a)). Moreover, individual demands larger than the vehicle capacity are allowed and there is always a feasible solution using the minimum number of vehicles, which does not happen in the CVRP. The SDVRP was introduced in Dror and Trudeau (1989) and Dror and Trudeau (1990) and it has been applied successfully to several real situations such as the management of a fleet of trucks in a feed distribution problem (Mullaseril et al. (1997)), the problem of determining the flight schedule for helicopters to off-shore platforms for exchanging crew people employed on these platforms (Sierksma and Tijssen (1998)), and a waste collection problem with vehicles with small capacities and customers with demands greater than the vehicle capacity (Archetti and Speranza (2004)), among others. A complete survey of the SDVRP can be found in Archetti and Speranza (2011). As the classical CVRP, the SDVRP is NP-hard (Dror and Trudeau (1990)). There are very few papers studying polyhedral properties of the SDVRP: some valid inequalities are derived in Dror et al. (1994) and a new family that defines facets of the polyhedron is introduced in Belenguer et al. (2000). Therefore, exact solution methods are few and cannot solve problems too large in size. In Jin et al. (2008), it is presented a column generation approach similar to the one in Sierksma and Tijssen (1998). They improve the error gap with respect to the results of Belenguer et al. (2000). In Jin et al. (2007), a two-stage algorithm is presented: the first stage creates clusters that cover all the demands and establishes a lower bound for the optimal value, and the second stage solves a traveling salesman problem at each cluster and establishes an upper bound. An exact algorithm for the SDVRP with Time Windows (where customers must be visited within a known time range) based on a set covering formulation and a column generation technique is proposed in Gendreau et al. (2006). In Desaulniers (2010), a branch-and-price-and-cut method for the variant of the SDVRP with time windows is proposed. Recently, it was developed in Archetti et al. (2011) the first branch-and-price-and-cut method for the simple SDVRP. They propose a decomposition of the problem similar to the one of Desaulniers (2010) but with two main differences: a) visited customers are always served, and b) the algorithm used to get the solution of the subproblem is different. The authors also provide a heuristic solution using the routes generated by the solution of the subproblem. Because of the already mentioned complexity of the problem, heuristic procedures are more often found in the literature. In Dror and Trudeau (1989) and Dror and Trudeau (1990), a first local search approach is presented. Tabu search heuristics are used in Ho and Haugland (2004), Archetti et al. (2006b) and Aleman and Hill (2010). In Ho and Haugland (2004) the authors compare the saving

A Randomized Granular Tabu Search Heuristic for the SDVRP

3

between solutions with and without split delivery in a model with time windows. A tabu search with only two procedures named Order Routes and Best Neighbor is implemented in Archetti et al. (2006b). They also add an improving phase using the GENIUS algorithm developed in Gendreau et al. (1992) and a k-split cycles elimination procedure. They compare their results with Dror and Trudeau (1989) and improve almost all the considered instances. Finally, in Aleman and Hill (2010) a tabu search and a learning procedure are developed. The method is based on a set of initial solutions that are used to build a new solution with higher quality. The generated solutions are improved with a variable neighborhood descendent procedure previously introduced in Aleman et al. (2008). Other metaheuristics based techniques can be found in Campos et al. (2008), where the first algorithm based on a scatter search methodology for this problem is proposed, and in Boudia et al. (2007), where a memetic algorithm with population management is implemented by combining a genetic algorithm with local search for intensification and diversification. In Archetti et al. (2008), the tabu search procedure of Archetti et al. (2006b) is used to identify which part of the solution space has a higher probability of containing a good solution. After the identification, an integer program is run to obtain improved feasible solutions. Another hybrid technique can be found in Chen et al. (2007), where a mixed integer program is combined with a record-to-record travel algorithm to produce high quality solutions. Finally, a review on different techniques for solving the SDVRP can be found in Aleman et al. (2009). This last paper also shows a new diversification methodology. In this paper, we present a tabu search heuristic that we call Randomized Granular Tabu Search (RGTS). The local search is done in a granular neighborhood based on the idea introduced in Toth and Vigo (2003), which seems that it has not been used too much later in local search techniques for vehicle routing problems. We use the granular neighborhood to reduce the computational time necessary to explore the neighborhoods of a solution. As far as the authors know, this paper is the first that uses this idea to solve the Split Delivery Vehicle Routing Problem. We define a granular threshold for each route by introducing information of the remaining vehicle capacity of the route. Our heuristic uses several local operators to improve a solution, two of them new in the literature. Some of these local operators focus on getting a good location of a node in a route and others focus on making “good” splits of the demands, whenever splitting is necessary. Also, local improving operators are used based on known properties of the problem and on the triangle inequality. An important idea is that this algorithm uses a random selection of the moves to be introduced at the current solution based on hierarchical probabilities. This concept is based on the idea that we can reach a better solution passing through non-best current solutions for some iterations of the granular tabu search. Moreover, unlike other algorithms, RGTS allows to search in neighborhoods of infeasible solutions in terms of the vehicle capacity. The infeasibility is fixed in a later corrective phase. Once the local search is over, a TSP is solved for each route to reorder the customers. Finally, we compare the performance of our method with the best heuristics found in the literature. The results show that RGTS is very competitive with the existing methods and improves many best-known solutions of the tested instances. This paper also presents a complete summary of the best known solutions for classic sets of instances. We put together and sort all these results to obtain the fairest possible comparison with our algorithm.

4

Berbotto, Garc´ıa, Nogales

The rest of the paper is structured as follows: Section 2 introduces the mathematical formulation for the SDVRP, Section 3 describes the RGTS and computational results are exposed in Section 4. Finally, some conclusions and further research directions are given in Section 5.

2 Mathematical Formulation

In this section, the notation and formulation of the SDVRP are introduced. Let C = {1, 2, . . . , n} be the set of customers, each customer i with a positive integer demand di , and let V = {0} ∪ C , where 0 is the depot (with no demand). Each arc (i, j ) ∀i, j ∈ V has a non-negative travel cost cij (usually distances). These costs are assumed to satisfy the triangle inequality. The SDVRP is defined on a directed graph G = (V, A), where A is the set of arcs {(i, j ) / i, j ∈ V, i 6= j}. A fleet of K identical vehicles, each with capacity Q ∈ Z+ , is available at the depot. We assume that K is equal to Kmin , the minimum number of vehicles needed to serve all the customers, which is determined as Kmin = dd(C )/Qe, where d(C ) is the sum of all the demands. Defining variables xijk as binary variables that take value one if vehicle k travels directly from i to j (and zero otherwise) and yik ≥ 0 as the demand of customer i ∈ C delivered by vehicle k = 1, 2, . . . , K , the vehicle flow formulation for the SDVRP (Archetti and Speranza (2008)) is: Min.

K X X

cij xijk

(i,j )∈A k=1

s.t.

K XX i∈I k=1 K XX

xijk ≥ 1, x0jk = K,

j∈I k=1 X

xipk −

i∈I X X

j ∈ C, (1)

X

(2)

xpjk = 0, p ∈ I, k = 1, 2, . . . , K, (3)

j∈I

xijk ≤ |S| − 1,

S ⊆ C, k = 1, 2, . . . , K, (4)

i∈S j∈S

yik ≤ di

X

xijk ,

i ∈ C, k = 1, 2, . . . , K, (5)

j∈I

X

yik ≤ Q,

k = 1, 2, . . . , K, (6)

yik = di ,

i ∈ C, (7)

i∈C K X k=1

xijk ∈ {0, 1}, yik ≥ 0,

(i, j ) ∈ A, k = 1, 2, . . . , K, i ∈ C, k = 1, 2, . . . , K.

The objective function is the total transportation cost of the K routes. Constraints (1) state that every customer must be visited by at least one vehicle while constraints (2) impose that exactly K routes are performed. Constraints (3)

A Randomized Granular Tabu Search Heuristic for the SDVRP

5

indicate that, if vehicle k visits node p, then it must leave it (flow conservation constraints) and (4) are the classical subtour elimination constraints. Constraints (5) guarantee that the customer can only be served if the vehicle visits it. Constraints (6) limit to Q the maximum load of each vehicle while constraints (7) ensure that the entire demand of each node is satisfied. Note that it is showed in Archetti et al. (2006a) that if the demands di and the capacity Q are integer values, then there exists always an optimal solution to the previous formulation where variables yik have integer values.

3 Solution method

In this section we present a tabu search heuristic with granular neighborhood that we call Randomized Granular Tabu Search (RGTS). First, some concepts about tabu search and local granular search are introduced. The neighborhood of a solution is a set of other solutions that can be obtained from a simple modification of the current solution. This modification (i.e., the transition from a solution to a new one) is called a move. In a simple local search, at each iteration, the best solution of the neighborhood that improves the current solution is considered. A weakness of this methodology is that, once the solution is close to a local optimum, it may be difficult to escape from it. The tabu search technique deals with this problem. A basic tabu search is based on operators designed to cross boundaries of local optimality or feasibility by guiding a local search procedure to explore the solution space beyond local optimality (Glover and Laguna (2001)). The main component of the tabu search is its use of adaptive memory: to escape from a local optimum and cycles, a tabu list or short term memory is used. This tabu list introduces the main attributes of a move and forbids new moves that imply at least one attribute from the list. Intensification and diversification components are also important in a tabu search. During the intensification stage, the algorithm searches on a neighborhood of the best solution by returning to attractive regions to search more thoroughly. On the diversification stage, the search is made on unvisited regions. After a given number of non-improving moves, the tabu search is stopped. Although this is an effective tool for avoiding local optima, the time required to explore all possible neighborhoods may be very large, especially in large instances. For this reason, a restricted neighborhood called granular neighborhood (Toth and Vigo (2003)) can be obtained from the standard one by removing the arcs that cannot belong to high quality feasible solutions, defining so an elite neighborhood. Because a granular neighborhood can be examined in much less time than the whole original neighborhood, the time required to reach a high quality solution is often much smaller. Considering the CVRP, Toth and Vigo (2003) observe that long arcs (i.e., expensive) have a low probability of being part of high quality solutions. In order to define what a “short” arc is, a granular threshold value g is defined as: g=β

z0

n+K

,

where β is a positive parameter and z 0 is the value of an initial CVRP solution obtained by a fast algorithm. The denominator represents the total number of arcs used in a solution with K vehicles. Note that the fractional term represents the

6

Berbotto, Garc´ıa, Nogales

average length of an arc in the initial solution. The granular neighborhood of a solution is defined on the arcs (i, j ) ∈ A such that cij ≤ g.

3.1 The Randomized Granular Tabu Search algorithm Our Randomized Granular Tabu Search heuristic has four phases: I) II) III) IV)

Initial solution. Granular tabu search. Capacity correction. Solution improvement.

3.1.1 Phase I: Initial solution

The initial solution is obtained with the well known Clarke and Wright savings algorithm (Clarke and Wright (1964)) based on the savings from merging two simple routes (0 − i − 0) and (0 − j − 0) into one route (0 − i − j − 0). The saving for nodes i, j ∈ C is defined as sij = c0i + c0j − cij . A decreasing ordered list of these savings is built by storing at each row the information of the nodes i and j and the saving sij . We use a parallel version of the Clarke and Wright algorithm for split delivery demands described in Aleman and Hill (2010), but slightly modified because we fix the number of vehicles to K = Kmin . We also evaluated the performance of RGTS using two different versions of the Clarke and Wright savings algorithm: a sequential version where the routes are built one by one, and a modified parallel version where splits at interior points are allowed (contrary to the parallel algorithm of Aleman and Hill (2010), which only allows splits at the extremes of the routes). We observed that the parallel version presented by Aleman and Hill (2010) gave better initial solutions for almost all the tested instances. Besides, in those cases where either the sequential algorithm or the modified parallel version improved the solution provided by Aleman and Hill (2010), the final solution of RGTS was not too much better. Therefore, for the sake of simplicity, we only use the parallel Clarke and Wright to obtain the initial solution. 3.1.2 Phase II: Granular tabu search

In this phase, seven different local search operators are considered and a tabu list is built. Each operator searches for an improving move in a granular neighborhood of the current solution and the best move (maximum saving) from among some promising moves randomly selected is used to obtain a new current solution. The move introduced at the current solution is stored in the tabu list and any other move that involves at least one attribute of the move in the list is tabu for the next t iterations. Thus, t is the maximum size of the tabu list. For the local search procedures we use a modified version of the granular threshold g that considers particular information regarding the remaining capacities of the vehicles. This granular threshold gk is defined for each vehicle k as: z C C gk = β 1 + k = gk0 + gk0 k , n+K +s C C

A Randomized Granular Tabu Search Heuristic for the SDVRP

7

where z is the value of the current solution, s is the number of splits of the current solution (which must be added to n + K to obtain the total number of arcs used in this solution), Ck is the remaining capacity of vehicle k, C represents the average of the remaining capacities of the routes that still have some unused capacity and gk0 = β z/(n + K + s). This new term in the definition of gk allows to create routes in a neighborhood of solutions that are not in the normal granular space but that are attractive because they have a relevant remaining capacity. Note that the granular threshold is defined for each route because it depends on the remaining capacity of each route k and it is calculated at each iteration. In the explanations of the local moves we use the following notation. Denoting with rk the route performed by vehicle k, we say that a node i ∈ rk stays in the 0 granular neighborhood of j ∈ rv if cij ≤ gv and we denote it as i ∈ gv (j ). Finally, yik 1 and yik are the demands of customer i served by vehicle k before introducing and after having introduced a move, respectively. The saving cost for each move is weighted by a factor that represents the remaining capacity of the route(s) where the node(s) will be inserted. Let D be the set of arcs to be deleted, F the set of arcs to be added, R the set of routes that receive new demand and H the set of nodes involved in the move. The weighted saving cost for move p is: ! X X X X ca − ca yj , sp = Ck − a∈D

a∈F

k∈R

j∈H

where Ck represents the remaining capacity of route k and yj represents the demand of node j to be relocated. As in the granular definition, the saving costs also favor those moves implying routes with a larger remaining capacity. The insertion cost of adding some nodes to a route is defined as the difference between the added arcs and the erased arcs. Similarly, the erasing cost of deleting some nodes from a route is the difference between the erased arcs and the added arcs. The operators used in our tabu search are the following: – Relocate Node (Relocate) : for customer nodes i ∈ rk and j ∈ rv , if either i ∈ gv (j ) or i ∈ gv (j + 1), and if the remaining capacity of rv allows the move (i.e., after the move, rv would deliver at most the allowed capacity limit), then move node i from rk to rv and locate it after node j on rv . Demands are set as 1 1 0 following: yik = 0 and yiv = yik (see Figure 1). – Exchange Node (Exchange) : for customer nodes i, t ∈ rk and j, p ∈ rv , if either i ∈ gv (p) or i ∈ gv (p +1), and either j ∈ gk (t) or j ∈ gk (t +1), and the remaining capacities of rk and rv allow the move, then move node i from rk to rv , locating it at after node p, and move j from rv to rk , locating it after node t. Demands 1 1 0 1 1 0 are set as following: yik = 0, yiv = yik and yjv = 0, yjk = yjv (see Figure 2). – 2-Opt (2Opt) : for customer nodes i ∈ rk and j ∈ rv , if i ∈ gv (j + 1) and j ∈ gk (i+1), switch the initial part of rk (rv ) until node i (node j ) and reconnect it with the part that starts at node j + 1 on rv (i + 1 on rk ). Demands are set 1 0 1 0 as following: yzv = yzk (yzk = yzv ) for every node z of rk (rv ) up to node i (j ). This move is done if the remaining capacities of vehicles k and v allow it (see

Figure 3). – Exchange Split (ExchSplit) : for a customer node i visited by routes rk and rv (i.e., i ∈ rk and i ∈ rv ) and customer nodes z, z + 1 ∈ rv and j ∈ rk , if either

8

Berbotto, Garc´ıa, Nogales 0 0 j ∈ gv (z ) or j ∈ gv (z + 1), and if yjk − yiv > 0, then erase i from rv and generate a new split at node j by inserting it after node z + 1 on rv . Demands are set as 1 1 0 0 1 0 0 1 0 following: yiv = 0, yik = yik + yiv , yjk = yjk − yiv and yjv = yiv (see Figure 4). – Delete Split and New Split (DelSplitNew) : for a customer node i visited by routes rk and rv and customer nodes j ∈ rk and z, z + 1 ∈ rt , if either j ∈ gt (z ) 0 0 or j ∈ gt (z + 1), yjk − yiv > 0 and the remaining capacity of rt allows the move, then erase i from rv and generate a new split on node j between route k and route t, locating it after node z on rt . Demands are set as following: 1 1 0 0 1 0 0 1 0 yiv = 0, yik = yik + yiv , yjk = yjk − yiv and yjt = yiv (see Figure 5). – Delete Split (DelSplit) : given a customer node i that is visited by more than one vehicle, let B be the set of routes that contain i in their paths and let B 0 be the set of routes without i in their paths. Let r∗ be the route in B 0 with the minimum insertion cost of node i. If for some node p ∈ r∗ , either i ∈ gr∗ (p) or i ∈ gr∗ (p + 1), and the remaining vehicle capacity of r∗ allows the move, then delete i from all routes in B and introduce the node in r∗ after node p. 1 1 Demands are set as following: yik = 0 ∀k ∈ B and yir ∗ = di (see Figure 6). – Delete Worst Split (DelWorstSplit) : let B be the set of routes that contain node i in their paths and let rk ∈ B be the route with the greatest saving cost of erasing i. Sort B by decreasing saving cost, delete i from the first route in B , 0 that is rk , and split yik among the other routes in B starting by the last one of

the list (the one with the lowest saving cost). When the vehicle is full, follow 0 with the next one and so on until having assigned the full amount yik (see Figure 7).

(a) Routes before Relocate.

(b) Routes after Relocate.

Fig. 1: Relocate operator.

The operators Relocate, Exchange and 2Opt are the classical ones developed for the traveling salesman problem and, jointly with with ExchSplit, are applied in Ho and Haugland (2004). DelWorstSplit is presented in Archetti et al. (2006b), but they use the insertion criterion in the inverse sense: starting by the route with the second greater saving cost. We use a different criterion because we consider that a route with a lower saving cost from erasing i has a smaller chance to delete this node in next moves and, thus, its demand must not be relocated. For example, think of a route k that visits the node i located on the segment that connect nodes j and z (that means, cji + ciz = cjz ) and assume that vehicle k serves a large demand yik . In this case, the saving cost from erasing i is equal to zero,

A Randomized Granular Tabu Search Heuristic for the SDVRP

(a) Routes before Exchange.

9

(b) Routes after Exchange.

Fig. 2: Exchange operator.

(a) Routes before 2Opt.

(b) Routes after 2Opt.

Fig. 3: 2Opt operator.

(a) Routes before ExchSplit.

(b) Routes after ExchSplit.

Fig. 4: ExchSplit operator.

the lowest possible because there is no extra cost for visiting i when the vehicle goes from j to z . This means that i is served by vehicle k with a cost of zero. On the other hand, if i presents a large saving cost from being erased from route k, it means that it is a promising node to be moved to other route in some future iteration and relocating yik might be complicated. The operators DelNewSplit and DelSplit are new and they focus on the relocation and deletion of splits. We introduce these operators because in a SDVRP with limited vehicle fleet size it is more important to optimize the split of the demands than when the fleet is unlimited. At every iteration, each operator searches for a move to improve the current solution based on the largest saving cost. For each operator, if there is more than one move with positive saving cost, they are sorted in a list in descendent saving

10

Berbotto, Garc´ıa, Nogales

(a) Routes before DelSplitNew.

(b) Routes after DelSplitNew.

Fig. 5: DelSplitNew operator.

(a) Routes before DelSplit.

(b) Routes after DelSplit.

Fig. 6: DelSplit operator, B={rk , rv }, B’={rt }.

(a) Routes before DelWorstSplit.

(b) Routes after DelWorstSplit.

Fig. 7: DelWorstSplit operator.

cost. Next, the algorithm uses randomness in the selection of the candidate moves for each operator. The idea is that, in order to obtain a better final solution, it is not necessary to pass through the best current solutions at each iteration. Thus, for each operator the algorithm selects randomly one move among the m best candidate moves by assigning higher probabilities to the moves with larger saving costs. The random selection of the candidate moves at each operator is defined as follows: 1. Fix the maximum candidate list size to m. 2. Obtain q ≤ m candidates.

A Randomized Granular Tabu Search Heuristic for the SDVRP

11

3. Assign the following probabilities to each candidate: – p1 = q/qˆ for the candidate with the largest saving cost, – p2 = (q − 1)/qˆ for the candidate with the second largest saving cost, – ... – pq = 1/qˆ for the candidate with the q -th largest saving cost, P where qˆ = qi=1 i. For example, if m = q = 5 for all the operators, the moves for each operator have the following probabilities: p1 = 5/15, p2 = 4/15, p3 = 3/15, p4 = 2/15, and p5 = 1/15. With these probabilities, one move is randomly selected for each operator. Once a candidate move has been selected for each operator, the one with the largest saving cost is implemented in the current solution. This move is stored in the tabu list and the reverse move is forbidden for t iterations. Table 1 shows an example of a possible iteration. Each operator presents a maximum of m = q = 5 candidate moves with positive saving costs and, with the probabilities explained above, the candidates in bold letter are selected. The best move among these selected candidates, ExchSplit, is implemented. Table 1: Tabu Phase with m=5 Moves Relocate Exchange 2opt ExchSplit DelSplitNew DelSplit DelWorstSplit ∗ By

saving cost;

10 7∗∗ 9 10 6 9 8 ∗∗ by

Candidates∗ 8 6∗∗ 4 6 3 2 8 7∗∗ 6 9∗∗ 3 2 5 2∗∗ 1 8∗∗ 7 4 7 3 2∗∗

I 2 1 3 1 0 2 0

⇒

Selected moves∗ 6 7 7 9 2 8 2

II

We choose

⇒

ExchSplit

random selection.

It must be noted that an alternative strategy consisting in selecting the best move for each process and then implementing randomly one of them was also tested but it was discarded because it performed worse. The intensification and diversification criterion are introduced as follows: if an improving move is found at an iteration, then the value of β in gk is reduced by a certain value φ% to intensify the search in this neighborhood. Otherwise, if there is no candidate move with positive saving cost at an iteration, the algorithm allows to introduce a non-positive move and a diversification strategy is applied: the value of β is increased by φ% for the next iteration to explore a new neighborhood. We do not fix a maximum value for β . Thus, at each iteration, either a diversification or an intensification strategy is applied. Finally, if a move is tabu but it leads to an improvement in the current best solution (incumbent solution), it is allowed to be a candidate move. Another important idea is that we allow the algorithm to explore neighborhoods of infeasible solutions in terms of vehicle capacity. The infeasibility in capacity is controlled by a parameter π such that the total demand served by a route cannot exceed the threshold (1 + π )Q, π ≥ 0. Thus, the tabu phase relaxes constraints (6) of the vehicle flow formulation to X yik ≤ (1 + π )Q, k = 1, . . . , K. i∈C

12

Berbotto, Garc´ıa, Nogales

The granular tabu phase stops when there are T consecutive iterations with no positive moves. Besides, this current solution may be infeasible in terms of the capacities of the vehicles, in which case this infeasibility is fixed in the next phase of the algorithm. 3.1.3 Phase III: Capacity correction

This phase is carried out only if there is some vehicle whose total served demand exceeds the capacity limit Q. It presents a set of operators combining up to three routes and up to three nodes that eliminate the excess of delivered product by combining different simple ways of making new splits or/and relocations. Given a route rk with excess of capacity and routes t and v with enough remaining capacity, the following operators are used to move product from the route with excess: 1. Two routes and one node: (a) For node i ∈ rk relocate i on rt . (b) For node i ∈ rk , rt relocate some demand of i from rk to rt . (c) For node i ∈ rk , introduce i in rt and split its demand. 2. Three routes and two nodes: (a) For nodes i, j ∈ rk , relocate them on rt and rv , respectively. (b) For nodes i, j ∈ rk , relocate i on rt , introduce j in rv and split yjk between rk and rv . If there is no operator among the previous ones that eliminates the whole excess of demand of the route, then we consider more complex operators that combine two or three of the ones in the previous list (considering a maximum of up to three nodes). If this new operator allows to remove the whole excess of demand, then it is implemented. The first of these operators (simple or complex) that eliminates totally the excess of delivered product is applied. If the excess of capacity of a route cannot be totally eliminated either with a simplex or a complex operator because there is no route with enough remaining capacity, then partial eliminations are allowed: the excess is reduced a ρ% (20% in our algorithm) and the procedure is restarted for the remaining excess. This guarantees that we will end with a solution where original capacity limits are not violated. 3.1.4 Phase IV: Solution improvement

After the capacity correction, the algorithm improves the obtained solution with three operators: 2-Split Cycles, Node-Position and Subtours-Elimination. The 2-Splits Cycle operator eliminates solutions where two routes share two split nodes. The motivation for this move is that it is shown in Dror and Trudeau (1989) that, if the costs satisfy the triangle inequality, then there is an optimal solution where any two routes share at most one customer. Subtour-Elimination is a new operator that deletes subtours in a given solution (based on constraints (4)) by erasing multiple visits of a vehicle to the same customer. In our method it is possible to obtain a solution with subtours in the previous phase because the local search is made in a bounded neighborhood. Thus, a relocation of node i in a route that already visits this node is something that may happen. Node-Position is another new operator where a node is moved to a different position if this leads to a reduction in the

A Randomized Granular Tabu Search Heuristic for the SDVRP

13

route cost. It tries to move a node only up to three positions backward or forward in the route. Next, the algorithm applies a granular local search with the same procedures of the tabu phase but setting π to 1 and selecting the best move at each step (i.e., the candidate list size is m = 1). Finally, a TSP is solved for every route looking for possible savings by reordering the nodes in the route. We use the standard TSP formulation with the MTZ subtour elimination constraints (Miller et al. (1960)). Only routes with at least four nodes are considered because for the routes with three nodes this is equivalent to the application of the procedure Node-position. We impose a time limit of 20 seconds per node for solving these problems. Because of the presence of randomness in the selection of the move at each local operator, the granular tabu search and capacity correction phases are executed M times. The tabu list is rebooted at each of these executions. To sum up, the structure of the algorithm is the following: 1. Phase I: Obtain an initial solution x0 with the Clarke and Wright algorithm adapted from Aleman and Hill (2010). Define BestSolution=IncumbentSolution:=x0 . 2. Repeat M times: (a) Define InitialSolution:=BestSolution. (b) Run Phase II with InitialSolution as starting solution. (c) Run Phase III and let BestSolution the final solution of this phase. (d) If value(BestSolution)