A multi-phase heuristic for the production routing problem

Accepted Manuscript

A multi-phase heuristic for the production routing problem Oguz ˘ Solyalı, Haldun Sural ¨ PII: DOI: Reference:

S0305-0548(17)30144-2 10.1016/j.cor.2017.06.007 CAOR 4261

To appear in:

Computers and Operations Research

Received date: Revised date: Accepted date:

14 August 2015 18 May 2017 1 June 2017

Please cite this article as: Oguz A multi-phase heuristic for the production routing ˘ Solyalı, Haldun Sural, ¨ problem, Computers and Operations Research (2017), doi: 10.1016/j.cor.2017.06.007

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 • The production routing problem is considered. • A multi-phase heuristic is proposed for the problem.

CR IP T

• Computational results show that the heuristic is very effective.

AC

CE

PT

ED

M

AN US

• The heuristic found new best solutions for the majority of the benchmark instances.

1

ACCEPTED MANUSCRIPT

A multi-phase heuristic for the production routing problem

Business Administration Program, Middle East Technical University, Northern Cyprus Campus, Kalkanli 99738, Mersin 10, Turkey, [email protected] b

Department of Industrial Engineering, Middle East Technical University, Ankara 06800, Turkey, [email protected]

AN US

June 7, 2017

Abstract

This study considers the production routing problem where a plant produces and distributes a single item to multiple retailers over a multi-period time horizon. The problem is to decide on when and how much to

M

produce and stock at the plant, when and how much to serve and stock at each retailer, and vehicle routes for shipments such that the sum of fixed production setup cost, variable production cost, distribution cost, and

ED

inventory carrying cost at the plant and retailers is minimized. A multi-phase heuristic is proposed for the problem. The proposed heuristic is a mathematical programming-based heuristic that relies on formulating and solving restricted versions of the problem as mixed integer programs. The computational experiments on

PT

benchmark instances show favorable results with regard to the quality of the solutions found at the expense of higher computing times on large instances. In particular, the heuristic managed to find new best solutions

CE

for the 65% of benchmark instances. Keywords: Supply chain, production routing, integer programming, heuristics

AC

a

CR IP T

O˘guz Solyalıa , Haldun S¨ uralb

1

Introduction

The integration of production, inventory and distribution routing decisions in supply chains offers tremendous cost saving opportunities to firms. In order to realize these savings, one has to solve a production routing

2

ACCEPTED MANUSCRIPT

problem (PRP), in which a plant produces single or multiple products and distributes them to multiple retailers over a multi-period time horizon by jointly considering all these decisions. The PRP typically involves production decisions at a plant, replenishment decisions at multiple retailers and routing decisions for vehicles; therefore, it can be seen as an integration of two-level lot sizing and vehicle routing problems, both of which are quite difficult to solve. Because the PRP is a very complex problem, there are only a few exact algorithms (see Adulyasak et al., 2014b; Archetti et al., 2011; Ruokokoski et al., 2010)

CR IP T

available in the literature and these algorithms can only solve small-scale instances to optimality. Therefore, many heuristics, e.g. an adaptive large neighborhood search (Adulyasak et al., 2014a), branch-and-price heuristics (Bard and Nananukul, 2009a, 2010), decomposition heuristics (Bertazzi et al., 2005; Boudia et al., 2008; Chandra, 1993; Chandra and Fisher, 1994; Lei et al., 2006), a GRASP (Boudia et al., 2007), Lagrangian heuristics (Fumero and Vercellis, 1999; Solyalı and S¨ ural, 2009), matheuristics (Absi et al., 2015; Archetti

AN US

et al., 2011), a memetic algorithm (Boudia and Prins, 2009), and tabu search algorithms (Armentano et al., 2011; Bard and Nananukul, 2009b) have been proposed for several variants of the PRP. The majority of these heuristics have been tested on the benchmark instances generated by Boudia et al. (2007) and Archetti et al. (2011) with the most successful ones being those proposed by Absi et al. (2015), Adulyasak et al. (2014a), and Archetti et al. (2011). For a recent review on the PRP, interested readers can refer to Adulyasak et al.

M

(2015). The PRP is also very related to the inventory routing problem (IRP) in that the PRP reduces to the IRP when there are no production decisions at the plant/supplier level (i.e., production amounts are assumed to be known or unlimited in IRP). See Coelho et al. (2014) for a recent review on the IRP.

ED

In this study, we introduce a multi-phase heuristic for the PRP. The first two phases of the proposed heuristic is based on the a priori tour idea, which first appeared in Pınar and S¨ ural (2006) and then in Solyalı

PT

and S¨ ural (2011) for solving the single-vehicle inventory routing problems. We here extend the a priori tour idea to multiple vehicles. The idea is to impose a fixed sequence of all retailers for all routes in which a vehicle follows the given sequence for the retailers to be visited and skips the non-visited ones. We obtain

CE

the fixed sequence of all retailers by solving a traveling salesman problem (TSP) in the first phase. In the second phase, using the obtained sequence, we formulate and solve a restricted version of the PRP, i.e., the

AC

PRP with the imposed a priori tour, as a mixed integer program (MIP) where the a priori tour is modelled as a multiple TSP with aggregate vehicle capacity constraints. Because individual vehicle capacities may be violated in the solution of the second phase, we solve a capacitated vehicle routing problem (CVRP) or a TSP for each period depending on the number of vehicles required in that period to obtain the vehicle routes

3

ACCEPTED MANUSCRIPT

in the third phase. Note that the routes obtained at the end of third phase may still be infeasible in the case of a limited number of vehicles and not allowing split deliveries. In the fourth phase, we try to obtain a feasible and improved solution exploiting an idea proposed in Archetti et al. (2012) and used in Coelho and Laporte (2013) for IRPs. Specifically, we formulate and solve a MIP for a restricted version of the PRP in which we impose the vehicle routes obtained in the third phase to the MIP and allow the insertion and

order to improve the solution obtained in the fourth phase.

CR IP T

removal of retailers to these vehicle routes. In the fifth phase, we solve a TSP for each period and vehicle in

Existing heuristics that involve solving a MIP like those of Absi et al. (2015), Adulyasak et al. (2014a), and Bard and Nananukul (2009b) consider approximating the distribution cost in the MIP by defining a visiting cost for retailers and associated variables. Unlike them, the proposed heuristic uses real distribution costs instead of approximating them by explicitly modeling routing decisions in the MIP in the second phase

AN US

at the expense of additional variables.

The proposed heuristic has been tested on the benchmark instances generated by Boudia et al. (2007) and Archetti et al. (2011) and compared with the most successful heuristics mentioned above. The computational experiments on benchmark instances show favorable results with regard to the quality of the solutions found at the expense of higher computing times on large instances. In particular, the heuristic managed to find new

M

best solutions for 60 of 90 instances generated by Boudia et al. (2007) and 623 of 960 instances generated by Archetti et al. (2011).

The rest of the manuscript is organized as follows. Section 2 describes the problem and provides a

ED

standard MIP formulation of the PRP considered in this paper. The multi-phase heuristic is presented in Section 3. Section 4 presents the computational results on benchmark instances. Finally, Section 5 concludes

Problem Formulation

CE

2

PT

the paper.

We consider a production-distribution system in which a plant denoted by 0 produces a single product and

AC

distributes to n retailers over a discrete time horizon with a length of T periods. Each retailer i faces external customer demand dit in period t and may keep inventory Iit to meet the demand without backordering. The plant decides on how much to produce in each period without exceeding its production capacity C, distributes to the retailers immediately, and keeps inventory for replenishing retailers in later periods. The plant also manages the inventories at the retailers by deciding on when and how much to ship to each retailer and

