# A Multiperiod Vehicle Lease Planning for Urban Freight Consolidation Network

Authors: Woosuk Yang; Taesu Cheong; Sang Hwa Song
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This paper considers a multiperiod vehicle lease planning problem for urban freight consolidation centers (UFCCs) in the urban freight transport network where short-term-leased and long-term-leased vehicles are hired together. The objective is to allocate the two kinds of leased vehicles optimally for direct transportation services from the associated origin node to the associated UFCC or from the associated UFCC to the associated destinations so as to satisfy a given set of period-to-period freight demands over a given planning horizon at total minimum vehicle allocation cost subject to demand-dependent transportation time restriction. The problem is formulated as an integer programming model and proven to be NP-hard in a strong sense. Thus, a Lagrangian heuristic is proposed to find a good solution efficiently. Numerical experiments show that the proposed algorithm finds good lower and upper bounds within reasonable time.


---

## Body

## 1. Introduction

Recent worldwide economic growth has accelerated urbanization, and people have been moving to cities for better jobs and more fulfilling lifestyle. In the late 2000s, it was reported that more people lived in cities than in rural areas, and it is expected that more people will move to urban areas in the near future [1, 2]. According to Blanco and Fransoo [3], the number of megacities with at least 10 million people is increasing, and in a decade, their contribution to world economy will steadily grow and constitute more than 20% of world GDP. As people start to live in concentrated and sometimes congested areas, traffic jams and air pollution continue to increase. As a result, many researchers have started to investigate the impacts of congestion due to urban transportation systems [2, 4].Freight transportation is one of the key contributors to traffic congestion and harmful pollutants in cities. Unlike typical passenger cars, freight vehicles are bigger in size and move more frequently due to the nature of their business. Furthermore, a recent trend toward just-in-time delivery makes the problem even worse [5]. To survive in competition and to meet increased customer expectation about timely delivery, companies have dispatched small packages more frequently to customers inside metropolitan areas. As average vehicle utilization went down and the number of vehicles dispatched increased, traffic congestion and the subsequent problems have been exacerbated to such an extent that they could cause serious social issues. To deal with the problems caused by urban freight transportation, integrated approaches for city logistics systems have been proposed. Crainic et al. [5] proposed an integrated short-term scheduling of operations management of logistics resources. Ehmke [1] and Ehmke et al. [6] proposed efficient routing systems integrating traffic information and logistics. Yang et al. [7] investigated a problem of designing a city logistics network considering green-house gas emissions. Thus, the research goal of urban freight transportation and city logistics should be aligned with congestion and pollution reduction with less degradation of the city center commercial activities. In this regard, multitier urban freight transportation system with consolidation and in-out synchronization, so-calledurban freight consolidation, has been proposed as an alternative transportation infrastructure for city logistics.Urban freight consolidation center (UFCC) is a logistics facility located at the boundary of urban areas to serve freight transportation to city centers including retailers, buildings, and construction sites [8]. Products are transported from their origins such as manufacturers to their destinations in cities through urban consolidation centers. They are first moved to urban consolidation centers, where incoming shipments are unloaded, sorted, and consolidated with other products from different origins. Sorted and consolidated products are then transferred directly to outbound vehicles without storage. In just-in-time delivery environment, it is quite important to send orders to customers when necessary, even though shipment size is not enough to fill a full truck load. By consolidating products in an intermediate point between origins and destinations, urban consolidation centers can contribute to an increase of vehicle utilization, an increase of average size of vehicles involved, and the reduction of delivery frequency. In addition to this, the usage of environmentally friendly vehicles such as electric and clean natural gas power vehicles could often contribute to the decrease of overall harmful emissions from freight transportation vehicles. By improving the overall loading factor of a vehicle destined for congested city centers, urban consolidation centers could effectively reduce the total travel distance and further reduce the impact of freight operations on traffic congestion [9]. In fact, the concept of UFCC has been tested with real business practice in European countries and Toronto, Canada [10–12].Research on urban freight consolidation centers focuses on the economic analysis of the consolidated freight transportation system. Su and Roorda [12] and Triantafyllou et al. [13] showed that urban consolidation centers can be successfully operated in real urban environments. It was reported that the trial systems were able to reduce harmful emissions and traffic congestion in city centers when properly managed and synchronized. Marcucci and Danielis [8] showed that, in their analysis, urban consolidation centers could attract a considerable amount of freight shipments bound to urban areas. Zhou and Wang [14] studied the issues related to development and construction of consolidation centers and showed that proper strategies based on public-private partnership can increase overall economic benefits of the system to participants in city logistics. While the feasibility of urban freight consolidation centers has been studied in depth, operation and network planning issues have been discussed (e.g., [15–17]) but still not been fully investigated yet. The detailed analysis of operation and planning of urban freight consolidation system should be done to be deployed in real business practice. In the supply chain context, an urban consolidation center is similar to the transportation system with cross-docking terminals. In cross-docking network, goods are moved through cross-docking terminals where shipments are sorted, consolidated, and transferred to outbound vehicles in a synchronized manner.Research on operation and management of the cross-docking based consolidation strategy has progressed in two directions. One is concerned with problems that are related to the internal operations at the consolidation center. Gue [18] analyzed the effects of scheduling trailers into doors on a layout of a freight consolidation center. Bartholdi III and Gue [19, 20] designed the layout of a freight consolidation center. Li et al. [21] studied a scheduling problem to minimize storage and order picking in a consolidation center. On the other hand, planning problems on a network level have been considered to locate consolidation centers, to allocate vehicles, and to make vehicle consolidation schedules. Ratliff et al. [22] and Chen et al. [23] considered a problem of making vehicle consolidation schedules for a transportation network. Donaldson et al. [24] considered a problem of allocating vehicles and making vehicle consolidation schedules. Sung and Song [25] and Sung and Yang [26] studied an integrated model of locating consolidation centers and allocating vehicles. All the works on a network level about consolidation based transportation have considered static freight demands, not varying with time. However, freight demands in various industries such as food, apparel, electronic goods, and logistics may be dynamic. Especially in urban freight transportation, demands tend to fluctuate over time. As a cost-efficient way of allocating vehicles to satisfy dynamic freight demands, the issue of vehicle supply on lease has received much research attention, in the situation where the unit-period vehicle lease cost depends on the lease term. The unit-period long-term vehicle lease cost is generally cheaper than the unit-period short-term vehicle lease cost [27–29]. In the case of static freight demands, the long-term vehicle lease is obviously better than the short-term vehicle lease, while in the case of dynamic freight demands, either one does not dominate over the other so that the two lease options need to be considered together. Furthermore, in an urban freight consolidation setting, we need to consider consolidation of shipments which makes the lease planning more difficult.Therefore, this paper considers a multiperiod vehicle lease planning problem in an urban freight consolidation network (MVLPUC). The problem, denoted byPMVLPUC, is concerned with optimally allocating the two kinds of leased vehicles for inbound and outbound transportation services so as to satisfy a given set of period-to-period freight demands over a given planning horizon at total minimum vehicle allocation cost subject to demand-dependent transportation time restriction. The planning horizon is divided into discrete time periods such as weeks or months. It is assumed that each freight demand is transported through a single path via one urban consolidation center (where the operations of sorting and consolidating are handled) located between origin and destination nodes, and each outbound vehicle at each urban consolidation center departs as soon as all the associated inbound vehicles arrive and the associated freight demands are sorted appropriately. It is also assumed that an unlimited number of homogeneous capacitated vehicles can be acquired through either long-term lease or short-term lease, with the unit-period short-term vehicle lease cost being greater than or equal to the unit-period long-term vehicle lease cost.This paper is organized as follows: Section2 introduces the problem formulation for PMVLPUC. In Section 3, we discuss the procedure for solving PMVLPUC based on Lagrangian relaxation and problem decompositions. We present the numerical experiments in Section 4 and then conclude our discussion in Section 5.


## 2. Model

In this section, we present the mathematical formulation for ProblemPMVLPUC. Before introducing the model formulation, all the parameters and decision variables used in this paper are given as follows.Sets and Parameters
I
,
J
,
K
,
V: sets of origin nodes representing manufacturers, destination nodes representing retailers, intermediate nodes representing UFCCs, and all nodes (i.e., V=I∪K∪J), respectively.
E
IN
,
E
OUT: sets of edges representing potential direct services from origin nodes to intermediate nodes and from intermediate nodes to destination nodes, respectively (EIN={e=(i,k)∣i∈I,k∈K} and EOUT={e=(k,j)∣k∈K,j∈J}).
E: set of all edges (i.e., E=EIN∪EOUT).
Q: set of freight demands, defined by an ordered pair of two nodes (i,j) for all i∈I and j∈J.
Q
k: subset of Q, each demand of which can be transported through the intermediated node k∈K within the associated transportation time restriction (to be explained later).
T: set of time periods.
d
i
j
p: quantity of freight demand (i,j)∈Q at time period p∈T.
h
i
k
L
,
h
k
j
L: unit-period long-term vehicle lease cost for (i,k)∈EIN and (k,j)∈EOUT, respectively.
h
i
k
p
S
,
h
k
j
p
S: unit-period short-term vehicle lease cost for (i,k)∈EIN and (k,j)∈EOUT, respectively, at time period p∈T (hikL≤hikpS and hkjL≤hkjpS).
Γ: vehicle capacity.
t
i
k
,
t
k
j: transportation time elapsed for each (i,k)∈EIN and (k,j)∈EOUT, respectively.
o
k: handling (sorting and consolidating) time at the intermediate node k∈K.
T
L
i
j: transportation time restriction required for freight demand (i,j)∈Q.
Decision Variables
X
i
j
k
p: 1 if freight demand (i,j)∈Q is transported in period p∈T through the intermediate node k∈K and 0 otherwise.
Y
i
k
,
Y
k
j: numbers of long-term-leased vehicles allocated for (i,k)∈EIN and (k,j)∈EOUT, respectively.
Z
i
k
p
,
Z
k
j
p: numbers of short-term-leased vehicles allocated for (i,k)∈EIN and (k,j)∈EOUT, respectively, in period p∈T.We note that, for notational simplicity, notations (or subscripts) e and (i,k) (or (k,j)) are used alternately to refer to freight demand (i,k) (or (k,j)) on the corresponding edge e. Furthermore, the two terms “edge” and “direct service” are used interchangeably in this paper. Then, we now present the problem formulation for Problem PMVLPUC as follows.Problem
P
M
V
L
P
U
C
(1a)
min
⁡
∑
i
,
k
∈
E
IN
T
h
i
k
L
Y
i
k
+
∑
p
∈
T
h
i
k
p
S
Z
i
k
p
+
∑
k
,
j
∈
E
OUT
T
h
k
j
L
Y
k
j
+
∑
p
∈
T
h
k
j
p
S
Z
k
j
p
(1b)
subject  to
∑
k
∈
K
∣
i
,
j
∈
Q
k
X
i
j
k
p
=
1
,
∀
i
,
j
∈
Q
,
∀
p
∈
T
(1c)
∑
j
∈
J
∣
i
,
j
∈
Q
k
d
i
j
p
X
i
j
k
p
≤
Γ
Y
i
k
+
Z
i
k
p
,
∀
i
,
k
∈
E
IN
,
∀
p
∈
T
(1d)
∑
i
∈
I
∣
i
,
j
∈
Q
k
d
i
j
p
X
i
j
k
p
≤
Γ
Y
k
j
+
Z
k
j
p
,
∀
k
,
j
∈
E
OUT
,
∀
p
∈
T
(1e)
Y
i
k
,
Z
i
k
p
,
Y
k
j
,
Z
k
j
p
∈
Z
+
,
X
i
j
k
p
∈
0,1
.The objective function (1a) represents the cost of allocating any long-term-leased and short-term-leased vehicles for edges. Constraints (1b) imply that all the freight demands have to be serviced for each period. Constraints (1c) and (1d) require that the total amount of demands transported through any edge should not exceed the total capacity of any allocated vehicles for each period. In regard to the set Qk in constraints (1b), (1c), and (1d), we assume that there exists the transportation time limit TLij for each freight demand (i,j) such that the sum of the transportation times between nodes and the handling time at an intermediate node k does not exceed TLij (i.e., tik+ok+tkj≤TLij). Thus, for each k∈K, the elements of set Qk can be identified as freight demands (i,j) in Q that satisfy the aforementioned condition with TLij at a preprocessing stage.


