Vehicle flow formulation for two-echelon time-constrained vehicle routing problem

Accepted Manuscript Vehicle flow formulation for two-echelon time-constrained vehicle routing problem Hongqi Li, Ming Bai, Yibin Zhao, Changzhi Dai

PII:

S2096-2320(19)30077-0

DOI:

https://doi.org/10.1016/j.jmse.2019.05.006

Reference:

JMSE 11

To appear in:

Journal of Management Science and Engineering

Please cite this article as: Li H., Bai M., Zhao Y. & Dai C., Vehicle flow formulation for two-echelon time-constrained vehicle routing problem, Journal of Management Science and Engineering, https:// doi.org/10.1016/j.jmse.2019.05.006. 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

Vehicle flow formulation for two-echelon time-constrained vehicle routing problem

[email protected]；[email protected]

SC

* Corresponding author: [email protected]

RI PT

Hongqi Li *, Ming Bai, Yibin Zhao, Changzhi Dai School of Transportation Science and Engineering, Beihang University. No. 37 Xueyuan Road, Haidian District, Beijing, 100191, China； [email protected]；

Abstract

EP

TE D

M AN U

Two-echelon routing problems, including variants such as the two-echelon vehicle routing problem (2E-VRP) and the two-echelon location routing problem (2E-LRP), involve assignment and location decisions. However, the two-echelon time-constrained vehicle routing problem (2E-TVRP) that caters to from-linehaul-to-delivery practices does not involve assignment decisions. This routing problem variant for networks with two echelons has not yet attracted enough research interest. Localized or long-distance services suffer from the lack of the assignment decisions between satellites and customers. Therefore, the 2E-TVRP, rather than using assignment decisions, adopts time constraints to decide the routes on each of the two interacting echelons: large-capacity vehicles transport cargoes among satellites on the first echelon, and small-capacity vehicles deliver cargoes from satellites to customers on the second echelon. This study introduces a mixed integer linear programming model for the 2E-TVRP and proposes a heuristic algorithm that incorporates the savings algorithm followed by a variable neighborhood search phase. Illustrative examples are used to test the mathematical formulation and the heuristic and a case study is used to demonstrate that the heuristic can effectively solve realistic-size instances of the 2E-TVRP.

Keywords: Vehicle routing; Two-echelon; Time constraints; Mixed integer linear programming;

AC C

Variable neighborhood search

1. Introduction

The loading capacities allow vehicles to be classified into two types: large-capacity vehicles that are fit for intercity linehaul transportation, and small-capacity vehicles that are fit for urban delivery. As large amounts of cargoes are generally transported among cities via large-capacity vehicles that are often forbidden inside cities, cargo transshipments are necessary between the large-capacity vehicles and small-capacity delivery vehicles. In the two-echelon city logistics systems reported in the literature (e.g., Crainic et al., 2009), the first echelon involves vehicles delivering cargoes from depots to satellites, and the second echelon involves vehicles delivering 1

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

cargoes from satellites to customers. At the satellites, the cargoes are transshipped between different vehicles. The two-echelon systems can be converted into cost-effective city logistics systems by integrating the underlying freight management and routing planning (Morais et al., 2014). Furthermore, varying the vehicle types in different echelons, allows the two-echelon logistics networks to save costs (Dondo et al., 2011). The family of routing problems for networks with two echelons, in the study that covers the two-echelon vehicle routing problem (2E-VRP) and the two-echelon location routing problem (2E-LRP) variants, has been referred to plan the routes of the two-echelon city logistics systems. In recent years, a considerable number of studies focusing on the family of routing problems for networks with two echelons have been published; however, Cuda et al. (2015) stated that the research is still relatively unexplored despite the important practical applications to city logistics systems. In the 2E-VRP, there is one depot and some capacitated satellites to serve customers. Decisions on the assignments between satellites and customers are used to formulate the 2E-VRP. In the 2E-LRP, there are alternative depots and satellites. The 2E-LRP makes location decisions in depots and satellites whereas the 2E-VRP only makes decisions in satellites. In addition, the truck and trailer routing problem (TTRP) is regarded as a special variant of the 2E-VRP in the literature (e.g., Cuda et al., 2015). Some exact algorithms were used to solve the 2E-VRP (e.g., Perboli and Tadei, 2010; Baldacci et al., 2013; Santos et al., 2013; Jepsen et al., 2013; Santos et al., 2014; Sitek and Wikarek, 2014; Perboli et al., 2018; Liu et al., 2018). In the literature, the 2E-VRP was solved mostly by heuristics, in which the assignment decisions between satellites and customers are known in advance and then the feasible solution is constructed. To create neighborhoods, the customer-to-satellite assignments are changed. Algorithms such as the multi-start heuristics (e.g., Crainic et al., 2011), adaptive large neighborhood search heuristic (e.g., Hemmelmayr et al., 2012), and variable neighborhood search (VNS) (e.g., Wang et al., 2017; Belgin et al., 2018) were adopted to solve the 2E-VRP. On the contrary, the 2E-LRP has seldom been solved by exact algorithms (e.g., Contardo et al., 2012; Diabat and Theodorou, 2015). Heuristics were generally used to solve the 2E-LRP (e.g., Boccia et al., 2010; Nguyen et al., 2010, 2012a, 2012b; Dalfard et al., 2013; Breunig et al., 2016; Pichka et al., 2018; Wang et al., 2018). In the literature, the assignment decisions between satellites and customers, which are denoted as binary variables, are the key to construct models for the routing problems of networks with two echelons. In practice, however, the assignment of customers to satellites is unnecessary when satellites have their own localized customers or the transport distance between any two satellites is too long. Typically, the logistics and supply-chain enterprises synchronously operate the intercity linehaul and urban delivery. Therefore, instead of the assignment decisions, other variables or constraints are necessary to model the routing problems of networks with two echelons. For example, Grangier et al. (2016) addressed the two-echelon multiple-trip vehicle routing problem with satellite synchronization that integrated synchronization constraints, time windows constraints, and multiple trips at the second echelon. Additionally, it was assumed that the first-echelon vehicles and the second-echelon vehicles must be at the satellite at the same time during the transshipment operation. To achieve this, an adaptive large neighborhood search was introduced. Li et al. (2016) proposed time constraints to decide routes on each of the two interacting echelons and to manage the dependence of the second echelon on the first one. For this, a savings-based algorithm(CW) improved by a local search phase was adopted. Furthermore, Li et

2

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

