Multiple neighborhood search, tabu search and ejection chains for the multi-depot open vehicle routing problem

Accepted Manuscript Multiple neighborhood search, tabu search and ejection chains for the multidepot open vehicle routing problem María Soto, Marc Sevaux, Andreas Reinholz, André Rossi PII: DOI: Reference:

S0360-8352(17)30118-3 http://dx.doi.org/10.1016/j.cie.2017.03.022 CAIE 4674

To appear in:

Computers & Industrial Engineering

Received Date: Accepted Date:

26 November 2015 17 March 2017

Please cite this article as: Soto, M., Sevaux, M., Reinholz, A., Rossi, A., Multiple neighborhood search, tabu search and ejection chains for the multi-depot open vehicle routing problem, Computers & Industrial Engineering (2017), doi: http://dx.doi.org/10.1016/j.cie.2017.03.022

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Multiple neighborhood search, tabu search and ejection chains for the multi-depot open vehicle routing problem

Abstract In this paper, we address the Multi-Depot Open Vehicle Routing Problem (MDOVRP), which is a generalization of the Capacitated Vehicle Routing Problem (CVRP) where vehicles start from different depots, visit customers, deliver goods and are not required to return to the depot at the end of their routes. The goal of this paper is two-fold. First, we have developed a general Multiple Neighborhood Search hybridized with a Tabu Search (MNSTS) strategy which is proved to be efficient and second, we have settled an unified view of ejection chains to be able to handle several neighborhoods in a simple manner. The neighborhoods in the proposed MNS-TS are generated from path moves and ejection chains. The numerical and statistical tests carried out over OVRP and MDOVRP problem instances from the literature show that MNS-TS outperforms the state-of-the-art methods. Keywords: CVRP, OVRP, MDOVRP, MNS, TS, ejection chains, path move 1. Introduction The Open Vehicle Routing Problem (OVRP) is a variant of the well known Capacitated Vehicle Routing Problem (CVRP) (Dantzig and Ramser, 1959), in which vehicles do not return to the depot. Each route ends at the last served customer. This problem was first introduce by Schrage (1981), then Sariklis and Powell (2000) introduced this problem as OVRP for the first time. One of the first applications of this problem is a special air cargo routing problem for FedEx introduced by Bodin et al. (1983). OVRP model is appropriate when the vehicles fleet is hired (Sariklis and Powell, 2000). For instance, in the transportation of handicapped persons (Soto et al., 2014), care centers for disabled people must provide daily transport to patients from their home to the care center as well as the return after treatment. These Preprint submitted to Computers & Industrial Engineering

March 18, 2017

care centers work with specialized services companies in the transportation of handicapped persons, which have adapted vehicles for the transportation of disabled people. For other real-life applications of this kind of problem, readers are referred to Li et al. (2007); Repoussis et al. (2010, 2007); Tarantilis et al. (2004, 2005). Furthermore, an excellent literature review about OVRP is provided in (Zachariadis and Kiranoudis, 2010) or in (Li et al., 2007). Some exact methods have been proposed to tackle OVRP. For instance, Letchford et al. (2007) propose a Branch-and-Cut algorithm and Pessoa et al. (2008) propose a more complex Branch-and-Cut-and-Price strategy to solve the problem. Exact methods reach optimal solutions at the cost of prohibitively large computational times for real-life problem instances. Therefore, metaheuristic approaches seem a good alternative to produce high quality solutions in a reasonable computational time. Readers interested to know more about the metaheuristic dedicated to OVRP are referred to (Reinholz and Schneider, 2013; Sariklis and Powell, 2000; MirHassani and Abolghasemi, 2011; Li et al., 2007; Repoussis et al., 2010; Fleszar et al., 2009; Zachariadis and Kiranoudis, 2010; Li et al., 2012; Yu et al., 2011; Derigs and Reuter, 2009). There are many real-life transportation problems, where the vehicles can depart from several depots. In the logistic field, this problem is well-known as the Multi-Depot Vehicle Routing Problem MDVRP (Wren and Holliday, 1972; Chao et al., 1993; Renaud et al., 1996; Cordeau et al., 1997a; Lim and Wang, 2005; Lau et al., 2010). Thus, the Multi-Depot Open Vehicle Routing Problem MDOVRP is a variant of the MDVRP, which is more complex than the single-depot version (OVRP). The MDOVRP has been introduced by Tarantilis and Kiranoudis (2002) to tackle a distribution problem of fresh meat. They proposed a list-based threshold-accepting algorithm. Also, Liu et al. (2014) present a mixed integer programming mathematical formulation and a hybrid genetic algorithm. Recently, Lalla-Ruiz et al. (2015) proposed a new mixed integer programming formulation for the MDOVRP, which improves the subtour elimination constraints proposed by Liu et al. (2014) and provides new sets of constraints. In this paper, we propose a Multiple Variable Neighborhood Search hybridized with a Tabu Search called MNS-TS to tackle the MDOVRP. The structure of this paper is organized as follows. The next section presents the MDOVRP. Then, the following sections describe the ingredients of the MNSTS. Section 3 introduces the Path Move, Section 4 explains the generation of our ejection chains and Section 5 presents the different neighborhoods we are 2

using. Section 6 presents the Tabu Search dedicated to the MDOVRP. Our greedy algorithm to generate initial solutions is presented in Section 7. Section 8 summarizes our MNS-TS approach. Extensive computational results as well as the statistical tests are given in Section 9. Section 10 concludes the paper. 2. Multi-depot open VRP This problem is defined on a complete directed graph G = (V, A), where M = {1, . . . , m} is the set of depots and N = {1, . . . , n} is the set of customers. Thus, V = M ∪ N is the set of all vertices in the graph G, with M ∩ N = ∅. The set of arcs is denoted by A. Each arc (i, j) in A represents connections between two vertices. ∀i ∈ N, qi is the demand of customer i. Depots have a zero demand, hence qi = 0, ∀i ∈ M. Depots have enough goods to deliver to all the customers, and an unlimited fleet of homogeneous vehicles with a capacity Q. Each arc (i, j) is associated with the distance di,j between the two vertices i and j. The traveling length of any vehicle route must be less than a given threshold D. The goal is to determine the vehicle routes which minimize the total traveling length satisfying the following constraints: • each route starts at any depot and finishes at the last visited customer, • each customer must be visited by exactly one vehicle, • the total demand of the customers on the route of any vehicle cannot exceed its capacity Q, • the total length of each vehicle route must not exceed D, the length limit. A solution S of the MDOVRP is represented as a quintuplet S = (X, C, E, R, O). Where matrix X represents the vehicle routes, vector Xk is the vehicle route k and xh,k = i means that vertex i is the h-th vertex visited by vehicle route k. Vector C is the used capacity, ck ∈ R∗+ is the capacity used by the vehicle route k. Vector E is the length of the vehicle routes, ek ∈ R∗+ is the traveled distance performed by the vehicle route k. Vector R is the route assigned to customers, ri = k indicates that customer i is served by the vehicle k. Finally, vector O is the position of the customers in the routes, oi = h means that customer i is the h-th customer in the vehicle route ri . Note that matrix 3

X, vectors R and O carry the same pieces of information, i.e. xh,k = i if and only if ri = k and oi = h. The objective function to minimize is the total traveling length and can be expressed as follows: f (S) =

X k

|Xk |−1

X

d(xh,k , xh+1,k )

(1)

h=1

3. Path and path moves This section introduces the path move, which is at the core of the ejection chains and the neighborhoods of MNS-TS. First, we define a path and a path move. Then, we presents the contribution length function, which computes the contribution length by moving a path. Finally, we define a function which allows us to determine when a path move is feasible or not. 3.1. Path A path is a sequence of customers composed by one or more consecutive customers. A path is denoted by Piα , ∀i ∈ N and α ∈ N∗ . It means that this path starts at customer i in the route ri and visits α customers. A path only composed by one customer i is denoted by Pi . For instance, in the following route: 5, 6, |{z} 7 , 1, 11, 3, 2, 10, 15, 8, 9 | {z } P7

P34

we have path P7 that is only composed by customer 7, and path P34 that starts at customers 3 and is composed by four consecutive customers: 3, 2, 10 and 15. We denote by `(Piα ) the length of path Piα : `(Piα )

=

oiX +α−2

d(xh,ri , xh+1,ri )

(2)

h=oi

