The cumulative capacitated vehicle routing problem: New formulations and iterated greedy algorithms

Accepted Manuscript

The Cumulative Capacitated Vehicle Routing Problem: New Formulations and Iterated Greedy algorithms ´ Samuel Nucamendi-Guillen, ´ Francisco Angel-Bello, Iris Mart´ınez-Salazar, Alvaro E. Cordero-Franco PII: DOI: Reference:

S0957-4174(18)30439-1 10.1016/j.eswa.2018.07.025 ESWA 12078

To appear in:

Expert Systems With Applications

Received date: Revised date: Accepted date:

14 November 2017 9 July 2018 10 July 2018

´ Please cite this article as: Samuel Nucamendi-Guillen, Iris Mart´ınez-Salazar, ´ Francisco Angel-Bello, Alvaro E. Cordero-Franco, The Cumulative Capacitated Vehicle Routing Problem: New Formulations and Iterated Greedy algorithms, Expert Systems With Applications (2018), doi: 10.1016/j.eswa.2018.07.025

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.

ACCEPTED MANUSCRIPT

Highlights

• Two Iterated greedy algorithms are developed.

CR IP T

• Two models for the Cumulative Capacitated Vehicle Routing Problem are proposed.

• The algorithms found new best known results for small instances.

• The algorithms found the best known results for medium size instances in less time.

AC

CE

PT

ED

M

AN US

• The algorithms found competitive results against the literature for large instances.

1

ACCEPTED MANUSCRIPT

The Cumulative Capacitated Vehicle Routing Problem: New Formulations and Iterated Greedy algorithms

CR IP T

b,∗ ´ Samuel Nucamendi-Guill´ena , Francisco Angel-Bello , Iris c d Mart´ınez-Salazar , Alvaro E. Cordero-Franco a

AN US

Universidad Panamericana Campus Guadalajara, Facultad de Ingenier´ıa, Zapopan, Jalisco CP 45010, M´exico b Tecnologico de Monterrey, Escuela de Ingenier´ıa y Ciencias. Ave. Eugenio Garza Sada 2501, Monterrey, Nuevo Le´ on CP 64849, M´exico c Universidad Aut´ onoma de Nuevo Le´ on, Graduated Program in Systems Engineering, Av. Universidad s/n. San Nicol´ as de los Garza, Nuevo Le´ on CP 66451, M´exico d Universidad Aut´ onoma de Nuevo Le´ on, Facultad de Ciencias F´ısico Matem´ aticas, Av. Universidad s/n. San Nicol´ as de los Garza, Nuevo Le´ on CP 66451, M´exico

Abstract

AC

CE

PT

ED

M

In this work, we address the Cumulative Capacitated Vehicle Routing Problem (CCVRP), a variant of the classical CVRP, which aims to minimize the sum of the arrival times to customers instead of the total traveled distance. This problem is relevant due to its applications in fields such as emergency logistics, transportation, wireless network and computing, among others. This paper presents the first two tractable integer formulations capable of solving instances with up to 44 nodes and two versions of an Iterated Greedy procedure for dealing with larger instances. Experimental results show that both formulations obtain optimal solutions in a reasonable amount of computational time. Regarding the two metaheuristic procedures, they were able to reach optimal and best known solutions and, in some cases, outperform the bounds obtained by exact methods for tested instances with up to 199 customers. Using larger instances, the proposed metaheuristics showed competitive results when compared to the best known values reported in the literature, with a significant reduction in computational time. ∗

Corresponding author Email addresses: [email protected] (Samuel Nucamendi-Guill´en), ´ [email protected] (Francisco Angel-Bello ), [email protected] (Iris Mart´ınez-Salazar), [email protected] (Alvaro E. Cordero-Franco)

Preprint submitted to Expert Systems with Applications

July 10, 2018

ACCEPTED MANUSCRIPT

Keywords: Vehicle routing, Integer programming, Metaheuristics, Latency, Cumulative costs

CR IP T

1. Introduction

AC

CE

PT

ED

M

AN US

The Cumulative Capacitated Vehicle Routing Problem (CCVRP) is a variation of the classical Capacitated Vehicle Routing Problem (CVRP) (Dantzig and Ramser, 1959), where, instead of minimizing the total travel cost or distance, its objective is to minimize the customers’ waiting time. This problem is categorized as a ”customer-centric” routing problem. Practical applications of the CCVRP can be found in several areas like disaster logistics and maintenance operations, or any situation where the arrival time to customers is crucial. CCVRP is a generalization of the Minimum Latency Problem (MLP) (Blum et al., 1994), where the former has gained special interest in recent years due to its complexity, variety and relevance of its applications. Some of the most important works in the MLP field are: reformulations of mathematical models (M´endez-D´ıaz et al., 2008, Ezzine et al., 2010, Angel-Bello et al., 2013, Bjeli´c et al., 2013); exact algorithms (M´endez-D´ıaz et al., 2008, Abeledo et al., 2013); and heuristic approaches (Salehipour et al., 2008, 2011, Ngueveu et al., 2010, Silva et al., 2012, Mladenovi´c et al., 2013, Bjeli´c et al., 2013). The CCVRP was first introduced by Ngueveu et al. (2010), who presented a mathematical model based on a vehicle routing problem with time windows (Kallehauge et al., 2006) and developed a memetic algorithm (MA) as a solution method. To validate their approach, the algorithm was evaluated using the Traveling Repairman Problem (TRP) instances (Salehipour et al., 2008) and modified CVRP instances (Christofides et al., 1979). Ribeiro and Laporte (2012) developed an adaptive large neighborhood search (ALNS) heuristic for the CCVRP. Their results outperformed the MA proposed by Ngueveu et al. (2010) in terms of computational time and quality of the solution. Later, Ke and Feng (2013) created a two-phase metaheuristic for the CCVRP based on perturbation and local search operators to improvement of initial solutions. This metaheuristic obtained better results than the MA proposed by Ngueveu et al. (2010) and the ALNS proposed by Ribeiro and Laporte (2012). In 2014, Lysgaard and Wøhlk (2014) presented the first exact algorithm for the CCVRP based on branch-and-cut-and-price proce-

3

ACCEPTED MANUSCRIPT

PT

ED

M

AN US

CR IP T

dure (BCP). According to the authors, their algorithm is capable of solving instances with up to 150 customers in reasonable computational time. Mart´ınez-Salazar et al. (2014) introduced the multitrip CCVRP (mtCCVRP), a variant of the CCVRP where a single vehicle has to minimize the total latency with multiple trips. They proposed two mixed integer formulations capable of solving instances up to 25 nodes and developed a GRASP metaheuristic, where high-quality solutions were obtained. Rivera et al. (2015) developed another mixed integer linear program for this problem and a multi-start iterated local search (MS-ILS) algorithm. They assessed their MS-ILS algorithm using CCVRP instances and compared its results with existing algorithms in the literature (Ngueveu et al., 2010, Ribeiro and Laporte, 2012, Ke and Feng, 2013). According to the authors, the MS-ILS found 7 out of 26 best known solutions, and for the rest of the instances, it obtained solution values with a deviation, on average, of 1.01%. Recently, Cinar et al. (2016) studied the cumulative vehicle routing problem with limited duration. In this case, the cumulative cost was considered as a step function of the load of the vehicle while customers are visited. The authors presented a mathematical model as an extension of the cumulative routing problem presented by Kara et al. (2008). A variant of this problem, considering stochastic demands, was studied by Gaur et al. (2016). In 2017, Sze et al. (2017) proposed a two-stage adaptive variable neighborhood search (AVNS) algorithm, capable of obtaining the best-known solutions so far, compared to the algorithms in the literature. In this paper we study the CCVRP, as previously described, where the main contributions are: • We propose two tractable mixed integer formulations for the problem.

AC

CE

• We develop two variants of an iterated greedy algorithm, capable of reaching optimal solutions for small and medium size instances, and obtaining competitive results for larger instances, regarding the BKS, consuming less computational time than methods from the literature.

The remainder of this paper is organized as follows: In Section 2 we define the problem and present the proposed mathematical formulations. Section 3 describes the design of the metaheuristic algorithm and its different components. In Section 4 we discuss results of computational experiments. Finally, in Section 5 we present our conclusions, major findings and key guidelines for future research. 4

ACCEPTED MANUSCRIPT

2. Problem Formulation

AN US

CR IP T

2.1. Problem Definition The CCVRP is defined formally through an undirected graph G = (V, A) where V = {0, 1, 2, . . . , n} is the node set and A is the arc set. Node 0 corresponds to the depot, and the remaining nodes are in the customers set V 0 = {1, 2, . . . , n}. Each arc (i, j) ∈ A has an associated travel time cij , whereas each node j ∈ V 0 has an associated service time sj and a known demand dj . There is a fixed fleet of k vehicles, each with capacity Q. The goal is to find a set of k disjoint tours, having in common only the node 0 in the first position, satisfying the demand of every customer while minimizing the sum of the customers waiting times (also called latencies), without exceeding the vehicles capacity.

AC

CE

PT

ED

M