al. (2018) introduced the two-echelon distribution system considering the variations in the real-time transshipment capacity, particularly for satellites having transshipment and consolidation operations. A mixed integer linear programming model was proposed. A savings-based algorithm followed by the variable neighborhood search phase was provided. This study caters to the practice that logistics and supply-chain enterprises synchronously operate the intercity linehaul and urban delivery, and focuses on the variant called the two-echelon time-constrained vehicle routing problem (2E-TVRP). The 2E-TVRP considers the synchronous routing problems on the linehaul echelon that involves large-capacity vehicles transporting cargoes among satellites and those on the delivery echelon that involves small-capacity vehicles delivering cargoes from satellites to customers. When introducing the 2E-TVRP, we also consider the background that the Ministry of Transport of China (MOTC) has been popularizing the tractor and semitrailer (T&St) transportation since 2010, because the T&St combination is being promoted as more energy-efficient than the single-unit trucks. Hundreds of logistics enterprises are presently encouraged to use the T&St to carry out intercity linehaul. However, the logistics enterprises face a problem: synchronously dispatching the linehaul vehicles (especially tractors) and the delivery vehicles (i.e., single-unit trucks) such that the timeliness of the urban delivery service is ensured. Besides, the Ministry of Commerce of China has issued the “National Plan on E-commerce (2016-2020)” that takes many measures, which includes synchronized dispatching of the linehaul vehicles and the delivery vehicles. Therefore, the 2E-TVRP has good practical background. In the 2E-TVRP, the first echelon involves large-capacity vehicles (i.e., T&St combinations in this study) that transport cargoes among satellites and the second echelon involves small-capacity vehicles (i.e., single-unit trucks in this study) that deliver the cargoes from the satellites to the customers. The 2E-TVRP involves two types of time constraints: the time windows constraints and synchronization constraints. Synchronization constraints are introduced on the meeting of the vehicles at satellites such that once the first-echelon vehicles have arrived, cargoes on the first-echelon vehicles are loaded onto the second-echelon vehicles as soon as possible. In practice, many enterprises use semitrailers as facilities for intermediate storage. Since satellites in the 2E-TVRP operate according to the vehicle synchronization, the 2E-TVRP does not particularly consider the intermediate storage constraints. In this study, we refer to Li et al. (2018), for the storage and transshipment capacity constraints in formulating the two-echelon scheme. The contributions of this study include: We first introduce the 2E-TVRP that uses time windows and synchronization constraints as the key for the mixed integer programming model. The nonlinear element in the objective function is linearized to put forward a formulation that makes it easier to find the optimal solution through some commercial IP solver software (e.g., CPLEX). The objective of the exact methods formulated by CPLEX, which may provide a benchmark for some small-scale test instances, is to find exact solutions to the test heuristics. We then provide a two-stage heuristic algorithm that includes a savings-based heuristic and a variable neighborhood search (VNS) phase. The formulations and the heuristic are tested by using some illustrative examples. The heuristic is then used to solve various realistic instances.

2. The 2E-TVRP formulation 2.1. Problem description As intermediary facilities on networks with two echelons, the satellites are not the origin of the

3

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

cargoes; the cargoes are transferred from one satellite to another (Shen and Honda, 2009). There is usually inter-satellite shipping where large-capacity vehicles transport cargoes between satellites located outside various cities: each satellite simultaneously acts as both the origin and the destination of the cargoes. The 2E-TVRP considers inter-satellite linehaul on the first echelon and delivery from satellites to customers on the second echelon. The 2E-TVRP has the following features. i) There are two types of vertices: satellites among which one serves as the central depot, and the customers served by satellites. For each satellite, any other satellite may be the cargo source. The inter-satellite shipping is full truckload (FTL) while the urban delivery is less-than-truckload (LTL). ii) Single-unit trucks and T&St combinations are respectively used for delivery and inter-satellite linehaul operations. A fleet of homogenous tractors is available at the central depot, and each satellite holds a fleet of homogenous trucks. The route duration constraints for tractors and trucks are respected. iii) On the loaded-semitrailer flow network of the first-echelon, all tractors are at the central depot in the beginning, and the tractor routes are closed. A tractor can run either with one loaded-semitrailer or alone. iv) On the second-echelon network, direct shipments from the central depot to its assigned customers are allowed; and each satellite is associated with a handling time for the unloading/loading operations. v) To manage the routing dependence between the two echelons, the 2E-TVRP makes use of the difference between the arrival time of the loaded-semitrailers and the departure time of the trucks at a satellite: At a satellite, the difference between the arrival time of the cargoes on the loaded-semitrailers and the departure time of the cargoes on the trucks is considered as cargo waiting time. It is encouraged to have the cargoes leave the satellites as soon as possible. The penalty cost of cargo waiting is divided into two parts by the earliest departing time of the used trucks. A two-echelon network with four satellites (one as the central depot) is illustrated in Figure 1. The colored squares denote the central depot (i.e., CD) and the remaining satellites (i.e., S1, S2, S3). The color of a shadow square denotes the type of cargo sources. Colored dots and serial numbers denote customers. The color of a shadow dot denotes the type of needed cargoes. The bold lines on the first echelon denote a tractor running with a loaded-semitrailer and the normal lines denote a tractor running alone. Furthermore, at is the arrival time of the loaded-semitrailers and Mindt is the earliest departure time of the trucks at a satellite. The planning period is 48 hours, and we use the time section [-24, 24]. The necessary time for cargo transshipment at satellites is assumed to be one hour. The feasible solution consists of two routes on the first echelon: CD→S2→S1→S2→S3→S2→CD, CD→S1→S3→S1→CD. The routes on the second echelon are as follows: i) CD→5→10→9→CD, CD→6→CD, CD→7→12→CD, CD→11→8→14→CD, CD→13→CD; ii) S1→16→20→S1, S1→17→S1, S1→18→S1, S1→19→S1, S1→21→15→S1, S1→22→S1, S1→23→24→S1; iii) S2→27→S2, S2→28→31→29→S2, S2→32→30→S2, S2→33→25→S2, S2→34→26→S2; iv) S3→35→44→S3, S3→37→38→40→S3, S3→39→S3, S3→41→42→S3, S3→43→36→S3.

4

M AN U

Figure 1. Example of a 2E-TVRP (with a feasible solution)

SC

RI PT

ACCEPTED MANUSCRIPT

2.2. Formulation

In this study, we aim to introduce a 2E-TVRP formulation that can be directly solved by CPLEX12.4. The definitions of the input parameters and variables are summarized in Table 1. We present the 2E-TVRP on a directed complete graph that includes the vertex set that consists of the satellite set (denoted as VD ) and the customer set (denoted as VC ), and the arc set that

TE D

includes the first-echelon arc set ( A1 ) and second-echelon arc set ( A2 ). On the first echelon, a tractor may visit the same satellite more than once. A set of dummy vertices of each satellite is added to ensure that a tractor visits a satellite vertex at most once. We use qi to denote the cargo

EP

volumes required by customer i, and Qij to denote the cargo flow between satellites i and j; Qij

AC C

is decided by qi , which is assumed to be known in advance. We use rij to denote the number of loaded-semitrailers flowing from satellite i to j; rij is estimated by Qij and the capacity and load factor of each semitrailer. Each customer or satellite has a time window that is denoted by the earliest and latest service-starting time. The service time of each customer, which is assumed to be known, is also considered. Table 1. Parameters and variables Sets

V = VD ∪ VC

V

Vertex set,

VD

Set of satellites,

VC

Set of customers,

VCm

Set of customers served by satellite m

v0 ∈VD

is also known as the central depot

VC = Um∈VD VCm

5

ACCEPTED MANUSCRIPT A = A1 ∪ A2

A

Arc set,

A1

First-echelon arc set,

A2

Second-echelon arc set,

A2m

Set of arcs in the service area of the satellite m,

A1 = {(i, j) | i, j ∈VD , i ≠ j}

A2 = Um∈VD A2m A2m = {(i, j ) | i, j ∈VCm ∪{m}, i ≠ j}

Set of dummy vertices of satellite m

m VDD VD′

Vertex set of the first echelon,

m VD′ = VD ∪ ( U VDD )

A1′

Arc set of the first echelon, A1′ = {(i,

RI PT

m∈VD

j) | i, j ∈VD′ , i ≠ j}

Input parameters Cargoes required by customer i. q i

Cargo flow between satellite i and j

Qst

δ1

Qt

& &

δ2

Capacities of each semitrailer and each single-unit truck, respectively

SC

Qij

Load factors of each semitrailer and each single-unit truck, respectively

rij

Number of loaded-semitrailers flowing from satellite i to j

[eim , lim ]

Time windows: If i ∈ VC , ei

m

and li

represent the earliest and latest service-starting time

of customer i. If i = m , e and l represent the earliest and latest opening time of satellite m Serving time of customer i served by satellite m m i

sim

m

M AN U

m

m i

(i, j ) ∈ A1′ , supposing d ij = d ji