Furthermore, we denote by `I (Piα) the length of the reverse path, expressed

4

by: `I (Piα )

=

oX i +1

d(xh,ri , xh−1,ri )

(3)

h=oi +α−1

Finally, Equation (4) explicits δ(Piα ) which is the total demand of the customers that are part of Piα : δ(Piα )

=

oiX +α−1

qxh,ri

(4)

h=oi

3.2. Path move A move consists in removing a path Piα from its route to reinsert it after some customers j in the same route or not. There are two different ways of reinserting a path, because the vertex sequence of a path may be reversed. Reversion is represented by ω in Ω = {1, 2}. Thus, ω is set to 2 if path Piα is reversed and inserted after customer j, otherwise ω is set to 1. A path move is denoted by a triplet (Piα , j, ω), which means that path Piα is inserted in the route of customer j right after it, with the reversion or not indicated by ω. With the goal of reducing the size of the neighborhoods presented in Section 5, for each customer i, we determine a neighborhood N (i) composed by the closest customers of i in G. For instance, the neighborhood of a vertex i can be composed by the vertices j such that the distance from vertex i to vertex j is less than a radius π, i.e., N (i) = {j ∈ V, d(i, j) ≤ π} Furthermore, we define by P(S), the set of all possible path moves. We also define Pα (S, L) as the set of path moves of size α that can be generated from a solution S and a list of customers L. Specifically, it is stated that Pα (S, L) = {(Piα, j, ω) ∈ P(S), i ∈ L, j ∈ N (i), ω ∈ Ω}. In the elementary situation defined by α = 1 and ω = 1, the set of possible path moves of size one is denoted by P(S, L). From the path moves, we can define intra and inter moves, because a path can involve either one or two vehicle routes. For instance, we can define the following well-known moves: • Relocate: It is composed by the path moves of the type (Pi , j, ω = 1), i.e. moving an unitary path either to the same or to a different vehicle route.

5

• 2-Opt: Path moves of type (Piα>1 , xoi −1,ri , ω = 2) constitute the 2-Opt move. It consists in moving a path of size α (greater than one) in the second way to the same position. • Route Splitting: It is composed by the path moves of the type (Piα , j, ω), where j is in M, i.e. moving a path to a new route which starts at the depot j. • Route Merging: It consists in moving a path composed by all customers of a route k next to customer j in the route rj , i.e. it is formed by the |X −1| path moves of type (Px0 ,kk , j, ω), where rj is different from the route k. 3.3. Contribution length function Considering a complete feasible solution S, we evaluate the length variation induced by performing a path move (Piα , j, ω). The function that computes this contribution length is denoted by g and is the sum of functions gR and gI . This function is expressed as follows: g(Piα, j, ω) = gR (Piα , j, ω) + gI (Piα , j, ω)

(5)

Function gR expressed by Equation (6) is the length involved by removing a path Piα from its route ri . If the route exists, i.e. ri 6= 0 and j 6= xoi −1,ri , then gR is equal to the length of the arc that connects the previous and next customers to the path minus the length of arcs that link the path with the route, minus `(Piα ). If the customer j is the previous customer of the path, xoi −1,ri , and ω = 2, then the length of the arc that connects the previous and next customers to the path is not added to function gR . Otherwise the function is equal to zero.  d(xoi −1,ri , xoi +α,ri ) − d(xoi −1,ri , i) − d(xoi +α−1,ri , xoi +α,ri ) − `(Piα), if    r 6= 0, j ∈ N \ {x i oi −1,ri } α gR (Pi , j, ω) =  −d(j, i) − d(xoi +α−1,ri , xoi +α,ri ) − `(Piα), if ri 6= 0, j = xoi −1,ri , ω = 2    0, otherwise. (6) Function gI stated by Equation (7) is the contribution length produced by inserting a path Piα in the same route as the one of customer j, right 6

after it. This function depends on ω to insert a path and the nature of the vertex j. When vertex j is a depot, it corresponds to splitting or creating a route. If ω = 1 the inserted length is equal to the length of the arc that connects the depot j and vertex i plus the length of the route. If ω = 2, then gI is equal to the length of the arc that links the depot and the last customer of the path plus `I (Piα ). When the customer j is the previous customer of the path, xoi −1,ri , and the way of inserting path is ω = 2, this is a 2-Opt move. Then, the length is equal to the distance between the customer j and the ending customer of the path plus the distance between the customer i and the subsequent customer to the path plus `I (Piα ). When the vertex j is a customer different from xoi −1,ri , it corresponds to reallocating a path or merging routes. Thus, the inserting length is equal to the length of the arc that connected the route with the path plus the length of the path either reversed or not, depending on ω, minus the length of the arc that link customer j with the following customer in the route rj .   d(j, i) + `(Piα ), if j ∈ M, ω = 1    α  d(j, xoi +α−1,ri ) + `I (Pi ), if j ∈ M, ω = 2      d(j, i) + d(xoi +α−1,ri , xoj +1,rj ) − d(j, xoj +1,rj ) + ` (Piα), if ω = 1,     j ∈ N \ {xoi −1,ri } α gI (Pi , j, ω) =  d(j, xoi +α−1,ri ) + d(i, xoj +1,rj ) − d(j, xoj +1,rj ) + `I (Piα ), if ω = 2,      j ∈ N \ {xoi −1,ri }    α   d(j, xoi +α−1,ri ) + d(i, xoi +α,ri ) + `I (Pi ), if j = xoi −1,ri , ω = 2    0, otherwise. (7) 3.4. Feasibility of path moves A path move is feasible if the capacity of the involved vehicle routes is respected as well as the length of the concerned vehicle routes. We defined the function h to determine if a path move (Piα , j, ω) can be performed in the feasible solution S = (X, C, E, R, O). This function is expressed by Equation (8). Two cases must be considered. The first one happens when the move involves only one route. In that case, checking the capacity limit of the route is not necessary, because no customer is added to the route. Only the length of the route has been modified. Thus, we need to check if the current length 7

of the route, erj , plus the contribution length of the path, g(Piα, j, ω), is less than the route length limit D. The second case is when the move involves two routes. The route where the path comes from, ri , is reduced in capacity and length, thus this route satisfies the feasibility condition, because the solution S is already feasible. We verify that the capacity of the route where the path is reinserted, crj , plus the demand of the path, δ(Piα), satisfies the vehicle capacity constraint. We also check that the current length of the route, erj , plus the length of the inserted path, gI (Piα , j, ω), is less than D. If the move is feasible, then the function h is equal to the contribution length of the path, otherwise h is set to infinity to indicate that the move is not feasible.  α α  g(Pi , j, ω), if rj = ri , erj + g(Pi , j, ω) ≤ D h(Piα , j, ω) = g(Piα, j, ω), if rj 6= ri , crj + δ(Piα ) ≤ Q, erj + gI (Piα , j, ω) ≤ D   +∞, otherwise. (8) 4. Ejection chains The neighborhoods used for the MDOVRP are based on the ejection chains introduced by Glover and Laguna (1997). Lin and Kernighan (1973) have been the pioneers using ejection chains by proposing a heuristic algorithm for the traveling salesman problem based of this new type of moves. Furthermore, ejection chains have been successfully applied to tackle the N P-hard problems presented in the electronics field (Sevaux et al., 2014; Soto et al., 2013, 2015). Ejection chains belong to variable depth procedures (Glover and Laguna, 1997). They are complex and powerful moves composed by a series of simple moves. The ejection chains used for the MDOVRP are the compound moves constituted by a series of paths moves that involves two different routes. In this paper, the ejection chain generated from the path moves is denoted by C. An ejection chain is a set of path moves C = {(Piα, j, ω), . . .}, where the first element is the path move which triggers the ejection chain. In fact, moving a path or “ejected from” its current state leads to trigger a cascading sequence of path moves, like a domino effect. Note that a path move can be seen as an unitary ejection chain. For instance, the ejection chain C = {(Piα , j, ω)} is only composed by the path move (Piα , j, ω). 8

Each ejection chain is associated to a contribution length that is the sum of the contribution length of each path move involved in the ejection chains. The contribution length function of an ejection chains is denoted by Γ(C) and is expressed as follows: X Γ(C) = g(Piα0 , j 0 , ω) (9) (Piα0 ,j 0 ,ω) ∈ C

