A branch-and-cut algorithm for the inventory routing problem with pickups and deliveries

A branch-and-cut algorithm for the inventory routing problem with pickups and deliveries

Journal Pre-proof

A branch-and-cut algorithm for the inventory routing problem with pickups and deliveries Claudia Archetti, M. Grazia Speranza, Maurizio Boccia, Antonio Sforza, Claudio Sterle PII: DOI: Reference:

S0377-2217(19)30813-6 https://doi.org/10.1016/j.ejor.2019.09.056 EOR 16084

To appear in:

European Journal of Operational Research

Received date: Accepted date:

29 January 2019 27 September 2019

Please cite this article as: Claudia Archetti, M. Grazia Speranza, Maurizio Boccia, Antonio Sforza, Claudio Sterle, A branch-and-cut algorithm for the inventory routing problem with pickups and deliveries, European Journal of Operational Research (2019), doi: https://doi.org/10.1016/j.ejor.2019.09.056

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2019 Published by Elsevier B.V.

Highlights • Formulation of the multi-vehicle inventory routing problem with pickup and delivery • Proposal of a branch-and-cut algorithm using known and new valid inequalities • Proposal of a new class of valid inequalities, called interval inequalities • Extension of single-item lot-sizing valid inequalities to the pickup nodes • Proposed method optimally solve most of the single/multi vehicle instances

1

A branch-and-cut algorithm for the inventory routing problem with pickups and deliveries Claudia Archetti Department of Information Systems, Decision Sciences and Statistics, ESSEC Business School, Paris, France [email protected]

M. Grazia Speranza Dept. of Economics and Management, University of Brescia, Italy [email protected]

Maurizio Boccia; Antonio Sforza; Claudio Sterle* Dept. of Electrical Engineering and Information Technology University of Naples, Italy * Istituto di Analisi dei Sistemi ed Informatica A. Ruberti, IASI-CNR, Rome, Italy [email protected]; [email protected]; [email protected]

Abstract This paper addresses the Inventory Routing Problem with Pickups and Deliveries (IRP-PD). A single commodity has to be picked up from several origins and distributed to several destinations. The commodity is made available at the supplier depot and at the pickup customers, with known amounts in each period of a discretized time horizon. The commodity is consumed by the delivery customers, whose demand, in each period of the time horizon, is known. The vehicles start and end their routes at the supplier depot. The objective is to determine a collection and distribution plan minimizing the sum of routing and inventory costs, satisfying inventory and capacity constraints. A mixed integer linear programming model is presented for the IRP-PD. Valid inequalities are proposed that are adapted from the inventory routing, the vehicle routing and the lot-sizing problems literature. Moreover, a new set of valid inequalities, called interval inequalities, is presented. A branch-and-cut algorithm is designed and tested on a large set of instances. Computational results show that the proposed approach is able to solve to Preprint submitted to EJOR

September 30, 2019

optimality 946 out of 1280 instances. Moreover, the approach outperforms the state-of-the art method for the single vehicle case, increasing the number of instances solved to optimality from 473 to 538 out of 640 instances tested, with 133 new best known solutions. Keywords: Routing, inventory management, pickups and deliveries, branch-and-cut, lot-sizing reformulation. 1. Introduction The most commonly investigated Inventory Routing Problem (IRP ) consists in the distribution of a commodity from a supplier to a set of customers over a discrete planning horizon with a homogeneous fleet of capacitated vehicles. The problem simultaneously tackles two different decision problems, i.e., vehicle routing and inventory management. The commodity is made available at the supplier, which could be a factory or a warehouse, and consumed by the customers. The great interest in this problem and its variants is motivated by applications arising in different contexts and mainly in supply chain management. The class of the IRP is extended to the class of the IRP with Pickups and Deliveries (IRP-PD) when vehicles pick up one or more products from multiple origins (pickup customers) and deliver it to multiple destinations (delivery customers). The IRP-PD has received little attention in the literature. A number of applications of the IRP-PD are mentioned in [2]. The most interesting ones are in rebalancing inventory problems in store chains, where some stores have an excess of inventory and others are running out of stock. Vehicles perform pickup and delivery operations aimed at rebalancing the inventory levels. In the sharing and collaborative economy, the IRP-PD finds also application in situations where a carrier, with a fleet of vehicles, picks up goods from producers and delivers goods to consumers over a period of time. Depending on the specific problem under investigation, the role of each customer, pickup or delivery, can be the same or change over the planning horizon. If the duration of the planning horizon is small, we can assume that the role of each customer does not change. If pickup and delivery customers correspond to plants and retailers, respectively, their role of pickup or delivery customers remains obviously unchanged during all the time horizon. On the contrary, if all the customers are retailers which can exchange a product depending on their stock levels, as it occurs in store chains, then 3

their role can change over the time horizon. In this paper, a problem belonging to the class of the IRP-PD is investigated. It extends to the case of multiple vehicles the problem introduced in [2], which considers an IRP-PD with no change of the customer role during the planning horizon, one commodity and one vehicle (1-1-IRP-PD). Hence, we tackle the single-commodity multi-vehicle version of the problem (1-MIRP-PD). A Mixed Integer Linear Programming (MILP ) model is presented together with an original branch-and-cut solution approach. More specifically, the main contributions of this paper can be summarized as follows: 1. The formulation proposed by [2] for the 1-1-IRP-PD is extended to the multiple vehicle case 1-M-IRP-PD in an aggregated form which does not suffer from the increasing dimension related to the number of vehicles. 2. The proposed formulation is improved by exploiting the single-item lotsizing substructure of the IRP. In particular, the capacitated lot-sizing formulation presented in [4] for delivery customers is extended to the case of pickup customers taking into account the different sequence of operations to be performed. Then, valid inequalities for the resulting polyhedron are proposed. 3. Different classes of valid inequalities are presented. Some of the valid inequalities are inherited from the 1-1-IRP-PD or the Vehicle Routing Problem (VRP) literature. Lot-sizing inequalities presented in [4] are adapted to the pickup and delivery case. In addition, a new class of inequalities, called interval inequalities, is proposed and its validity of is demonstrated. These inequalities are defined considering separately pickup and delivery customers. Thus, they can be used also in the classical IRP. 4. The formulation and the valid inequalities are tested on a large set of benchmark instances. The results show that the valid inequalities improve the lower bound at the root node of the branch-and-bound tree by 24.84% on average. For the case with multiple vehicles, the branchand-cut algorithm solves to optimality 946 out of 1280 tested instances with up to 50 customers and a horizon of 3 periods or 30 customers and a horizon of 6 periods. For the single vehicle case, the branch-and-cut algorithm solves to optimality 538 out of 640 instances, substantially increasing the number of 473 instances solved to optimality in [2]. The 4