d ij

Distance of arc

d ijm

Distance of arc (i, j ) ∈ A2 , supposing d ij = d ji

vl & v t

Velocity of a tractor running with a loaded-semitrailer or a tractor running alone, respectively

v

Velocity of a single-unit truck

m

m

(i, j ) ∈ A1′ tijl = dij vl

Running time of a tractor with a loaded-semitrailer covering the arc

TE D

t

l ij

m

(i, j ) ∈ A1′ , tijt = dij v t

tijt

Running time of a tractor covering alone the arc

tijm

Running time of a single-unit truck covering the arc (i, j ) ∈ A2 ,

cijl

Running cost of a tractor with a loaded-semitrailer covering the arc

η1

EP

where t ij

c

m

tijm = dijm v (i, j ) ∈ A1′ .

is the fuel consumption coefficient of a tractor pulling on a loaded-semitrailer

Running cost of a tractor covering alone the arc

(i, j ) ∈ A1′ . cijt = η2 dij

where

η2

is the fuel

consumption coefficient of a tractor running alone Running cost of a single-unit truck covering the arc (i, j ) ∈ A2 . cij = η3 d ij where

AC C

cijm

H Km maxT1

m

η3

is

Route duration for a tractor Route duration for a truck

αm β1 & β 2

Handling cost of a unit of cargo at satellite m

M

m

the fuel consumption coefficient of a single-unit truck running Set of tractors held by the central depot Set of trucks held by satellite m

maxT2 P Tope

γ

m

Duration of the planning period Time for cargo transshipment at satellites

Waiting cost of a loaded-semitrailer for a unit of time at satellites, and of a unit of cargo at satellites after the earliest departing time of trucks, respectively Waiting cost of a truck for a unit of time at customers A large enough positive integer

6

ACCEPTED MANUSCRIPT Variables

dt h

Departure time of the h-th tractor ( h ∈ H )

athi

Arrival time of the h-th tractor at satellite i ( i ∈VD′ )

wrhi

Waiting time of the loaded-semitrailer pulled-on by the h-th tractor at satellite i ( i ∈VD′ ). If the h-th tractor does not arrive at satellite i or arrives alone at satellite i, wrhi = 0

dt

Departure time of the k-th truck ( k ∈ K ). It is assumed that dt1 ≤ dt2 ≤ ... ≤ dtk

m k

m

m

m

m

atkim

Arrival time of the k-th truck at customer i ( i ∈ VC )

wtkim

Waiting time of the k-th truck at customer i ( i ∈ VC ) The binary decision-variable, which is 1 only if the h-th tractor with a loaded-semitrailer covers m

l xijh

(i, j ) ∈ A1′

the arc

RI PT

m

(i, j ) ∈ A1′

t xijh

The binary decision-variable, which is 1 only if the h-th tractor alone covers the arc

m yijk

The binary decision-variable, which is 1 only if the k-th truck covers the arc (i, j ) ∈ A2

ξh

The binary variable, which is 1 only if the h-th tractor is used

ωkm

The binary variable, which is 1 only if the k-th truck is used

M AN U

SC

m

The vehicle flow formulation for the 2E-TVRP is as follows: The objective function is defined as equation(1).

∑ (β

m∈VD k∈K m

2

⋅ (dtkm − dt1m )⋅

∑q ⋅ ∑ j

j∈VCm

TE D

+∑

yijkm )

(1)

i∈{vm ∪VCm }

EP

The objective function includes six parts: the tractor running cost on the first echelon; the cargo handling cost at the satellites; the penalty cost of loaded-semitrailers waiting at satellites (before the earliest departure time of single-unit trucks); the truck running cost on the second echelon; the penalty cost of single-unit trucks waiting at customers; and the penalty cost of cargoes on single-unit trucks waiting at satellites (after the earliest departure time of single-unit trucks). All parts of the objective function are estimated by the fuel consumption (unit: L).

AC C

i) Constraints for the second echelon:

∑y

m mjk

= ωkm ,

∀m ∈ VD , ∀k ∈ K m

(2)

j∈VCm

m y mjk =

∑

j∈VCm

∑

∑∑y k ∈K

i∈VCm

∑

(4)

∀m ∈ VD

(5)

m

∑ ∑ k∈ K

(3)

j∈VCm

m ymjk ≤| K m |,

∑∑ m

∀ m ∈ VD , ∀ k ∈ K m

≥ ( ∑ q j ) / (Qt ⋅ δ 2 ), ∀m ∈ VD

m mjk

j∈VCm k ∈K m

j∈VCm

y mjmk ,

j∈VCm

m yijk = 1,

∀ m ∈ VD , ∀ j ∈ VCm

(6)

∪{ m }

m yijk ≤ ω km ,

∀ m ∈ VD , ∀ j ∈ VCm , ∀ k ∈ K m

(7)

∀ m ∈ VD , ∀ j ∈ VCm , ∀k ∈ K m

(8)

i∈VCm ∪{ m }

∑ i∈VCm ∪{ m }

m yijk =

∑

y mjik ,

i∈VCm ∪{ m}

7

ACCEPTED MANUSCRIPT m yijk ⋅ qi ≤ Q t ⋅ δ 2 ,

∀ m ∈ VD , ∀ k ∈ K m

(9)

m dtkm + tmjm − M ⋅ (1 − ymjk ) ≤ atkjm ,

∀m ∈VD , ∀j ∈VCm , ∀k ∈ K m

(10)

m m dtkm + tmj +M ⋅ (1 − ymjk ) ≥ atkjm ,

∀m ∈ VD , ∀j ∈ VCm , ∀k ∈ K m

(11)

∑ ∑ i∈VCm

j∈VCm ∪{ m }

atkim + wtkim + sim + tijm − M ⋅ (1 − yijkm ) ≤ atkjm , ∀m ∈VD , ∀i ∈VCm , ∀j ∈VCm ∪ {m}, ∀k ∈ K m

(12)

at + wt + s + t + M ⋅ (1 − y ) ≥ at , m ki

m i

m ij

m ijk

m kj

RI PT

m ki

∀m ∈ VD , ∀i ∈ VCm , ∀j ∈ VCm ∪ {m}, ∀k ∈ K m

(13)

dt ≥ e ,

∀m ∈ VD , k ∈ K

m

(14)

dt ≤ l ,

∀m ∈VD , k ∈ K

m

m m

m k

m m

e ⋅

y ≤ at + wt ≤ l ⋅

∑

m i

m ijk

m ki

m ki

∑

m i

j∈VCm ∪{m}

(15)

m ijk

y ,

j∈VCm ∪{ m}

∀m ∈ VD , ∀i ∈ VCm , ∀k ∈ K m

∑

∑

y ⋅t + ∑ s ⋅ m ijk

m ij

i∈VCm ∪{m} j∈VCm ∪{m}

i∈VCm

(16)

y + ∑ wt ≤ maxT2 ,

∑

m i

SC

m k

m ijk

m ki

j∈VCm ∪{m}

i∈VCm

∀m ∈ VD , k ∈ K (17) Constraints (2) and (3) guarantee that if the k-th truck held by satellite m is used, its route should be closed. Constraints (4) and (5) initialize the number of trucks held by satellite m. Constraints (6)-(8) ensure that each customer is visited once by the assigned truck. Constraint (9) restricts the capacity of the used trucks. The relationship between the departure time and the arrival time of the trucks at customers is indicated by constraints (10)-(13). Constraints (14) and (15) ensure that the trucks depart in the opening time of satellite m. Constraint (16) indicates that the trucks respect the time windows of customers; i.e., if the truck arrives at the customer before the time window opening, it should wait. Constraint (17) respects the route duration of the trucks.

