The Flexible Periodic Vehicle Routing Problem

Accepted Manuscript

The Flexible Periodic Vehicle Routing Problem Claudia Archetti, Elena Fernandez, Diana L. Huerta-Munoz ´ ˜ PII: DOI: Reference:

S0305-0548(17)30073-4 10.1016/j.cor.2017.03.008 CAOR 4216

To appear in:

Computers and Operations Research

Received date: Revised date: Accepted date:

28 June 2016 16 March 2017 17 March 2017

Please cite this article as: Claudia Archetti, Elena Fernandez, Diana L. Huerta-Munoz, The Flex´ ˜ ible Periodic Vehicle Routing Problem, Computers and Operations Research (2017), doi: 10.1016/j.cor.2017.03.008

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 • A Mathematical model for the Flexible Periodic Vehicle Routing Problem (FPVRP) isproposed. • A Load-based formulation is developed to solve the FPVRP. • A set of several inequalities are proposed to reinforce the load-based FPVRPformulation.

CR IP T

• A Load-based formulation for the Periodic Vehicle Routing Problem (PVRP) and a load-based formulation for the Flexible Periodic Vehicle Routing Problem with InventoryConstraints (FPVRP-IC) are proposed for a fair comparison in the computationalexperience.

• A computational experience was carried out in order to show that FPVRP formulationoutper-

AC

CE

PT

ED

M

AN US

forms PVRP and FPVRP-IC formulations in terms of solution quality.

1

CR IP T

ACCEPTED MANUSCRIPT

The Flexible Periodic Vehicle Routing Problem

AN US

Claudia Archetti, Elena Fernández, Diana L. Huerta-Muñoz

Periodic Problem

IRP

PVRP

- No limited by a schedule

AC

CE

PT

ED

M

- Limited by a schedule

- No fixed # of visits

- Fixed # of visits

- Fixed quantity to deliver

+ Flexibilty

- Quantity to deliver limited by an inventory level

- No limited by a schedule

FPVRP

- No fixed # of visits - No inventory levels

$$$$$

Optimal Solution

CR IP T

ACCEPTED MANUSCRIPT

The Flexible Periodic Vehicle Routing Problem Claudia Archetti

Department of Economics and Management

AN US

Università di Brescia Brescia, Italy

[email protected] Elena Fernández

M

Department of Statistics and Operations Research Universitat Politècnica de Catalunya Barcelona, Spain

ED

[email protected]

Diana L. Huerta-Muñoz

1

PT

Department of Statistics and Operations Research

AC

CE

Universitat Politècnica de Catalunya

1

Corresponding author

Barcelona, Spain [email protected] Submitted: June 2016 Accepted: March 2017

Abstract

CR IP T

ACCEPTED MANUSCRIPT

This paper introduces the Flexible Periodic Vehicle Routing Problem (FPVRP) where a carrier has to establish a distribution plan to serve his customers over a planning horizon. Each customer has a total demand that must be served within the horizon and a limit on the maximum quantity

AN US

that can be delivered at each visit. A fleet of homogeneous capacitated vehicles is available to perform the services and the objective is to minimize the total routing cost. The FPVRP can be seen as a generalization of the Periodic Vehicle Routing Problem (PVRP) which instead has fixed service frequencies and schedules and where the quantity delivered at each visit is fixed. Moreover, the FPVRP shares some common characteristics with the Inventory Routing Problem (IRP) where inventory levels are considered at each time period and, typically, an inventory cost is involved in

M

the objective function. We present a worst-case analysis which shows the advantages of the FPVRP with respect to both PVRP and IRP. Moreover, we propose a mathematical formulation for the

ED

problem, together with some valid inequalities. Computational results show that adding flexibility

PT

improves meaningfully the routing costs in comparison with both PVRP and IRP.

AC

CE

Keywords: Periodic Vehicle Routing; Inventory Routing; Service Frequency; Flexibility.

ACCEPTED MANUSCRIPT

1

Introduction

Vehicle Routing Problems (VRPs) (Toth and Vigo, 2014) deal with the definition of a distribution plan where a fleet of vehicles is used to deliver (or pick-up) goods to a set of customers. These problems involve a single-period plan, meaning that the delivery operations start at the beginning of the period and should terminate within the end of it when all customers are served. However,

CR IP T

several real–world applications deal with periodic delivery operations over a given time horizon. For example, periodic delivery service may be found in delivery of groceries or collection of waste. We refer to Francis et al. (2008) for a more exhaustive list of periodic service applications. The Periodic Vehicle Routing Problem (PVRP) (Campbell and Wilson, 2014) is a vehicle routing problem where the service to the customers has to be provided over multiple periods. It is assumed that customers must be served with a certain frequency and according to a given schedule and they should receive

AN US

a fixed quantity at each visit. The problem is to choose the visit schedule for each customer and to organize vehicle routes. Another broad class of routing problems dealing with multi-period plans is the class of Inventory Routing Problems (IRPs) (Bertazzi and Speranza, 2012, 2013; Coelho et al., 2013). IRPs differ from PVRPs in that they do not fix the service frequency nor the quantity delivered to each customer at each visit. Instead, each customer has a predefined rate of goods consumption per period, and the inventory level is calculated at each time period. The distribution

M

plan has to be such that each customer must be able to satisfy his consumption rate at each time period.

ED

In this paper, we present a new problem which deals with periodic demand without resorting to inventories. We call this problem the Flexible Periodic Vehicle Routing Problem (FPVRP). In the FPVRP, customers are assigned to a total demand that has to be satisfied within the end of

PT

the planning horizon. The quantity delivered at each visit should not exceed the customer capacity, which is typically lower than the total demand. Thus, multiple visits to each customer are performed over the planning horizon. The FPVRP can be seen as an extension of the PVRP where no fixed

CE

frequency schedule is set. Instead, if the customer capacity corresponds to the ratio between the total demand and the frequency established in the PVRP, then the FPVRP becomes a generalization

AC

of the PVRP where each customer must be visited a number of times which is at least equal to the PVRP frequency. Moreover, the quantity delivered at each visit has to be defined. This clearly increases the flexibility with respect to the PVRP setting and may produce cost savings. The FPVRP is also related to the IRP, where no fixed frequency is set for the customer visits

and the delivered quantity is a decision variable. However, contrary to the IRP, no inventory level is considered in the FPVRP. Again, this gives the FPVRP a higher flexibility with respect to the IRP which may result in cost savings. The main contributions of the paper are the following. We introduce the FPVRP, which is a new and challenging problem dealing with flexibility in periodic delivery operations. We show that, 1

ACCEPTED MANUSCRIPT

theoretically, the FPVRP can produce arbitrarily large improvements with respect to both the PVRP and the IRP when considering the routing cost. Indeed, such improvements come at the expenses of introducing additional decisions to the problem related to defining when the customers should be visited and how much should be delivered at each visit. This clearly increases the complexity of the problem, as will become evident in the computational results. We then propose a mathematical formulation of the FPVRP based on commodity flow variables, which we call load-based formulation.

CR IP T

This choice is consistent with recent works on other VRP variants for which load-based formulations have been proposed (Letchford and Salazar-González, 2015; Archetti et al., 2014). We also present a load-based formulation for the PVRP and for the IRP, which we use for the computational experiments, as they are quite effective and they do not require the use of sophisticated techniques (like branch-and-cut or column generation) for their implementation.

AN US

The paper is organized as follows. In Section 2, a review of the related literature is presented. Then, Section 3 provides a formal definition of the FPVRP. The worst case ratio between the optimal values of the PVRP and the FPVRP and the optimal values of the IRP and the FPVRP are analyzed, respectively, in Section 4. Load-based mathematical programming formulations for the FPRVP, the IRP, and for the PVRP are proposed in Section 5. In order to strengthen the formulation and improve the computing times, several families of valid inequalities are introduced

M

for the FPVRP in Section 6. Computational experiments and the analysis of numerical results are

2

ED

shown in Section 7. Finally, some concluding remarks and future work are given in Section 8.

Literature Review

PT

Vehicle Routing Problems (Toth and Vigo, 2014) are focused on finding an optimal design of routes in which a fleet of vehicles departs from a given depot to satisfy the demand of a set of customers.

CE

VRPs were introduced by Dantzig and Ramser (1959) and are NP–Hard (Garey and Johnson, 1979) because they generalize the well-known Traveling Salesman Problem (TSP) (Shmoys et al., 1985). Recently, different variants have been proposed to solve real–world applications. The book by Toth

AC

and Vigo (2014) overviews a good number of them. Specific applications require to handle periodic delivery operations and service flexibility to attend the customers demands over a given time horizon.

The two main classes of routing problems dealing with periodic delivery operations are PVRPs and IRPs.

PVRPs were introduced by Beltrami and Bodin (1974), a formal definition and the first mathematical formulation were given by Russell and Igo (1979) and Christofides and Beasley (1984), respectively. Different studies in which these problems are addressed can be found in Francis et al. (2006), Baldacci et al. (2011), and Hemmelmayr et al. (2009). Recent studies on the PVRP are mainly devoted to the development of heuristic algorithms (Vidal et al., 2012; Cacchiani et al., 2014; 2

ACCEPTED MANUSCRIPT

Michallet et al., 2014; Escobar et al., 2014; Rahimi-Vahed et al., 2015). For a recent survey on the PVRP the reader is refereed to Campbell and Wilson (2014). The IRP was introduced by Bell et al. (1983). An extensive review of the IRP literature, including different inventory policies, is given in Coelho et al. (2013) and in two recent tutorials (Bertazzi and Speranza, 2012, 2013). Apart from the two main classes of problems discussed above, there are many recent contributions