4

ACCEPTED MANUSCRIPT

guarantees that neither retailers nor itself will stock-out in any period. When a production occurs at the plant in a period, a fixed setup cost f , independent of the size of production, and a variable production cost b per unit produced are incurred. For each unit kept in inventory at a retailer and plant, a holding cost hi is incurred at the end of a period. The amount of inventory kept at a retailer (plant) cannot exceed the storage capacity Ui (U0 ) of retailer i (plant). A fleet of m homogeneous vehicles based at the plant is available for distribution of the product to retailers. Each vehicle, departing from and returning back to the

CR IP T

plant, can visit several retailers in a route without exceeding its capacity Q. A vehicle traveling from facility i to facility j incurs a distribution cost cij . We assume that a vehicle can only perform a single tour in every period and the splitted deliveries to retailers are not allowed. The PRP is to decide on when and how much to produce and stock at the plant, when and how much to serve and stock at each retailer, and vehicle routes for shipments such that the sum of fixed production setup cost, variable production cost, distribution cost,

AN US

and inventory carrying cost at the plant and retailers is minimized. Let pt be the quantity produced at the plant in period t, qit the quantity shipped to retailer i in period t, Iit (I0t ) the inventory level of retailer i (plant) at the end of period t, yt 1 if a production occurs at the plant in period t and 0 otherwise, z0t the number of vehicles dispatched from the plant in period t, zit 1 if retailer i is visited in period t and 0 otherwise, and xijt 1 if node j is visited immediately after node i by a vehicle in period t and 0 otherwise.

F : min

M

Then, we can formulate the PRP as follows:

T X n n n X X X [ hi Iit + f yt + bpt + cij xijt ] t=1 i=0

n X

ED

s.t. I0,t−1 + pt =

qit + I0t

i=1

PT

Ii,t−1 + qit = dit + Iit

CE

pt ≤ min{C,

n X T X i=1 k=t

dik }yt

AC

Iit ≤ Ui

qit ≤ min{Ui + dit , Q, n X i=1

(1)

i=0 j=0,j6=i

T X k=t

dik }zit

qit ≤ λQz0t

z0t ≤ m

5

1≤t≤T

(2)

1 ≤ i ≤ n, 1 ≤ t ≤ T

(3)

1≤t≤T

(4)

0 ≤ i ≤ n, 1 ≤ t ≤ T

(5)

1 ≤ i ≤ n, 1 ≤ t ≤ T

(6)

1≤t≤T

(7)

1≤t≤T

(8)

ACCEPTED MANUSCRIPT

n X

xijt = zit

0 ≤ i ≤ n, 1 ≤ t ≤ T

(9)

xjit = zit

0 ≤ i ≤ n, 1 ≤ t ≤ T

(10)

S ⊆ {1, ..., n} : |S| ≥ 2, 1 ≤ t ≤ T

(11)

xijt ∈ {0, 1}

0 ≤ i ≤ n, 0 ≤ j ≤ n, i 6= j, 1 ≤ t ≤ T

(12)

yt ∈ {0, 1}

1≤t≤T

(13)

zit ∈ {0, 1}

1 ≤ i ≤ n, 1 ≤ t ≤ T

(14)

z0t ∈ Z+

1≤t≤T

(15)

pt , qit , I0t , Iit ≥ 0

1 ≤ i ≤ n, 1 ≤ t ≤ T.

(16)

j=0,j6=i n X

j=0,j6=i

X i∈S

zit −

X

qit /Q

i∈S

AN US

i∈S j∈S,i6=j

xijt ≤

CR IP T

X X

The objective function (1) in F minimizes the sum of inventory holding costs at the retailers and plant, fixed and variable production costs at the plant and distribution costs. Constraints (2) and (3) are inventory balance equations for the plant and retailers, respectively. Constraints (4) ensure that the total production at the plant cannot exceed the available capacity and the production is possible only if the necessary fixed cost is paid. Constraints (5) are for restricting the storage capacity while constraints (6) are for bounding the

M

shipment quantity to retailers. Constraints (6) also ensure that a retailer is visited if that retailer receives some quantity in a period. Constraints (7) and (8) limit the total shipment quantity to retailers by the

ED

fleet capacity. Note that we introduce the parameter λ in constraints (7). When the PRP is solved to optimality, one must set λ to 1. However, in the heuristic we propose in Section 3, the parameter λ < 1

PT

reduces the available vehicle capacity in the second phase of the multi-phase heuristic in order to facilitate finding a feasible solution at the third phase and acts as a factor to remove the impact of no split delivery. Constraints (9) and (10) are degree constraints and impose that each retailer is visited once in a period if

CE

necessary and that a particular vehicle enters and leaves a particular retailer. Constraints (9) and (10) also impose that the number of vehicles dispatched in a period enters and leaves the plant. Constraints (11) are

AC

the generalized fractional subtour elimination constraints (see Adulyasak et al., 2014b), which do not only prevent subtours but also ensure that individual vehicle capacities are not exceeded. Constraints (12)–(14) are for binary restrictions, (15) are for integrality restrictions, and (16) are for nonnegativity of variables.

6

ACCEPTED MANUSCRIPT

3

The Multi-Phase Heuristic

The multi-phase heuristic decomposes the PRP into five phases. We solve a giant TSP over all retailers and plant in the first phase. The second phase decides on the production quantities at the plant, the retailers to be visited in each period and the delivery quantities to these retailers by formulating and solving a restricted version of the problem as a MIP in which the giant tour obtained in the first phase is imposed to the MIP.

CR IP T

Then, we solve a CVRP or a TSP depending on the number of vehicles needed in each period in the third phase so as to obtain the vehicle routes. In the fourth phase, we improve the solution obtained in the earlier phases by formulating and solving a restricted version of the problem as a MIP, in which the vehicle routes obtained in the third phase are imposed and the insertion and removal of retailers to these routes are allowed. This phase decides on the production and delivery quantities in each period. Finally, we solve a TSP for each period and vehicle in the last phase in order to obtain the best vehicle routes for the deliveries decided

AN US

in the fourth phase.

Note that the F formulation in Section 2 is defined on a directed graph in order to facilitate the introduction of the MIP in the second phase. However, our multi-phase heuristic and TSP and CVRP tools we use are all assuming an undirected underlying graph.

3.1

The First Phase

M

Next we describe the five phases of our heuristic in detail.

ED

The first phase is to solve the TSP over all nodes which yields a giant tour. The giant tour specifies the sequence of retailers to be visited in a period. In order to store this sequence, once the giant TSP is solved, we form two sets αi and βi for each node i (0 ≤ i ≤ n). βi (αi ) denotes the set of nodes that can be visited

PT

before (after) node i according to the sequence of node i on the giant tour. Note that α0 and β0 contain all retailers (i.e., α0 = {1, ..., n} and β0 = {1, ..., n}). Also, note that αi and βi involve plant 0 for each retailer

CE

i where 1 ≤ i ≤ n. In order to clarify the sets αi and βi , we present an example with a three-retailer PRP instance in Table 1.

The Second Phase

AC

3.2

The second phase is to solve the PRP on a restricted tour that is the same as the one defined in Section 2 with the exception of imposing the fixed sequence of retailers to visit according to the sequence specified by the a priori tour obtained in the first phase. We formulate and solve the following MIP, F 0 , for the second

7

ACCEPTED MANUSCRIPT

Table 1: All αi and βi values for a 3-retailer PRP instance on a given giant tour 0–2–3–1–0 αi {1,2,3} {0} {0,1,3} {0,1}

βi {1,2,3} {0,2,3} {0} {0,2}