TE D

M AN U

m

ii) Constraints for the routing interaction:

m ∀i ∈ VD′ , ∀j ∈ {m} ∪ VDD , ∀m ∈ VD , ∀h ∈ H

(18)

wrhj ≥ dt − athj − M ⋅ (1 − x ) m 1

l ijh

AC C

EP

m (19) ∀i ∈ VD′ , ∀j ∈ {m} ∪ VDD , ∀m ∈ VD , ∀h ∈ H As the key of the 2E-TVRP formulation, constraints (18) and (19) are used to decide routes on each of the two interacting echelons, which guarantee that the loaded-semitrailers arrive at the satellite before the trucks depart, while ensuring that the necessary handling time for cargo transshipments is respected.

iii) Constraints for the first echelon:

∑ (x

l 0 jh

j∈VD′

+ x0t jh ) =

(20)

∑ (x

j∈VD′

l j 0h

+ xtj 0 h ),

∀h ∈ H

(21) (22) (23)

∑ (x

i∈VD′

l ijh

t + xijh )=

∑ (x

i∈VD′

l jih

+ xtjih ),

∀h ∈ H , ∀j ∈ VD′ \ {v0 }

8

(24)

ACCEPTED MANUSCRIPT

∑ ∑

h∈H

l xijh = rm1m2 ,

∑

m1 i∈{m1}∪VDD

m2 j∈{m2 }∪VDD

∀m1 ∈ VD , ∀m2 ∈ VD

(25)

∀h ∈ H , ∀j ∈VD′ \{v0 }

(26)

dth + t0l j + M ⋅ (1 − x0l jh ) ≥ athj ,

∀h ∈ H , ∀j ∈VD′ \{v0 }

(27)

dth + t0t j − M ⋅ (1 − x0t jh ) ≤ athj ,

∀h ∈ H , ∀j ∈VD′ \{v0 }

(28)

dth + t0t j + M ⋅ (1 − x0t jh ) ≥ athj ,

∀h ∈ H , ∀j ∈VD′ \{v0}

(29)

athi + t − M ⋅ (1 − x ) ≤ athj ,

∀h ∈ H , ∀i ∈VD′ \{v0 }, ∀j ∈VD′

(30)

l athi + tijl + M ⋅ (1 − xijh ) ≥ athj ,

∀h ∈ H , ∀i ∈VD′ \{v0 }, ∀j ∈VD′

(31)

l ij

l ijh

RI PT

dth + t0l j − M ⋅ (1 − x0l jh ) ≤ athj ,

(32)

∑ ∑ (x

l ijh

i∈VD′ j∈VD′

⋅ t + x ⋅ t ) ≤ maxT1 , l ij

t ijh

∀i ∈ VD′ , ∀j ∈ VD′ , ∀h ∈ H

l t xijh + xijh ≤ 1,

+∑x i∈VD′

t jih

∀j ∈ VD′ \ {v0 }, ∀h ∈ H

≤ 1,

(33)

(34)

(35) (36)

M AN U

∑x

i∈VD′

t ijh

∀h ∈ H

t ij

SC

∀h ∈ H , ∀i ∈VD′ \{v0 }, ∀j ∈VD′

t athi + tijt + M ⋅ (1 − xijh ) ≥ athj ,

TE D

Constraints (20) and (21) ensure that each used tractor departs from and finally returns to the central depot. Constraint (22) guarantees that each satellite vertex (including the dummies) is visited at most once. Constraints (23) and (24) guarantee that a used tractor may visit a satellite vertex at most once. Constraint (25) ensures that the loaded-semitrailer demand is served. Constraints (26)-(33), which indicate the relationship between the departure time and the arrival time of the tractors at satellites, ensure the continuity of routes. Constraint (34) restricts the route duration of the tractors. Constraint (35) ensures that a used tractor may cover an arc either by running with a loaded-semitrailer or by running alone. Constraint (36) ensures that a used tractor cannot cover two adjacent arcs alone.

l xijh ∈{0,1},

ω ∈{0,1}, m k

x = x = 0, l ijh

t ijh

(37) (38)

∀i, j ∈VD′ , ∀h ∈ H

(39)

∀i, j ∈VD′ , ∀h ∈ H

(40)

AC C

t xijh ∈{0,1},

EP

iv) Some auxiliary constraints: P P − ⋅ ξh ≤ dth ≤ ⋅ ξh , ∀h ∈ H 2 2

∀m ∈VD , ∀k ∈ K

(41) m

(42) (43) (44)

∀m ∈ VD , i, j ∈ {m} ∪ V , ∀h ∈ H m DD

(45)

Constraints (37) and (38) restrict the planning period, which is assumed to be P hours. Constraints (39)-(43) specify the domain of some variables. Constraints (44) and (45) initialize some variables.

2.3. Reformulation In Section 2.2 we introduced a mixed integer nonlinear programming model for the 2E-TVRP. The commercial IP solver software (i.e., CPLEX12.4 used in this study) probably fails to provide a benchmark for the small-scale test instances. To overcome the computational difficulties, the 9

ACCEPTED MANUSCRIPT following linearization method that equivalently transforms the mixed integer nonlinear programming model into a mixed integer linear programming (MILP) model that can be solved directly by CPLEX12.4 is proposed. The nonlinear element in the 2E-TVRP formulation is involved by the following objective function:

∑ ∑ (β

2

⋅ (dtkm − dt1m )⋅

m∈VD k∈K m

∑q j∈VCm

j

⋅

∑

yijkm )

i∈{vm ∪VCm }

RI PT

A new variable ( ct kim ) is used to denote the departure time of the k-th truck at satellite m where

M AN U

SC

the k-th truck serves a customer set including customer i. The new variable makes it easier to calculate the departure time difference between the earliest departing truck and other trucks used at a satellite. The calculated time difference corresponds to the waiting time of the cargoes. The penalty cost of cargoes on single-unit trucks waiting at satellites (after the earliest departure of single-unit trucks) is equivalently transformed into the penalty cost of cargoes of each customer waiting at satellites. Constraints including the new variables are as follows: (L-1)

(L-2)

TE D

Constraints (L-1) and (L-2) ensure that the departure time of the k-th truck is equal to that of the cargoes of customers served by the k-th truck. The nonlinear element in the objective function is rewritten as follows:

EP

Besides, the following constraint that ensures that all kinds of departure times of cargoes are not earlier than the departure time of the earliest departing truck is added:

ctkim ≥ dt1m ,

∀m ∈ VD , ∀i ∈ VCm , ∀k ∈ K m

（L-3）

AC C

The linearization provides the MILP formulation that can be directly solved by CPLEX12.4. However, compared to the mixed integer nonlinear programming model, the MILP formulation introduces extra variables and constraints. For instance, if we consider including four satellites (among which one is used as the central depot) and 40 customers as an example, the mixed integer nonlinear programming model involves 10128 variables and 38231 constraints while the MILP model involves 10528 variables and 39431 constraints. If the last part of the objective function in Section 2.2 is abandoned, i.e., if the penalty cost of cargoes on single-unit trucks waiting at the satellite is ignored, CPLEX12.4 can solve the formulation. We denote the 2E-TVRP formulation ignoring the penalty cost of cargoes on single-unit trucks waiting at satellites in the objective function as the following model (called BM model). (BM model) The objective function is defined as (46).

10

ACCEPTED MANUSCRIPT

(46) Subject to constraints (2)-(45).

2.4. Illustrative examples

M AN U

SC

RI PT