2.2. Mathematical Formulations We propose two mixed integer formulations for the CCVRP. The first one is developed from a single-flow formulation for the Multiple Traveling Salesman Problem (Gavish and Graves, 1978) and its modification to the Multiple Minimum Latency Problem (mMLP) (Angel-Bello et al., 2017). The second one is based on time-dependent mathematical models for the k-Travelling Repairmen Problem (k-TRP) (Nucamendi-Guill´en et al., 2016) and for the mMLP (Angel-Bello et al., 2017). Both formulations are generalized to the CCVRP by adding the capacity constraints. In the first model, let xij be a binary variable equal to 1 if and only if arc (i, j) is used in any route, yij is an integer variable that represents the number of remaining nodes in a route after visiting customer i when xij = 1, and it is equal to zero when xij = 0. Also, let vij be a variable that denotes the sum of the remaining demands after node i in a route when xij = 1, and it is equal to zero when xij = 0. Thus, CCVRP is formulated as:

min z =

n X n X i=0 j=1 j6=i

5

cij ∗ yij

(1)

ACCEPTED MANUSCRIPT

subject to: n X

x0i = k

(2)

j=0 j6=i n X

xi0 = k

(3)

xij = 1

i = 1, 2, . . . , n

xij = 1

j = 1, 2, . . . , n

i=0 i6=j

y0i +

X j=1 j6=i

yji −

X

yij = 1

j=1 j6=i

i, j = 1, 2, . . . , n;

yij ≤ (n − k)xij

i, j = 1, 2, . . . , n;

v0j ≤ Q ∗ x0j

ED

vij ≥ dj xij

i = 1, 2, . . . , n

M

yij ≥ 2xij − xj0

(5)

i = 1, 2, . . . , n

yij ≥ xij

y0i ≤ (n − k + 1)x0i

(4)

AN US

i=1 n X

CR IP T

i=1

n X

i 6= j

(7)

i 6= j

(9)

i = 0, 1, . . . , n; j = 1, 2, . . . , n;

i, j = 1, 2, . . . , n; j = 1, 2, . . . , n

(6)

i 6= j

i 6= j

(10) (11) (12)

vij ≤ (Q − di )xij X X v0i + vji − vij = di

i, j = 1, 2, . . . , n; i = 1, 2, . . . , n

(14)

y0i ≥ 0

i = 1, 2, . . . , n

(15)

i, j = 1, 2, . . . , n;

vij ≥ 0

i, j = 1, 2, . . . , n;

CE

yij ≥ 0

v0i ≥ 0

AC

xij ∈ {0, 1}

(13)

j=1 j6=i

PT

j=1 j6=i

i 6= j

(8)

i = 1, 2, . . . , n i, j = 0, 1, . . . , n;

i 6= j

(16)

i 6= j

(18)

j 6= i

(17) (19)

In this formulation, the objective function (1) aims to minimize the total latency of the system. Constraints (2) and (3) ensure that exactly k arcs leave and return to the depot. Constraints (4) and (5) impose that exactly one arc enters and leaves each node associated with a customer. Constraints 6

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 1: Representation of the multilevel network

AC

CE

PT

ED

M

(7) and (10) are valid inequalities proposed by Godinho et al. (2009) for the unit demand VRP. Constraints (8) and (9) force variables yij to be zero when xij = 0. Constraints (6), (8) and (9) avoid sub-tours and allow us to calculate a customer position on the routes. Constraints (11) are derived from a generalization of restrictions (7). Constraints (12) and (13) force variables vij to be zero when xij = 0, and, in conjunction with (14), they estimate the load per vehicle. Finally, constraints (15) - (19) establish the nature of the variables. This formulation consists of n2 N − n binary variables, n2 + nN real variables, and 4n2 + 2n + 4 constraints, where N = n − k + 1. To better understand the second proposed mathematical model, we use the multilevel network presented in Figure 1 (Angel-Bello et al., 2017). In this network, N = n − k + 1 denotes the maximum number of customers that might belong to a route. Level 1 contains a copy of every node associated with a customer, while levels 2 to N have a copy for every node. Finally, level N+1 consists of a copy of the depot. Each path must start with a copy of the depot and must visit nodes at lower levels, ending at a node in level 1. If the starting node is at level m, then, the corresponding path contains m − 1 nodes associated with clients. Two paths cannot include a node associated with the same client, even in different levels. We say that a node is active in a certain level when it is visited by a path at that same level. (r) For this model, the following decision variables are defined. Let si be a binary variable equal to 1 if and only if node i is active at level r in the (r) multilevel network and let tij be a binary variable equal to 1 if and only if arc (i, j) is used to link the node i in level r + 1 with node j in level r in the

7

ACCEPTED MANUSCRIPT

(r)

min z =

n X

c0j

j=1

N X

(1)

=k

si

i=1 N X n X

(r)

PT j=1 j6=i

n X

CE

+

i = 1, 2, . . . , n

ED

(N )

t0i = si n X (r) (r+1) tij = si (r) t0j

(r)

r=1

i=1 j=1 j6=i

(r)

rtij

t0j = k

r=1 j=1

(N )

cij

N X

AN US

=1

n X n X

M

r=1 n X

(r)

si

+

r=1

subject to: N X

(r) rt0j

CR IP T

multilevel network. Note that si = 1 means that there are (r -1) customers after customer i in the path where customer i belongs, and trij = 1 means that the arc (i, j) belongs to a path in which, after customer i, the number of remaining customers is exactly r. Finally, let vij be a variable that denotes P (r) the sum of node demands after node i in a path when N r=1 tij = 1 and it is equal to 0 otherwise. The resulting formulation is described below:

(r)

tij = sj

(20)

(21) (22)

(23)

i = 1, 2, . . . , n

(24)

i = 1, 2, . . . , n; r = 1, 2 . . . , N − 1

(25)

j = 1, 2, . . . , n; r = 1, 2 . . . , N

(26)

AC

i=1 i6=j

8

ACCEPTED MANUSCRIPT

N X

(r)

t0j

j=1 j6=i

(r)

si

vji −

N X

N −1 X r=1 n X

(r)

tij

i, j = 1, 2, . . . , n;

vij = di

i = 1, 2, . . . , n

j=1 j6=i

(r)

tij

i, j = 1, 2, . . . , n;

r=1

≥0

v0j ≥ 0 vij ≥ 0 (r)

(29)

i 6= j

(30)

(31)

j = 1, 2, . . . , n

(32)

j 6= i

j = 1, 2, . . . , n; r = 1, 2, . . . , N

(r)

i, j = 1, . . . , n;

j 6= i;

(33) (34) (35)

r = 1, 2, . . . , N − 1

M

tij ∈ {0, 1}

(28)

i = 1, 2, . . . , n : r = 1, 2, . . . , N

i, j = 1, 2, . . . , n;

t0j ∈ {0, 1}

i 6= j

CR IP T

n X

vij ≥ dj

(27)

r=1

vij ≤ (Q − di ) v0i +

j = 1, 2, . . . , n

AN US

v0j ≤ Q

AC

CE

PT

ED

In this second formulation, the objective function (20) minimizes the sum of the latencies for every customer. Constraints (21) ensure that node i is active at only one level. Constraints (22) ensure that exactly k edges will arrive at level 1, while constraint (23) ensures that exactly k edges will depart from node 0 copies. Constraints (24) impose that node 0 at level N + 1 must be linked with nodes that are active at level N . Constraints (25) and (26) ensure that arcs departing from active nodes in level r + 1 can arrive only to active nodes in level r. Constraints (27) and (28) force variables to be zero when the arc (i, j) is not used in the solution. Constraints (27) guarantee that vehicle capacity is not exceeded. Constraints (27), in conjunction with restrictions (28), (29) and (30), are the sub-tour elimination constraints and allow us to estimate the load of each vehicle. Finally, constraints from (31) to (35) establish the nature of the variables. Note that because of binary (r) (r) variables tij , constraints (24), (25) and (26) force variables si to take binary values.

9

ACCEPTED MANUSCRIPT

CR IP T

This model has n2 N −n2 +n binary variables, n2 +nN real variables, and 3n2 − 2nr + 3n + 2 constraints. It can be seen that the second formulation requires fewer binary variables with respect to the first model and fewer constraints. 3. Metaheuristic approach

AC

CE

PT

ED

M

AN US

In this section, we describe two versions of a metaheuristic algorithm using iterated greedy (IG) strategies, which has shown to be very efficient for solving the uncapacitated problem variant (Nucamendi-Guill´en et al., 2016). Each algorithm version is composed of a constructive procedure, an improvement procedure and a destructive procedure. The improvement and destructive procedures are the same for both versions, while the constructive strategies are different in each of them. In our algorithm, a solution is represented by S = {S1 , S2 , · · · , Sk }, where Sl represents a route with nl customers (l = 1, 2, · · · , k). The first version of the constructive procedure is based on a parallelbuilding constructive strategy, while the other version is based on a clustering strategy. The destructive procedure receives a feasible solution and applies a controlled randomized strategy to eliminate customers from routes. As an output of the destructive procedure we obtain a set of unrouted customers Sa and a partial destroyed solution Sp = {Sp1 , Sp2 , · · · , Spk } with nl , (l = 1, 2, · · · , k) customers per route. Concerning the improvement procedure, it is based on a variable neighborhood descent with random neighborhood ordering (RVND) mechanism. A pseudo-code of the general IG algorithm is shown in Figure 2. The algorithm inputs are: the number of customers (n), number of vehicles (k ), set of unrouted customers (Sa) and cost matrix (C ). Step 1 initializes the partial solution (Sp), the number of customers per route (nl ), the iteration counter (iter ) and the maximum number of iterations without solution improvement (MaxIter ). Due to vehicle capacity constraints, the applied constructive procedure could create a solution with unrouted customers (infeasible solution). For this reason, steps 3 and 4 (later steps 9 and 10) are repeated until a feasible solution (S ), with objective value (L), is found. The improvement procedure is applied over every single feasible solution, saving the newly obtained solution bS and its objective value bL (step 6). Then, the algorithm goes into a loop where an iterative process