## 3. Solution Approach

In this section, we propose a heuristic approach based on Lagrangian relaxation. ProblemPMVLPUC can be proven as NP-hard in a strong sense in the same manner as shown in Sung and Song [25]. Moreover, PMVLPUC considers time-varying demands, so it may be too complex to derive an exact algorithm to find optimal solutions even for small-sized problem instances. Therefore, we here propose an efficient heuristic method based on the Lagrangian relaxation.


### 3.1. Lagrangian Relaxation

Before relaxing constraints inPMVLPUC for Lagrangian relaxation heuristic, the following problem PMVLPUC+ is derived by introducing another decision variable X¯ijkp (equivalent to variables Xijkp) and adding a set of constraints Xijkp≤X¯ijkp to PMVLPUC.Problem
P
M
V
L
P
U
C
+
(2a)
υ
OPT
=
min
⁡
∑
i
,
k
∈
E
IN
T
h
i
k
L
Y
i
k
+
∑
p
∈
T
h
i
k
p
S
Z
i
k
p
+
∑
k
,
j
∈
E
OUT
T
h
k
j
L
Y
k
j
+
∑
p
∈
T
h
k
j
p
S
Z
k
j
p
(2b)
subject  to
∑
k
∈
K
∣
i
,
j
∈
Q
k
X
i
j
k
p
=
1
,
∀
i
,
j
∈
Q
,
∀
p
∈
T
(2c)
∑
j
∈
J
∣
i
,
j
∈
Q
k
d
i
j
p
X
i
j
k
p
≤
Γ
Y
i
k
+
Z
i
k
p
,
∀
i
,
k
∈
E
IN
,
∀
p
∈
T
(2d)
X
i
j
k
p
≤
X
¯
i
j
k
p
,
∀
k
∈
K
,
∀
i
,
j
∈
Q
k
,
∀
p
∈
T
(2e)
∑
i
∈
I
∣
i
,
j
∈
Q
k
d
i
j
p
X
¯
i
j
k
p
≤
Γ
Y
k
j
+
Z
k
j
p
,
∀
k
,
j
∈
E
OUT
,
∀
p
∈
T
(2f)
Y
i
k
,
Z
i
k
p
,
Y
k
j
,
Z
k
j
p
∈
Z
+
,
X
i
j
k
p
,
X
¯
i
j
k
p
∈
0,1
.Problem PMVLPUC+ has some obviously redundant variables X¯ijkp and the associated constraints (2d) while it yields the interesting problem structure which is good to apply the Lagrangian relaxation method. That is, if constraints (2b) and (2d) in PMVLPUC+ are Lagrangian relaxed, then the resulting problem can be decomposed into single-edge problems. Thus, constraints (2b) and (2d) are Lagrangian relaxed with Lagrange multipliers μijp′ (unrestricted) and μijkp′′ (≥0), respectively. For given μ≡(μijp′,μijkp′′)∀i,j,k,p, the resulting problem PLRμ can be derived as follows.Problem
P
L
R
μ
(3)
υ
μ
=
min
⁡
∑
i
,
j
∈
Q
∑
p
∈
T
μ
i
j
p
′
+
∑
i
,
k
∈
E
IN
T
h
i
k
L
Y
i
k
+
∑
p
∈
T
h
i
k
p
S
Z
i
k
p
-
∑
j
∈
J
∣
i
,
j
∈
Q
k
θ
i
j
k
p
X
i
j
k
p
+
∑
k
,
j
∈
E
OUT
T
h
k
j
L
Y
k
j
+
∑
p
∈
T
h
k
j
p
S
Z
k
j
p
-
∑
i
∈
I
∣
i
,
j
∈
Q
k
δ
i
j
k
p
X
¯
i
j
k
p
subject  to
2c
,
2e
,
2f
,
where
θ
i
j
k
p
=
μ
i
j
p
′
-
μ
i
j
k
p
′
′
,
δ
i
j
k
p
=
μ
i
j
k
p
′
′
.It is evident that υ(μ) is a lower bound on υOPT for any given μ. Let υ∗(μ)=maxμυ(μ), and let υLP be a lower bound (on υOPT) obtained by solving the linear programming (LP) relaxation problem of PMVLPUC+. Then, the proposed Lagrangian relaxation PLRμ provides a good lower bound as indicated in Proposition 1. Since the proposition below is obvious to show, we omit the proof.Proposition 1.
ConsiderυOPT≥υ∗(μ)≥υLP.Furthermore, as mentioned earlier, ProblemPLRμ can be decomposed into single-edge problems of selecting demand and allocating vehicles (SEPDV), resulting in PSEPDVik,μ and PSEPDVkj,μ for each (i,k)∈EIN and (k,j)∈EOUT, respectively, as follows.Problem
P
S
E
P
D
V
i
k
,
μ
(4a)
υ
SEPDV
i
k
μ
=
min
⁡
T
h
i
k
L
Y
i
k
+
∑
p
∈
T
h
i
k
p
S
Z
i
k
p
-
∑
j
∈
J
∣
i
,
j
∈
Q
k
θ
i
j
k
p
X
i
j
k
p
(4b)
subject  to
∑
j
∈
J
∣
i
,
j
∈
Q
k
d
i
j
p
X
i
j
k
p
≤
Γ
Y
i
k
+
Z
i
k
p
,
∀
p
∈
T
(4c)
Y
i
k
,
Z
i
k
p
∈
Z
+
,
X
i
j
k
p
∈
0,1
,
∀
p
∈
T
.
Problem
P
S
E
P
D
V
k
j
,
μ
(5)
υ
SEPDV
k
j
μ
=
m
i
n
T
h
k
j
L
Y
k
j
+
∑
p
∈
T
h
k
j
p
S
Z
k
j
p
-
∑
i
∈
I
∣
i
,
j
∈
Q
k
δ
i
j
k
p
X
¯
i
j
k
p
subject  to
∑
i
∈
I
∣
i
,
j
∈
Q
k
d
i
j
p
X
¯
i
j
k
p
≤
Γ
Y
k
j
+
Z
k
j
p
,
∀
p
∈
T
Y
k
j
,
Z
k
j
p
∈
Z
+
,
X
¯
i
j
k
p
∈
0,1
,
∀
p
∈
T
.It is straightforward to show that, for given μ,(6)υμ=∑i,j∈Q∑p∈Tμijp′+∑i,k∈EINυSEPDVikμ+∑k,j∈EOUTυSEPDVkjμand hence, υ(μ) can be evaluated through υSEPDVik(μ) and υSEPDVkj(μ). Since PSEPDVik,μ and PSEPDVkj,μ have the same problem structure, we only discuss the procedure for solving PSEPDVik,μ in the next section.


### 3.2. SolvingPSEPDVik,μ