phase decisions: T X n n X X X [ hi Iit + f yt + bpt + cij xijt ] t=1 i=0

i=0 j∈αi

s.t. (2) − (8), (13) − (16), X xijt = zit j∈αi

X

xjit = zit

j∈βi

0 ≤ xijt ≤ 1

(17)

0 ≤ i ≤ n, 1 ≤ t ≤ T

(18)

0 ≤ i ≤ n, 1 ≤ t ≤ T

(19)

0 ≤ i ≤ n, j ∈ αi , 1 ≤ t ≤ T.

(20)

AN US

F 0 : min

CR IP T

i 0 1 2 3

M

Besides a slight change in the last term of the objective function, the difference of F 0 from F is that constraints (18)–(20) replaced constraints (9)–(12) of F . The constraints (18)–(19) impose the a priori tour

ED

obtained in the first phase to the formulation so that the retailers are visited in the order dictated by the a priori tour by skipping the non-visited retailers. For instance, let us consider the example in Table 1 and assume that retailers 1 and 2 are to be visited in period t with q1t + q2t < C. In this case, the a priori tour

PT

would impose a single vehicle to follow the route 0–2–1–0. There are two important computational advantages of using F 0 instead of F in solving PRP. The first one

CE

is that we do not need to tackle with the subtour elimination constraints (11) whose numbers are increasing exponentially. The second one is that we do not need to deal with the difficult decision of routing vehicles because the routes will be automatically derived given the optimal zit values. Note that integral xijt values

AC

are obtained when F 0 is solved because given the integral zit values, the constraint matrix regarding xijt ,

i.e., constraints (18)–(20), corresponds a transportation problem for each period which is well-known to be totally unimodular. On the other hand, due to only having the aggregate vehicle capacity constraints (7), individual vehicle capacities may be exceeded and the obtained solution may be infeasible. Note that the

8

ACCEPTED MANUSCRIPT

solution yielded by F 0 is feasible in the case of a single vehicle. At the end of the second phase, we decide on when and how much to produce at the plant as well as when and how much to ship to retailers. Note that in the presence of multiple vehicles, we model the a priori tour as a multiple TSP in F 0 without using an additional index to denote different vehicles in x and z variables. Indeed, using an additional index for different vehicles is straightforward to do and it enables to easily enforce the feasibility of individual vehicle capacities but the resulting formulation is computationally very demanding and we could not even

3.3

CR IP T

solve the resulting formulation for the relatively small-size instances (i.e., when n = 50) in reasonable times.

The Third Phase

Given the solution of the second phase, we decide on the routes of the vehicles in the third phase. Specifically, using the values of qit and z0t variables in the solution of F 0 , we solve a CVRP if z0t > 1 or a TSP if z0t = 1

AN US

in each period t. In the CVRPs the demands of retailers are equal to the qit values while the retailers to be visited in a period are those with qit > 0 in the TSPs. Note that the solution obtained in the second and third phases may be infeasible with regard to the number of vehicles needed when the number of vehicles is limited and the split delivery is not allowed. Due to the aggregate capacity constraints (7) and not allowing

3.4

The Fourth Phase

M

split delivery, the number of vehicles needed can be at most m + 1.

ED

In order to obtain a feasible and improved solution to the PRP, in the fourth phase we formulate a restricted version of the problem as a MIP in which we dictate the routes found in the third phase to the PRP but allow insertion and removal of retailers in these routes. Specifically, we define sitk (resp. ritk ) as a binary

PT

variable to account for the insertion (resp. removal) of retailer i into (resp. from) route k in period t with the associated cost Γitk (resp. ∆itk ). Based on the solution of third phase, we let aitk be a parameter being

CE

equal to 1 if retailer i is visited in period t by vehicle k, and 0 otherwise. Note that if the solution of the third phase involves usage of a fleet size greater than m in a period, then we take only m routes and disregard the route with the least number of visited retailers in that period. As a result, we formulate and solve the

AC

following MIP, IF , in the fourth phase:

IF : min

T X n n X m X X [ hi Iit + f yt + bpt + ((1 − aitk )Γitk sitk − aitk ∆itk ritk )] t=1 i=0

i=1 k=1

9

(21)

ACCEPTED MANUSCRIPT

s.t. (4), (5), (13), n X m X

0 qitk + I0t

(22)

1 ≤ i ≤ n, 1 ≤ t ≤ T

(23)

1 ≤ i ≤ n, 1 ≤ t ≤ T, 1 ≤ k ≤ m

(24)

0 qitk ≤Q

1 ≤ t ≤ T, 1 ≤ k ≤ m

(25)

(aitk − aitk ritk + (1 − aitk )sitk ) ≤ 1

1 ≤ i ≤ n, 1 ≤ t ≤ T

(26)

(aitk ritk + (1 − aitk )sitk ) ≤ L

1 ≤ t ≤ T, 1 ≤ k ≤ m

(27)

1 ≤ i ≤ n, 1 ≤ t ≤ T, 1 ≤ k ≤ m

(28)

i=1 k=1

Ii,t−1 +

m X

0 qitk = dit + Iit

k=1 0 qitk ≤ Mit (aitk − aitk ritk + (1 − aitk )sitk ) n X

i=1 m X

k=1 n X i=1

ritk ∈ {0, 1}, sitk ∈ {0, 1}

AN US

0 pt , qitk , I0t , Iit ≥ 0

CR IP T

1≤t≤T

I0,t−1 + pt =

1 ≤ i ≤ n, 1 ≤ t ≤ T, 1 ≤ k ≤ m,

(29)

0 where qitk be the quantity shipped to retailer i in period t by vehicle k, L is the number of allowable insertion PT and removal moves, and Mit = min{Ui + dit , Q, k=t dik }.

The objective function (21) in IF is the sum of inventory holding costs at plant and retailers, fixed and

M

variable production costs, and the costs of inserting and removing retailers to/from vehicle routes. Equations (22) and (23) are the inventory balance equations for the plant and retailers, respectively. Constraints (24)

ED

ensure that a retailer i can only be shipped a positive quantity in a period t by a vehicle k if that retailer is visited in t by k. Constraints (25) guarantee that individual vehicle capacities are not exceeded. Constraints (26) do not allow split deliveries by ensuring that each retailer can be visited by at most one vehicle in a

PT

period. Constraints (27) limit the total number of insertions and removals by L for each period and each route. Note that setting a larger L would be better for eliminating infeasible routes that may be found at

CE

the end of the third phase but the computed insertion and removal costs would have a higher chance to be inaccurate. For example, if only one insertion or removal occurs in a route, the change in the distribution

AC

cost of the route is accurately computed by s and r variables, whereas if two consecutive retailers’ insertion or removal is concerned, the change in the distribution cost of the route may not match with the one computed by s and r variables. Constraints (28) are for nonnegativity of variables and (29) are for binary restrictions.

10

ACCEPTED MANUSCRIPT

3.5

The Fifth Phase

0 Using the values of qitk variables found in the fourth phase, the routing of vehicles are determined in the 0 fifth phase by solving a TSP for each vehicle k in each period t. Note that only the retailers with qitk >0

are to be visited by vehicle k in period t.

Computational Results

CR IP T

4

We have performed our computational experiments on the benchmark instances introduced by Archetti et al. (2011) and Boudia et al. (2007) to assess the computational performance of the multi-phase heuristic (abbreviated as 5P). We have solved the mathematical formulations in the second and fourth phases (i.e., F 0 and IF ) using Cplex 12.5 with its default settings and a single thread, the CVRPs in the third phase

AN US