average optimality gap is 0.89%, which improves the previous 1.22%. Moreover, 133 new best known solutions are obtained. The rest of the paper is organized as follows. Section 2 provides a literature review of the IRP-PD. The 1-M-IRP-PD is described and formulated in Section 3, where valid inequalities are also presented. In Section 4 we describe the branch-and-cut algorithm, while Section 5 is devoted to computational tests and results. Some conclusions are drawn in Section 6. 2. Literature review The class of the IRP, where only delivery operations are performed, has been widely studied in the literature. A review on the IRP is out of the scope of this work. Hence, we just address the interested reader to the following surveys and recent works on the topic: [9, 10, 14, 5]. The same observation is valid also for pickup and delivery problems PDP, for which we just refer to the surveys [8, 7] and the recent works published in [20, 29]. On the contrary, the literature on the IRP-PD that considers both pickup and delivery operations is relatively scarce. Most of the contributions on the IRP-PD are related to maritime applications, in particular to the transportation of liquid and bulk materials. Given the size and the capacity of the ships, the quantities to be transported and the extent of the distances and of the traveling times, planning routes performing both pickup and delivery operations at the visited ports is much more efficient than separately planning the two kinds of operations. Recent contributions on the IRP-PD arising in maritime transportation can be found in [17, 16, 1, 15]. The reader is referred to [13] and [12] for an overview of maritime inventory routing problems with pickup and delivery operations. Although relevant applications of the IRP-PD can be found in road transportation, to the best of our knowledge, the first contribution on the IRP-PD appeared in [6]. In this paper the authors considered an IRP characterized by the presence of satellite facilities where vehicles can be reloaded to continue the delivery operations to the customers until a route duration constraint is reached. They proposed three heuristic decomposition approaches for the problem. Almost 10 years later, in [26] the authors tackled an IRP-PD defined on wide geographic areas and for a long planning horizon, with delivery routes including also pickup points. The studied problem models the distribution 5

of argon where customers have to be served by vehicles which continuously move on the network and pickup products at different facilities with limited production capacity. The authors proposed a randomized greedy algorithm, later refined in [27], where an integer multi-commodity flow formulation on a time-expanded network and a branch-and-cut algorithm were proposed. In [25] an IRP-PD was considered where pickup operations of different commodities have to be performed before the delivery to the final customers. The possibility of delivery splitting is also considered. An integer linear programming formulation was proposed and solved on small size instances. The recent contribution [28], motivated by a real-life application to the replenishment of automated teller machines, studies a problem where each node can be both a pickup and a delivery node, depending on the specific period of the planning horizon, as money may have to be delivered to or picked up from an automated teller machine. Finally, IRP-PD applications have been studied in [19], [23] and [22]. In particular, in [19], the authors tackle a particular IRP-PD problem arising in the management of the returnable transport items in a two stage supply chain system. The authors provide a MILP model, taking into account operational constraints and peculiarities, and propose a clustering based heuristic to solve realistic size instances and perform a scenario analysis. In [22] and [23], the authors tackle an IRP-PD arising in the industrial gas supply chain and related reverse logistics. Pick-up and delivery operations involve liquefied industrial gases and empty tanks backhauling, respectively. Both works provide a deterministic mathematical programming formulation for the related IRP-PD and exploit the Economic Model Predictive Control (EMPC) to deal with demand uncertainties. As already mentioned, the single-commodity and single-vehicle problem (1-1-IRP-PD) was introduced in [2], where a commodity has to be picked up from pickup customers and delivered to delivery customers over a given planning horizon with a single capacitated vehicle. The vehicle starts and ends its routes at the supplier depot, where the commodity can be stocked and the role of each customer does not change over all the planning horizon. The objective is to determine a distribution and collection plan minimizing the sum of routing cost and inventory cost at the customers. In the following sections of this paper we study the extension of this problem to the multiple vehicle case.

6

3. Problem description and formulation The 1-M-IRP-PD is defined on a network G = (N 0 , A), where N 0 and A are the node and the arc sets. More precisely, N 0 = {0} ∪ N , where 0 is the depot, and N = NP ∪ ND is the set of nodes representing pickup customers (NP ) and delivery customers (ND ). A discrete planning horizon T = {1, ..., T } is considered. Let us also define T 0 = T ∪ {0}. Five quantities are associated with each node i ∈ N : dit , hi , Ii0 , Li and Ui , where dit is the constant quantity made available (if i ∈ NP ) or demanded (if i ∈ ND ) at t ∈ T , hi is the unitary holding cost, Ii0 , Li and Ui are the initial inventory level, the lower and upper inventory limits, respectively. An initial inventory level I00 is associated with the depot, without lower and upper limits. Finally, a transportation cost cij is associated with each arc (i, j) ∈ A and a fleet of m homogeneous vehicles with capacity Q is available. Each vehicle can perform at most one route starting and ending at the depot in each time period t, t ∈ T . The quantities to be picked up or delivered cannot be split. The objective is to determine the collection/distribution plan that minimizes the total cost, which is given by the sum of the routing cost and the inventory cost. No stockout or backlog is allowed. The proposed formulation for the 1-M-IRP-PD makes use of the following decision variables: 1. continuous variables qit identifying the quantity picked up (delivered) at node i ∈ NP (i ∈ ND ), at time t ∈ T ; 2. binary variables yit taking value 1 if the node i ∈ N is visited at time t ∈ T , 0 otherwise; 3. integer variables y0t counting the number of vehicles used (routes performed) at time t ∈ T ; 4. continuous variables Iit representing the inventory level at node i ∈ N 0 at time t ∈ T 0 ; 5. binary variables xijt taking value 1 if arc (i, j) ∈ A is traversed at time t ∈ T , 0 otherwise; 6. continuous variables lijt representing the quantity on load of vehicles traversing arc (i, j) ∈ A at time t ∈ T . The flow formulation for the 1-M-IRP-PD, based on the notation summarized in Table 1, is as follows:

7

Table 1: IRP-PD notation Sets N N0 NP , ND A T T0

Set of customer nodes Set of all nodes, including the depot, N 0 = N ∪ {0} Set of pickup and delivery customers, respectively, NP ∪ ND = N , NP ∩ ND = ∅ Set of arcs (i, j), i, j ∈ N 0 Set of time periods T ∪ {0}

Parameters Demand of customer i ∈ N at time period t ∈ T Unitary holding cost at node i ∈ N 0 Lower and upper inventory limit at customer i ∈ N , respectively Initial inventory level at node i ∈ N 0 Transportation cost associated with arc (i, j) ∈ A Fleet size Vehicle capacity

dit hi Li , Ui Ii0 cij m Q Variables

Quantity picked up (delivered) at customer i ∈ N at time t Equals to 1 if node i ∈ N is visited at time t ∈ T , 0 otherwise Equals to the number of vehicles leaving the depot at time t ∈ T Inventory level at node i ∈ N 0 at time t ∈ T Equals to 1 if arc (i, j) ∈ A is used at time t ∈ T , 0 otherwise Quantity loaded on the vehicle traversing arc (i, j) ∈ A at time t ∈ T

qit yit y0t Iit xijt lijt

min z =

X X

XX

cij xijt +

X

j∈N

(1a)

x0jt − y0t = 0,

X

j∈N 0

X

hi Iit ,

i∈N 0 t∈T

(i,j)∈A t∈T

j∈N 0

xijt −

X

xjit = 0,

j∈N 0

xijt − yit = 0,

XX i∈S j∈S

xijt ≤

X i∈S

yit − ymt ,

Iit − Ii,t−1 − dit + qit = 0,

Iit − Ii,t−1 + dit − qit = 0, X X I0t − I0,t−1 + qit − qit = 0, i∈N D

i∈N P

8

t∈T,

(1b)

i ∈ N 0, t ∈ T ,

(1c)

i ∈ N,t ∈ T ,

(1d)

S ⊆ N , m ∈ S, t ∈ T ,

(1e)

i ∈ N P, t ∈ T ,

