The open vehicle routing problem with decoupling points

Accepted Manuscript

The Open Vehicle Routing Problem with Decoupling Points Reza Atefi , Majid Salari , Leandro C. Coelho , Jacques Renaud PII: DOI: Reference:

S0377-2217(17)30660-4 10.1016/j.ejor.2017.07.033 EOR 14583

To appear in:

European Journal of Operational Research

Received date: Revised date: Accepted date:

29 March 2017 5 July 2017 7 July 2017

Please cite this article as: Reza Atefi , Majid Salari , Leandro C. Coelho , Jacques Renaud , The Open Vehicle Routing Problem with Decoupling Points, European Journal of Operational Research (2017), doi: 10.1016/j.ejor.2017.07.033

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

CR IP T

We describe a practical routing problem arising in the distribution of cookies The problem is a generalization of the open vehicle routing problem At some locations the load is dropped off and delivered by another carrier Real data were provided by our industrial partner and its carriers The proposed algorithm leads to the best known solutions for the classical problem

AC

CE

PT

ED

M

AN US

    

1

ACCEPTED MANUSCRIPT

The Open Vehicle Routing Problem with Decoupling Points Reza Atefi1, Majid Salari1, Leandro C. Coelho2,3,4, Jacques Renaud2,3,* 1

AN US

July, 2017

CR IP T

Department of Industrial Engineering, Ferdowsi University of Mashhad, P.O. Box 91779-48951, Mashhad, Iran 2 Département opérations et systèmes de décision, Faculté des sciences de l’administration, Université Laval, Québec, G1V 0A6, Canada 3 Centre interuniversitaire de recherche sur les réseaux d’entreprise, la logistique et le transport (CIRRELT), Québec, Canada 4 Canada Research Chair in Integrated Logistics, Québec, Canada * Corresponding author: [email protected]

Abstract

AC

CE

PT

ED

M

In this paper we introduce the open vehicle routing problem with decoupling points (OVRP-DP). This practical problem is faced by companies dealing with carriers to ship their goods over large territories. In this case it may be profitable to use more than one carrier to perform a specific expedition: the first one leaves the depot and performs part of the deliveries, drops off all remaining load, and the second carrier continues from that point onwards. This drop off location is called the decoupling point of the route. This problem generalises the classical OVRP in which each route must be performed by only one carrier. We model this problem using a realistic multi-drop less-than-truckload cost function composed of a non-linear transportation cost, a detour cost and a drop cost. We have developed a tailored Iterated Local Search (ILS) algorithm which handles the special features of the problem. The efficiency of the ILS was demonstrated by obtaining all best known solutions on a set of classical OVRP instances and improving it for one instance. Then, using real orders and transportation costs obtained from industrial partners, we clearly show the benefit of using decoupling points to optimize transportation costs. The performance of the ILS is analysed and shown to be very robust and superior to what can be obtained with a commercial solver.

Keywords : Open vehicle routing, common carriers, iterated local search, decoupling points

2

ACCEPTED MANUSCRIPT

1

Introduction

Many companies face practical transportation situations that could be modeled and solved as variants of the well-known vehicle routing problem (VRP). The objective of the classical VRP is to find the best set of routes starting and ending at a depot for a fleet of identical capacitated vehicles such that each customer is visited exactly once (Laporte

CR IP T

2009, Toth and Vigo 2014). Customers’ orders can be defined by a weight and a volume, and the vehicle capacity can be measured in terms of these two quantities. Extensions of the problem also consider a maximum route duration, which is often used when planning home deliveries. Coelho et al. (2016) reviewed practical applications of the VRP for oil, gas, fuel, retail applications, waste management, mail and small package delivery, and

AN US

food distribution.

In other contexts, vehicles need not return to their depot at the end of the route, as it is the case with external carriers, for example. When one hires such a carrier, it is paid for the delivery based on distance and other business factors, but one does not (explicitly) pay for the vehicle return trip, which may be empty or not, depending on the carrier ability to

M

consolidate its operations. This situation was first described by Schrage (1981). When vehicles are not required to return to the depot, they simply finish their routes at the last

ED

visited customer; in this case, we refer to the open vehicle routing problem (OVRP). A solution to the OVRP is then a set of Hamiltonian paths connected to the depot, as

PT

opposed to Hamiltonian cycles in the VRP. The objective of the OVRP is generally to minimize the number of vehicles and/or a

CE

measure of cost incurred by the trips of the vehicles (Sariklis and Powell 2000, Brandão 2004, Li et al. 2007). For both the VRP and the OVRP a cost cij is paid when vehicles

AC

travel from node i to j. Generally, costs cij are considered as a function of the distance between i and j, the most used ones being the Euclidean or the road distance. Sariklis and Powell (2000) were among the first to propose a heuristic for the OVRP. The method was a sequential two-phase algorithm starting with a clustering phase and followed by a routing procedure based on the use of minimum spanning trees. Brandão (2004) developed a tabu search algorithm based on a nearest neighbor method followed by an improvement phase based on the unstringing and stringing procedure of Gendreau 3

ACCEPTED MANUSCRIPT

et al. (1992) to Hamiltonian paths. Tarantilis et al. (2004a) developed a decision support system based on a metaheuristic algorithm and a geographical information system to solve the OVRP. They later proposed a threshold accepting approach (Tarantilis et al. 2004b) and a single parameter metaheuristic (Tarantilis et al. 2005). Fu et al. (2005) developed a tabu search algorithm in which the initial solution is based on an

CR IP T

optimization-based farthest first heuristic. They used vertex reassignment, vertex swap, 2-opt and a tails swap to improve the solution. Li et al. (2007) developed a set of largescale test problems and applied a record-to-record algorithm. A branch-and-cut algorithm for the minimum cost capacitated OVRP was developed by Letchford et al. (2007). Repoussis et al. (2007) developed a route-construction insertion-based sequential

AN US

approach for solving the OVRP with time windows. Derigs and Reuter (2009) solved the OVRP with an attribute-based hill climber heuristic and Fleszar et al. (2009) used a variable neighborhood search algorithm. Repoussis et al. (2010) developed a hybrid evolution strategy for the OVRP. They reported results against the algorithm of Tarantilis et al. (2004a, 2004b) and of Tarantilis et al. (2005). Salari et al. (2010) developed a

M

destruct-and-repair method in which customers are removed in a random way and reinserted by solving an integer linear model. Zachariadis and Kiranoudis (2010)

ED

developed a local search metaheuristic which examines wide solution neighborhoods. An algorithm based on the iterative variation and selection concept of Reinholz (2008) that integrates a stochastic multiple neighborhood search in a (1+1)-evolutionary strategy was

PT

proposed by Reinholz and Schmeider (2013). Finally, a multiple neighborhood search algorithm hybridized with a tabu search strategy for the multi-depot OVRP was proposed

CE

by Soto et al. (2017), who also presented results for classical OVRP instances. Of relevance to the context of this paper are two variants of the OVRP. Aksen et al.

AC