CR IP T

dealing with problems which extend or combine the IRP and the PVRP to deal with flexibility issues. Some related literature is discussed in the following section.

Finally, we note that the FPVRP, the PVRP and the IRP share the following common characteristic: customers can be visited more than once in the optimal solution. In the three above mentioned problems, multiple visits to a customer are due to the need of visiting a customer more than once in the considered planning horizon. However, it may be convenient to have multiple visits

AN US

to customers also when considering a single period, as shown by the studies on the Split Delivery Vehicle Routing Problem (SDVRP). SDVRP is a generalization of the classical VRP where multiple visits to customers are allowed. It has attracted a lot of attention in the recent years. For a recent survey on the problem the reader is refereed to Archetti and Speranza (2012).

Flexibility in Periodic Routing Problems

M

2.1

There are several works where service flexibility is studied. The first one we report is strictly related to the PVRP and is called the PVRP with Service Choice (PVRP–SC). The PVRP–SC is a variant

ED

of the PVRP in which the visit frequency is a decision variable. In particular, it allows to visit customers more often than their predefined frequency or, at least, at the same frequency. The amount of product to deliver at each visit depends on the assigned schedule and the selected delivery strategy.

PT

A profit that affects the objective function is associated with the visit frequency decision variable. The aim is to minimize the total travel cost minus the service benefit. An Integer Programming (IP)

CE

formulation for this problem was given by Francis et al. (2006) and a continuous approximation model was provided by Francis and Smilowitz (2006). Table 1 shows the main differences among the

AC

PVRP, the PVRP–SC and the IRP.

3

ACCEPTED MANUSCRIPT

Problem

Periodicity

PVRP

Predefined schedule

Same quantity at each Depends on the selected

frequency modeled as a

schedule and the delivery

decision variable

strategy at each customer

No predefined schedule

Modeled as a decision

and unconstrained

variable. Depends on the

number of visits

replenishment policy

IRP

Minimize routing cost

visit

Predefined schedule, visit

SC

Objective

Minimize the total travel

CR IP T

PVRP-

Delivered Quantity

cost minus service benefit Minimize holding and routing costs

AN US

Table 1: Comparisons among PVRP, PVRP–SC and IRP

As can be seen in the table and as stated by Francis et al. (2006) and Francis and Smilowitz (2006), the PVRP-SC relaxes the fixed frequency constraint of the PVRP and thus allows a higher degree of flexibility. The FPVRP, which is the focus of this study, further increases the degree of flexibility by relaxing the fixed schedule and the fixed quantity constraint. They can thus be

M

seen as different strategies to handle periodic routing problems. Concerning the IRP, the degree of flexibility is similar to the one of the FPVRP, i.e., no fixed frequency or schedule or quantity.

ED

However, inventory constraints have to be taken into account. Additional studies related to periodicity in routing problem are analyzed in the following. Rusdiansyah and Tsao (2005) integrate IRP and PVRP with Time Windows (PVRPTW)

PT

features to solve a real–world application in the field of delivery of products in vending–machine supply chains. They develop a mathematical formulation that combines both inventory and periodic

CE

routing, the flexibility criterion is applied by considering service frequency as a decision variable. The objective function attempts to minimize routing, holding and visit frequency costs. They propose different heuristic algorithms for the solution of the problem. Francis et al. (2007) develop

AC

a Tabu Search (TS) heuristic for the PVRP–SC that incorporates different measures to analyze the trade-off between operational flexibility and operational complexity (the ability to modify operating conditions). Hemmelmayr et al. (2009) propose several solution approaches based on IP and Variable Neighborhood Search (VNS) to evaluate delivery strategies for blood products supplies

which combine IRP and PVRP features. Pacheco et al. (2012) present a Mixed–Integer Linear Programming (MILP) formulation and a two–phase solution method, which allow flexibility in the delivery dates, for a real–world problem of a bakery company in Spain. Aksen et al. (2012) propose two different MILP formulations that combine PVRP and IRP characteristics for a waste vegetable oil collection problem. The aim is to minimize the sum of collection, inventory and purchasing costs. 4

ACCEPTED MANUSCRIPT

Archetti et al. (2015) study a routing problem with due dates where flexibility is related to the time at which each customer can be visited. Flexibility in routing problems has been analyzed also concerning different issues with respect to frequency or delivery schedule. Some examples are provided by the works in Hashimoto et al. (2006), where delivery in time windows is analyzed, García Calvillo (2010), where the focus is on the analysis of the flexibility on delivery dates, and Há et al. (2014) who study the benefits of flexibility

3

CR IP T

on fleet size.

Problem Definition

In this section we first give a formal definition of the FPVRP as described in the previous sections. For the sake of completeness, we recall afterwards the definition of the PVRP and of the IRP,

AN US

respectively. These definitions will be used in Section 4 to show the potential savings of the FPVRP with respect to them.

3.1

The Flexible Periodic Vehicle Routing Problem

Consider a complete and directed graph G = (N, A) with set of nodes N = {0} ∪ C where 0 denotes

M

the depot and C = {1, . . . , n} the set of customers. Let T = {1, . . . , H} be a discrete set of time

periods. Each customer i ∈ C has a total demand Wi over T and a storage capacity wi . The quantity qit delivered at each visit cannot be greater than wi and the sum of the quantities qit over T must

ED

be equal to Wi . A fleet of vehicles K = {1, . . . , m}, with a homogeneous capacity Q is available to

perform the service. A cost cij ≥ 0 is associated with each arc (i, j) ∈ A and is paid every time a

vehicle traverses the arc. The FPVRP is defined as the problem of finding the set of routes which

The Periodic Vehicle Routing Problem

CE

3.2

PT

minimizes the total routing costs satisfying customer demands.

In the PVRP, as it is defined in Christofides and Beasley (1984), a set of schedules, S, is given, where each schedule consists of a set of days in which customers receive service. Allocating a customer

AC

to a schedule implies that the customer will be visited, receiving the same amount of product wi , in every day of the schedule, i.e., Si = {s ∈ S :

P

t∈T

ast = fi } where Si is the schedule chosen for

customer i, fi is the visit frequency for customer i and   1 if day t ∈ T belongs to schedule s ∈ S, ast =  0 otherwise.

Three decisions have to be taken: select a schedule for each customer, assign customers to vehicles and generate the routes to be performed in each period of the time horizon. 5

ACCEPTED MANUSCRIPT

The PVRP is the problem of finding the set of routes that satisfy the periodic demand of customers, given a selected schedule for each customer, with a minimum total routing cost. A link between the FPVRP and the PVRP can be established by defining the storage capacity wi as the ratio between the total demand Wi and the expected frequency of visits fi as defined in the PVRP. Thus, wi =

Wi fi .

This way, in the FPVRP customer i has to be visited at least fi times,

the frequency defined in the PVRP. The customer may be visited more frequently if this leads to

3.3

CR IP T

cost savings.

The Inventory Routing Problem

The IRP is defined on the same graph G used to define the FPVRP. The only change with respect to the FPVRP setting is that customers are no more associated with a total demand Wi . Instead,

AN US

a demand dti is defined for each customer i ∈ C at each time period t ∈ T . Moreover, a starting

inventory level Ii0 is associated with each customer, together with a capacity wi . The distribution plan has to be such that each customer is able, at each time period, to satisfy the demand dti , thus he/she must have a sufficient quantity Iit in inventory. Moreover, the quantity delivered at each visit plus the inventory available when the visit is performed, should not exceed the capacity wi . As done in Archetti et al. (2014), we assume that, at each customer, the inventory level at time t ∈ T is the time t. No shortages are allowed.

M

inventory level at time t − 1 plus the amount delivered at time t minus the amount consumed at The aim of the IRP is to determine the quantity of product to deliver to each customer and the

ED

corresponding service routes of minimum total routing cost guaranteeing that there is no shortage at each customer in each time period.

Note that the IRP, as defined in Bertazzi and Speranza (2012, 2013) and Coelho et al. (2013),

PT

includes inventory holding costs in the objective function and inventory constraints at the supplier. In the following we do not consider any of these elements, in order to have a fair comparison with

CE

the FPVRP. Henceforth, such a version of the IRP will be referred to as the FPVRP with Inventory Constraints (FPVRP-IC).

FPVRP Savings: Worst-case Analysis

AC

4

In this section we illustrate the potential savings that may be achieved when comparing FPVRP

with PVRP and FPVRP–IC, respectively. In particular, we show theoretically the maximum savings that can be obtained in each case. Let z(P ) denote the optimal value to problem P .

Theorem 1. There exists no finite bound for the ratio

z(P V RP ) z(F P V RP ) .

Proof. Consider the following instance of the PVRP in which |T | = 2, |K| = Q and fi = 1 for