using the heuristic of Groer et al. (2010), and the TSPs in the first, third and fifth phases using Concorde (Applegate et al., 2007). We have performed the experiments on a 2.40 GHz Workstation with 4x12 GB RAM and 12 cores operating on Windows 7 (64-bit).

The test instances as generated by Archetti et al. (2011) involve six periods and 14 (A1), 50 (A2), and 100 (A3) retailers in three sets, each of which has 480 test instances. Being larger scale instances than those

M

of Archetti et al. (2011), the instances generated by Boudia et al. (2007) involve 20 periods and respectively 50, 100, and 200 retailers in three sets of instances (B1, B2, and B3) with each set having 30 test instances. A summary of the benchmark instances are presented in Table 2. We have compared our results with the

ED

results of the heuristics (IM-VRP and IM-MultiTSP) proposed by Absi et al. (2015) and the adaptive large neighborhood search heuristic (ALNS) proposed by Adulyasak et al. (2014a) on these benchmark instances.

PT

The results with the best parameter settings for ALNS have been used (i.e., 500 iterations for the set B3 and 1000 iterations for the sets A1–A3 and B1–B2). The local search heuristic (H) proposed by Archetti et al.

CE

(2011) have not been included into the comparison, because H requires a correction in its implementation (Archetti, 2015). ALNS has been implemented on a 2.10 GHz Workstation with 2 GB RAM while IMVRP and IM-MultiTSP have been implemented on a 2.67 GHz Workstation with 3.48 GB RAM. Although

AC

different computational platforms were used for the experiments, they have similar efficiencies.

11

ACCEPTED MANUSCRIPT

Table 2: Summary of Benchmark Instances∗

4.1

B1 30 20 50 5 V C C C V 0 C

B2 30 20 100 9 V C C C V 0 C

B3 30 20 200 13 V C C C V 0 C

CR IP T

A3 480 6 100 ∞ C ∞ ∞ C 0 V C

AN US

Instance set A1 A2 Number of test instances 480 480 Number of periods (T ) 6 6 Number of retailers (n) 14 50 Number of trucks (m) 1 ∞ Demand (dit ) C C Production capacity (C) ∞ ∞ Plant inventory capacity (U0 ) ∞ ∞ Retailer inventory capacity (Ui ) C C Initial inventory at plant (I00 ) 0 0 Initial inventory at retailers (Ii0 ) V V Vehicle capacity (Q) C C ∗ V: Varying, C: Constant, ∞: Unlimited. (Table from Adulyasak et al., 2014a)

Results on Archetti et al. (2011) Instances

Each of the instance sets A1–A3 is composed of four classes. The 120 instances in the first class are the original instances and the rest of the instances are generated using these first class instances. The instances

M

in the second class have 10 times larger variable production cost b than the first class instances. Note that the fixed setup cost f is also affected by b as f = 10b in these instances. The third class instances have five

ED

times larger distribution cost than the first class instances. Finally, the fourth class instances do not have any inventory holding cost at the retailers.

Because the fleet size is one in the set A1 and unlimited in the sets A2 and A3, which means that there

PT

will be no feasibility issue for the solution obtained at the end of third phase, we set λ = 1. In order to determine the value of L, we have performed a preliminary experiment in which different values of L from

CE

three to nine as well as infinity (L = ∞ means constraints (27) are removed from the IF formulation) are tried in order to assess their impact on the solution quality. We have randomly selected 20 instances from each of the four classes in A2 for experimentation; thus, 80 instances have been selected in total among 480

AC

instances in A2. The results revealed that the average objective function value over these instances ranged from 448,790.7 to 448,808.0 (see Table 18 in appendix for details). These results are quite robust and the objective function values are close to each other for different values of L. We use the best result, L = 6, for the IF formulation in the fourth phase. In order to limit the computational time, we set a time limit of 1800

12

ACCEPTED MANUSCRIPT

Table 3: Summary of Computational Results for Sets A1–A3

CR IP T

A1 (n = 14) A2 (n = 50) A3 (n = 100) Heuristic Ave. Cost # Best Ave. Cost # Best Ave. Cost # Best 5P 179,313.6 290 587,602.1 309 1,084,957.3 317 IM-VRP – – 588,182.7 130 1,086,168.9 153 IM-MultiTSP 179,423.6 307 589,774.5 31 1,090,640.8 10 ALNS 181,802.9 2 590,210.0 13 1,089,588.6 0

seconds and a 0.01% integrality gap % (IG%) between the best upper (UB) and lower bounds (LB) found (i.e., IG% = 100x(U B − LB)/U B) for the solution of the IF formulation in the sets A2 and A3. The IF formulation is solved to optimality in the set A1. The F 0 formulation is solved to optimality in all instances. A summary of the computational results on sets A1–A3 is presented in Table 3, where the first column

AN US

indicates the name of the heuristic, the columns named as ‘Ave. Cost’ give the average of the objective function values obtained over 480 instances, and the columns named as ‘# Best’ give the number of the times a heuristic finds the smallest objective function value. Note that columns 2–3, 4–5, 6–7 show the results for the instances with n = 14, 50, and 100, respectively. The third column in Table 3 shows the number of the times a heuristic finds the optimal solution. The optimal solutions for the instances in set

M

A1 (i.e., single vehicle instances) have been found using the branch-and-cut algorithm in Ruokokoski et al. (2010). Note that IM-VRP and IM-MultiTSP are equivalent in the case of a single vehicle and their best

ED

performing variant is called the iterative heuristic with a multi-start scheme (IM-MS) in Absi et al. (2015). In other words, the results of IM-MultiTSP for the set A1 in Table 3 are obtained by IM-MS. The results presented in Table 3 reveal that 5P yielded superior solutions on average than the other heuristics in all

PT

sets. 5P also found significantly more best solutions than the other heuristics in all sets except the set A1. In particular, 5P managed to find new best solutions in 623 instances out of 960 in sets A2 and A3. Note

CE

that 5P found the same best solution with others in three instances of set A2. The detailed results of IM-MS along with ALNS and 5P on the set A1 have been reported in Tables 4 and 5. The computational results reveal that 5P provided the best results in terms of percent average gap

AC

between the heuristic solution value and the optimal solution value and computational times for all classes except class IV with a slight margin. The average gap percentages with respect to the best known solutions and the average computational

time for the sets A2 and A3 are presented in Tables 6 and 7, respectively. It can be seen in Table 6 that the

13

ACCEPTED MANUSCRIPT

Table 4: Average Gap% with respect to the Optimal Solutions for Set A1 ALNS 1.70 0.36 8.43 0.93 2.85

IM-MS 0.09 0.01 0.57 0.02 0.17

5P 0.03 0.00 0.18 0.03 0.06

CR IP T

Classes Class I Class II Class III Class IV All

Table 5: Average Computational Times (seconds) for Set A1 ALNS 9.2 8.9 7.6 8.7 8.6

IM-MS 251.0 214.2 237.2 216.9 229.8

5P 5.0 4.9 4.7 5.2 5.0

AN US

Classes Class I Class II Class III Class IV All

percent average gaps of 5P with respect to the best known solution values are the lowest compared to those of three competitor heuristics. In terms of the computational time, 5P is the fastest heuristic in the set A2

M

whereas it is the second fastest heuristic together with ALNS after IM-VRP in the set A3 according to Table 7. Although 5P needed on average 2.2 more minutes than IM-VRP for A3 instances, solving large-scale PRP

ED

instances in less than four minutes is quite acceptable. The average computational times in seconds required by each phase of 5P are given in Table 8, where the first column indicates the instance set and the rest of the columns designates the computational times

PT

required by each phase of 5P. Note that the IF formulation in the fourth phase of 5P could not be solved to optimality for 18 instances in the set A3 within the time limit with the largest IG% = 0.12%. Apart from