In this section, we present the solution procedure forPSEPDVik,μ so that we can eventually evaluate the value υ(μ) of PLRμ for each μ. First, for given μ, Proposition 2 and Corollary 3 are characterized so as to reduce the solution space of X and (Y,Z) variables, respectively, in PSEPDVik,μ.Proposition 2.
For givenμ, the following statements hold in PSEPDVik,μ: (a)
Ifθijkp≤0 for some j∈J and p∈T, then there is an optimal solution with Xijkp=0.(b)
Ifθijkp≥hikpS⌈dijp/Γ⌉ for some j∈J and p∈T, then there is an optimal solution with Xijkp=1.Proof.
(a) Let us compare the situations ofXijkp=1 and Xijkp=0 for some j∈J and p∈T such that θijkp≤0. The value of the objective function (4a) for the former case cannot be smaller than that for the latter case because the coefficient of the X variable in (4a) is nonnegative. Also, the former case makes the constraint (4b) tighter than the latter case. Therefore, the former case (i.e., Xijkp=1) cannot give a better solution than the latter case (i.e., Xijkp=0).
(b) This can be shown in the same manner as in the proof of (a), and hence we omit its proof.LetQik be a set of freight demands that can be transported through (i,k) where i∈I and k∈K. Note that Qk=⋃i∈I‍Qik. For p∈T, let Q~ikp0 and Q~ikp1 be subsets of Qik satisfying Propositions 2(a) and 2(b) in PSEPDVik,μ, respectively, and let Q~ikp=Q~ikp0∪Q~ikp1.Corollary 3.
For givenμ and {Xijkp} fixed at the associated values according to Proposition 2, there is an optimal solution in PSEPDVik,μ with (Yik,Zikp)∈(y,zp)∣Y~ikmin⁡≤y≤Y~ikmax⁡,0≤zp≤Z~ikp,ymax⁡,∀p∈T, where (7)Y~ikmin⁡=minp∈T⁡∑j∈J∣i,j∈Q~ikp1dijpΓ,Y~ikmax⁡=maxp∈T⁡∑j∈J∣i,j∈Qik∖Q~ikp0dijpΓ,Z~ikp,ymax⁡=∑j∈J∣i,j∈Qik∖Q~ikp0dijpΓ-y.Proof.
The constraint (4b) can be rewritten as(8)∑j∈J∣i,j∈Q~ikp1dijp+∑j∈J∣i,j∈Qik∖Q~ikpdijpXijkp≤ΓYik+Zikpfor all p∈T. Since all the variables are binary variables, Y~ikmin⁡≤∑j∈J∣i,j∈Q~ikp1‍dijp/Γ≤Yik+Zikp holds for all p∈T. Since it is assumed that 0<hikL≤hikpS for all p∈T, any feasible solution with (Xijkp,Yik,Zikp)=(xjp,y,zp), where y<Y~ikmin⁡ cannot yield a better solution value than the feasible solution with (Xijkp,Yik,Zikp)=xjp,Y~ikmin⁡,max0,zp-Y~ikmin⁡+y. Moveover, any feasible solution with (Xijkp,Yik,Zikp)=(xjp,y,zp), where y>Y~ikmax⁡ cannot yield a better solution value than the feasible solution with (Xijkp,Yik,Zikp)=xjp,Y~ikmax⁡,0. Therefore, there is an optimal solution with Y~ikmin⁡≤Yik≤Y~ikmax⁡. Also, with Yik fixed at y∈y′∈Z+∣Y~ikmin⁡≤y′≤Y~ikmax⁡, it is obvious that there is an optimal solution with 0≤Zikp≤Z~ikp,ymax⁡.Then, the restricted problem withYik fixed at y∈y′∈Z+∣Y~ikmin⁡≤y′≤Y~ikmax⁡ in PSEPDVik,μ, denoted by PSEPDVik,μ,y, can be decomposed into single-period bounded knapsack problems (BKP), PBKPikp,μ,y, as follows.Problem
P
B
K
P
i
k
p
,
μ
,
y
(9)
υ
BKP
i
k
p
μ
,
y
=
max
⁡
h
i
k
p
S
Z
¯
i
k
p
+
∑
j
∈
J
∣
i
,
j
∈
Q
i
k
∖
Q
~
i
k
p
0
θ
i
j
k
p
X
i
j
k
p
subject  to
∑
j
∈
J
∣
i
,
j
∈
Q
i
k
∖
Q
~
i
k
p
d
i
j
p
X
i
j
k
p
+
Γ
Z
¯
i
k
p
≤
Γ
∑
j
∈
J
∣
i
,
j
∈
Q
i
k
∖
Q
~
i
k
p
0
d
i
j
p
Γ
-
∑
j
∈
J
∣
i
,
j
∈
Q
~
i
k
p
1
d
i
j
p
Z
¯
i
k
p
∈
0,1
,
…
,
∑
j
∈
J
∣
i
,
j
∈
Q
i
k
∖
Q
~
i
k
p
0
d
i
j
p
Γ
-
y
,
X
i
j
k
p
∈
0,1
.Note that Z¯ikp=Z~ikp,ymax⁡-Zikp. We next discuss how to solve PBKPikp,μ,y. Let υSEPDVik(μ,y) denote the optimal value of PSEPDVik,μ,y. We first transform PBKPikp,μ,y into the associated 0-1 knapsack problem, and the knapsack problem can be solved by the dynamic programming algorithm proposed by Toth [30]. We remark that Z¯ikp takes the finite number of integer values in the set Ψ≡0,1,…,∑j∈J∣i,j∈Qik∖Q~ikp0‍dijp/Γ-y. Furthermore, for each value of Z¯ikp in Ψ, it is straightforward to show that the problem PBKPikp,μ,y reduces to a 0-1 knapsack problem with decision variables Xijkp for all j. Therefore, we iteratively solve a 0-1 knapsack problem for each value of Z¯ikp in Ψ, compare the objective function values for each case, and then identify the optimal solution for PBKPikp,μ,y. Once all the problems of PBKPikp,μ,y are solved for each value of Yik as mentioned above, the optimal value of the problem PSEPDVik,μ, υSEPDVik(μ), is computed as(10)υSEPDVikμ=minY~ikmin⁡≤y≤Y~ikmax⁡⁡υSEPDVikμ,y=minY~ikmin⁡≤y≤Y~ikmax⁡⁡ThikLy+∑p∈ThikpSZ~ikp,ymax⁡-υBKPikpμ,y-∑j∈J∣i,j∈Q~ikp1θijkp.The number of problems PBKPikp,μ,y and PSEPDVik,μ,y to be solved can be further reduced by Propositions 4(a) and 4(b), respectively. Let Xijkp,y∗,Z¯ikp,y∗ be the optimal solution of PBKPikp,μ,y with explicit dependence on y, and furthermore, let us denote Try as Try=p∈T∣Z¯ikp,y∗≤Z~ikp,ymax⁡-r=Z~ikp,(y+r)max⁡ for Y~ikmin⁡≤y≤Y~ikmax⁡ and 0<r≤Z~ikp,ymax⁡. Note that Try⊂Tly⊂T for 0≤l<r≤Z~ikp,ymax⁡.Proposition 4.
(a) GivenXijkp,y∗,Z¯ikp,y∗ for some y∈y′∈Z+∣Y~ikmin⁡≤y′≤Y~ikmax⁡ and some p∈T, if p∈Try for some r∈r′∈Z+∣0<r′≤Z~ikp,ymax⁡, then Xijkp,(y+r)∗,Z¯ikp,(y+r)∗=Xijkp,y∗,Z¯ikp,y∗.
(b) GivenXijkp,y∗,Z¯ikp,y∗ for some y∈y′∈Z+∣Y~ikmin⁡≤y′≤Y~ikmax⁡ and all p∈T, if ∑p∈Try‍hikpS≥ThikL for some r∈r′∈Z+∣0<r′≤Z~ikp,ymax⁡, then υSEPDVik(μ,y+l)≥υSEPDVik(μ,y+r) for 0≤l<r.In summary, we now present the solution procedure forPSEPDVik,μ,PROC(SEPDVik,μ), in Algorithm 1 based on the whole discussion above. For reference, Figure 1 illustrates the relationship among all the problems discussed in this section.Algorithm 1:Solution procedure for PSEPDVik,μ.
P
ROC(SEPDVik,μ)
begin
(1) LetYik∗,Zikp∗ for all p∈T and Xijkp∗ for all j∈{j∈J∣(i,j)∈Qk} be Y, Z and X-variables of
the optimal solution ofPSEPDVik,μ respectively;
(2) SetυSEPDVik(μ)=∞;
(3) Fix someXijkps at the associated values according to Proposition 2 and determine Y~ikmin,
Y~ikmax and Z~ikp,ymax for all y∈{y′∈Z+∣Y~ikmin≤y′≤Y~ikmax} by Corollary 3;
(4) SetTry=ϕ for all r∈{r′∈Z+∣0<r′≤Z~ikp,ymax⁡} and y∈{y′∈Z+∣Y~ikmin⁡≤y′≤Y~ikmax⁡};
(5) Setαy=0 and βpy=0 for all y∈{y′∈Z+∣Y~ikmin⁡≤y′≤Y~ikmax⁡} and p∈T where αy and
βpy are binary variables having 1 only if PSEPDVik,μ and PBKPikp,μ,y are not necessary
to be solved respectively, and 0 otherwise;
(6) Sety=Y~ikmin⁡;
(7)whiley≤Y~ikmaxdo
(8)ifαy=0then
(9)     Setp=1;
(10)whilep≤Tdo
(11)ifβpy=0then
(12)       SolvePBKPikp,μ,y and set the resulting solution as (Xijkp,y∗,Z¯ikp,y∗);
(13)       SetTr-l,y+l=Tr-l,y+1∪{p} for all r∈{r′∈Z+∖{0}∣Z¯ikp,y∗≤Z~ikp,ymax-r′} and l∈{l′∈Z+∣0≤l′<r};
end
(14)     Set(Xijkp,y+r∗,Z¯ikp,y+r∗)=(Xijkp,y∗,Z¯ikp,y∗) and βp,y+r=1 for all
r∈{r′∈Z+∖{0}∣Z¯ikp,y∗≤Z~ikp,ymax⁡-r′} by Proposition 4(a);
(15)     Setp=p+1;
end
(16)     Setαy+l=1 for all l∈l′∈Z+∣0<l′<r∗ by Proposition 4(b) where
r∗=maxr∈Rr such that R={r∈Z+∣0<r≤Z~ikp,ymax,∑p∈TryhikpS≥ThikL};
(17)ifυSEPDVik(μ,y)<υSEPDVik(μ)then
(18)     SetυSEPDVik(μ)=υSEPDVik(μ,y) (refer to (10));
(19)     SetYik∗=y and (Xijkp∗,Z¯ikp∗)=(Xijkp,y∗,Z~ikp,ymax-Z¯ikp,y∗) for all p∈T;
end
end
(20)   Sety=y+1;
end
endFigure 1
Problem decomposition ofPMVLPUC.


### 3.3. Finding the Lagrange Multipliers

So far, we examined how to address the problemPSEPDVik,μ for each (i,k)∈EIN when the Lagrangian multiplier μ is given. We now discuss the procedure for finding the Lagrangian multiplier μ. For the purpose, we adapt the subgradient optimization procedure which is one of the most popular methods to find a good set of Lagrange multipliers [31]. Specifically, the Lagrange multipliers μ≡(μijp′,μijkp′′)∀i,j,k,p are generated iteratively as(i)
μijp,s+1′=μijp,s′+γsGijp,s′ for all (i,j)∈Qk, p∈T,(ii)
μijkp,s+1′′=max0,μijkp,s′′+γs+γsGijp,s′ for all k∈T, (i,j)∈Qk, p∈T,where Gijp,s′ and Gijkp,s′′ are the subgradients for constraints (2b) and (2d), respectively, and γs is the stepsize at iteration s.When it comes to the subgradientsGijp,s′ and Gijkp,s′′ in the formulas above, they are determined by the solution of PLRμ given μ as follows: Gijp,s′=1-∑{k∈K∣(i,j)∈Qk}‍Xijkp for all (i,j)∈Qk, p∈T, and Gijkp,s′′=Xijkp-X¯ijkp for all k∈K,(i,j)∈Qk, p∈T. Moreover, the stepsize γs at iteration s is determined by(11)γs=ηUB-υμs∑i,j∈Q∑p∈TGijp,s′2+∑k∈K∣i,j∈QkGijkp,s′′2,where UB represents the best upper bound of υOPT found up to iteration s and η is a control parameter for γs. In this paper, η is set to 2 at the beginning and then halved if the lower bound is not improved in a predetermined number of consecutive iterations.


### 3.4. Finding a Feasible Solution

A solution obtained by solvingPLRμ based on the discussion in Section 3.2 may be infeasible to PMVLPUC+, and hence we propose a two-phase heuristic to deal with infeasibility issues as follows: In the first phase, a construction heuristic (CH) is derived to modify any infeasible solutions to PMVLPUC+ into feasible ones, and then, in the second phase, a tabu search-based heuristic (TSH) is utilized to improve the solution of (CH). The details are presented in the subsequent sections.


#### 3.4.1. Construction Heuristic (CH)

