# Multimodal Green Logistics Network Design of Urban Agglomeration with Stochastic Demand

Authors: Jiehui Jiang; Dezhi Zhang; Shuangyan Li; Yajie Liu
Journal: Journal of Advanced Transportation
 (2019)

---

## Abstract

This study investigates a multimodal green logistics network design problem of urban agglomeration with stochastic demand, in which different logistics authorities among the different cities jointly optimize the logistics node configurations and uniform carbon taxes over logistics transport modes to maximize the total social welfare of urban agglomeration and consider logistics users’ choice behaviors. The users’ choice behaviors are captured by a logit-based stochastic equilibrium model. To describe the game behaviors of logistics authorities in urban agglomeration, the problem is formulated as two nonlinear bilevel programming models, namely, independent and centralized decision models. Next, a quantum-behaved particle swarm optimization (QPSO) embedded with a Method of Successive Averages (MSA) is presented to solve the proposed models. Simulation results show that to achieve the overall optimization layout of the green logistics network in urban agglomeration the logistics authorities should adopt centralized decisions, construct a multimode logistics network, and make a reasonable carbon tax.


---

## Body

## 1. Introduction

The logistics network design problem (LNDP) has been widely studied [1–3]. Generally, this problem comprises two subproblems: (i) a location problem, namely, how to decide the locations of logistics nodes (such as logistics parks, distribution centers, and logistics terminals), and (ii) an allocation problem, namely, how to route the flow of goods to origin–destination (O-D). Since the construction of urban logistics nodes has the characteristics of a large investment, long period, and great risk, the rationality and feasibility of the plan need to be ensured with the objective guidance of theory and method. Meanwhile, with the development of society and the improvement of people’s living standards, traffic congestion and transportation-related environmental issues have become the focus of attention for scholars and governments [4–6]. Against this unsettling backdrop, the traditional single mode of transport does not meet the needs of modern cities. Therefore, how to design a multimodal logistics network with high efficiency, safety, and environmental friendliness is an important issue to be solved urgently in the development of urban agglomeration logistics.To locate logistics nodes, as addressed in the conventional facility location problem (FLP) [7], the multimodal logistics network design problem (MLNDP) of urban agglomeration needs to address the rational allocation of limited resources between different cities in urban agglomerations. Because of the limited resources in urban agglomerations, there is a great game among each local government in the resources (e.g., politics, economy, and population), and this game will derive many game-disordered problems under the restriction of the government’s political and economic behavior. Local governments usually establish the logistics network based on the situation of their own precinct to fit the regional development with little coordination and communication with other surrounding regions, resulting in the lack of effective integration of the basic logistics resources [8]. Therefore, game theory is used to analyze the game relationship among city governments in urban agglomeration, which is of great significance to the rational allocation of resources.To route the flow of goods, multimodal transportation is a hybrid of transportation that involves two classic transportation services: the single mode of transport and combined transport [9]. Compared with the single transportation mode, the multimodal transportation can make full use of existing logistics networks of infrastructure and the advantage of various modes of transport, integrating transport capacity resources, meeting the transport demand, and achieving the target (e.g., time, cost, and profit) [10, 11]. In the multimodal transportation, there is no effective connection between each mode of transportation, and the design of the multimodal transportation plan needs to be further improved and optimized [12]. However, urban traffic congestion and transportation-related environmental issues (e.g., carbon dioxide (CO2), nitrogen oxide, and sulfur oxide) are the bottlenecks that restrict the sustainable development of cities [13, 14]. Consequently, there is an inevitable trend toward the future development of urban agglomeration green logistics systems to develop multimodal and integrative transportation networks.


### 1.1. Literature Review

We classify the literature related to our research into three branches: (1) urban agglomeration logistics, (2) green logistics network design for facility location, and (3) logistics network design with choice behaviors.


#### 1.1.1. Urban Agglomeration Logistics

The urban agglomeration logistics system is an important part of the urban agglomeration economy as it plays an important and supporting role in the sustainable development and evolution of the urban agglomeration. Mullen and Marsden [15] explored the adaptive adjustment mechanism of regional logistics and regional economy. Li et al. [16] emphasized the relationship between traffic investment and traffic efficiency in urban agglomeration. Jiang et al. [17] analyzed the bidirectional relationship between multimodal transportation investment and economic development. However, some scholars studied the spatial and temporal distribution characteristics of urban agglomeration logistics [18, 19]. Lindsey et al. [20] studied the relationships among freight transport, economic market drivers, and industrial space demand. Kumar et al. [21] investigated spatial patterns of transportation and logistics cluster in the US regions, applying spatial cluster and econometric analyses. The above studies mainly focus on the linkage mechanism between logistics and economy, demand and spatial and temporal distribution, and few focus on the green logistics system of urban agglomeration. Table 1 summarizes the main features of the proposed models in the reviewed studies.Table 1
Characteristics of the reviewed models proposed for the green logistics network design problem.
Study
Model characteristics
Objective functions
Decision variables
Solution procedure
Research areas
decision makers
Transport modes
Demands
Choice behaviors
Number
Type
Capacity
Routing
Carbon tax
Exact
Heuristic
S
M
S
M
S
M
De
St
Type
Utility
S
M
Xiong and Wang [12]
✓
✓
✓
✓
-
-
✓
Cost, time
✓
✓
MOGA
Harris et al. [25]
✓
✓
✓
✓
-
-
✓
Cost, CO2 emissions
✓
✓
SEAMO2
Yang et al. [26]
✓
✓
✓
✓
-
-
✓
Cost
✓
✓
Cplex
Turken et al. [24]
✓
✓
✓
✓
-
-
✓
Cost
✓
✓
FLAA
Rezaee et al. [2]
✓
✓
✓
✓
-
-
✓
Cost
✓
✓
Cplex
Gao et al. [28]
✓
✓
✓
✓
-
-
✓
Cost
✓
✓
Cplex
Yamada et al. [37]
✓
✓
✓
✓
UE
Cost, time
✓
Benefit–cost ratio
✓
✓
GA,GLS
Meng and Wang [38]
✓
✓
✓
✓
UE
Cost, time
✓
Cost
✓
✓
GA
Wang and Meng [39]
✓
✓
✓
✓
UE
Cost, time
✓
Cost
✓
✓
B&B
Zhang et al. [35]
✓
✓
✓
✓
UE
Cost, time
✓
Benefit–cost ratio
✓
✓
GA+ FWA
This study
✓
✓
✓
✓
SUE
Cost, time, CO2 emission
✓
Social benefit
✓
✓
✓
QPSO+ MSA
S: single, M: multiple, De: deterministic, St: stochastic, UE: user equilibrium, SUE: stochastic user equilibrium, GA: Genetic Algorithm, and MOGA: Multiobjective Genetic Algorithm.
SEAMO2: Simple Evolutionary Multiobjective Optimization 2, FLAA: Facility Linear Approximation Algorithm, GLS: Genetic Local Search, and B&B: Branch-and-Bound.
FW: Frank–Wolfe Algorithm, QPSO: Quantum behaved Particle Swarm Optimization, and MSA: Method of Successive Averages.


#### 1.1.2. Green Logistics Network Design for Facility Location

Green logistics aims to restrain the environmental damage (e.g., greenhouse gas emissions, noise, and accidents) caused by logistics activity and develop a sustainable balance between economic, social, and environmental objectives [22]. Facility location modeling is a strategic planning design approach that selects the optimal set of facilities from a set of potential facilities [23]. As a branch of network design, the green logistics network design problem originates from the classic facility location problem (FLP) [7, 24] and then extends it to the location-allocation problem (LAP) [25–27].In the background of the vigorously developing green economy, Turken et al. [24] investigated the impact of carbon tax on plant capacity and location decisions of a firm. Yang et al. [26] studied the low-carbon network design problem of the third-party logistics distribution system. To solve location-allocation problem with a carbon emissions constraint, Rezaee et al. [2] developed a two-stage stochastic programming model, and Gao et al. [28] proposed five multidimensional mixed-integer nonlinear programming models. In practice, the urban agglomeration logistics network may comprise multiple urban logistics networks, which are independently and separately managed by local governments (or logistics authorities) with different targets. Nevertheless, the literature always supposed that the network is totally designed by a single management body, which may not suit the design of the logistics network in urban agglomeration.


#### 1.1.3. Green Logistics Network Design under Route Choice Models

In general, route choice models should evaluate the utility of each route and route the flow of passengers or goods from origin to destination. Although they have been widely used to solve urban passenger transportation problems by researchers and engineers [29–33], using these models to solve urban cargo transportation problems has received little attention in the literature.The logistics network design problem based on route choice behavior can be characterized as a bilevel programming model or an equivalent mathematical program with equilibrium constraints (MPEC) model. At the upper level, the logistics authority determines the number, location, and capacity of logistics facilities (nodes). At the lower level, shippers and carriers follow the user equilibrium (UE) principle or stochastic user equilibrium (SUE) principle [34–36]. Yamada et al. [37] presented a bilevel programming model to simulate the multimodal freight transport network design problem, where the lower-level model is a multimodal multiclass user equilibrium model. To address the intermodal hub-and-spoke network design problem for multiple stakeholders and multitype containers, Meng and Wang [38] established an MPEC model and used a variational inequality to describe the operators’ route choice behavior. Wang and Meng [39] extended the research work of Meng and Wang [38] by considering congestion effects and piecewise linear cost functions and solved it with a solution algorithm based on nonlinear optimization and branch-and-bound global optimization.However, the route choice models of these studies did not consider the carbon emissions cost, perception error, and elastic demand. In addition, they mainly focus on static, determined network design issues. As strategic decision planning, the design of a logistics network in urban agglomeration involves long-term planning of a between 5- and 20-year timeline, where the logistics demand is often uncertain with the development of urban evolution and industrial structure optimization. Hence, it is more realistic and reasonable to design the logistics network with demand uncertainty and stochastic route choice.


### 1.2. Objectives and Contributions

In view of the abovementioned realistic problems and specific characteristics of logistics development in urban agglomerations, the objective of this paper is to coordinate the distribution of various types of logistics nodes (e.g., logistics parks, distribution centers and logistics terminals) in the urban agglomeration logistics network from the overall perspective of urban agglomeration and fully consider the low carbon requirement together with demand uncertainties. The main contributions of this paper are as follows. First, to characterize different decision-making behaviors among multiple local authorities, the multimodal logistics network design problem is formulated as two nonlinear bilevel programming models. At the upper level, the logistics authority of each city attempts to maximize the social welfare including producer surplus and consumer surplus. At the lower level, the logistics users’ route choice decisions follow logit-based stochastic user equilibrium (SUE) with elastic demand under logistics demand scenario. Second, a heuristic solution algorithm that is a combination of quantum behaved particle swarm optimization (QPSO) and Method of Successive Averages (MSA) is developed to solve the proposed bilevel programming model. Third, the optimal number, scale, and location of logistics nodes, analysis of the impact on the distribution of city group logistics network of infrastructure investment budget, and carbon tax are illustrated with an example.The remainder of this article is arranged as follows. The basic considerations of this paper, including general assumptions and network representation, are described in Section2. A nonlinear bilevel programming formulation is proposed in Section 3, and the solution methods are presented in Section 4. A numerical example is presented to illustrate the availability of previous models in Section 5. Finally, the conclusion and future studies are discussed in Section 6.


## 1.1. Literature Review

We classify the literature related to our research into three branches: (1) urban agglomeration logistics, (2) green logistics network design for facility location, and (3) logistics network design with choice behaviors.


### 1.1.1. Urban Agglomeration Logistics

The urban agglomeration logistics system is an important part of the urban agglomeration economy as it plays an important and supporting role in the sustainable development and evolution of the urban agglomeration. Mullen and Marsden [15] explored the adaptive adjustment mechanism of regional logistics and regional economy. Li et al. [16] emphasized the relationship between traffic investment and traffic efficiency in urban agglomeration. Jiang et al. [17] analyzed the bidirectional relationship between multimodal transportation investment and economic development. However, some scholars studied the spatial and temporal distribution characteristics of urban agglomeration logistics [18, 19]. Lindsey et al. [20] studied the relationships among freight transport, economic market drivers, and industrial space demand. Kumar et al. [21] investigated spatial patterns of transportation and logistics cluster in the US regions, applying spatial cluster and econometric analyses. The above studies mainly focus on the linkage mechanism between logistics and economy, demand and spatial and temporal distribution, and few focus on the green logistics system of urban agglomeration. Table 1 summarizes the main features of the proposed models in the reviewed studies.Table 1
Characteristics of the reviewed models proposed for the green logistics network design problem.
Study
Model characteristics
Objective functions
Decision variables
Solution procedure
Research areas
decision makers
Transport modes
Demands
Choice behaviors
Number
Type
Capacity
Routing
Carbon tax
Exact
Heuristic
S
M
S
M
S
M
De
St
Type
Utility
S
M
Xiong and Wang [12]
✓
✓
✓
✓
-
-
✓
Cost, time
✓
✓
MOGA
Harris et al. [25]
✓
✓
✓
✓
-
-
✓
Cost, CO2 emissions
✓
✓
SEAMO2
Yang et al. [26]
✓
✓
✓
✓
-
-
✓
Cost
✓
✓
Cplex
Turken et al. [24]
✓
✓
✓
✓
-
-
✓
Cost
✓
✓
FLAA
Rezaee et al. [2]
✓
✓
✓
✓
-
-
✓
Cost
✓
✓
Cplex
Gao et al. [28]
✓
✓
✓
✓
-
-
✓
Cost
✓
✓
Cplex
Yamada et al. [37]
✓
✓
✓
✓
UE
Cost, time
✓
Benefit–cost ratio
✓
✓
GA,GLS
Meng and Wang [38]
✓
✓
✓
✓
UE
Cost, time
✓
Cost
✓
✓
GA
Wang and Meng [39]
✓
✓
✓
✓
UE
Cost, time
✓
Cost
✓
✓
B&B
Zhang et al. [35]
✓
✓
✓
✓
UE
Cost, time
✓
Benefit–cost ratio
✓
✓
GA+ FWA
This study
✓
✓
✓
✓
SUE
Cost, time, CO2 emission
✓
Social benefit
✓
✓
✓
QPSO+ MSA
S: single, M: multiple, De: deterministic, St: stochastic, UE: user equilibrium, SUE: stochastic user equilibrium, GA: Genetic Algorithm, and MOGA: Multiobjective Genetic Algorithm.
SEAMO2: Simple Evolutionary Multiobjective Optimization 2, FLAA: Facility Linear Approximation Algorithm, GLS: Genetic Local Search, and B&B: Branch-and-Bound.
FW: Frank–Wolfe Algorithm, QPSO: Quantum behaved Particle Swarm Optimization, and MSA: Method of Successive Averages.