We have studied two ways to generate ejection chains. The first one consists in producing ejection chains from an infeasible path move. The other one consists in generating ejection chains only from the feasible moves. In the following two paragraphs, we explain how to generate ejection chains from infeasible and feasible path moves. 4.1. Ejection chains from infeasible path moves The operator which produces ejection chains from an infeasible path move, (Piα , j, ω) is named ECβ (Piα , j, ω) and Figure 1 presents its pseudocode. This operator is only applied to path moves that do not respect the constraints of capacity or length i.e., when h(Piα , j, ω) = +∞. The operator ECβ () receives an infeasible path move that triggers other path moves of size β that make the ejection chain returning a feasible solution. Thus, this operator finishes when all the routes satisfy the capacity and length constraints. Initially, the infeasible route k is the route where the path Piα is moved, rj . The set L of customers to scan is equal to the customers visited by route k, and the set of forbidden routes I is empty. The set I avoids the infinite cycles between routes when building an ejection chain. In fact, infeasible path moves can be accepted in the ejection chain generating cycles between routes. For instance, path Paα from the infeasible route ri is moved to route rj , thus ri becomes feasible but not rj . Then, infeasible path Pbα from route rj is moved to route ri , causing ri infeasible again. And then, path Pcα from ri is moved to rj , letting rj infeasible in its turn. This can continue infinitely. In order to avoid cycles, problematic routes that have become feasible are stored in the set I. Also, we define the auxiliary set A which is constituted by the path moves (Paβ , b, ω) in Pβ (S, L), where a path move involving a forbidden route is accepted only if it is feasible (h(Paβ , b, ω) 6= +∞). This forces the operator to search for other path moves that either do not involve 9

forbidden routes or involve forbidden routes that will be combined to yield feasible solutions. At each iteration the operator searches for a path move (Paβ , b, ω) in A which produces the minimum contribution length g. Then, the path move is added in the ejection chain C. If the current problematic route k is feasible, then we search for another problematic route k 0 and update the set of forbidden routes I and the set of customers L. By construction, ECβ guarantees the generation of finite ejection chains that produce feasible solutions. There are two aspects that ensure the finite feasible ejection chains. The first one is the use of sets I and A, which breaks the cycles. The other aspect is the use of the path moves which split the routes. This kind of path moves always allow an infeasible route to become feasible by splitting it in several feasible routes.

1

2 3 4 5 6 7 8

9 10 11 12 13 14 15

Input: Solution S, size of path moves, β, and trigger infeasible path move (Piα , j, ω), such that h(Piα , j, ω) = +∞ Output: Ejection chain C C ← {(Piα , j, ω)} Infeasible route k ← rj L = {i ∈ N, ri = k} Forbidden set of routes I = ∅ while k 6= 0 do (Piβ0 , j 0 , ω 0 ) ← arg min g(Paβ , b, ω), where

∀(Paβ ,b,ω)∈A β A = {(Pa , b, ω) ∈ Pβ (S, L), h(Paβ , b, ω) 6= +∞ ∀rb ∈ C ← C ∪ {(Piβ0 , j 0 , ω 0 )} if ck + δ(Piβ0 ) ≤ Q and erj + gI (Piβ0 , j 0 , ω) < D then

I}

I ← I ∪ {k} Search for a new infeasible route k 0 k ← k 0 L = {i ∈ N, ri = k} end end

Figure 1: Operator ECβ (Piα , j, ω) to generate ejection chains from infeasible path moves

10

4.2. Ejection chains from feasible path moves The operator which produces ejection chains from feasible path moves is named ECFβ and Figure 2 presents its pseudo-code. This operator is only applied to path moves that respect the constraints of capacity and length. It receives a path move (Piα , j, ω) and searches for the feasible path move in A with the minimum contribution. The set A contains the path moves of size β, (Paβ , b, ω), in Pβ (S, L) that only involve the routes concerned in the trigger path move, i.e. ra = rj and rb = ri . Note that ECFβ also generates path exchange like “Path Exchanging Operator” introduced by Reinholz and Schneider (2013). For instance, when β and α are equal to 1, the operator ECFβ generates a swap between the involved routes. When the path moves include the last customer of a route, ECFβ is equivalent to “End Restricted Path Exchange” (Reinholz and Schneider, 2013) which consists in cutting two routes into two paths and exchanging their ends.

1

2 3 4

Input: Solution S, size of path moves, β, and feasible path move (Piα , j, ω), such that h(Piα , j, ω) 6= +∞ Output: Ejection chain C L = {i ∈ N, ri = rj } (Piβ0 , j 0 , ω 0) ← arg min h(Paβ , b, ω), A=

5

C←

∀(Paβ ,b,ω)∈A β {(Pa , b, ω) ∈ Pβ (S, L), ra {(Piα , j, ω), (Piβ0 , j 0 , ω 0)}

= rj , rb = ri }

Figure 2: Operator ECFβ (Piα , j, ω) to generate ejection chains from feasible path moves

5. Neighborhoods From the path moves and ejection chains we are able to define many neighborhoods. In fact, we can generate different neighborhoods by adapting the parameters of the path moves. However, this section defines only the neighborhoods used by our metaheuristic. We determine our neighborhoods based on a given feasible solution S = (X, C, E, R, O). Thus, we propose neighborhoods generated from the path moves which are only constituted 11

by intra route moves. We also propose neighborhoods generated from the ejection chains which are composed by either inter or intra route moves. These neighborhoods are detailed in the following paragraphs. 5.1. Intra route neighborhoods Relocate neighborhood NR (S). is composed by the path moves of size one that are moved to another position in the route. It can be expressed as follows: NR (S) = {(Pi , j, ω) ∈ P(S), ω = 1, i ∈ N, j ∈ N (i), ri = rj } 2-Opt neighborhood N2opt (S). is constituted by the path moves that are reversed and moved to the same position. It is expressed as follows: N2opt (S) = {(Piα , j, ω) ∈ P(S), ∀α, ω = 2, i ∈ N, j ∈ N (i), j = xoi −1,ri } 5.2. Inter route neighborhoods Inter route relocation neighborhood NIR (S). moves a path of size one to another route provided that the path move is feasible: NRI (S) = {(Pi , j, ω) ∈ P(S), ω = 1, i ∈ N, j ∈ N (i), ri 6= rj , h(Pi , j, ω) 6= +∞} Inter route large relocation neighborhood NBR (S). is composed by the feasible path moves that relocate the paths containing the last customer of routes: NBR (S) = {(Piα, j, ω) ∈ P(S), α = |Xri | − ri + 1, i ∈ N, j ∈ N (i), ri 6= rj , h(Pi , j, ω) 6= +∞} 5.3. Ejection chains from infeasible moves We now present the neighborhoods composed by ejections chains generated from infeasible path moves applying the operator ECβ . These neighborhoods are inter route. Neighborhood NαE (S). for α ∈ {1, 2, 3} are constituted by the ejection chains generated from the infeasible path moves of size α whose contribution length is minimum: NαE (S) = {EC(Piα0 , j 0 , ω 0 ), (Piα0 , j 0 , ω 0 ) = arg min g(Piα, j, ω)} (Piα ,j,ω)∈Pα (S,N )

12

Neighborhood NBE (S). is composed by the ejection chain generated from the infeasible path move containing the last customer whose contribution length is minimum: NBE (S) = {EC(Piα0 , j 0 , ω 0 ), (Piα0 , j 0 , ω 0) = arg min g(Piα, j, ω), α = |Xri |− (Piα ,j,ω)∈Pα (S,N )

ri + 1} 5.4. Ejection chains from feasible moves We describe the neighborhoods generated from the operator ECFβ applied to feasible path moves. Swap neighborhood Nswap (S). is a special case of ejection chains. It is composed by the ejections chains that swap two paths of size one. Nswap (S) = {ECFβ (Pi , j, ω), β = 1, (Pi , j, ω) ∈ P(S, N), h(Pi , j, ω) 6= +∞} End restricted path exchange neighborhood Nex (S). is composed by the ejection chains that exchange two paths containing the last customers of the routes. Nex (S) = {ECFβ (Piα , j, ω), (Piα, j, ω) ∈ Pα (S, N), h(Piα , j, ω) 6= +∞, α = |Xri | − ri + 1, β = |Xrj | − rj + 1} 5.5. Other neighborhoods Finally, we propose two composed neighborhoods. The first one called N intra which is composed by the intra route neighborhoods. The other one denoted by N inter composed by almost all inter route neighborhoods except Neighborhoods NαE (S) for α equals to 2 and 3. 6. Tabu search Tabu search (Glover and Laguna, 1997) is a local search procedure which iteratively moves from the current solution to another one. To escape from the local optimum, tabu search employs a tabu list which can be seen as the memory of the metaheuristic. The tabu list keeps the solutions that have been explorer in the recent past. A solution presented in the tabu list is a forbidden one called “tabu solution”. At each iteration, tabu search accepts a solution which is the best non-tabu solution in the neighborhood but might be worse than the current solution. It allows the metaheuristic to visit unexplored regions.