Fifteen small example networks including 3, 4, or 5 satellites (among which one is used as the central depot) are randomly generated to illustrate the performance of the proposed vehicle flow formulation for the 2E-TVRP. It is assumed that the number of customers served by each satellite is the same, and the cargo demand of all customers served by one satellite is from the same origin satellite. The time windows of customers are designed by the method proposed by Solomon (1987). An example is denoted by “RAND n-num_C-SL_type_num-i”, where n is the number of satellites, num_C is the number of customers served by a satellite, SL_type is the number of cargo origins of a satellite, and i is the example sequence in the same n, num_C, and SL_type_num. In our computational experiments, a type of T&St combination and a single-unit truck that satisfy the fuel-efficiency requirements of the “Regulation of Supervising Vehicle Fuel Consumption” (No. 11/2009 Decree of the MOTC) are used. The values of the input parameters included in the model are shown in Table 2. Table 2. The values of input parameters included in the computational experiments

Value (unit)

Parameter

Value (unit)

Qst Qt δ1 δ2 vl vt V η1 η2

17 (t) 6 (t) 0.7 0.6 50 (km/h) 50 (km/h) 40 (km/h) 0.29 (L/km) 0.16 (L/km)

η3 maxT1 maxT2 Tope αm β1 β2 γ P

0.14 (L/km) 16 (h) 10 (h) 1 (h) 0.4 (L/t) 2.0 (L/h)

EP

TE D

Parameter

0.17 (L/ t ⋅ h ) 2.75 (L/h) 48 (h)

Note: The MOTC issues (on the website (http://atestsc.mot.gov.cn/pub/publish/)) hundreds of vehicle types that

AC C

satisfy the No. 11/2009 Decree. Once the vehicle type is selected, the fuel consumption coefficient and the velocity are obtained from the information on the website.

The MILP and the BM models are solved directly by CPLEX12.4 on several computers with the same configurations (i.e., an Intel(R) Core(TM) 3.2 GHz processor with 8 GB memory using Windows XP). For each example, CPLEX12.4 runs with default settings until it finds the optimal solution. The computational results are shown in Table 3: Column 1 indicates the examples; Columns 2-3 show the results of the objective function (Obj) of the MILP and BM models; Columns 4-5 show the computation time (seconds) of the MILP and BM formulations, respectively; and Column 6 compares the routes in the solutions provided by the MILP and BM formulations. We observe that for eight of the fifteen examples, CPLEX12.4 can provide the optimal solutions.

11

ACCEPTED MANUSCRIPT Table 3. Results of illustrative examples solved optimally by CPLEX12.4 Example RAND3-5-1-1

Obj

Time (s)

Change in the routes of MILP solution with

MILP

BM

MILP

BM

respect to those of BM solution

448.2

439.3

328

45

Departure time of one second-echelon route is advanced

RAND3-5-1-2

449.7

436.9

273

13

For two satellites, the number of delivery routes decreases by 1

428.8

421.8

85

51

The departure time of one first-echelon route is postponed;

RI PT

RAND3-5-1-3

the

departure

time

of

two

second-echelon routes is postponed; and the departure time of one second-echelon route is advanced 379.7

368.3

5896

26

For one satellite, the number of delivery routes

SC

RAND3-5-1-4

decreases by 1 RAND3-5-1-5

381.2

376.2

30645

1602

The departure time of one second-echelon route is postponed; and the departure time of

M AN U

three second-echelon routes is advanced

RAND4-5-1-1

544.3

531.1

16288

233

No change

RAND4-5-1-2

491.6

480.4

11119

699

The departure time of one first-echelon route is postponed;

the

departure

time

of

one

second-echelon route is postponed; and the departure time of two second-echelon routes is

Advanced

511.4

493.2

34514

205

TE D

RAND4-5-1-4

The departure time of one second-echelon route is advanced

Note: CPLEX12.4 cannot output the optimal solutions for the randomly-generated examples: “RAND4-5-1-3”, “RAND4-5-1-5”, and “RAND5-5-1-i” (i = 1, 2, 3, 4, 5).

AC C

EP

The following observations are made from Table 3: i) From the angle of global optimization of the two-echelon logistics systems, the truck departure times at the satellites are different, due to which the waiting time of cargoes of each customer at the satellites is different. Including the penalty cost of the cargoes on single-unit trucks waiting at satellites leads to larger objective values. ii) The MILP formulation involves more variables and constraints than the BM model. Therefore, the computation time of CPLEX12.4 for the MILP formulation is obviously longer than that for the BM formulation. Moreover, for the same number of satellites and customers, the computation time of CPLEX12.4 for each example is different. The location of the vertices and time windows are important factors that influence the computation time. iii) The penalty cost of cargoes on single-unit trucks waiting at satellites involved by the objective may have effects such as the following: the departure time of single-unit trucks is either advanced or postponed; the arrival time of the tractors is postponed; and the customers included in the truck routes are changed. iv) CPLEX 12.4 could only provide the optimal solutions for eight of the fifteen illustrative examples. Therefore, it is necessary to design heuristic algorithms that can solve the 2E-TVRP 12

ACCEPTED MANUSCRIPT formulation to attain suboptimal solutions of realistic-size instances.

3. A two-stage heuristic

RI PT

A two-stage heuristic algorithm is proposed in this study to solve the 2E-TVRP formulation. The first stage aims to construct an initial solution. Through the CW algorithm, the second-echelon vehicle routes are constructed first and then the first-echelon vehicle routes are constructed. The second stage aims to improve the initial solution constructed by the CW algorithm through the VNS algorithm. To construct the initial solution, the second-echelon routes for each satellite are constructed by the CW on condition that constraints (2)-(17) are respected. Furthermore, the arrival time of the first-echelon vehicles for each satellite is restricted by constraints (18) and (19); and the first-echelon routes are constructed through the CW by respecting constraints (20)-(36).

3.1. The initial solution

TE D

M AN U

SC

3.1.1. Second-echelon routes The routing problem for each satellite in the 2E-TVRP is the VRP with time windows (VRPTW). To construct a new truck route, a customer with the minimum of latest-service-starting time of customers is inserted into the route, and the customer with the maximum savings is inserted step-wise. The savings, which are expressed as time savings, include distance savings and waiting-time savings. 3.1.2. First-echelon routes Once the second-echelon routes of the u-th satellite are initialized, the earliest departure time of trucks (denoted as E_DT(u)) at the u-th satellite is restricted. The loaded-semitrailers that carry cargoes for the earliest truck should arrive at the u-th satellite before E_DT(u)－Tope (where Tope is the time required for cargo transshipments). As mentioned in Section 2, the T&St combination of vehicles is chosen for the inter-satellite linehaul. The first-echelon VRP is abstracted on a loaded-semitrailer flow network. A set of FTLs, which are actually origin-destination (O-D) pairs with loaded-semitrailer demand, waits for transport. In the beginning, all tractors remain at the central depot. The number of loaded-semitrailers arriving at a satellite is determined by the delivery demand of the satellite. The loaded-semitrailers of the same O-D pair are denoted as order. Table 4 shows the information of the order.

EP

Table 4. The information of the Order

O (start)

D (end)

num

1 2 ⁞ m

i1 i2 ⁞ im

j1 j2 ⁞ jm

ri1j1 ri2j2 ⁞ rimjm

AC C

Order index

The initial solution is made up of direct trips. In a direct trip, a tractor departs from the central depot to satellite i, pulls on a loaded-semitrailer from satellite i to satellite j, drops off the loaded-semitrailer at satellite j and then returns directly to the central depot. To construct a new tractor route, an order with the minimum latest-arriving-time of a loaded-semitrailer is selected from the order set. Then, a new trip is merged by referring to the maximum savings. The details are obtained from Li et al. (2015).

3.2. The VNS The VNS, which uses the systematic change of neighborhoods both within a descent phase to find a local optimum and in a perturbation phase to get out of the corresponding valley, is a simple and effective metaheuristic (Mladenović and Hansen, 1997). The VNS systematically exploits the 13

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