As mentioned above, a construction heuristic (CH) intends to deal with possible infeasibility issues of solutions obtained by solving PLRμ, and (CH) is performed at each subgradient iteration. Let Xijkp∗ and X¯ijkp∗ be the X- and X¯-variable values of the optimal solution of PSEPDVik,μ and PSEPDVkj,μ, respectively. First, the corresponding demand for all (i,j)∈Q and p∈T is reassigned to the UFCC with the lexicographical maximum value of the three-tuple relation Xijkp∗+X¯ijkp∗,θijkpXijkp∗+δijkpX¯ijkp∗,-hikL+hkjL. Then, from the reassigned demands, demand quantity Dep for all e∈E and p∈T is computed as Dep=∑{j∈J∣(i,j)∈Qk}‍dijpXijkp or Dep=∑{i∈I∣(i,j)∈Qk}‍dijpX¯ijkp. Finally, to allocate vehicles for each e∈E at the minimum vehicle cost, we solve the following vehicle allocation problem (VAP), PVAPe.Problem
P
V
A
P
e
(12a)
υ
VAP
e
=
min
⁡
T
h
e
L
Y
e
+
∑
p
∈
T
h
e
p
S
Z
e
p
(12b)
subject  to
D
e
p
≤
Γ
Y
e
+
Z
e
p
,
∀
p
∈
T
(12c)
Y
e
,
Z
e
p
∈
Z
+
,
∀
p
∈
T
.Proposition 5 characterizes the optimal solutions of PVAPe.Proposition 5.
There is an optimal solution ofPVAPe with (Ye,Zep)∈(y,Zep,y)∣y=y[1],…,y[|T|], where Zep,y=max0,⌈Dep/Γ⌉-y and y[i] is the ith largest y.Proof.
As in Corollary3, it is straightforward to show that there is an optimal solution of PVAPe with (Ye,Zep)∈(y,Zep,y)∣y[1]≤y≤y[|T|]. For an integer m such that y[n]≤m<y[n+1] and for n∈1,…,|T|-1,(13)Zep,m=0ifp≤nyp-motherwise,Zep,m+1=0ifp≤nyp-m+1otherwise.Let υVAPe(y) be the value of the objective function (12a) associated with Ye fixed at y∈[y1,y[|T|]] in PVAPe. Then, υVAPe(m+1)-υVAPe(m)=(m+1)|T|heL-∑p=1|T|‍Ze[p],m+1 − m|T|heL-∑p=1|T|‍Ze[p],m=|T|heL-∑p=n+1|T|‍he[p]S. It implies that υVAPe(m) is a linear function in m over the range y[n]≤m<y[n+1]. Therefore, there is an optimal solution of PVAPe with (Ye,Zep)∈(y,Zep,y)∣y=y[1],…,y[|T|].By Proposition5, PVAPe can be solved by finding a solution of having the minimum cost among |T| number of solutions which are associated with Ye fixed at y[p] for all p∈T.


#### 3.4.2. Tabu Search-Based Heuristic (TSH)

After the first phase (CH), a tabu search-based heuristic (TSH) is performed to improve the solution obtained from (CH). We note that (TSH) is expected to require much more computational time than (CH), and hence, to save the elapsed time of the overall algorithm, (TSH) is executed only at the subgradient iterations satisfying the following two conditions ((a) and (b)) together: (a) the best lower bound is improved and (b) the solution obtained by (CH) at the current subgradient iteration is different from the solutions obtained by (CH) so far. In this paper, we call the conditions (a) and (b) as “(TSH)-activating conditions.”Tabu search is one of the popular metaheuristics designed to escape the trap of local optimality [32], and it has been used frequently in various network-related problems as in Sung and Song [25], Crevier et al. [33], and Üster and Maheshwari [34]. We now explain the main features of the tabu search used in (TSH), includingmove operation,tabu restriction andaspiration criterion,move evaluation, andtermination condition.Themove operation used in this paper is to modify assignment of some freight demand, say q=(i,j)∈Q, from its current UFCC, say l∈K, to the alternative UFCC, say k∈K, for some period, say p∈T. It corresponds to modifying the value of Xijlp from 1 to 0 and the value of Xijkp from 0 to 1 in Problem PMVLPUC+. Thetabu restriction is employed such that the UFCC assignment for any freight demand should not be changed duringtabu_size number of iterations except for the case that satisfies the aspiration criterion (to be explained later). To implement the tabu restriction, as a new solution is selected, the information of the freight demand associated with the move operation is recorded in a list (calledtabu list) and kept in the list over thetabu_size number of iterations. In this paper, thetabu_size is chosen randomly from the discrete uniform distribution in the closed interval [10,20]. Theaspiration criterion adapted here is as follows: a solution from the move operation for the freight demands recorded in thetabu list is eligible to be chosen as a new solution only if the move improves the best upper bound. Moreover, when we choose a new solution among neighbor solutions at each iteration, themove evaluation function used in this paper is to maximize {mv-A·n} where mv, n, and A represent the decrease amount of the objective function value, the frequency that the corresponding move is used to produce a new solution so far, and the associated penalty with the frequency, respectively. When it comes to mv, it is computed simply by subtracting the objective value of the new solution from the objective value of the current solution, and if mv is positive, A is set to 0 and 50 otherwise. We note that this move evaluation function is known to have an advantage in diversifying the search direction efficiently [32]. Lastly, if the total number of performed iterations amounts tomax_iter or the best upper bound has not been improved duringmax_no_imp consecutive iterations, (TSH) isterminated. In this paper,max_iter andmax_no_imp are set at 500 and 30, respectively.


### 3.5. Overall Solution Algorithm

We now summarize all the discussion and present the algorithmPROC(MVLPUC) for solving PMVLPUC in Algorithm 2. It can be terminated either after a predetermined number (1000 in this paper) of iterations or when η gets too small (i.e., η≤0.005 in this paper). It can also be terminated when UB-LB/LB×100 is less than 0.1%, where UB and LB represent the best upper and lower bounds searched so far, respectively.Algorithm 2:Overall solution procedure for PMVLPUC.
P
ROC(MVLPUC)
begin
(1) Initialize the necessary parameters and Lagrange multipliers (i.e.,s=1, γs=2, η=2,
LB=-∞, UB=∞, μijp,s′=1 and μijkp,s′′=1);
(2)   InitializeNILB=0 where NILB represents the number of consecutive iterations for
which the best lower bound has not been improved;
repeat
(3)     Solve the Lagrangian problemPL(μs), and let υ(μs) be the resulting objective value;
ifLB<υ(μs)then
(4)      SetLB=υ(μs) and NILB=0;
else
(5)     SetNILB=NILB+1
end
(6)    Find a feasible solution(X,X¯,Y,Z)UB,s by (CH) where (X,X¯,Y,Z)UB,s represents the
feasible solution toPMVLPUC+ at iteration s;
if (TSH)-activating conditions are satisfiedthen
(7)     Update(X,X¯,Y,Z)UB,s by (TSH);
end
(8)      LetυUB,s be the objective value of the solution (X,X¯,Y,Z)UB,s for PMVLPUC+;
(9)ifUB>υUB,sthen
(10)    SetUB=υUB,s, and set (X,X¯,Y,Z)∗=(X,X¯,Y,Z)UB,s where (X,X¯,Y,Z)∗ represents
the best feasible solution searched until iterations;
end
ifNILB=30then
(11)    Setη=η/2 and NILB=0;
end
(12)    Updateμ by the subgradient optimization procedure (refer to Section 3.3);
(13)    Sets=s+1;
untilthe aforementioned termination conditions (i.e., s>1000, η≤0.005 or (UB-LB/LB)×100<0.1) are met;
end


## 3.1. Lagrangian Relaxation

Before relaxing constraints inPMVLPUC for Lagrangian relaxation heuristic, the following problem PMVLPUC+ is derived by introducing another decision variable X¯ijkp (equivalent to variables Xijkp) and adding a set of constraints Xijkp≤X¯ijkp to PMVLPUC.Problem
P
M
V
L
P
U
C
+
(2a)
υ
OPT
=
min
⁡
∑
i
,
k
∈
E
IN
T
h
i
k
L
Y
i
k
+
∑
p
∈
T
h
i
k
p
S
Z
i
k
p
+
∑
k
,
j
∈
E
OUT
T
h
k
j
L
Y
k
j
+
∑
p
∈
T
h
k
j
p
S
Z
k
j
p
(2b)
subject  to
∑
k
∈
K
∣
i
,
j
∈
Q
k
X
i
j
k
p
=
1
,
∀
i
,
j
∈
Q
,
∀
p
∈
T
(2c)
∑
j
∈
J
∣
i
,
j
∈
Q
k
d
i
j
p
X
i
j
k
p
≤
Γ
Y
i
k
+
Z
i
k
p
,
∀
i
,
k
∈
E
IN
,
∀
p
∈
T
(2d)
X
i
j
k
p
≤
X
¯
i
j
k
p
,
∀
k
∈
K
,
∀
i
,
j
∈
Q
k
,
∀
p
∈
T
(2e)
∑
i
∈
I
∣
i
,
j
∈
Q
k
d
i
j
p
X
¯
i
j
k
p
≤
Γ
Y
k
j
+
Z
k
j
p
,
∀
k
,
j
∈
E
OUT
,
∀
p
∈
T
(2f)
Y
i
k
,
Z
i
k
p
,
Y
k
j
,
Z
k
j
p
∈
Z
+
,
X
i
j
k
p
,
X
¯
i
j
k
p
∈
0,1
.Problem PMVLPUC+ has some obviously redundant variables X¯ijkp and the associated constraints (2d) while it yields the interesting problem structure which is good to apply the Lagrangian relaxation method. That is, if constraints (2b) and (2d) in PMVLPUC+ are Lagrangian relaxed, then the resulting problem can be decomposed into single-edge problems. Thus, constraints (2b) and (2d) are Lagrangian relaxed with Lagrange multipliers μijp′ (unrestricted) and μijkp′′ (≥0), respectively. For given μ≡(μijp′,μijkp′′)∀i,j,k,p, the resulting problem PLRμ can be derived as follows.Problem
P
L
R
μ
(3)
υ
μ
=
min
⁡
∑
i
,
j
∈
Q
∑
p
∈
T
μ
i
j
p
′
+
∑
i
,
k
∈
E
IN
T
h
i
k
L
Y
i
k
+
∑
p
∈
T
h
i
k
p
S
Z
i
k
p
-
∑
j
∈
J
∣
i
,
j
∈
Q
k
θ
i
j
k
p
X
i
j
k
p
+
∑
k
,
j
∈
E
OUT
T
h
k
j
L
Y
k
j
+
∑
p
∈
T
h
k
j
p
S
Z
k
j
p
-
∑
i
∈
I
∣
i
,
j
∈
Q
k
δ
i
j
k
p
X
¯
i
j
k
p
subject  to
2c
,
2e
,
2f
,
where
θ
i
j
k
p
=
μ
i
j
p
′
-
μ
i
j
k
p
′
′
,
δ
i
j
k
p
=
μ
i
j
k
p
′
′
.It is evident that υ(μ) is a lower bound on υOPT for any given μ. Let υ∗(μ)=maxμυ(μ), and let υLP be a lower bound (on υOPT) obtained by solving the linear programming (LP) relaxation problem of PMVLPUC+. Then, the proposed Lagrangian relaxation PLRμ provides a good lower bound as indicated in Proposition 1. Since the proposition below is obvious to show, we omit the proof.Proposition 1.
ConsiderυOPT≥υ∗(μ)≥υLP.Furthermore, as mentioned earlier, ProblemPLRμ can be decomposed into single-edge problems of selecting demand and allocating vehicles (SEPDV), resulting in PSEPDVik,μ and PSEPDVkj,μ for each (i,k)∈EIN and (k,j)∈EOUT, respectively, as follows.Problem
P
S
E
P
D
V
i
k
,
μ
(4a)
υ
SEPDV
i
k
μ
=
min
⁡
T
h
i
k
L
Y
i
k
+
∑
p
∈
T
h
i
k
p
S
Z
i
k
p
-
∑
j
∈
J
∣
i
,
j
∈
Q
k
θ
i
j
k
p
X
i
j
k
p
(4b)
subject  to
∑
j
∈
J
∣
i
,
j
∈
Q
k
d
i
j
p
X
i
j
k
p
≤
Γ
Y
i
k
+
Z
i
k
p
,
∀
p
∈
T
(4c)
Y
i
k
,
Z
i
k
p
∈
Z
+
,
X
i
j
k
p
∈
0,1
,
∀
p
∈
T
.
Problem
P
S
E
P
D
V
k
j
,
μ
(5)
υ
SEPDV
k
j
μ
=
m
i
n
T
h
k
j
L
Y
k
j
+
∑
p
∈
T
h
k
j
p
S
Z
k
j
p
-
∑
i
∈
I
∣
i
,
j
∈
Q
k
δ
i
j
k
p
X
¯
i
j
k
p
subject  to
∑
i
∈
I
∣
i
,
j
∈
Q
k
d
i
j
p
X
¯
i
j
k
p
≤
Γ
Y
k
j
+
Z
k
j
p
,
∀
p
∈
T
Y
k
j
,
Z
k
j
p
∈
Z
+
,
X
¯
i
j
k
p
∈
0,1
,
∀
p
∈
T
.It is straightforward to show that, for given μ,(6)υμ=∑i,j∈Q∑p∈Tμijp′+∑i,k∈EINυSEPDVikμ+∑k,j∈EOUTυSEPDVkjμand hence, υ(μ) can be evaluated through υSEPDVik(μ) and υSEPDVkj(μ). Since PSEPDVik,μ and PSEPDVkj,μ have the same problem structure, we only discuss the procedure for solving PSEPDVik,μ in the next section.