### 1.1.2. Green Logistics Network Design for Facility Location

Green logistics aims to restrain the environmental damage (e.g., greenhouse gas emissions, noise, and accidents) caused by logistics activity and develop a sustainable balance between economic, social, and environmental objectives [22]. Facility location modeling is a strategic planning design approach that selects the optimal set of facilities from a set of potential facilities [23]. As a branch of network design, the green logistics network design problem originates from the classic facility location problem (FLP) [7, 24] and then extends it to the location-allocation problem (LAP) [25–27].In the background of the vigorously developing green economy, Turken et al. [24] investigated the impact of carbon tax on plant capacity and location decisions of a firm. Yang et al. [26] studied the low-carbon network design problem of the third-party logistics distribution system. To solve location-allocation problem with a carbon emissions constraint, Rezaee et al. [2] developed a two-stage stochastic programming model, and Gao et al. [28] proposed five multidimensional mixed-integer nonlinear programming models. In practice, the urban agglomeration logistics network may comprise multiple urban logistics networks, which are independently and separately managed by local governments (or logistics authorities) with different targets. Nevertheless, the literature always supposed that the network is totally designed by a single management body, which may not suit the design of the logistics network in urban agglomeration.


### 1.1.3. Green Logistics Network Design under Route Choice Models

In general, route choice models should evaluate the utility of each route and route the flow of passengers or goods from origin to destination. Although they have been widely used to solve urban passenger transportation problems by researchers and engineers [29–33], using these models to solve urban cargo transportation problems has received little attention in the literature.The logistics network design problem based on route choice behavior can be characterized as a bilevel programming model or an equivalent mathematical program with equilibrium constraints (MPEC) model. At the upper level, the logistics authority determines the number, location, and capacity of logistics facilities (nodes). At the lower level, shippers and carriers follow the user equilibrium (UE) principle or stochastic user equilibrium (SUE) principle [34–36]. Yamada et al. [37] presented a bilevel programming model to simulate the multimodal freight transport network design problem, where the lower-level model is a multimodal multiclass user equilibrium model. To address the intermodal hub-and-spoke network design problem for multiple stakeholders and multitype containers, Meng and Wang [38] established an MPEC model and used a variational inequality to describe the operators’ route choice behavior. Wang and Meng [39] extended the research work of Meng and Wang [38] by considering congestion effects and piecewise linear cost functions and solved it with a solution algorithm based on nonlinear optimization and branch-and-bound global optimization.However, the route choice models of these studies did not consider the carbon emissions cost, perception error, and elastic demand. In addition, they mainly focus on static, determined network design issues. As strategic decision planning, the design of a logistics network in urban agglomeration involves long-term planning of a between 5- and 20-year timeline, where the logistics demand is often uncertain with the development of urban evolution and industrial structure optimization. Hence, it is more realistic and reasonable to design the logistics network with demand uncertainty and stochastic route choice.


## 1.1.1. Urban Agglomeration Logistics

The urban agglomeration logistics system is an important part of the urban agglomeration economy as it plays an important and supporting role in the sustainable development and evolution of the urban agglomeration. Mullen and Marsden [15] explored the adaptive adjustment mechanism of regional logistics and regional economy. Li et al. [16] emphasized the relationship between traffic investment and traffic efficiency in urban agglomeration. Jiang et al. [17] analyzed the bidirectional relationship between multimodal transportation investment and economic development. However, some scholars studied the spatial and temporal distribution characteristics of urban agglomeration logistics [18, 19]. Lindsey et al. [20] studied the relationships among freight transport, economic market drivers, and industrial space demand. Kumar et al. [21] investigated spatial patterns of transportation and logistics cluster in the US regions, applying spatial cluster and econometric analyses. The above studies mainly focus on the linkage mechanism between logistics and economy, demand and spatial and temporal distribution, and few focus on the green logistics system of urban agglomeration. Table 1 summarizes the main features of the proposed models in the reviewed studies.Table 1
Characteristics of the reviewed models proposed for the green logistics network design problem.
Study
Model characteristics
Objective functions
Decision variables
Solution procedure
Research areas
decision makers
Transport modes
Demands
Choice behaviors
Number
Type
Capacity
Routing
Carbon tax
Exact
Heuristic
S
M
S
M
S
M
De
St
Type
Utility
S
M
Xiong and Wang [12]
✓
✓
✓
✓
-
-
✓
Cost, time
✓
✓
MOGA
Harris et al. [25]
✓
✓
✓
✓
-
-
✓
Cost, CO2 emissions
✓
✓
SEAMO2
Yang et al. [26]
✓
✓
✓
✓
-
-
✓
Cost
✓
✓
Cplex
Turken et al. [24]
✓
✓
✓
✓
-
-
✓
Cost
✓
✓
FLAA
Rezaee et al. [2]
✓
✓
✓
✓
-
-
✓
Cost
✓
✓
Cplex
Gao et al. [28]
✓
✓
✓
✓
-
-
✓
Cost
✓
✓
Cplex
Yamada et al. [37]
✓
✓
✓
✓
UE
Cost, time
✓
Benefit–cost ratio
✓
✓
GA,GLS
Meng and Wang [38]
✓
✓
✓
✓
UE
Cost, time
✓
Cost
✓
✓
GA
Wang and Meng [39]
✓
✓
✓
✓
UE
Cost, time
✓
Cost
✓
✓
B&B
Zhang et al. [35]
✓
✓
✓
✓
UE
Cost, time
✓
Benefit–cost ratio
✓
✓
GA+ FWA
This study
✓
✓
✓
✓
SUE
Cost, time, CO2 emission
✓
Social benefit
✓
✓
✓
QPSO+ MSA
S: single, M: multiple, De: deterministic, St: stochastic, UE: user equilibrium, SUE: stochastic user equilibrium, GA: Genetic Algorithm, and MOGA: Multiobjective Genetic Algorithm.
SEAMO2: Simple Evolutionary Multiobjective Optimization 2, FLAA: Facility Linear Approximation Algorithm, GLS: Genetic Local Search, and B&B: Branch-and-Bound.
FW: Frank–Wolfe Algorithm, QPSO: Quantum behaved Particle Swarm Optimization, and MSA: Method of Successive Averages.


## 1.1.2. Green Logistics Network Design for Facility Location

Green logistics aims to restrain the environmental damage (e.g., greenhouse gas emissions, noise, and accidents) caused by logistics activity and develop a sustainable balance between economic, social, and environmental objectives [22]. Facility location modeling is a strategic planning design approach that selects the optimal set of facilities from a set of potential facilities [23]. As a branch of network design, the green logistics network design problem originates from the classic facility location problem (FLP) [7, 24] and then extends it to the location-allocation problem (LAP) [25–27].In the background of the vigorously developing green economy, Turken et al. [24] investigated the impact of carbon tax on plant capacity and location decisions of a firm. Yang et al. [26] studied the low-carbon network design problem of the third-party logistics distribution system. To solve location-allocation problem with a carbon emissions constraint, Rezaee et al. [2] developed a two-stage stochastic programming model, and Gao et al. [28] proposed five multidimensional mixed-integer nonlinear programming models. In practice, the urban agglomeration logistics network may comprise multiple urban logistics networks, which are independently and separately managed by local governments (or logistics authorities) with different targets. Nevertheless, the literature always supposed that the network is totally designed by a single management body, which may not suit the design of the logistics network in urban agglomeration.


## 1.1.3. Green Logistics Network Design under Route Choice Models

In general, route choice models should evaluate the utility of each route and route the flow of passengers or goods from origin to destination. Although they have been widely used to solve urban passenger transportation problems by researchers and engineers [29–33], using these models to solve urban cargo transportation problems has received little attention in the literature.The logistics network design problem based on route choice behavior can be characterized as a bilevel programming model or an equivalent mathematical program with equilibrium constraints (MPEC) model. At the upper level, the logistics authority determines the number, location, and capacity of logistics facilities (nodes). At the lower level, shippers and carriers follow the user equilibrium (UE) principle or stochastic user equilibrium (SUE) principle [34–36]. Yamada et al. [37] presented a bilevel programming model to simulate the multimodal freight transport network design problem, where the lower-level model is a multimodal multiclass user equilibrium model. To address the intermodal hub-and-spoke network design problem for multiple stakeholders and multitype containers, Meng and Wang [38] established an MPEC model and used a variational inequality to describe the operators’ route choice behavior. Wang and Meng [39] extended the research work of Meng and Wang [38] by considering congestion effects and piecewise linear cost functions and solved it with a solution algorithm based on nonlinear optimization and branch-and-bound global optimization.However, the route choice models of these studies did not consider the carbon emissions cost, perception error, and elastic demand. In addition, they mainly focus on static, determined network design issues. As strategic decision planning, the design of a logistics network in urban agglomeration involves long-term planning of a between 5- and 20-year timeline, where the logistics demand is often uncertain with the development of urban evolution and industrial structure optimization. Hence, it is more realistic and reasonable to design the logistics network with demand uncertainty and stochastic route choice.


## 1.2. Objectives and Contributions

In view of the abovementioned realistic problems and specific characteristics of logistics development in urban agglomerations, the objective of this paper is to coordinate the distribution of various types of logistics nodes (e.g., logistics parks, distribution centers and logistics terminals) in the urban agglomeration logistics network from the overall perspective of urban agglomeration and fully consider the low carbon requirement together with demand uncertainties. The main contributions of this paper are as follows. First, to characterize different decision-making behaviors among multiple local authorities, the multimodal logistics network design problem is formulated as two nonlinear bilevel programming models. At the upper level, the logistics authority of each city attempts to maximize the social welfare including producer surplus and consumer surplus. At the lower level, the logistics users’ route choice decisions follow logit-based stochastic user equilibrium (SUE) with elastic demand under logistics demand scenario. Second, a heuristic solution algorithm that is a combination of quantum behaved particle swarm optimization (QPSO) and Method of Successive Averages (MSA) is developed to solve the proposed bilevel programming model. Third, the optimal number, scale, and location of logistics nodes, analysis of the impact on the distribution of city group logistics network of infrastructure investment budget, and carbon tax are illustrated with an example.The remainder of this article is arranged as follows. The basic considerations of this paper, including general assumptions and network representation, are described in Section2. A nonlinear bilevel programming formulation is proposed in Section 3, and the solution methods are presented in Section 4. A numerical example is presented to illustrate the availability of previous models in Section 5. Finally, the conclusion and future studies are discussed in Section 6.


## 2. Basic Considerations

### 2.1. General Assumptions

To facilitate the presentation of essential ideas without the loss of generality, the following basic assumptions are made.A1. For the simplicity of expression, the urban agglomeration is assumed to be a city that is a set of all single cities, and the planning period is assumed to be one week.A2. In the urban logistics system, the logistics nodes investment and subsidy are determined by the logistics authority.A3. The disutility of each service route is measured by transport time, transport cost, and CO2 emission taxes (if any). Logistics users select their logistics service routes which are associated with their own perceptions of service disutility.A4. The urban agglomeration comprises several cities, and each city has one logistics authority. In decentralized decision-making, the decision-making of logistics departments in urban agglomeration is completely independent.


### 2.2. Network Representation