idea of neighborhood change, both in descent to local minima and in escape from the valleys that contain them (Hansen et al., 2006). By increasingly exploring the distant neighborhoods of the current incumbent solution, a new solution may be obtained only if there is an improvement. Thus, most variables are at their optimal values and this good character of the optimal solution will be retained and used to obtain promising neighboring solutions. Moreover, a local search is repeatedly applied on the neighboring solutions to obtain a local optimum. The VNS scheme is used to solve the 2E-TVRP. The local search operators adopt the first improvement, and the stop condition adopts a maximum number of iterations. For each satellite, the routing problem for the second-echelon is actually the VRPTW; similarly, for the central depot, the routing problem for the first-echelon is the VRPTW. There are two types of operators called relocate and λ-interchange that are generally used in the literature for optimizing the VRPTW. In this study, the λ-interchange operator has been used for searching neighborhood solutions. The implementation strategy for neighborhood operators is selected through computational trials. The process of the heuristic based on the VNS is shown in Figure 2. The initial solution constructed by the CW is denoted as best_S. Besides, best_S is renewed in course of the VNS computation. The iterations of the complete VNS phase are denoted as Global_iteration. The iterations of the VNS for the second echelon and those for the first echelon are denoted as sl_iteration and fl_iteration, respectively. During the optimization of the first-echelon routes through the VNS, fMID_S represents the feasible solution obtained in the shaking phase, and fMID_SS represents the optimized solution obtained in the local search phase. Similarly, during the optimization of the second-echelon routes through the VNS, sMID_S represents the feasible solution obtained in the shaking phase, and sMID_SS represents the optimized solution obtained in the local search phase; and Obj(s) is the objective value of the attained solution s.

14

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 2. VNS algorithm for the 2E-TVRP

AC C

3.2.1 Shaking operators The shaking operators are the stochastic components of the VNS, which can select neighbors randomly. It is difficult to balance the computation effectiveness and the chance of getting out of local optima. Neighborhood structures should allow to sufficiently perturb the incumbent solution and ensure the new solution while retaining the important parts of the incumbent. In the literature, the λ-interchange is one of the best neighborhood structures for optimizing the VRPTW. We consider three types of λ-interchange, which are denoted as (λ1, λ2)-interchange, as the shaking operators. Two routes are randomly selected; and then, λ1 continuous segments of one route and λ2 continuous segments of the other route are selected for exchange. The study adopts (1,0)-interchange, (1,1)-interchange, and (1,2)-interchange. Figure 3 shows several examples of the λ-interchange. Figure 3(a) shows an example of the (1,0)-interchange operator: Two routes (i.e., “0-1-2-3-4-0” and “0-5-6-7-8-0”) are selected. A vertex or an order of the former route is removed and inserted into the latter route. Figure 3(b) shows an example of the (1,1)-interchange

15

ACCEPTED MANUSCRIPT

M AN U

SC

RI PT

operator. A vertex or an order of the former route and a vertex or an order of the latter route are exchanged. Figure 3(c) shows an example of the (1,2)-interchange operator. A vertex or an order of the former route and two continuous segments of the latter route are exchanged.

Figure 3. Examples of the λ-interchange operators

EP

TE D

3.2.2 Local search operators The local search operators are the deterministic components of the VNS, which starts immediately after a new solution is found through shaking operators. We used five types of neighborhood structures as local search operators: (1,0)-interchange, (2,0)-interchange, (1,1)-interchange, (1,2)-interchange, and (2,2)-interchange. Examples for the five neighborhood structures are shown in Figure 3. Figure 3(d) shows an example of the (2,2)-interchange operator: Two continuous segments of the former route and two continuous segments of the latter route are exchanged. Figure 3(e) shows an example of the (2,0)-interchange operator: Two continuous segments of the former route are removed and inserted into the latter route.

3.3. Computational results of the illustrative examples

AC C

The two-stage heuristic is used to solve the eight illustrative examples in Section 2.4. The performance results of the heuristic are summarized in Table 5. The results of the heuristic and those of the MILP solved by CPLEX are also compared: Column 1 indicates the examples; Columns 2-3 show the results of the objective function (Obj), i.e., the results of the heuristic (HA) and the MILP formulations, respectively; Column 4 shows the percentage gap between the results of the heuristic and MILP; and Columns 5-6 show the computation time of the heuristic and MILP, respectively. Table 5. Results of illustrative examples solved by the heuristic Example

Obj

Time (s)

HA

MILP

Gap (%)

HA

MILP

RAND3-5-1-1

450.9

448.2

0.6

29.5

328

RAND3-5-1-2

451.6

449.7

0.4

27.0

273

RAND3-5-1-3

433.4

428.8

1.1

29.3

85

16

ACCEPTED MANUSCRIPT RAND3-5-1-4

388.4

379.7

2.3

28.0

5896

RAND3-5-1-5

389.7

381.2

2.2

50.3

30645

RAND4-5-1-1

550.3

544.3

1.1

95.3

16288

RAND4-5-1-2

502.2

491.6

2.2

65.2

11119

RAND4-5-1-4

515.3

511.4

0.8

101.5

34514

SC

RI PT

The heuristic proves to be a general algorithmic structure and can produce the results in less than two minutes for the eight examples. The experimental results of the eight randomly generated examples show that the two-stage heuristic performs well in terms of the computation time and solution quality. Of the eight examples for which optimal solutions have been provided by CPLEX12.4, the percentage gap between the results of the heuristic and that of CPLEX12.4 is less than 2.3%, and the average percentage gap between the results of the heuristic and that of CPLEX12.4 equals 1.3%. Therefore, the experimental results prove that the proposed two-stage heuristic can effectively solve randomly generated small-scale instances, and the heuristic solutions are close to the optimal results.

M AN U

4. Case study 4.1. The realistic instances

AC C

EP

TE D

The purpose of the case study on realistic instances is to evaluate the applicability of the introduced model and the two-stage heuristic for solving realistic-size problems. A trucking company, simply named SDHLG, is the object of our case study. SDHLG has 23 satellites distributed in 23 prefecture-level or county-level cities in the Shandong Province of China. SDHLG uses hundreds of single-unit trucks for from-linehaul-to-delivery service. The transport network of SDHLG includes two echelons: On the first echelon, large-capacity vehicles are used for intercity linehaul among the 23 satellites. On the second echelon, small-capacity trucks that depart from satellites deliver cargoes to customers. In addition to adhering with the policy on popularizing T&St combinations issued by the MOTC, SDHLG is substituting T&St combinations for single-unit trucks on the inter-satellite linehaul echelon. The data on the loaded-semitrailer demand for the inter-satellite echelon and a small proportion of data on customers are obtained. The data on the first echelon (e.g., loaded-semitrailer demand, distances among satellites) is known in advance while the data on the second echelon (e.g., customer locations, time windows) is imperfect. We attempted to supplement the data on the second echelon according to SDHLG’s practice. SDHLG uses the satellite in Jinan City (JNA) as the central depot. In a planning period, 56 tractors are used for serving the inter-satellite linehaul echelon. The case study assumes that any satellite can be regarded as the central depot. There are 23 instances that are classified by the central depot location. Some data on the realistic instances are shown in Table 6 and Figure 4. The case study adopts 23 instances with up to 23 satellites (among which one is the central depot) and 1008 customers. Table 6. Some data on the realistic instances Central

Arrived

Loaded-semitrailer’s

Customers of the

depot*

loaded-semitrailer

origins

satellite

PI-1

JNA

21

14

272

PI-2

QD

8

5

97

PI-3

HD

1

1

15

Instance