CE

computational times of phases 2 and 4 in 5P, other phases are comparable in terms of computational times for all sets.

Results on Boudia et al. (2007) Instances

AC

4.2

We have performed a sensitivity analysis on randomly chosen six instances among 30 instances in B1 to determine the best λ and L pair. We have tested different pairs of λ and L where we set λ = 0.95 + 0.01Λ, Λ = 0, 1, ..., 5 and L = 3, 4..., 9 and L = ∞. The average objective function values of six instances for

14

ACCEPTED MANUSCRIPT

Table 6: Average Gap% with respect to the Best Known Solutions for Sets A2 and A3 ALNS 1.13 0.17 3.52 0.19 1.25 0.96 0.13 3.56 0.29 1.23

IM-MultiTSP 1.09 0.09 2.73 0.46 1.09 1.82 0.11 3.72 0.66 1.58

IM-VRP 0.17 0.05 0.99 0.10 0.33 0.18 0.03 1.14 0.21 0.39

5P 0.05 0.01 0.09 0.06 0.05 0.03 0.01 0.06 0.01 0.03

CR IP T

A3

Classes Class I Class II Class III Class IV All Class I Class II Class III Class IV All

AN US

Sets A2

Table 7: Average Computational Times (seconds) for Sets A2 and A3 ALNS 50.2 49.5 42.7 44.1 46.6 228.5 217.6 197.8 206.0 212.5

CE

PT

ED

A3

Classes Class I Class II Class III Class IV All Class I Class II Class III Class IV All

IM-MultiTSP 338.5 235.7 317.9 375.8 317.0 514.2 497.4 509.3 507.0 507.0

M

Sets A2

IM-VRP 25.6 21.7 22.6 27.7 24.4 85.5 76.1 75.1 86.6 80.8

5P 16.4 13.8 15.8 25.4 17.8 323.5 50.6 349.6 125.1 212.2

AC

Table 8: Average Computational Times (seconds) of Each Phase for Sets A1–A3 Instance Set A1 A2 A3

Phase 1 0.6 0.6 1.4

Phase 2 0.9 4.7 33.0

15

Phase 3 1.6 1.0 2.8

Phase 4 0.3 6.3 166.0

Phase 5 1.6 5.2 9.1

ACCEPTED MANUSCRIPT

each λ and L pair are close to each other and contained in the interval [339,556.3, 343,975.7] (see Table 19 in appendix for details). Based on the best result in these experiments, we set λ = 1 for the formulation F 0 in the second phase and L = ∞ for the formulation IF in the fourth phase. In order to control the computational time, we set a time limit of 3600 (1), 7200 (2), and 14400 (4) seconds (hours) for the solution of the F 0 formulation in the sets B1, B2, and B3, respectively. Likewise, we set a 3600 seconds (1 hour) time limit and a 0.1% integrality gap between the best upper and lower bounds found for the solution of the IF

CR IP T

formulation.

The summary of the computational results on the instances introduced by Boudia et al. (2007) are presented in Table 9, where the first column indicates the name of the heuristic, the columns named as ‘Ave. Cost’ give the average of the objective function values obtained over 30 instances, and the columns named as ‘# Best’ give the number of the times a heuristic finds the smallest objective function value. Table

AN US

10 shows the average gap percentage over 30 instances with respect to the best known solutions for each instance set and for each heuristic. Based on the results given in Tables 9 and 10, it is easy to observe that 5P significantly outperformed the other three heuristics in sets B1 and B2 in terms of both the number of best solutions found and the average gap percentage with respect to the best known solutions. However, on the largest set B3, 5P yielded the poorest average performance over 30 instances among these heuristics

M

despite the fact that it found the best solution in seven of these instances.

Next we compare the average computational times of the heuristics in Table 11, where the first column indicates the instance set and the rest of the columns the computational times (in seconds) required by each

ED

heuristic. As seen in Table 11, the computational times of 5P and IM-MultiTSP are quite close to each other and they are substantially higher than those of the rest. Thus, one can say that 5P performs well

PT

in sets B1 and B2 compared to existing heuristics at the expense of higher computing times. Considering the complexity of the PRP and the sizes of instances solved, solving instances in sets B1 and B2 within 2.1 hours (i.e., 7487.5 seconds) on average to obtain good-quality solutions is acceptable. On the other hand,

CE

5P needs much more time than IM-VRP in B3 and gives worse results than IM-VRP. The reasons of the obtained results in B3 will be elaborated in the sequel.

AC

In order to shed light on the performance and computational time requirements of 5P, the average computational times in seconds for each phase of 5P and the average integrality gap % for the MIPs in the second and fourth phases in sets B1–B3 are presented in Table 12, where the first column indicates the instance set, columns 2–6 the computational times required by the phases 1–5 of 5P, respectively, and columns

16

ACCEPTED MANUSCRIPT

Table 9: Summary of Computational Results for Sets B1–B3

AN US

CR IP T

B1 (n = 50) B2 (n = 100) B3 (n = 200) Ave. Cost # Best Ave. Cost # Best Ave. Cost # Best 5P 343,439.1 27 631,351.4 26 882,128.5 7 IM-VRP 348,990.2 0 636,294.5 2 865,905.5 18 IM-MultiTSP 348,175.2 1 639,271.7 1 879,447.1 4 ALNS 346,878.3 2 636,962.1 1 876,761.4 1

Table 10: Average Gap% with respect to the Best Known Solutions for Sets B1–B3 IM-MultiTSP 1.41 1.32 1.93

M

ALNS 1.03 0.96 1.61

IM-VRP 1.64 0.85 0.36

5P 0.02 0.07 2.24

PT

ED

Instance Set B1 B2 B3

AC

CE

Table 11: Average Computational Times (seconds) for Sets B1–B3 Instance Set B1 B2 B3

ALNS 481.3 1569.9 5794.2

IM-MultiTSP 1653.2 9483.5 19270.1

17

IM-VRP 550.9 2054.2 4197.0

5P 2464.3 7487.5 16365.3

ACCEPTED MANUSCRIPT

Table 12: Average Integrality Gaps% and Computational Times (seconds) of Each Phase for Sets B1–B3 Instance Set B1 B2 B3

Phase 1 0.6 1.0 2.4

Phase 2 2423.9 (10) 7200.2 (30) 14400.3 (30)

Phase 3 5.5 12.9 24.6

Phase 4 19.1 242.1 1899.0

Phase 5 15.2 31.4 38.9

IG-F 0 % 0.08 2.42 6.95

IG-IF % 0.08 0.09 0.09

CR IP T

7–8 the average integrality gap % for the MIPs in second and fourth phases, respectively. The numbers in parenthesis in the second column show the number of instances that F 0 could not be solved to optimality within the time limit. As indicated by Table 12, in terms of computational time, the most demanding phase of 5P is the second phase. A significant amount of time is needed to obtain a good solution for the F 0 formulation in the second phase as the size of the instances gets larger. Furthermore, the F 0 formulation

AN US

could not be solved to optimality within the time limits in all instances of sets B2 and B3, and the average percent gap between the best found upper and lower bounds for F 0 increases significantly with an increase in the size of the instances. The computational requirements in phases 1, 3, and 5 are negligible compared to those in phases 2 and 4. The IF formulation was solved within the one-hour time limit with a very small integrality gap even for the largest instances. The seventh column of Table 12 reveals that the poor average performance of 5P in set B3 may be due to not being able to solve the F 0 formulation near-optimality as

M

the remaining IG-F 0 % is not small. Therefore, we performed further computational experiments on the first six instances of set B3 with greater time limits set for the solution of F 0 .