10

ACCEPTED MANUSCRIPT

4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

CR IP T

3

AN US

2

M

1

// Iterated Greedy Procedure Data: n, k, Sa = {1, 2, · · · , n}, C Initialization: Spl = {0}, nl ← 0(l = 1, 2, · · · , k), iter ← 0, M axIter ← T repeat (S, L) ← ConstructiveP rocedure(n, k, Sa, Sp, nl , C) if (Sa 6= ∅) then (Sp, Sa, nl ) ← DestructiveP rocedure(S) until (Sa = ∅); (bS, bL) ← ImprovementP rocedure(S, L, C) while (iter < M axIter) do repeat (Sp, Sa, nl ) ← DestructiveP rocedure(bS) (S, L) ← ConstructiveP rocedure(n, k, Sa, Sp, nl , C) until (Sa = ∅) (S 0 , L0 ) ← ImprovementP rocedure(S, L, C) if (bL > L0 ) then (bS, bL) ← (S 0 , L0 ) iter ← 0 end iter ← iter + 1 end Result: (bS, bL)

ED

Figure 2: Pseudo-code for the general Iterated Greedy Procedure

CE

PT

of destruction-construction and improvement is applied, saving the best solution so far (steps 8-17). The stopping criterion is a predefined number of iterations without improvement (step 7). To avoid the undesired termination of positive progression, the iteration counter is reset every time a better solution is found (step 15). In next subsections, we present a detailed description of each procedure.

AC

3.1. Constructive Procedure For the metaheuristic procedure, two different constructive strategies were developed, one based on a Parallel Building (PRB) strategy and the other based on a Clustering strategy (CE). For the PRB strategy, we adapted the generalized regret criterion used by Potvin and Rousseau (1993). The second strategy (CE) consists of grouping customers into clusters. In this strategy, for each customer in Sa, a measure of proximity from each cluster 11

ACCEPTED MANUSCRIPT

CR IP T

is computed and the customer with the lowest value is selected to be inserted into its nearest cluster. Both constructive procedures require as input the number of customers (n), the number of vehicles (k ), the set of unrouted customers (Sa), the current partial solution (Sp), the cost matrix C, the number of customers per route (nl ) in Sp and the remaining vehicle capacities (RQl , l = 1, 2, · · · , k).

AC

CE

PT

ED

M

AN US

Parallel Route Building Strategy This constructive procedure takes the ideas from a parallel route constructive approach, originally proposed by Potvin and Rousseau (1993), for the Vehicle Routing Problem with Time Windows. The main idea of this approach is to use the generalized regret criterion to include unrouted customers in a solution under construction Sp, that is, the procedure calculates for each unrouted customer the k best feasible insertion points over all partial routes, computes the differences between the values of all insertion points and the value of the best insertion point and sums the values of these differences. Then, the customer with the greater sum is inserted into its best place in Sp. An outline of the PRB procedure is shown in Figure 3. The first time the constructive procedure is applied to build an initial solution, all routes in Sp are empty and RQl = Q (l = 1, 2, · · · , k). For successive constructions, partial empty routes could be obtained as a result of applying the destructive procedure. In any of these cases, the procedure selects, for empty routes, customers in Sa with the greatest cost value respect to the depot (step 3). Steps 6 to 11 correspond to our implementation of the generalized regret measure. In addition to the solution and its objective value, the procedure returns Sa to control the feasibility of the obtained solution. The value of an insertion point for a customer j is calculated as the increment that the objective function would have as result of inserting customer j at that position into the partial route. Since the searching process for the best insertion point in a sequence in problems with cumulative costs is very time-consuming, we developed a new efficient insertion strategy to achieve the same purpose. Given a route Sr with m customers, its total latency is calculated as L(Sr ) = mc[0][1] + (m − 1)c[1][2] + · · · + 2c[m−2][m−1] + c[m−1][m]

(36)

As a result, using expression (36), we developed the formula (37) that 12

ACCEPTED MANUSCRIPT

3 4 5 6 7 8 9 10 11 12 13

M

14

CR IP T

2

AN US

1

// PRB ConstructiveProcedure() Data: n, k, Sa, Sp, C, nl , RQl (l = 1, 2, · · · , k) while ((Sa 6= {∅}) and (there is enough capacity in any vehicle)) do if (there are empty routes in Sp) then Initialize them with customers j ∈ Sa that have greater values of C0j Update Sa and RQl end foreach costumer ∈ Sa do Determine the k best feasible insertion points over all partial routes Compute the differences between the values of all insertion points and the value of the best insertion point Sum the values of these differences end Insert the customer with the greater sum into its best place in Sp Update Sa and RQl end S ← Sp and compute L for the current solution S Result: (S, L) and Sa Figure 3: Outline of the PRB Constructive Procedure

CE

PT

ED

enables us to use the calculations made in the insertion point q to evaluate the q + 1 insertion point.  (m + 1)c0j + m(cj[1] − c0[1] ), if q = 1    Pq−1 c l=1 [l−1][l] + (m − q + 2)c[q−1]j (37) ∆rjq =  +(m − q + 1)(c − c ), if 2 ≤ q ≤ m j[q] [q−1][q]   Pm if q = m + 1 l=1 c[l−1][l] + c[m]j ,

AC

In formula (37), ∆rqj is the increment that would have the total latency of partial route Sr as a result of inserting client j at position q in route Sr . [q] denotes the client who occupies the position q in the route and m is the current number of customers in route r. Using the formula (37) in an iterative way starting from q = 1, it is guaranteed to perform O(1) elementary operations to evaluate each insertion point.

13

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

Clustering Strategy The developed clustering procedure follows the principles of centered clustering presented by Mulvey and Beck (1984). These principles help to reduce the dispersion among the customers in each cluster. This version generates as many clusters as vehicles are required, and it initializes and updates cluster’s central positions in the same way. The process to initialize clusters is the following: the first customer is randomly selected from Sa and is assigned as the center of the first cluster. Then, the most distant customer from the previously chosen one is selected as the center of the second cluster. For the third cluster, the farthest customer from any of the two centers is chosen, and so on. The process finishes after every cluster in Sp has a defined initial center. The coordinates of the initial C centers (aC r , br ) correspond to the coordinates of assigned customers, who, at the same time, are removed from the set of unrouted customers Sa. The assignment of the remaining (n-k) unrouted customers considers the proximity criterion concerning the centers of partial clusters. Every time the procedure adds a new customer to a cluster r, the coordinates of its center are updated using the expression (38): P P i∈Cr bi ∗ di i∈Cr ai ∗ di C C ; br = P (38) ar = P i∈Cr di i∈Cr di

AC

CE

PT

ED

where (ai , bi ) and di denote the coordinates and the demands, respectively, for the customers belonging to cluster r. In addition, each cluster has an associated limited capacity corresponding to vehicle capacity Q. Therefore, every time that a new customer is added to a cluster, its remaining capacity is updated by subtracting the demand of the added customer. Then, we say that a cluster is feasible for a given customer if the corresponding demand does not exceed the remaining cluster capacity. Now, we describe how the assignment of customers to clusters is performed. This strategy calculates, for each customer in Sa, the Euclidean distance from the center of each feasible cluster and inserts the customer with the shortest distance into its nearest cluster. We named this procedure as complete evaluation (CE) clustering procedure. Figure 4 shows the pseudocode for this procedure. The process of assigning customers to clusters stops when either Sa is empty or the procedure cannot insert any customer into any cluster due to 14

ACCEPTED MANUSCRIPT

3 4 5 6 7 8 9 10 11 12 13 14

CR IP T

2

AN US

1

// CE Clustering Procedure() Data: n, k, Sa, Sp, C, nl , RQl (l = 1, 2, · · · , k) while ((Sa 6= {∅}) and (there is enough capacity in any cluster)) do if (there are empty clusters) then Select initial centers for them Update Sa and RQl end forall the customer ∈ Sa do Compute the distances to the center of each feasible cluster end Insert the customer with the shortest distance into its closest cluster Update Sa and RQl Update cluster centres end Build a solution S by assigning customers in clusters to routes Compute L for the current solution S Result: (S, L) and Sa

M

Figure 4: Clustering procedure based on complete evaluation

PT

ED

capacity constraint. If we have a feasible solution, the route that the cluster follows is designed by routing customers in the same order they were selected when first included in the cluster. The resulting set of routes is defined as the initial solution (S ). Otherwise, the loop described in steps 3 and 4 from figure 2 is applied.

AC

CE

3.2. Destructive Procedure The destructive-constructive process implemented in our IG algorithm works as a diversification generator phase. Its purpose is to produce diverse initial solutions for the improvement phase. The proposed destructive procedure randomly eliminates a group of customers in S and assigns them to Sa. The partially destroyed solution S is then assigned to a partial solution Sp. In order to enhance the diversity among future initial solutions, we control the removing process by avoiding recently removed customers being removed again within a given number of future iterations (NIL). If a customer that occupies position p in S is selected, it is blocked for NIL successive iterations, and the blocked positions are saved in a list (LBP ). The values of the parameters number of customers to remove (ncr) and N IL are described in 15