each customer i. There are three sets of customers. The first set is composed by Q customers with 6

ACCEPTED MANUSCRIPT

wi = Wi = Q. In the second set there are Q customers with wi = Wi = Q − 1 and the third set have Q customers with wi = Wi = 1. Moreover, all customers can be visited either in t = 1 or in

t = 2. Customers and depot locations are as follows. Each customer in the first set is co-located with a customer of the second set and they are spread around a circle centered at the depot with radius δ < 1 and a distance  apart. Customers of the third set are spread around a circle centered at the depot with a radius 1, perfectly aligned along the radius with the customers of the first and apart (Figure 1a). Note that, as the fleet is composed by Q vehicles

 δ

CR IP T

second set, with a distance

and |T | = 2, the maximum number of routes during T is 2Q. Moreover, as the total demand of customers is 2Q2 , all vehicles must be used and fully loaded in both periods.

The optimal and only solution of the PVRP is the one where Q routes are used to serve the customers in the first set with direct trips to the depot and the other Q routes serve one customer in the second set and one customer in the third set each. The cost of this solution is 2Qδ + 2Q

AN US

(Figure 1b).

The optimal solution of the FPVRP is the following (Figure 1c). In t = 1 one route is used to serve all customers in the third set. There remains Q − 1 routes in t = 1 and Q routes in t = 2.

They are constructed as follows. Without loss of generality, choose one customer of the second set as the first customer. Number all customers in a clockwise direction as i1 , ..., i|Q| if they belong to

the second set and j1 , ..., j|Q| if they belong to the first set. The first route delivers Q − 1 units

M

to customer i1 and 1 unit to customer j1 . The second route delivers Q − 1 units to customer j1 and 1 unit to customer i2 . The third route delivers delivers Q − 2 units to customer i2 and 2 units

ED

to customer j2 , and so on. We obtain Q + Q − 1 routes where the last route delivers Q units to

customer j|Q| . The odd routes are performed in t = 2 while the even routes are performed in t = 1

(as no customer can be served more than once in the same period). The cost of this solution is

PT

2 + (Q − 1)( δ ) + 2Qδ + (Q − 1)(2δ + ).

The ratio between z(P V RP ) and z(F P V RP ) is therefore

go to 0 this ratio tends to infinity.

AC

CE

goes to infinity and , δ and

 δ

2Qδ+2Q 2+(Q−1)( δ )+2Qδ+(Q−1)(2δ+) .

7

When Q

CR IP T

ACCEPTED MANUSCRIPT

(b) PVRP Solution

M

AN US

(a) Original Graph

t=2

ED

t=1

(c) FPVRP Solution

PT

Figure 1: Worst–case Analysis 

CE

The proof presented above is quite similar to the one proposed by Gueguen (1999) for the analysis of the benefits of the SDVRP with respect to the VRP.

AC

We now compare the FPVRP with the FPVRP-IC as defined in Section 3.3. The following result

holds.

Theorem 2. There exists no finite bound for the ratio

z(F P V RP −IC) z(F P V RP ) .

Proof. Let us consider the instance introduced in the proof of Theorem 1. In the FPVRP-IC, initial inventory levels at customers must be defined. We define them as follows: the initial inventory level is equal to 0 for the customers in the first set, to Q − 1 for the customers in the second set and to 1 for the customers in the third set. The customers demands are equal to Q − 1 and 1, at time t = 1 and t = 2, for customers of the second and third set, respectively. For the customers in 8

ACCEPTED MANUSCRIPT

the first set, demand at time t = 1 is equal to Q and demand at t = 2 is equal to 0. Note that the customers of the second and third set cannot be served at time t = 1 as their initial inventory level is equal to the storage capacity. On the other side, customers of the first set have to be served at time t = 1 as their initial inventory level is 0 and the demand at t = 1 is positive. Thus, the only feasible solution for the FPVRP-IC corresponds to the solution of the PVRP shown in the proof of Theorem 1, i.e., Q routes are used to serve the customers in the first set with direct trips to the

CR IP T

depot at time t = 1 and Q routes serve one customer in the second set and one customer in the third set each at time t = 2. The cost of this solution is 2Qδ + 2Q. The solution of the FPVRP does not change. Thus the ratio is the same and tends to infinity.

Mathematical Formulations

AN US

5



In this section we propose a load-based MILP formulation for the FPVRP. We also present the load-based MILP formulations for the FPVRP-IC and the PVRP, which have been used in the computational experiments.

Traditional VRP formulations with multiple vehicles use decision variables with a vehicle index that shows which vehicle traverses each arc. This results in a high number of decision variables,

M

particularly in problems where decisions must be made at the different periods of a given time horizon, like the ones we study in this paper. In order to mitigate this difficulty, recent works on

ED

several VRP variants have proposed the use of the so-called load-based formulations, in which decision variables identify the arcs used in the solutions without making explicit the vehicles that traverse them (Letchford and Salazar-González, 2015; Archetti et al., 2014). For this, an additional

PT

set of continuous commodity flow variables is needed to guarantee that routes are properly defined. Such formulations tend to be quite effective in practice, although their linear programming (LP)

CE

relaxations are usually weaker than their traditional counterpart. On the one hand they have a smaller number of variables. On the other hand their implementation does not require the use of

AC

sophisticated techniques, like branch-and-cut or column generation.

5.1

MILP Formulation of the FPVRP

For the load-based formulation for the FPVRP we introduce two sets of binary decision variables,

which identify the arcs that are traversed and the customers that are visited at each time period, respectively. We also define a set of integer decision variables that indicate the number of vehicles that are used at each time period. This can be easily computed by counting the number of arcs leaving the depot at each time period. We finally define two additional sets of continuous variables. The first one indicates the load of the vehicles when they traverse the arcs. The second one, gives the amount of product that is delivered to each visited customer. The definition of the decision 9

ACCEPTED MANUSCRIPT

variables is the following:

  1 t = • yij  0

If customer i is visited at time period t, otherwise. If arc (i, j) is traversed at time period t, otherwise.

• z0t : Number of vehicles used at time period t ∈ T .

(i ∈ C, t ∈ T ) ((i, j) ∈ A, t ∈ T )

CR IP T

  1 • zit =  0

t : Load of the vehicle traversing arc (i, j) ∈ A, at time period t ∈ T . • lij

• qit : Quantity delivered to customer i ∈ C at time period t ∈ T .

X X

t∈T (i,j)∈A

s.t.

t cij yij

qit ≤ wi zit X

i∈C

qit ≤ Qz0t

t yij = zit

ED

X

j|(i,j)∈A

X

t yij =

PT

j|(i,j)∈A

X

AC

CE

j|(i,j)∈A

i ∈ C, t ∈ T

M

min

AN US

The formulation of the FPVRP is as follows.

X

t yji

j|(j,i)∈A

t lij −

X

j|(j,i)∈A

t t lij ≤ Qyij

X

(2)

t∈T

(3)

i ∈ C, t ∈ T

(4)

i ∈ N, t ∈ T

(5)

   −qit , i ∈ C t P t lji = i ∈ N, t ∈ T  qi0 , i = 0  i0 ∈C

(1)

(i, j) ∈ A, t ∈ T

(6) (7)

t∈T

(8)

t y0i

t∈T

(9)

qit = Wi

i∈C

(10)

t∈T

qit ≥ 0

i ∈ C, t ∈ T

(11)

j|(0,j)∈A

z0t = X

z0t zit

t yij

t y0j ≤ |K|

X

i∈C

∈Z

∈ {0, 1} ∈

t {0, 1}, lij

t∈T i ∈ C, t ∈ T

≥0

(i, j) ∈ A, t ∈ T

10

(12) (13) (14)

ACCEPTED MANUSCRIPT

Objective function (1) minimizes the routing costs. Constraints (2) impose that the quantity delivered to each customer does not exceed wi . Constraints (3) establish that the total quantity delivered at time t does not exceed the total capacity of the vehicles used at time t. Constraints (4) state that one arc has to exit from a node of a visited customer. Constraints (5) are the flow conservation constraints. Constraints (6) are the load conservation constraints. Constraints (7) impose that the vehicle load does not exceed its capacity and link y and l variables. Constraints

CR IP T

(8) ensure that the number of vehicles used is at most |K|. Constraints (9) ensure that variables z0t take the appropriate value. Constraints (10) impose that the total quantity delivered at each customer at the end of the time horizon is equal to Wi . Finally, Constraints (11)–(14) define the domain of the variables.

The above formulation has |T | × (|C| + |A|) binary variables, |T | general integer variables, and

5.2

AN US

|T | × (|C| + |A|) continuous variables. The number of constraints is |T | × (6|C| + 2|A| + 6) + |C|.

Mathematical Formulation of the FPVRP-IC

Below we present a load-based formulation for the FPVRP-IC, largely based on the formulation proposed by Archetti et al. (2014). In addition to the l, q, y, and z variables used in the FPVRP it uses another set of continuous variables that represent the customers inventory level at each time

M

period. As done in Archetti et al. (2014), the calculation of the inventory levels occurs after delivery of qit and consumption of dti are carried out, and these levels cannot exceed wi − qit .

ED

Decision Variables:

PT

• Iit : Inventory level at customer i ∈ C at the end of time period t ∈ T .