(1f)

D

i ∈ N ,t ∈ T ,

(1g)

t∈T,

(1h)

i ∈ N0

Ii0 = Ii0 ,

i ∈ N D, t ∈ T ,

Ii,t−1 + qit ≤ Ui , Iit ≥ Li ,

X

j∈N 0

X

j∈N

X

i ∈ N ,t ∈ T ,

(1l)

i ∈ N P, t ∈ T ,

(1q)

i ∈ N D, t ∈ T ,

(1r)

j ∈ N,t ∈ T ,

(1s)

i ∈ N , t ∈ T , (1m) i ∈ N , t ∈ T , (1n) i ∈ N , t ∈ T , (1o) (i, j) ∈ A, t ∈ T , (1p)

j∈N 0

ljit − qit −

(1k)

P

P

Ii,t−1 − qit ≥ Li , qit ≤ min{Q, Ui − Li }yit , yit ≤ y0t ≤ m, lijt ≤ Qxijt , X X ljit + qit − lijt = 0, j∈N 0

(1j)

D

i ∈ N ,t ∈ T ,

Iit ≤ Ui ,

(1i)

lijt = 0,

j∈N 0

l0jt ≤ I0,t−1 ,

xijt ∈ {0, 1}, lijt ≥ 0, yit ∈ {0, 1}, qit ≥ 0, Iit ≥ 0, y0t ∈ N,

(i, j) ∈ A, t ∈ T , (1t) i ∈ N , t ∈ T , (1u) i ∈ N 0 , t ∈ T , (1v) t ∈ T . (1w)

The objective function (1a) is the sum of inventory and routing costs. Constraints (1b) compute the number of vehicles (routes) leaving the depot at each time period, while constraints (1c) and (1d) are related to route construction. In particular, (1c) are flow conservation constraints, while (1d) impose that the number of arcs leaving from a customer node is equal to 1 if the node is visited, 0 otherwise. As yit is binary for i ∈ N , this means that at most one visit is performed at each customer in each time period. Constraints (1e) are subtour elimination constraints. Constraints (1f), (1g) and (1h) are inventory balance constraints for pickup nodes, delivery nodes and the depot, respectively. Constraints (1i) impose the initial inventory conditions. Constraints (1j) and (1k), and constraints (1l) and (1m), define lower and upper inventory limits for delivery and pickup nodes, respectively. They are consistent with the following sequence of operations: first the node is served by a vehicle performing either delivery or pickup operations without exceeding the upper and lower inventory limits; then, either the demand or the supply is satisfied; finally, the inventory level is updated. Constraints 9

(1n) impose the limits on the quantity that can be picked up or delivered to a node. Constraints (1o) link variables yit , i ∈ N and t ∈ T , with y0t . Moreover, they impose that the number of routes in each time period does not exceed the number of available vehicles. Constraints (1p) are vehicle capacity constraints. Constraints (1q) and (1r) are load balance constraints, while constraints (1s) limit the quantity distributed from the depot with respect to the available inventory level. Finally, constraints (1t), (1u), (1v) and (1w) define the domain of the variables. Note that formulation (1) differs from the one proposed in [2] for the single vehicle case in constraints (1o) and (1s), which handle the multiple vehicle case. Different valid inequalities may be added to enhance formulation (1). In Section 3.1 the inequalities proposed in [2] for the 1-1-IRP-PD, which remain valid for the multiple vehicle case, are reported. In Sections 3.2 and 3.3 inequalities derived from the lot-sizing reformulation of the problem for delivery and pickup nodes, respectively, are described. In Section 3.4 Rcut inequalities are presented which have been proposed in the literature for different routing problems. Finally, in Section 3.5 a new class of valid inequalities is proposed. 3.1. Valid inequalities from 1-1-IRP-PD The following inequalities were proposed in [2] for the single vehicle problem 1-1-IRP-PD and are valid for the multiple vehicle case too. Ii(t−k−1) t X

j=1

yij ≥

Ii(t−k−1)



≥

k X

j=0

& Pt−1

j=1



di(t−j)  1 − dij − Ii0 + Li

'

k X

j=0

yij ≥

yi(t−j)  + Li

min{Q, Ui − Li }    k k X X ≤ Ui −  di(t−j)  1 − yi(t−j)  j=0

t X



& Pt−1

j=1

dij + Ii0 − Ui

j=0

'

min{Q, Ui − Li } P X X j∈N D qjt xijt ≥ Q i∈N P ∪{0} j∈N D P X X i∈N D qit xijt ≥ Q i∈N D j∈N P ∪{0} P X X j∈N P qjt xijt ≥ Q D P

j=1

i∈N

∪{0} j∈N

10

i ∈ N D , t ∈ T , k = 0, . . . , t − 1

(2a)

i ∈ N D, t ∈ T

(2b)

i ∈ N P , t ∈ T , k = 0, . . . , t − 1,

(2c)

i ∈ NP , t ∈ T ,

(2d)

t∈T,

(2e)

t∈T,

(2f)

t∈T,

(2g)

X

X

i∈N P j∈N D ∪{0}

xijt ≥

P

i∈N P

qit

t∈T.

Q

(2h)

Inequalities (2a)-(2b) and (2c)-(2d) are related to the inventory level and the number of visits to pickup and delivery nodes, respectively. Constraints (2e)-(2h), performing a sort of shrinking of all the pickup/delivery nodes in a single node [18], link the variables representing the quantities picked up/delivered to the arc variables. 3.2. Lot-sizing inequalities for constant-demand delivery nodes In [4] the authors addressed the single vehicle IRP and presented a tight reformulation for the lot-sizing subproblems arising for each delivery node. In addition, they presented a set of valid inequalities exploiting the lot-sizing structure of the subproblem. This reformulation, together with the valid inequalities, can be applied to the delivery nodes in the 1-M-IRP-PD. We now report the valid inequalities for completeness. They are defined for any delivery node i, i ∈ ND , that satisfies the following conditions: i) node i has constant demand, i.e, dit = di , t ∈ T ; ii) the difference between the upper bound Ui and the lower bound Li on the inventory level is a multiple of di , i.e., Ui − Li = V di . In the following, we show the inequalities for the cases V ∈ {2, 3}. The inequalities can be generalized to any value of V . We note that, given the sequence of the operations at the delivery nodes, the following condition on the inventory level is verified: Li ≤ Iit ≤ Li + (V − 1)di , i ∈ ND , t ∈ T . The inequalities for each node i, i ∈ ND are the following. Case V = 3

Iit ≤ Li + di (1 + yit ), Iit−1 + di yit ≥ Li + di , Iit−1 + 3di yit ≥ di + Iit , di (yit−1 + 2yit ) + Li ≥ Iit , yit−2 + yit−1 + yit ≥ 1, Iit−2 + di (2yit−1 + yit ) ≥ Li + 2di , Iit+1 + di (yit−1 + yit ) ≥ Iit , 11

t∈T, t∈T, t∈T, t ∈ T , t ≥ 2, t ∈ T , t ≥ 3, t ∈ T , t ≥ 2, t ∈ T , 2 ≤ t ≤ T − 1,

(3a) (3b) (3c) (3d) (3e) (3f) (3g)

Iit−1 + Iit+1 + 2di yit ≥ di + Iit + Li , di + Iit + 2di yit+1 ≥ Iit−1 + Iit+1 − Li ,