13

In the tabu search developed for the MDOVRP, the tabu list keeps the set of customers visited by a route. Thus, a path move is set tabu if the set of customers visited by at least one of the concerned routes is in the tabu list. An ejection chain is tabu if one of its path moves is tabu. A tabu ejection chain can be performed if it produces a solution which is better than the current best solution. This is known as the aspiration criterion (Glover and Laguna, 1997). Figure 3 is the pseudo-code for the tabu search, T S(N ), developed to explore our neighborhoods. It is applied to a neighborhood N and a feasible initial solution S, and returns the best solution reached. At each iteration, it searches for the non tabu ejection chains associated with the minimum contribution length to perform it and obtain a new solution S 0 . The aspiration criterion is applied too. The solution S 0 is accepted as the current one even if it is not better than the best solution. Then, the best solution between S ∗ and S is kept, i.e the solution whose objective value f is the minimum will be the best one. And the tabu list is updated on the FIFO basis (First In First Out) with the routes concerned by path moves in the performed ejection chain. The search finishes when the stopping criterion is met, i.e. when the best solution is not improved after a fixed number of iterations.

1 2 3 4 5 6 7

Input: Valid solution S and Neighborhood N Output: Best solution found, S ∗ Initially f (S ∗ ) ← +∞ Iterative phase while stopping condition is not reached do C ← arg minΓ(C 0 ), where C 0 is not tabu or f (C 0 ) < f (S ∗ ) C 0 ∈N

8 9 10 11 12

Perform ejection chain C to obtain S 0 Do S ← S 0 Keep best solution between S and S ∗ Update the tabu list and the stopping criterion end Figure 3: Tabu Search, T S(N )

14

7. Initial Solution We propose an insertion greedy algorithm to generate initial solutions for the MDOVRP. Figure 4 presents the pseudo-code for this algorithm. Initially, the solution S is empty, vectors C, E, R, O involved in the solution are set to 0. For each depot, we set rk different to zero to trigger the algorithm. And, the list of customers to scan, L, is composed by all customers. At each iteration, we search for the path move, associated with the minimum contribution length. Then, the path move is performed in the solution and the inserted customers is removed from the list L. The algorithm terminates when all customers have been assigned to a route.

1 2 3 4 5 6 7 8 9 10

Initially S = (X, 0, 0, 0, 0), X = 0, C = 0, E = 0, R = 0 and C = 0 rj 6= 0 for all k in M L←N Insertion phase repeat (Pi0 , j 0 , ω 0) ← arg min h(Pi , Pj , ω) Perform the L ← L \ {i0 } until L = ∅;

(Pi ,j,ω)∈P (S,L),rj 6=0 path move (Pi0 , j 0 , ω 0)

in solution S

Figure 4: Greedy insertion algorithm for the MDOVRP

8. Multiple neighborhood search hybridized with tabu search (MNSTS) Multiple neighborhood search MNS is based on the principle of systematically taking advantage of several optimization approaches. It has been successfully applied for addressing logistic problems such as Vehicle Routing Problems (S¨orensen et al., 2008; Reinholz and Schneider, 2013; Jin et al., 2012). Furthermore, because of its adaptability, MNS is very popular in the commercial VRP packages, where it constitutes the main mechanism of intensification and diversification in the search process (S¨orensen et al., 2008).

15

The MNS-TS developed for the MDOVRP applies a sequential strategy for combining neighborhoods called token-ring search (Gaspero and Schaerf, 2003, 2006). This strategy changes the neighborhoods in a circular way, always starting from the best solution found by the previous one. The token-ring search has shown good results addressing various problems such as timetabling problems (Gaspero and Schaerf, 2006, 2003), heterogeneous vehicle routing problem with time windows and carrier-dependent costs (Ceschia et al., 2011), or balanced graph partitioning problem (Benlic and Hao, 2011) among others. Figure 5 presents the pseudo-code for this hybrid metaheuristic. This framework returns the best solution reached. It starts with the initial solution S0 generated by our greedy insertion algorithm. Then, each route of the solution S0 is improved by the neighborhood N intra to obtain the solution S. Initially, the best solution is equal to the current solution. At each iteration, we explore the two largest neighborhoods N inter and N intra with the tabu search. Then the best solution is kept. To diversify the search, a neighborhood based on ejection chains is explored by the tabu search. These neighborhoods change in token-ring strategy, at each iteration. The hybrid approach terminates when the best solution can not be improved after a fixed number of iterations. 9. Numerical experiments This section presents the computational experiments carried out for assessing the performance of the MNS-TS proposed in this work. The proposal has been implemented in C++ and compiled with gcc 4.11 in Linux OS 12.10 using an Intel Core i7-2600 processor system at 3.40 GHz with 16 GB RAM. The computational experiments have been conducted over the set of well-known instances for the OVRP (single-depot case) and MDOVRP. 9.1. OVRP results To assess the performance of the MNS-TS for the single-depot case, i.e. for the OVRP, we have employed the so-called A, B, E, F, M, and P benchmark CVRP instances1 . Also, we have done our experiments on the benchmark instances from Christofides et al. (1979) and Fisher and Jaikumar (1981). 1

CVRP instances are available at www.branchandcut.org

16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

Output: Best solution found, S ∗ Initially Initial Solution S0 generated by the greedy insertion algorithm S ← T S(N intra(S0 )) S∗ ← S k←1 Iterative phase while stopping condition is not reached do S ← T S(N inter(S)) S ← T S(N intra(S)) Update S ∗ if k = 1 then S ← T S(NBE (S)) if k = 2 then S ← T S(N2E (S)) if k = 3 then S ← T S(N3E (S)) Update S ∗ if k = 3 then k ← k + 1 else k ← 1 Update the stopping criterion end Figure 5: M N S − T S(N )

Our MNS-TS have been tested over 110 mainstream instances for the OVRP. In all of these instances, a customer is represented as a vertex located in the Euclidean plane. The cost of an arc is equal to the Euclidean distance between its end-vertices. The instance name indicates the source of the instance, the number of vertices n, the number of vehicles k which is fixed at the minimum possible, and the length limit, D. For example, E-n51-k5 is taken from the E benchmark (Christofides and Eilon, 1969) has 51 vertices and 4 vehicles and does not have a limit for the length of the routes. C06-n50-d200 is taken from Christofides et al. (1979) has 50 vertices and the length of the route must not exceed 200 units. A Branch-and-Cut (Letchford et al., 2007) and Branch-Cut-and-Price (Pessoa et al., 2008) have been proposed to exactly tackle the OVRP. These two methods will be used for comparison to either the optimal solution when 17

known or to the best lower bound. To our knowledge, the best metaheuristics developed for the OVRP are the following ones: a Hybrid (1+1)-Evolutionary Strategy proposed by Reinholz and Schneider (2013), an ILP improvement procedure proposed Salari et al. (2010), and a Tabu Search Heuristic (ABHC) proposed by Derigs and Reuter (2009). For each instance, we have chosen to compare our approach to the metaheuristic which reaches the best function value for the OVRP. This best value will be named “best metaheuristic” in the following tables. We also compare the performance of the MNS-TS is compared to the best exact method. After some computational tests, we have set the parameters involved in the execution of our approach as follows. The stopping criterion is set to 10 iterations without any improvement in the best solution. The size of the tabu list is set to 20. Thus, Table 1 displays the numerical results from our computational experiments. The first column is the instance name. The second column is the lower bound reported by the exact methods; the values with an asterisk (*) are the known optimal solutions. The upper bound has not been reported here because the metaheuristics from the literature reach better solutions when optimal solution is not known. The following four columns are the information about the best metaheuristics for the OVRP such as the name of the best metaheuristic, the value of the objective function f . We also report here the improvement of the best metaheuristic over the lower bound or the optimal solution provided by the exact methods, which is denoted by ML and expressed in percentage. Since the different methods have been executed on different computers, cpu* times have been scaled according to Linpack benchmarks2 . The last four columns show the information on the MNS-TS such as objective value reached for each instance and the computational time spent. Also, we compute the distance of MNS-TS to the lower bound or to the optimal solution, ML , and the improvement over the best metaheuristic of the literature, MB , both of them are expressed in percentage. The best solution values are displayed in bold. In the last column, the values in bold mean the MNS-TS is better than the best metaheuristic, and values in italics mean the best metaheuristic is better than the MNS-TS. 2