## 3.2. SolvingPSEPDVik,μ

In this section, we present the solution procedure forPSEPDVik,μ so that we can eventually evaluate the value υ(μ) of PLRμ for each μ. First, for given μ, Proposition 2 and Corollary 3 are characterized so as to reduce the solution space of X and (Y,Z) variables, respectively, in PSEPDVik,μ.Proposition 2.
For givenμ, the following statements hold in PSEPDVik,μ: (a)
Ifθijkp≤0 for some j∈J and p∈T, then there is an optimal solution with Xijkp=0.(b)
Ifθijkp≥hikpS⌈dijp/Γ⌉ for some j∈J and p∈T, then there is an optimal solution with Xijkp=1.Proof.
(a) Let us compare the situations ofXijkp=1 and Xijkp=0 for some j∈J and p∈T such that θijkp≤0. The value of the objective function (4a) for the former case cannot be smaller than that for the latter case because the coefficient of the X variable in (4a) is nonnegative. Also, the former case makes the constraint (4b) tighter than the latter case. Therefore, the former case (i.e., Xijkp=1) cannot give a better solution than the latter case (i.e., Xijkp=0).
(b) This can be shown in the same manner as in the proof of (a), and hence we omit its proof.LetQik be a set of freight demands that can be transported through (i,k) where i∈I and k∈K. Note that Qk=⋃i∈I‍Qik. For p∈T, let Q~ikp0 and Q~ikp1 be subsets of Qik satisfying Propositions 2(a) and 2(b) in PSEPDVik,μ, respectively, and let Q~ikp=Q~ikp0∪Q~ikp1.Corollary 3.
For givenμ and {Xijkp} fixed at the associated values according to Proposition 2, there is an optimal solution in PSEPDVik,μ with (Yik,Zikp)∈(y,zp)∣Y~ikmin⁡≤y≤Y~ikmax⁡,0≤zp≤Z~ikp,ymax⁡,∀p∈T, where (7)Y~ikmin⁡=minp∈T⁡∑j∈J∣i,j∈Q~ikp1dijpΓ,Y~ikmax⁡=maxp∈T⁡∑j∈J∣i,j∈Qik∖Q~ikp0dijpΓ,Z~ikp,ymax⁡=∑j∈J∣i,j∈Qik∖Q~ikp0dijpΓ-y.Proof.
The constraint (4b) can be rewritten as(8)∑j∈J∣i,j∈Q~ikp1dijp+∑j∈J∣i,j∈Qik∖Q~ikpdijpXijkp≤ΓYik+Zikpfor all p∈T. Since all the variables are binary variables, Y~ikmin⁡≤∑j∈J∣i,j∈Q~ikp1‍dijp/Γ≤Yik+Zikp holds for all p∈T. Since it is assumed that 0<hikL≤hikpS for all p∈T, any feasible solution with (Xijkp,Yik,Zikp)=(xjp,y,zp), where y<Y~ikmin⁡ cannot yield a better solution value than the feasible solution with (Xijkp,Yik,Zikp)=xjp,Y~ikmin⁡,max0,zp-Y~ikmin⁡+y. Moveover, any feasible solution with (Xijkp,Yik,Zikp)=(xjp,y,zp), where y>Y~ikmax⁡ cannot yield a better solution value than the feasible solution with (Xijkp,Yik,Zikp)=xjp,Y~ikmax⁡,0. Therefore, there is an optimal solution with Y~ikmin⁡≤Yik≤Y~ikmax⁡. Also, with Yik fixed at y∈y′∈Z+∣Y~ikmin⁡≤y′≤Y~ikmax⁡, it is obvious that there is an optimal solution with 0≤Zikp≤Z~ikp,ymax⁡.Then, the restricted problem withYik fixed at y∈y′∈Z+∣Y~ikmin⁡≤y′≤Y~ikmax⁡ in PSEPDVik,μ, denoted by PSEPDVik,μ,y, can be decomposed into single-period bounded knapsack problems (BKP), PBKPikp,μ,y, as follows.Problem
P
B
K
P
i
k
p
,
μ
,
y
(9)
υ
BKP
i
k
p
μ
,
y
=
max
⁡
h
i
k
p
S
Z
¯
i
k
p
+
∑
j
∈
J
∣
i
,
j
∈
Q
i
k
∖
Q
~
i
k
p
0
θ
i
j
k
p
X
i
j
k
p
subject  to
∑
j
∈
J
∣
i
,
j
∈
Q
i
k
∖
Q
~
i
k
p
d
i
j
p
X
i
j
k
p
+
Γ
Z
¯
i
k
p
≤
Γ
∑
j
∈
J
∣
i
,
j
∈
Q
i
k
∖
Q
~
i
k
p
0
d
i
j
p
Γ
-
∑
j
∈
J
∣
i
,
j
∈
Q
~
i
k
p
1
d
i
j
p
Z
¯
i
k
p
∈
0,1
,
…
,
∑
j
∈
J
∣
i
,
j
∈
Q
i
k
∖
Q
~
i
k
p
0
d
i
j
p
Γ
-
y
,
X
i
j
k
p
∈
0,1
.Note that Z¯ikp=Z~ikp,ymax⁡-Zikp. We next discuss how to solve PBKPikp,μ,y. Let υSEPDVik(μ,y) denote the optimal value of PSEPDVik,μ,y. We first transform PBKPikp,μ,y into the associated 0-1 knapsack problem, and the knapsack problem can be solved by the dynamic programming algorithm proposed by Toth [30]. We remark that Z¯ikp takes the finite number of integer values in the set Ψ≡0,1,…,∑j∈J∣i,j∈Qik∖Q~ikp0‍dijp/Γ-y. Furthermore, for each value of Z¯ikp in Ψ, it is straightforward to show that the problem PBKPikp,μ,y reduces to a 0-1 knapsack problem with decision variables Xijkp for all j. Therefore, we iteratively solve a 0-1 knapsack problem for each value of Z¯ikp in Ψ, compare the objective function values for each case, and then identify the optimal solution for PBKPikp,μ,y. Once all the problems of PBKPikp,μ,y are solved for each value of Yik as mentioned above, the optimal value of the problem PSEPDVik,μ, υSEPDVik(μ), is computed as(10)υSEPDVikμ=minY~ikmin⁡≤y≤Y~ikmax⁡⁡υSEPDVikμ,y=minY~ikmin⁡≤y≤Y~ikmax⁡⁡ThikLy+∑p∈ThikpSZ~ikp,ymax⁡-υBKPikpμ,y-∑j∈J∣i,j∈Q~ikp1θijkp.The number of problems PBKPikp,μ,y and PSEPDVik,μ,y to be solved can be further reduced by Propositions 4(a) and 4(b), respectively. Let Xijkp,y∗,Z¯ikp,y∗ be the optimal solution of PBKPikp,μ,y with explicit dependence on y, and furthermore, let us denote Try as Try=p∈T∣Z¯ikp,y∗≤Z~ikp,ymax⁡-r=Z~ikp,(y+r)max⁡ for Y~ikmin⁡≤y≤Y~ikmax⁡ and 0<r≤Z~ikp,ymax⁡. Note that Try⊂Tly⊂T for 0≤l<r≤Z~ikp,ymax⁡.Proposition 4.
(a) GivenXijkp,y∗,Z¯ikp,y∗ for some y∈y′∈Z+∣Y~ikmin⁡≤y′≤Y~ikmax⁡ and some p∈T, if p∈Try for some r∈r′∈Z+∣0<r′≤Z~ikp,ymax⁡, then Xijkp,(y+r)∗,Z¯ikp,(y+r)∗=Xijkp,y∗,Z¯ikp,y∗.
(b) GivenXijkp,y∗,Z¯ikp,y∗ for some y∈y′∈Z+∣Y~ikmin⁡≤y′≤Y~ikmax⁡ and all p∈T, if ∑p∈Try‍hikpS≥ThikL for some r∈r′∈Z+∣0<r′≤Z~ikp,ymax⁡, then υSEPDVik(μ,y+l)≥υSEPDVik(μ,y+r) for 0≤l<r.In summary, we now present the solution procedure forPSEPDVik,μ,PROC(SEPDVik,μ), in Algorithm 1 based on the whole discussion above. For reference, Figure 1 illustrates the relationship among all the problems discussed in this section.Algorithm 1:Solution procedure for PSEPDVik,μ.
P
ROC(SEPDVik,μ)
begin
(1) LetYik∗,Zikp∗ for all p∈T and Xijkp∗ for all j∈{j∈J∣(i,j)∈Qk} be Y, Z and X-variables of
the optimal solution ofPSEPDVik,μ respectively;
(2) SetυSEPDVik(μ)=∞;
(3) Fix someXijkps at the associated values according to Proposition 2 and determine Y~ikmin,
Y~ikmax and Z~ikp,ymax for all y∈{y′∈Z+∣Y~ikmin≤y′≤Y~ikmax} by Corollary 3;
(4) SetTry=ϕ for all r∈{r′∈Z+∣0<r′≤Z~ikp,ymax⁡} and y∈{y′∈Z+∣Y~ikmin⁡≤y′≤Y~ikmax⁡};
(5) Setαy=0 and βpy=0 for all y∈{y′∈Z+∣Y~ikmin⁡≤y′≤Y~ikmax⁡} and p∈T where αy and
βpy are binary variables having 1 only if PSEPDVik,μ and PBKPikp,μ,y are not necessary
to be solved respectively, and 0 otherwise;
(6) Sety=Y~ikmin⁡;
(7)whiley≤Y~ikmaxdo
(8)ifαy=0then
(9)     Setp=1;
(10)whilep≤Tdo
(11)ifβpy=0then
(12)       SolvePBKPikp,μ,y and set the resulting solution as (Xijkp,y∗,Z¯ikp,y∗);
(13)       SetTr-l,y+l=Tr-l,y+1∪{p} for all r∈{r′∈Z+∖{0}∣Z¯ikp,y∗≤Z~ikp,ymax-r′} and l∈{l′∈Z+∣0≤l′<r};
end
(14)     Set(Xijkp,y+r∗,Z¯ikp,y+r∗)=(Xijkp,y∗,Z¯ikp,y∗) and βp,y+r=1 for all
r∈{r′∈Z+∖{0}∣Z¯ikp,y∗≤Z~ikp,ymax⁡-r′} by Proposition 4(a);
(15)     Setp=p+1;
end
(16)     Setαy+l=1 for all l∈l′∈Z+∣0<l′<r∗ by Proposition 4(b) where
r∗=maxr∈Rr such that R={r∈Z+∣0<r≤Z~ikp,ymax,∑p∈TryhikpS≥ThikL};
(17)ifυSEPDVik(μ,y)<υSEPDVik(μ)then
(18)     SetυSEPDVik(μ)=υSEPDVik(μ,y) (refer to (10));
(19)     SetYik∗=y and (Xijkp∗,Z¯ikp∗)=(Xijkp,y∗,Z~ikp,ymax-Z¯ikp,y∗) for all p∈T;
end
end
(20)   Sety=y+1;
end
endFigure 1
Problem decomposition ofPMVLPUC.