t ∈ T , t ≤ T − 1, t ∈ T , t ≤ T − 1.

(3h) (3i)

Case V = 2 di yit + Iit−1 ≥ Li + di , di yit + Li ≥ Iit ,

t∈T, t∈T.

(4a) (4b)

3.3. Lot-sizing inequalities for constant-demand pickup nodes The inequalities presented in Section 3.2 for delivery customers can be adapted to the case of pickup customers. To this aim, let us consider the difference between the operations performed at pickup and delivery nodes, respectively. At delivery nodes, the demand consumes the quantity in inventory. This quantity is then refilled by the deliveries which are such that the inventory level never exceeds the maximum level. On the contrary, at pickup nodes, the production increases the quantity in inventory. This quantity is collected by the vehicles before exceeding the maximum inventory level and the quantity left in inventory is at least equal to the minimum level. Let us consider a pickup node i with constant production rate, i.e., dit = di , t ∈ T . For such node, we have the lot-sizing subproblem defined by inequalities (1f), (1l), (1m) and (1n). Under the assumption Ui − Li = V di , this subproblem can be rewritten for each i, i ∈ NP , as follows: t∈T, t∈T, t∈T.

Ii,t−1 + di = qit + Ii,t Ii,t ≤ V di + Li qit ≤ V di yit

(5a) (5b) (5c)

On the basis of the lot-sizing reformulation for pickup nodes, we obtain valid inequalities for pickup customers by replacing (Iit −Li ) with (V di −Iit + Li ) in constraints (3a) - (4b). This replacement can be explained by observing that (Iit −Li ) represents the amount available in the delivery node i to satisfy its own demand in the successive periods. Likewise, the pickup node i can supply at most the quantity (Ui −Iit ) in the successive periods. Recalling that Ui = V di +Li , the replacement can be obtained by substitution. We also note that, given the sequence of the operations at the pickup nodes, the following condition on the inventory level is verified: Li + di ≤ Iit ≤ Li + V · di , i ∈ NP , t∈T. 12

Case V = 3 Iit + di yit ≥ Li + 2di , Iit−1 ≤ Li + 2di + di yit , Iit + 3di yit ≥ di + Iit−1 , di (yit−1 + 2yit ) + Iit ≥ 3di + Li , yt−2 + yt−1 + yt ≥ 1, Iit−2 ≤ Li + di (1 + 2yt−1 + yt ), Iit + di (yit−1 + yit ) ≥ Iit+1 , It−1 + It+1 ≤ 2di (1 + yit ) + Iit + Li , Iit−1 + Iit+1 + 2di yit+1 ≥ 2di + Iit + Li ,

t∈T, t∈T, t∈T, t ∈ T , t ≥ 2, t ∈ T , t ≥ 3, t ∈ T , t ≥ 2, t ∈ T , t ≤ T − 1, t ∈ T , 2 ≤ t ≤ T − 1, t ∈ T , t ≤ T − 1.

(6a) (6b) (6c) (6d) (6e) (6f) (6g) (6h) (6i)

t∈T t∈T

(7a) (7b)

Case V = 2 Iit−1 ≤ Li + di yit + di , di yit + Iit ≥ Li + 2di ,

3.4. R-cut inequalities In this section valid inequalities are proposed that are derived from the R-cut constraints already used in the literature for different routing problems [21, 3, 11]. Let us consider a time period t ∈ T and two nodes i, j ∈ N . Let us denote by C(i, j) a cut between i and j, i.e., C(i, j) = {(u, v) ∈ A|i, u ∈ S, j, v ∈ S 0 , S ∪ S 0 = N , S ∩ S 0 = ∅}. Moreover, let C(0, {i, j}) be a cut between the depot and nodes {i, j}, i.e., C(0, {i, j}) = {(u, v) ∈ A|0, u ∈P S, i, j, v ∈ S 0 , S∪ S 0 = N 0 , S ∩ S 0 = ∅}. Finally, let us define X t (C(S, S 0 )) = (u,v)∈C(S,S 0 ) xuvt . Then, the following inequality: X t (C(0, {i, j})) + X t (C(i, j)) + X t (C(j, i)) ≥ yit + yjt

(8)

is valid for the 1-M-IRP-PD. Indeed, when yit + yjt ≤ 1, then X t (C(0, {i, j})) ≥ yit + yjt and the inequality is satisfied. If yit + yjt = 2, then i and j may be served by two vehicles (X t (C(0, {i, j})) ≥ 2), or by a single vehicle (X t (C(0, {i, j})) ≥ 1 and X t (C(i, j)) + X t (C(j, i)) ≥ 1).

13

We may couple inequalities (8) with the lot-sizing structure of constantdemand nodes presented in Sections 3.2 and 3.3. Let us consider two customers i and j with constant demand and for which Ui − Li = 2di and Uj − Lj = 2dj . Then, the inequality: X t (C(0, {i, j})) + X t (C(i, j)) + X t (C(j, i)) ≥ 2 − yi(t−1) − yj(t−1)

(9)

is valid for the 1-M-IRP-PD. In fact, if yi(t−1) + yj(t−1) = 0, then yit + yjt = 2. If instead only one among yi(t−1) and yj(t−1) is equal to 0, say yi(t−1) , then yit = 1 and X t (C(0, {i, j})) ≥ 1. Note that inequalities (9) can be easily generalized to the case where Ui − Li = V di and Uj − Lj = V dj , V > 2. However, in the experiments presented in Section 5, we used only the inequalities for the case V = 2. Inequalities (8) can be strengthened for the single vehicle case (1-1-IRPPD) as follows: X t (C(i, j)) + X t (C(j, i)) ≥ yit + yjt − y0t

(10)

for each i, j ∈ N . Note that, if y0t = 0 then yit = yjt = 0 and the inequality is satisfied. If y0t = 1 and one of the two nodes is not served in t then yit + yjt − y0t ≤ 0. If y0t = 1 and the two nodes are served in t, at least one arc in C(i, j) ∪ C(j, i) must be traversed. As for inequalities (8), we may couple inequalities (10) with the lot-sizing structure of constant-demand nodes. In particular, consider two customers i and j with constant demand and for which Ui − Li = 2di and Uj − Lj = 2dj . Then, the inequality: X t (C(i, j)) + X t (C(j, i)) ≥ 1 − yi(t−1) − yj(t−1)

(11)

is valid for the 1-1-IRP-PD. Indeed, if yi(t−1) −yj(t−1) = 0, then the two nodes must be served in t and at least one arc in C(i, j) ∪ C(j, i) must be traversed. As for inequalities (9), also (11) can be easily generalized to the case where Ui − Li = V di and Uj − Lj = V dj , V > 2. 3.5. Interval inequalities In this section we present a new class of valid inequalities, called interval inequalities, derived specifically for the 1-M-IRP-PD. The interval inequalities are a generalization of inequalities (2b) and (2d) presented above. While 14