AC

CE

The FPVRP-IC formulation is the following:

min

X X

t∈T (i,j)∈A

s.t.

(15)

cij yijt

Iit = Iit−1 − dti + qit

qit ≤ wi − Iit−1

Equations (2) – (9)

i ∈ C, t ∈ T i ∈ C, t ∈ T

(16) (17)

Equations (11) – (14) (18) The objective function (15) is the minimization of the total routing cost. Constraints (16) and (??) determine the inventory level over time and avoid stock–out situations. Note that, 11

ACCEPTED MANUSCRIPT

since T = {1, . . . , H}, when t = 1, Iit−1 = Ii0 corresponds to the initial inventory at customer

i. Constraints (17) ensure that delivered quantity does not exceed the maximum quantity need at each customer and at each time period. The remaining constraints have the same meaning as in the FPVRP model. Mathematical Formulation of the PVRP

CR IP T

5.3

Below we present the formulation for the PVRP we have used in our experiments. To the best of our knowledge this is the first load-based formulation in the literature for the classical PVRP. It uses the same y, z and l variables as the FPVRP formulation. In addition, the following set of binary decision variables is defined to determine the schedule that is chosen for each customer:  

AN US

1 If customer i ∈ C is visited according to schedule s ∈ Si , • vis =  0 otherwise. Then, the load-based formulation for the PVRP is:

t∈T (i,j)∈A

s.t.

X

vis = 1

s∈Si

zit = yijt

cij yijt

M

X X

X

ED

min

vis ast

s∈Si zit +

2

X

t lij

PT

≤

zjt

CE

j|(i,j)∈A

−

X

t lji

j|(j,i)∈A

AC

Equations (3) – (5)

=

    

−wi zit , P

j∈C

wj zjt ,

i∈C

i=0

(19)

i∈C

(20)

t ∈ T, i ∈ C

(21)

t ∈ T ; i, j ∈ C(i 6= j)

(22)

i ∈ N, t ∈ T

(23)

i ∈ C, s ∈ Si

(24)

Equations (7) – (8) Equations (13) – (14) vis ∈ {0, 1}

The objective function (19) minimizes the total routing costs. Constraints (20) ensure that a schedule is assigned to each customer. Constraints (21) relate the selected schedules with customer visits. Constraints (22) allow to use edges connecting two customers only if both customers are assigned to the same period. Constraints (23) are flow conservation constraints. Finally, Constraints (24) define the domain of the new set of variables vis . 12

ACCEPTED MANUSCRIPT

6

Inequalities

In order to strengthen the formulation and improve the computing times, several families of inequalities have been tested for the FPVRP formulation presented above. Note that classical inequalities that can be used when the amount of product delivered to each customer at each time period is known, like the ones described in Letchford and Salazar-González (2015),

CR IP T

cannot be applied in our case because such an amount is a decision variable. The inequalities that have been considered are the following:

Inequalities I1. Sum of final loads: All vehicles return to the depot with an empty load. Indeed, these inequalities are not valid, as there are feasible FPVRP solutions that do not satisfy them. Instead, they are optimality cuts, since there is at least one optimal solution

AN US

that satisfies them. Therefore, they can be used to reduce the domain of the solutions that are explored.

XX

t lj0 =0

(25)

t∈T j∈C

Inequalities I2. Symmetry-Breaking of Routes: The following inequalities partially

M

break the symmetry of the routes by exploiting the fact that arc costs are symmetric. According to them only the routes with a certain orientation will be considered (as the

ED

same route in the opposite orientation will have the same cost). Among the two possible orientations for a route we chose the one that starts with the customer with the lowest index.

PT

For this, we impose

CE

t yi0 ≤

X

t y0r ,

∀i ∈ C, ∀t ∈ T

r≤i

(26)

That is, if arc (i, 0) enters the depot at time period t, then there must be a lowest index

AC

arc (0, r) with r ≤ i which exits the depot at this time period. The symmetry-breaking inequalities are also optimality cuts, but not valid inequalities. Valid Inequalities I3: The relation between variables y and z can be imposed in several

ways. In our formulation this is done in a two-step fashion. On the one hand, the flow balance constraints (6) relate the load variables l to the q’s, which, in turn, activate the z’s. On the other hand, constraints (7) relate the load variables l to the arc variables y. Nevertheless, this relation can be stated in a more direct way similarly to the constraints imposed in

13

ACCEPTED MANUSCRIPT

Christofides and Beasley (1984) for the PVRP. In particular, no arc yij can be used at time period t unless customers i and j are visited at the same time period. Therefore, the following inequalities are valid: yijt

zit + zjt ≤ , 2

∀i, j ∈ C, ∀t ∈ T

(27)

period the above inequalities can be reinforced to: t yijt + yji ≤

∀i < j ∈ C, ∀t ∈ T

Computational Experiments

(28)

AN US

7

zit + zjt , 2

CR IP T

Taking into account that no arc will be traversed in both directions in the same time

In this section we present the numerical results of the computational experiments we have conducted. The aim of these experiments was twofold: on the one hand to analyze the computational difficulty of the FPVRP and the effectiveness of the proposed load-based formulation. On the other hand, to highlight the benefits derived from allowing flexibility in the PVRP and the IRP, by comparing the solutions produced by the three models. For the

M

computational experiments, all three formulations were implemented in C++ with ILOG Concert Technology API and CPLEX 12.5.0.0, running on a HP Intel(R)-Xeon(R) 2.4GHz

ED

Workstation with 32GB RAM (Win Server 2012, 64 bits). Default parameters were used. All computing times were limited to 14400 seconds. Before presenting the numerical results we describe the sets of benchmark instances we

PT

have used. Data instances and complementary tests can be found at http://or-brescia.

7.1

CE

unibs.it/instances.

Benchmark instances

AC

Three different sets of benchmark instances have been used: • Set 1 (S1): The IRP instances proposed by Archetti et al. (2014). This set consists of 40 benchmark instances with |C| ∈ {5, 10, 15, 20} and a 3–period time horizon, i.e. |T | = 3. Note that larger instances were considered in Archetti et al. (2014) but we do

not consider them in our experiments as they turned out to be too hard to solve for the FPVRP. • Set 2 (S2): The PVRP instances proposed by Francis et al. (2006). This set consists of 24 instances with |C| ∈ {7, 9, 11, 15, 49} and time horizon |T | = 5. Instances were 14

ACCEPTED MANUSCRIPT

generated in a similar way to those used by Francis et al. (2006) and according to the database provided by one of the authors. Schedules for each customer are assigned according to their frequency of visits as explained in the following. • Set 3 (S3): A newly generated set of 35 PVRP instances with clustered customers. The time horizon used in the whole set is |T | = 5. Other parameters are: the number

CR IP T

of customers |C|, the vehicle capacity Q, the number of clusters p, a radius r which determines the coverage area of each cluster, and a parameter β which together with r

determines the minimum distance β × r among the centers of the clusters. Five instances

were generated for each combination of |C| ∈ {10, 15, 20} and r ∈ {0.15, 0.30}, plus five

more instances with |C| = 20 and r = 0.50. Vehicles capacities, Q, have been set to 200,

250 and 300 for 10, 15 and 20 customers, respectively. For |C| = 10, we set the number

AN US

of clusters to p = 2, whereas for instances with |C| ∈ {15, 20} we set p = 3. When

r ∈ {0.15, 0.30}, the value of β has been set to 2, which avoids clusters overlapping. Instead, for r = 0.50 we have set β = 1, allowing that a customer can belong to more than one cluster.

Instances have been generated as follows:

M

– The depot is located at the center of the unit square.

– p centers (one for each cluster) are randomly generated from a uniform distribution

ED

in the unit square. Each generated center must satisfy the minimum distance condition β × r with respect to the others.

PT

– Once all the centers have been fixed, the customers are generated in such a way that clusters are balanced, i.e. each cluster contains up to

l

|C| p

m

customers. Clusters

are progressively filled: first, one customer is generated for each cluster; then, a

CE

second one; and so on, until |C| customers have been randomly generated using an uniform distribution in the circles around the cluster’s centers of radius r.

AC

– A demand wi , i ∈ C, is randomly generated from an integer uniform distribution in [1, Q].

– A number of visits is randomly associated to each customer according to the following options: fi = 5, 3 or 2.

– A set of schedules Si is assigned to each customer according to its frequency fi . Figure 2 shows an example of an instance with |C| = 10, Q = 10 and p = 4. As can be

seen, two clusters have three customers each and two other clusters have two customers

15

ACCEPTED MANUSCRIPT

each. Nodes with the same color represent customers with the same assigned schedule. The number inside each node represents the assigned demand wi . ൒ ߚ‫ݎ‬

2

‫ݎ‬

5

9

7 ‫ݎ‬

3 0.5

CR IP T

1.0

‫ݕ‬

3 6

‫ݎ‬

AN US

1 ‫ݎ‬

4

4

0

0.5

1.0

‫ݔ‬

M

Figure 2: Clustered instance with |C| = 10, Q = 10, p = 4, and |S| = 3.

• Schedules for S2 and S3: Different schedule options are considered for the tests.

ED