## 3.3. Finding the Lagrange Multipliers

So far, we examined how to address the problemPSEPDVik,μ for each (i,k)∈EIN when the Lagrangian multiplier μ is given. We now discuss the procedure for finding the Lagrangian multiplier μ. For the purpose, we adapt the subgradient optimization procedure which is one of the most popular methods to find a good set of Lagrange multipliers [31]. Specifically, the Lagrange multipliers μ≡(μijp′,μijkp′′)∀i,j,k,p are generated iteratively as(i)
μijp,s+1′=μijp,s′+γsGijp,s′ for all (i,j)∈Qk, p∈T,(ii)
μijkp,s+1′′=max0,μijkp,s′′+γs+γsGijp,s′ for all k∈T, (i,j)∈Qk, p∈T,where Gijp,s′ and Gijkp,s′′ are the subgradients for constraints (2b) and (2d), respectively, and γs is the stepsize at iteration s.When it comes to the subgradientsGijp,s′ and Gijkp,s′′ in the formulas above, they are determined by the solution of PLRμ given μ as follows: Gijp,s′=1-∑{k∈K∣(i,j)∈Qk}‍Xijkp for all (i,j)∈Qk, p∈T, and Gijkp,s′′=Xijkp-X¯ijkp for all k∈K,(i,j)∈Qk, p∈T. Moreover, the stepsize γs at iteration s is determined by(11)γs=ηUB-υμs∑i,j∈Q∑p∈TGijp,s′2+∑k∈K∣i,j∈QkGijkp,s′′2,where UB represents the best upper bound of υOPT found up to iteration s and η is a control parameter for γs. In this paper, η is set to 2 at the beginning and then halved if the lower bound is not improved in a predetermined number of consecutive iterations.


## 3.4. Finding a Feasible Solution

A solution obtained by solvingPLRμ based on the discussion in Section 3.2 may be infeasible to PMVLPUC+, and hence we propose a two-phase heuristic to deal with infeasibility issues as follows: In the first phase, a construction heuristic (CH) is derived to modify any infeasible solutions to PMVLPUC+ into feasible ones, and then, in the second phase, a tabu search-based heuristic (TSH) is utilized to improve the solution of (CH). The details are presented in the subsequent sections.


### 3.4.1. Construction Heuristic (CH)

As mentioned above, a construction heuristic (CH) intends to deal with possible infeasibility issues of solutions obtained by solving PLRμ, and (CH) is performed at each subgradient iteration. Let Xijkp∗ and X¯ijkp∗ be the X- and X¯-variable values of the optimal solution of PSEPDVik,μ and PSEPDVkj,μ, respectively. First, the corresponding demand for all (i,j)∈Q and p∈T is reassigned to the UFCC with the lexicographical maximum value of the three-tuple relation Xijkp∗+X¯ijkp∗,θijkpXijkp∗+δijkpX¯ijkp∗,-hikL+hkjL. Then, from the reassigned demands, demand quantity Dep for all e∈E and p∈T is computed as Dep=∑{j∈J∣(i,j)∈Qk}‍dijpXijkp or Dep=∑{i∈I∣(i,j)∈Qk}‍dijpX¯ijkp. Finally, to allocate vehicles for each e∈E at the minimum vehicle cost, we solve the following vehicle allocation problem (VAP), PVAPe.Problem
P
V
A
P
e
(12a)
υ
VAP
e
=
min
⁡
T
h
e
L
Y
e
+
∑
p
∈
T
h
e
p
S
Z
e
p
(12b)
subject  to
D
e
p
≤
Γ
Y
e
+
Z
e
p
,
∀
p
∈
T
(12c)
Y
e
,
Z
e
p
∈
Z
+
,
∀
p
∈
T
.Proposition 5 characterizes the optimal solutions of PVAPe.Proposition 5.
There is an optimal solution ofPVAPe with (Ye,Zep)∈(y,Zep,y)∣y=y[1],…,y[|T|], where Zep,y=max0,⌈Dep/Γ⌉-y and y[i] is the ith largest y.Proof.
As in Corollary3, it is straightforward to show that there is an optimal solution of PVAPe with (Ye,Zep)∈(y,Zep,y)∣y[1]≤y≤y[|T|]. For an integer m such that y[n]≤m<y[n+1] and for n∈1,…,|T|-1,(13)Zep,m=0ifp≤nyp-motherwise,Zep,m+1=0ifp≤nyp-m+1otherwise.Let υVAPe(y) be the value of the objective function (12a) associated with Ye fixed at y∈[y1,y[|T|]] in PVAPe. Then, υVAPe(m+1)-υVAPe(m)=(m+1)|T|heL-∑p=1|T|‍Ze[p],m+1 − m|T|heL-∑p=1|T|‍Ze[p],m=|T|heL-∑p=n+1|T|‍he[p]S. It implies that υVAPe(m) is a linear function in m over the range y[n]≤m<y[n+1]. Therefore, there is an optimal solution of PVAPe with (Ye,Zep)∈(y,Zep,y)∣y=y[1],…,y[|T|].By Proposition5, PVAPe can be solved by finding a solution of having the minimum cost among |T| number of solutions which are associated with Ye fixed at y[p] for all p∈T.


### 3.4.2. Tabu Search-Based Heuristic (TSH)

After the first phase (CH), a tabu search-based heuristic (TSH) is performed to improve the solution obtained from (CH). We note that (TSH) is expected to require much more computational time than (CH), and hence, to save the elapsed time of the overall algorithm, (TSH) is executed only at the subgradient iterations satisfying the following two conditions ((a) and (b)) together: (a) the best lower bound is improved and (b) the solution obtained by (CH) at the current subgradient iteration is different from the solutions obtained by (CH) so far. In this paper, we call the conditions (a) and (b) as “(TSH)-activating conditions.”Tabu search is one of the popular metaheuristics designed to escape the trap of local optimality [32], and it has been used frequently in various network-related problems as in Sung and Song [25], Crevier et al. [33], and Üster and Maheshwari [34]. We now explain the main features of the tabu search used in (TSH), includingmove operation,tabu restriction andaspiration criterion,move evaluation, andtermination condition.Themove operation used in this paper is to modify assignment of some freight demand, say q=(i,j)∈Q, from its current UFCC, say l∈K, to the alternative UFCC, say k∈K, for some period, say p∈T. It corresponds to modifying the value of Xijlp from 1 to 0 and the value of Xijkp from 0 to 1 in Problem PMVLPUC+. Thetabu restriction is employed such that the UFCC assignment for any freight demand should not be changed duringtabu_size number of iterations except for the case that satisfies the aspiration criterion (to be explained later). To implement the tabu restriction, as a new solution is selected, the information of the freight demand associated with the move operation is recorded in a list (calledtabu list) and kept in the list over thetabu_size number of iterations. In this paper, thetabu_size is chosen randomly from the discrete uniform distribution in the closed interval [10,20]. Theaspiration criterion adapted here is as follows: a solution from the move operation for the freight demands recorded in thetabu list is eligible to be chosen as a new solution only if the move improves the best upper bound. Moreover, when we choose a new solution among neighbor solutions at each iteration, themove evaluation function used in this paper is to maximize {mv-A·n} where mv, n, and A represent the decrease amount of the objective function value, the frequency that the corresponding move is used to produce a new solution so far, and the associated penalty with the frequency, respectively. When it comes to mv, it is computed simply by subtracting the objective value of the new solution from the objective value of the current solution, and if mv is positive, A is set to 0 and 50 otherwise. We note that this move evaluation function is known to have an advantage in diversifying the search direction efficiently [32]. Lastly, if the total number of performed iterations amounts tomax_iter or the best upper bound has not been improved duringmax_no_imp consecutive iterations, (TSH) isterminated. In this paper,max_iter andmax_no_imp are set at 500 and 30, respectively.


## 3.4.1. Construction Heuristic (CH)

As mentioned above, a construction heuristic (CH) intends to deal with possible infeasibility issues of solutions obtained by solving PLRμ, and (CH) is performed at each subgradient iteration. Let Xijkp∗ and X¯ijkp∗ be the X- and X¯-variable values of the optimal solution of PSEPDVik,μ and PSEPDVkj,μ, respectively. First, the corresponding demand for all (i,j)∈Q and p∈T is reassigned to the UFCC with the lexicographical maximum value of the three-tuple relation Xijkp∗+X¯ijkp∗,θijkpXijkp∗+δijkpX¯ijkp∗,-hikL+hkjL. Then, from the reassigned demands, demand quantity Dep for all e∈E and p∈T is computed as Dep=∑{j∈J∣(i,j)∈Qk}‍dijpXijkp or Dep=∑{i∈I∣(i,j)∈Qk}‍dijpX¯ijkp. Finally, to allocate vehicles for each e∈E at the minimum vehicle cost, we solve the following vehicle allocation problem (VAP), PVAPe.Problem
P
V
A
P
e
(12a)
υ
VAP
e
=
min
⁡
T
h
e
L
Y
e
+
∑
p
∈
T
h
e
p
S
Z
e
p
(12b)
subject  to
D
e
p
≤
Γ
Y
e
+
Z
e
p
,
∀
p
∈
T
(12c)
Y
e
,
Z
e
p
∈
Z
+
,
∀
p
∈
T
.Proposition 5 characterizes the optimal solutions of PVAPe.Proposition 5.
There is an optimal solution ofPVAPe with (Ye,Zep)∈(y,Zep,y)∣y=y[1],…,y[|T|], where Zep,y=max0,⌈Dep/Γ⌉-y and y[i] is the ith largest y.Proof.
As in Corollary3, it is straightforward to show that there is an optimal solution of PVAPe with (Ye,Zep)∈(y,Zep,y)∣y[1]≤y≤y[|T|]. For an integer m such that y[n]≤m<y[n+1] and for n∈1,…,|T|-1,(13)Zep,m=0ifp≤nyp-motherwise,Zep,m+1=0ifp≤nyp-m+1otherwise.Let υVAPe(y) be the value of the objective function (12a) associated with Ye fixed at y∈[y1,y[|T|]] in PVAPe. Then, υVAPe(m+1)-υVAPe(m)=(m+1)|T|heL-∑p=1|T|‍Ze[p],m+1 − m|T|heL-∑p=1|T|‍Ze[p],m=|T|heL-∑p=n+1|T|‍he[p]S. It implies that υVAPe(m) is a linear function in m over the range y[n]≤m<y[n+1]. Therefore, there is an optimal solution of PVAPe with (Ye,Zep)∈(y,Zep,y)∣y=y[1],…,y[|T|].By Proposition5, PVAPe can be solved by finding a solution of having the minimum cost among |T| number of solutions which are associated with Ye fixed at y[p] for all p∈T.