(2b) and (2d) count the number of times a delivery and a pickup node have to be visited within the first time period and any period t of the time horizon, interval inequalities consider any pair of time periods. We first present them for delivery nodes and then we discuss their adaptation to pickup nodes. 3.5.1. Interval inequalities for delivery nodes Let us consider a delivery node i and two periods t1 and t2 , with 1 ≤ t1 < t2 ≤ T . The following step function &
' Dtt12 − Ii(t1 −1) − Li i , (12) ηt1 t2 (Ii(t1 −1) ) = K P2 where Dtt12 = tt=t dit and K = min{Q, Ui − Li }, computes the minimum 1 number of times the delivery node i has to be visited in the time interval [t1 , t2 ] if the inventory level at time t1 − 1 is equal to Ii(t1 −1) . Let X = {¯ x1 , x¯2 , ..., x¯|X| }, with x¯h < x¯h+1 for each h = 1, ..., |X| − 1, be the set of the step points of function ηt1 t2 . More formally, X = {Ii(t1 −1) |Li ≤ Ii(t1 −1) ≤ Ui ∧ ηti1 t2 (Ii(t1 −1) ) = ηti1 t2 (Ii(t1 −1) − ) + 1, 0 <  ≤ K}. Proposition 3.1. The step function ηti1 t2 for delivery nodes satisfies the following properties: i) x¯1 − Li = Dtt12 − K(ηti1 t2 (Li ) − 1) ≤ K, ii) x¯h = K(h − 1) + x¯1 for each h ∈ {1, ..., |X|}, iii) ηti1 t2 (¯ xh ) = ηti1 t2 (Li ) − h for each h ∈ {1, ..., |X|}. Proof. i) ηti1 t2 (Li ) is the maximum value of function ηti1 t2 which corresponds to the case where Ii(t1 −1) = Li and, thus, a quantity at least equal to Dtt12 has to be transported to i in [t1 , t2 ]. x¯1 corresponds to the inventory level for which one vehicle less than ηti1 t2 (Li ) is needed to serve i in [t1 , t2 ] and is equal to Dtt12 − K(ηti1 t2 (Li ) − 1) + Li . By construction, the quantity x¯1 − Li is lower than K. ii) Given the value x¯1 , the x¯h values, h > 1, are obtained by adding K(h − 1) to x¯1 . 15

iii) By construction, starting from the initial value ηti1 t2 (Li ), the value of function ηti1 t2 decreases by 1 unit every time a new step value x¯h is reached, h ∈ {1, ..., |X|}. Theorem 3.1. The inequality t2 X t=t1

yit ≥ ηti1 t2 (Li ) − γ Ii(t1 −1) − Li




i ∈ N D , 1 ≤ t1 < t2 ≤ T,

(13)

where γ = 1/(¯ x1 − Li ), is valid for the 1-M-IRP-PD. Proof. We will prove that ηti1 t2 (Ii(t1 −1) ) ≥ ηti1 t2 (Li ) − γ Ii(t1 −1) − Li




(14)

from which inequality (13) follows. Let us denote by X0 = {Li } ∪ X the set obtained by adding the element x¯0 = Li to the set X. Note that the inventory level Ii(t1 −1) can be defined as Ii(t1 −1) = x¯h¯ +
θ, θ ≥ 0, where x¯h¯ , x¯h ∈ X0 , corresponds to the value such that ηti1 t2 Ii(t1 −1) = ηti1 t2 (¯ xh ). Let us consider the right-hand side of inequality (14). We have that: ηti1 t2 (Li ) −

Ii(t1 −1) − Li x¯¯ + θ − Li = ηti1 t2 (Li ) − h = (¯ x1 − Li ) (¯ x1 − Li )

= ηti1 t2 (Li ) −

¯ − 1) + θ (¯ x1 − Li ) + K(h . (¯ x1 − Li )

As (¯ x1 − Li ) ≤ K and θ ≥ 0, we have: ηti1 t2 (Li ) −

¯ − 1) + θ
(¯ x1 − Li ) + K(h ¯ = ηi ≤ ηti1 t2 (Li ) − h t1 t2 Ii(t1 −1) (¯ x1 − Li )

and the inequality holds.

16

Figure 1: Step function ηti1 t2 (Ii(t1 −1) ) for delivery nodes

In Figure 1 an example of the function ηti1 t2 and its linear approximation given by ηti1 t2 (Li )−γ(Ii(t1 −1) −Li ) is shown. Note that there are several linear approximations which can be used as lower bounds of function ηti1 t2 (Ii(t1 −1) ), namely all linear equations corresponding to straight lines lying below the step function. We chose to use inequality (14) because of the setting of the benchmark instances (described in Section 5), for which the interval inequalities are more effective for low values of xh . Finally, note that inequality (13) holds also for the classical IRP. 3.5.2. Interval inequalities for pickup nodes The step function ηti1 t2 (Ii(t1 −1) ) is adapted to pickup nodes as follows: &
' Dtt12 − Ui − Ii(t1 −1) i ηt1 t2 (Ii(t1 −1) ) = , (15) K where Dtt12 and K have the meaning as above. Proposition 3.2. The step function ηti1 t2 (Ii(t1 −1) ) for pickup nodes has the following properties: i) Ui − x¯1 = Dtt12 − K(ηti1 t2 (Ui ) − 1) ≤ K, ii) x¯h = x¯1 − K(h − 1) for each h ∈ {1, ..., |X|}, 17

Figure 2: Step Function ηti1 t2 (Ii(t1 −1) ) for pickup nodes

iii) ηti1 t2 (¯ xh ) = ηti1 t2 (Ui ) − h for each h ∈ {1, ..., |X|}. Theorem 3.2. The inequality t2 X t=t1

yi ≥ ηti1 t2 (Ui ) − γ(Ui − Ii(t1 −1) )

i ∈ N P , 1 ≤ t1 < t2 ≤ T,

(16)

where γ = 1/(Ui − x¯1 ), is valid for the 1-M-IRP-PD. For the sake of the brevity, we did not report the proof of the step function properties and of the resulting theorem, which can be easily derived from the one provided for the delivery nodes. In Figure 2 we show a representation of the step function ηti1 t2 (Ii(t1 −1) ) for pickup nodes and its linear approximation given by ηti1 t2 (Ui ) − γ(Ui − Ii(t1 −1) ). 4. Branch-and-cut algorithm The branch-and-cut algorithm adopted for the solution of the 1-M-IRPPD solves the proposed model strengthened through the valid inequalities discussed in Sections 3.1, 3.2, 3.3 and 3.5, which are all in polynomial number. More precisely, at each node of the branch-and-bound tree, we solve the linear relaxation of formulation (1a)-(1w), without subtour elimination constraints (1e), integrated by valid inequalities (2), (3) - (4), (6) - (7), (13) and (16). 18

Subtour elimination constraints and R-cut constraints, (1e) and (8) respectively, are added dynamically once violated. The separation of (1e) is performed using the max flow - min cut procedure described in [24]. Inequalities (8), as shown in Section 3.4, involve three different cuts, i.e. C(0, {i, j}), C(i, j), and C(j, i). Hence, their separation is obtained by repeating, for each of them, the same max flow-min cut procedure used for the subtour elimination constraints [11]. Both classes of constraints are added at each node of the branch-and-bound tree. 5. Computational results In this section the results of the computational experiments, aimed at assessing the effectiveness of the valid inequalities and the performance of the branch-and-cut algorithm, are presented. The experiments have been carried out on the benchmark instances proposed in [2] for the single vehicle case. The instances were generated by adapting IRP instances. For the sake of completeness, we recall the main features of the instances. • Routing costs. Instances provide the Euclidean coordinates of the depot and the customers. Costs cij are calculated as Euclidean distances. • Inventory costs. Two sets are considered: in the first, inventory costs are randomly generated in the interval [0.1, 0.5], while in the second in the interval [0.01, 0.05]. In the following, we will use ‘high inventory cost’ and ‘low inventory cost’ to denote the instances in the first and in the second set, respectively. • Time horizon. Two lengths for the time horizons are considered, T = 3, 6. • Number of customers. The number of customers is |N | = 5t, with t = 1, ..., 6 when T = 6 and t = 1, ..., 10 when T = 3. • Customer demands. The demands of all customers are constant over the planning horizon, i.e., dit = di , t ∈ T, i ∈ N . • Vehicle capacity. Four different capacity values are considered. Being Q the vehicle capacity in the original IRP instances, the four values of vehicle capacity for the 1-1-IRP-PD are obtained by considering dQ · αe, where α takes values: 1, 1/2, 1/4, 1/8. 19