ED

We present the results of these additional computational experiments in Table 13, where the first column indicates the instance number, column 2 the best known upper bounds (obtained by IM-VRP) before this

PT

study, columns 3–5 the upper bounds obtained by 5P when F 0 is solved with the time limits of four, six, and eight hours, respectively, and columns 6–8 the average integrality gap % for F 0 . According to the results in Table 13, the solution quality of 5P is the best in five of these six instances when eight hours of time is

CE

allowed for the solution of F 0 . Thus, we can see that better upper bounds could be obtained by 5P if it was possible to solve F 0 near-optimality in reasonable times. When the time limit increased, the upper bounds

AC

obtained by 5P became better on average and the remaining integrality gaps for F 0 usually decreased. Note that there are some instances in which the quality of the solution given by 5P decreased as the time limit for F 0 increased. This is due to the fact that the TSP solver we used (i.e. Concorde) can find alterative optimal tours with the same optimal tour length in different runs and using a different giant tour in F 0 can yield a

different solution.

18

ACCEPTED MANUSCRIPT

Table 13: Computational Results for the First Six Instances of Set B3

5

Conclusion

5P-8 hours 4 hours 861,423 4.81 853,406 9.95 854,293 3.98 849,737 6.26 850,380 4.36 869,635 8.32 856,479 6.28

IG-F 0 % 6 hours 3.82 5.25 6.37 4.84 1.12 6.23 4.60

8 hours 2.37 1.98 1.91 1.67 1.47 4.23 2.27

CR IP T

Instance IM-VRP 1 875,836 2 859,356 3 860,738 4 860,695 5 861,275 6 872,453 Average 865,059

Upper Bounds 5P-4 hours 5P-6 hours 866,373 874,580 907,265 873,904 866,063 872,673 865,499 866,189 862,217 847,752 887,197 883,171 875,769 869,712

AN US

We have proposed an effective multi-phase heuristic for the production routing problem. The computational results have shown that our heuristic yielded favorable results compared to the existing ones in the literature on the benchmark instances. In particular, our heuristic managed to improve the best known results for the majority of the benchmark instances available in the literature at the expense of higher computing times on large instances.

A good feature of our heuristic is that because it is a mathematical programming-based heuristic, it can

M

easily be adapted for different variants of the problem (e.g., the IRP, the PRP with multiple products, the PRP with order-up-to level policy at retailers, the PRP with time-varying costs, etc.).

ED

As a new perspective, tailored algorithms to solve the PRP with the imposed a priori tour in the second phase and the restricted version of the PRP in the fourth phase can improve the quality of the obtained

PT

solution and can reduce the solution time.

CE

Appendix. Detailed Computational Results In this section, we provide the objective function values found by our heuristic on the instances introduced

AC

by Archetti et al. (2011) and Boudia et al. (2007). We also present the results of the sensitivity analysis which enabled us to fine-tune important parameters (i.e., λ and L) of our heuristic. The objective function values found by 5P for the sets A1–A3 are presented in Tables 14–16, where N

and I denote the instance group and instance number, respectively. Table 17 provides the objective function values found by 5P for the sets B1–B3. All instances, the optimal objective function values for the set

19

ACCEPTED MANUSCRIPT

A1, and the results obtained by all heuristics in this study are provided at http://tol2.metu.edu.tr/testinstances-production-routing-problem. The results of the sensitivity analysis performed on the 80 instances of A2 are presented in Table 18. Similarly, Table 19 provides the results of the sensitivity analysis performed on six instances of B1.

CR IP T

References Absi, N., C. Archetti, S. Dauzere-Peres, D. Feillet. 2015. A two-phase iterative heuristic approach for the production routing problem. Transportation Science 49(4) 784–795.

Adulyasak, Y., J. F. Cordeau, R. Jans. 2014a. Optimization-based adaptive large neighborhood search for

AN US

the production routing problem. Transportation Science 48(1) 20–45.

Adulyasak, Y., J. F. Cordeau, R. Jans. 2014b. Formulations and branch-and-cut algorithms for multivehicle production and inventory routing problems. INFORMS Journal on Computing 26(1) 103–120. Adulyasak, Y., J. F. Cordeau, R. Jans. 2015. The production routing problem: A review of formulations and solution algorithms. Computers & Operations Reseach 55 141–152.

M

Applegate, D. L., R. E. Bixby, V. Chavatal, W. J. Cook. 2007. The traveling salesman problem: A computational study. Princeton University Press, New Jersey.

ED

Archetti, C. 2015. Private Communication.

Archetti, C., L. Bertazzi, A. Hertz, M. G. Speranza. 2012. A hybrid heuristic for an inventory routing

PT

problem. INFORMS Journal on Computing 24(1) 101–116. Archetti, C., L. Bertazzi, G. Paletta, M. G. Speranza. 2011. Analysis of the maximum level policy in a

CE

production-distribution system. Computers & Operations Reseach 38(12) 1731–1746. Armentano, V. A., L. Shiguemoto A, A. Lokketangen. 2011. Tabu search with path relinking for an integrated

AC

production-distribution problem. Computers & Operations Reseach 38(8) 1199–1209.

Bard, J. F., N. Nananukul. 2009a. Heuristics for a multiperiod inventory routing problem with production decisions. Computers & Industrial Engineering 57(3) 713–723.

20

ACCEPTED MANUSCRIPT

Table 14: Objective Function Values Found by 5P for the Set A1 1 66907 69184 73768 108475 110504 115574 57425 59454 64524 97225 99435 106294 95311 99540 108800 136879 140856 150606 85829 89806 99556 125629 129856 141326 26543 26979 28558 65518 66822 70328 29556 30020 32059 68736 69920 73829 226343 226779 228358 598318 599634 603128 229356 229820 231859 601536 602720 606629

2 62818 63767 67283 97698 98575 102278 54703 55580 59012 87462 89045 93702 93545 95210 101114 128435 130028 137216 85440 87033 92447 118189 120758 127137 22703 23158 23862 55110 55802 57582 25767 26341 27494 58438 59383 61214 189012 189478 190182 498630 499322 501102 192110 192661 193814 501935 502903 504734

21

I 3 52770 53656 55794 83394 84601 86328 43909 45116 46843 74179 74956 77243 75239 77815 81492 105851 108435 111784 66366 68950 72299 96636 98790 102699 20125 20190 20844 49111 49590 51209 22580 22675 23333 51848 52159 53832 172135 172200 172854 454471 454950 456569 174590 174685 175477 457208 457519 459192

4 64097 66719 72026 102010 104142 110054 55535 57667 63579 91381 94017 100884 93894 98664 109297 131728 135933 147325 85253 89582 100850 121099 125808 137840 24411 24787 26119 59662 60477 63129 27457 28053 29915 62938 63937 67041 205311 205713 207019 542061 542877 545529 208323 208953 210815 545338 546337 549441

CR IP T

N 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96

AN US

5 34814 35164 35806 70799 71156 72233 27444 27801 28878 60580 61301 63943 38112 38810 40387 73722 74352 76313 30367 30997 32958 63860 64497 68023 204644 204989 205636 523679 524036 525113 197274 197631 198708 513468 514181 516823 207942 208640 210217 526602 527232 529193 200197 200827 202788 516740 517377 520903

M

4 36970 37185 37785 75284 75601 76346 28809 29126 29871 64453 65128 67025 40068 40917 41685 78265 78809 80074 31790 32334 33599 67597 68461 70964 218020 218286 218685 557684 558001 558746 209709 210026 210771 546853 547528 549425 221116 221817 222585 560665 561209 562474 212690 213234 214499 549997 550861 553364

ED