ACCEPTED MANUSCRIPT

subsection 4.3. Figure 5 shows an outline of the destructive procedure.

2 3 4 5 6

CR IP T

1

// Destructive Procedure() Data: S, n, k, nl , Sa, ncr, LBP Compute NIL Randomly select ncr unblocked positions in S Remove from S the customers in the selected positions and add them to Sa Block the selected positions for N IL iterations in LBP Assign the partially destroyed solution S to Sp Update nl (l = 1, 2, · · · , k) Result: nl , Sp, Sa

AN US

Figure 5: Outline of the destructive procedure

PT

ED

M

3.3. Improvement Procedure The improvement procedure is based on a Variable Neighborhood Descent with Random Neighborhood Ordering (RVND) (Mladenovi´c and Hansen, 1997, Subramanian et al., 2010). The RVND operates on a set of structures of neighborhoods N = {N1 , N2 , · · · , Nr }. It randomly chooses a neighborhood from N to be explored. If a selected neighborhood fails to improve the current best solution, it is deleted from the set N, and RVND randomly selects another neighborhood from N to continue the search. When a better solution is found, the set N is restored and the search process is restarted. The RVND procedure ends when the set N is empty. In our implementation of the RVND procedure, we define the following five different neighborhoods:

CE

• (N1 ) Intra-route interchange of two customers: Two customers from a given route interchange their positions.

AC

• (N2 ) Intra-route insertion: A customer is removed from its current position and it is reinserted in another position towards the end of the route.

• (N3 ) Intra-route 2-opt: Two arcs are removed from a route and the two resulting subsequences are reconnected in a different way to obtain a new route.

16

ACCEPTED MANUSCRIPT

• (N4 ) Inter-route interchange of two customers: Two customers from two different routes interchange their positions.

CR IP T

• (N5 ) Inter-route insertion: A customer is removed from a route and it is inserted into another route.

AN US

When the objective function is to minimize the total latency, local moves impact the latency values of a large number of customers. As a consequence, direct cost evaluation of a simple move is, in general, significantly timeconsuming. For this reason, we developed effective strategies to explore the defined neighborhoods. Using these strategies, the number of operations to explore a neighborhood in the CCVRP is equivalent to the number of operations that are needed to explore the same neighborhood in the CVRP. In (N1 ) only two customers change their positions, then using the formula (36) it is easy to obtain the following expression (39) to evaluate the total latency variation in a route when customers [i] and [j] exchanged their positions

M

∆ij = (n − i + 1)(c[i−1][j] − c[i−1][i] ) + (n − i)(c[j][i+1] − c[i][i+1] ) +(n − j + 1)(c[j−1][i] − c[j−1][j] ) + (n − j)(c[i][j+1] − c[j][j+1] )

(39)

PT

ED

In (N2 ) the intra-route insertion movements are evaluated through swap movements of consecutive positions and the movement value is calculated in an incremental way. To evaluate the swap movements between customers in adjacent positions, we derive the expression (40) ∆i,i+1 = (n − i + 1)(c[i−1][i+1] − c[i−1][i] ) + (n − i)(c[i+1][i] − c[i][i+1] ) +(n − i − 1)(c[i][i+2] − c[i+1][i+2] )

(40)

AC

CE

In (N3 ) two non-adjacent arcs linking positions i − 1 and i and positions j −1 and j are deleted and the chain of customers in positions i, i+1, · · · , j − 1 is inverted. To evaluate the 2-opt moves, we follow the ideas presented in Ngueveu et al. (2010). In this kind of movement the customers in the positions up to i − 1 and in the positions from j to n do not change their positions, so the arcs linking them do not contribute to the variation of the total latency. For that reason, it is only necessary to focus on the arcs belonging to the inverted chain. If the neighborhood exploration is viewed as an iterative process starting with j = i + 2, it is easy to obtain a procedure 17

ACCEPTED MANUSCRIPT

CR IP T

with O(1) operations to evaluate a movement using the calculations of the previous movement. In (N4 ) only customers [i] and [j] change their positions, then using the formula (36) it is easy to obtain the expression (41) to evaluate the total latency variation ∆i,j = (m1 − i + 1)(c[i−1][j] − c[i−1][i] ) + (m1 − i)(c[j][i+1] − c[i][i+1] ) +(m2 − j − 1)(c[j−1][i] − c[j−1][j] ) + (m1 − j)(c[i][j+1] − c[j][j+1] )

(41)

M

AN US

where m1 and m2 are the numbers of customers in the selected routes. In (N5 ) a customer is removed from its position i in a route and it is inserted into position j in another route. A similar expression to (37) can be obtained to evaluate the total latency variation when a customer is removed from its current position in a route. To evaluate its insertion in another route, the expression (37) is used. Using these two expressions iteratively, it is possible to evaluate moves in (N5 ) with O(1) elementary operations. Despite all the strategies used to reduce computing time, exploring neighborhoods in the presence of cumulative costs is still time consuming. For this reason, we decided to use first improvement (F I) criterion to streamline the neighborhood exploration process.

ED

4. Computational experiments

CE

PT

In this section, we describe the characteristics of the data instances and the computational experiments used to evaluate the proposed models and the metaheuristic procedures. The developed mixed integer formulations were solved using solver CPLEX 12.4 under Concert technology with CPU time limit of 2 hours per instance. The metaheuristic procedures were implemented using C++. All of the experiments were performed on an Intel Core 2 Duo CPU at 2.4 GHz and 3.21 GB RAM processor, under Windows OS.

AC

4.1. Description of the used instances We use four sets of instances to conduct the experiments. The first group corresponds to instances (MT-DMP) proposed by Angel-Bello et al. (2013) and ranging from 10 to 40 customers. The second group contains an extended group of CVRP instances (A, B, E, M, P), proposed by Christofides et al. (1979) and used by Lysgaard and Wøhlk (2014). In these instances, each 18

ACCEPTED MANUSCRIPT

AN US

CR IP T

distance measure is rounded to integer values and the size of the fleet k is set as the minimum required to carry out all customer demands. The third group involves 7 instances (CMT) also proposed by Christofides et al. (1979) and ranging from 50 to 199 nodes. They are a subset of the E and M instances, but the distances are calculated using double precision. In fact, instances CMT1, CMT2, CMT3, CMT4, CMT11 and CMT12 correspond to instances E-n51-k5, E-n76-k10, E-n101-k8, M-n151-k12, M-n121-k7 and M-n101-k10, respectively. The CMT instances were used by Ngueveu et al. (2010), Ribeiro and Laporte (2012), Ke and Feng (2013) and Lysgaard and Wøhlk (2014). Finally, the last group involves instances (GKWC) ranging from 200 to 483 customers proposed by Golden et al. (1998) and used in Ribeiro and Laporte (2012), Ke and Feng (2013).

AC

CE

PT

ED

M

4.2. Comparison between proposed formulations The proposed formulations are compared regarding the number of instances solved up to optimality, and the elapsed CPU time required to reach optimal solutions. The single-flow based formulation will be referred as Formulation 1, while the time-dependent based formulation will be referred as Formulation 2. We also implemented formulations from literature (Ngueveu et al., 2010, Ribeiro and Laporte, 2012). We do not report results for these formulations because using them we could only solve instances of up to 10 clients. For 15-client instances, the optimizer always reported out of memory. Experimental results are reported in Table 1, where instances are grouped in subsets of 25 instances according to the number of nodes and vehicles. Column 1 shows the name of the group of tested instances (25 instances were evaluated per group) whereas columns 2 and 3 display the number of customers (n) and vehicles (k), respectively. Columns 4 and 8 report the number of instances solved to optimality (NSI) per group of 25 tested instances for Formulations 1 and 2, respectively. Columns 5 and 9 report the maximum CPU time required to solve to optimality any instance, while columns 6 and 10 report the minimum CPU time required in seconds. Finally, columns 7 and 11 present the average CPU time spent by the solver over each instance using Formulations 1 and 2, respectively. CPU times reported in columns 6-10 were calculated considering, only, instances solved to optimality. As can be seen from table 1, Formulation 2 is faster and solves larger size instances when compared with Formulation 1. Furthermore, Formulation 2 also presents a superior performance in terms of quality. It enables the 19

ACCEPTED MANUSCRIPT

solved instances and CPU Formulation 2 Max Min Ave 0.750 0.078 0.238 5.063 0.204 1.081 70.406 0.844 11.637 2050.170 4.782 157.405 3095.080 26.000 495.027 7428.560 34.719 1650.746

CR IP T

Table 1: Comparison of formulations regarding the number of times for MT-DMP instances Instance Formulation 1 name n k NSI Max Min Ave NSI 10 2 25 3.953 0.265 0.846 25 15 3 25 26.156 0.156 7.738 25 MT20 3 25 1569.160 14.250 414.990 25 DMP 25 3 5 4755.060 1235.500 3008.610 25 30 3 – – – – 25 40 4 – – – – 19

AC

CE

PT

ED

M

AN US