http://www.roylongbottom.org.uk/linpack%20results.htm

18

At the bottom of Table 1, we summarize the main results about computational comparison such as the number of best solutions returned by each method, the average computational time and the average improvement brought by MNS-TS over its best competitors. The best metaheuristics in the literature reach the best solution for 95 instances. MNS-TS returns the optimal solutions for all instances except two when it is known. Furthermore, MNS-TS reaches the best solution for 104 instances out of 110 and improves the best solution of the literature by an average 0.01%. The competitors are 0.11% over the lower bound while MNSTS is 0.10% over the lower bound. It means that our MNS-TS has improved the difference between the lower bound in 11.78% on average. Note that all methods are able to reach a large number of optimal solution, introducing a bias in the average results. Furthermore, on average the computational time spent by the MNS-TS is much less than the one spent by the best metaheuristics from the state-of-the-art. Table 1: Comparison of best literature metaheuristic of the on and MNS-TS Instance Name P-n16-k8 P-n19-k2 P-n20-k2 P-n21-k2 E-n22-k4 P-n22-k2 P-n22-k8 E-n23-k3 P-n23-k8 E-n30-k3 B-n31-k5 A-n32-k5 A-n33-k5 A-n33-k6 E-n33-k4 A-n34-k5 B-n34-k5 B-n35-k5 A-n36-k5 A-n37-k5 A-n37-k6 A-n38-k5 B-n38-k6 A-n39-k5 A-n39-k6 B-n39-k5

Lower bound

Name

235.06* 168.57* 170.28* 163.88* 252.61* 167.19* 345.87* 442.98* 302.87* 393.51* 362.73* 487.31* 424.54* 462.43* 511.26* 508.17* 458.76* 557.33* 519.46* 486.24* 581.07* 498* 445.63* 549.68* 533.07* 322.54*

ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC

Best metaheuristic f cpu* 235.06 168.57 170.28 163.88 252.61 167.19 345.87 442.98 302.87 393.51 362.73 487.31 424.54 462.43 511.26 508.17 458.76 557.33 519.46 486.24 581.07 498.00 445.63 549.68 533.07 322.54

19

0.70 1.10 1.20 1.30 1.20 1.50 1.10 1.70 0.90 2.00 2.50 2.60 2.60 2.40 2.40 2.50 2.70 2.70 2.80 3.10 2.30 2.50 3.70 2.80 3.30 3.80

ML % 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

f 235.06 168.57 170.28 163.88 252.62 167.19 345.87 442.99 302.87 393.51 362.73 487.31 424.54 462.43 511.26 508.17 458.76 557.33 519.46 486.24 581.07 498.00 445.63 549.69 533.07 322.54

MNS-TS cpu ML %

Mf %

0.51 0.00 0.18 0.00 0.20 0.00 0.16 0.00 0.24 0.00 0.23 0.00 1.11 0.00 0.28 0.00 12.92 0.00 0.55 0.00 8.46 0.00 0.33 0.00 0.40 0.00 0.47 0.00 0.54 0.00 0.69 0.00 0.48 0.00 1.32 0.00 0.92 0.00 1.42 0.00 0.92 0.00 0.64 0.00 0.81 0.00 1.24 0.00 0.56 0.00 0.81 0.00 continued on next

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 page

Table 1 – continued from previous page Instance Name P-n40-k5 B-n41-k6 B-n43-k6 A-n44-k6 B-n44-k7 A-n45-k6 A-n45-k7 B-n45-k5 B-n45-k6 F-n45-k4 P-n45-k5 A-n46-k7 A-n48-k7 B-n50-k7 B-n50-k8 C01-n50 C06-n50-d200 P-n50-k10 P-n50-k7 P-n50-k8 B-n51-k7 E-n51-k5 P-n51-k10 B-n52-k7 A-n53-k7 A-n54-k7 A-n55-k9 P-n55-k10 P-n55-k15 P-n55-k7 P-n55-k8 B-n56-k7 B-n57-k7 B-n57-k9 A-n60-k9 P-n60-k10 P-n60-k15 A-n61-k9 A-n62-k8 A-n63-k10 A-n63-k9 B-n63-k10 A-n64-k9 B-n64-k9 A-n65-k9 P-n65-k10 B-n66-k9 B-n67-k10 B-n68-k9 A-n69-k9 P-n70-k10 F11-n72

Lower bound

name

349.55* 483.07* 428.17* 617.385* 501.31* 648.67* 685.16* 488.07* 403.81* 391.81* 583.54* 669.83* 437.15* 720.43 440.44* 397.38* 422.99 625.14* 416.06* 441.19* 655.18* 709.27* 669.06* 444.31* 411.58* 412.55* 420.49* 646.36* 482.09* 783.18* 775.26 837.07* 839.25 520.47* 728.59* 522.5* 755.27* 701.71* 757.76* 548.34 -

ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC Hybrid Hybrid ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC Hybrid

Best metaheuristic f cpu* 349.55 483.07 428.17 617.39 501.31 654.01 685.16 488.07 403.81 463.90 391.81 583.54 669.83 437.15 720.79 412.95 412.95 440.44 397.38 417.75 625.14 416.06 483.43 441.19 655.18 709.27 669.06 444.31 541.64 411.58 412.55 420.48 617.97 869.32 798.01 482.09 569.43 657.81 783.18 778.94 941.53 837.07 848.16 520.47 728.59 522.50 755.27 618.92 701.72 757.76 552.65 177.00

20

3.30 3.10 4.20 3.30 3.70 4.20 3.60 3.50 3.10 4.60 4.00 4.40 4.20 5.40 4.30 3.62 2.42 3.60 4.60 4.50 4.50 4.00 3.40 5.50 4.70 4.40 4.90 5.00 4.60 5.80 5.20 6.70 7.60 5.60 5.90 5.10 5.00 7.00 5.70 6.10 5.10 6.30 5.70 5.60 5.70 6.30 5.80 7.80 7.10 6.90 5.80 261.18

ML % 0.00 0.00 0.00 0.00 0.00 0.82 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.05 0.00 0.00 -1.24 0.00 0.00 2.24 0.00 0.00 0.00 0.00 0.00 -0.50 0.00 0.00 0.00 -4.39 0.15 2.26 0.00 1.62 -1.01 0.00 0.47 0.62 0.00 1.06 0.00 0.00 0.00 0.00 0.46 0.00 0.00 0.79 0.00

f 349.55 483.07 428.17 617.39 501.31 648.67 685.16 488.07 403.81 463.89 391.81 583.54 669.83 437.15 722.04 412.95 412.95 440.44 397.38 417.75 625.14 416.06 482.59 441.19 655.19 709.27 669.06 444.31 541.64 411.58 412.55 420.48 617.97 869.31 798.01 482.09 569.43 657.81 783.18 778.94 941.53 837.07 848.15 520.46 728.59 522.50 755.27 617.43 701.71 757.76 552.65 177.00

MNS-TS cpu ML %

MB %

1.91 0.00 0.00 2.27 0.00 0.00 2.53 0.00 0.00 0.77 0.00 0.00 1.53 0.00 0.00 15.36 0.00 -0.82 7.53 0.00 0.00 1.10 0.00 0.00 25.73 0.00 0.00 0.36 0.00 0.00 0.43 0.00 0.00 0.77 0.00 0.00 1.71 0.00 0.00 1.55 0.00 0.00 3.67 0.22 0.17 0.19 0.00 1.01 0.00 10.80 0.00 0.00 3.15 0.00 0.00 1.85 -1.24 0.00 1.50 0.00 0.00 1.06 0.00 0.00 2.93 2.06 -0.17 1.33 0.00 0.00 2.30 0.00 0.00 2.18 0.00 0.00 1.83 0.00 0.00 8.16 0.00 0.00 6.78 -0.50 0.00 1.77 0.00 0.00 0.70 0.00 0.00 1.93 0.00 0.00 0.79 -4.39 0.00 20.84 0.15 0.00 3.72 2.26 0.00 6.39 0.00 0.00 23.25 1.62 0.00 11.67 -1.01 0.00 8.51 0.00 0.00 2.00 0.47 0.00 22.62 0.62 0.00 25.71 0.00 0.00 14.01 1.06 0.00 1.74 0.00 0.00 2.57 0.00 0.00 1.36 0.00 0.00 245.01 0.00 0.00 2.80 0.22 -0.24 18.14 0.00 0.00 13.69 0.00 0.00 2.90 0.79 0.00 24.62 0.00 0.00 continued on next page