To model logistics service of the urban agglomeration, we first represent the demand network, multimodal logistics physical network, and multimodal logistics service network.Demand Network. Let O
∈
N
0 be the set of logistics demand origin nodes and let D
∈
N
0 be the set of logistics demand destination nodes, where N
0 is the set of logistics nodes (or transfer nodes) including existing and potential logistics parks, distribution centers, and logistics terminals. Denote by W
⊂
O
×
D the set of logistics demand origin-destination (O–D) pairs. For a given logistics O-D pair, Figure 1(a) shows the different types of logistics demands such as industrial demand (K1), commercial demand (K2), and agricultural demand (K3). These demands are served by the multimodal logistics network, as shown in Figures 1(b)–1(e).Figure 1
Multimodal logistics network representation.
(a)
Logistics user demand networkG
d
=
N
d
,
A
d
(b)
Multimodal logistics physical networkG
0
=
N
0
,
A
0
(c)
Representation of the virtual subnetwork of logistics transfer node 1
(d)
Multimodal logistics service supernetworkG
=
M
,
N
,
A
g
,
A
t
(e)
Logistics service routeMultimodal Logistics Physical Network. Suppose that a multimodal logistics physical network needs to be designed, which is made up of a set of logistics nodes and a set of logistics links or arcs of different transport modes. Each logistics node or link is provided by a logistics operator.We denote the multimodal logistics physical network as a directed networkG
0
=
N
0
,
A
0, where N
0 is the set of logistics nodes and A
0 is the set of logistics links. Let N
t
⊂
N
0 denote the set of logistics transfer nodes, and let M
0 be the set of modes. Each link a
∈
A
0 is represented by a triplet s
-
a
,
e
-
a
,
m
a, where s
-
a
,
e
-
a
∈
N
0 are the starting point and ending point of link a, respectively, and m
a
∈
M
0 is the transportation mode on the link.Analogously, for any logistics transfer nodei
∈
N
t, let A
i
n
i denote the set of logistics links pointing into transfer node i, and let A
o
u
t
i denote the set of logistics links stemming out of transfer node i. For each logistics link a
∈
A
i
n
i, we make a copy of transfer node i as its ending point. Let N
i
n
i denote the set of heads (copies) of all logistics arcs in A
i
n
i after copying. Similarly, For each logistics link a
∈
A
o
u
t
i, we make a copy of transfer node i as its ending point. Let N
o
u
t
i denote the set of tails (copies) of all logistics arcs in A
o
u
t
i after copying. Thus, the set of virtual arcs at transfer node i is denoted by A
t
i
:
=
j
,
k
,
j
∈
N
i
n
i
,
k
∈
N
o
u
t
i, and the virtual subnetwork of logistics transfer node i can be represented by a directed graph G
i
=
N
i
n
i
∪
N
o
u
t
i
,
A
t
i.Figure1(b) shows an example of a multimodal logistics physical network G
0
=
N
0
,
A
0, where N
0
=
O
,
D
,
1
,
2
,
3 and A
0
=
a
i
∣
i
=
1,2
,
…
,
6
. We suppose N
t
=
1
,
2
,
3 and M
=
R
a
i
l
w
a
y
,
E
x
p
r
e
s
s
w
a
y
. Then, the virtual subnetworks at each logistics transfer node i
∈
N
tcan be generated. Figure 1(c) shows the example of the virtual subnetwork G
1
=
N
t
1
,
A
t
1at logistics transfer node 1, where N
t
1
=
1
a
,
1
b
,
1
c and A
t
1
=
a
7
,
a
8.Multimodal Logistics Service Network. A hypernetwork G
=
M
0
,
N
,
A
g
,
A
t is used to construct a multimodal logistics service network, where M
0 represents the type of route transport modes, N
:
=
⋃
i
∈
N
t
N
i
n
i
∪
N
o
u
t
i
∪
N
i
∈
N
0
∖
N
t denotes the set of nodes after transfer node copies, A
g
:
=
A
p denotes the set of logistics links, and A
t
:
=
⋃
i
∈
N
t
A
t
i denotes the set of transfer arcs. Based on the example of multimodal logistics physical network in Figure 1(b) and the virtual transfer subnetworks, the logistics service supernetwork G
=
M
0
,
N
,
A
g
,
A
t is then generated and contains highway subnetworks and railway subnetworks (Figure 1(d)).Feasible Sets of Route. Although in the multimodal logistics service supernetwork there are some routes from the logistics demand origin to the destination, only the routes that satisfy certain conditions, such as cost, time, or CO2 emissions, are called feasible routes. Thus, for any logistics O-D pair w
∈
W, let R
w denote the set of logistics feasible service routes connecting w in the multimodal logistics service supernetwork, and let δ
w
r
,
m be one if service route r
∈
R
w uses transport mode m
∈
M and zero otherwise.


## 2.1. General Assumptions

To facilitate the presentation of essential ideas without the loss of generality, the following basic assumptions are made.A1. For the simplicity of expression, the urban agglomeration is assumed to be a city that is a set of all single cities, and the planning period is assumed to be one week.A2. In the urban logistics system, the logistics nodes investment and subsidy are determined by the logistics authority.A3. The disutility of each service route is measured by transport time, transport cost, and CO2 emission taxes (if any). Logistics users select their logistics service routes which are associated with their own perceptions of service disutility.A4. The urban agglomeration comprises several cities, and each city has one logistics authority. In decentralized decision-making, the decision-making of logistics departments in urban agglomeration is completely independent.


## 2.2. Network Representation

To model logistics service of the urban agglomeration, we first represent the demand network, multimodal logistics physical network, and multimodal logistics service network.Demand Network. Let O
∈
N
0 be the set of logistics demand origin nodes and let D
∈
N
0 be the set of logistics demand destination nodes, where N
0 is the set of logistics nodes (or transfer nodes) including existing and potential logistics parks, distribution centers, and logistics terminals. Denote by W
⊂
O
×
D the set of logistics demand origin-destination (O–D) pairs. For a given logistics O-D pair, Figure 1(a) shows the different types of logistics demands such as industrial demand (K1), commercial demand (K2), and agricultural demand (K3). These demands are served by the multimodal logistics network, as shown in Figures 1(b)–1(e).Figure 1
Multimodal logistics network representation.
(a)
Logistics user demand networkG
d
=
N
d
,
A
d
(b)
Multimodal logistics physical networkG
0
=
N
0
,
A
0
(c)
Representation of the virtual subnetwork of logistics transfer node 1
(d)
Multimodal logistics service supernetworkG
=
M
,
N
,
A
g
,
A
t
(e)
Logistics service routeMultimodal Logistics Physical Network. Suppose that a multimodal logistics physical network needs to be designed, which is made up of a set of logistics nodes and a set of logistics links or arcs of different transport modes. Each logistics node or link is provided by a logistics operator.We denote the multimodal logistics physical network as a directed networkG
0
=
N
0
,
A
0, where N
0 is the set of logistics nodes and A
0 is the set of logistics links. Let N
t
⊂
N
0 denote the set of logistics transfer nodes, and let M
0 be the set of modes. Each link a
∈
A
0 is represented by a triplet s
-
a
,
e
-
a
,
m
a, where s
-
a
,
e
-
a
∈
N
0 are the starting point and ending point of link a, respectively, and m
a
∈
M
0 is the transportation mode on the link.Analogously, for any logistics transfer nodei
∈
N
t, let A
i
n
i denote the set of logistics links pointing into transfer node i, and let A
o
u
t
i denote the set of logistics links stemming out of transfer node i. For each logistics link a
∈
A
i
n
i, we make a copy of transfer node i as its ending point. Let N
i
n
i denote the set of heads (copies) of all logistics arcs in A
i
n
i after copying. Similarly, For each logistics link a
∈
A
o
u
t
i, we make a copy of transfer node i as its ending point. Let N
o
u
t
i denote the set of tails (copies) of all logistics arcs in A
o
u
t
i after copying. Thus, the set of virtual arcs at transfer node i is denoted by A
t
i
:
=
j
,
k
,
j
∈
N
i
n
i
,
k
∈
N
o
u
t
i, and the virtual subnetwork of logistics transfer node i can be represented by a directed graph G
i
=
N
i
n
i
∪
N
o
u
t
i
,
A
t
i.Figure1(b) shows an example of a multimodal logistics physical network G
0
=
N
0
,
A
0, where N
0
=
O
,
D
,
1
,
2
,
3 and A
0
=
a
i
∣
i
=
1,2
,
…
,
6
. We suppose N
t
=
1
,
2
,
3 and M
=
R
a
i
l
w
a
y
,
E
x
p
r
e
s
s
w
a
y
. Then, the virtual subnetworks at each logistics transfer node i
∈
N
tcan be generated. Figure 1(c) shows the example of the virtual subnetwork G
1
=
N
t
1
,
A
t
1at logistics transfer node 1, where N
t
1
=
1
a
,
1
b
,
1
c and A
t
1
=
a
7
,
a
8.Multimodal Logistics Service Network. A hypernetwork G
=
M
0
,
N
,
A
g
,
A
t is used to construct a multimodal logistics service network, where M
0 represents the type of route transport modes, N
:
=
⋃
i
∈
N
t
N
i
n
i
∪
N
o
u
t
i
∪
N
i
∈
N
0
∖
N
t denotes the set of nodes after transfer node copies, A
g
:
=
A
p denotes the set of logistics links, and A
t
:
=
⋃
i
∈
N
t
A
t
i denotes the set of transfer arcs. Based on the example of multimodal logistics physical network in Figure 1(b) and the virtual transfer subnetworks, the logistics service supernetwork G
=
M
0
,
N
,
A
g
,
A
t is then generated and contains highway subnetworks and railway subnetworks (Figure 1(d)).Feasible Sets of Route. Although in the multimodal logistics service supernetwork there are some routes from the logistics demand origin to the destination, only the routes that satisfy certain conditions, such as cost, time, or CO2 emissions, are called feasible routes. Thus, for any logistics O-D pair w
∈
W, let R
w denote the set of logistics feasible service routes connecting w in the multimodal logistics service supernetwork, and let δ
w
r
,
m be one if service route r
∈
R
w uses transport mode m
∈
M and zero otherwise.


## 3. Model Formulation

### 3.1. The General Stochastic Bilevel Programming Model

Stochastic bilevel programming combines the characteristics of stochastic programming and bilevel programming and introduces random scenarios to describe the uncertainties involved in the model. The most commonly used stochastic bilevel programming model is the expected bilevel programming model, which is described in the following mathematical form [40]:(1)
U
0
m
i
n
x
E
F
x
,
y
ξ
,
ξ
s
.
t
.
E
G
x
,
y
ξ
,
ξ
≤
0where response function y
=
y
x is implicitly defined by(2)
L
0
m
i
n
y
f
x
,
y
ξ
,
ξ
s
.
t
.
g
x
,
y
ξ
,
ξ
≤
0Obviously, the expected bilevel programming model comprises two submodels, (U0), which is defined as an upper level and (L0), which is a lower level.x, F
x
,
y
ξ
,
ξ and G
x
,
y
ξ
,
ξdenote the decision vector, objective function, and constraint set of the upper-level decision-makers or system managers, respectively. E
· is the expectation operator with respect to random scenarios ξ. y
ξ, f
x
,
y
ξ
,
ξ, and g
x
,
y
ξ
,
ξdenote the decision vector, objective function, and constraint set of the lower-level decision-makers or users, respectively.Assume that the random scenariosξ have a finite number of states,Ψ
=
σ
s
,
s
∈
Ω. Let s
=
1,2
,
…
,
Ω index be its possible realizations and let p
s be their respective probabilities. Then, we can now express the extensive form of the stochastic bilevel program as follows:(3)
U
1
m
i
n
x
∑
s
∈
Ω
p
s
F
x
,
y
s
σ
s
,
σ
s
s
.
t
.
∑
s
∈
Ω
p
s
G
x
,
y
s
σ
s
,
σ
s
≤
0where response function y
=
y
x is implicitly defined by(4)
L
1
m
i
n
y
s
f
x
,
y
s
σ
s
,
σ
s
s
.
t
.
g
x
,
y
s
σ
s
,
σ
s
≤
0


### 3.2. Decentralized Decision in Multiple Cities

Since the urban logistics system is exclusively managed by a logistics authority but serves for all logistics users (i.e., shippers and carriers) equally, the urban logistics network design is typically formulated as a Stackelberg game that can well characterize the interactions between the logistics authority and users.As shown in Figure2, a bilevel program is used to model the leader-follower behaviors between logistics authority and users. At the upper level, the logistics authority attempts to maximize the social welfare by planning the investment capacity of logistics nodes and determining the carbon tax on transportation service links. At the lower level, the logistics users’ reactions and choice decisions to the urban logistics network design scheme will be assumed to follow logit-based stochastic user equilibrium (SUE) under the logistics demand scenario. Moreover, the main notations are provided in Appendix A.Figure 2
Decision framework of green logistics system.


#### 3.2.1. Upper-Level Model of Green Logistics Network Design for a Single City

As the upper-level decision maker, the logistics authority aims to maximize the total social welfare of its own city. It is well known that the total social welfare comprises the consumer surplus and producer surplus in scenarios
∈
Ω, which can be expressed as(5)
S
W
k
s
X
k
,
y
k
,
V
k
=
C
S
k
s
+
P
S
k
,
1
s
+
P
S
k
,
2
swhere(6)
C
S
k
s
=
∑
w
∈
W
k
∫
0
q
k
,
w
s
D
w
-
1
w
d
w
-
∑
w
∈
W
k
q
k
,
w
s
λ
k
,
w
s
,
∀
k
∈
K
,
s
∈
Ω
(7)
P
S
k
,
1
s
=
∑
i
∈
I
k
C
k
,
i
-
η
k
,
i
f
k
,
i
s
+
∑
i
∈
P
k
m
a
x
x
k
,
i
-
z
k
,
i
,
0
S
u
b
k
,
i
-
∑
i
∈
I
k
C
k
,
i
0
x
k
,
i
ρ
,
∀
k
∈
K
,
s
∈
Ω
(8)
P
S
k
,
2
s
=
∑
m
∈
M
0
∑
a
∈
A
k
g
v
k
,
a
m
,
s
l
k
,
a
c
k
,
a
m
-
τ
k
,
a
m
+
∑
m
∈
M
0
∑
a
∈
A
k
g
v
k
,
a
m
,
s
l
k
,
a
e
a
m
y
k
,
∀
k
∈
K
,
s
∈
ΩEq. (6) formulates the consumer surplus, and D
w
-
1
w is the inverse function of the logistics demand function. Consumer surplus is the extra benefit logistics users’ gains when the costs they actually pay are less than what they would be prepared to pay. Eqs. (7) and (8) formulate the producer surplus of all logistics nodes and logistics arcs, respectively. Producer surplus is a measure of producer welfare. Taking (7) as an example, the first term is transfer profits, the second term is subsidy revenues, and the last is park construction costs. Then, the upper-level model of green logistics network design in single city k
∈
K is formulated as(9)
m
a
x
X
k
,
y
k
S
W
X
k
,
y
k
,
V
k
=
∑
s
∈
Ω
p
k
s
C
S
k
s
+
P
S
k
,
1
s
+
P
S
k
,
2
ssubject to(10)
0
≤
x
k
,
i
≤
S
k
,
i
m
a
x
,
∀
i
∈
P
k
(11)
0
≤
y
k
≤
C
h
k
m
a
x
(12)
∑
i
∈
I
k
C
k
,
i
0
g
x
k
,
i
ρ
≤
B
kwhere V
k
X
k
,
y
k can be obtained by solving the upper-level model.The objective function (9) formulates the expected social welfares. Constraint (10) represents the establishment or expansion capacity restraint of the logistics nodes; constraint (11) denotes the constraint of carbon tax; and constraint (12) represents the construction investment constraint of the logistics nodes.