solver to reach optimal results with instances up to 40 nodes. However, with Formulation 1, no instance larger than 25 nodes was solved. The formulations performance was also evaluated over different instances used by Lysgaard and Wøhlk (2014). Results are shown in table 2. Column 1 indicates the instance name, while columns 2 and 3 display the number of customers (n) and number of vehicles (k), respectively. Columns 4 and 6 report the optimal objective function value or the objective function value for the best integer solution, while columns 5 and 7 present the CPU time (in seconds) using Formulations 1 and 2, respectively. The symbol “-” means that the solver stopped reporting out of memory without finding any feasible solution. Again, Formulation 2 outperforms Formulation 1, allowing to solve instances with up to 44 nodes (A-n45-k7). As shown in table 2, only 5 instances were solved using Formulation 1. In instance E-n23-k3, using both formulations the solver reported an out of memory situation, but it was able to identify an integer solution (in bold). This solution coincides with the optimal solution value reported by Lysgaard and Wøhlk (2014). When dealing with instances larger than 44 nodes, no optimal solution was found using any of the two formulations. Formulation 2 permits to solve instances 1.7 times larger than Formulation 1, and it confirms that a time-dependent based formulation outperforms one based on a single-flow formulation. 4.3. Performance evaluation of the metaheuristic algorithm In this subsection we describe the computational experiments carried out to evaluate the performance of the implemented metaheuristic algorithms. The IG version using the PRB constructive procedure is referred to as IGPRB, while the constructive version using a clustering procedure with complete evaluation is referred to as IG-CE. The two proposed versions were 20

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

Table 2: Evaluation of the performance of formulations on A, B, E, and P instances Instance Formulation 1 Formulation 2 name n k Obj. Value CPU Time Obj. Value CPU Time A-n32-k5 31 5 – – 2192 113.141 A-n33-k5 32 6 – – 1725 456.625 A-n33-k6 32 6 – – 1612 802.672 A-n34-k5 33 5 – – 2104 3026.860 A-n36-k5 35 6 – – 2279 103.985 A-n37-k5 36 6 – – 1747 12.297 A-n37-k6 36 6 – – 2241 14459.8 A-n45-k7 44 7 – – 2831 8803.44 B-n31-k5 30 5 – – 1830 1210.840 B-n34-k5 33 5 – – 2271 7318.73 E-n22-k4 21 4 845 717.844 845 36.297 E-n23-k3 22 3 1908 1908 E-n33-k4 32 6 – – 2852 424.516 P-n16-k8 15 8 396 0.171 396 0.297 P-n19-k2 18 2 849 3912.520 849 54.938 P-n20-k2 19 2 – – 924 26.625 P-n21-k2 20 2 – – 928 6.031 P-n22-k8 21 8 681 25.890 681 14.282 P-n23-k8 22 8 616 92.953 616 430.922 P-n40-k5 39 5 – – 1541 140.000

AC

CE

PT

ED

used to solve the test instances previously described. In order to validate the robustness of the metaheuristic, the IG-PRB and IG-CE algorithms were run 30 times per instance. After performing empirical experimentation, the values of the parameters werenadjusted. j ko The number of customers to remove (ncr ) is initialized as max 1, n8 and it is increased by 1 every 20 consecutive iterations without improvement from the best solution obtained so far. The value of ncr is reset every j k time a new best solution is found or when it reaches a value equal to n . 2 In order to tune the value of ncr we conducted the following empirical experimentation. First, we randomly selected 5 instances from each group of instances: MT-DMP, Augerat (A-P), Augerat (E&M) and GKWC. From MTDMP set, we selected the instances: MT-DMP15s1-04, MT-DMP20s1-18, MT-DMP25s1-06, MT-DMP30s1-22 and MT-DMP40s1-10. The instances taken from the other sets are marked with the symbol (•) in tables 5, 7, A.9, A.10 and A.11. Then, these 20 instances were solved 20 times for each of the 21

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 6: Average percentage gaps for each tested ncr value

AC

CE

PT

ED

M

n n n 10 values of ncr = { 11 , 10 , 9 , · · · , n2 }. For each value of ncr, we determine in each group of 5 instances the gap regarding the best solution within the group. This process is repeated 20 times and the average gap is calculated. In figure 6 are displayed the average percentages of gaps for each tested ncr value. n to n9 , the percentages of gaps are small and It can be seen that from 11 with a difference of less than 0.1 between these ncr values. For values from n8 to n2 we can observe a small increase. When ncr reaches its maximum value (ncr = n2 ), larger deviations are obtained, allowing a greater diversification in the search space. Due to this fact, we decided to start the perturbation procedure using the value of n8 and increasing it in order to escape from local minima. The value of NIL is computed as: j n k N IL = −1 (42) ncr

This way of calculating the NIL value decreases the number of blocking iterations as the value of ncr increases. The maximum value of ncr was set in such a way as to accelerate the constructive procedure to build a new initial solution. 22

ACCEPTED MANUSCRIPT

AN US

CR IP T

4.3.1. Experimental results using MT-DMP, A, B, E, M and P instances The MT-DMP instances provide travel time matrix instead of node coordinates. Therefore these instances were solved only by the IG-PRB. In order to assess the performance of the metaheuristic algorithm, optimal solutions provided by the Formulation 2 were used. Table 3 summarizes the obtained results. As in table 1, instances are grouped in subsets of 25 according to number of nodes and number of vehicles. Column 1 denotes the name of the group of instances, and columns 2 and 3 refer to the number of customers and number of vehicles, respectively. Column 4 reports the number of optimal solutions found by the metaheuristic algorithm. Columns 5, 6 and 7 report the maximum, minimum and average CPU time in seconds, respectively. Table 3: Comparison of the IG-PRB algorithm results with optimal values for MT-DMP instances

NSI 24 25 25 25 24 18

ED

MTDMP

n k 10 2 15 3 20 3 25 3 30 3 40 4

Elapsed CPU Time (seconds) Maximum Minimum Average 0.008 0.003 0.006 0.015 0.010 0.012 0.139 0.071 0.087 0.156 0.072 0.094 1.399 0.321 0.678 3.128 1.641 2.211

M

Group name

AC

CE

PT