17

ACCEPTED MANUSCRIPT ZB

4

2

55

PI-5

ZZH

2

2

23

PI-6

TZH

1

1

15

PI-7

DY

1

1

11

PI-8

GR

1

1

11

PI-9

YT

3

2

37

PI-10

LZH

1

1

13 36

LYA

3

2

LK

1

1

PI-13

WF

21

11

PI-14

JNI

2

2

PI-15

TA

1

1

PI-16

WH

2

2

PI-17

RZH

1

PI-18

BZH

1

PI-19

DZH

1

PI-20

LCH

1

PI-21

LYI

PI-22

HZ

PI-23

LW

RI PT

PI-11 PI-12

SC

PI-4

11

265 26

14

25 11

1

12

1

13

1

15

1

1

10

1

1

10

1

1

11

M AN U

1

Note: * For an instance, one satellite is regarded as the central depot, and the central depots of all other instances

AC C

EP

TE D

are regarded as satellites.

Figure 4. Location of satellites and some customers on SDHLG two-echelon logistics network

4.2. Computational results The parameter values presented in Table 2 are adopted in the realistic instances. The heuristic results for the 23 instances are shown in Table 7. The values of the six parts of the objective function are also recorded. We denote the tractor running cost on the first echelon as Obj1, the penalty cost of loaded-semitrailers waiting at satellites as Obj2, the penalty cost of cargoes on single-unit trucks waiting at satellites as Obj3, the handling cost at satellite as Obj4, the truck running cost on the second echelon as Obj5, and the penalty cost of single-unit trucks waiting at

18

ACCEPTED MANUSCRIPT customers as Obj6. In Table 7, Column 1 indicates the realistic instance, Columns 2-7 indicate the values of the six parts of the objective function, Column 8 shows the objective values, Columns 9-10 show the number of used tractors (TaN) and trucks (TkN), respectively, and Column 11 shows the computation time (in seconds). The computation time of the heuristic for any realistic instance is less than 100 minutes. This shows that the two-stage heuristic can efficiently solve realistic-size instances.

RI PT

Table 7. Computational results of the realistic instances Obj1

Obj2

Obj3

Obj4

Obj5

Obj6

Obj

TaN

TkN

Time (s)

PI-1

6282

928

808

381

4414

186

12999

28

319

5756

PI-2

7145

728

807

381

4457

174

13692

32

317

5706

PI-3

7840

545

799

381

4473

156

14194

37

318

5608

PI-4

6365

792

800

381

4466

163

12967

29

316

5934

PI-5

11188

405

802

381

4451

152

17379

53

318

5870

PI-6

10222

464

808

381

4441

149

16465

48

318

5825

PI-7

7490

672

798

381

4484

159

13984

34

317

5804

PI-8

6893

704

820

381

4398

160

13356

33

317

5918

PI-9

10494

482

805

381

4438

158

16758

49

317

5632

PI-10

7762

601

808

381

4411

164

14127

37

318

5810

PI-11

8171

557

799

381

4435

163

14506

39

318

5667

PI-12

9346

464

801

381

4433

139

15564

45

318

5721

PI-13

5811

1021

808

381

4406

160

12587

26

320

5652

PI-14

10094

447

796

381

4452

162

16332

47

316

5752

PI-15

7308

801

797

381

4419

167

13873

33

315

5716

PI-16

13471

324

806

381

4445

160

19587

63

317

5885

PI-17

8596

458

801

381

4430

189

14855

41

318

5837

PI-18

7385

645

800

381

4433

151

13795

35

317

5652

PI-19

8611

688

800

381

4447

154

15081

38

319

5621

PI-20

8517

583

805

381

4399

181

14866

39

315

5692

PI-21

8959

458

802

381

4429

159

15188

44

316

5538

M AN U

TE D

12406

375

798

381

4409

155

18524

55

318

5723

6934

706

799

381

4408

181

13409

32

316

5732

AC C

PI-23

EP

PI-22

SC

Instance

Although the case study results may be affected by the evaluated values of input parameters and the company practice, several specific conclusions are summarized as follows: i) The 2E-TVRP focuses on the synchronous routing problems on the intercity linehaul echelon that uses tractor and semitrailer combinations and on the intercity delivery echelon that uses single-unit trucks. The introduced formulation for the 2E-TVRP can provide routes for enterprises to synchronously dispatch the tractors and the single-unit trucks so as to ensure the timeliness of the from-linehaul-to-delivery service. ii) The Obj of the 23 instances varies from 12587 to 19587, and its average is 14960. The number of used tractors of the 23 instances varies from 26 to 63, and the average is 40. The number of used single-unit trucks of the 23 instances does not show an obvious variation. When the satellite with different locations and various loaded-semitrailer flows are respectively regarded as the central depot, the objective value of the heuristic solution obviously varies. Both the central 19

ACCEPTED MANUSCRIPT

5. Conclusion

TE D

M AN U

SC

RI PT

depot location and the loaded-semitrailer flow are important for SDHLG to save costs. Instances where the central depots are located in cities such as JNA, ZB, WF, TA, or LW have relatively good solutions. iii) The average running cost of the tractors is 215.6 (L per tractor), while the maximum is 226.6 (L per tractor) and the minimum is 203.6 (L per tractor). Similarly, the average running cost of single-unit trucks is 14.0 (L per truck), while the maximum is 14.1 (L per truck), and the minimum is 13.8 (L per truck); the average running cost of all vehicles is 36.3 (L per vehicle), while the maximum is 47.1 (L per vehicle), and the minimum is 29.5 (L per vehicle). Moreover, for all the realistic instances, the running cost of trucks is observed to have a slight variation. Therefore, it is advisable to seek more savings from routings of the inter-satellite linehaul echelon. For all the instances, about 57.3% of the total cost is consumed in the running course of tractors while 29.6% of the total cost is consumed in the running course of trucks. Moreover, there exists some relationship between the running cost of tractors and the total cost: using straight lines to fit the relationships, we find y = 0.9200x + 7069.2 (where y denotes Obj and x denotes Obj1) with a correlation coefficient of 0.9989. However, there is no obvious relationship between the running cost of trucks and the total cost. iv) At present, SDHLG uses the satellite in JNA as the central depot. However, the computational results show that the satellite in JNA is not the most satisfactory location for the central depot. The results of the PI-1 and PI-13 instances are more satisfactory, which indicates that using the satellite in WF as the central depot decreases the total cost by 3.2% and the number of used tractors by 7.1%. SDHLG practically uses 56 tractors to carry out the inter-satellite linehaul using the satellite in JNA as the central depot, while the optimized solution from the computational results adopts only 28 tractors. Therefore, even though the 2E-TVRP model has certain assumptions and limitations, there still exists the opportunity for SDHLG to adjust tractor routes to decrease the number of used tractors.

AC C

EP

In this study, we introduced a variant (called the 2E-TVRP) of some practical applicability that has not yet attracted enough research interest. Considering the lack of the assignment decisions between satellites and customers because of localized or too-long distance services, the 2E-TVRP adopts time constraints to make routes on each of the two interacting echelons, instead of using assignment decisions. The time constraints include time windows constraints and synchronization constraints between vehicles of the two echelons: large-capacity vehicles transport cargoes among satellites on the first echelon, and the small-capacity vehicles deliver cargoes from satellites to customers on the second echelon. An MILP model for the 2E-TVRP is constructed and the MILP formulation for the 2E-TVRP is directly solved by CPLEX12.4 to find an optimal solution that can provide a benchmark for some small-scale instances to calibrate and test heuristics. We provided a heuristic algorithm that incorporates savings-based heuristics followed by a VNS phase. Both the vehicle flow formulations and the two-stage heuristic are tested by solving illustrative examples. In the case study, the heuristic effectively solved some realistic instances with up to 23 satellites and 1008 customers. Future research can focus on more effective exact or heuristic algorithms for the 2E-TVRP.