#### 3.2.2. Lower-Level Model of Green Logistics Network Design

Travel/Transfer Time. To capture the difference in attributes of different modes of transport, for each transport link a
∈
A
k
g, we consider the following service time function:(13)
t
k
,
a
m
,
s
v
k
,
a
m
,
s
=
t
k
,
a
m
,
0
1
+
0.15
v
k
,
a
m
,
s
C
a
p
k
,
a
m
4
,
m
=
1,2
t
k
,
a
m
,
0
+
t
k
,
m
d
m
a
x
v
k
,
a
m
,
s
-
C
a
p
k
,
a
m
,
0
C
a
p
k
,
a
m
,
m
=
3,4
∀
k
∈
K
,
a
∈
A
k
g
,
s
∈
Ωwhere t
a
m
,
0 is the link free-flow transport service time, and t
k
,
m
d is the average transport time interval. For HGVs or LGVs, the Bureau of Public Roads-type (BPR) function can be adopted to estimate transport service time. For railways or waterways, we consider service time function as a function of link free-flow transport service time and departure interval time [11, 35]. Similarly, for each virtual transfer arc, we use the following service time function:(14)
t
k
,
a
i
,
s
f
k
,
i
s
=
t
k
,
a
i
,
0
1
+
α
0
f
k
,
i
s
x
k
,
i
β
0
,
x
k
,
i
≠
0
+
∞
,
,
x
k
,
i
=
0
∀
k
∈
K
,
i
∈
I
k
,
a
∈
A
k
t
,
s
∈
Ωwhere t
k
,
a
i
,
0 is the arc free-flow transfer service time and α
0 and β
0 are impedance parameters.Route Utility and Flow. According to A3, each logistics service route is associated with a given actual cost (disutility), which can be expressed as(15)
u
k
,
w
r
,
s
=
C
k
,
w
r
,
s
+
τ
v
o
t
T
k
,
w
r
,
s
+
G
k
,
w
r
,
s
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
s
∈
Ωwhere τ
v
o
t is the value of time. C
k
,
w
r
,
s, T
k
,
w
r
,
s, and G
k
,
w
r
,
s represent the transportation cost, logistics service time, and CO2 emission cost on service route r between O-D pair w in city k
∈
K
- under logistics demand scenarios s, respectively, which are expressed as(16)
C
k
,
w
r
,
s
=
∑
m
∈
M
0
∑
a
∈
A
k
g
c
k
,
a
m
l
k
,
a
δ
a
,
r
m
+
∑
i
∈
I
k
∑
a
∈
A
k
t
c
k
,
a
i
δ
i
,
r
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
s
∈
Ω
(17)
T
k
,
w
r
,
s
=
∑
m
∈
M
0
∑
a
∈
A
k
g
t
k
,
a
m
,
s
v
k
,
a
m
,
s
δ
a
,
r
m
+
∑
i
∈
I
k
∑
a
∈
A
k
t
t
k
,
a
i
,
s
f
k
,
i
s
δ
i
,
r
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
m
∈
M
,
s
∈
Ω
(18)
G
k
,
w
r
,
s
=
∑
m
∈
M
0
∑
a
∈
A
k
g
e
m
l
k
,
a
y
k
δ
a
,
r
m
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
m
∈
M
,
s
∈
ΩIt is worth noting that the transportation cost (time) of each route includes the transport cost (time) on links and the transfer cost (time) at parks. Due to variations in perception, the route service disutility is perceived differently by each logistics user, and thus, the perceived disutility of each route is treated as a random variable. If the random variable can be considered to obeyGumbel distribution[41], then the path flow f
k
,
w
r
,
s on route r
∈
R
w between O-D pair w
∈
W
k in single city k
∈
K and scenario s
∈
Ω can be given by(19)
f
k
,
w
r
,
s
=
q
k
,
w
s
e
x
p
-
θ
u
k
,
w
r
,
s
∑
r
∈
R
w
e
x
p
-
θ
u
k
,
w
r
,
s
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
s
∈
Ωwhere θ represents the sensitivity of route selection disutility and q
k
,
w
s denotes the logistics demand function. Eq. (19) is the most widely used flow assignment method in traffic planning and represents the flow assigned to each feasible route for each logistic demand [41]. To capture the logistics users’ responses to logistics service disutility, we assume that the elastic demand function between a generic O-D pair w
∈
W
k is a monotonically decreasing function of the O-D service cost u
k
,
w
r
,
s between this O-D pair. The generic elastic demand function [42] is expressed as follows:(20)
q
k
,
w
s
=
D
w
λ
k
,
w
s
=
q
-
k
,
w
s
exp
⁡
-
β
λ
k
,
w
s
,
∀
k
∈
K
,
w
∈
W
k
,
s
∈
Ωwhere β represents the sensitivity to the expected service cost and λ
k
,
w
s is the expected minimum perceived service cost between O-D pair w
∈
W
k in single city k
∈
K and scenario s
∈
Ω. In the case of stochastic user equilibrium (SUE) assignment with elastic demand [42], the expected minimum service cost between an O-D pair under logistics demand scenario s
∈
Ω in city k
∈
K could be expressed as(21)
λ
k
,
w
s
=
-
1
θ
ln
⁡
∑
r
∈
R
w
exp
⁡
-
θ
u
k
,
w
r
,
sEquivalence Model. According to the above analysis, for any given logistics authority decision X
k
,
y
k, the lower-level model of green logistics network design in a single city k can be formulated as the following equivalent minimization program:(22)
m
i
n
Z
F
=
1
θ
∑
w
∈
W
k
∑
r
∈
R
w
f
k
,
w
r
,
s
ln
⁡
f
k
,
w
r
,
s
-
1
-
1
θ
∑
w
∈
W
k
q
k
,
w
s
ln
⁡
q
k
,
w
s
-
1
-
∑
w
∈
W
k
∫
0
q
k
,
w
s
D
w
-
1
w
d
w
+
∑
m
∈
M
0
∑
a
∈
A
k
g
τ
v
o
t
t
k
,
a
m
,
s
v
k
,
a
s
+
c
k
,
a
m
l
k
,
a
+
e
m
l
k
,
a
y
k
v
k
,
a
s
+
∑
i
∈
I
k
∑
a
∈
A
k
t
τ
v
o
t
t
k
,
a
i
f
k
,
i
s
+
c
k
,
a
i
f
k
,
i
ssubject to(23)
∑
r
∈
R
w
f
k
,
w
r
,
s
=
q
k
,
w
s
,
∀
k
∈
K
,
w
∈
W
k
,
s
∈
Ω
(24)
f
k
,
w
r
,
s
>
0
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
s
∈
Ωwhere(25)
v
k
,
a
s
=
∑
w
∈
W
∑
r
∈
R
w
f
k
,
w
r
,
s
δ
a
,
r
m
,
∀
k
∈
K
,
a
∈
A
k
g
,
s
∈
Ω
(26)
f
k
,
i
s
=
∑
w
∈
W
∑
r
∈
R
w
f
k
,
w
r
,
s
δ
i
,
r
,
∀
k
∈
K
,
i
∈
I
k
,
s
∈
ΩConstraint (23) is the flow conservation constraint. Constraint (24) is the nonnegativity constraints of service route flows. Constraint (25) defines the relationship between the route flow and path flow. Constraint (26) defines the relationship between the route flow and logistics nodes flow, respectively.Proposition 1.
The minimization program (22)-(26) is equivalent to the logit-based stochastic user equilibrium (SUE) assignment with elastic demand under logistics demand scenario s
∈
Ω in city k
∈
K. The proof of Proposition 1 is given in Appendix B.


### 3.3. Centralized Decision in Urban Agglomeration

In the previous subsection, we study the problem that several logistics authorities perform logistics network design decisions independently. However, such independence might be adverse to the overall benefit of urban agglomeration, since it does not explore the full potential of the investment to improve the need of all logistics users on the entire logistics network. We propose a centralized decision model in which logistics authorities cooperate with each other to maximize the total social welfare in urban agglomeration.The upper-level model of green logistics network design in urban agglomeration is formulated as(27)
m
a
x
X
K
+
1
,
y
K
+
1
S
W
X
K
+
1
,
y
K
+
1
,
V
K
+
1
=
∑
s
∈
Ω
p
s
C
S
K
+
1
s
+
P
S
K
+
1,1
s
+
P
S
K
+
1,2
ssubject to(28)
0
≤
x
i
≤
S
i
m
a
x
,
∀
i
∈
P
K
+
1
(29)
0
≤
y
K
+
1
≤
C
h
K
+
1
m
a
x
(30)
∑
i
∈
I
K
+
1
C
i
0
g
x
i
ρ
≤
B
K
+
1The lower level program is still the traffic assignment problem described in Section3.2.2.


## 3.1. The General Stochastic Bilevel Programming Model

Stochastic bilevel programming combines the characteristics of stochastic programming and bilevel programming and introduces random scenarios to describe the uncertainties involved in the model. The most commonly used stochastic bilevel programming model is the expected bilevel programming model, which is described in the following mathematical form [40]:(1)
U
0
m
i
n
x
E
F
x
,
y
ξ
,
ξ
s
.
t
.
E
G
x
,
y
ξ
,
ξ
≤
0where response function y
=
y
x is implicitly defined by(2)
L
0
m
i
n
y
f
x
,
y
ξ
,
ξ
s
.
t
.
g
x
,
y
ξ
,
ξ
≤
0Obviously, the expected bilevel programming model comprises two submodels, (U0), which is defined as an upper level and (L0), which is a lower level.x, F
x
,
y
ξ
,
ξ and G
x
,
y
ξ
,
ξdenote the decision vector, objective function, and constraint set of the upper-level decision-makers or system managers, respectively. E
· is the expectation operator with respect to random scenarios ξ. y
ξ, f
x
,
y
ξ
,
ξ, and g
x
,
y
ξ
,
ξdenote the decision vector, objective function, and constraint set of the lower-level decision-makers or users, respectively.Assume that the random scenariosξ have a finite number of states,Ψ
=
σ
s
,
s
∈
Ω. Let s
=
1,2
,
…
,
Ω index be its possible realizations and let p
s be their respective probabilities. Then, we can now express the extensive form of the stochastic bilevel program as follows:(3)
U
1
m
i
n
x
∑
s
∈
Ω
p
s
F
x
,
y
s
σ
s
,
σ
s
s
.
t
.
∑
s
∈
Ω
p
s
G
x
,
y
s
σ
s
,
σ
s
≤
0where response function y
=
y
x is implicitly defined by(4)
L
1
m
i
n
y
s
f
x
,
y
s
σ
s
,
σ
s
s
.
t
.
g
x
,
y
s
σ
s
,
σ
s
≤
0


## 3.2. Decentralized Decision in Multiple Cities

Since the urban logistics system is exclusively managed by a logistics authority but serves for all logistics users (i.e., shippers and carriers) equally, the urban logistics network design is typically formulated as a Stackelberg game that can well characterize the interactions between the logistics authority and users.As shown in Figure2, a bilevel program is used to model the leader-follower behaviors between logistics authority and users. At the upper level, the logistics authority attempts to maximize the social welfare by planning the investment capacity of logistics nodes and determining the carbon tax on transportation service links. At the lower level, the logistics users’ reactions and choice decisions to the urban logistics network design scheme will be assumed to follow logit-based stochastic user equilibrium (SUE) under the logistics demand scenario. Moreover, the main notations are provided in Appendix A.Figure 2
Decision framework of green logistics system.


### 3.2.1. Upper-Level Model of Green Logistics Network Design for a Single City

As the upper-level decision maker, the logistics authority aims to maximize the total social welfare of its own city. It is well known that the total social welfare comprises the consumer surplus and producer surplus in scenarios
∈
Ω, which can be expressed as(5)
S
W
k
s
X
k
,
y
k
,
V
k
=
C
S
k
s
+
P
S
k
,
1
s
+
P
S
k
,
2
swhere(6)
C
S
k
s
=
∑
w
∈
W
k
∫
0
q
k
,
w
s
D
w
-
1
w
d
w
-
∑
w
∈
W
k
q
k
,
w
s
λ
k
,
w
s
,
∀
k
∈
K
,
s
∈
Ω
(7)
P
S
k
,
1
s
=
∑
i
∈
I
k
C
k
,
i
-
η
k
,
i
f
k
,
i
s
+
∑
i
∈
P
k
m
a
x
x
k
,
i
-
z
k
,
i
,
0
S
u
b
k
,
i
-
∑
i
∈
I
k
C
k
,
i
0
x
k
,
i
ρ
,
∀
k
∈
K
,
s
∈
Ω
(8)
P
S
k
,
2
s
=
∑
m
∈
M
0
∑
a
∈
A
k
g
v
k
,
a
m
,
s
l
k
,
a
c
k
,
a
m
-
τ
k
,
a
m
+
∑
m
∈
M
0
∑
a
∈
A
k
g
v
k
,
a
m
,
s
l
k
,
a
e
a
m
y
k
,
∀
k
∈
K
,
s
∈
ΩEq. (6) formulates the consumer surplus, and D
w
-
1
w is the inverse function of the logistics demand function. Consumer surplus is the extra benefit logistics users’ gains when the costs they actually pay are less than what they would be prepared to pay. Eqs. (7) and (8) formulate the producer surplus of all logistics nodes and logistics arcs, respectively. Producer surplus is a measure of producer welfare. Taking (7) as an example, the first term is transfer profits, the second term is subsidy revenues, and the last is park construction costs. Then, the upper-level model of green logistics network design in single city k
∈
K is formulated as(9)
m
a
x
X
k
,
y
k
S
W
X
k
,
y
k
,
V
k
=
∑
s
∈
Ω
p
k
s
C
S
k
s
+
P
S
k
,
1
s
+
P
S
k
,
2
ssubject to(10)
0
≤
x
k
,
i
≤
S
k
,
i
m
a
x
,
∀
i
∈
P
k
(11)
0
≤
y
k
≤
C
h
k
m
a
x
(12)
∑
i
∈
I
k
C
k
,
i
0
g
x
k
,
i
ρ
≤
B
kwhere V
k
X
k
,
y
k can be obtained by solving the upper-level model.The objective function (9) formulates the expected social welfares. Constraint (10) represents the establishment or expansion capacity restraint of the logistics nodes; constraint (11) denotes the constraint of carbon tax; and constraint (12) represents the construction investment constraint of the logistics nodes.