For each combination of the parameters described above, five instances are generated by randomly changing customer demands. For the single vehicle case, we tested all the 640 generated instances (see [2] for more details). We adapted the single vehicle instances to the multiple vehicle case using the values dQ/me as vehicle capacity, where m is the number of vehicles varying between 2 and 5. However, we considered only single vehicle instances obtained with α = 1 and α = 1/2. The reason is that for α = 1/4, 1/8 many multiple vehicle instances turned out to be infeasible. For each combination of the parameters described above, five instances are generated by randomly changing customer demands. Thus, we obtained 1280 instances for the multiple vehicle case. In the following, the results obtained for the multiple vehicle case are described first. Then, the results obtained for the single vehicle case are presented, for which the proposed branch-and-cut algorithm is compared with the algorithm of [2]. The experiments have been performed on an Intel(R) Core(TM) i7-8700k, 3.70 GHz, 16.00 GB of RAM. The branch-and-cut algorithm has been coded in C language using Cplex 12.7 with default setting as MIP solver, with a maximum computing time of one hour. 5.1. Results for the multiple vehicle case We first show the results of the analysis of the effectiveness of the valid inequalities proposed in Sections 3.1–3.5. Then, we analyze the performance of the branch-and-cut algorithm. Table 2 shows the effect of the valid inequalities on the lower bound at the root node of the branch-and-bound tree, whereas Table 3 shows the performance of two different versions of the proposed branch-and-cut algorithm. The results are presented in different aggregated ways, grouping the instances on the basis of the different values of: α (vehicle capacity); T (horizon); inventory cost (high or low); m (number of vehicles; n (number of customers). The detailed results are available at the following url: http://www.ing.unisannio.it/boccia/IRP_PD.html. Table 2 reports the information about the percentage improvement of the lower bound at the root node obtained by adding the different classes of valid inequalities, in an additive way. All improvements are computed with respect to the linear programming relaxation of formulation (1), called LP from now on. More precisely, for each instance group (whose cardinality is reported in the second column) the percentage improvement of the following combination 20

of cuts and inequalities are given: LP plus the subtour elimination constraints (LP + SUB); LP + SUB plus the valid inequalities presented in Sections 3.1 and 3.5 (LP + SUB + VI); LP + SUB + VI plus the lot-sizing reformulation valid inequalities described in Sections 3.2 and 3.3 (LP + SUB + VI + LS); LP + SUB + VI + LS plus the R-cuts (LP + SUB + VI + LS + RCUTS); and LP + SUB + VI + LS + RCUTS plus Cplex presolving and cuts (LP + SUB + VI + LS + RCUTS + Cplex). All the considered variants provide an appreciable percentage improvement of the lower bound with the only exception of the variants which adds the R-cuts. The introduction of the subtour elimination constraints and of the lot-sizing inequalities provide the greatest improvement. When analyzing the impact of the parameters that define the instances, we notice that there is no remarkable impact apart for the case of instances with low inventory cost, for which the percentage improvement of each variant is generally much higher, and the lot-sizing inequalities are more effective. Even though the R-cut inequalities do not seem to be very effective from the results of Table 2, we decided to test them in the following set of experiments to analyze their impact when running the branch-and-cut algorithm up to termination. In Table 3, the results obtained by the two variants of branch-and-cut algorithm, with or without R-cut inequalities, are shown. For each group of instances we report the following information: number of instances in the set (# inst); number of instances for which a feasible solution has been found or the infeasibility is proved (# feas/infeas); number of instances for which the optimal solution was found or the infeasibility is proved (# opt/infeas); average percentage optimality gap (Av % opt. gap); maximum percentage optimality gap (Max % opt. gap); average CPU time overall the instances, in seconds (Av CPU time); average number of branchand-bound nodes (Av # B&B nodes). A feasible solution is found for 1264 out of 1280 instances and an optimal solution for 946 instances using Rcut inequalities and 1261 and 884, respectively, without R-cut inequalities. Hence, even though R-cut inequalities do not provide significant improvement of the lower bound at the root node, they allow to optimally solve a larger number of instances and contribute to a significant reduction of both the computational time and the average number of branch-and-bound nodes. Moreover, the average optimality gap is 0.87 with the R-cut inequalities and 1.16 without.

21

Table 2: Average % increase of the LP bound at the root node

α=1 α = 1/2 T=3 T=6 high low m=2 m=3 m=4 m=5 n=5 n = 10 n= 15 n= 20 n = 25 n = 30 n = 35 n = 40 n = 45 n = 50 All

# inst

LP +SUB

LP+SUB +VI

LP+SUB+VI +LS

LP+SUB+VI LS+RCUTS

640 640 800 480 640 640 320 320 320 320 160 160 160 160 160 160 80 80 80 80 1280

12.39 10.99 11.95 11.25 6.34 17.04 12.85 13.33 10.60 9.97 9.37 8.22 15.28 14.63 12.60 9.65 13.10 14.01 11.98 8.42 11.69

14.69 13.51 14.68 13.14 8.28 19.92 14.64 14.71 13.50 13.54 17.13 13.27 17.99 16.88 13.57 10.15 13.07 13.90 11.78 8.85 14.10

25.16 20.38 22.34 23.49 12.93 32.62 25.86 23.90 21.33 19.99 20.65 20.81 27.67 26.84 24.37 21.25 21.81 22.93 19.68 16.76 22.77

25.18 20.39 22.34 23.52 12.94 32.63 25.88 23.92 21.34 20.00 20.65 20.81 27.69 26.86 24.39 21.28 21.81 22.93 19.68 16.76 22.78

22

LP+SUB+VI LS+RCUTS +Cplex 25.18 23.11 24.23 25.87 14.18 35.51 27.06 25.67 23.76 22.88 22.37 25.14 30.11 29.60 26.57 22.72 22.96 23.60 20.41 17.51 24.84

Table 3: Results of the B&C algorithm with and without R-cut inequalities

23

α=1 α = 1/2 T=3 T=6 high low m=2 m=3 m=4 m=5 n=5 n = 10 n= 15 n= 20 n = 25 n = 30 n = 35 n = 40 n = 45 n = 50 All

# inst

# feas /infeas

640 640 800 480 640 640 320 320 320 320 160 160 160 160 160 160 80 80 80 80 1280

640 624 793 471 637 631 320 318 315 311 160 160 160 160 157 154 80 80 79 74 1264

B&C Algorithm - full version # opt Av % Max % Av /infeas opt. opt. CPU gap gap time 570 0.22 16.07 549.05 376 1.53 44.25 1673.30 647 0.59 21.61 856.08 299 1.34 44.25 1536.33 481 0.45 16.48 1070.03 465 1.29 44.25 1152.32 296 0.12 4.92 380.97 256 0.47 XX 890.32 210 1.25 43.85 1452.73 184 1.66 44.25 1720.67 160 0.00 0.00 2.51 137 0.26 5.30 716.87 131 0.29 3.95 839.12 102 1.22 16.23 1458.46 108 1.45 44.25 1407.41 95 1.73 43.85 1664.31 60 0.59 7.46 1155.55 53 0.53 5.15 929.42 51 1.46 21.61 1444.06 39 1.60 19.51 2072.39 946 0.87 44.25 1111.17

Av # B&B nodes

# feas /infeas

11527.50 28388.83 11070.26 34784.67 19228.42 20675.86 9078.99 17999.09 24953.45 27799.29 394.82 42259.50 27475.84 29256.14 21589.53 15529.93 14856.79 8354.34 9294.73 13463.33 19951.57

637 624 794 467 631 630 320 318 314 309 160 160 160 160 154 151 80 80 80 74 1261

B&C Algorithm - No R-cuts # opt Av % Max % Av /infeas opt. opt. CPU gap gap time 507 0.76 29.54 895.53 377 1.57 49.49 1676.78 611 0.80 19.25 1039.28 273 1.77 49.49 1697.61 451 0.62 49.49 1237.52 433 1.70 32.86 1334.80 297 0.12 4.73 379.96 255 0.54 18.27 887.67 144 2.41 29.54 2175.41 188 1.60 49.49 1701.58 160 0.00 0.00 1.93 132 0.29 4.51 793.35 127 0.39 4.06 950.11 96 1.48 19.27 1592.61 94 2.14 29.54 1655.37 85 2.19 49.49 1919.18 54 0.72 9.39 1461.99 57 0.79 6.02 1219.57 45 1.95 12.59 1766.04 34 2.36 19.29 2305.79 884 1.16 49.49 1286.16

Av # B&B nodes 16262.03 28571.68 13586.99 37133.29 21374.31 23459.40 8501.14 17546.84 34412.67 29206.76 305.03 46694.53 30998.54 30566.48 24721.61 18140.97 17761.66 11439.85 12377.61 14236.23 22416.85

5.2. Results for the single vehicle case In Tables 4 and 5 the results obtained on the 640 single vehicle instances proposed in [2] are shown. The instances are grouped in different ways, according to the different values of: α (vehicle capacity); T (horizon); inventory cost (high or low); n (number of customers). In the tables information is reported about the performance of the proposed branch-and-cut algorithm (referred to as B&C) and of the branch-and-cut algorithm presented in [2] (referred to as ACS). The information about the ACS is directly taken from the paper, as the workstations on which the experiments have been performed are comparable. For the B&C algorithm, we used the version tested for the multiple vehicle case which provided the best performance, i.e., the one with all valid inequalities and the R-cuts inequalities. In Table 4 the performances of the two branch-and-cut algorithms are compared in terms of solution quality. For each group of instances the following information is reported: number of instances in the set (# inst); number of instances for which a feasible solution is found or the infeasibility is proved (# feas/infeas); number of instances for which the optimal solution was found or the infeasibility is proved (# opt/infeas); average percentage optimality gap (Av % opt. gap); maximum percentage optimality gap (Max % opt. gap); number of instances for which an algorithm provides better solutions than the other (column ‘> ACS’ when B&C is better than ACS, and ‘> B&C’ in the opposite case). In Table 5 the performance of the two algorithms is compared considering the computational aspects, i.e. the average CPU time (Av CPU time) and average number of branch-and-bound nodes (Av # B&B nodes). From Table 4 it can be seen that the B&C algorithm finds a feasible solution for 627 out of 640 instances and an optimal solution for 538 instances, against 612 feasible and 473 optimal solutions of the ACS algorithm. The average optimality gap of the B&C algorithm is 0.89, which improves the 1.22 gap of the ACS algorithm. Moreover, B&C improves the solution obtained by ACS on 133 instances. The optimality gaps are higher for the instances with very low capacity (α = 1/8), longer horizon and low inventory. In particular, the instances with α = 1/8 are characterized by a narrow solution space. The computational performance of the two algorithms (see Table 5) is similar, with a slight improvement for the B&C algorithm computational time, despite the higher number of explored nodes.

24

Table 4: Comparison of B&C and ACS: solution quality # inst

25

α=1 α = 1/2 α = 1/4 α = 1/8 T=3 T=6 high low n=5 n = 10 n = 15 n = 20 n = 25 n = 30 n = 35 n = 40 n = 45 n = 50 All

160 160 160 160 400 240 320 320 80 80 80 80 80 80 40 40 40 40 640

# feas /infeas 160 160 160 147 396 231 316 311 80 80 80 80 74 77 40 40 39 37 627

# opt /infeas 160 160 142 76 351 187 273 265 80 75 72 63 61 60 34 34 29 30 538

B&C Av % opt. gap 0.00 0.00 0.29 3.49 0.62 1.35 0.53 1.26 0.00 0.20 0.29 1.69 1.07 1.76 0.61 0.66 1.60 1.55 0.89

Max % opt. gap 0.00 0.00 5.09 59.10 29.31 59.10 21.88 59.10 0.00 4.63 5.51 47.06 21.88 59.10 8.80 8.36 23.61 29.31 59.10

> ACS 8 8 43 74 56 77 66 67 0 9 15 19 16 32 7 8 10 17 133

# feas /infeas 158 159 158 137 392 220 306 306 80 80 80 77 76 67 40 39 38 35 612

# opt /infeas 152 152 111 58 331 142 235 238 80 68 64 57 52 40 31 32 27 22 473

ACS Av % opt. gap 0.08 0.08 0.73 4.45 1.12 1.42 0.75 1.69 0.00 0.39 0.50 1.13 1.65 1.95 0.96 1.14 4.94 1.81 1.22

Max % opt. gap 5.23 5.23 7.87 56.91 56.91 30.39 47.11 56.91 0.00 5.23 5.94 12.11 21.22 30.39 10.78 17.27 56.91 19.52 56.91

> B&C 0 0 3 14 7 10 9 8 0 1 1 4 6 3 1 0 1 0 17

Table 5: Comparison of B&C and ACS: computational time and nodes α=1 α = 1/2 α = 1/4 α = 1/8 T=3 T=6 high low n=5 n = 10 n = 15 n = 20 n = 25 n = 30 n = 35 n = 40 n = 45 n = 50 All

# inst 160 160 160 160 400 240 320 320 80 80 80 80 80 80 40 40 40 40 640

B&C Av CPU time Av # of nodes 31.56 147.18 27.73 310.72 885.54 12300.69 2197.71 23258.59 610.89 4897.7 1078.21 15848.62 757.37 8635.26 813.89 9373.34 3.15 277.89 479.28 18520.13 637.68 14104.65 860.71 10800.9 986.45 9264.39 1237.93 7823.69 860.52 7790.38 702.87 4108.35 1183.26 4012.25 1413.12 6574.5 785.63 9004.30

ACS Av CPU time Av # of nodes 281.99 321.11 326.24 667.21 1263.09 7702.24 2427.28 14224.24 744.35 2838.79 1625.16 10545.21 1074.15 5354 1075.1563 6103.39 1.93 524.59 599.28 20949.03 753.61 7392.35 1304.45 7002.66 1384.86 4349.46 2023.51 2749.46 952.6 2125.9 1005.9 1439.7 1341.13 1137.55 1759.53 1020.58 1074.65 5728.70

6. Conclusions In this paper we have studied the inventory routing problem where pickup and delivery customers have to be served with multiple vehicles. A mixed integer linear programming formulation is presented that uses aggregated variables to represent the amount of commodity traversing each arc of the graph. Moreover, different classes of valid inequalities are proposed, some of which adapted from the literature and one new class, the interval inequalities. These latter inequalities contribute to improving the lower bound at the root node. Computational results on benchmark multiple vehicle instances show that the proposed branch-and-cut algorithm optimally solves, in one hour of computational time, 946 out of 1280 instances, and it finds feasible solutions for 1264 instances with an average optimality gap of 0.87%. On the single vehicle instances, the algorithm was shown to outperform the one proposed in [2] with 538 out of 640 instances solved to optimality, against the previous 473, and an average optimality gap of 0.89%, against the previous 1.22%. Moreover, the proposed algorithm finds 133 new best known solutions. Future research directions include investigating the impact on the interval inequalities on other problems, for example on the classical IRP. Heuristics 26

for larger instances should also be designed. The problem studied in this paper is the first and probably most essential routing problem with pickup and delivery customers to be served with a fleet of vehicles over a period of time. Several variants of this problem would find application in real-life applications. Acknowledgements The research activity of the third author was partially supported by the Italian Ministry of Education, National Research Program (PRIN) 2015, contract no. 20153TXRX9. The research activity of the forth and fifth author was partially funded by the University “Federico II” of Naples within the research project MOSTOLOG project (A multi-objective approach for Sustainable Logistic Systems). References References [1] A. Agra, M. Christiansen, and A. Delgado. Mixed integer formulations for a short sea fuel oil distribution problem. Transportation Science, 47:108–124, 2013. [2] C. Archetti, M. Christiansen, and M. G. Speranza. Inventory routing with pickups and deliveries. European Journal of Operations Research, 268:314–324, 2018. [3] P. Avella, M. Boccia, and I. Vasilyev. Lifted and local reachability cuts for the vehicle routing problem with time windows. Computers & Operations Research, 40:2004–2010, 2013. [4] P. Avella, M. Boccia, and L. A. Wolsey. Single-item reformulations for a vendor managed inventory routing problem: computational experience with benchmark instances. Networks, 64:129–138, 2014. [5] P. Avella, M. Boccia, and L. A. Wolsey. Single-period cutting planes for inventory routing problems. Transportation Science, 52:497–508, 2018.

27

[6] J. F. Bard, L. Huang, P. Jaillet, and M. Dror. A decomposition approach to the inventory routing problem with satellite facilities. Transportation Science, 32:189–203, 1998. [7] M. Battarra, J.-F. Cordeau, and M. Iori. Pickup-and-delivery problems for goods transportation. vehicle routing: Problems, methods, and applications. MOS/SIAM series on optimization, 2:161–192, 2014. [8] G. Berbeglia, J.-F. Cordeau, I. Gribkovskaia, and G. Laporte. Static pickup and delivery problems: A classification scheme and survey. TOP, 15:1–31, 2007. [9] L. Bertazzi and M. G. Speranza. Inventory routing problems: an introduction. EURO Journal on Transportation and Logistics, 1:307–326, 2012. [10] L. Bertazzi and M. G. Speranza. Inventory routing problemswith multiple customers. EURO Journal on Transportation and Logistics, 2:255– 275, 2012. [11] M. Boccia, T. G. Crainic, A. Sforza, and C. Sterle. Multi-commodity location-routing: Flow intercepting formulation and branch-and-cut algorithm. Computers & Operations Research, 89:94–112, 2018. [12] M. Christiansen and K. Fagerholt. Ship routing and scheduling in industrial and tramp shipping. Vehicle routing: problems, methods, and applications. MOS/SIAM series on optimization, 2e, pages 381–408, 2014. [13] M. Christiansen, K. Fagerholt, B. Nygreen, and D. Ronen. Ship routing and scheduling in the new millennium. European Journal of Operational Research, 228:467–483, 2013. [14] L. C. Coelho, J.-F. Cordeau, and G. Laporte. Thirty years of inventory routing. Transportation Science, 48:1–19, 2013. [15] N. C. P. Edirisinghe and R. J. W. James. Fleet routing position-based model for inventory pickup under production shutdown. European Journal of Operational Research, 236:736–747, 2014. [16] F. G. Engineer, K. C. Furman, G. L. Nemhauser, M.W.P. Savelsbergh, and J.-H. Song. A branch-price-and cut algorithm for single-product maritime inventory routing. Operations Research, 60:106–122, 2012. 28

[17] R. Grønhaug, M. Christiansen, G. Desaulniers, and J. Desrosiers. A branch-and-price method for a liquefied natural gas inventory routing problem. Transportation Science, 44:400–415, 2010. [18] H. Hern´andez-P´erez and J. J. Salazar-Gonz´alez. A branch-and-cut algorithm for a traveling salesman problem with pickup and delivery. Discrete Applied Mathematics, 145:126–139, 2004. [19] G. Iassinovskaia, S. Limbourg, and F. Riane. The inventory-routing problem of returnable transport items with time windows and simultaneous pickup and delivery in closed-loop supply chains. International Journal of Production Economics, 183:570–582, 2017. [20] C. Li, L. Gong, Z. Luo, and A. Lim. A branch-and-price-and-cut algorithm for a pickup and delivery problem in retailing. Omega, 89:71 – 91, 2019. [21] J. Lysgaard. Reachability cuts for the vehicle routing problem with time windows. European Journal of Operational Research, 175:210–223, 2006. [22] A. Nikolakopoulos. An economic model predictive approach for inventory routing and control with time windows constraints: Application in the distribution of industrial gases. Chemical Engineering, 52:925–930, 2016. [23] A. Nikolakopoulos and I. Ganas. Economic model predictive inventory routing and control. Central European Journal of Operations Research, 25:587–609, 2017. [24] M. Padberg and G. Rinaldi. A branch-and-cut algorithm for the resolution of large-scale symmetric traveling salesman problems. SIAM Rev., 33:60–100, 1991. [25] N. Ramkumar, P. Subramanian, T. T. Narendran, and K. Ganesh. Mixed integer linear programming model for multi-commodity multidepot inventory routing problem. OPSEARCH, 49:413 – 429, 2012. [26] M. Savelsbergh and J.-H. Song. Inventory routing with continuous moves. Computers & Operations Research, 34:1744–1763, 2007.

29

[27] M. Savelsbergh and J.-H. Song. An optimization algorithm for the inventory routing problem with continuous moves. Computers & Operations Research, 35:2266–2282, 2008. [28] R. G. van Anholt, L. C. Coelho, G. Laporte, and I. F. A. Vis. An inventory-routing problem with pickups and deliveries arising in the replenishment of automated teller machines. Transportation Science, 50:1077–1091, 2016. [29] Y. Wang, J. Zhang, K. Assogba, Y. Liu, M. Xu, and Y. Wang. Collaboration and transportation resource sharing in multiple centers vehicle routing optimization with delivery and pickup. Knowledge-Based Systems, 160:296–310, 2018.

30