– fi = 2 : Si = {(0, 0, 1, 0, 1), (0, 1, 0, 0, 1), (0, 1, 0, 1, 0), (1, 0, 0, 1, 0), (1, 0, 1, 0, 0)}, – fi = 3 : Si = {(0, 1, 0, 1, 1), (0, 1, 1, 0, 1), (1, 0, 1, 0, 1), (1, 0, 1, 1, 0), (1, 1, 0, 1, 0)},

7.2

PT

– fi = 5 : Si = {(1, 1, 1, 1, 1)}.

Evaluation of Alternative Formulations for the FPVRP

CE

Our first set of experiments focused on the evaluation of the different formulations for the FPVRP and their effectiveness. Preliminary tests showed that load-based formulations for the

AC

FPVRP greatly outperform vehicle index formulations. Thus, this second type was excluded from any further consideration in our experiments. Initially, we tested several variations of the formulation presented in Section 5.1 and the

obtained results indicated that the best performance was attained with formulation (1)-(14).

Then, we compared this base formulation against alternative reinforcements resulting from the addition of different combinations of inequalities presented in Section 6. For this test, we selected a subset of benchmark instances whose optimality was particularly difficult to prove with the base formulation. This subset consists of 16 S2 instances (with |K| ∈ {3, 4}) plus 16

ACCEPTED MANUSCRIPT

the 35 new S3 instances. For these experiments we considered the Case A of the schedule assignment mentioned in Section 7.1. Average results of the percentage optimality gaps at termination of the computation (over all instances in each group) are summarized in Table 2. These optimality gaps have been computed as 100 × BestSol−LB , where BestSol is the best-known feasible solution value for each BestSol instance (among all tested versions), and LB is the lower bound at termination of the tested

CR IP T

version. Each row corresponds to a group of instances. The first column indicates the set to which the group of instances belongs. The next three columns give some characteristics of the instances of each group: r (for instances in S3 ), |C| and Q (an average value is reported for S2). The number of instances in each group is given in column Inst. Average results of the

combinations tested are given in the following six columns: entries in column Base correspond

AN US

to the base formulation (1)-(14), whereas under I1, I2, I3, we show the results of the base formulation reinforced respectively with (i) the empty load return inequalities (25); (ii) the symmetry breaking inequalities (26); (iii) the valid inequalities (28). The last two columns correspond to the combinations of (i) with (ii) and of (i) with (iii). In each case, the number of instances optimally solved, within the allowed computing time, is given in parenthesis next to the percentage of optimality gap. These values are only given for the groups of instances

M

of S2 because it was never possible to prove the optimality for any of the instances in S3. Other combinations of inequalities were also tested but did not give significantly different r

S2 S3

0.15

15 15 10 15 20 10 15 20 20

Q 999.38 200 250 300 200 250 300 300

CE

0.30

|C|

AC

0.50

Inst.

Base

I1

I2

I3

I1+I2

I1+ I3

8 8 5 5 5 5 5 5 5

0.43 (3) 0.32 (3) 2.58 3.83 5.71 1.51 1.97 2.09 2.59

0.30 (3) 0.29 (4) 2.78 3.82 5.71 1.56 1.97 2.15 2.50

0.38 (5) 0.31 (5) 2.65 3.76 5.67 1.54 1.94 2.16 2.55

0.45 (3) 0.52 (3) 3.16 4.15 5.68 1.92 2.15 2.35 2.97

0.26 (5) 0.26 (5) 2.49 3.55 5.65 1.51 1.89 2.15 2.49

0.36 (3) 0.35 (3) 2.93 3.90 5.57 1.71 2.05 2.26 2.85

PT

Set

ED

results.

Table 2: Summary of results of FPVRP reinforced with different combinations of inequalities.

Results show that strengthened formulations reduced the percentage gap in nearly all cases. Broadly speaking, both the empty load return inequalities (25) and the symmetry breaking inequalities (26) were effective. The effectiveness of inequalities (28) that relate y and z variables is not so clear and in some cases the base formulation alone gives better results than if they are used. We attribute this behavior to the high number of inequalities of 17

ACCEPTED MANUSCRIPT

this family (O(|C|2 |T |)). As can be seen, the best results are those of column I1+I2, which

corresponds to the base formulation reinforced with (25) and (26). Therefore, all subsequent tests used the corresponding strengthened FPVRP formulation. The results obtained with the strengthened FPVRP formulation for the full set of benchmark instances are summarized in Table 3. Again, rows correspond to groups of instances with similar characteristics. The first five columns describe the characteristics of the instances: r

CR IP T

(for instances in S3 ), |C|, |K|, Q, and Dt (the total demand that must be distributed at each

time period). When not all the instances in the group have the same parameters, minimum and maximum values are displayed. In the following columns, Inst. gives the number of instances in each group, and O/F gives the number of instances of the group that terminated

with a solution with proven optimality (O) and with only proven feasibility (F). The average

AN US

percentage optimality gaps at termination are given in column %Gap. Finally, the last two columns refer to the computing times: avrg. for the average times and range for the minimum and maximum of such values. r

|K|

|C| 5 10 15 20

S2

7 9 11 15 49

3 3 3 3 4

10 15 20 10 15 20 20

4-8 6 - 10 10 - 12 5-8 6-9 10 - 13 7 - 14

0.15 0.30

3 3 3 3

%Gap

79 - 175 272 - 480 340 - 619 512 - 867

193 - 304 458 - 640 681 - 845 999- 1156

10 10 10 10

10/0 10/0 9/1 6/4

0.00 0.01 0.34 2.32

1.40 202.90 3086.80 7214.50

0-7 7 - 942 101 - 14400 147 - 14400

496 - 871 1037 - 1208 546 - 1521 757 - 1240 1111

699 1206 851 1360 4443

4 4 8 8 1

4/0 4/0 7/1 5/3 0/1

0.00 0.00 0.26 0.24 49.76

0.25 27.25 1965.50 6313.88 14400

0-1 9 - 57 22 - 14400 17 - 14400 14400

200 250 300 200 250 300 300

664 - 1401 1495 - 2484 2813 - 3408 820 - 1436 1488 - 2225 2824 - 3691 2059 - 4076

5 5 5 5 5 5 5

0/5 0/5 0/5 0/5 0/5 0/5 0/5

2.49 3.55 5.65 1.51 1.89 2.15 2.49

14400 14400 14400 14400 14400 14400 14400

14400 14400 14400 14400 14400 14400 14400

avrg.

range

CE

0.50

-

O/F

PT

S3

2 2 2 2

Inst.

ED

S1

Time

Dt

Q

M

Set

Table 3: Summary of results of selected FPVRP formulation for the complete set of instances.

AC

The obtained results highlight the difficulty of the FPVRP. This will become more evident

when the results of the FPVRP are compared with the results of the FPVRP-IC and the PVRP. Still, it was possible to optimally solve 35 out of the 40 S1 instances and 20 out of the 24 S2 instances. For the set S1, the percentage optimality gap at termination was always below 5%, except for one 20 customer instance (abs1n20_2). Optimality gaps of unsolved S2 instances were always below 1%, except for a 11 customer instance (Instn12t5k3_1), with a 2.11% gap, and the largest 49 customer instance for which the relative percentage deviation between the upper and lower bounds at termination was of nearly 25%. The results of the individual instances in S1 and S2 (see Tables 4, 5, and 6 below) show that, when all other 18

ACCEPTED MANUSCRIPT

parameters are similar, a clear indicator of the difficulty of an instance is the fleet size. As can be seen, the new instances of set S3 are considerably harder to solve than those of sets S1 and S2 of similar sizes. We believe this is due to two factors. The first is the fleet size which is much larger in S3 than in S1 and S2. The second is the proximity of the clustered customers (particularly for the small radius 0.15), which makes it particularly difficult to discriminate among solutions that permute the order in which neighboring customers are

7.3

CR IP T

visited.

Comparison of the FPVRP with other VRPs with periodic demand

The final series of experiments were oriented to analyze the potential advantage of the FPVRP relative to the FPVRP-IC and the PVRP. Since all three models focus on the overall routing

AN US

costs throughout the time horizon, potential advantages can be quantified in terms of the percentage relative reduction in the objective function value of the compared models. In particular, we will use the value 



ZMod − ZFPVRP %Imp = max 0, × 100, ZMod

M

where ZFPVRP and ZMod are the best-known values for the FPVRP and the compared model, respectively.

Comparison between the FPVRP and the FPVRP-IC

ED

7.3.1

For the comparison between the FPVRP and the FPVRP-IC we used the set of benchmark

PT

instances S1. Each instance was run with both the strengthened FPVRP formulation and the FPVRP-IC formulation with a time limit of 14400 seconds. The results are shown in Table 4 where the first five columns show the name of the instances and their characteristics. The next

CE

two blocks of four columns each correspond to the FPVRP and to the FPVRP-IC results, respectively. Column Status, indicates whether the instance was solved to prove optimality

AC

(O), or a feasible solution was found but its optimality was not proven (F); BestSol gives the value of the best solution obtained in the run, BestLB is for the lower bound at termination, and Time is the computing time consumed (in seconds). The last column of the table, %Imp, gives the percentage relative improvement obtained with the FPVRP with respect to the FPVRP-IC. The results of Table 4 show that, except for the small instances with only five customers, in general, the computing times required by the FPVRP are substantially greater than those of the FPVRP-IC. This becomes more evident as the size of the instances increases and is reflected by the fact that five instances were not optimally solved for the FPVRP within the 19