### 3.2.2. Lower-Level Model of Green Logistics Network Design

Travel/Transfer Time. To capture the difference in attributes of different modes of transport, for each transport link a
∈
A
k
g, we consider the following service time function:(13)
t
k
,
a
m
,
s
v
k
,
a
m
,
s
=
t
k
,
a
m
,
0
1
+
0.15
v
k
,
a
m
,
s
C
a
p
k
,
a
m
4
,
m
=
1,2
t
k
,
a
m
,
0
+
t
k
,
m
d
m
a
x
v
k
,
a
m
,
s
-
C
a
p
k
,
a
m
,
0
C
a
p
k
,
a
m
,
m
=
3,4
∀
k
∈
K
,
a
∈
A
k
g
,
s
∈
Ωwhere t
a
m
,
0 is the link free-flow transport service time, and t
k
,
m
d is the average transport time interval. For HGVs or LGVs, the Bureau of Public Roads-type (BPR) function can be adopted to estimate transport service time. For railways or waterways, we consider service time function as a function of link free-flow transport service time and departure interval time [11, 35]. Similarly, for each virtual transfer arc, we use the following service time function:(14)
t
k
,
a
i
,
s
f
k
,
i
s
=
t
k
,
a
i
,
0
1
+
α
0
f
k
,
i
s
x
k
,
i
β
0
,
x
k
,
i
≠
0
+
∞
,
,
x
k
,
i
=
0
∀
k
∈
K
,
i
∈
I
k
,
a
∈
A
k
t
,
s
∈
Ωwhere t
k
,
a
i
,
0 is the arc free-flow transfer service time and α
0 and β
0 are impedance parameters.Route Utility and Flow. According to A3, each logistics service route is associated with a given actual cost (disutility), which can be expressed as(15)
u
k
,
w
r
,
s
=
C
k
,
w
r
,
s
+
τ
v
o
t
T
k
,
w
r
,
s
+
G
k
,
w
r
,
s
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
s
∈
Ωwhere τ
v
o
t is the value of time. C
k
,
w
r
,
s, T
k
,
w
r
,
s, and G
k
,
w
r
,
s represent the transportation cost, logistics service time, and CO2 emission cost on service route r between O-D pair w in city k
∈
K
- under logistics demand scenarios s, respectively, which are expressed as(16)
C
k
,
w
r
,
s
=
∑
m
∈
M
0
∑
a
∈
A
k
g
c
k
,
a
m
l
k
,
a
δ
a
,
r
m
+
∑
i
∈
I
k
∑
a
∈
A
k
t
c
k
,
a
i
δ
i
,
r
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
s
∈
Ω
(17)
T
k
,
w
r
,
s
=
∑
m
∈
M
0
∑
a
∈
A
k
g
t
k
,
a
m
,
s
v
k
,
a
m
,
s
δ
a
,
r
m
+
∑
i
∈
I
k
∑
a
∈
A
k
t
t
k
,
a
i
,
s
f
k
,
i
s
δ
i
,
r
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
m
∈
M
,
s
∈
Ω
(18)
G
k
,
w
r
,
s
=
∑
m
∈
M
0
∑
a
∈
A
k
g
e
m
l
k
,
a
y
k
δ
a
,
r
m
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
m
∈
M
,
s
∈
ΩIt is worth noting that the transportation cost (time) of each route includes the transport cost (time) on links and the transfer cost (time) at parks. Due to variations in perception, the route service disutility is perceived differently by each logistics user, and thus, the perceived disutility of each route is treated as a random variable. If the random variable can be considered to obeyGumbel distribution[41], then the path flow f
k
,
w
r
,
s on route r
∈
R
w between O-D pair w
∈
W
k in single city k
∈
K and scenario s
∈
Ω can be given by(19)
f
k
,
w
r
,
s
=
q
k
,
w
s
e
x
p
-
θ
u
k
,
w
r
,
s
∑
r
∈
R
w
e
x
p
-
θ
u
k
,
w
r
,
s
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
s
∈
Ωwhere θ represents the sensitivity of route selection disutility and q
k
,
w
s denotes the logistics demand function. Eq. (19) is the most widely used flow assignment method in traffic planning and represents the flow assigned to each feasible route for each logistic demand [41]. To capture the logistics users’ responses to logistics service disutility, we assume that the elastic demand function between a generic O-D pair w
∈
W
k is a monotonically decreasing function of the O-D service cost u
k
,
w
r
,
s between this O-D pair. The generic elastic demand function [42] is expressed as follows:(20)
q
k
,
w
s
=
D
w
λ
k
,
w
s
=
q
-
k
,
w
s
exp
⁡
-
β
λ
k
,
w
s
,
∀
k
∈
K
,
w
∈
W
k
,
s
∈
Ωwhere β represents the sensitivity to the expected service cost and λ
k
,
w
s is the expected minimum perceived service cost between O-D pair w
∈
W
k in single city k
∈
K and scenario s
∈
Ω. In the case of stochastic user equilibrium (SUE) assignment with elastic demand [42], the expected minimum service cost between an O-D pair under logistics demand scenario s
∈
Ω in city k
∈
K could be expressed as(21)
λ
k
,
w
s
=
-
1
θ
ln
⁡
∑
r
∈
R
w
exp
⁡
-
θ
u
k
,
w
r
,
sEquivalence Model. According to the above analysis, for any given logistics authority decision X
k
,
y
k, the lower-level model of green logistics network design in a single city k can be formulated as the following equivalent minimization program:(22)
m
i
n
Z
F
=
1
θ
∑
w
∈
W
k
∑
r
∈
R
w
f
k
,
w
r
,
s
ln
⁡
f
k
,
w
r
,
s
-
1
-
1
θ
∑
w
∈
W
k
q
k
,
w
s
ln
⁡
q
k
,
w
s
-
1
-
∑
w
∈
W
k
∫
0
q
k
,
w
s
D
w
-
1
w
d
w
+
∑
m
∈
M
0
∑
a
∈
A
k
g
τ
v
o
t
t
k
,
a
m
,
s
v
k
,
a
s
+
c
k
,
a
m
l
k
,
a
+
e
m
l
k
,
a
y
k
v
k
,
a
s
+
∑
i
∈
I
k
∑
a
∈
A
k
t
τ
v
o
t
t
k
,
a
i
f
k
,
i
s
+
c
k
,
a
i
f
k
,
i
ssubject to(23)
∑
r
∈
R
w
f
k
,
w
r
,
s
=
q
k
,
w
s
,
∀
k
∈
K
,
w
∈
W
k
,
s
∈
Ω
(24)
f
k
,
w
r
,
s
>
0
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
s
∈
Ωwhere(25)
v
k
,
a
s
=
∑
w
∈
W
∑
r
∈
R
w
f
k
,
w
r
,
s
δ
a
,
r
m
,
∀
k
∈
K
,
a
∈
A
k
g
,
s
∈
Ω
(26)
f
k
,
i
s
=
∑
w
∈
W
∑
r
∈
R
w
f
k
,
w
r
,
s
δ
i
,
r
,
∀
k
∈
K
,
i
∈
I
k
,
s
∈
ΩConstraint (23) is the flow conservation constraint. Constraint (24) is the nonnegativity constraints of service route flows. Constraint (25) defines the relationship between the route flow and path flow. Constraint (26) defines the relationship between the route flow and logistics nodes flow, respectively.Proposition 1.
The minimization program (22)-(26) is equivalent to the logit-based stochastic user equilibrium (SUE) assignment with elastic demand under logistics demand scenario s
∈
Ω in city k
∈
K. The proof of Proposition 1 is given in Appendix B.


## 3.2.1. Upper-Level Model of Green Logistics Network Design for a Single City

As the upper-level decision maker, the logistics authority aims to maximize the total social welfare of its own city. It is well known that the total social welfare comprises the consumer surplus and producer surplus in scenarios
∈
Ω, which can be expressed as(5)
S
W
k
s
X
k
,
y
k
,
V
k
=
C
S
k
s
+
P
S
k
,
1
s
+
P
S
k
,
2
swhere(6)
C
S
k
s
=
∑
w
∈
W
k
∫
0
q
k
,
w
s
D
w
-
1
w
d
w
-
∑
w
∈
W
k
q
k
,
w
s
λ
k
,
w
s
,
∀
k
∈
K
,
s
∈
Ω
(7)
P
S
k
,
1
s
=
∑
i
∈
I
k
C
k
,
i
-
η
k
,
i
f
k
,
i
s
+
∑
i
∈
P
k
m
a
x
x
k
,
i
-
z
k
,
i
,
0
S
u
b
k
,
i
-
∑
i
∈
I
k
C
k
,
i
0
x
k
,
i
ρ
,
∀
k
∈
K
,
s
∈
Ω
(8)
P
S
k
,
2
s
=
∑
m
∈
M
0
∑
a
∈
A
k
g
v
k
,
a
m
,
s
l
k
,
a
c
k
,
a
m
-
τ
k
,
a
m
+
∑
m
∈
M
0
∑
a
∈
A
k
g
v
k
,
a
m
,
s
l
k
,
a
e
a
m
y
k
,
∀
k
∈
K
,
s
∈
ΩEq. (6) formulates the consumer surplus, and D
w
-
1
w is the inverse function of the logistics demand function. Consumer surplus is the extra benefit logistics users’ gains when the costs they actually pay are less than what they would be prepared to pay. Eqs. (7) and (8) formulate the producer surplus of all logistics nodes and logistics arcs, respectively. Producer surplus is a measure of producer welfare. Taking (7) as an example, the first term is transfer profits, the second term is subsidy revenues, and the last is park construction costs. Then, the upper-level model of green logistics network design in single city k
∈
K is formulated as(9)
m
a
x
X
k
,
y
k
S
W
X
k
,
y
k
,
V
k
=
∑
s
∈
Ω
p
k
s
C
S
k
s
+
P
S
k
,
1
s
+
P
S
k
,
2
ssubject to(10)
0
≤
x
k
,
i
≤
S
k
,
i
m
a
x
,
∀
i
∈
P
k
(11)
0
≤
y
k
≤
C
h
k
m
a
x
(12)
∑
i
∈
I
k
C
k
,
i
0
g
x
k
,
i
ρ
≤
B
kwhere V
k
X
k
,
y
k can be obtained by solving the upper-level model.The objective function (9) formulates the expected social welfares. Constraint (10) represents the establishment or expansion capacity restraint of the logistics nodes; constraint (11) denotes the constraint of carbon tax; and constraint (12) represents the construction investment constraint of the logistics nodes.


## 3.2.2. Lower-Level Model of Green Logistics Network Design