(2007) studied a version where the vehicles end their routes at special driver nodes. They also investigated the problem under a maximum route duration and with time deadlines for visiting the customers. They developed a tabu search based on the parallel Clarke and Wright (1964) algorithm followed by node exchanges. Another variant of the OVRP in which the maximum time spent on the vehicle by one person must be minimized was studied by López-Sánchez (2014). This version was motivated by a real world application in Spain where employees must be picked up in several locations and taken to their 4

ACCEPTED MANUSCRIPT

working place. We finally note that in the relay network design in freight transportation (Tarab et al., 2002, Uster and Kewcharoenwong, 2011) one must identify relay nodes where drivers can rest while performing routes that last several days. Obviously, this problem is relevant when incorporating work hours regulation into the optimization, which is not the case of the OVRP or of our context.

CR IP T

As it can be seen, the OVRP has been studied as a direct extension of the classical VRP, at the exception of an objective function which is generally composed of minimizing the number of vehicles and the traveling distance. We note, however, that some papers focus on the single objective version of the problem, minimizing only transportation costs, such as Soto et al. (2017). In our paper, we present a new practical extension of the OVRP

AN US

which is based on a collaboration with one of the most important Canadian producers of cookies and snacks. In their daily operations, they assign less-than-truckload deliveries to their carriers. This leads to two practical and important generalizations of the OVRP. The first one is to consider the real and more complex cost structure of carriers. The second one is to analyze the profitability of using a carrier up to a given point (called a

M

decoupling point), and continuing the deliveries with another carrier which may offer better tariffs for the last portion of the route. To the best of our knowledge it is the first

ED

time that this problem, that we call the open vehicle routing problem with decoupling points (OVRP-DP), is studied.

PT

In Section 2 we describe the practical motivation and context of the OVRP-DP. Section 3 formally describes the OVRP-DP, its cost components and presents an integer linear

CE

model for solving it. An Iterated Local Search (ILS) solving method is presented in Section 4. Practical context and data are presented in Section 5 along with extensive

AC

computational experiments. Results are also reported on standard OVRP instances and show that the developed method consistently obtains the best known results, leading also to a new best known solution. Our conclusions follow in Section 5.

2

Practical motivation and context

When companies need to transport their products over very long distances, they often require service of external carriers. Broadly speaking, external carriers offer three kinds

5

ACCEPTED MANUSCRIPT

of services: parcels, truckload (TL) and less-than-truckload (LTL) (Lapierre et al. 2004). Parcels are generally associated with courier companies delivering letters and small boxes. At the opposite, a TL is used when an expeditor requires a large portion of a vehicle capacity such that the carrier will not need to consolidate and add more cargo. For TL, shipping rates are rather simple and mainly depend on the origin-destination. In LTL

CR IP T

shipments, the carrier consolidates loads from different expeditors and performs the transportation to various destinations. Shipping rates are much more complicated and may depend on the products’ characteristics. In practice, rates may be based on dollars per hundred pounds (CTW tariff), on the number of pallets or on the linear feet requirement. Two other components also complicate the shipping rate structure. First, the

AN US

unit rates decrease when the quantity shipped increases. This means that a heavier shipment will have smaller rates per unit. Second, carriers may modulate their rates depending on the final destination. Carriers will more likely offer good rates for destinations where they can easily find return freight. Alternatively, they may charge higher rates for destinations where they cannot find cargo to bring back to the area of

M

their depot, or which are very far away from their original hub.

Our problem is motivated by one of the most important Canadian manufacturer of

ED

cookies, cereal bars and other light snacks. Our partner recently signed contracts with two major American retailers. As a result, every week the company needs to supply its

PT

clients’ distribution centers (DCs) disseminated throughout all the United States. Figure 1 shows the location of our partner distribution center as a diamond and those of both new

CE

customers correspond to circles and squares. In order to optimize its distribution costs, our partner negotiated multi-drop routes with

AC

its carriers. This means that each truck departing from our partner’s DC will be almost full and will visit some destinations, generally between two and five. This transportation scheme is a variant of the TL and the LTL, and its cost is adjusted (penalized) for each additional destination (drop cost) and the additional mileage incurred in this detour.

6

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 1. Locations of our partner distribution center (diamond) and all distribution centers of both new customers (circles and squares)

When analyzing the transportation cost structure, it became apparent that sending a

M

carrier too far away from its operational base significantly increased the costs. As an example, a carrier located in Quebec may offer a good rates for covering the Quebec-

ED

Alabama segment. However, covering the Alabama-San Diego route may be significantly more costly as it brings the Canadian carrier very far away from its operating base in

PT

Quebec. In this case, some carriers operating around Alabama may offer much better tariffs to San Diego. By knowing the tariffs of different carriers located in different cities,

CE

it is possible to create routes combining them, such as shipping a large load from Quebec to Alabama, and in Alabama hiring a local carrier to continue the deliveries up to San

AC

Diego. In this case, Alabama is called the decoupling point of the route. Based on this real application, we introduce, model and solve the problem of determining the best decoupling points for OVRP routes, such as to take advantage of the scenario just described.

7

ACCEPTED MANUSCRIPT

3

Problem definition and mathematical formulation

The OVRP-DP is defined over a graph G = (V, A) where V is the node set and A is the arc set. In V, node 0 corresponds to our partner DC and N = {1, …, n} is the set of the n cities where the client’s DC are located, V = N  {0}. We consider qi as the demand of DC i (this can be a weight, a number of pallets or of linear-feet). The arc set A represents the

CR IP T

links between the nodes and dij is the distance between i and j. Without loss of generality, we consider that we have access to an unlimited number of a single type of vehicle of capacity Q, as long-distance transportation is almost all performed by 53 linear-feet trailer trucks (thus, Q = 53). In the following we define a route as a vehicle starting from the DC 0, visiting some clients and finishing its trip at the last location i. The last location

AN US

of a route is called the decoupling point. In Figure 2, route k=2 visits nodes 1, 2 and 3 which is the decoupling point. A subroute starts from a last client location which has been visited by a route, i.e., the decoupling point i, and visits another set of clients, finishing its subroute at its last client location. In Figure 2, subroute k=1 strarts from decoupling point 9 and visits node 10. Subroute k=2 starts from decoupling point 3 and

M

visits nodes 6 and 7. A subroute always follows a route starting at the decoupling point. A route may be followed by a subroute or not. In Figure 2, routes k=3 and k=4 are not

ED

followed by a subroute. A vehicle performing a route must carry the total load of its route and of any subsequent subroute. Hence, a route ends at a decoupling point, only one

locations.

PT

subroute follows a decoupling point, and decoupling points correspond to customer

CE

The objective of the OVRP-DP is to determine the set of routes and subroutes visiting all customers and minimizing the transportation cost function. The cost function is described

AC

next.

8

ACCEPTED MANUSCRIPT

0

k=1

k=4

k=2 9

1

10

Depot

8

k=3

Clients visited before the decoupling point

4 2

Clients visited after the decoupling point

5 3

CR IP T

Decoupling point

6

Arcs before the decoupling point 7

Arcs after the decoupling point

AN US

Figure 2. Different types of routes and subroutes

In order to state the cost function let us also define

 V as the subset of customers

included in a route (subroute) starting at i and ending at j. Similarly, let ∑

. Node j is included in

M

load of the route (subroute),

i. As each customer must be delivered by only one vehicle, are in at least one route, ⋃