ACCEPTED MANUSCRIPT

maximum computing time, whereas this number reduces to three for the FPVRP-IC. The entries in column %Impr show that there are instances where both models have the same optimal value, although in most cases the FPVRP produces an improvement (up to 9.48%) with respect to the FPVRP-IC which increases with the size of the instances. It is worth noting that, even for the two instances that were not optimally solved with the FPVRP but were optimally solved with the FPVRP-IC, the best FPVRP solution at termination

CR IP T

was considerably better than the optimal FPVRP-IC solution. These results show that, if possible, it is worthwhile to increase flexibility and not consider inventory levels, as done in the FPVRP, as this may lead to remarkable savings. FPVRP

FPVRP-IC

Dt

Status

BestSol

BestLB

Time

Status

2 3 2 3 2 3 2 3 2 3

144 96 118 79 228 152 134 89 175 117

193 193 158 158 304 304 179 179 234 234

O O O O O O O O O O

1301.85 1335.88 1088.72 1494.37 2109.51 2864.95 1504.27 2224.13 1091.97 1386.18

1301.85 1335.86 1088.72 1494.37 2109.51 2864.87 1504.27 2223.94 1091.97 1386.18

0 0 0 2 1 1 7 3 0 0

O O O O O O O O O O

10

2 3 2 3 2 3 2 3 2 3

476 317 408 272 343 229 411 274 480 320

635 635 545 545 458 458 548 548 640 640

O O O O O O O O O O

1936.15 2369.40 2491.71 3194.02 1980.71 2372.91 2115.97 2756.47 1746.02 2014.42

1936.01 2369.16 2491.52 3193.71 1980.71 2372.73 2115.78 2756.20 1746.02 2014.42

96 737 7 185 17 13 9 942 11 12

abs1n15_1 abs1n15_2 abs2n15_1 abs2n15_2 abs3n15_1 abs3n15_2 abs4n15_1 abs4n15_2 abs5n15_1 abs5n15_2

15

2 3 2 3 2 3 2 3 2 3

619 413 592 395 633 422 538 359 510 340

826 826 790 790 845 845 718 718 681 681

O O O O O O O O O F

1915.91 2349.28 2161.09 2388.97 2373.10 2646.11 2064.15 2403.11 2192.45 2678.85

1915.77 2349.05 2160.89 2388.73 2372.90 2645.86 2063.95 2402.87 2192.23 2634.11

abs1n20_1 abs1n20_2 abs2n20_1 abs2n20_2 abs3n20_1 abs3n20_2 abs4n20_1 abs4n20_2 abs5n20_1 abs5n20_2

20

799 533 782 521 768 512 749 499 867 578

1066 1066 1043 1043 1024 1024 999 999 1156 1156

O F O O O O F F O F

2345.27 3004.62 2148.82 2365.64 2283.53 2529.42 2888.29 3408.32 2854.24 3562.11

2345.04 2741.40 2148.63 2365.40 2283.31 2529.19 2858.71 3287.19 2853.96 3403.85

AC

CE

5

2 3 2 3 2 3 2 2 3 3

M

abs1n10_1 abs1n10_2 abs2n10_1 abs2n10_2 abs3n10_1 abs3n10_2 abs4n10_1 abs4n10_2 abs5n10_1 abs5n10_2

abs1n5_1 abs1n5_2 abs2n5_1 abs2n5_2 abs3n5_1 abs3n5_2 abs4n5_1 abs4n5_2 abs5n5_1 abs5n5_2

PT

|K|

BestLB

Time

1301.85 1335.88 1088.72 1494.37 2302.82 2864.95 1650.73 2224.13 1091.97 1386.18

1301.85 1335.88 1088.72 1494.23 2302.82 2864.95 1650.59 2224.13 1091.97 1386.18

0 0 0 4 2 0 0 1 0 4

0.00 0.00 0.00 0.00 8.39 0.00 8.87 0.00 0.00 0.00

O O O O O O O O O O

1960.99 2429.55 2554.79 3214.05 1980.71 2410.50 2240.93 2943.14 1848.20 2151.45

1960.82 2429.55 2554.79 3213.88 1980.71 2410.50 2240.73 2942.87 1848.20 2151.45

18 31 14 19 6 39 35 279 14 19

1.27 2.48 2.47 0.62 0.00 1.56 5.58 6.34 5.53 6.37

175 5437 898 719 101 116 455 7297 1271 14399

O O O O O O O O O O

1915.91 2402.36 2185.68 2388.97 2373.10 2646.11 2199.78 2572.55 2309.75 2959.31

1915.89 2402.14 2185.46 2388.97 2373.10 2646.11 2199.57 2572.30 2309.53 2959.02

21 295 194 39 11 20 188 705 88 5829

0.00 2.21 1.13 0.00 0.00 0.00 6.17 6.59 5.08 9.48

6768 14399 156 4511 615 147 14399 14399 2352 14399

O F O O O O O F O F

2410.91 3103.40 2148.82 2393.13 2283.53 2529.42 3136.22 3664.52 2859.60 3567.47

2410.70 2867.75 2148.82 2392.90 2283.53 2529.28 3135.91 3602.53 2859.35 3507.53

6343 14399 18 658 23 16 2782 14399 379 14399

2.72 3.18 0.00 1.15 0.00 0.00 7.91 6.99 0.19 0.15

Average Improvement

Table 4: Comparison between FPVRP and FPVRP-IC using S1 instances.

20

%Imp

BestSol

AN US

Q

|C|

ED

Instance

2.56

ACCEPTED MANUSCRIPT

7.3.2

Comparison between the FPVRP and the PVRP

For the comparison between the FPVRP and the PVRP we have used instances of both S2 and S3. Each instance in these two sets was run with the strengthened FPVRP formulation (1)-(14) and with the PVRP formulation (19)-(24). Again, the time limit for each run was 14400 seconds. The results for the instances of set S2 are presented in Table 5, where instances

CR IP T

with the same number of customers are ordered by decreasing capacity. The columns labeled as %Gap show the relative CPLEX MIP gap after the termination of each optimization. Most of the instances with |C| = 7, |C| = 9 and |C| = 11 were solved to optimality with

both the FPVRP and the PVRP formulations. Once more we can observe that, in terms of the computing times, the FPVRP formulation is more demanding than that for the PVRP. This could be expected as the FPVRP incorporates additional decisions to the ones of the

AN US

PVRP. While the computing times to optimally solve the PVRP instances remain negligible, they become significant for the FPVRP as the size of the instances increases. In particular, for five of the instances the optimality of the best solution found could not be proven within the allowed computing time. Concerning the improvement of the FPVRP with respect to PVRP, we notice that no improvement can be perceived related to the flexibility in the delivered quantities since the total vehicles capacities are enough to cover the demand of all customers

M

in the same day. Nevertheless, the FPVRP still shows slight improvements related to the non-dependency of schedules. Moreover, although the FPVRP produced no improvement in

ED

the large 49 customers instance, it was able to find a feasible solution within the time limit.

AC

CE

PT

On the contrary, the PVRP produced no solution for that instance.

21

ACCEPTED MANUSCRIPT

%Imp

Q

Status

BestSol

BestLB

%Gap

Time

Status

%Gap

Time

3

699

871 765 705 496

O O O O

5.82 5.82 5.82 6.30

5.82 5.82 5.82 6.30

0.00 0.00 0.00 0.00

0 1 0 0

O O O O

0.00 0.00 0.00 0.00

0 0 0 0

0.00 0.00 0.00 0.00

9

3

1206

1208 1058 1045 1037

O O O O

6.76 6.86 6.90 6.90

6.76 6.86 6.90 6.90

0.00 0.01 0.01 0.01

9 17 26 57

O O O O

0.00 0.01 0.00 0.00

6 7 6 7

0.05 0.57 0.00 0.00

Instn12t5k3_6 Instn12t5k3_5 Instn12t5k3_7 Instn12t5k3_8 Instn12t5k3_4 Instn12t5k3_2 Instn12t5k3_3 Instn12t5k3_1

11

3

851

1521 1491 1399 1146 925 802 748 546

O O O O O O O F

4.47 4.47 4.47 4.47 4.47 4.50 4.54 4.83

4.47 4.47 4.47 4.47 4.47 4.50 4.54 4.73

0.01 0.01 0.01 0.01 0.01 0.01 0.01 2.01

22 31 37 26 44 102 1062 14400

O O O O O O O O

0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01

6 7 9 9 8 30 29 1404

0.00 0.00 0.00 0.00 0.00 0.67 0.47 0.59

Instn16t5k3_3 Instn16t5k3_4 Instn16t5k3_1 Instn16t5k3_2 Instn16t5k3_8 Instn16t5k3_7 Instn16t5k3_5 Instn16t5k3_6

15

3

1360

1240 1232 1056 1030 1027 851 802 757

O O O O O F F F

5.62 5.62 5.62 5.63 5.63 5.74 5.76 5.78

5.62 5.62 5.62 5.63 5.63 5.70 5.72 5.75