Travel/Transfer Time. To capture the difference in attributes of different modes of transport, for each transport link a
∈
A
k
g, we consider the following service time function:(13)
t
k
,
a
m
,
s
v
k
,
a
m
,
s
=
t
k
,
a
m
,
0
1
+
0.15
v
k
,
a
m
,
s
C
a
p
k
,
a
m
4
,
m
=
1,2
t
k
,
a
m
,
0
+
t
k
,
m
d
m
a
x
v
k
,
a
m
,
s
-
C
a
p
k
,
a
m
,
0
C
a
p
k
,
a
m
,
m
=
3,4
∀
k
∈
K
,
a
∈
A
k
g
,
s
∈
Ωwhere t
a
m
,
0 is the link free-flow transport service time, and t
k
,
m
d is the average transport time interval. For HGVs or LGVs, the Bureau of Public Roads-type (BPR) function can be adopted to estimate transport service time. For railways or waterways, we consider service time function as a function of link free-flow transport service time and departure interval time [11, 35]. Similarly, for each virtual transfer arc, we use the following service time function:(14)
t
k
,
a
i
,
s
f
k
,
i
s
=
t
k
,
a
i
,
0
1
+
α
0
f
k
,
i
s
x
k
,
i
β
0
,
x
k
,
i
≠
0
+
∞
,
,
x
k
,
i
=
0
∀
k
∈
K
,
i
∈
I
k
,
a
∈
A
k
t
,
s
∈
Ωwhere t
k
,
a
i
,
0 is the arc free-flow transfer service time and α
0 and β
0 are impedance parameters.Route Utility and Flow. According to A3, each logistics service route is associated with a given actual cost (disutility), which can be expressed as(15)
u
k
,
w
r
,
s
=
C
k
,
w
r
,
s
+
τ
v
o
t
T
k
,
w
r
,
s
+
G
k
,
w
r
,
s
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
s
∈
Ωwhere τ
v
o
t is the value of time. C
k
,
w
r
,
s, T
k
,
w
r
,
s, and G
k
,
w
r
,
s represent the transportation cost, logistics service time, and CO2 emission cost on service route r between O-D pair w in city k
∈
K
- under logistics demand scenarios s, respectively, which are expressed as(16)
C
k
,
w
r
,
s
=
∑
m
∈
M
0
∑
a
∈
A
k
g
c
k
,
a
m
l
k
,
a
δ
a
,
r
m
+
∑
i
∈
I
k
∑
a
∈
A
k
t
c
k
,
a
i
δ
i
,
r
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
s
∈
Ω
(17)
T
k
,
w
r
,
s
=
∑
m
∈
M
0
∑
a
∈
A
k
g
t
k
,
a
m
,
s
v
k
,
a
m
,
s
δ
a
,
r
m
+
∑
i
∈
I
k
∑
a
∈
A
k
t
t
k
,
a
i
,
s
f
k
,
i
s
δ
i
,
r
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
m
∈
M
,
s
∈
Ω
(18)
G
k
,
w
r
,
s
=
∑
m
∈
M
0
∑
a
∈
A
k
g
e
m
l
k
,
a
y
k
δ
a
,
r
m
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
m
∈
M
,
s
∈
ΩIt is worth noting that the transportation cost (time) of each route includes the transport cost (time) on links and the transfer cost (time) at parks. Due to variations in perception, the route service disutility is perceived differently by each logistics user, and thus, the perceived disutility of each route is treated as a random variable. If the random variable can be considered to obeyGumbel distribution[41], then the path flow f
k
,
w
r
,
s on route r
∈
R
w between O-D pair w
∈
W
k in single city k
∈
K and scenario s
∈
Ω can be given by(19)
f
k
,
w
r
,
s
=
q
k
,
w
s
e
x
p
-
θ
u
k
,
w
r
,
s
∑
r
∈
R
w
e
x
p
-
θ
u
k
,
w
r
,
s
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
s
∈
Ωwhere θ represents the sensitivity of route selection disutility and q
k
,
w
s denotes the logistics demand function. Eq. (19) is the most widely used flow assignment method in traffic planning and represents the flow assigned to each feasible route for each logistic demand [41]. To capture the logistics users’ responses to logistics service disutility, we assume that the elastic demand function between a generic O-D pair w
∈
W
k is a monotonically decreasing function of the O-D service cost u
k
,
w
r
,
s between this O-D pair. The generic elastic demand function [42] is expressed as follows:(20)
q
k
,
w
s
=
D
w
λ
k
,
w
s
=
q
-
k
,
w
s
exp
⁡
-
β
λ
k
,
w
s
,
∀
k
∈
K
,
w
∈
W
k
,
s
∈
Ωwhere β represents the sensitivity to the expected service cost and λ
k
,
w
s is the expected minimum perceived service cost between O-D pair w
∈
W
k in single city k
∈
K and scenario s
∈
Ω. In the case of stochastic user equilibrium (SUE) assignment with elastic demand [42], the expected minimum service cost between an O-D pair under logistics demand scenario s
∈
Ω in city k
∈
K could be expressed as(21)
λ
k
,
w
s
=
-
1
θ
ln
⁡
∑
r
∈
R
w
exp
⁡
-
θ
u
k
,
w
r
,
sEquivalence Model. According to the above analysis, for any given logistics authority decision X
k
,
y
k, the lower-level model of green logistics network design in a single city k can be formulated as the following equivalent minimization program:(22)
m
i
n
Z
F
=
1
θ
∑
w
∈
W
k
∑
r
∈
R
w
f
k
,
w
r
,
s
ln
⁡
f
k
,
w
r
,
s
-
1
-
1
θ
∑
w
∈
W
k
q
k
,
w
s
ln
⁡
q
k
,
w
s
-
1
-
∑
w
∈
W
k
∫
0
q
k
,
w
s
D
w
-
1
w
d
w
+
∑
m
∈
M
0
∑
a
∈
A
k
g
τ
v
o
t
t
k
,
a
m
,
s
v
k
,
a
s
+
c
k
,
a
m
l
k
,
a
+
e
m
l
k
,
a
y
k
v
k
,
a
s
+
∑
i
∈
I
k
∑
a
∈
A
k
t
τ
v
o
t
t
k
,
a
i
f
k
,
i
s
+
c
k
,
a
i
f
k
,
i
ssubject to(23)
∑
r
∈
R
w
f
k
,
w
r
,
s
=
q
k
,
w
s
,
∀
k
∈
K
,
w
∈
W
k
,
s
∈
Ω
(24)
f
k
,
w
r
,
s
>
0
,
∀
k
∈
K
,
w
∈
W
k
,
r
∈
R
w
,
s
∈
Ωwhere(25)
v
k
,
a
s
=
∑
w
∈
W
∑
r
∈
R
w
f
k
,
w
r
,
s
δ
a
,
r
m
,
∀
k
∈
K
,
a
∈
A
k
g
,
s
∈
Ω
(26)
f
k
,
i
s
=
∑
w
∈
W
∑
r
∈
R
w
f
k
,
w
r
,
s
δ
i
,
r
,
∀
k
∈
K
,
i
∈
I
k
,
s
∈
ΩConstraint (23) is the flow conservation constraint. Constraint (24) is the nonnegativity constraints of service route flows. Constraint (25) defines the relationship between the route flow and path flow. Constraint (26) defines the relationship between the route flow and logistics nodes flow, respectively.Proposition 1.
The minimization program (22)-(26) is equivalent to the logit-based stochastic user equilibrium (SUE) assignment with elastic demand under logistics demand scenario s
∈
Ω in city k
∈
K. The proof of Proposition 1 is given in Appendix B.


## 3.3. Centralized Decision in Urban Agglomeration

In the previous subsection, we study the problem that several logistics authorities perform logistics network design decisions independently. However, such independence might be adverse to the overall benefit of urban agglomeration, since it does not explore the full potential of the investment to improve the need of all logistics users on the entire logistics network. We propose a centralized decision model in which logistics authorities cooperate with each other to maximize the total social welfare in urban agglomeration.The upper-level model of green logistics network design in urban agglomeration is formulated as(27)
m
a
x
X
K
+
1
,
y
K
+
1
S
W
X
K
+
1
,
y
K
+
1
,
V
K
+
1
=
∑
s
∈
Ω
p
s
C
S
K
+
1
s
+
P
S
K
+
1,1
s
+
P
S
K
+
1,2
ssubject to(28)
0
≤
x
i
≤
S
i
m
a
x
,
∀
i
∈
P
K
+
1
(29)
0
≤
y
K
+
1
≤
C
h
K
+
1
m
a
x
(30)
∑
i
∈
I
K
+
1
C
i
0
g
x
i
ρ
≤
B
K
+
1The lower level program is still the traffic assignment problem described in Section3.2.2.


## 4. Solution Methods

### 4.1. The Quantum Behaved Particle Swarm Optimization (QPSO)

The QPSO was introduced by Sun et al. [43], combining quantum theory based on the particle swarm optimization (PSO). Different from that in PSO, particles in QPSO have no velocity vectors; in addition, parameters that need to be adjusted are far fewer. It has been widely used in theory and engineering practice because of its global convergence together with comparative simplicity. In the QPSO, the particles are updated with the four following equations:(31)
m
b
e
s
t
t
u
=
m
b
e
s
t
,
1
t
u
,
m
b
e
s
t
,
2
t
u
,
…
,
m
b
e
s
t
,
D
t
u
=
1
M
Q
∑
j
∈
M
Q
P
j
t
u
=
1
M
Q
∑
j
∈
M
Q
P
j
,
1
t
u
,
1
M
Q
∑
j
∈
M
Q
P
j
,
2
t
u
,
…
,
1
M
Q
∑
j
∈
M
Q
P
j
,
D
t
u
(32)
p
j
,
d
t
u
=
φ
j
,
d
·
P
j
,
d
t
u
+
1
-
φ
j
,
d
·
P
g
,
d
t
u
,
∀
j
∈
M
Q
,
d
∈
D
,
φ
j
,
d
~
U
0,1
(33)
α
t
u
=
α
0
-
α
0
-
α
1
·
t
u
T
(34)
X
-
j
,
d
t
u
+
1
=
p
j
,
d
t
u
±
α
t
u
·
m
b
e
s
t
,
d
t
u
-
X
-
j
,
d
t
u
·
l
n
1
μ
j
,
d
t
u
,
∀
j
∈
M
Q
,
d
∈
D
,
μ
j
,
d
~
U
0,1where the X
- denotes the particle’s position. m
b
e
s
t denotes the mean best position of all the particles’ best positions. The P
j and P
g are the particle’s personal best position and the global best position, respectively. α
t
u is the expansion-contraction coefficient, where t is the current iteration number and T is the maximum number of iterations; in general, u
0
=
1, u
1
=
0.5.


### 4.2. Hybrid QPSO Algorithm for Bilevel Decision Problem

To penalize the candidate solutions violating construction investment constraints (12), we first define the evaluation function as the sum of the objective function and penalty terms as follows:(35)
F
X
k
,
y
k
,
V
k
=
S
W
X
k
,
y
k
,
V
k
-
r
e
m
a
x
∑
i
∈
I
k
C
k
,
i
0
x
k
,
i
ρ
-
B
k
,
0
2subject to (10),(11),(13)-(26), and k
∈
K
-wherer
e is the positive variable penalty coefficient.The above model can then be solved by a solution algorithm based on the QPSO and Method of Successive Averages (MSA) hybrid algorithm (Figure3). The detailed process of the above hybrid algorithm is described as follows.Figure 3
Flow chart of QPSO and MSA hybrid algorithm.Step 1.
External Loop Initialization.
(1)
Set the QPSO parameters: population size (M
Q) and related parameters(α
0
,
α
1).
(2)
Set the termination criterion for the external loop operation: maximum iteration (T).
(3)
Set the external loop iteration countert
u
=
0 and randomly generate particle’s initial position X
-
k
,
j
t
u
=
X
k
,
j
t
u
,
y
k
,
j
t
u
,
∀
j
∈
M
Q
,
k
∈
K
-, in which X
-
k
,
j
t
u is sampled randomly in the feasible space. If t
u
=
0, let the personal best (pbest) position of each particle P
k
,
j
t
u
=
X
-
k
,
j
t
u
.Step 2.
Inner Loop Operation (Solve the Lower-Level Flow Assignment).
For a given logistics authority decisionX
-
k
,
j
t
u, use the method of successive averages (MSA) to solve the lower-level flow assignment problem for each scenarios
∈
Ω and obtain optimal solutions V
k
X
-
k
,
j
t
u
.Step 3.
Calculate the mean of the personal best (pbest) positions of all particles m
b
e
s
t
,
k
t
u by using (28).Step 4.
Update the expansion-contraction coefficientα
k
t
u by using (33).Step 5.
For each particle, calculate the personal fitness valueF
X
-
k
,
j
t
u
,
V
k
X
-
k
,
j
t
u andpbest position P
k
,
j
t
u.
IfF
X
-
k
,
j
t
u
,
V
k
X
-
k
,
j
t
u
>
F
P
k
,
j
t
u
,
V
k
X
-
k
,
j
t
u, update P
k
,
j
t
u
=
X
-
k
,
j
t
u
.Step 6.
Update the global best (gbest) position P
k
,
g
t
u
,
g
=
a
r
g
max
j
∈
M
Q
⁡
F
P
k
,
j
t
u
,
V
k
X
-
k
,
j
t
u
.Step 7.
For each particle, generate the new local attractor pointp
k
,
j
,
d
t
u by using (32).Step 8.
For each particle, update the particle’s positionX
-
k
,
j
,
d
t
u
+
1 by using (34).Step 9.
Check the termination criterion for the external loop operation. If (t
u
≥
T), then stop the calculation and go to Step 10; otherwise, let t
u
=
t
u
+
1, and go to Step 2.Step 10.
Output the final optimal solutions.In the upper model, as shown in Figure4, each particle is represented as one string comprising P
k
+
1 continuous numbers, where P
k is the number of entities vector X
k
=
x
1
,
x
2
,
…
,
x
P
k; that is, D
=
P
k
+
1, where D is a set of D consecutive positive integers.Figure 4
Representation of individual particle.The method of successive averages used in Step2 is described in the following substeps.Step
2.1.Inner Loop Initialization.
(1)
Set the termination criterion for the external loop operation (ζ
l).
(2)
Set the initial link flowv
k
,
a
s
0
=
0 and calculate the route flow f
k
,
w
r
,
s
0 by using (19).
(3)
Determine effective service route setR
w
,
w
∈
W
k
.
(4)
Set the inner loop iteration counter (t
l
=
1).Step
2.2. Calculate the link flow v
k
,
a
s
t
l and node flow f
k
,
w
r
,
s
t
l by using (25) and (26), respectively.Step
2.3. Calculate the link service time t
k
,
a
m
,
s
v
k
,
a
m
,
s
t
l and node service time t
k
,
a
i
,
s
f
k
,
i
s
t
l by using (13) and (14), respectively.Step
2.4. Calculate the disutility u
w
r
m
s
t
l by using (15).Step
2.5. Calculate the auxiliary route flow f
¯
w
r
m
s
t
l by using (19).Step
2.6. Update route flow by using method of successive average (MSA).Letf
w
r
m
s
t
l
+
1
=
f
w
r
m
s
t
l
+
1
/
t
l
f
¯
w
r
m
s
t
l
-
f
w
r
m
s
t
lStep
2.7. Check the termination criterion for the inner loop operation.Let(36)
g
a
p
s
t
l
=
∑
m
∈
M
∑
w
∈
W
k
∑
r
∈
R
w
f
w
r
m
s
t
l
+
1
-
f
w
r
m
s
t
l
2
∑
m
∈
M
∑
w
∈
W
k
∑
r
∈
R
w
f
w
r
m
s
t
l
,if g
a
p
s
t
l
≤
ζ
l, then stop the calculation and go to Step 2.8; otherwise, let t
l
=
t
l
+
1 and go to Step 2.2.Step
2.8. Return to the inner loop optimal solutions V
k
X
-
k
,
j
t
u
.


## 4.1. The Quantum Behaved Particle Swarm Optimization (QPSO)