Table 1 – continued from previous page Instance Name

Lower bound

name

C02-n75 C07-n75-d160 E-n76-k10 E-n76-k14 E-n76-k7 E-n76-k8 P-n76-k4 P-n76-k5 B-n78-k10 A-n80-k10 C03-n100 C08-n100-D230 C12-n100 C14-n100-D1040 E-n101-k14 E-n101-k8 P-n101-k4 C11-n120 C13-n120-D720 F12-n135 C04-n150 C09-n150-D200 C05-n200 C10-n200-D200 G5-n201-D1800 G1-n241-D650 G6-n281-D1500 G2-n321-D900 G7-n361-D1300 G3-n401-D1200 G8-n441-D1200 G4-n481-D1600

559.62 621.27 530.02* 537.24* 522.95* 525.64* 1,061.90 534.24 704.69 639.74* 621.75* 762.89 -

Hybrid Hybrid ABHC ABHC ABHC ABHC ABHC ABHC ABHC ABHC Hybrid Hybrid Hybrid Hybrid ABHC ABHC ABHC Hybrid Hybrid Hybrid Hybrid Hybrid Hybrid Hybrid ABHC ILP-Im ABHC ILP-Im ILP-Im ABHC ILP-Im ABHC

# best solutions Average

Best metaheuristic f cpu*

ML %

f

564.06 566.93 567.14 625.36 530.02 537.24 523.21 525.64 725.79 1,071.01 639.25 642.11 534.24 552.64 713.20 640.42 621.75 678.54 836.55 769.55 733.13 741.44 868.44 871.58 6,018.52 4,573.53 7,731.46 7,251.74 9,156.74 9,803.80 10,344.37 12,420.16

11.16 35.35 6.70 6.70 9.80 8.20 9.60 8.20 8.60 8.50 748.73 301.86 3.67 19.97 14.00 15.80 18.60 55.70 1,642.66 1,414.14 401.79 797.79 3,640.47 3,655.04 36.50 31.88 404.70 42.20 50.45 97.70 66.19 883.70

1.34 0.66 0.00 0.00 0.05 0.00 1.45 0.86 1.21 0.11 0.00 0.87 -

564.06 566.93 567.14 624.62 530.02 537.24 523.21 525.63 724.55 1,069.29 639.25 642.11 534.24 552.64 711.80 639.74 625.97 685.08 837.01 769.55 733.13 741.43 868.11 869.94 6,018.50 4,574.12 7,705.52 7,269.38 9,138.82 9,801.86 10,338.42 12,390.15

136.61

0.11

95

MNS-TS cpu ML %

MB %

0.74 0.86 7.75 50.16 2.09 3.23 1.83 2.51 22.88 17.10 7.53 1.47 45.85 9.50 56.61 2.95 1.49 348.48 479.27 47.40 37.36 13.46 974.74 952.18 34.87 583.03 281.96 144.45 350.16 316.36 934.53 834.91

1.34 0.54 0.00 0.00 0.05 0.00 1.28 0.70 1.01 0.00 0.68 0.87 -

0.00 0.00 0.00 -0.12 0.00 0.00 0.00 0.00 -0.17 -0.16 0.00 0.00 0.00 0.00 -0.20 -0.11 0.68 0.96 0.05 0.00 0.00 0.00 -0.04 -0.19 0.00 0.01 -0.34 0.24 -0.20 -0.02 -0.06 -0.24

65.21

0.10

-0.01

104

9.2. MDOVRP results To assess the performance of the MNS-TS for the MDOVRP, we conducted computational experiments over the set of instances introduced by Cordeau et al. (1997b).

21

Table 2: Comparison of MNS-TS, MIP, HGA

Instances Name n m

UB

MIP LB

4 4 5 2 2 3 4 2 3 4 5 2 4 6 4 4 4 4 4 4 6 6 6 6

386.18 375.93 474.57 662.22 608.73 611.99 613.96 3,156.89 2,795.43 2,704.64 2,844.66 953.26 1,883.53 2,818.36 647.03 978.82 1,423.48 1,523.28 1,817.52 2,221.94 821.25 1,254.45 1,659.39 2,295.01

* * * 617.04 599.56 591.73 582.84 2,397.05 2,262.03 2,231.19 2,180.01 * * * * * 1,421.06 1,468.56 1,606.47 1,905.89 * 1,252.40 1,548.17 1,884.18

# best solutions Average

15

50 50 75 100 100 100 100 249 249 249 249 80 160 240 48 96 144 192 240 288 72 144 216 288

22

p01 p02 p03 p04 p05 p06 p07 p08 p09 p10 p11 p12 p15 p18 pr01 pr02 pr03 pr04 pr05 pr06 pr07 pr08 pr09 pr10

cpu*

f

HGA f

∗

44.23 386.18 0.65 375.93 44.22 474.57 4,703.94 667.02 4,703.94 612.30 4,703.94 614.93 4,703.94 615.24 4,703.94 2,913.63 4,703.94 2,738.54 4,703.94 2,577.48 4,703.94 2,561.85 0.48 953.26 27.90 1,888.67 245.58 2,830.78 0.67 647.03 14.33 981.63 4,703.94 1,445.65 4,703.94 1,548.93 4,703.94 1,746.57 4,703.94 2,069.01 2.46 821.25 4,703.94 1,266.87 4,703.94 1,643.56 4,703.94 2,037.22 7 2,955.82

cpu*

cpu

f

MNS-TS MU % ML %

388.82 9.02 386.18 1.46 0.00 376.90 9.47 376.78 1.23 0.23 477.69 26.66 474.57 6.04 0.00 669.49 56.19 670.13 6.00 1.19 616.35 59.45 609.81 14.82 0.18 621.13 59.71 616.90 26.24 0.80 623.91 57.30 615.98 21.48 0.33 2,911.51 711.73 2,852.56 55.39 -9.64 2,755.28 720.16 2,625.04 59.81 -6.10 2,616.02 730.94 2,539.60 48.61 -6.10 2,598.56 727.67 2,508.45 85.04 -11.82 955.80 32.21 953.25 1.21 0.00 1,899.98 220.17 1,885.80 9.18 0.12 2,854.81 678.28 2,818.34 6.65 0.00 647.03 8.56 647.03 0.56 0.00 985.69 55.47 979.82 9.07 0.10 1,458.75 197.50 1,426.95 28.35 0.24 1,586.67 240.42 1,517.80 14.09 -0.36 1,755.26 669.27 1,758.51 95.94 -3.25 2,085.62 1,154.75 2,009.97 68.08 -9.54 824.09 26.59 821.25 1.18 0.00 1,277.79 164.44 1,266.17 5.63 0.93 1,666.02 508.61 1,615.42 17.35 -2.65 2,076.08 1,137.37 2,013.86 175.46 -12.25 2

0.00 0.23 0.00 8.60 1.71 4.25 5.69 19.00 16.05 13.82 15.07 0.00 0.12 0.00 0.00 0.10 0.41 3.35 9.46 5.46 0.00 1.10 4.34 6.88

Mf ∗ %

Mf %

0.00 0.23 0.00 0.47 -0.41 0.32 0.12 -2.10 -4.14 -1.47 -2.08 0.00 -0.15 -0.44 0.00 -0.18 -1.29 -2.01 0.68 -2.85 0.00 -0.06 -1.71 -1.15

-0.68 -0.03 -0.65 0.10 -1.06 -0.68 -1.27 -2.02 -4.73 -2.92 -3.47 -0.27 -0.75 -1.28 0.00 -0.60 -2.18 -4.34 0.19 -3.63 -0.34 -0.91 -3.04 -3.00