20

ACCEPTED MANUSCRIPT Besides, some variants of the 2E-TVRP that include, for instance, heterogeneous fleet or multi-depots on the first echelon are expected. The MILP model and the practical instances presented in this study serve as baselines for future investigations.

Acknowledgements：This research work was supported by the Research Grant from the

RI PT

National Natural Science Foundation of China (grant number 71672005).

Conflicts of Interest: The authors declare no conflict of interest.

SC

References

Baldacci, R., Mingozzi, A., Roberti, R., & Wolfler, C. R. (2013). An exact algorithm for the two-echelon capacitated vehicle routing problem. Operations Research, 61(2), 298-314.

M AN U

Belgin, O., Karaoglan, I., & Altiparmak, F. (2018). Two-echelon vehicle routing problem with simultaneous pickup and delivery: Mathematical model and heuristic approach. Computers & Industrial Engineering, 115, 1-16. Boccia, M., Crainic, T. G., Sforza, A., & Sterle, C. (2010). A metaheuristic for a two echelon location-routing problem. Lecture Notes in Computer Science, 6049, 288-301.

Breunig, U., Schmid, V., Hartl, R. F., & Vidal, T. (2016). A large neighbourhood based heuristic for two-echelon routing problems. Computers & Operations Research, 76, 208-225.

Contardo, C., Hemmelmayr, V., & Crainic, T. G. (2012). Lower and upper bounds for the two-echelon capacitated location-routing problem. Computers & Operations Research, 39(12), 3185-3199.

TE D

Crainic, T. G., Mancini, S., Perboli, G., & Tadei, R. (2011). Multi-start heuristics for the two-echelon vehicle routing problem. Lecture Notes in Computer Science, 6622, 179-190. Crainic, T. G., Ricciardi, N., & Storchi, G. (2009). Models for evaluating and planning city logistics systems. Transportation Science, 43(4), 432-454.

Cuda, R., Guastaroba, G., & Speranza, M. G. (2015). A survey on two-echelon routing problems. Computers &

EP

Operations Research, 55, 185-199.

Dalfard, V. M., Kaveh, M., & Nosratian, N. E. (2013). Two meta-heuristic algorithms for two-echelon location-routing problem with vehicle fleet capacity and maximum route length constraints. Neural

AC C

Computing and Applications, 23(7-8), 2341-2349.

Diabat, A., & Theodorou, E. (2015). A location-inventory supply chain problem: reformulation and piecewise linearization. Computers & Industrial Engineering, 90, 381-389.

Dondo, R., Méndez, C. A., & Cerdá, J. (2011). The multi-echelon vehicle routing problem with cross docking in supply chain management. Computers & Chemical Engineering, 35(12), 3002-3024.

Grangier, P., Gendreau, M., Lehuédé, F., & Rousseau, L. M. (2016). An adaptive large neighborhood search for the two-echelon multiple-trip vehicle routing problem with satellite synchronization. European Journal of Operational Research, 254(1), 80-91. Hansen, P., Mladenović, N., & Urošević, D. (2006). Variable neighborhood search and local branching. Computers & Operations Research, 33(10), 3034-3045. Hemmelmayr, V. C., Cordeau, J. F., & Crainic, T. G. (2012). An adaptive large neighborhood search heuristic for two-echelon vehicle routing problems arising in city logistics. Computers & Operations Research, 39(12), 3215-3228. 21

ACCEPTED MANUSCRIPT Jepsen, M., Spoorendonk, S., & Ropke, S. (2013). A branch-and-cut algorithm for the symmetric two-echelon capacitated vehicle routing problem. Transportation Science, 47(1), 23-37. Li, H., Lv, T., & Li, Y. (2015). The tractor and semitrailer routing problem with many-to-many demand considering carbon dioxide emissions. Transportation Research Part D, 34, 68-82. Li, H., Zhang, L., Lv, T., & Chang, X. (2016). The two-echelon time-constrained vehicle routing problem in linehaul-delivery systems. Transportation Research Part B, 94, 169-188. Li, H., Liu, Y., Jian, X., & Lu, Y. (2018). The two-echelon distribution system considering the real-time

RI PT

transshipment capacity varying. Transportation Research Part B, 110, 239-260.

Liu, T., Luo, Z., Qin, H., & Lim, A. (2018). A branch-and-cut algorithm for the two-echelon capacitated vehicle routing problem with grouping constraints. European Journal of Operational Research, 266(2), 487-497.

Mladenović, N., & Hansen, P. (1997). Variable neighborhood search. Computers & Operations Research, 24(11), 1097-1100.

SC

Morais, V. W., Mateus, G. R., & Noronha, T. F. (2014). Iterated local search heuristics for the vehicle routing problem with cross-docking. Expert Systems with Applications, 41(16), 7495-7506.

Nguyen, V. P., Prins, C., & Prodhon, C. (2010). A multi-start evolutionary local search for the two-echelon location routing problem. Hybrid Metaheuristics Lecture Notes in Computer Science, 6373, 88-102.

M AN U

Nguyen, V. P., Prins, C., & Prodhon, C. (2012a). Solving the two-echelon location routing problem by a GRASP reinforced by a learning process and path relinking. European Journal of Operational Research, 216(1), 113-126.

Nguyen, V. P., Prins, C., & Prodhon, C. (2012b). A multi-start iterated local search with tabu list and path relinking for the two-echelon location-routing problem. Engineering Applications of Artificial Intelligence, 25(1), 56-71.

Perboli, G., & Tadei, R. (2010). New families of valid inequalities for the two-echelon vehicle routing problem.

TE D

Electronic notes in discrete mathematics, 36, 639-646.

Perboli, G., Tadei, R., & Fadda, E. (2018). New valid inequalities for the two-echelon capacitated vehicle routing problem. Electronic Notes in Discrete Mathematics, 64(1), 75-84. Pichka, K., Bajgiran, A. H., Petering, M. E., Jang, J., & Yue, X. (2018). The two echelon open location routing problem: Mathematical model and hybrid heuristic. Computers & Industrial Engineering, 121, 97-112.

EP

Santos, F. A., da Cunha, A. S., & Mateus, G. R. (2013). Branch-and-price algorithms for the two-echelon capacitated vehicle routing problem. Optimization Letters, 7(7), 1537-1547. Santos, F. A., Mateus, G. R., & da Cunha, A. S. (2014). A branch-and-cut-and-price algorithm for the two-echelon

AC C

capacitated vehicle routing problem. Transportation Science, 49(2), 355-368. Shen, S. Y., & Honda, M. (2009). Incorporating lateral transfers of vehicles and inventory into an integrated replenishment and routing plan for a three-echelon supply chain. Computers & Industrial Engineering, 56(2),

754-775.

Sitek, P., & Wikarek, J. (2014). A novel integrated approach to the modelling and solving of the Two-Echelon Capacitated Vehicle Routing Problem. Production & Manufacturing Research, 2(1), 326-340.

Solomon, M. M. (1987). Algorithms for the Vehicle Routing and Scheduling Problems with Time Window Constraints. Operations Research, 35(2), 254-265. Wang, Y., Assogba, K., Liu, Y., Ma, X., Xu, M., & Wang, Y. (2018). Two-echelon location-routing optimization with time windows based on customer clustering. Expert Systems with Applications, 104, 244-260. Wang, K., Shao, Y., & Zhou, W. (2017). Matheuristic for a two-echelon capacitated vehicle routing problem with environmental considerations in city logistics service. Transportation Research Part D, 57, 262-276.

22