The QPSO was introduced by Sun et al. [43], combining quantum theory based on the particle swarm optimization (PSO). Different from that in PSO, particles in QPSO have no velocity vectors; in addition, parameters that need to be adjusted are far fewer. It has been widely used in theory and engineering practice because of its global convergence together with comparative simplicity. In the QPSO, the particles are updated with the four following equations:(31)
m
b
e
s
t
t
u
=
m
b
e
s
t
,
1
t
u
,
m
b
e
s
t
,
2
t
u
,
…
,
m
b
e
s
t
,
D
t
u
=
1
M
Q
∑
j
∈
M
Q
P
j
t
u
=
1
M
Q
∑
j
∈
M
Q
P
j
,
1
t
u
,
1
M
Q
∑
j
∈
M
Q
P
j
,
2
t
u
,
…
,
1
M
Q
∑
j
∈
M
Q
P
j
,
D
t
u
(32)
p
j
,
d
t
u
=
φ
j
,
d
·
P
j
,
d
t
u
+
1
-
φ
j
,
d
·
P
g
,
d
t
u
,
∀
j
∈
M
Q
,
d
∈
D
,
φ
j
,
d
~
U
0,1
(33)
α
t
u
=
α
0
-
α
0
-
α
1
·
t
u
T
(34)
X
-
j
,
d
t
u
+
1
=
p
j
,
d
t
u
±
α
t
u
·
m
b
e
s
t
,
d
t
u
-
X
-
j
,
d
t
u
·
l
n
1
μ
j
,
d
t
u
,
∀
j
∈
M
Q
,
d
∈
D
,
μ
j
,
d
~
U
0,1where the X
- denotes the particle’s position. m
b
e
s
t denotes the mean best position of all the particles’ best positions. The P
j and P
g are the particle’s personal best position and the global best position, respectively. α
t
u is the expansion-contraction coefficient, where t is the current iteration number and T is the maximum number of iterations; in general, u
0
=
1, u
1
=
0.5.


## 4.2. Hybrid QPSO Algorithm for Bilevel Decision Problem

To penalize the candidate solutions violating construction investment constraints (12), we first define the evaluation function as the sum of the objective function and penalty terms as follows:(35)
F
X
k
,
y
k
,
V
k
=
S
W
X
k
,
y
k
,
V
k
-
r
e
m
a
x
∑
i
∈
I
k
C
k
,
i
0
x
k
,
i
ρ
-
B
k
,
0
2subject to (10),(11),(13)-(26), and k
∈
K
-wherer
e is the positive variable penalty coefficient.The above model can then be solved by a solution algorithm based on the QPSO and Method of Successive Averages (MSA) hybrid algorithm (Figure3). The detailed process of the above hybrid algorithm is described as follows.Figure 3
Flow chart of QPSO and MSA hybrid algorithm.Step 1.
External Loop Initialization.
(1)
Set the QPSO parameters: population size (M
Q) and related parameters(α
0
,
α
1).
(2)
Set the termination criterion for the external loop operation: maximum iteration (T).
(3)
Set the external loop iteration countert
u
=
0 and randomly generate particle’s initial position X
-
k
,
j
t
u
=
X
k
,
j
t
u
,
y
k
,
j
t
u
,
∀
j
∈
M
Q
,
k
∈
K
-, in which X
-
k
,
j
t
u is sampled randomly in the feasible space. If t
u
=
0, let the personal best (pbest) position of each particle P
k
,
j
t
u
=
X
-
k
,
j
t
u
.Step 2.
Inner Loop Operation (Solve the Lower-Level Flow Assignment).
For a given logistics authority decisionX
-
k
,
j
t
u, use the method of successive averages (MSA) to solve the lower-level flow assignment problem for each scenarios
∈
Ω and obtain optimal solutions V
k
X
-
k
,
j
t
u
.Step 3.
Calculate the mean of the personal best (pbest) positions of all particles m
b
e
s
t
,
k
t
u by using (28).Step 4.
Update the expansion-contraction coefficientα
k
t
u by using (33).Step 5.
For each particle, calculate the personal fitness valueF
X
-
k
,
j
t
u
,
V
k
X
-
k
,
j
t
u andpbest position P
k
,
j
t
u.
IfF
X
-
k
,
j
t
u
,
V
k
X
-
k
,
j
t
u
>
F
P
k
,
j
t
u
,
V
k
X
-
k
,
j
t
u, update P
k
,
j
t
u
=
X
-
k
,
j
t
u
.Step 6.
Update the global best (gbest) position P
k
,
g
t
u
,
g
=
a
r
g
max
j
∈
M
Q
⁡
F
P
k
,
j
t
u
,
V
k
X
-
k
,
j
t
u
.Step 7.
For each particle, generate the new local attractor pointp
k
,
j
,
d
t
u by using (32).Step 8.
For each particle, update the particle’s positionX
-
k
,
j
,
d
t
u
+
1 by using (34).Step 9.
Check the termination criterion for the external loop operation. If (t
u
≥
T), then stop the calculation and go to Step 10; otherwise, let t
u
=
t
u
+
1, and go to Step 2.Step 10.
Output the final optimal solutions.In the upper model, as shown in Figure4, each particle is represented as one string comprising P
k
+
1 continuous numbers, where P
k is the number of entities vector X
k
=
x
1
,
x
2
,
…
,
x
P
k; that is, D
=
P
k
+
1, where D is a set of D consecutive positive integers.Figure 4
Representation of individual particle.The method of successive averages used in Step2 is described in the following substeps.Step
2.1.Inner Loop Initialization.
(1)
Set the termination criterion for the external loop operation (ζ
l).
(2)
Set the initial link flowv
k
,
a
s
0
=
0 and calculate the route flow f
k
,
w
r
,
s
0 by using (19).
(3)
Determine effective service route setR
w
,
w
∈
W
k
.
(4)
Set the inner loop iteration counter (t
l
=
1).Step
2.2. Calculate the link flow v
k
,
a
s
t
l and node flow f
k
,
w
r
,
s
t
l by using (25) and (26), respectively.Step
2.3. Calculate the link service time t
k
,
a
m
,
s
v
k
,
a
m
,
s
t
l and node service time t
k
,
a
i
,
s
f
k
,
i
s
t
l by using (13) and (14), respectively.Step
2.4. Calculate the disutility u
w
r
m
s
t
l by using (15).Step
2.5. Calculate the auxiliary route flow f
¯
w
r
m
s
t
l by using (19).Step
2.6. Update route flow by using method of successive average (MSA).Letf
w
r
m
s
t
l
+
1
=
f
w
r
m
s
t
l
+
1
/
t
l
f
¯
w
r
m
s
t
l
-
f
w
r
m
s
t
lStep
2.7. Check the termination criterion for the inner loop operation.Let(36)
g
a
p
s
t
l
=
∑
m
∈
M
∑
w
∈
W
k
∑
r
∈
R
w
f
w
r
m
s
t
l
+
1
-
f
w
r
m
s
t
l
2
∑
m
∈
M
∑
w
∈
W
k
∑
r
∈
R
w
f
w
r
m
s
t
l
,if g
a
p
s
t
l
≤
ζ
l, then stop the calculation and go to Step 2.8; otherwise, let t
l
=
t
l
+
1 and go to Step 2.2.Step
2.8. Return to the inner loop optimal solutions V
k
X
-
k
,
j
t
u
.


## 5. Numerical Examples

As shown in Figure5, an example multimodal logistics network of urban agglomeration is used to illustrate the availability of the proposed bilevel model and QPSO algorithm. The solution procedure is coded by MATLAB R2012a and runs on a desktop Lenovo G3240 with an Intel Pentium 3.10 GHz and 4 GB RAM.Figure 5
An example multimodal logistics network of urban agglomeration.


### 5.1. Main Data and Parameter Values

The network comprises 19 nodes, 42 arcs, and six logistics demand O-D pairs. The logistics demand for each O-D pair, the basic data of each logistics node, and each arc in the network are provided in AppendixC. To examine the effect of logistics demand uncertainty on logistics network design decisions, we consider three logistics demand scenarios. The total construction budget for urban agglomeration is 15,000 $/week. The sensitivity parameters θ and β are 0.8 and 0.001, respectively. The value of time is 8 $/h. The unit investment subsidy for each logistics node is 0.5 $/week. The average emissions of expressway, railway, and waterway are 0.132, 0.022, and 0.016 kg/ton-km, respectively. The maximum carbon tax on arc is 1 $/kg. Unless otherwise specified, these input data are considered unchanged in the following analysis.


### 5.2. Numerical Results and Discussions

For comparison, the do-nothing model and the decentralized and centralized decision models are considered in this paper. In the do-nothing model and decentralized decision model, the logistics nodes budget for each single city is set as 0 and 5,000 $/week, respectively. However, in the centralized decision model, the budget is designed under a total budget constraint, 15,000 $/week. The optimal solutions and performance comparison for these three models are given in Tables2 and 3, respectively.Table 2
Optimal solutions for three logistics network design models.
Models
City 1
City 2
City 3
X1
X2
X7
Y1
X3
X4
X8
X9
Y2
X5
X6
X10
Y3
Do-nothing
0
0
0
0
0
0
0
0
0
0
0
0
0
Decentralized
Tax
662
498
0
0.22
602
585
0
0
0.17
205
1000
0
0.24
Nontax
273
939
0
0
720
468
0
0
0
182
988
0
0
Centralized
Tax
630
495
0
0.21
829
704
0
0
0.21
0
1000
0
0.21
Nontax
229
994
0
0
750
742
0
0
0
0
952
0
0Table 3
Comparison of the three logistics network design model performances.
Models
Total social welfare ($)
Consumer surplus ($)
Producer surplus ($)
Percentage of combined transport
Average carbon emission rate (kg/ton-km)
Do-nothing
2,449,637
2,345,459
104,178
0%
0.1320
Decentralized
Tax
3,889,751
3,707,380
182,371
41.88%
0.0739
Nontax
3,879,730
3,698,102
181,628
36.08%
0.0752
Centralized
Tax
3,950,308
3,765,081
185,227
42.00 %
0.0730
Nontax
3,927,741
3,743,486
184,255
36.36%
0.0741Table2 clearly shows that under the given investment budget five logistics parks located at (X1-X4, X6) are constructed in the centralized decision model, and the total processing capacity is 3,658 tons/week. While the number of logistics parks increases to six (X1-X6) in the decentralized decision model, the total processing capacity decreases to 3,552 tons/week. By comparing the centralized decision model with the decentralized decision model, it shows that in the decentralized decision model the investment budget for the city 1 is reasonable, while the investment budget for the city 2 is insufficient, and the investment budget for the city 3 is in surplus. Similarly, as shown in Table 3, among the three design models, the centralized decision model has the best performance, followed by the decentralized decision model. Compared with the do-nothing model, the centralized decision model and the decentralized decision model can reduce the average carbon emission rate by approximately 45% after the construction of a multimode network.Given the above, to achieve overall optimality in the layout of urban agglomeration logistics networks, each logistics authority in the urban agglomeration should make joint decisions, allocate the investment budget, and plan the quantity and scale of the logistics nodes rationally. At the same time, the carbon emissions generated during the logistics and transportation can be greatly reduced by constructing a multimodal network and setting a reasonable carbon tax.


### 5.3. Impact Analysis of Budget Allocation on Network Performance

As seen from Table2, the investment budget of city one is reasonable. To analyze the impact of budget distribution on network performance in the decentralized decision model, we set B
1
=
5,000 $/week and carried out numerical experiments of nine budget allocation scenarios (with the same total investment budget, B
2
+
B
3
=
10,000 $/week). Numerical results are presented in Table 4 and Figure 6.Table 4
Optimal solutions for decentralized decision model under different budget allocation.
Scenario
City 1
City 2
City 3
($/week)
X1
X2
X7
Y1
X3
X4
X8
X9
Y2
X5
X6
X10
Y3
(5-1-9)∗103
662
498
0
0.22
199
0
0
0
0.54
1000
1000
10
0.03
(5-2-8)∗103
662
498
0
0.22
450
0
0
0
0.48
1000
1000
0
0.04
(5-3-7)∗103
662
498
0
0.22
579
114
0
0
0.42
737
1000
0
0.08
(5-4-6)∗103
662
498
0
0.22
590
341
0
0
0.33
470
1000
0
0.14
(5-5-5)∗103
662
498
0
0.22
602
585
0
0
0.17
205
1000
0
0.24
(5-6-4)∗103
662
498
0
0.22
806
654
0
0
0.16
0
999
0
0.37
(5-7-3)∗103
662
498
0
0.22
994
738
0
0
0.12
0
728
0
0.51
(5-8-2)∗103
662
498
0
0.22
1000
1000
0
0
0.05
0
465
0
0.63
(5-9-1)∗103
662
498
0
0.22
1000
1000
0
74
0.04
0
216
0
0.84Figure 6
Social welfare for each city under different budget allocation.From Table4, first, we can see that the construction of the logistics park takes priority over the general logistics node. Second, with the increase in the investment budget in city 2, the number of logistics parks and logistics nodes will increase, and processing capacity will be enhanced gradually. To be specific, when the investment budget is equal to 1,000-2,000$/week, one logistics park is constructed; when the investment budget is equal to 3,000-8,000$/week, two logistics parks are constructed; and when the investment budget is more than 9,000$/week, two logistics parks and one logistics node are constructed. Finally, the optimal carbon tax gradually declines with the increase in the investment volume. This is because before the construction of the logistics park, there is a single highway with high carbon emissions in the logistics network, and after the construction of the logistics park there will be a single highway and multimodal transportation such as highway-railway and highway-waterway in the logistics network, which effectively reduces the carbon tax in the network. As shown in Figure 6, the social welfare increases with the increase in the investment budget for a single city. However, from the perspective of total social welfare, when the investment budget is allocated in the scenario of 5,000-6,000-4,000 $/week, the total social welfare is the maximum for the urban agglomeration.


### 5.4. Impact Analysis of Mode Choices on Network Performance