0.01 0.00 0.00 0.01 0.01 0.73 0.76 0.45

52 17 31 2100 5113 14400 14400 14399

O O O O O O O O

0.00 0.00 0.00 0.01 0.01 0.01 0.00 0.00

22 26 33 1743 1356 8054 37 26

0.00 0.00 0.00 0.40 0.40 0.38 0.07 0.00

Instn50t5k4

49

4

4443

1111

F

18.69

12.48

33.25

14400 I — — Average Improvement

|K|

Instn10t5k3_3 Instn10t5k3_1 Instn10t5k3_2 Instn10t5k3_4

Instn8t5k3_1 Instn8t5k3_3 Instn8t5k3_2 Instn8t5k3_4

7

AN US

Dt

|C|

M

Instance

PVRP

CR IP T

FPVRP

0.00 0.14

ED

Table 5: Comparison between FPVRP and PVRP formulations using S2 instances.

On the other hand, Table 6 shows the results of the comparison between the FPVRP

PT

and the PVRP with the S3 instances. For the PVRP instances, computing times increased significantly in some of them. The exception was instances FPVRP_n10k5t5_1, for which no feasible PVRP solution could be found in the 14400 seconds allowed. The strengthened

CE

FPVRP formulation was not able to prove the optimality of the best solution found for any of the instances in the set, although it was always able to find a feasible solution. In

AC

general, the percentage optimality gaps at termination are small, with higher values (up to 5.51%) for the instances with the smallest radius r = 0.15. Despite not knowing whether or not the best-known FPVRP solutions are optimal, in all cases they produce substantial improvements in the routing costs with respect to the optimal PVRP solutions. In these instances, the benefit of using FPVRP is more noticeable for instances of S2 because of their structure. In fact, the improvement is significant since it is allowed to partition customer demands into several periods. The range of such improvements is 3.03%- 12.27%, with an average improvement of 7.43% in comparison to the PVRP.

22

ACCEPTED MANUSCRIPT

FPVRP 0.15

FPVRP_n10k6t5_1 FPVRP_n10k6t5_2 FPVRP_n10k5t5_3 FPVRP_n10k5t5_4 FPVRP_n10k8t5_5 FPVRP_n15k9t5_1 FPVRP_n15k9t5_2 FPVRP_n15k7t5_3 FPVRP_n15k7t5_4 FPVRP_n15k6t5_5 FPVRP_n20k10t5_1 FPVRP_n20k12t5_2 FPVRP_n20k10t5_3 FPVRP_n20k13t5_4 FPVRP_n20k12t5_5

0.30

FPVRP_n20k14t5_1 FPVRP_n20k10t5_2 FPVRP_n20k7t5_3 FPVRP_n20k10t5_4 FPVRP_n20k11t5_5

0.50

10

15

20

10

15

20

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

14 10 7 10 11

PVRP

%Imp

Q

Dt

Status

BestSol

BestLB

% Gap

Time

Status

% Gap

Time

200

938 664 874 738 1401 2317 1495 2484 1767 1606 2813 3408 3138 2957 2874

F F F F F F F F F F F F F F F

20.78 12.44 13.23 13.53 25.91 34.26 17.42 25.25 32.11 23.89 24.58 36.10 23.70 35.35 29.45

20.30 12.12 12.87 13.05 25.58 33.39 16.51 24.45 31.25 22.84 23.58 35.45 22.79 34.53 28.89

2.36 2.64 2.8 3.68 1.29 2.61 5.51 3.27 2.75 4.60 4.24 1.83 3.99 2.37 1.94

14400 14400 14400 14400 14400 14403 14400 14400 14400 14400 14402 14400 14404 14401 14399

I O O O O O F O O F F F F F F

— 0.01 0.00 0.01 0.00 0.01 1.69 0.01 0.01 1.54 2.33 0.48 2.47 4.81 4.31

14399 145 11 35 1 2445 14399 1213 4587 14399 14399 14399 14399 14399 14399

— 4.37 4.10 10.81 6.60 5.55 12.27 10.12 5.81 7.28 5.25 7.71 11.03 10.34 6.36

1028 1176 922 820 1436 2116 2225 1692 1530 1488 2824 3537 2849 3691 3308

F F F F F F F F F F F F F F F

19.04 13.89 14.50 14.38 19.40 27.16 29.72 27.84 17.43 20.59 29.57 31.55 26.07 42.61 34.05

18.97 13.59 14.16 14.04 19.34 26.64 29.19 27.38 17.06 20.18 28.97 30.78 25.18 41.62 33.76

0.37 2.21 2.4 2.42 0.31 1.95 1.82 1.68 2.17 2.03 2.07 2.50 3.53 2.38 0.86

14400 14399 14400 14400 14400 14400 14401 14400 14401 14400 14400 14400 14400 14407 14400

O O O O O O O F O F F F F F F

0.01 0.01 0.01 0.00 0.01 0.01 0.01 1.68 0.01 0.88 3.44 0.58 4.27 1.47 2.3

79 122 269 19 8 105 161 14399 8992 14399 14399 14399 14399 14399 14403

5.58 11.05 7.03 4.39 10.08 3.03 9.60 4.57 4.58 7.70 9.41 5.94 8.89 8.07 8.41

4076 2786 2059 2825 3043

F F F F F

32.30 29.33 23.25 24.80 36.45

31.57 28.77 22.68 24.02 35.50

2.31 1.95 2.51 3.25 2.68

0.69 14400 5.46 14399 1.61 14399 4.09 14399 2.38 14399 Improvement

6.32 9.55 5.12 11.45 4.22 7.43

250

300

200

250

300

300

ED

20

|K|

CR IP T

FPVRP_n10k5t5_1 FPVRP_n10k4t5_2 FPVRP_n10k5t5_3 FPVRP_n10k4t5_4 FPVRP_n10k8t5_5 FPVRP_n15k10t5_1 FPVRP_n15k6t5_2 FPVRP_n15k10t5_3 FPVRP_n15k8t5_4 FPVRP_n15k7t5_5 FPVRP_n20k10t5_1 FPVRP_n20k12t5_2 FPVRP_n20k11t5_3 FPVRP_n20k10t5_4 FPVRP_n20k10t5_5

|C|

AN US

r

M

Instance

14400 14400 14400 14400 14400

F F F F F Average

Analysis of the FPVRP benefits

CE

7.3.3

PT

Table 6: Comparison between FPVRP and PVRP formulations using S3 instances.

In this section we focus on the analysis of the benefits of the FPVRP with respect to the FPVRP-IC and to the PVRP that emerged in the computational experiments described in

AC

the previous section. To this aim, we present a summary of the improvements of the objective function achieved by the FPVRP relative to different problem parameters. We start by focusing on the comparison between FPVRP and FPVRP-IC on the instances

of set S1. Results are shown in Table 7 where instances are clustered on the basis of the number of customers, first, and fleet size, second. For each cluster of instances we report the average and maximum improvement achieved by the FPVRP.

23

ACCEPTED MANUSCRIPT

Avrg. Improvement 1.73% 3.22% 3.07% 2.23% Avrg. Improvement 3.11% 2.02%

Max. Improvement 8.87% 6.37% 9.48% 7.91% Max. Improvement 8.87% 9.48%

CR IP T

|C| 5 10 15 20 |K| 2 3

Table 7: Improvement of FPVRP versus FPVRP-IC on S1 instances.

The table shows that the improvements increase with the number of customers except for the case with |C| = 20 while they slightly decrease when the number of vehicles increases.

The maximum improvement is 9.48%.

AN US

We now analyze the benefits of the FPVRP with respect to the PVRP. Table 8 summarizes the results for the instances of set S2. Instances are clustered by number of customers only as all instances have the same number of vehicles. We do not consider instance Instn50t5k4 in this analysis.

M

Avrg. Improvement 0.00% 0.16% 0.43% 0.21%

Max. Improvement 0.00% 0.57% 0.67% 0.40%

ED

|C| 7 9 11 15

Table 8: Improvement of FPVRP versus PVRP on S2 instances.

PT

The results show that the FPVRP benefits are quite fluctuating with respect to the number of customers. Nevertheless, FPVRP still shows slight improvements with respect to the PVRP.

CE

Finally, Table 9 focuses on the comparison between FPVRP and PVRP on the instances of Set S3. Instances are clustered by the number of customers and radius. We did not consider

AC

the number of vehicles for these instances as it varies with the number of customers. |C| 10 15 20 r 0.15 0.30 0.50

Avrg. Improvement 7.11% 7.05% 7.87% Avrg. Improvement 7.69% 7.22% 7.33%

Max. Improvement 11.05% 12.27% 11.45% Max. Improvement 12.27% 11.05% 11.45%

Table 9: Improvement of FPVRP versus PVRP on S3 instances. 24

ACCEPTED MANUSCRIPT

We can notice that, as for instances of set S2, the benefits are quite fluctuating with respect to the number of customers and the radius. The maximum savings are 12.27%. When comparing the benefits of the FPVRP versus the FPVRP-IC and PVRP, respectively, we notice that much larger improvements are achieved by the FPVRP with respect to PVRP than with respect to FPVRP-IC. This is clearly due to the way the instances are constructed and also due to the fact that PVRP is suffering from the rigidity derived from the fixed