14 344.25

31.62

-2.40% 4.82% -0.76% -1.56%

We compare the performance of our method with the exact method (MIP) proposed by Lalla-Ruiz et al. (2015) and with the Hybrid Genetic algorithm (HGA) proposed by Liu et al. (2014). Table 2 displays the results for the MDOVRP. The first three columns characterize the instances with the name of the instance, the number of customers n and the number of depots m. The three following columns are the results of MIP (upper and lower bounds denoted by UB and LB respectively and the scaled computational time cpu*). For HGA, we present the best value, f ∗ , and the average value, f, reached for each instance. Also we report the scaled computational time spent by HGA. The last six columns report the information about the MNS-TS. We first report the value of the objective function and the computational time spent for each instance. We also report the improvement of MNS-TS over the upper and the distance to the lower bound (MU and ML ) and the improvement of MNS-TS over the best value and average value of HGA (Mf ∗ and Mf ). At the bottom of Table 2, we summarize the main results of this numerical comparison such as the number of best solutions returned by each method, the average computational time and the average improvement brought by MNS-TS over its competitors. Again, the best solution values are displayed in bold. 9.3. Statistical test The results of the metaheuristics assessed over both set of instances are mutually independent. Thus we applied the non parametric statistical test for the solution value and computational time for separately comparing (univariate model (Chiarandini et al., 2007)) the performance in terms of solution quality and computational time. The goal of the statistical test is to detect significant differences between the performance of different approaches. Results are detailed in Appendix A. The statistical tests demonstrate that MNS-TS is the fastest approach while HGA and MIP are equivalent in terms of computational time for this set of instances. Moreover, MNS-TS is better than the average performance of HGA but there is not a significant difference in the performance in terms of the best solution value. From the statistical test and visual analysis, we can conclude that MNS-TS outperforms the state-of-the-art method for the OVRP and MDOVRP in terms computational time. And it outperforms the state-of-the-art methods for the OVRP in terms of solution value while it remains similar to the MIP for the MDOVRP. 23

10. Conclusions The success of the MNS-TS is due to the successive exploration of neighborhoods. This allows the search to obtain a good balance between intensification and diversification. In fact, each neighborhood ensures a good harmony between the intensity and diversity of the search. The composed neighborhoods ensure the intensity by exploring several neighborhoods and choose the best move. While the diversity is guaranteed by the tabu search which allows us to accept bad solutions and so exploring new parts of the solution space. The different neighborhoods are presented in a unified view thanks to the Path Moves and the Ejection Chains. In this paper, we propose a Multiple Neighborhood Search hybridized with a Tabu Search to address the Multiple Depot Open Vehicle Routing Problem. This approach have been tested on the mainstream instances of the literature for the OVRP and MDOVRP. From the statistical tests, we can conclude that our procedure outperforms the state-of-the-art methods for the OVRP and MDOVRP in terms of computational time. MNS-TS outperforms the approaches devised for the OVRP in terms of solution value. For the MDOVRP, the performance MNS-TS and MIP in terms of the solution value are comparable. References References Benlic, U., J-K. Hao. 2011. An effective multilevel tabu search approach for balanced graph partitioning. Computers & Operations Research 38(7) 1066–1075. Bodin, L., B. Golden, A. Assad, M Ball. 1983. Routing and scheduling of vehicles and crews: The state of the art. Computers & Operations Research 10(2) 63–211. Ceschia, S., L. Gaspero, A. Schaerf. 2011. Tabu search techniques for the heterogeneous vehicle routing problem with time windows and carrierdependent costs. Journal of Scheduling 14(6) 601–615. Chao, I-M., B.L. Golden, E. Wasil. 1993. A new heuristic for the multi-depot vehicle routing problem that improves upon best-known solutions. Am. J. Math. Manage. Sci. 13(3-4) 371–406. 24

Chiarandini, M., L. Paquete, M. Preuss, E. Ridge. 2007. Experiments on metaheuristics: Methodological overview and open issues. Tech. Rep. DMF-2007-03-003, The Danish Mathematical Society, Denmark. Christofides, N., S. Eilon. 1969. An algorithm for the vehicle-dispatching problem. OR 20(3) 309–318. Christofides, N., A. Mingozzi, P. Toth. 1979. Combinatorial optimization, chap. The vehicle routing problem. John Wiley, Chichester, 315–338. Conover, W.J. 1999. Practical nonparametric statistic. 3rd ed. Wiley, New York, USA. Cordeau, J-F., M. Gendreau, G. Laporte. 1997a. A tabu search heuristic for periodic and multi-depot vehicle routing problems. Networks 30(2) 105–119. Cordeau, J-F., M. Gendreau, G. Laporte. 1997b. A tabu search heuristic for periodic and multi-depot vehicle routing problems. Networks 30(2) 105–119. Dantzig, G. B., J. H. Ramser. 1959. The truck dispatching problem. Management Science 6(1) 80–91. Derigs, U., K. Reuter. 2009. A simple and efficient tabu search heuristic for solving the open vehicle routing problem. The Journal of the Operational Research Society 60(12) 1658–1669. Fisher, M.L., R. Jaikumar. 1981. A generalized assignment heuristic for vehicle routing. Networks 11(2) 109–124. Fleszar, K., I.H. Osman, K.S. Hindi. 2009. A variable neighbourhood search algorithm for the open vehicle routing problem. European Journal of Operational Research 195(3) 803–809. Friedman, M. 1937. The use of ranks to avoid the assumption of normality implicit in the analysis of variance. Journal of the American Statistical Association 32 675–701. Gaspero, L., A. Schaerf. 2006. Neighborhood portfolio approach for local search applied to timetabling problems. Journal of Mathematical Modelling and Algorithms 5(1) 65–89. 25

Gaspero, L. Di, A. Schaerf. 2003. Multi-neighbourhood local search with application to course timetabling. Proc. of the 4th Int. Conf. on the Practice and Theory of Automated Timetabling (PATAT-2002), number 2740 in Lecture Notes in Computer Science. Springer-Verlag, 262–275. Glover, F., M. Laguna. 1997. Tabu Search. Kluwer Academic Publisher, Dordrecht, The Netherlands. Jin, J., T.G. Crainic, A. Lokketangen. 2012. A parallel multi-neighborhood cooperative tabu search for capacitated vehicle routing problems. European Journal of Operational Research 222(3) 441–451. Lalla-Ruiz, E., C. Exp´osito-Izquierdo, S. Taheripour, Voss S. 2015. An improved formulation for the multi-depot open vehicle routing problem. Operations Research & Decision Theory, OR Spectrum . Lau, H.C.W., T. Chan, W.T. Tsui, W.K. Pang. 2010. Application of genetic algorithms to solve the multidepot vehicle routing problem. Automation Science and Engineering, IEEE Transactions on 7(2) 383–392. Letchford, A.N., J. Lysgaard, R.W. Eglese. 2007. A branch-and-cut algorithm for the capacitated open vehicle routing problem. Journal of the Operational Research Society 58(12) 1642–1651. Li, F., B. Golden, E. Wasil. 2007. The open vehicle routing problem: Algorithms, large-scale test problems, and computational results. Computers & Operations Research 34(10) 2918–2930. Li, X., S. Leung, P. Tian. 2012. A multistart adaptive memory-based tabu search algorithm for the heterogeneous fixed fleet open vehicle routing problem. Expert Systems with Applications 39(1) 365–374. Lim, A., F. Wang. 2005. Multi-depot vehicle routing problem: a one-stage approach. Automation Science and Engineering, IEEE Transactions on 2(4) 397–402. Lin, S., B. Kernighan. 1973. An effective heuristic algorithm for the traveling salesman problem. Operations Research 21 498–516. Liu, R., Z. Jiang, N. Geng. 2014. A hybrid genetic algorithm for the multidepot open vehicle routing problem. OR Spectrum 36(2) 401–421. 26