To fully explore the impact of mode choices on network performance, four mode combination strategies are considered: (I) Expressway, (II) Expressway-Railway, (III) Expressway-Waterway, and (IV) Expressway-Railway-Waterway. The first strategy is to simulate the traditional single-road transport structure. Considering the incomplete accessibility of other transportation modes, this paper reasonably constructs the latter three multimodal strategies based on road transportation. The network performance of different mode strategies under different budget allocation scenarios is shown in Table5.Table 5
Network performance of different mode strategies under different budget allocation.
Scenario
Total social welfare ($)
Average carbon emission rate (kg/ton-km)
($/week)
I
II
III
IV
I
II
III
IV
(5-1-9)∗103
2,463,707
3,292,622
3,263,794
3,558,279
0.1320
0.0880
0.0900
0.0797
(5-2-8)∗103
2,463,707
3,420,352
3,372,572
3,615,201
0.1320
0.0841
0.0870
0.0778
(5-3-7)∗103
2,463,707
3,568,935
3,511,567
3,777,362
0.1320
0.0805
0.0839
0.0745
(5-4-6)∗103
2,463,707
3,666,169
3,634,317
3,848,109
0.1320
0.0782
0.0824
0.0751
(5-5-5)∗103
2,463,707
3,666,167
3,634,048
3,889,751
0.1320
0.0782
0.0815
0.0739
(5-6-4)∗103
2,463,707
3,666,229
3,633,969
3,914,493
0.1320
0.0782
0.0815
0.0738
(5-7-3)∗103
2,463,707
3,579,201
3,571,180
3,912,885
0.1320
0.0805
0.0832
0.0740
(5-8-2)∗103
2,463,707
3,472,760
3,473,270
3,866,673
0.1320
0.0835
0.0858
0.0752
(5-9-1)∗103
2,463,707
3,349,296
3,357,314
3,756,729
0.1320
0.0870
0.0891
0.0776As can be seen from Table5, in all budget allocation scenarios, the total social welfares under multimodal transport strategies II, III, and IV are significantly better than the traditional single-road strategy I, with an average increase of 42.88%, 41.85%, and 53.97%, respectively. In this context, the total social welfare under the strategy IV is better than the strategies II and III, with an average increase of 7.81% and 8.58%, respectively. Furthermore, strategy III is slightly inferior to strategy II in most scenarios. To reduce the impact on the environment, the logistics authority is also concerned about the average carbon emission rate of the urban agglomeration. Table 5 also shows that the average carbon emission rate under different multimodal transport strategies can effectively reduce CO2 emissions. Compared with the strategy I, adopting multimodal transport strategy II, III, and IV results in an average decrease of 61.23%, 55.61%, and 74.41%, respectively. Therefore, the development of multimodal transport network can not only improve social welfare but also effectively reduce CO2 emissions. In addition, other realistic conditions of each strategy should be considered when selecting a multimodal strategy. For example, taking strategy I as the benchmark, there is an expected target to reduce the average carbon emission rate in the regional multimodal network design. If the target is set to 50%, all multimodal strategies are feasible. If set to 60%, strategies II and IV are feasible. If set to 70%, only strategy IV is feasible.


## 5.1. Main Data and Parameter Values

The network comprises 19 nodes, 42 arcs, and six logistics demand O-D pairs. The logistics demand for each O-D pair, the basic data of each logistics node, and each arc in the network are provided in AppendixC. To examine the effect of logistics demand uncertainty on logistics network design decisions, we consider three logistics demand scenarios. The total construction budget for urban agglomeration is 15,000 $/week. The sensitivity parameters θ and β are 0.8 and 0.001, respectively. The value of time is 8 $/h. The unit investment subsidy for each logistics node is 0.5 $/week. The average emissions of expressway, railway, and waterway are 0.132, 0.022, and 0.016 kg/ton-km, respectively. The maximum carbon tax on arc is 1 $/kg. Unless otherwise specified, these input data are considered unchanged in the following analysis.


## 5.2. Numerical Results and Discussions

For comparison, the do-nothing model and the decentralized and centralized decision models are considered in this paper. In the do-nothing model and decentralized decision model, the logistics nodes budget for each single city is set as 0 and 5,000 $/week, respectively. However, in the centralized decision model, the budget is designed under a total budget constraint, 15,000 $/week. The optimal solutions and performance comparison for these three models are given in Tables2 and 3, respectively.Table 2
Optimal solutions for three logistics network design models.
Models
City 1
City 2
City 3
X1
X2
X7
Y1
X3
X4
X8
X9
Y2
X5
X6
X10
Y3
Do-nothing
0
0
0
0
0
0
0
0
0
0
0
0
0
Decentralized
Tax
662
498
0
0.22
602
585
0
0
0.17
205
1000
0
0.24
Nontax
273
939
0
0
720
468
0
0
0
182
988
0
0
Centralized
Tax
630
495
0
0.21
829
704
0
0
0.21
0
1000
0
0.21
Nontax
229
994
0
0
750
742
0
0
0
0
952
0
0Table 3
Comparison of the three logistics network design model performances.
Models
Total social welfare ($)
Consumer surplus ($)
Producer surplus ($)
Percentage of combined transport
Average carbon emission rate (kg/ton-km)
Do-nothing
2,449,637
2,345,459
104,178
0%
0.1320
Decentralized
Tax
3,889,751
3,707,380
182,371
41.88%
0.0739
Nontax
3,879,730
3,698,102
181,628
36.08%
0.0752
Centralized
Tax
3,950,308
3,765,081
185,227
42.00 %
0.0730
Nontax
3,927,741
3,743,486
184,255
36.36%
0.0741Table2 clearly shows that under the given investment budget five logistics parks located at (X1-X4, X6) are constructed in the centralized decision model, and the total processing capacity is 3,658 tons/week. While the number of logistics parks increases to six (X1-X6) in the decentralized decision model, the total processing capacity decreases to 3,552 tons/week. By comparing the centralized decision model with the decentralized decision model, it shows that in the decentralized decision model the investment budget for the city 1 is reasonable, while the investment budget for the city 2 is insufficient, and the investment budget for the city 3 is in surplus. Similarly, as shown in Table 3, among the three design models, the centralized decision model has the best performance, followed by the decentralized decision model. Compared with the do-nothing model, the centralized decision model and the decentralized decision model can reduce the average carbon emission rate by approximately 45% after the construction of a multimode network.Given the above, to achieve overall optimality in the layout of urban agglomeration logistics networks, each logistics authority in the urban agglomeration should make joint decisions, allocate the investment budget, and plan the quantity and scale of the logistics nodes rationally. At the same time, the carbon emissions generated during the logistics and transportation can be greatly reduced by constructing a multimodal network and setting a reasonable carbon tax.


## 5.3. Impact Analysis of Budget Allocation on Network Performance

As seen from Table2, the investment budget of city one is reasonable. To analyze the impact of budget distribution on network performance in the decentralized decision model, we set B
1
=
5,000 $/week and carried out numerical experiments of nine budget allocation scenarios (with the same total investment budget, B
2
+
B
3
=
10,000 $/week). Numerical results are presented in Table 4 and Figure 6.Table 4
Optimal solutions for decentralized decision model under different budget allocation.
Scenario
City 1
City 2
City 3
($/week)
X1
X2
X7
Y1
X3
X4
X8
X9
Y2
X5
X6
X10
Y3
(5-1-9)∗103
662
498
0
0.22
199
0
0
0
0.54
1000
1000
10
0.03
(5-2-8)∗103
662
498
0
0.22
450
0
0
0
0.48
1000
1000
0
0.04
(5-3-7)∗103
662
498
0
0.22
579
114
0
0
0.42
737
1000
0
0.08
(5-4-6)∗103
662
498
0
0.22
590
341
0
0
0.33
470
1000
0
0.14
(5-5-5)∗103
662
498
0
0.22
602
585
0
0
0.17
205
1000
0
0.24
(5-6-4)∗103
662
498
0
0.22
806
654
0
0
0.16
0
999
0
0.37
(5-7-3)∗103
662
498
0
0.22
994
738
0
0
0.12
0
728
0
0.51
(5-8-2)∗103
662
498
0
0.22
1000
1000
0
0
0.05
0
465
0
0.63
(5-9-1)∗103
662
498
0
0.22
1000
1000
0
74
0.04
0
216
0
0.84Figure 6
Social welfare for each city under different budget allocation.From Table4, first, we can see that the construction of the logistics park takes priority over the general logistics node. Second, with the increase in the investment budget in city 2, the number of logistics parks and logistics nodes will increase, and processing capacity will be enhanced gradually. To be specific, when the investment budget is equal to 1,000-2,000$/week, one logistics park is constructed; when the investment budget is equal to 3,000-8,000$/week, two logistics parks are constructed; and when the investment budget is more than 9,000$/week, two logistics parks and one logistics node are constructed. Finally, the optimal carbon tax gradually declines with the increase in the investment volume. This is because before the construction of the logistics park, there is a single highway with high carbon emissions in the logistics network, and after the construction of the logistics park there will be a single highway and multimodal transportation such as highway-railway and highway-waterway in the logistics network, which effectively reduces the carbon tax in the network. As shown in Figure 6, the social welfare increases with the increase in the investment budget for a single city. However, from the perspective of total social welfare, when the investment budget is allocated in the scenario of 5,000-6,000-4,000 $/week, the total social welfare is the maximum for the urban agglomeration.


## 5.4. Impact Analysis of Mode Choices on Network Performance

To fully explore the impact of mode choices on network performance, four mode combination strategies are considered: (I) Expressway, (II) Expressway-Railway, (III) Expressway-Waterway, and (IV) Expressway-Railway-Waterway. The first strategy is to simulate the traditional single-road transport structure. Considering the incomplete accessibility of other transportation modes, this paper reasonably constructs the latter three multimodal strategies based on road transportation. The network performance of different mode strategies under different budget allocation scenarios is shown in Table5.Table 5
Network performance of different mode strategies under different budget allocation.
Scenario
Total social welfare ($)
Average carbon emission rate (kg/ton-km)
($/week)
I
II
III
IV
I
II
III
IV
(5-1-9)∗103
2,463,707
3,292,622
3,263,794
3,558,279
0.1320
0.0880
0.0900
0.0797
(5-2-8)∗103
2,463,707
3,420,352
3,372,572
3,615,201
0.1320
0.0841
0.0870
0.0778
(5-3-7)∗103
2,463,707
3,568,935
3,511,567
3,777,362
0.1320
0.0805
0.0839
0.0745
(5-4-6)∗103
2,463,707
3,666,169
3,634,317
3,848,109
0.1320
0.0782
0.0824
0.0751
(5-5-5)∗103
2,463,707
3,666,167
3,634,048
3,889,751
0.1320
0.0782
0.0815
0.0739
(5-6-4)∗103
2,463,707
3,666,229
3,633,969
3,914,493
0.1320
0.0782
0.0815
0.0738
(5-7-3)∗103
2,463,707
3,579,201
3,571,180
3,912,885
0.1320
0.0805
0.0832
0.0740
(5-8-2)∗103
2,463,707
3,472,760
3,473,270
3,866,673
0.1320
0.0835
0.0858
0.0752
(5-9-1)∗103
2,463,707
3,349,296
3,357,314
3,756,729
0.1320
0.0870
0.0891
0.0776As can be seen from Table5, in all budget allocation scenarios, the total social welfares under multimodal transport strategies II, III, and IV are significantly better than the traditional single-road strategy I, with an average increase of 42.88%, 41.85%, and 53.97%, respectively. In this context, the total social welfare under the strategy IV is better than the strategies II and III, with an average increase of 7.81% and 8.58%, respectively. Furthermore, strategy III is slightly inferior to strategy II in most scenarios. To reduce the impact on the environment, the logistics authority is also concerned about the average carbon emission rate of the urban agglomeration. Table 5 also shows that the average carbon emission rate under different multimodal transport strategies can effectively reduce CO2 emissions. Compared with the strategy I, adopting multimodal transport strategy II, III, and IV results in an average decrease of 61.23%, 55.61%, and 74.41%, respectively. Therefore, the development of multimodal transport network can not only improve social welfare but also effectively reduce CO2 emissions. In addition, other realistic conditions of each strategy should be considered when selecting a multimodal strategy. For example, taking strategy I as the benchmark, there is an expected target to reduce the average carbon emission rate in the regional multimodal network design. If the target is set to 50%, all multimodal strategies are feasible. If set to 60%, strategies II and IV are feasible. If set to 70%, only strategy IV is feasible.


## 6. Conclusion and Future Studies

This paper studied a multimodal green logistics network design problem of urban agglomeration with four marked features: stochastic demand, congestion effects, game behaviors of multiple logistics authorities, and route choice behaviors of logistics users. To describe the game behaviors of logistics authorities in urban agglomeration, two nonlinear bilevel programming models were proposed, decentralized, and centralized decision models. Compared with the centralized decision model, the decentralized decision model is mainly used to analyze the impact of budget allocation of urban agglomeration. At the upper level, each logistics authority aims to maximize the total social welfare of its own city by planning the investment capacity of logistics nodes and determining the carbon tax on transportation service links. The lower-level subproblem is a logit-based stochastic user equilibrium (SUE) problem, regarding the logistics demand scenario and carbon tax. Next, to solve the nonlinear bilevel programming models, a quantum-behaved particle swarm optimization (QPSO) embedded with a Method of Successive Averages (MSA) for stochastic user equilibrium flow assignment was proposed.Numerical examples show that centralized decisions are the optimal decisions for improving the entire logistics network performance, followed by decentralized decisions among multiple logistics authorities under a budget allocation constraint. Meanwhile, compared with the single mode logistics network, the constructed multimode logistics network can reduce the average carbon emission rate by approximately 45%. In addition, by adopting centralized decisions, constructing a multimode logistics network and formulating a reasonable carbon tax, the logistics authority can achieve the overall optimization layout of the green logistics network in urban agglomeration. In short, the proposed models and algorithms can help the logistics authorities make scientific decisions on the multimodal green logistics network design problem of urban agglomeration.On the basis of this study, some main extensions can be made in the future, such as applying our proposed models and algorithms to a large and realistic logistics network and establishing multiobjective optimization model or robust optimization model.



---

*Source: 4165942-2019-08-07.xml*