## 3.4.2. Tabu Search-Based Heuristic (TSH)

After the first phase (CH), a tabu search-based heuristic (TSH) is performed to improve the solution obtained from (CH). We note that (TSH) is expected to require much more computational time than (CH), and hence, to save the elapsed time of the overall algorithm, (TSH) is executed only at the subgradient iterations satisfying the following two conditions ((a) and (b)) together: (a) the best lower bound is improved and (b) the solution obtained by (CH) at the current subgradient iteration is different from the solutions obtained by (CH) so far. In this paper, we call the conditions (a) and (b) as “(TSH)-activating conditions.”Tabu search is one of the popular metaheuristics designed to escape the trap of local optimality [32], and it has been used frequently in various network-related problems as in Sung and Song [25], Crevier et al. [33], and Üster and Maheshwari [34]. We now explain the main features of the tabu search used in (TSH), includingmove operation,tabu restriction andaspiration criterion,move evaluation, andtermination condition.Themove operation used in this paper is to modify assignment of some freight demand, say q=(i,j)∈Q, from its current UFCC, say l∈K, to the alternative UFCC, say k∈K, for some period, say p∈T. It corresponds to modifying the value of Xijlp from 1 to 0 and the value of Xijkp from 0 to 1 in Problem PMVLPUC+. Thetabu restriction is employed such that the UFCC assignment for any freight demand should not be changed duringtabu_size number of iterations except for the case that satisfies the aspiration criterion (to be explained later). To implement the tabu restriction, as a new solution is selected, the information of the freight demand associated with the move operation is recorded in a list (calledtabu list) and kept in the list over thetabu_size number of iterations. In this paper, thetabu_size is chosen randomly from the discrete uniform distribution in the closed interval [10,20]. Theaspiration criterion adapted here is as follows: a solution from the move operation for the freight demands recorded in thetabu list is eligible to be chosen as a new solution only if the move improves the best upper bound. Moreover, when we choose a new solution among neighbor solutions at each iteration, themove evaluation function used in this paper is to maximize {mv-A·n} where mv, n, and A represent the decrease amount of the objective function value, the frequency that the corresponding move is used to produce a new solution so far, and the associated penalty with the frequency, respectively. When it comes to mv, it is computed simply by subtracting the objective value of the new solution from the objective value of the current solution, and if mv is positive, A is set to 0 and 50 otherwise. We note that this move evaluation function is known to have an advantage in diversifying the search direction efficiently [32]. Lastly, if the total number of performed iterations amounts tomax_iter or the best upper bound has not been improved duringmax_no_imp consecutive iterations, (TSH) isterminated. In this paper,max_iter andmax_no_imp are set at 500 and 30, respectively.


## 3.5. Overall Solution Algorithm

We now summarize all the discussion and present the algorithmPROC(MVLPUC) for solving PMVLPUC in Algorithm 2. It can be terminated either after a predetermined number (1000 in this paper) of iterations or when η gets too small (i.e., η≤0.005 in this paper). It can also be terminated when UB-LB/LB×100 is less than 0.1%, where UB and LB represent the best upper and lower bounds searched so far, respectively.Algorithm 2:Overall solution procedure for PMVLPUC.
P
ROC(MVLPUC)
begin
(1) Initialize the necessary parameters and Lagrange multipliers (i.e.,s=1, γs=2, η=2,
LB=-∞, UB=∞, μijp,s′=1 and μijkp,s′′=1);
(2)   InitializeNILB=0 where NILB represents the number of consecutive iterations for
which the best lower bound has not been improved;
repeat
(3)     Solve the Lagrangian problemPL(μs), and let υ(μs) be the resulting objective value;
ifLB<υ(μs)then
(4)      SetLB=υ(μs) and NILB=0;
else
(5)     SetNILB=NILB+1
end
(6)    Find a feasible solution(X,X¯,Y,Z)UB,s by (CH) where (X,X¯,Y,Z)UB,s represents the
feasible solution toPMVLPUC+ at iteration s;
if (TSH)-activating conditions are satisfiedthen
(7)     Update(X,X¯,Y,Z)UB,s by (TSH);
end
(8)      LetυUB,s be the objective value of the solution (X,X¯,Y,Z)UB,s for PMVLPUC+;
(9)ifUB>υUB,sthen
(10)    SetUB=υUB,s, and set (X,X¯,Y,Z)∗=(X,X¯,Y,Z)UB,s where (X,X¯,Y,Z)∗ represents
the best feasible solution searched until iterations;
end
ifNILB=30then
(11)    Setη=η/2 and NILB=0;
end
(12)    Updateμ by the subgradient optimization procedure (refer to Section 3.3);
(13)    Sets=s+1;
untilthe aforementioned termination conditions (i.e., s>1000, η≤0.005 or (UB-LB/LB)×100<0.1) are met;
end


## 4. Numerical Experiments