MirHassani, S.A., N. Abolghasemi. 2011. A particle swarm optimization algorithm for open vehicle routing problem. Expert Systems with Applications 38(9) 11547–11551. Pessoa, A., M.P. de Arag˜ao, E. Uchoa. 2008. Robust branch-cut-and-price algorithms for vehicle routing problems. B. Golden, S. Raghavan, E. Wasil, eds., The Vehicle Routing Problem: Latest Advances and New Challenges, Operations Research/Computer Science Interfaces, vol. 43. Springer US, 297–325. Pratt, J.W. 1959. Remarks on zeros and ties in the wilcoxon signed rank procedures. Journal of the American Statistical Association 54(287) 655– 667. Reinholz, A., H. Schneider. 2013. A hybrid (1+1)-evolutionary strategy for the open vehicle routing problem. L. Di Gaspero, A. Schaerf, T. St¨ utzle, eds., Advances in Metaheuristics, Operations Research/Computer Science Interfaces Series, vol. 53. Springer New York, 127–141. Renaud, J., G. Laporte, F. Boctor. 1996. A tabu search heuristic for the multi-depot vehicle routing problem. Computers & Operations Research 23(3) 229–235. Repoussis, P.P., C.D. Tarantilis, O. Br¨aysy, G. Ioannou. 2010. A hybrid evolution strategy for the open vehicle routing problem. Computers & Operations Research 37(3) 443–455. Repoussis, P.P., C.D. Tarantilis, G. Ioannou. 2007. The open vehicle routing problem with time windows. The Journal of the Operational Research Society 58(3) 355–367. Salari, M., P. Toth, A. Tramontani. 2010. An ILP improvement procedure for the open vehicle routing problem. Computers & Operations Research 37(12) 2106–2120. Sariklis, D., S. Powell. 2000. A heuristic method for the open vehicle routing problem. The Journal of the Operational Research Society 51(5) 564–573. Schrage, L. 1981. Formulation and structure of more complex/realistic routing and scheduling problems. Networks 11(2) 229–232. 27

Sevaux, M., A. Rossi, M. Soto, A. Duarte, R. Mart´ı. 2014. Grasp with ejection chains for the dynamic memory allocation in embedded systems. Soft Computing 18(8) 1515–1527. S¨orensen, K., M. Sevaux, P. Schittekat. 2008. ”multiple neighbourhood” search in commercial vrp packages: Evolving towards self-adaptive methods. Carlos Cotta, Marc Sevaux, Kenneth Srensen, eds., Adaptive and Multilevel Metaheuristics, Studies in Computational Intelligence, vol. 136. Springer Berlin Heidelberg, 239–253. Soto, M., A. Rossi, M. Sevaux. 2013. Iterative approaches for a dynamic memory allocation problem in embedded systems. European Journal of Operational Research 231(1) 34–42. Soto, M., A. Rossi, M. Sevaux. 2015. A multiple neighborhood search for dynamic memory allocation in embedded systems. Journal of Heuristics 21(6) 719–749. Soto, M., M. Sevaux, A. Rossi. 2014. Optimisation du transport de personnes handicap´ees. ROADEF, 15`eme congr`es de la Soci´et´e Fran¸caise de Recherche Op´erationnelle et d’Aide a` la D´ecision. Bordeaux, France. Tarantilis, C.D., G. Ioannou, C.T. Kiranoudis, G.P. Prastacos. 2005. Solving the open vehicle routeing problem via a single parameter metaheuristic algorithm. The Journal of the Operational Research Society 56(5) 588– 596. Tarantilis, C.D., C.T. Kiranoudis. 2002. Distribution of fresh meat. Journal of Food Engineering 51(1) 85–91. Tarantilis, C.D, C.T Kiranoudis, V.S Vassiliadis. 2004. A threshold accepting metaheuristic for the heterogeneous fixed fleet vehicle routing problem. European Journal of Operational Research 152(1) 148–158. Wilcoxon, F. 1945. Individual Comparisons by Ranking Methods. Biometrics Bulletin 1(6) 80–83. Wren, A., A. Holliday. 1972. Computer scheduling of vehicles from one or more depots to a number of delivery points. Operational Research Quarterly (1970-1977) 23(3) 333–344. 28

Yu, S., C. Ding, K. Zhu. 2011. A hybrid ga-ts algorithm for open vehicle routing optimization of coal mines material. Expert Systems with Applications 38(8) 10568–10573. Zachariadis, E., C. Kiranoudis. 2010. An open vehicle routing problem metaheuristic for examining wide solution neighborhoods. Computers & Operations Research 37(4) 712–723.

29

Appendix A. Statistical tests To statistically evaluate the result displayed by Table 1 on the OVRP instances, we apply the well-known Wilcoxon test (Wilcoxon, 1945) because there are two methods to compare. While, we apply the non parametric Friedman test (Friedman, 1937) to evaluate the result of Table 2 for the MDOVRP instances. Appendix A.1. Statistically assessing OVRP results. There are two methods to calculate the signed-rank test. The first one is the classic Wilcoxon test (Wilcoxon, 1945), which discards any tied data and then calculates the signed ranks. The second method incorporates tied values in the ranking procedure. It was proposed by Pratt (1959). We have computed the two non parametric tests. The p-value for the classic Wilcoxon test is 0.0323 and for the Wilcoxon test proposed by Pratt is 0.0235. Both p-values indicate that there is a difference between the two methods using the typical significance level of 0.05 as the threshold between rejecting or not rejecting the null hypothesis. For the time, the p-value for the Wilcoxon test is 0.0163. This result means that there is a significant difference between the computational time of MNS-TS and the best metaheuristics. We can conclude that MNS-TS is better in terms of the objective function and computational time than the best metaheuristics from the state-of-the-art (ILP Improvement, ABHC and the Hybrid strategy). Appendix A.2. Statistically assessing MDOVRP results. For HGA, we know the best value and the average value for the solution value. First, we apply the Friedman test using the best value reached by HGA, by MIP and by MNS-TS. The p-value for this data is 0.1530 which means there is not a significant difference in terms of objective function for these three methods. Second, we apply the Friedman test using the average value for HGA, MIP and MNS-TS. The p-value returned by the test is 7.98 × 10−4 . Then, the null hypothesis is rejected at the level of significance 0.05. Therefore, there exists at least one approach whose performance is different from at least one of the others. Finally, we use a Post-Hoc paired comparison to know which metaheuristics are different (Conover, 1999). Table A.3 summarizes p-values for the paired comparisons for the cost. The bold values means the approaches are different in terms of solution value. 30

Table A.3: p-value for Post-Hoc paired comparison for the solution value

MIP

MNS-TS

HGA 2.570 × 10−2 MIP -

8.288 × 10−4 5.530 × 10−1

The results returned by the Post-Hoc comparison show that the solution value reached by MNS-TS and MIP are different from the average quality of the solutions reached by HGA. Moreover, there does not exist a significant difference between the performance of MNS-TS and MIP in terms of the solution value. The p-value for the Friedman test for the computational time is 2.379 × −7 10 , then the null hypothesis is rejected at the level of significance 0.05. Then, we can assume that there exists at least one approach whose performance in terms of CPU time is different from at least one of the others. Table A.4 summarizes p-values for the paired Post-Hoc comparison for the computational time. Once again, bold values mean the approaches are different in terms of CPU time. Table A.4: p-value for Post-Hoc paired comparison for the CPU time

HGA MIP

MIP

MNS-TS

4.805 × 10−1 -

8.725 × 10−5 7.537 × 10−7

This test means that there is a significant difference between the computational time of MNS-TS and one two other methods. However, there does not exist a significant difference between the CPU time of HGA and MIP. MNS-TS is the fastest approach while HGA and MIP are equivalent in terms of computational time for this set of instances. Summarizing the experimentation above, we can say that MNS-TS has the best performance in terms of the computational time while in terms of the solution value it is better than the average performance of the HGA but there is not a significant difference in the performance in terms of the best solution value. From the statistical test and visual analysis, we can therefore conclude that MNS-TS outperforms the state-of-the-art method for the OVRP and MDOVRP in terms computational time. And it outperforms the state-of31

the-art methods for the OVRP in terms of solution value while it remains similar to the MIP for the MDOVRP.

32

Multiple neighborhood search, tabu search and ejection chains for the multi-depot open vehicle routing problem María Soto, Marc Sevaux1 Université de Bretagne-Sud Lab-STICC – UMR 6285 – CNRS 2 rue de Saint Maudé F-56321 Lorient, France [email protected] [email protected] André Rossi Université d’Angers LERIA – EA 2645 – CNRS UFR Sciences 2 Boulevard Lavoisier 49045 Angers, France [email protected] Andreas Reinholz DLR, Institute of Solar Research Linder Höhe D-51147 Köln, Germany [email protected]

1 Corresponding

author