=  and all clients

ED

PT

(

)

( )

is defined as:

(| |

)

(1)

CE

( )

, but not node

. When operating under multiple-drop tariffs, the total

cost of a route/subroute delivering customers ( )



and

be the total

AC

These costs are detailed next. 1. Transportation cost: This term consists of the cost to hire a carrier to leave its depot with a given load and to drive up to a last node, as is the current practice and normal in OVRP settings. We define the transportation cost

as a function

of the initial node i (starting location), of the load and final destination determined from

, hence

(

).

9

ACCEPTED MANUSCRIPT

2. Detour cost: The second term is the detour cost, defined

. It represents the cost

of the extra miles driven to reach other nodes different from the origin i and the final destination j, i.e., all intermediary nodes in {

}, the detour cost is computed as ) where

( ). Given that

, hence (

is unit cost per extra mile, defined by the carrier.

CR IP T

3. Drop cost: The third and last term of the total cost function is associated with each and is paid (| |

additional stop (the first drop is free). This cost is called times, hence 3.1

(| |

).

Numerical example

)

AN US

Figure 3 shows an example of a multi-drop route starting from Quebec City (node A) to Baltimore (node B), moving to Charlotte (node C) and ending at Dallas (node D). Table 1 presents the road distances and the demand of each destination in linear-feet. In this example the cost per mile of a detour is 1.50 $ and the cost per additional drop is

.

Table 2 gives real shipment cost (in $) for these destinations and different number of

M

linear-feet dispatched.

CE

PT

ED

Table 1. Road distances (in miles) and demand (in linear-feet) Cities (vi) Demand qi Baltimore Charlotte Dallas Quebec City (vA) 703 1 141 1 902 Baltimore (vB) 9 455 1 360 Charlotte (vC) 14 1 032 Dallas (vD) 16

AC

From Quebec to Baltimore (9) Charlotte (14) Dallas (16) From Baltimore to Charlotte (14) Dallas (16) From Charlotte to Dallas (16)

Table 2. Real shipment costs (in $) Number of linear-feet 9 14 16 30 511

1 335 1 454 2 762

659 1 311 394 818 619 10

39

684 1302

AN US

CR IP T

ACCEPTED MANUSCRIPT

ED

M

Figure 3. A Quebec – Baltimore – Charlotte – Dallas route

Table 3 details costs function (1) for four different shipment options:

2 481 $.

PT

 Option 1: Three individual trips out of Quebec for a total cost of 511 + 659 + 1 311 =  Option 2: One trip from Quebec to Baltimore with 39 linear-feet, and a second trip

CE

starting in Baltimore to Charlotte and Dallas. The total cost is then composed of 1 335$ (39 linear-feet from Quebec to Baltimore), plus 1 302$ (30 linear-feet from

AC

Baltimore to Charlotte and Dallas), plus a detour cost of 190.5$ (127  1.50) plus 100 $ as a drop cost at Charlotte. The total cost is then 2 925.50$.

 Option 3: Similar to Option 2 but using Charlotte as a decoupling point for a total cost of 2 201.50$.  Option 4: A single long trip starting in Quebec and ending in Dallas. In this case, transportation, detour and drop cost are respectively 2 762, 28.5 and 200 for a total cost of 2 990.50$. 11

ACCEPTED MANUSCRIPT

Table 3. Cost calculation for four combinations of shipments

{vB } {vC} {vD}

( ) ( ) ( )

Option 2 {vA, vB } {vB, vC, vD}

{vB } {vC, vD}

( ) ( )

)

( ) 0 0 0

0 (

{vC, vD}

{vD}

Option 4 {vA, vB, vC, vD}

( ) 1 454 ( ) 619 ( ) 2 762

(

)

19×1.50 0

AN US

{vB, vC}

)

0 0 0 Total = )

127×1.50 Option 3 {vA, vB, vC}