This section presents the numerical results on the performance of the proposed Lagrangian heuristic algorithmPROC(MVLPUC) for PMVLPUC. The algorithmPROC(MVLPUC) was implemented in C/C++ language, and the performance ofPROC(MVLPUC) is compared with those of the following three algorithms as benchmarks: (a) the Lagrangian heuristic where only (CH) is used inPROC(MVLPUC) while (TSH) is excluded, (b) a commercial optimization softwareIBM ILOG CPLEX Optimization Studio, and (c) the algorithm where (CH) is applied at every node of the branch-and-bound tree forIBM ILOG CPLEX. The three algorithms are denoted by (LCH), (CX), and (CXCH), respectively.For the numerical experiments, a variety of problem instances were randomly generated according to the underlying network structure, the number of nodes (i.e.,|V|), and the number of time periods (i.e., |T|). Three types of underlying networks were considered in the experiments including origin-centric, equal, and destination-centric networks according to the number of origin and destination nodes. Origin-centric networks (destination-centric networks) have larger numbers of origin (destination) nodes than those of destination (origin) nodes, and equal networks have almost equivalent numbers of origin nodes as those of destination nodes. Moreover, three |V| values including 30, 40, and 50 and three |T| including 3, 6, and 12 were considered in the experiment.To generate problem instances, each node was randomly selected on a100×100 grid. Some of the nodes were connected by arcs, which were also chosen randomly. Unit-period long-term vehicle lease cost heL was calculated based on Euclidean distance between nodes, and unit-period short-term vehicle lease cost hepS was randomly chosen in the range from 1.5 to 2 times of heL. Such a structure of vehicle lease costs reflects the actual practice in Korean logistics industry. Moreover, transportation time te was randomly chosen in the range from 0.8 to 1.2 times of heL, and handling time at k∈K, ok, was set as ok=minmini∈Itik,minj∈Jtkj. Transportation time restriction TLij was set as TLij=0.1mink∈Ktik+ok+tkj+0.9maxk∈Ktik+ok+tkj. Demand quantity dijp was randomly selected in the range from 1 to 90, and vehicle capacity Γ was set at 100. For each combination of underlying network structure, number of nodes, and number of periods, four problem instances were generated, and hence, in total, 108 problem instances were generated.Tables1–3 show the experimental results for the problem instances with |V|=30,40 and 50 nodes, respectively. The tables present the best lower bound (LB), the best upper bound (UB), the percentage deviation between the two (%Gap), and the elapsed time in seconds (Time[s]). The percentage deviation, %Gap, is computed as UB-LB/LB×100.Table 1
Test results of the proposed algorithm for the problem instances with|V|=30.
I
K
J
T
=
3
T
=
6
T
=
12
LB
UB
%Gap
Time [s]
LB
UB
%Gap
Time [s]
LB
UB
%Gap
Time [s]
20
5
5
212996.81
224714
5.50
53.16
461155.84
489788
6.21
98.66
929519.50
1002248
7.82
230.75
20
5
5
216811.19
241662
11.46
65.11
460583.37
504220
9.47
107.86
948275.80
1025174
8.11
255.70
22
3
5
279876.66
298848
6.78
240.95
573183.77
611568
6.70
461.84
1192304.31
1289542
8.16
691.83
22
3
5
272613.72
294054
7.86
212.24
554945.75
600348
8.18
432.34
1170739.75
1265118
8.06
419.22
13
5
12
332774.01
355556
6.85
57.92
683250.91
728424
6.61
127.08
1434397.96
1536928
7.15
272.66
13
5
12
304624.33
340794
11.87
48.02
650095.80
714146
9.85
112.13
1343684.15
1470888
9.47
296.05
13
3
14
413042.22
440862
6.74
50.27
872661.27
943590
8.13
109.41
1819085.70
1956382
7.55
325.81
13
3
14
408811.15
439939
7.61
58.99
850207.47
916393
7.78
132.24
1768298.40
1911291
8.09
302.02
5
5
20
243255.92
257018
5.66
42.53
532275.53
565122
6.17
75.63
1088310.70
1168786
7.39
191.48
5
5
20
211330.59
238778
12.99
41.20
510637.65
562726
10.20
96.13
1066249.72
1139144
6.84
211.16
5
3
22
268067.28
282174
5.26
33.77
547072.94
578438
5.73
75.20
1129752.50
1206642
6.81
120.70
5
3
22
287200.52
320211
11.49
56.61
563613.65
613835
8.91
68.11
1188680.85
1282241
7.87
129.75
Average
8.34
80.06
7.83
158.05
7.78
287.26Table 2
Test results of the proposed algorithm for the problem instances with|V|=40.
I
K
J
T
=
3
T
=
6
T
=
12
LB
UB
%Gap
Time [s]
LB
UB
%Gap
Time [s]
LB
UB
%Gap
Time [s]
26
7
7
424766.97
471054
10.90
183.97
1414851.73
1522332
7.60
409.52
2868060.92
3075664
7.24
1393.03
26
7
7
360749.62
399556
10.76
215.84
751912.45
836474
11.25
387.16
1539449.29
1722836
11.91
952.08
29
4
7
418927.58
441254
5.33
105.27
865273.19
934006
7.94
209.19
1771774.46
1933194
9.11
485.31
29
4
7
485306.04
515006
6.12
170.86
1016416.93
1075553
5.82
263.59
2071972.70
2216043
6.95
625.59
17
7
16
550512.84
608016
10.45
200.25
1180080.13
1298552
10.04
655.23
2426278.35
2681688
10.53
2074.97
17
7
16
528692.03
585014
10.65
190.52
1099526.11
1211485
10.18
502.03
2236024.28
2462131
10.11
1370.17
18
4
18
734759.20
783204
6.59
412.23
1587132.99
1749738
10.25
1085.08
3176697.85
3560520
12.08
3414.08
18
4
18
647147.85
700289
8.21
312.63
1383698.67
1485373
7.35
665.97
2769306.69
2975719
7.45
1468.78
7
7
26
376082.83
409574
8.91
190.22
825644.77
894690
8.36
477.67
1678824.58
1874644
11.66
1249.97
7
7
26
387078.61
419447
8.36
146.84
784096.19
854760
9.01
495.31
1594309.42
1748548
9.67
1147.27
7
4
29
493481.06
515596
4.48
129.53
1026893.95
1077740
4.95
281.89
2091129.61
2201314
5.27
644.69
7
4
29
527232.35
568513
7.83
121.63
1103782.36
1189023
7.72
259.09
2240109.02
2419115
7.99
667.61
Average
8.22
198.32
8.37
474.31
9.17
1291.13Table 3
Test results of the proposed algorithm for the problem instances withV=50.
I
K
J
T
=
3
T
=
6
T
=
12
LB
UB
%Gap
Time [s]
LB
UB
%Gap
Time [s]
LB
UB
%Gap
Time [s]
30
10
10
504429.05
547776
8.59
603.09
1044268.14
1151346
10.25
2322.48
2138696.38
2348434
9.81
3356.77
30
10
10
601784.87
679393
12.90
588.67
1083093.64
1227322
13.32
1893.45
2069257.08
2307938
11.53
3370.17
34
6
10
619464.52
669098
8.01
531.14
1308349.43
1417378
8.33
1819.80
2636765.90
2854142
8.24
4683.91
34
6
10
603913.87
664708
10.07
485.77
1267603.33
1396176
10.14
1293.69
2626944.97
2887162
9.91
5836.92
20
10
20
810314.00
875940
8.10
1041.22
1671447.84
1809086
8.23
2197.13
3393505.82
3669690
8.14
5413.17
20
10
20
822869.60
892081
8.41
933.94
1669648.98
1843641
10.42
2719.97
3427822.66
3727153
8.73
5259.75
22
6
22
1056175.36
1146126
8.52
722.27
2185505.69
2372490
8.56
2221.61
4404884.25
4825116
9.54
4340.27
22
6
22
937230.54
1019646
8.79
938.88
2215333.07
2404507
8.54
1831.72
3958931.06
4289352
8.35
6382.42
10
10
30
651651.93
701858
7.70
839.84
1322310.58
1428274
8.01
1842.22
2753031.46
2960746
7.54
4175.56
10
10
30
669202.65
737403
10.19
932.56
1381848.96
1505385
8.94
2470.44
2829709.60
3071734
8.55
4343.27
10
6
34
685536.17
755326
10.18
592.84
1447676.89
1591850
9.96
1656.83
2938744.04
3228866
9.87
6059.58
10
6
34
685081.79
764624
11.61
917.94
1431087.35
1595460
11.49
3033.92
2912678.41
3215532
10.40
5632.03
Average
9.42
760.68
9.68
2108.60
9.22
4904.48For the problem instances with|V|=30,PROC(MVLPUC) gave solutions with the following performance: (a) the duality gap (shown in the “%Gap” column) ranges from 5.26% to 12.99% with the average of 7.98%, and (b) the elapsed time ranges from 33.77 s to 691.83 s with the average of 175.12 s. The duality gap does not seem to be dependent on both |T| and the network structure, as observed that the duality gap does not change consistently according to them. On the other hand, the elapsed time seems to increase not exponentially but linearly as |T| increases while it does not seem to be dependent on the network structure. For example, the average elapsed times for the problem instances with |T|=3,6, and 12 are 80.06 s, 158.05 s, and 287.26 s, respectively. The aforementioned trends can be observed in the problem instances with |V|=40 and 50 as well.From all the experiments, the performance changes of the proposed algorithm according to the number of nodes (i.e.,|V|) can be summarized as follows. As |V| increases, the duality gap seems to marginally increase. For example, the average duality gaps for the problem instances with |V|=30,40, and 50 are 7.98%, 8.58%, and 9.44%, respectively, and this may be due to the fact that, as |V| increases, the number of the associated Lagrange multipliers greatly increases as shown in μijkp′′, and hence, the associated lower bound may become looser accordingly. We also observe that, as |V| increases, the elapsed time increases rapidly, but not exponentially. For example, the average elapsed times for the problem instances with |V|=30,40, and 50 are 175.12 s, 654.59 s, and 2591.26 s, respectively.We now compare the performance ofPROC(MVLPUC) with those of the aforementioned three algorithms including (LCH), (CX), and (CXCH), respectively. They were tested on the small-sized problem instances (or instances with |V|=30). For (LCH), the same termination conditions as inPROC(MVLPUC) were used, while for (CX) and (CXCH), an execution time limit of five hours was applied as a termination condition.Table4 shows the numerical results of the compared algorithms for the problem instances with |V|=30. The table gives the percentage deviations between the best lower bounds obtained from each of the compared algorithms and the best lower bound obtained by the proposed algorithm (LB in Table 1) and the percentage deviations between the best upper bounds obtained by each of the compared algorithms and the best upper bound obtained by the proposed algorithm (UB in Table 1). For (LCH) and both of (CX) and (CXCH), the percentage deviations of the lower bounds are denoted by ILBLCH and ILBCX, respectively. Note that (CX) and (CXCH) give the same lower bound obtained from the LP relaxation of Problem PMVLPUC. For (LCH), (CX), and (CXCH), the percentage deviations of the upper bounds are denoted by IUBLCH, IUBCX, and IUBCXCH, respectively. In detail, the percentage deviations of the lower bounds and the percentage deviations of the upper bounds were computed by LB-LBC/LBC×100 and UBC-UB/UBC×100, respectively, where LBC and UBC are the best lower bound and the best upper bound obtained by the corresponding compared algorithm, respectively. For example, if IUBLCH is positive, it means that the best upper bound of the proposed algorithm is better than that of (LCH) at the associated rate. From the table, the proposed algorithm may give the almost same lower bounds as (LCH) and better upper bounds than (LCH) from the fact that ILBLCH is almost zero and ULBLCH ranges from 0.36% to 11.88%. This may imply that (TSH) embedded inPROC(MVLPUC) improves the upper bound significantly. Moreover,PROC(MVLPUC) may give better lower bounds than (CX) and (CXCH). For the long planning horizon,PROC(MVLPUC) may give better upper bounds than (CX) and (CXCH), while for the short planning horizon, it may not. For example, the associated average values of IUBCX with |T|=3,6, and 12 are -2.35%, 2.56%, and 6.44%, respectively. This may imply that the proposed algorithm is suitable for the large-sized problem instances. Also, it can be observed that (CXCH) outperforms (CX) marginally.Table 4
Test results of the compared algorithms for the problem instances withV=30.
I
K
J
T
=
3
T
=
6
T
=
12
IL
B
L
C
H
I
L
B
C
X
I
U
B
L
C
H
I
U
B
C
X
I
U
B
C
X
C
H
I
L
B
L
C
H
I
L
B
C
X
I
U
B
L
C
H
I
U
B
C
X
I
U
B
C
X
C
H
I
L
B
L
C
H
I
L
B
C
X
I
U
B
L
C
H
I
U
B
C
X
I
U
B
C
X
C
H
20
5
5
0.01
10.58
5.52
-
0.17
-
1.17
-
0.19
11.45
3.73
1.38
1.16
-
0.09
9.37
3.53
4.14
0.26
20
5
5
0.00
11.14
3.90
-
6.46
-
6.35
-
0.14
10.38
4.68
-
2.77
-
2.27
-
0.08
10.21
4.66
3.61
-
0.32
22
3
5
0.45
9.55
11.88
2.30
-
0.43
0.00
9.18
8.11
10.27
2.05
0.09
9.22
4.08
7.54
4.11
22
3
5
-
0.01
8.48
4.00
-
3.38
-
0.76
0.02
8.88
6.29
-
0.47
-
1.72
0.04
8.65
3.97
9.44
4.98
13
5
12
-
0.14
10.05
2.80
1.12
-
1.10
-
0.13
9.25
7.05
3.18
1.54
-
0.71
8.34
3.92
6.77
3.88
13
5
12
-
0.31
9.50
5.78
-
2.84
-
2.72
-
0.56
8.76
6.76
1.07
1.87
0.11
9.65
6.44
7.89
6.85
13
3
14
-
0.26
5.17
5.96
1.13
0.01
-
0.17
5.52
4.07
8.08
3.13
-
0.09
5.60
5.75
11.74
7.45
13
3
14
-
0.20
6.19
5.09
-
1.86
-
1.45
-
0.07
6.33
5.61
8.67
5.70
-
0.08
6.29
6.29
9.87
8.04
5
5
20
-
0.05
15.65
4.44
-
1.48
-
1.26
-
0.56
16.84
5.09
1.00
0.12
-
0.26
16.59
3.19
0.36
1.32
5
5
20
-
0.02
12.99
0.36
-
8.46
-
8.34
0.00
16.94
3.87
-
5.58
-
5.31
-
0.16
16.57
4.13
1.42
0.17
5
3
22
-
0.06
13.93
2.91
-
1.57
-
1.87
0.02
13.92
3.82
0.99
-
0.39
-
0.58
12.42
3.86
4.14
6.70
5
3
22
-
0.01
13.45
3.48
-
6.51
-
6.97
-
0.16
12.66
3.17
-
2.62
-
2.75
-
0.68
11.69
2.96
5.66
3.84
Average
-
0.05
10.56
4.68
-
2.35
-
2.70
-
0.18
10.39
5.13
2.56
0.69
-
0.16
9.97
4.47
6.44
4.32In summary, the proposed algorithmPROC(MVLPUC) may give good lower and upper bounds within reasonable time, and the duality gap seems not to be dependent on the network structure and the number of time periods, while it seems to increase a little, as |V| increases. Furthermore, as presented in Table 3, the average duality gap for the problem instances with |V|=50 is still less than 10%, and the elapsed time seems not to be dependent on the network structure, while it increases as |V| or |T| increases. However, the increase in the elapsed time is not exponential in either case.


## 5. Concluding Remarks

In this paper, a multiperiod vehicle lease planning problem in an urban freight distribution network with UFCCs is considered where vehicles for each direct service are leased in short term or long term. The problem under consideration is motivated by the realistic situation where the freight demands fluctuate dynamically in time, and urban freight distribution can be made through UFCCs. The problem is formulated as an integer programming model and can be proven to be NP-hard in a strong sense. Thus, we propose a Lagrangian-based heuristic to address the concomitant problem. Various numerical experiments are performed to evaluate the performance of the proposed algorithm, and they indicate that the proposed algorithm works well. The proposed algorithm may be applied to vehicle lease planning for industrial logistic services for handling various products including food, apparel, electronic goods, and logistics companies.As extended research, several models derived from the proposed model in this paper can be considered. For example, a problem where each freight demand can be delivered through more than a single UFCC may be interesting. Moreover, it may be necessary to deal with vehicle specific lease periods. Lastly, proposing the efficient heuristic-based approach for solving large-scale problems would be the topic for the future research.



---

*Source: 921482-2015-10-19.xml*