The metaheuristic reached the optimal value in 141 out of 144 tested instances. For the 3 remaining instances, the gaps from the optimal solution are 0.007% for the 10-node instances, 0.07% for the 30-node instances, and 0.4% for the 40-node instances. In the next experiments, we compare the two versions IG-PRB and IGCE with the results of the BCP algorithm reported by Lysgaard and Wøhlk (2014). Detailed results of these experiments are reported in tables A.9 A.11 in the appendix section. Table 4 summarizes results of BCP, IG-PRB and IG-CE on A, B, E, M and P instances. In this table, column 1 indicates the name of the instance group, columns 2, 4 and 6 display the number of the best solutions (#BKS) found by each method, while columns 3, 5 and 7 show the average CPU time (Tav ) (in seconds) spent by the BCP, IG-PRB and IG-CE respectively. From table 4 it can be observed that IG-PRB is the winner regarding the number of best solutions. It found 84 optimal solutions, while BCP was capable of finding 61, and IG-CE only 59. 23

ACCEPTED MANUSCRIPT

Table 4: Summary of results of BCP, IG-PRB and IG-CE on A, B, E, M and P instances BCP #BKS Tav 26 401.56 15 1420.78 4 2575 16 759.55

IG-PRB #BKS CPU time 27 13.56 22 14.22 14 92.55 21 33.41

IG-CE #BKS Tav 19 21.68 19 40.25 7 192.93 14 30.35

CR IP T

Instance Group A B E&M P

AC

CE

PT

ED

M

AN US

4.3.2. Experimental results using CMT instances We carried out the following experiments to assess the performance of the proposed metaheuristic procedures by comparing the results obtained using state-of-the-art metaheuristic methodologies. These are memetic algorithms MA1 and MA2 (Ngueveu et al., 2010), adaptive large neighborhood search (ALNS) (Ribeiro and Laporte, 2012),Two-Phase metaheuristic (Ke and Feng, 2013) and adaptive variable neighborhood search (AVNS) (Sze et al., 2017). To do that we use CMT (Christofides et al., 1979) and GWKC (Golden et al., 1998) instances. Table 5 summarizes the best results obtained by each of the algorithms, while table 6 shows the objective, average values and the associated gaps on CMT instances for our IG procedures. In Table 5 columns 1, 2 and 3 denote the name of the instance, number of customers and number of vehicles, respectively. Columns 4 - 10 indicate the best value obtained for each of the procedures presented in the literature. In Table 6, columns 1, 2 and 3 denote the name of the instance, number of customers and number of vehicles, respectively. Column 4 displays the best known solution values (BKV) according to Sze et al. (2017). Columns 5 and 11 display the best solution values obtained by IG-PRB and IG-CE, respectively, while columns 6 and 12 show the relative gaps associated with columns 5 and 11, respectively. Columns 7 and 13 display the average objective value of the two versions of the proposed IG algorithm, columns 8 and 14 report the coefficient of variation related to the average values reported in the columns 7 and 13, while columns 9 and 15 show the relative gaps associated with these columns. Finally, columns 10 and 16 display the average CPU time in seconds. The values in columns 7, 8, 9, 10, 13, 14, 15 and 16 are averaged over 30 runs of the algorithms. Gaps are computed as: gap(ObjV al) = 100 ∗ (ObjV al − BKV )/BKV 24

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Table 5: Best known results obtained by each procedure for CMT instances Instance n k MA1 MA2 ALNS twoAVNS IGIG-CE phase PRB CMT1* 50 5 2230.35 2230.35 2230.35 2230.35 2230.35 2230.35 2230.35 CMT2* 75 10 2421.90 2429.18 2391.63 2391.63 2391.63 2391.63 2391.63 CMT3* 100 8 4073.12 4073.12 4045.42 4045.42 4045.42 4045.42 4045.42 CMT4* 150 12 4987.52 4987.52 4987.52 4987.52 4987.52 4987.52 4987.52 CMT5 • 199 17 5810.12 5810.12 5838.32 5809.59 5806.02 5826.42 5809.59 CMT11 120 7 7317.98 7347.49 7315.87 7314.55 7314.55 7320.61 7314.55 CMT12* • 100 10 3558.92 3559.23 3558.92 3558.92 3558.92 3558.92 3558.92

25

Instance Name CMT1* CMT2* CMT3* CMT4* CMT5 CMT11 CMT12* Average

n 50 75 100 150 199 120 100

k 5 10 8 12 17 7 10

AC BKV 2230.35 2391.63 4045.42 4987.52 5806.02 7314.55 3558.92

ObjVal 2230.35 2431.19 4068.27 4992.69 5840.38 7365.86 3558.92

M

Time 10.48 38.81 77.92 48.77 72.88 40.90 80.90 52.95

Gap 0.00 0.00 0.00 0.00 0.06 0.00 0.00 0.01

ObjVal 2230.35 2406.21 4045.42 4987.52 5824.93 7342.33 3558.92

Average CoefVar Gap 0.00 0.00 0.17 0.61 0.00 0.00 0.00 0.00 1.17 0.33 0.42 0.38 0.00 0.00 0.25 0.19

IG-CE

CR IP T

Best ObjVal 2230.35 2391.63 4045.42 4987.52 5809.59 7314.55 3558.92

AN US

Average CoefVar Gap 0.00 0.00 0.16 1.65 0.11 0.56 0.19 0.10 1.25 0.59 0.18 0.70 0.00 0.00 0.27 0.51

IG-PRB

ED

Gap 0.00 0.00 0.00 0.00 0.35 0.08 0.00 0.06

PT

Best ObjVal 2230.35 2391.63 4045.42 4987.52 5826.42 7320.61 3558.92

CE Table 6: Comparison with the best known solution values for CMT instances Time 7.36 21.65 64.79 48.25 62.29 22.17 61.61 41.16

ACCEPTED MANUSCRIPT

26

ACCEPTED MANUSCRIPT

CR IP T

As it can be noticed in Table 6, the two IG versions present an effective performance in terms of solution quality. For these instances, the IG-CE algorithm was able to find 6 out of 7 best known results while the IG-PRB obtained the best known solutions for 5 out of 7 instances. Regarding averages results, IG-CE obtained solutions with an average gap of less than 0.19%, while IG-PRB obtained average gap less than 0.51%.

AN US

4.3.3. Experimental results using GKWC instances Tables 7 and 8 summarize the results obtained for the GKWC instances. These tables have the same structure as the Tables 5 and 6, respectively. Solution values used in these tables for MA1, MA2, ALNS, two-phase and AVNS algorithms were also taken from Sze et al. (2017). The values in columns 7, 8, 9, 10, 13, 14, 15 and 16 in table 8 are averaged over 30 runs of the algorithms.

AC

CE

PT

ED

M

Table 7: Best known results obtained by each procedure for GKWC instances Instance n k MA1 MA2 ALNS two-phase AVNS IG-PRB GKWC1 240 9 54815.17 54826.53 54786.92 54683.36 54749.30 54881.10 GKWC2 320 10 100836.90 100800.33 100662.53 100666.06 100560.00 101368.00 GKWC3• 400 10 171277.26 171311.81 171613.59 171040.04 170924.00 172579.00 GKWC4 480 10 262584.23 262646.37 263433.03 262145.12 261993.00 266620.00 GKWC5 200 5 114163.64 114163.64 114494.66 114163.64 114163.64 114211.00 GKWC6• 280 7 140430.09 140463.67 140804.64 140430.09 140430.09 141229.00 GKWC7 360 8 183282.64 190371.53 180481.56 183509.62 179388.00 185166.00 GKWC8 440 10 194312.60 194273.58 194988.74 193855.58 193698.00 196198.00 GKWC9 255 14 4730.70 4737.32 4725.58 4723.77 4723.95 4758.02 GKWC10 323 16 6732.36 6721.16 6713.92 6714.94 6712.53 6731.82 GKWC11 399 18 9243.05 9240.37 9214.07 9217.63 9210.45 9276.19 GKWC12 483 19 12629.37 12659.15 12526.17 12549.73 12507.50 12755.90 GKWC13 252 26 3653.07 3686.62 3628.30 3635.44 3632.61 3639.14 GKWC14• 320 29 5770.02 5826.43 5216.80 5290.05 5205.75 5376.90 GKWC15 396 33 7077.48 8576.61 7010.41 7026.96 7010.13 7218.81 GKWC16 480 37 9300.74 9347.02 9250.98 9271.40 9248.37 9436.26 GKWC17 240 22 3089.99 3095.32 3065.46 3064.47 3063.02 3076.06 GKWC18• 300 27 4528.16 5083.70 4221.14 4311.06 4216.01 4315.89 GKWC19 360 33 5570.35 5600.05 5523.38 5513.72 5502.08 5547.39 GKWC20• 420 38 7413.58 8453.96 7223.08 7237.40 7209.31 7305.50

27

IG-CE 55008.40 101135.20 172774.80 266843.40 114603.40 141682.80 184514.00 196367.00 4732.72 6729.21 9230.44 12669.16 3650.09 5280.07 7091.78 9346.82 3074.99 4283.61 5546.40 7337.28

Instance Name GWKC1 GWKC2 GWKC3 GWKC4 GWKC5 GWKC6 GWKC7 GWKC8 GWKC9 GWKC10 GWKC11 GWKC12 GWKC13 GWKC14 GWKC15 GWKC16 GWKC17 GWKC18 GWKC19 GWKC20 Average

n 240 320 400 480 200 280 360 440 255 323 399 483 252 320 396 480 240 300 360 420

k 9 10 10 10 5 7 8 10 14 16 18 19 26 29 33 37 22 27 33 38

AC BKV 54683.36 100560 170924 261993 114163.64 140430.09 179388 193698 4723.77 6712.53 9210.45 12507.5 3628.3 5205.75 7010.13 9248.37 3063.02 4216.01 5502.08 7209.31

M

Time 84.29 200.12 334.86 322.19 65.36 164.62 115.14 254.99 86.91 140.19 202.66 271.40 271.40 83.43 141.65 236.24 67.90 86.90 142.28 166.92 162.00

28

ObjVal 55121.11 101409 173006.45 267263.6 115104.65 142409.25 185097.65 196275.7 4739.09 6736.76 9243.18 12716.83 3663.39 5304.59 7119.69 9343.65 3088.34 4313.51 5563.11 7343.45

Average CoefVar Gap 0.01 0.8 0.02 0.84 0.03 1.22 0.03 2.01 0.01 0.82 0.01 1.41 0.01 3.18 0.03 1.33 0.18 0.32 0.26 0.36 0.42 0.36 0.32 1.67 0.15 0.97 0.34 1.9 0.32 1.56 1.4 1.03 0.14 0.83 0.22 2.31 0.22 1.11 0.23 1.86 0.22 1.29

IG-CE

CR IP T

Best ObjVal Gap 54937.6 0.46 100962 0.4 172651 1.01 266184 1.6 114435 0.24 141313 0.63 184190 2.68 196070 1.22 4730.52 0.14 6723.56 0.16 9226.73 0.18 12662.6 1.24 3645.77 0.48 5267.54 1.19 7069.27 0.84 9326.18 0.84 3069.57 0.21 4276.49 1.43 5542.26 0.73 7314.08 1.45 0.86

AN US

Average CoefVar Gap 0.02 0.61 0.05 0.84 0.04 1.28 0.01 3.06 0.01 0.27 0.02 0.73 0.01 3.34 0.01 1.35 0.32 0.96 0.19 1.6 0.31 2.55 0.41 2.42 0.39 1.07 0.26 3.83 0.3 5.18 1.55 3.19 0.84 0.62 0.07 3.48 0.23 1.51 0.34 2.4 0.27 2.02

IG-PRB ObjVal 55016.64 101404.83 173108.27 270006.47 114473.39 141460.57 185384.38 196310.05 4768.93 6820.25 9445.40 12810.13 3667.14 5404.87 7372.96 9543.19 3082.15 4362.60 5584.98 7382.38

ED

Best ObjVal Gap 54881.1 0.36 101368 0.82 172579 1.09 266620 2.78 114211 0.04 141229 0.57 185166 3.22 196198 1.29 4758.02 0.73 6731.82 0.41 9276.19 1.27 12755.9 1.99 3639.14 0.88 5376.9 3.29 7218.81 3.87 9436.26 2.03 3076.06 0.63 4315.89 2.37 5547.39 1.41 7305.5 2.21 1.56

PT

CE Table 8: Comparison with the best known solution values for GKWC instances Time 220.13 604.78 1151.41 1564.36 80.80 304.43 481.88 1260.61 293.59 632.18 1242.19 1895.16 261.83 465.30 931.60 2036.26 161.43 282.77 557.14 874.85 765.14

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

CR IP T

The results obtained on GKWC instances confirm the ones obtained with CMT instances. The IG-CE algorithm obtained better objective values than the IG-PBR algorithm. From Table 8, we can see that IG-CE obtained solutions with an average gap of 0.86% regarding the BKV, while the average gap obtained by IG-PBR algorithm was 1.56%. In summary, these results show that the IG-PRB algorithm performed better when dealing with small-size instances (MT-DMP, A, B, E, M and P), while IG-CE algorithm produces better results for medium and large-size instances (CMT and GKWC).

AN US

5. Conclusions

AC

CE

PT

ED

M

In this paper, the first two tractable models for the CCVRP are presented. The first formulation follows a single-flow approach; the second formulation is based in the time-dependent TSP. The time-dependent formulation proved to be more efficient than the single-flow formulation, allowing CPLEX to solve larger instances and spending a significant less CPU time for solving instances of the same size. Particularly, the time-dependent formulation was able to solve instances up to 45 nodes, while the single-flow based formulation solved instances up to 25 nodes. Two versions of an iterated greedy algorithm are presented. One version with a constructive procedure based on a PRB strategy and another one with a procedure based on clustering. The Iterated greedy algorithm based on a PRB procedure showed a good performance for solving small instances, finding the optimal solution for 141 out of 144 instances. The Iterated greedy algorithm with a complete evaluation constructive procedure showed a competitive performance for solving larger instances, obtaining results with less than 1% gap when compared with the best solution values from the literature. Moreover, for A, B, E, M and P instances where the optimal solution is not known, the metaheuristics were able to find feasible solutions with objective values that improve the best known results. In particular, the solution for the instance E-n30-k3 was proved to be optimal, and a new feasible solution for the instance Pn55k8 is presented, which improves the best result known for this instance, reported as optimal in the literature. Future works involve extensions of the CCVRP considering multiple trips per vehicle, multiple depots and location-routing problems. Other interesting research avenues could be to include multiple objectives or to consider priority indexes in the customers. 29

ACCEPTED MANUSCRIPT

Acknowledgements

CR IP T

This work was partially supported by Universidad Panamericana through the grant ”Fomento a la Investigaci´on UP 2017”, under project code UPCI-2017-ING-GDL-07 and by the Mexican National Council of Science and Technology (Proyecto CONACyT 280081). These supports are gratefully acknowledged. References

AN US

Abeledo, H., Fukasawa, R., Pessoa, A., and Uchoa, E. (2013). The time dependent traveling salesman problem: polyhedra and algorithm. Mathematical Programming Computation, 5(1):27–55.

Angel-Bello, F., Alvarez, A., and Garc´ıa, I. (2013). Two improved formulations for the minimum latency problem. Applied Mathematical Modelling, 37(4):2257–2266.

M

´ Angel-Bello, F., Cardona-Vald´es, Y., and Alvarez, A. (2017). Mixed integer formulations for the multiple minimum latency problem. Operational Research, pages 1–30.

ED

Bjeli´c, N., Vidovi´c, M., and Popovi´c, D. (2013). Variable neighborhood search algorithm for heterogeneous traveling repairmen problem with time windows. Expert Systems with Applications, 40(15):5997–6006.

CE

PT

Blum, A., Chalasani, P., Coppersmith, D., Pulleyblank, B., Raghavan, P., and Sudan, M. (1994). The minimum latency problem. In Proceedings of the twenty-sixth annual ACM symposium on Theory of computing, pages 163–171. ACM.

AC

Christofides, N., Mingozzi, A., and Toth, P. (1979). The vehicle routing problem. In N. Christofides, A. Mingozzi, P. Toth and C. Sandi, eds., Combinatorial optimization. John Wiley and Sons, London. Cinar, D., Gakis, K., and Pardalos, P. M. (2016). A 2-phase constructive algorithm for cumulative vehicle routing problems with limited duration. Expert Systems with Applications, 56:48–58. Dantzig, G. B. and Ramser, J. H. (1959). The truck dispatching problem. Management science, 6(1):80–91. 30

ACCEPTED MANUSCRIPT

Ezzine, I. O., Semet, F., and Chabchoub, H. (2010). New formulations for the traveling repairman problem. In Proceedings of the 8th International Conference of Modeling and Simulation, pages 1889–1894.

CR IP T

Gaur, D. R., Mudgal, A., and Singh, R. R. (2016). Approximation algorithms for cumulative vrp with stochastic demands. In Conference on Algorithms and Discrete Applied Mathematics, pages 176–189. Springer.

Gavish, B. and Graves, S. C. (1978). The travelling salesman problem and related problems. Working paper, Massachusetts Institute of Technology, Operations Research Center.

AN US

Godinho, M. T., Gouveia, L., Magnanti, T. L., Pesneau, P., and Pires, J. (2009). On the unit demand vehicle routing problem: Flow based inequalities implied by a time-dependent formulation. In Proceedings of International Network Optimization Conference (INOC).

M

Golden, B. L., Wasil, E. A., Kelly, J. P., and Chao, I. M. (1998). Metaheuristics in vehicle routing. In Fleet management and logistics, pages 33–56. Kluwer, Boston.

ED

Kallehauge, B., Larsen, J., and Madsen, O. B. (2006). Lagrangian duality applied to the vehicle routing problem with time windows. Computers & Operations Research, 33(5):1464–1487.

PT

Kara, I., Kara, B. Y., and Kadri Yeti, M. (2008). Cumulative vehicle routing problems. In T. Caric and H. Gold eds., Vehicle Routing Problem, pages 85–98. InTech.

CE

Ke, L. and Feng, Z. (2013). A two-phase metaheuristic for the cumulative capacitated vehicle routing problem. Computers & Operations Research, 40(2):633–638.

AC

Lysgaard, J. and Wøhlk, S. (2014). A branch-and-cut-and-price algorithm for the cumulative capacitated vehicle routing problem. European Journal of Operational Research, 236(3):800–810. Mart´ınez-Salazar, I., Angel-Bello, F., and Alvarez, A. (2014). A customercentric routing problem with multiple trips of a single vehicle. Journal of the Operational Research Society, 66(8):1312–1323. 31

ACCEPTED MANUSCRIPT

M´endez-D´ıaz, I., Zabala, P., and Lucena, A. (2008). A new formulation for the traveling deliveryman problem. Discrete Applied Mathematics, 156(17):3223–3237.

CR IP T

Mladenovi´c, N. and Hansen, P. (1997). Variable neighborhood search. Computers & Operations Research, 24(11):1097–1100.

Mladenovi´c, N., Uroˇsevi´c, D., and Hanafi, S. (2013). Variable neighborhood search for the travelling deliveryman problem. 4OR, 11(1):57–73.

AN US

Mulvey, J. M. and Beck, M. P. (1984). Solving capacitated clustering problems. European Journal of Operational Research, 18(3):339–348. Ngueveu, S. U., Prins, C., and Calvo, R. W. (2010). An effective memetic algorithm for the cumulative capacitated vehicle routing problem. Computers & Operations Research, 37(11):1877–1885.

M

Nucamendi-Guill´en, S., Mart´ınez-Salazar, I., Angel-Bello, F., and MorenoVega, J. M. (2016). A mixed integer formulation and an efficient metaheuristic procedure for the k-travelling repairmen problem. Journal of the Operational Research Society.

ED

Potvin, J.-Y. and Rousseau, J.-M. (1993). A parallel route building algorithm for the vehicle routing and scheduling problem with time windows. European Journal of Operational Research, 66(3):331–340.

PT

Ribeiro, G. M. and Laporte, G. (2012). An adaptive large neighborhood search heuristic for the cumulative capacitated vehicle routing problem. Computers & Operations Research, 39(3):728–735.

CE

Rivera, J. C., Afsar, H. M., and Prins, C. (2015). A multistart iterated local search for the multitrip cumulative capacitated vehicle routing problem. Computational Optimization and Applications, 61(1):159–187.

AC

¨ Salehipour, A., S¨orensen, K., Goos, P., and BRAYSY, O. (2008). An efficient grasp+ vnd metaheuristic for the traveling repairman problem. Technical report.

Salehipour, A., S¨orensen, K., Goos, P., and Br¨aysy, O. (2011). Efficient grasp+ vnd and grasp+ vns metaheuristics for the traveling repairman problem. 4OR, 9(2):189–209. 32

ACCEPTED MANUSCRIPT

Silva, M. M., Subramanian, A., Vidal, T., and Ochi, L. S. (2012). A simple and effective metaheuristic for the minimum latency problem. European Journal of Operational Research, 221(3):513–520.

CR IP T

Subramanian, A., Drummond, L. M. d. A., Bentes, C., Ochi, L. S., and Farias, R. (2010). A parallel heuristic for the vehicle routing problem with simultaneous pickup and delivery. Computers & Operations Research, 37(11):1899–1911.

AN US

Sze, J. F., Salhi, S., and Wassan, N. (2017). The cumulative capacitated vehicle routing problem with min-sum and min-max objectives: An effective hybridisation of adaptive variable neighbourhood search and large neighbourhood search. Transportation Research Part B: Methodological, 101:162 – 184. Appendix A. Detailed results of PRB and CE procedures over A, B, E, M and P instances

AC

CE

PT

ED

M

In tables A.9, A.10, and A.11, column 1 refers to the instance name, columns 2 and 3 refer to the instance size in terms of number of customers and number of vehicles, respectively. Column 4 reports the best known solution values (BKV) reported by Lysgaard and Wøhlk (2014). Columns 5 and 8 indicate the best objective value (bestVal ), while columns 6 and 9 report the average objective value (averVal ) for each IG procedure. Finally, columns 7 and 10 report the coefficient of variation related to the average values reported in the columns 6 and 9. Instance name marked with ’*’ indicates that the solution optimality is not proven. To determine the average values the algorithms were run 30 times. From tables A.9, A.10 and A.11, it can be observed that the two versions performed well over every tested A, B, E, M and P instance. In particular, the IG-PRB algorithm outperformed IG-CE algorithm obtaining equal or better solutions for every tested instance. However, for instances with proved optimal solutions, the IG-PRB algorithm found 61 out of 62 instances, while IG-CE algorithm found only 58. In addition, in instances that the IG-PRB algorithm did not reach optimal solutions, the maximum gap was around 2.90% from the best found solution, and increased to 2.99% for the average solution values.

33

CR IP T

ACCEPTED MANUSCRIPT

Table A.9: Comparison of the metaheuristic results with the best known values on A instances IG-PRB averVal CoefVar(%) 2192 0.00 1725 0.00 1612 0.00 2104.95 0.20 2279 0.00 1970 0.00 2241 0.00 2084 0.00 2312 0.00 2216 0.00 2563 0.00 2883.05 1.12 2831 0.00 2373 0.00 3101 0.00 3118.4 0.13 3364.45 0.39 2606.32 0.29 3447 0.02 2905.15 1.05 3926.9 0.07 4633.8 0.25 3260.45 0.27 4138 0.13 3488.25 0.16 3538.55 0.39 5936.8 0.16

AC

CE

bestVal 2192 1725 1612 2104 2279 1970 2241 2084 2312 2216 2563 2848 2831 2373 3101 3115 3357 2588 3446 2868 3925 4630 3256 4135 3487 3528 5929

bestVal 2192 1725 1612 2104 2279 1970 2241 2084 2312 2216 2563 2848 2831 2373 3101 3122 3382 2588 3446 2869 3925 4681 3256 4137 3549 3542 5944

AN US

BKV 2192 1725 1612 2104 2279 1970 2241 2084 2312 2216 2563 2848 2831 2373 3101 3115 3357 2588 3446 3652 3925 4630 3256 4135 3487 3528 7174

M

k 5 5 6 5 5 5 6 5 5 6 6 6 7 7 7 7 7 9 9 9 8 9 10 9 9 9 10

ED

n 31 32 32 33 35 36 36 37 38 38 43 44 44 45 47 52 53 54 59 60 61 62 62 63 64 68 79

PT

Instance Name A-n32-k5 • A-n33-k5 A-n33-k6 A-n34-k5 A-n36-k5 A-n37-k5 A-n37-k6 A-n38-k5 A-n39-k5 A-n39-k6 • A-n44-k6 A-n45-k6 A-n45-k7 A-n46-k7• A-n48-k7 A-n53-k7 A-n54-k7 A-n55-k9• A-n60-k9 A-n61-k9* A-n62-k8 A-n63-k9 A-n63-k10 A-n64-k9 A-n65-k9 A-n69-k9 A-n80-k10*•

34

IG-CE averVal CoefVar(%) 2192 0.00 1725 0.00 1612 0.00 2104 0.00 2279 0.00 1970 0.00 2241 0.00 2091.05 0.45 2317.70 0.06 2216 0.00 2583.10 0.10 3002.60 0.25 2831 0.00 2373 0.00 3101 0.00 3128.95 0.02 3416.20 0.09 2606.44 0.16 3448.05 0.07 2955.45 0.10 3930.40 0.05 4720.55 0.53 3282.95 0.31 4148.00 0.15 3657.80 1.16 3571.25 0.38 5966.45 0.21

ACCEPTED MANUSCRIPT

CE AC

bestVal 1830 2271 2846 2164 1960 2329 2123 2295 2386 2057 2261 2953 3133 2573 2358 3885 4500 4379 2608 4132 2868 4058 845 1908 1984 2852 2213 2920 2686 2366 2093 3954 2955 3565 7230 4917

IG-PRB averVal CoefVar(%) 1830 0.00 2271 0.00 2846 0.00 2164 0.00 1960 0.00 2329 0.00 2123 0.00 2296.6 0.14 2386 0.35 2067.7 0.41 2261 0.00 2953 0.00 3140.9 0.21 2573 0.00 2358 0.00 3910.4 0.23 4500 0.00 4380 0.05 2618.4 0.24 4142.9 0.26 2868 0.00 4058 0.00 845 0.00 1908 0.00 1984 0.00 2852 0.00 2243.9 0.77 2928.25 0.86 2691 0.50 2403.55 0.35 2105.55 0.29 3984.95 0.22 2963.25 0.52 3570.85 0.17 7273.85 0.48 4991.5 1.02

bestVal 1830 2271 2846 2164 1960 2329 2123 2295 2386 2110 2261 2953 3133 2573 2358 4340 4500 4379 2621 4175 2868 4058 845 1908 1984 2852 2213 2920 2703 2557 2093 3968 2976 3574 7359 4980

AN US

BKV 1830 2271 2846 2103 1968 2329 2123 2295 2386 – 2293 2953 3133 2573 2358 – 4500 4379 3222 4872 2871 4058 845 1908 1987 2852 2213 3503 3674 3364 – 5166 4022 3566 8539 6314

M

k 5 5 5 6 5 6 6 7 5 6 7 8 7 7 7 7 9 10 9 9 10 9 4 3 3 4 5 7 8 10 14 8 14 10 7 12

ED

n 30 33 34 37 38 40 42 43 44 44 49 49 50 51 55 56 56 62 63 65 66 67 21 22 29 32 50 75 75 75 75 100 100 100 120 150

PT

Instance Name B-n31-k5 B-n34-k5• B-n35-k5 B-n38-k6 B-n39-k5* B-n41-k6 B-n43-k6 B-n44-k7 B-n45-k5 B-n45-k6*• B-n50-k7* B-n50-k8 B-n51-k7 B-n52-k7 B-n56-k7• B-n57-k7* B-n57-k9 B-n63-k10 B-n64-k9* • B-n66-k9* B-n67-k10* B-n68-k9 E-n22-k4 E-n23-k3 E-n30-k3* E-n33-k4 E-n51-k5 E-n76-k7* E-n76-k8* E-n76-k10* E-n76-k14* E-n101-k8* E-n101-k14* M-n101-k10* M-n121-k7* M-n151-k12*

CR IP T

Table A.10: Comparison of the metaheuristic results with the best known values on B, E and M instances

35

IG-CE averVal CoefVar(%) 1830 0.00 2271 0.00 2846 0.00 2164 0.00 1960 0.00 2330.05 0.03 2123 0.00 2295.6 0.03 2413.7 0.13 2181.3 1.60 2261.1 0.03 2954.4 0.06 3631.15 0.23 2574.7 0.05 2358 0.00 4521.5 1.80 4500.6 0.02 4394.8 0.21 2700.8 1.09 4191.65 0.35 2873.4 0.12 4071.3 0.16 845 0.00 1908 0.00 1984 0.00 2852 0.00 2294.15 1.93 2403.55 1.11 2101.52 3.10 2960.47 0.32 2725.42 0.53 2997.55 0.29 4004.05 0.58 3583.5 0.21 7345.49 0.96 4987.52 0.79

CR IP T

ACCEPTED MANUSCRIPT

Table A.11: Comparison of the metaheuristic results with the best known values on P instances IG-PRB averVal CoefVar(%) 396 0.00 849 0.00 924 0.00 928 0.00 681 0.00 616 0.00 1542.6 0.21 1901.95 0.89 1554 0.00 1552.95 0.88 1347.15 0.03 1496.25 0.74 1764 0.00 1620 0.00 1463 0.00 1704 0.00 1510.1 0.07 1956.3 0.26 2147.8 1.22 4645.85 0.75 3843.15 1.21 7010.55 0.99

AC

CE

bestVal 396 849 924 928 681 616 1541 1894 1554 1533 1347 1489 1764 1620 1463 1704 1509 1948 2121 4610 3795 6943

bestVal 396 849 924 928 681 616 1541 1894 1554 1682 1347 1489 1764 1620 1463 1705 1517 1949 2264 4677 3916 6943

AN US

BKV 396 849 924 928 681 616 1541 1894 1554 – 1347 1487 1764 1768 1463 1704 1509 1948 2121 6191 5463 8707

M

k 8 2 2 2 8 8 5 5 7 8 10 10 7 8 10 10 15 10 10 4 5 4

ED

n 15 18 19 20 21 22 39 44 49 49 49 50 54 54 54 59 59 64 69 75 75 100

PT

Instance Name P-n16-k8 P-n19-k2• P-n20-k2 P-n21-k2 P-n22-k8 P-n23-k8 P-n40-k5 • P-n45-k5 P-n50-k7 P-n50-k8* P-n50-k10 P-n51-k10 P-n55-k7 P-n55-k8*• P-n55-k10 P-n60-k10 P-n60-k15 P-n65-k10 P-n70-k10• P-n76-k4* P-n76-k5* P-n101-k4* •

36

IG-CE averVal CoefVar(%) 396 0.00 849 0.00 925.55 0.30 928 0.00 681 0.00 616 0.00 1541.6 0.13 1894 0.00 1554.6 0.06 1770.75 3.00 1355.75 0.38 1500.3 1.26 1764 0.00 1620 0.00 1463.3 0.00 1719.2 0.34 1526.45 0.40 1967.35 0.43 2308.4 1.10 4756.85 0.95 3972.75 0.73 7020.1 1.14