(| |

( ) 511 659 1 311 2 481 $

CR IP T

(

Shipments Options Option 1 {vA, vB} {vA, vC} {vA, vD}

(

)

19×1.50

0 ( ) 100 Total = ( ) 100 0

Total = ( ) 200 Total =

1 335 1 592.50

2 927.50 $

1 582.5 619

2 201.5 $ 2 990.5 2 990.5 $

M

As it can be seen from this small real example, the use of decoupling points can decrease

3.2

ED

the cost from 2 481$ down to 2 201.50$, an 11.2% decrease.

Mathematical formulation

In order to present the mathematical model, we define the following parameters and

K clj

AC

clij

: Set of vehicles

CE

Parameters

PT

decision variables.

: Cost of an l-foot long load leaving the depot to decoupling point j : Cost of an l-foot long load leaving the decoupling point i and finishing its route at node j : Cost of each additional drop on a route/subroute

Q

: Maximum capacity of a vehicle, in linear-foot

dij

: Distance, in miles, between node i and node j



: The number of levels in the transportation cost function

12

ACCEPTED MANUSCRIPT

Decision variables xijk

: 1 if vehicle k travels from i to j where j may be a decoupling point i  N  {0}

yjk

jN

k K

: 1 if j is the decoupling point used by vehicle k. If route k is not followed jN

zijk

CR IP T

by a subroute, the last point of route k is a decoupling point. k K

: 1 if vehicle k travels from i to j

i  N is the decoupling point or follows it j  N follows the decoupling point jk

AN US

k K

: Binary variable set to 1 if j is the endpoint of subroute k, 0 otherwise. : Total demand of route k when it leaves the depot

: Demand of route k after the decoupling point (excluding the demand of the decoupling point)

: Used capacity on the vehicle k when leaving i to j. For this variable, n+1

M

uijk

is a copy of the depot.

j  N {n+1}

ED

i  N  {0} wlkj

k K

: Binary variable equal to 1 if the price interval l is used by route k

PT

starting from the depot and ending at j l = 1, …, 

AC

ijk

: Binary variable equal to 1 if the price interval l is used by subroute k

CE

̅̅̅̅̅̅

kK

starting at the decoupling point i and ending at j l = 1, …, 

kK

: Binary variable equal to 1 if subroute k starts at i and ends at j. i, j N

k K

Hence, variables yik are used to identify whether node i N is used as a decoupling point for vehicle k; arcs belonging to a route and up to a decoupling point are identified by the

13

ACCEPTED MANUSCRIPT

xijk variables, while arcs of the subroutes correspond to the zijk variables. The OVRP-DP can then be formulated as follows:

Minimize:

∑



(1a)

∑∑ ∑∑∑

CR IP T



∑

̅̅̅̅̅̅

(1b)

(∑ ∑(

)

∑

∑

(∑ ∑(

)

∑∑

∑∑∑

(

∑ (∑ ∑

∑

ED

∑

M

)

∑ ∑ *

+

AC ∑

(3)

∑∑

(6)

∑ ∑

(2)

(5)

* +

∑

)

∑

CE

∑ ∑

(1c)

(4)

PT

∑

)

∑

Subject to: ∑(

))

AN US

∑

(7a) ∑

(7b)

∑

(8a)

∑

(8b) (8c) 14

ACCEPTED MANUSCRIPT

∑ ∑

∑∑

(9)

* +

∑∑

(10) (

)

(11)

CR IP T

∑ ∑ ∑ ∑ ̅̅̅̅̅̅

(12)

∑ ∑ ̅̅̅̅̅̅

(13a)

∑ ∑ ̅̅̅̅̅̅

AN US

(13b)

∑ ̅̅̅̅̅̅ ∑ ∑

M

∑

∑

∑

∑

(15) (16)

(18)

+

(19a) (19b)

CE

*

(14)

(17)

PT

∑∑

ED

∑

(13c)

∑

(19c)

∑∑

AC

(20)

The first part of the objective function (1a) represents the load-based cost associated with each vehicle k leaving the depot with a load lying in the interval wlkj and ending at decoupling point j. Its second part is related to the load-based cost of the subroute of vehicle k starting from decoupling point i to endpoint j. Note that these two points are needed to identify the applicable transportation cost function. Part (1b) of the objective 15

ACCEPTED MANUSCRIPT

function computes the detour costs for each route or subroute k. In the first term the distance of all used arcs minus d0j is computed, and the second term has a similar function but starting from the decoupling point instead of the depot. The starting and ending points of subroute k are identified using variable

. Part (1c) of the objective

function computes the drop cost related to each additional stop. It does so by calculating

CR IP T

the number of arcs activated on the routes. The second part of (1c) computes the number of arcs in a subroute arcs and then subtracts the number of ending points.

Constraints (2) enforce that a vehicle travels from the depot to at most one node. A route visits a decoupling point only if the node is set as such by constraints (3). Constraints (4) impose capacity on the vehicles. Constraints (5) link the zijk with the yik variables. Each

AN US

customer j is visited by one vehicle as per constraints (6). Constraints (7a) and (7b) identify the decoupling point j on each route, yjk, and link the arcs before the decoupling point with variables xijk, and after the decoupling point with variables zjik. Also, constraints (7b) impose that if j is a decoupling point, it is set as the starting point of a subroute (

as a linearization of the product

). Constraints (9) compute the total demand of a route (with variables xijk)

M

(

). Constraints (8a)–(8c) set the value of

including the load of any subsequent subroutes (with variables zijk). Likewise, constraints

ED

(10) compute the total demand of a subroute. Constraints (11) identify the price interval l associated with the total demand of a route k ending at node j. If route k does not end at j

PT

(yjk=0), all wlkj will be set to zero in order to minimize their value in the objective function (1a). Constraints (12) determine the price interval of subroute k starting at a

CE

decoupling point i and ending at j via variables ̅̅̅̅̅̅. These variables are defined and set using constraints (13a)–(13c). Constraints (13b) determine the last node of subroute by identifying a node that is entered but never left. A combination of (13b) and (13c) allows

AC

the use of the right cost function for subroute k from i to j if i is a decoupling point (13a), j is the end of a subroute (13b) and both i and j are on the same subroute (13c). Constraints (14) identify a decoupling point as a node entered by a vehicle but never left. Constraints (15) manage the end of subroute k. Constraints (16)—(20) enforce the MillerTucker-Zemlin subtour elimination. Constraints (16) set the load when departing from the depot. By (17) uijk can have a non-zero value only if the related link xijk is used. Constraints (18) decrease the flow exiting j by its demand qj. By constraints (19a) to 16

ACCEPTED MANUSCRIPT

(19c) we ensure that we can go to the copy of the depot from i if and only if i is the end of a route (in case we have no subroute) or a subroute. Finally, by (20) vehicles returning to the copy of the depot are empty.

4

An Iterated Local Search algorithm

CR IP T

In this section we propose an Iterated Local Search algorithm tailored for the OVRP-DP. We first present the main procedures and phases of this algorithm in Sections 4.1–4.5, and we describe how they are used together within an ILS in Section 4.6. 4.1

Initialization

The Initialization phase is a random-based insertion procedure and is described in

AN US

Algorithm 1. It starts by randomly selecting the first client (line 3). At each step, a client will be randomly selected and inserted into its best position in the route or subroute of the corresponding vehicle, i.e., the position that leads to the minimum extra insertion cost (line 5). After each client insertion, the Decoupling point update procedure (line 6) is used to determine which client should be used as decoupling point in the considered

Decoupling point update procedure

ED

4.2

M

route/subroute. This procedure is explained next.

In a given route/subroute k having nk customers, each of these customers can be selected as the decoupling point leading to as many routes/subroutes combinations. This

CE

PT

procedure determines the best decoupling position to minimize function ( ). Algorithm 1: The Initialization Procedure

AC

0: Route index ; 1: While (there is an un-routed client) { 2: ; 3: Select a client randomly and visit it as the first client of route k; 4: while ((there is an un-routed client) and (the capacity of the working route is not violated)) { 5: Select a client randomly and insert that in the best feasible position of the kth route or subroute; 6: ( ) 7: } 8: }

17

ACCEPTED MANUSCRIPT

4.3

Swap procedure

We implement a classical Swap procedure which may be executed under a first or a best improvement strategy. In this procedure the goal is to improve the cost of the solution by exchanging two strings of consecutive customers. Let n1 and n2 be the number of

CR IP T

customers visited on two routes/subroutes. The procedure generates all possible combinations of k1 < n1 and k2 < n2 between each pair of routes/subroutes, and for each values of k1 and k2, identifies each possible string of k1 consecutive customers on the first route/subroute and of k2 consecutive customers on the second route/subroute, and swaps

AN US

them.

While evaluating these swaps, we also consider both orientations of each string, leading to four swapping possibilities. In this procedure, decoupling points may appear in any of the two strings. In order to evaluate the cost of the resulting routes/subroutes, we also reoptimize the selection of the decoupling points with the Decoupling point update

M

procedure. Note that when the procedure is applied to a single route/subroute, the length of the two strings are k1 and (n1- k1) with

. The procedure is repeated

4.4

ED

until no more improvements are found.

Extraction-Insertion procedure

PT

In this procedure the goal is to improve the solution by repositioning a sequence of consecutive customers. The length of this sequence is variable. The following

CE

possibilities are evaluated:

1- Inserting the extracted sequence between two visited customers of any route or

AC

subroute.

2- In case a route does not have a subroute, the sequence can be considered as a subroute, hence being inserted in the end of the route.

3- In case a route does not have a subroute, the sequence can be inserted at the beginning of the route, converting the original route into a subroute. 4- The extracted sequence can be used to start a new route.

18

ACCEPTED MANUSCRIPT

The original sequence and its opposite direction are considered, and all feasible insertion positions are evaluated and the best choice is implemented. If an extracted sequence cannot be inserted elsewhere improving the solution, it is inserted back into its original position and a new sequence is evaluated. This procedure iterates as long as the solution improves by extracting and inserting all

CR IP T

possible sequences. The pseudocode of this procedure is presented in Algorithm 2. Algorithm 2: The Extraction-Insertion procedure

AC

CE

PT

ED

M

AN US

1. Repeat { 2. ; (vehicle index) 3. While ( ){ 4. For ( ){ 5. If ( ) 6. 7. Else 8. 9. 10. 11. While ( ) { 12. While ( ) { 13. : Extract the sequence ( ) from the solution ( is the customer visited in the position). Let S’ be the reverted sequence ( ); 14. Try to insert and S’ back into the solution: 15. ● Allocate and S’ between two visited customers of the available routes or sub-routes 16. ● Allocate and S’ as the sub-route of an available route 17. ● Allocate and S’ as a route and transform the initial route as a subroute 18. ● Create a new route. 19. If (improvement found) 20. Implement the move leading to the greatest cost improvement and update the solution; 21. ( ) 22. Else 23. Insert S back in its original position; 24. } 25. 26. ;+ 27. } 28. } 29.Until (no more improvements are found) } 19

ACCEPTED MANUSCRIPT

4.5

Shaking procedures

In the following we present three shaking procedures used to perturb the solution and thus avoid local optima and explore large neighborhoods.

CR IP T

Shaking procedure 1: Two routes are randomly selected and all their clients are extracted and put into a set A of free clients. Then a client is randomly selected from A and inserted into its best position, i.e., the position with the minimum cost increase. If no feasible insertion position can be found, the customer is assigned to a new route. When all the

Decoupling Point Update procedure.

AN US

customers of A are reinserted, the best decoupling points are updated using the

Shaking procedure 2: Here, each sequence containing up to three consecutive customers is extracted from the solution. To do so, we first compute the cost savings obtained by extracting any sequence from the solution. As an example, suppose that we have the following route (

) and that we consider the 2-customer sequence (

). The

M

cost savings obtained by extracting the sequence (k, l) is equal to

ED

. After having calculated the cost savings for all possible sequences, the algorithm uses a roulette wheel mechanism to select which ones to extract. Thus, the probability of selecting sequence (k, l) is

PT

Sequences are selected until at least

⁄ where

is the sum of all savings.

percent of the customers have been removed. Their

CE

reinsertion is performed as in Shaking procedure 1. Shaking procedure 3: This procedure is based on the concept of regions such as those used in the sweep or petal algorithms (Renaud et al. 1996). For each customer i, we first

AC

calculate its polar angle

with respect to the depot considered as the center of the plane.

Then we compute what we call the vehicle angle

which corresponds to the average of

each customer’s polar angle visited by vehicle k. Without loss of generality, assume vehicles are relabeled in increasing order of their angle as follows.

20

. Then customers are removed

ACCEPTED MANUSCRIPT

We first remove all customers from vehicle (route) k, selected randomly, and put them in the set of extracted customers A. Then, more similar routes are removed: either route k-1 or k+1, or two other routes, either (k-2;k-1), (k-1;k+1) or (k+1;k+2) with the same probability. Extracted customers are added to A. Finally, all remaining customers, if any, which are located between the angles of the selected vehicles are also extracted from the

CR IP T

solution and added to the set A. The reinsertion of extracted customers is performed as previously described. 4.6

The overall ILS algorithm

The previous procedures are organized as follows to build up our overall ILS algorithm. It begins by generating a number  of solutions by using the Initialization, Swap and

AN US

Extraction – Insertion procedures. The  obtained solutions are stored in decreasing order of their cost in the ordered set BKS (lines 1-10 of Algorithm 3). In the second phase, the algorithm improves each of the  solutions from the BKS set. This is done repeatedly by randomly selecting a shaking procedure from Section 4.5 followed by the Swap and the

M

Extraction – Insertion procedures from Section 4.3 and 4.4. The search is stopped after MRep repetitions without improvement (lines 13 to 26). If a new best solution is found,

ED

BKS is updated and the search restarts with the best solution (lines 19 to 25). Otherwise the search is performed on the next best known solution. The complete search is stopped after having tried to improve the th best solution. After extensive preliminary testing,

PT

 = 5 and MRep = 100 were retained for leading to high quality solutions in short

5

CE

computing times. The ILS pseudocode is presented in Algorithm 3.

Computational results

AC

In this section we report the results of computational experiments run to evaluate the performance of the algorithm and to assess the value of decoupling points. The mathematical model implemented in Section 3 was coded in Visual Studio C++ within CPLEX version 12.3. The ILS algorithm described in Section 4 was coded in Visual Studio C++. All computations were executed on a machine equipped with a Dual core 2.6 GHz, having 2 GB of RAM. In Section 5.1 we evaluate the performance of our ILS on

21

ACCEPTED MANUSCRIPT

classical OVRP instances without decoupling points. The OVRP-DP instances are

AC

CE

PT

ED

M

AN US

CR IP T

described and evaluated in Section 5.2.

22

ACCEPTED MANUSCRIPT

Algorithm 3: The ILS algorithm best known solutions) (); Repeat {

(

)

CR IP T

(

Until (no improvements are obtained) } Update best known solution; (

){

( ( (

(

)

); ( )

(

;

); (

)

)) {

)

ED

M

} Update (

AN US

Repeat {

);

Until (no improvements are obtained for MRep consecutive iterations)}

Results on classical OVRP instances

CE

5.1

(the set of )*

(

PT

1. 2. 3. 4. 5. 6. 7. 8. 9. 10. } 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. }

The performance of the ILS was evaluated on classical OVRP benchmark instances taken from Christofides et al. (1979) and Fischer (1994). We use instances C1-C5, C11, C12,

AC

F11 and F12 as these do not contain route duration restrictions. Best known solutions for these instances are reported in Salari et al. (2010, Table 5). Details about these instances are given in Zachariadis and Kiranoudis (2010, Table 2). We report results for the hierarchical OVRP minimizing first the number of vehicles and then the routing cost. In Table 4 we first present for each instance (Inst.) the number of customers n, the number of vehicles, |K|, and the best known solution from the literature BKS. As we run

23

ACCEPTED MANUSCRIPT

the ILS five times, we present the Best solution, the Average solution, its Standard deviation and the Average computing time in seconds. When analyzing these solutions, one must carefully consider the number of vehicles used. As mentioned in Section 1, some papers minimize only the distance, while others minimize first the number of vehicles, and then the distance. We use the second approach. When comparing the results

CR IP T

summarized by Reinholz and Schneider (2013, Table 1) and by Salari et al. (2010, Table 5), and considering the minimum number of vehicles for each instance, our algorithm always provides the best solution, also improving the BKS for instance C5. The robustness of our algorithm can also be seen from the fact that its average solution was equal to the BKS for six instances, meaning that the optimum was found on each of the

AN US

five runs. Clearly, even if not originally designed for solving the OVRP, the proposed ILS performs as well as the best available algorithms for this problem. We note that Soto et al. (2017, Table 1) reported lower solution values for instances C1 (412.95), C2 (564.06), C3 (639.25) and C5 (868.11). However these results were obtained with more vehicles, respectively 6, 11, 9 and 17. For instances with the same number of vehicles, we

M

obtained the same results for C4, F11 and F12 and a better solution for C11.

C1 C2 C3 C4 C5 C11 C12 F11 F12

50 75 100 150 199 120 100 71 134

|K| 5 10 8 12 16 7 10 4 7

AC

ILS results

BKS

Best solution

Average solution

Standard deviation

Average time (s)

416.06 567.14 639.74 733.13 879.37 682.12 534.24 177.00 769.55

416.06 567.14 639.74 733.13 878.94* 682.12 534.24 177.00 769.55

416.58 567.14 639.88 733.13 880.46 682.12 534.24 177.00 769.55

0.71 0.00 0.30 0.00 2.12 0.00 0.00 0.00 0.00

16.4 13.6 193.6 233.4 178.0 165.0 153.6 135.3 212.9

PT

n

CE

Inst.

ED

Table 4. Performance of the ILS on classical hierarchical OVRP instances

*New improved solution

24

ACCEPTED MANUSCRIPT

5.2

Experiments on the OVRP-DP

Figure 1 showed the distribution network of our partner for which we have all the relevant real data. The main distribution center (which is also a production facility) is located in Québec city, eastern Canada (in green). From this DC, 32 distribution centers for two major retailers (in blue and red) spread throughout all the USA are served. Each

CR IP T

location is identified by its latitude and longitude coordinates, and road distances computed in miles. An additional drop cost is set to 160 $. Transportation costs from Quebec to each city and for each load quantity (in linear-foot) are available in Canadian Dollars. These costs were effective July 2016, obtained from the private carrier of our partner. We obtained from Canadian and American carriers working with our partner

AN US

their complete cost structure for each combination of cities. When such information was not completely available, we requested the price of a full load, which was then broken down into one-feet intervals to obtain a complete cost structure.

From this real commercial supply chain, we have generated 77 instances. Each instance is

M

labeled "Ins-x-y" in which "x" represents the number of customers and "y" the instance number. We generated 20 instances having 8, 10 and 12 customers where for each

ED

instance the customer locations are selected randomly within the set of 32 clients’ DC. We also generated 12 instances containing 16 customers and 5 instances with 20 customers. As an example of the required data, for a 16 customers instance we have 18

PT

different matrices: one 17×17 global distance matrix (represented by the dij parameter), one 16×50 (as we have 50 cost levels) cost matrix from the depot to all customers

CE

(represented by the clj parameter), and 16 other cost matrices of dimension 15×50 for the transportation costs between the 16 customers (represented by the clij parameter). These

AC

instances are available upon request. Table 5 presents the results for 40 small OVRP-DP instances containing up to 10 customers. Column OVRP was obtained by using the following procedures: Initialization, Swap and Extraction and Insertion. These procedures were applied without the Decoupling point update. In this case, the last visited customer is set as the decoupling point. The OVRP+ column was obtained as the OVRP one but the Swap and the Extraction and Insertion procedures now use the Decoupling point update procedure. 25

ACCEPTED MANUSCRIPT

Then each OVRP-DP model was solved two times with CPLEX. CPLEX was first run using its default parameter settings, and then using the relaxation induced neighborhood search more often than the default value (the RINS parameter was set to 1 instead of -1). A maximum computing time of 10 800 seconds was allowed for each run. We use superscripts D,

R

or

=

to indicate whether the best solution was found with the default

CR IP T

parameters, the RINS parameters, or if they provide the same solution, respectively. The last four columns show the results of the complete ILS which was run five times on each instance. ILS results are presented under the columns Average, Standard deviation, Best solution and Worst solution. Comparing the OVRP solution with the OVRP+ one shows the cost reduction achievable by considering the use of decoupling points. For the 40

AN US

instances, the average cost of the OVRP is 12 109.57 $ while the average of the best ILS solution is 8 882.11 $, a 36.33% reduction. The average OVRP+ solution is 11 325.89 $, meaning that optimizing the location of the decoupling points (without the ILS) leads to an average reduction of about 2 444 $. Finally, the average best ILS solution is 8 882.10 $, a further improvement of 27.50% over the OVRP+.

M

Computing times for CPLEX are not reported as they always reached the maximum allotted time of 10 800 seconds and none of the 40 solutions were proven optimal. On

ED

these instances, a better solution was obtained two times with the default parameters, six times with the RINS parameters, and other 12 solutions were the same. We also note that

PT

the best solution reported by the ILS was always better or equal to that of CPLEX (28 equal solutions, 12 better solutions provided by the ILS algorithm). On these 40 instances

CE

CPLEX reported solution costs on average 0.53% higher than those of the ILS. These results show that the performance of the ILS is very stable on these instances as the same solution was obtained on each of the five runs leading to a standard deviation of 0.00.

AC

Computing times for the ILS are not reported as the longest run was under 0.2 seconds. Table 6 shows the results for larger instances. On these 37 instances the ILS yields an average cost reduction of 25.37% over the OVRP results and of 18.20% over the OVRP+ ones. This clearly shows the importance of considering the decoupling points at each stage of the solution process within an ILS. Again, our ILS algorithm proves its stability by producing the same solution over the five runs for 30 of the 37 instances. As expected,

26

ACCEPTED MANUSCRIPT

the quality of CPLEX solutions quickly deteriorates, with an average deviation over the ILS of 1.98% for instances with 12 customers and 10.66% for those with 16 customers. For instances with 20 customers, CPLEX, with both sets of parameters, always ran out of memory before the time limit of 10 800 seconds. The average gap obtained by CPLEX was about 127% with the default parameters and 134% with the RINS parameters. The

CR IP T

ILS used less than 0.5 seconds for the 16 customer instances, and up to 1.12 seconds for the 20 customer instances.

6

Conclusion

In this paper we have described a new practical vehicle routing problem where multiple carriers are used to deliver goods over large territories. Based on a collaboration with

AN US

industrial partners, we have modeled this problem as an OVRP and we have introduced the notion of decoupling points, which generalizes the classical OVRP. Here, a route can be performed by more than one carrier. We have modeled this problem mathematically and designed an efficient ILS algorithm for its solution. Even if not specifically designed

M

to solve the OVRP, our ILS produces the best known results on a set of classical OVRP instances, also giving one new best solution. Based on real orders and real transportation costs obtained from our industrial partners we generated several real-life instances. By

ED

comparing results between the OVRP and the OVRP-DP we have evaluated the benefits of using more than only one carrier per route. Extensive computational results

PT

demonstrate the efficiency of the developed ILS and the advantages of using the new

AC

CE

decoupling point concept.

27

ACCEPTED MANUSCRIPT

Table 5. Performance of the ILS on small OVRP-DP instances ILS

AC

9658.58 8189.45 7269.21 10084.77 11328.56 9448.72 12192.61 10353.58 9931.62 9914.25 11576.31 11610.43 10636.99 9590.19 8544.13 5746.46 10203.25 8925.86 11120.30 9618.79 9797.20 12750.36 15442.80 16624.93 11246.99 11632.79 12504.06 11768.90 11778.94 14607.08 13059.14 13521.55 12966.32 11402.01 11537.24 12533.98 13461.52 14314.01 12040.92 9367.92 14530.05 12854.57 11325.89

7217.75R 5906.50= 6317.91= 6085.95D 7247.23= 8033.49R 9678.30= 7033.87= 6777.05= 8534.46R 9712.13R 10433.48= 8387.98= 6954.41= 6068.78D 5036.01= 7985.56R 7775.44= 8561.79R 7731.92= 7 594.28 9948.11= 11238.99= 12471.18D 7740.60= 9482.24D 10693.58R 8699.47D 10473.72D 11678.41D 10791.52= 10917.41R 9493.00D 8262.46D 8970.02D 10176.23D 11236.34D 12110.43R 10791.19D 8030.36D 12571.16R 10307.30 8950.78

Standard deviation 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 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

AN US

10398.80 8189.45 7545.80 10617.80 12006.30 11471.45 13895.95 10806.35 10464.65 10230.65 13279.65 13568.70 11392.65 9980.15 9058.65 6348.25 10492.35 9771.70 11120.30 9695.15 10516.73 12796.25 15442.80 16680.10 11409.90 13202.90 15121.80 12874.55 14065.35 15967.20 14817.65 13521.55 14309.65 11893.80 12004.80 13093.05 14567.65 14836.55 12545.65 9626.00 15270.75 13702.40 12109.57

Average solution 7217.75 5906.50 6317.91 6085.95 7247.23 8033.49 9678.30 6710.94 6777.05 8534.46 9707.34 10311.63 8387.98 6954.41 6068.78 5036.01 7985.56 7775.44 8545.15 7731.92 7 550.69 9948.11 11238.99 12471.18 7740.60 9480.21 10199.66 8699.47 10473.72 11598.58 10791.52 10592.16 9182.21 8114.91 8970.02 10171.08 11236.34 11969.05 10791.19 8030.36 12571.16 10213.53 8882.11

ED

28

Best solution 7217.75 5906.50 6317.91 6085.95 7247.23 8033.49 9678.30 6710.94 6777.05 8534.46 9707.34 10311.63 8387.98 6954.41 6068.78 5036.01 7985.56 7775.44 8545.15 7731.92 7 550.69 9948.11 11238.99 12471.18 7740.60 9480.21 10199.66 8699.47 10473.72 11598.58 10791.52 10592.16 9182.21 8114.91 8970.02 10171.08 11236.34 11969.05 10791.19 8030.36 12571.16 10213.53 8882.11

Worst solution 7217.75 5906.50 6317.91 6085.95 7247.23 8033.49 9678.30 6710.94 6777.05 8534.46 9707.34 10311.63 8387.98 6954.41 6068.78 5036.01 7985.56 7775.44 8545.15 7731.92 7 550.69 9948.11 11238.99 12471.18 7740.60 9480.21 10199.66 8699.47 10473.72 11598.58 10791.52 10592.16 9182.21 8114.91 8970.02 10171.08 11236.34 11969.05 10791.19 8030.36 12571.16 10213.53 8882.11

CR IP T

CPLEX

M

OVRP

CE

Data-8-1 Data-8-2 Data-8-3 Data-8-4 Data-8-5 Data-8-6 Data-8-7 Data-8-8 Data-8-9 Data-8-10 Data-8-11 Data-8-12 Data-8-13 Data-8-14 Data-8-15 Data-8-16 Data-8-17 Data-8-18 Data-8-19 Data-8-20 Avg Data-10-1 Data-10-2 Data-10-3 Data-10-4 Data-10-5 Data-10-6 Data-10-7 Data-10-8 Data-10-9 Data-10-10 Data-10-11 Data-10-12 Data-10-13 Data-10-14 Data-10-15 Data-10-16 Data-10-17 Data-10-18 Data-10-19 Data-10-20 Avg Global Avg

OVRP

PT

Inst.

+

ACCEPTED MANUSCRIPT

Table 6. Performance of the ILS on large OVRP-DP instances ILS

CE

AC

16097.29D 12702.81D 10598.61D 13318.16R 11128.72R 12635.48D 12752.08D 11172.02D 12019.15D 10317.21D 14292.42R 12967.57R 10444.49D 12934.68D 14849.24R 11072.78D 10483.78D 12925.04D 14121.05D 10881.32D 12 395.98 18707.45= 14900.38R 18263.45= 16089.00= 17705.85= 16581.83D 18892.14D 16875.45= 19031.67D 18598.65D 18310.10= 16585.30= 17 829.24 OOM= OOM= OOM= OOM= OOM=

Average solution 15958.22 12532.59 10471.45 13233.88 10995.80 12328.10 12535.36 10839.77 11859.65 9752.90 14062.00 12113.29 10169.29 12778.30 14725.29 11007.38 10434.67 12771.30 14010.53 10690.38 12163.51 17142.08 13844.20 16658.20 13852.37 16490.70 15069.06 16701.91 15417.20 16216.24 16854.63 17110.19 15485.83 15903.55 20210.52 20953.63 18299.96 20145.29 18575.93 19637.06 14386.43

Std. deviation 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 6.72 57.78 0.00 0.00 0.00 0.00 3.22 0.00 0.00 0.00 0.00 6.82 0.00 0.00 197.28 72.00 120.00 0.00 0.00 33.00 0.00 100.04 0.00 0.00 0.00 20.00 15.15

29

Best solution 15958.22 12532.59 10471.45 13233.88 10995.80 12328.10 12535.36 10839.77 11859.65 9752.90 14062.00 12113.29 10169.29 12778.30 14713.95 10838.42 10434.67 12771.30 14010.53 10690.38 12 154.49 17142.08 13844.20 16658.20 13852.37 16461.71 15069.06 16701.91 15085.13 16192.85 16737.67 17110.19 15485.83 15 861.76 20210.52 20788.10 18299.96 20145.29 18575.93 19 603.96 14 636.54

Worst solution 15958.22 12532.59 10471.45 13233.88 10995.80 12328.10 12535.36 10839.77 11859.65 9752.90 14062.00 12113.29 10169.29 12778.30 14729.07 11026.15 10434.67 12771.30 14010.53 10690.38 12164.64 17142.08 13844.20 16658.20 13852.37 16492.23 15069.06 16701.91 15573.32 16426.78 16971.60 17110.19 15485.83 15 943.98 20210.52 21022.03 18299.96 20145.29 18575.93 19 650.75 14 402.00

CR IP T

CPLEX

AN US

19237.86 13686.37 13172.17 16425.25 14476.48 14953.07 14691.00 12556.79 14526.52 11301.03 18230.91 15871.86 12394.71 14614.81 16457.37 13420.39 13009.79 15306.81 18290.82 13622.09 14812.30 20241.02 17491.64 21073.84 14992.30 19714.34 18087.67 19926.10 17502.95 20800.72 19390.42 20878.39 18635.45 19061.23 24632.43 22787.44 22338.66 23725.63 21667.80 23030.39 17300.89

PT

20623.65 Data-12-1 Data-12-2 13686.85 Data-12-3 13943.00 Data-12-4 17165.95 Data-12-5 17002.35 Data-12-6 15571.00 Data-12-7 15241.35 Data-12-8 12716.50 Data-12-9 15017.90 Data-12-10 11650.20 Data-12-11 18459.95 Data-12-12 16046.30 Data-12-13 13007.60 Data-12-14 16318.15 Data-12-15 18398.25 Data-12-16 13959.25 Data-12-17 13530.00 Data-12-18 17170.15 Data-12-19 20579.70 16230.40 Data-12-20 15815.92 Avg 22331.45 Data-16-1 Data-16-2 18332.60 Data-16-3 21444.20 Data-16-4 15665.40 Data-16-5 20284.80 Data-16-6 18739.05 Data-16-7 21135.95 Data-16-8 19013.20 Data-16-9 21718.70 Data-16-10 21040.75 Data-16-11 21675.90 Data-16-12 18951.85 20027.82 Avg 26482.70 Data-20-1 Data-20-2 24588.35 Data-20-3 23348.55 Data-20-4 24656.15 23248.10 Data-20-5 24464.77 Avg 18350.70 Global Avg OOM : Out of memory

OVRP

M

OVRP

ED

Inst.

+

ACCEPTED MANUSCRIPT

Acknowledgements This research was partly supported by grant numbers 2014-05764 and 0172633 from the Canadian Natural Sciences and Engineering Research Council. We thank our industrial partner for providing us with their orders and customer locations and to two Canadian carriers for giving us access to their LTL transportation costs through major Canadian

CR IP T

and US cities. We thank Marc Sevaux and Maria Soto for their information about their results. We also thank an associate editor and two anonymous referees for their valuable comments on an earlier version of this paper.

References

AN US

D. Aksen, Z. Özyurt & N. Aras, Open vehicle routing problem with driver nodes and time deadlines. Journal of the Operational Research Society, 2007, 58, 1223-1234. T. H. Ali, S. Radhakrishnan, S. Pulat & N. C. Gaddipati, Relay network design in freight transportation

systems.

Transportation

Research

Part

E:

Logistics

and

Transportation Review, 2002, 38, 405-422.

M

J. Brandão, A tabu search for the open vehicle routing problem. European Journal of

ED

Operational Research, 2004, 157, 552-564. N. Christofides, A. Mingozzi & P. Toth, The vehicle routing problem. In: Christofides N., Mingozzi A., Toth P., Sandi C., editors. Combinatorial Optimization. Chichester:

PT

Wiley: 1979, 313-338.

CE

G. Clarke & J. W. Wright, Scheduling of vehicle from a central depot to a number of delivery points. Operations Research, 1964, 12, 568-581.

AC

L. C. Coelho, J. Renaud & G. Laporte, Road-based goods transportation: a survey of realworld logistics applications from 2000 to 2015. INFOR: Information Systems and Operations Research, 2016, 54, 79-96.

U. Derigs & K. Reuter, A simple and efficient tabu search heuristic for solving the open vehicle routing problem, Journal of the Operational Research Society, 2009, 60, 1658-1669.

30

ACCEPTED MANUSCRIPT

M. Fisher, Optimal solutions of vehicle routing problems using minimum k-trees. Operations Research, 1994, 42, 626-642. K. Fleszar, I. H. Osman & K. S. Hindi A variable neighbourhood search algorithm for the open vehicle routing problem, European Journal of Operational Research, 2009,

CR IP T

195, 803-809. Z. Fu, R. W. Eglese & L. Li, A new tabu search heuristic for the open vehicle routing problem. Journal of the Operational Research Society, 2005, 56, 267-274.

M. Gendreau, A. Hertz & G. Laporte, New insertion and post-optimization procedures for the traveling salesman problem. Operations Research, 1992, 40, 1086-1094.

AN US

S. Lapierre, A. B. Ruiz, & P. Soriano, Designing distribution networks: Formulations and solution heuristic. Transportation Science, 2004, 38, 174-187.

G. Laporte, Fifty years of vehicle routing. Transportation Science, 2009, 43, 408-416. A. N. Letchford, J. Lysgaard & R. W. Eglese, A branch-and-cut algorithm for the

M

capacitated open vehicle routing problem. Journal of the Operational Research Society, 2007, 58, 1642-1651.

ED

F. Li, B. L. Golden & E. A. Wasil, The open vehicle routing problem: Algorithms, large scale test problems, and computational results. Computers & Operations Research,

PT

2007, 34, 2918-2930.

A. D. López-Sánchez, A. G. Hernández-Díaz, D. Vigo, R. Caballero & J. Molina, A

CE

multi-start algorithm for a balanced real-world open vehicle routing problem. European Journal of Operational Research, 2014, 238, 104-113.

AC

A. Reinholz, Optimizing logistics networks with metaheuristics. In Optimization of logistics systems, Methods and Experiences, 85-98. Symposium of the Collaborative Research Center 559 ―Modeling of large logistics networks‖, Verlag Praxiswissen, Fraunhofer IML, Dortmund, DE, 2008. A. Reinholz & H. Schneider, A hybrid (1 + 1)-evolutionary strategy for the open vehicle routing problem. In L. Di Gaspero, A. Schaerf, T. Stützle (Eds.), Advances in 31

ACCEPTED MANUSCRIPT

metaheuristics, Operations Research/Computer Science Interfaces Series 53, Springer, New York, 2013, 127-141. J. Renaud, F. F. Boctor & G. Laporte, An improved petal heuristic for the vehicle routeing problem. Journal of the Operational Research Society, 1996, 47, 329-336.

CR IP T

P. P. Repoussis, C. D. Tarantilis & G. Ioannou, The open vehicle routing problem with time windows. Journal of the Operational Research Society, 2007, 58, 355-367.

P. P. Repoussis, C. D. Tarantilis, O. Bräysy & G. Ioannou A hybrid evolution strategy for the open vehicle routing problem. Computers & Operations Research, 2010, 37, 443-455.

AN US

M. Salari, P. Toth, & A. Tramontani, An ILP improvement procedure for the open vehicle routing problem. Computers & Operations Research, 2010, 51, 2106-2120. D. Sariklis & S. Powell, A heuristic method for the open vehicle routing problem. Journal of the Operational Research Society, 2000, 51, 564-573.

M

L. Schrage L., Formulation and structure of more complex/realistic routing and scheduling problems. Networks, 1981, 11, 229-232.

ED

M. Soto, M. Sevaux, A. Rossi & A. Reinholz, Multiple neighborhood search, tabu search and ejection chains for the multi-depot open vehicle routing problem. Computers &

PT

Industrial Engineering, 2017, 107, 211-222. C. D. Tarantilis, D. Diakoulaki & C. T. Kiranoudis, Combination of geographical

CE

information system and efficient routing algorithm for real life distribution operations. European Journal of Operational Research, 2004a, 152, 437-453.

AC

C. D. Tarantilis, G. Ioannou, C. T. Kiranoudis, & G. P. Prastacos, A threshold accepting approach to the open vehicle routing problem. RAIRO, 2004b, 38, 345-360.

C. D. Tarantilis, C. T. Kiranoudis, G. Ioannou, & G. P. Prastacos, Solving the open vehicle routing problem via a single parameter metaheuristic algorithm. Journal of the Operational Research Society, 2005, 56, 588-596.

32

ACCEPTED MANUSCRIPT

P. Toth & D. Vigo, Vehicle routing. Problems, methods and applications. P. Toth & D. Vigo editors, MOS-SIAM Series on Optimization, 2014. H. Üster H. & P. Kewcharoenwong, Strategic design and analysis of relay network in truckload transportation. Transportation Science, 2011, 45, 505-523.

CR IP T

C. T. Zachariadis & C. T. Kiranoudis, An open vehicle routing problem metaheuristic for examining wide solution neighborhoods. Computers & Operations Research, 2010,

AC

CE

PT

ED

M

AN US

37, 716-723.

33