CR IP T

schedules and fixed delivery quantities. Part of this rigidity is overcome with the FPVRP-IC. However, the improvements of the FPVRP versus FPVRP-IC show that inventory constraints may remarkably affect solution costs.

Finally, we made additional computational experiments in order to investigate the impact of customer demands on the benefits produced by the PVRP. The results confirm that the highest

AN US

savings can be attained when the capacity is tight and demands are rather homogeneous, whereas the smallest savings are obtained with larger capacities and high variations in the demands. Details about these tests can be found at http://or-brescia.unibs.it/instances. A similar observation was made in Archetti et al. (2008) concerning the benefits of the SDVRP with respect to the VRP.

Conclusions and Future Research

M

8

In this paper we have introduced the Flexible Periodic Vehicle Routing Problem (FPVRP)

ED

where a carrier has to establish a distribution plan to serve his customers over a planning horizon, using a fleet of homogeneous capacitated vehicles. Each customer has a total demand

PT

that must be served within the horizon, although the quantity delivered at each visit cannot exceed a fixed quantity. The objective is to minimize the total routing cost over the time horizon. The FPVRP is related to other VRPs in which customers have to be served with

CE

a certain regularity, like the Periodic Vehicle Routing Problem (PVRP) and the Inventory Routing Problem (IRP). We have presented a worst-case analysis which shows the theoretical

AC

advantages of the FPVRP with respect to both PVRP and IRP. Furthermore, we have presented a load-based mixed integer formulation for the FPVRP, together with some families of valid inequalities. The results of the computational experiments that have been carried out show that the FPVRP may produce substantial improvements in the routing costs in comparison with both PVRP and IRP. Given the difficulty of optimally solving FPVRP instances of relatively small sizes, a promising avenue for future research is to develop an efficient heuristic able to produce high quality solutions in small computing times.

25

ACCEPTED MANUSCRIPT

Acknowledgements The research of the last two authors has been partially supported by the Spanish Ministry of Economy and Competitiveness end EDRF funds through grant MTM2015-63779-R (MINECO/FEDER). The research of the third author has been partially supported by the acknowledged.

CR IP T

Mexican National Council for Science and Technology (CONACyT). This support is gratefully The authors are grateful to Karen Smilowitz for providing the necessary information to reproduce the set of instances S2 used in this paper.

Thanks are due to both referees for their constructive comments, which have led to a considerable improvement of the manuscript.

AN US

References

Aksen, D., Kaya, O., Salman, F. S., and Akça, Y. (2012). Selective and periodic inventory routing problem for waste vegetable oil collection. Optimization Letters, 6(6):1063–1080. Archetti, C., Bianchessi, N., Irnich, S., and Speranza, M. G. (2014). Formulations for an

M

inventory routing problem. International Transactions in Operational Research, 21(3):353– 374.

Archetti, C., Jabali, O., and Speranza, M. G. (2015). Multi–period vehicle routing problem

ED

with due dates. Computers and Operations Research, 61:122–134. Archetti, C., Savelsbergh, M., and Speranza, M. G. (2008). To split or not to split: That is

PT

the question. Transportation Research E, 44:114–123. Archetti, C. and Speranza, M. (2012). Vehicle routing problems with split deliveries. Inter-

CE

national Transactions in Operational Research, 39:3–22. Baldacci, R., Bartolini, E., Mingozzi, A., and Valleta, A. (2011). An exact algorithm for the period routing problem. Operations Research, 59(1):228–241.

AC

Bell, W. J., Dalberto, L. M., Fisher, M. L., Greenfield, A. J., Jaikumar, R., Kedia, P., Mack, R. G., and Prutzman, P. J. (1983). Improving the distribution of industrial gases with an on-line computerized routing and scheduling optimizer. Interfaces, 13(6):4–23.

Beltrami, E. J. and Bodin, L. D. (1974). Networks and vehicle routing for municipal waste collection. Networks, 4(1):65–94. Bertazzi, L. and Speranza, M. (2012). Inventory routing problems: an introduction. EURO Journal on Transportation and Logistics, 1:307–326.

26

ACCEPTED MANUSCRIPT

Bertazzi, L. and Speranza, M. (2013). Inventory routing problems with multiple customers. EURO Journal on Transportation and Logistics, 2:255–275. Cacchiani, V., Hemmelmayr, V. C., and Tricoire, F. (2014). A set–covering based heuristic algorithm for the periodic vehicle routing problem. Discrete Applied Mathematics, 163, Part 1:53–64. Matheuristics 2010.

CR IP T

Campbell, A. M. and Wilson, J. H. (2014). Forty years of periodic vehicle routing. Networks, 63(1):2–15.

Christofides, N. and Beasley, J. E. (1984). The period routing problem. Networks, 14(2):237– 256.

Coelho, L. C., Cordeau, J.-F., and Laporte, G. (2013). Thirty years of inventory routing. Transportation Science, 48(1):1–19.

AN US

Dantzig, G. B. and Ramser, J. H. (1959). The truck dispatching problem. Managment Science, 6(1):80–91.

Escobar, J. W., Linfati, R., Toth, P., and Baldoquin, M. G. (2014). A hybrid granular tabu search algorithm for the multi-depot vehicle routing problem. Journal of Heuristics, 20(5):483–509.

M

Francis, M. P., Smilowitz, K., and Tzur, M. (2008). The period vehicle routing problem and its extensions. In Golden, B., raghavan, S., and Wasil, E., editors, The Vehicle Routing

ED

Problem: Latest Advances and New Challenges, pages 73–102. Springer, New York. Francis, P. and Smilowitz, K. (2006). Modeling techniques for periodic vehicle routing problems. Transportation Research Part B: Methodological, 40(10):872–884.

PT

Francis, P., Smilowitz, K., and Tzur, M. (2006). The period vehicle routing problem with service choice. Transportation Science, 40(4):439–454.

CE

Francis, P. M., Smilowitz, K. R., and Tzur, M. (2007). Flexibility and complexity in periodic distribution problems. Naval Research Logistics, 2(54):136–150. García Calvillo, I. D. (2010). Un Enfoque Metaheurístico para un Problema de Ruteo con

AC

Flexibilidad en las Fechas de Entrega. PhD thesis, Universidad Autónoma de Nuevo León, México.

Garey, M. R. and Johnson, D. S. (1979). Computers and Intractability: A Guide to the Theory of NP-Completeness. W. H. Freeman & Co., New York, NY, USA. Gueguen, C. (1999). ’Méthodes de résolution exacte pour les problèmes de tournées de véhicules. PhD thesis, Ecole Centrale Paris. Há, M. H., Bostel, N., Langevin, A., and Rousseau, L.-M. (2014). An exact algorithm and a

27

ACCEPTED MANUSCRIPT

metaheuristic for the generalized vehicle routing problem with flexible fleet size. Computers & Operations Research, 43:9–19. Hashimoto, H., Ibaraki, T., Imahori, S., and Yagiura, M. (2006). The vehicle routing problem with flexible time windows and traveling times. Discrete Applied Mathematics, 154(16):2271–2290.

CR IP T

Hemmelmayr, V., Doerner, K., Hartl, R., and Savelsbergh, M. W. P. (2009). Delivery strategies for blood products supplies. OR Spectrum, 31(4):707–725.

Letchford, A. N. and Salazar-González, J. J. (2015). Stronger multi-commodity flow formulations of the capacitated vehicle routing problem. European Journal of Operational Research, 244(3):730 – 738.

Michallet, J., Prins, C., Amodeo, L., Yalaoui, F., and Vitry, G. (2014). Multi-start iterated

AN US

local search for the periodic vehicle routing problem with time windows and time spread constraints on services. Computers & Operations Research, 41:196–207. Pacheco, J., Álvarez, A., García, I., and Ángel Bello, F. (2012). Optimizing vehicle routes in a bakery company allowing flexibility in delivery dates. Journal of the Operational Research Society, 63(5):569–581.

M

Rahimi-Vahed, A., Crainic, T. G., Gendreau, M., and Rei, W. (2015). Fleet-sizing for multi-depot and periodic vehicle routing problems using a modular heuristic algorithm.

ED

Computers & Operations Research, 53:9 – 23. Rusdiansyah, A. and Tsao, D. (2005). An integrated model of the periodic delivery problems for vending–machine supply chains. Journal of Food Engineering, 70(3):421–434.

PT

Russell, R. and Igo, W. (1979). An assignment routing problem. Networks, 9(1):1–17. Shmoys, D. B., Rinnooy, K. A. H. G., Lenstra, J. K., and Lawler, E. L., editors (1985). The

CE

Traveling Salesman Problem: A Guide Tour of Combinatorial Optimization. John Wiley & Sons Ltd., New York, USA. Toth, P. and Vigo, D., editors (2014). Vehicle Routing: Problems, Methods, and Applications.

AC

Society for Industrial and Applied Mathematics, Philadelphia, 2 edition.

Vidal, T., Crainic, T. G., Gendreau, M., Lahrichi, N., and Rei, W. (2012). A hybrid genetic algorithm for multidepot and periodic vehicle routing problems. Operations Research, 60(3):611–624.

28