I 3 31107 31328 31476 63090 63124 63403 23605 23639 23918 52836 53113 54383 33992 34169 34486 65387 65509 65949 25902 26024 26464 55611 55706 56864 183117 183338 183486 468450 468484 468763 175615 175649 175928 458196 458473 459743 186002 186179 186496 470747 470869 471309 177912 178034 178474 460971 461066 462224

PT

2 35183 35428 35691 70056 70214 70789 27061 27219 27794 59725 60290 61534 38260 38643 39355 73133 73395 74345 30138 30400 31350 62900 63774 65169 201503 201748 202011 513576 513734 514309 193381 193539 194114 503245 503810 505054 204580 204963 205675 516653 516915 517865 196458 196720 197670 506420 507294 508689

CE

1 40390 40600 41278 82787 83212 83978 31893 32162 32928 71289 72096 74698 43735 44506 45510 85781 86277 87516 34731 35227 36466 74460 75137 78236 240619 240835 241078 615743 616012 616778 231693 231962 232728 604089 604896 607498 243854 244442 245310 618581 619077 620316 234531 235027 236266 607260 607937 611036

AC

N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

5 61636 64411 73168 97082 99907 108570 53727 56522 65215 87587 90026 100446 90854 96318 113469 126300 131668 148871 82945 88313 105516 116805 121873 140821 22893 23280 25090 55613 56776 60155 26027 26464 29146 58879 59958 64211 192723 193110 194920 508493 509660 513035 195857 196288 198976 511759 512838 517091

ACCEPTED MANUSCRIPT

Table 15: Objective Function Values Found by 5P for the Set A2 1 147171 152422 163676 263721 268972 280283 120133 125146 137354 236683 241696 252908 183179 191622 213301 299729 310144 332610 151899 166314 187538 272169 282738 306368 76883 78213 80080 191981 193097 194871 84061 86130 90223 200032 201987 206518 679253 680466 682450 1798301 1799417 1801257 686715 688490 692518 1806352 1808329 1812758

2 137695 143352 154660 252295 257956 268079 112295 117415 128139 226895 232015 241758 168192 179877 200548 282792 294477 315948 142438 153528 173210 257038 266805 288107 74864 75452 77537 188840 189960 191717 80812 82910 86690 195578 196901 201293 666916 668192 670016 1767080 1768231 1769957 672949 675058 679155 1774349 1776627 1779526

22

I 3 133299 138562 147222 246399 251662 260735 110085 115305 123989 223184 228405 237170 166487 176645 194956 279587 289745 308800 141282 150316 169529 256356 265204 283562 74362 75377 77276 186380 187337 189391 80851 82688 86326 194076 195778 199287 658293 659270 661191 1742995 1743977 1745989 664879 666615 670328 1750764 1752858 1756654

4 152425 161048 174895 281475 290098 302896 127592 135153 147706 256596 263635 276341 191874 207165 235288 320924 336787 362396 165760 179479 205461 294810 309835 333253 84765 86529 89150 214137 215285 217830 92433 95445 100809 222004 224053 230087 755709 757308 760001 2001175 2002229 2005100 763308 766221 772411 2010357 2012170 2018552

CR IP T

N 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96

AN US

5 103318 104414 107437 214266 215928 219353 83660 85222 88225 191980 193472 196629 113182 116538 122292 222305 225601 232604 91770 94881 101456 200025 203131 209259 668352 670014 672986 1720028 1721538 1724670 646785 648267 651214 1683738 1685308 1688485 678299 682120 687587 1727479 1730819 1736774 654652 657720 663527 1692406 1695619 1702106

M

4 117830 119181 121555 251180 252872 255419 96508 98212 100759 225561 227262 229809 129059 131707 136279 258755 261966 267259 104394 107456 112599 233444 236583 241649 796767 798221 800603 2059798 2060963 2063974 773433 774911 777689 2020509 2021469 2024517 807914 810188 814806 2067046 2069937 2075322 780992 783870 789125 2028431 2031567 2036973

ED

I 3 103963 104743 106541 219951 220937 222760 84070 85053 86880 197171 198153 199980 112645 114865 118688 226522 228553 232311 90606 92605 96158 203706 205705 209248 695550 696374 698270 1795108 1795901 1797765 673038 673877 675791 1759143 1760142 1761999 704360 706201 710579 1801745 1803222 1806954 679671 680808 684931 1766118 1767802 1771589

PT

2 110290 110978 112414 228197 229303 231284 87502 88579 90589 202102 203179 205170 118480 120930 124256 234476 236487 240463 93644 95658 99651 208244 210258 214251 710197 710982 712390 1825789 1827069 1828904 685263 686458 688371 1787084 1788113 1789740 718372 720780 724104 1831910 1833809 1838039 691395 693273 697497 1793374 1795615 1799536

CE

1 114470 115536 117416 234959 236005 238288 91249 92295 94578 207799 208845 211128 124249 126223 129939 242244 244240 248802 98429 100530 105092 215084 217080 221642 725286 726552 728238 1862175 1863260 1865547 700244 701374 703287 1820666 1821672 1823892 734830 736496 741292 1868968 1871259 1876172 707047 709367 713612 1828211 1830646 1834391

AC

N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

5 138702 146345 163383 247182 254595 270650 115842 123915 138518 224491 232165 246750 179322 194687 228099 287434 302937 335965 155877 170995 201313 263989 277943 310462 72709 74361 77517 178998 180733 183465 80426 84350 90538 187847 191052 197133 630259 631908 635067 1665798 1667620 1670353 637827 641768 648088 1674769 1677695 1683933

ACCEPTED MANUSCRIPT

Table 16: Objective Function Values Found by 5P for the Set A3 1 263087 274450 300807 481311 493983 518043 217044 229871 253905 435847 447521 472673 326696 350036 401246 542879 569116 617676 279847 304990 354495 497675 520436 571836 143280 145880 150662 362525 364578 369460 155874 160173 169573 374668 379024 388794 1301758 1303998 1308936 3446939 3449003 3453920 1314231 1318993 1328737 3461148 3465568 3475430

2 252218 264421 287056 456972 468183 491038 204441 215858 238334 408889 419708 442244 311937 332297 377874 515703 538667 582969 262161 283760 331524 466011 487610 532623 133806 135807 139972 337715 339529 343941 144762 149109 158010 348603 353235 361778 1209250 1211520 1216143 3203450 3205539 3210083 1221835 1225782 1234904 3216525 3220313 3229582

23

I 3 262468 277672 313147 476782 491843 522675 216157 232872 267032 429429 444778 474807 330259 361247 426489 544003 576545 635645 286129 313161 379198 495425 526193 589810 140998 143750 149684 354236 357373 362875 154227 160057 171720 367621 373334 385082 1268208 1270924 1276647 3355918 3358580 3364664 1282374 1288430 1300961 3371355 3376551 3388252

4 249333 257846 281010 446408 455391 479049 196643 207353 229167 395056 403784 426708 301195 323008 366901 500907 521108 562565 250692 268901 315574 450689 466512 513340 128966 130904 135523 326600 328726 332977 139398 143915 152509 337406 341070 350232 1171057 1172993 1177294 3102183 3104370 3108656 1181686 1186172 1194871 3114292 3118274 3126829

CR IP T

N 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96

AN US

5 205786 208410 213719 435946 438642 444237 169771 172955 178540 391121 394241 400063 224280 230614 240491 451699 458097 468539 182994 189276 200570 404344 410626 421994 1395912 1398343 1403403 3601116 3603968 3609604 1353464 1356133 1361850 3526437 3528941 3534650 1413672 1419551 1430654 3614386 3619868 3631729 1366805 1372126 1383818 3541422 3547189 3558327

M

4 195940 197335 201798 402684 404501 408397 153929 155746 160211 351529 353346 357811 211773 214927 223279 413435 417366 426405 164680 168711 177771 362280 366174 374908 1256179 1257679 1261727 3221967 3224261 3228250 1208619 1210517 1214784 3141888 3143905 3148072 1271317 1275136 1283161 3232529 3236420 3245026 1219429 1223082 1231966 3154099 3157884 3167025

ED

I 3 198475 200937 206927 420746 423451 429251 163263 166427 172334 376648 379627 385738 216689 222231 233494 437250 442831 453723 176681 182720 195376 389824 395760 408328 1343153 1346127 1351536 3466384 3469768 3475739 1301793 1305130 1311091 3392522 3395127 3401082 1361345 1367334 1378537 3480016 3486440 3498131 1315396 1321788 1333704 3407585 3412872 3425079

PT

2 195207 197046 201709 407705 409918 414616 158721 160836 165378 362571 364779 369080 211186 215568 224718 422281 426713 435389 170054 174318 183639 373957 378313 387232 1289167 1290969 1295175 3318374 3320805 3324683 1246776 1249103 1253459 3243201 3245307 3249873 1304263 1309319 1317427 3329777 3334140 3343395 1258218 1262456 1271980 3256612 3260661 3269148

CE

1 200990 203275 207989 429194 431328 436030 167393 169704 174814 386193 388539 393597 219814 223756 234049 444080 448728 459095 179918 184839 194308 398742 403390 413140 1376231 1378352 1382908 3556480 3558928 3563908 1337298 1340020 1344362 3485382 3487784 3492521 1393687 1398581 1408311 3568751 3574334 3583346 1349592 1354819 1364192 3500325 3504614 3514359

AC

N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

5 269561 284911 312635 491264 507084 535303 222443 237954 265534 443792 458701 487001 336457 366924 429151 557809 586644 643973 287726 319288 379699 508494 539400 595743 145894 148669 154498 367063 370072 375517 158771 164510 176190 380447 386324 397056 1317663 1320516 1325955 3487870 3490667 3496244 1331304 1336780 1348698 3503317 3508755 3519495

ACCEPTED MANUSCRIPT

Table 17: Objective Function Values Found by 5P for Sets B1–B3

PT

ED

B3 866,373 907,265 866,063 865,499 862,217 887,197 883,939 859,674 930,653 896,542 864,602 886,261 881,990 868,718 846,854 935,476 854,055 857,877 869,641 873,565 880,292 867,850 875,685 871,828 875,332 924,333 875,154 919,585 863,724 945,610

CR IP T

B2 627,481 636,179 626,542 620,637 640,475 625,173 618,610 614,107 629,671 631,606 633,159 639,171 628,002 630,926 629,486 635,671 648,858 639,500 638,363 646,467 632,601 639,488 624,745 642,607 634,006 624,516 622,326 620,148 637,569 622,453

AN US

B1 340,313 357,608 344,605 340,510 348,207 348,338 349,067 347,139 345,314 335,588 332,822 337,257 343,915 347,021 347,469 338,888 346,411 353,899 347,945 337,715 332,891 329,150 344,988 337,110 349,216 345,900 337,618 337,063 361,210 337,995

M

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

AC

CE

Table 18: The average objective function values of 80 instances from A2 found by 5P for different values of L and λ = 1 L 3 4 5 6 7 8 9 ∞

Average objective function value 448,806.9 448,794.3 448,797.2 448,790.7 448,793.8 448,808.0 448,804.9 448,799.6

24

ACCEPTED MANUSCRIPT

Table 19: The average objective function values of six instances from B1 found by 5P for different values of L and λ λ 0.95 341,588.0 342,023.0 341,816.5 341,857.3 341,899.8 341,988.5 341,686.2 341,828.0

0.96 341,944.7 341,215.8 341,557.5 341,515.3 341,888.7 341,340.0 341,999.2 341,960.5

0.97 341,684.7 341,284.8 341,159.2 340,892.0 340,956.3 341,477.3 341,288.3 341,472.0

0.98 343,499.8 343,975.7 343,800.2 341,578.2 341,393.7 341,481.0 341,235.0 341,138.7

0.99 342,187.2 340,562.0 341,574.2 340,794.8 340,523.2 340,956.8 341,124.0 340,871.3

1.00 343,240.2 340,792.4 339,626.2 339,871.3 340,281.3 340,087.7 339,741.8 339,556.3

CR IP T

L 3 4 5 6 7 8 9 ∞

Bard, J. F., N. Nananukul. 2009b. The integrated production–inventory–distribution-routing problem. Jour-

AN US

nal of Scheduling 12(3) 257–280.

Bard, J. F., N. Nananukul. 2010. A branch-and-price algorithm for an integrated production and inventory routing problem. Computers & Operations Reseach 37(12) 2202–2017.

Bertazzi, L., G. Paletta, M. G. Speranza. 2005. Minimizing the total cost in an integrated vendor-managed inventory system. Journal of Heuristics 11(5) 393–419.

M

Boudia, M., M.A. O. Louly, C. Prins. 2007. A reactive grasp and path relinking for a combined production— distribution problem. Computers & Operations Reseach 34(11) 3402–3419.

ED

Boudia, M., M.A. O. Louly, C. Prins. 2008. Fast heuristics for a combined production planning and vehicle routing problem. Production Planning & Control 19(2) 85–96.

PT

Boudia, M., C. Prins. 2009. A memetic algorithm with dynamic population management for an integrated production—distribution problem. European Journal of Operational Research 195(3) 703–715.

CE

Chandra, P. 1993. A dynamic distribution model with warehouse and customer replenishment requirements. Journal of the Operational Research Society 44(7) 681–692.

AC

Chandra, P., M. L. Fisher. 1994. Coordination of production and distribution planning. European Journal of Operational Research 72(3) 503–517.

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

25

ACCEPTED MANUSCRIPT

Coelho, L. C., G. Laporte. 2013. A branch-and-cut algorithm for the multi-product multi-vehicle inventoryrouting problem. International Journal of Production Research 51(23-24) 7156–7169. Fumero, F., C. Vercellis. 1999. Synchronized development of production, inventory and distribution schedules. Transportation Science 33(3) 330–340.

Mathematical Programming Computation 2(2) 79–101.

CR IP T

Groer, C., B. Golden, E. Wasil. 2010. A library of local search heuristics for the vehicle routing problem.

Lei, L., S. Liu, A. Ruszczynski, S. Park. 2006. On the integrated production, inventory and distribution routing problem. IIE Transactions 38(11) 955–970.

¨ H. S¨ Pınar, O., ural. 2006. Coordinating inventory and transportation in vendor managed systems. R. Meller

AN US

et al., ed., Proceedings of the Material Handling Research Colloquium 2006 . 459–474.

Ruokokoski, M., O. Solyalı, J. F. Cordeau, R. Jans, H. S¨ ural. 2010. Efficient formulations and a branch-andcut algorithm for a production-routing problem. Tech. Rep. G-2010-66, HEC Montreal, Canada. Solyalı, O., H. S¨ ural. 2009. A relaxation based solution approach for the inventory control and vehicle routing

and Optimization. World Scientific.

M

problem in vendor managed systems. S.K. Neogy, A.K. Das, R.B. Bapat, eds., Modeling, Computation

ED

Solyalı, O., H. S¨ ural. 2011. A branch-and-cut algorithm using a strong formulation and an a priori tour

AC

CE

PT

based heuristic for an inventory-routing problem. Transportation Science 45(3) 335–345.

26