# Retrofitting Transportation Network Using a Fuzzy Random Multiobjective Bilevel Model to Hedge against Seismic Risk

Authors: Lu Gan; Jiuping Xu
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This paper focuses on the problem of hedging against seismic risk through the retrofit of transportation systems in large-scale construction projects (LSCP). A fuzzy random multiobjective bilevel programming model is formulated with the objectives of the retrofit costs and the benefits on two separate levels. After establishing the model, a fuzzy random variable transformation approach and fuzzy variable approximation decomposition are used to deal with the uncertainty. An approximation decomposition-based multi-objective AGLNPSO is developed to solve the model. The results of a case study validate the efficiency of the proposed approach.


---

## Body

## 1. Introduction

Transportation networks play a very important role in both urban and rural areas, as well as in industrial sites such as large-scale construction sites. Liu et al. [1] stated that transportation networks are critical infrastructure and their smooth operation is important for maintaining the normal functions of society. However, disasters, especially earthquakes, cause not only tremendous economic losses and social chaos but also enormous damage to infrastructure (e.g., 2008 Wenchuan Earthquake, 2010 Chile Earthquake, and 2011 Japan Earthquake). Thus, as Liu et al. [1] pointed out, seismic risk control should also consider the effect that damaged or destroyed transportation networks have on the effectiveness of postdisaster rescue and repair activities and the associated socioeconomic losses. Under a seismic risk threat, retrofit decisions are considered to be an effective protective measure and can have a significant impact on these systems [1–3]. Therefore, promoting retrofit decisions for transportation networks is necessary to hedge against seismic risk.The research in this area has mainly focused on the retrofitting of bridges for transportation networks [4–6]. Werner et al. [2] extended seismic retrofits to highway systems. Afterwards, Liu et al. [1] established a two-stage stochastic programming model for retrofit decisions for transportation network protection. This previous research, however, has primarily focused on urban transportation, but it is essential that transportation networks in large-scale construction projects (LSCP) also be considered. As a critical infrastructure, the smooth operation of these networks is important for maintaining the normal progress of these projects. Therefore, it is necessary to control the seismic risk for LSCP transportation networks to mitigate losses. When considering LSCP transportation network retrofits, there are significant challenges. First, these transportation systems have not only permanent links and temporary links to consider but must also assess the critical links (i.e., bridge, tunnel, etc.) and the noncritical links. Secondly, the retrofit decision making environment is a mutual environment involving an investor who pays for the retrofit and an administrator who controls the transportation systems. Thirdly, a consideration of the environmental costs for the investor has increasingly become necessary for social and economic development. Lastly, a majority of the previous research has assumed that seismic damage is classified into five categories and there is a set of discrete probabilities associated with each of the five damage categories. In practice, however, the situation is often not that simple, and the description of the possible result of seismic damage is vague and uncertain. In this case, this needs to be qualified with a vague perception of a crisp but unobservable random variable. Hence, due to the complexity of assessing the seismic risk to property, seismic damage is subject to uncertainty with both fuzziness and randomness, that is, fuzzy random in nature. More recently, since Kwakernaak [7] proposed the concept of the fuzzy random variable, considerable research has been done, which has allowed for its application in many areas [7–13]. Unfortunately, there has been little research which has discussed a mixture of fuzziness and randomness in a transportation network retrofit problem. Therefore, the uncertainty with hybrid fuzziness and randomness induced by the seismic damage risk to property needs to be further studied and elaborated.The fuzzy random variable was proposed by Kwakernaak [7] who regarded it as “random variables whose values are not real, but fuzzy numbers.” From another view, Puri and Ralescu [14] and Klement et al. [15] regarded a fuzzy random variable as a random fuzzy set. Fuzzy random variables represent a well-formalized concept which has underlain many recent probabilistic and statistical studies involving data obtained from a random experiment when these data are assumed to be fuzzy set valued [16]. Therefore, in a transportation network retrofit problem, the description of seismic damage is considered a fuzzy random variable, that is, a discrete distribution variable with a vague perception (i.e., triangular fuzzy number). Several research works have demonstrated how these fuzzy random coefficients can be converted into crisp values. Usually, at first, the fuzzy random variables are transformed into fuzzy numbers using the fuzzy expected values [17] or transformed into (
α
1
,
σ
)-level trapezoidal fuzzy variables through an approach proposed in Xu and Liu [12]. Then, these fuzzy numbers are transformed into deterministic values using their expected value [18] or (
α
,
β
)-satisfactory solution to the programming is determined using fuzzy coefficients [12]. In this case, based on the properties of the fuzzy random seismic damage in this study, the theorem and the proof presented in Xu and Liu [12] are adjusted to allow for a discrete random distribution to obtain the equivalent fuzzy bilevel programming model. Then, using the theorem proposed by Zhang et al. [19], decomposition is used on these fuzzy variables to derive an approximate solution to the model.Under these emerging challenges, this paper formulates a fuzzy random multiobjective bilevel programming model for a transportation network retrofit decision to hedge against seismic risk in an LSCP. The distinctions in the link types allow for the recognition of the retrofit and reconstruction costs. The investor and the administrator are the decision-makers on two separate levels. Retrofit costs which include the environmental costs and the retrofit benefits are the two objectives of the investor, and the retrofit benefits are the objective of the administrator. In order to describe the hybrid uncertainty of possible seismic damage, fuzzy random variables are introduced in the programming model, the use of which has been applied in many areas [10, 18]. To cope with the proposed fuzzy random multiobjective bilevel programming model, a transformation approach is used to obtain an equivalent fuzzy bilevel programming model. This approach transforms the fuzzy random variables in the model into fuzzy variables which are similar to trapezoidal fuzzy variables. Then, decomposition is utilized on these fuzzy variables using a fuzzy number decomposition theorem [19]. To solve the model, an approximation decomposition-based multiobjective AGLNPSO is developed in this paper. Through the decomposition of the fuzzy variables, the models are successively solved until termination, and the approximation solutions are obtained. The multiobjective AGLNPSO is a combination of the Pareto Archived Evolution Strategy (PAES) [20] and the AGLNPSO [21] which is developed by incorporating an adaptive particle swarm optimization (APSO) [22] with a GNLPSO [23] and a multiple objectives particle swarm optimization (MOPSO) [24].This study contributes to the literature by adopting the work of Liu et al. [1] to an LSCP and describing the complex uncertain seismic damage scenario using fuzzy random variables. Bilevel decisions involving the investor and the administrator, distinctions between the various link types, and the specification of the retrofit decisions into several ranks according to the seismic damage scenario provide a more reasonable and practical description of the problem. The consideration of the environmental costs in the transportation network in an LSCP enhances the focus for management. To the best of our knowledge, an integrated approach to deal with fuzzy random variables has not been previously comprehensively studied. The approximation decomposition-based multiobjective AGLNPSO is developed as a useful tool to solve the problem, in which both the bilevel and multiobjective environments are considered.The remainder of this paper is as follows. The problem description, the fuzzy random multiobjective bilevel programming model, the transformation approach, and the approximation decomposition are given in Section2. An approximation decomposition-based multiobjective AGLNPSO is developed in Section 3. A case study is presented in Section 4. Finally, advantages, limitations, and possible future extensions of this work are presented in Section 5.


## 2. Modeling

In this section, the concepts for the LSCP transportation network, the bilevel decision framework, the environmental costs, and the fuzzy random seismic damage scenario are introduced. A multiobjective bilevel programming model for the problem considering fuzziness and randomness is established. See in the following the notations used to describe the model.Index
a
:
Link in transportation network,a
∈
A
b:
Node in transportation network,b
∈
B
v:
Variable environment cost,v
∈
V
f:
Fixed environmental cost,f
∈
F
i:
Retrofit output,i
∈
I
j:
Retrofit activity,j
∈
J
k:
Origin-destination pair considered as commodity,k
∈
{
1
,
…
,
K
}.Variables
m
a
=
{
1
,
Permanent
link
0
,
Temporary
link
n
a
=
{
1
,
Critical
link
0
,
Noncritical
link
c
v
a
p: Increased variable retrofit costs for permanent link by basic rank (i.e., rank 1)
c
v
a
t: Variable retrofit costs for temporary link by basic rank (i.e., rank 1)
c
f
i
p: Increased fixed retrofit costs for permanent link
c
f
i
t: Fixed retrofit costs for temporary link
ρ: Weight of environmental costs
c
e
v
p: Increased variable environmental costs for permanent link by basic rank (i.e., rank 1)
c
e
v
t: Variable environmental costs for temporary link by basic rank (i.e., rank 1)
c
e
f
f: Fixed environmental costs
p
e
j
v
v: Percent of activity cost center j in variable environment cost v
p
e
i
f
f: Percent of output i in fixed environment cost f
c
e
j
c: Variable environmental costs of activity cost center j
a
m
j: Cost of driver at activity cost center j
r
a
j: Driver rate at activity cost centerj
a
m
i
j: Cost of driver for output i in activity cost center j
C: Retrofit costs including environmental costs
Q: Retrofit benefits
ξ
¯
~
a: Preretrofit link damage state for link a
Ξ
¯
~
a: Postretrofit link damage state for linka
c
r
v
a
p: Increased variable reconstruction cost for permanent link by basic rank (i.e., rank 1)
c
r
v
a
t: Variable reconstruction cost for temporary link by basic rank (i.e., rank 1)
c
r
f
i
p: Increased fixed reconstruction cost for permanent link
c
r
f
i
t: Fixed reconstruction cost for temporary link
γ: Weight coefficient conversion time to monetary value
t
i
a
0: Free flow travel time and link a
α: Coefficient of BPR function
f
l
a: The total flow on link a
c
a
a
′: Practical capacity of link a is set at 90% of the design capacity
β: BPR function coefficient
c
a
b: Capacity of node b
W: Node-commodity adjacency matrix
M: Link-commodity adjacency matrix.Decision Variables
(1)
u
a
:
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
x
k
:
≤
0
,
∀
k
=
1
,
…
,
K
.


### 2.1. LSCP Transportation Network

The LSCP transportation network is composed of an internal road system and an external road system and is always built based on the existing links around the site which are connected with the newly built links according to transportation need. There are two types of links (i.e., permanent links and temporary links) which vary considerably in terms of quality. In addition, according to the different functions, the links are divided into critical links and noncritical links. Critical links are those which have vital transport functions such as bridges and tunnels and which should be preferentially taken into account [1]. The retrofit decisions for the different link types vary. That is, the links being considered for the retrofit are either considered to be permanent or critical. The retrofit and reconstruction costs for the temporary links are lower than those for the permanent links. Further, the retrofit decision is specific with 0 (i.e., no retrofit) and there are several ranks according to the seismic damage scenarios.


### 2.2. Bilevel Decision Framework

In this paper, the seismic hazard retrofit decision for an LSCP transportation network involves two participants (i.e., the investor who pays for the retrofit and the administrator who controls the transportation). Therefore, these two participants are the decision-makers on two levels, both of whom successively make the retrofit decisions. The investor on upper level decides which retrofit rank should be taken for each link within the range and, therefore, the two objectives on this level are the retrofit costs including the environmental costs and the retrofit benefits. The administrator decides on the commodity flow (i.e., the transportation network flow once seismic damage has occurred) on the lower level according to the decision results of the upper level. On this level, the retrofit benefits are the primary objective. In this paper, the retrofit benefits are quantified as reconstruction and travel delay cost savings. The investor on the upper level affects the decisions of the administrator on lower level, but does not fully control them. The administrator makes their decision autonomously based on the scope of the decision of the upper level.


### 2.3. Environmental Costs

In recent years, more attention has been paid to environmental problems as these have begun to seriously affect both local communities and the economy. Thus, it is essential to consider the environmental costs in the LSCP. In addition, as environmental costs affect the overall project cost, it is necessary to effectively record and calculate environmental costs in the LSCP. Generally, many environmental costs are not usually tracked systematically or attributed to the related processes and outputs but are simply summed and added to total cost [25]. The fact that environmental costs are not fully recorded often leads to distorted calculations [25]. Activity based costing (ABC) is an effective method to record and calculate environmental costs [26]. Cooper provided a comprehensive discussion of ABC [27–30], following the pioneering work of [31, 32]. This method treats activities as accounting objects, and identifies and measures the amount of activities using cost drivers.The environmental costs for retrofitting LSCP transportation networks based on the ABC are as shown in Figure1. As can be seen, the environmental costs are fully recorded according to the environmental cost categories proposed in [25]. Note that some environmental cost categories are related to the work processes, and others are not. Then, all the outputs involved in retrofitting the LSCP transportation network are determined and an analysis of the activity processes and a definition of the activities are prepared to determine the environmental costs. Every activity corresponds to an activity cost center. Jasch [25] proposed that the costs should be more precisely allocated to cost centers. Therefore, environmental costs are either directly allocated to each activity cost center or systematically traced to the responsible environmental media. Of course, if those costs are attributed to outputs directly (i.e., they are not related to the work processes), it is not necessary to allocate them to activity cost centers. Then, the cost drivers for the activity cost centers are determined and the cost driver amounts measured to calculate the cost driver rates. Finally, the environmental costs of each output are determined.Figure 1
Obtain environmental costs for retrofitting transportation network in LSCP based on ABC.Therefore, by using a complete recording method, distortion in the environmental costs can be avoided and through this precise allocation it is easier to effectively manage these costs as it is possible to systematically trace them to the related processes and outputs.


### 2.4. Fuzzy Random Seismic Damage Scenario

To better understand the concepts for the fuzzy random seismic damage scenario, this subsection gives some basic knowledge for the definition and properties of the fuzzy random variables. After Zadeh [33] proposed the concept of fuzzy sets, many scholars have usually tied fuzziness to randomness as possible random outcomes have to be described using fuzzy sets. To describe this fuzziness and randomness, Kwakernaak [7] proposed the concept of fuzzy random variables in 1978. Kruse and Meyer [17] then worked on an expanded version of a similar model. In addition, Puri and Ralescu [14] and Klement et al. [15] also defined fuzzy random variables from other angles. In this paper, the fuzzy random variables are defined in the real number set. This makes the above definitions equivalent [34]. Here, the definition proposed by [14] is utilized.In the following,ℝ is denoted as the set of all real numbers, ℱ
c
(
ℝ
) is denoted as the set of all fuzzy variables, and 𝒦
c
(
ℝ
) is denoted as all of the nonempty bounded close intervals.Definition 1 (see [14]).
In a given probability space(
Ω
,
𝒜
,
Pr
), a mapping ξ
¯
~
:
Ω
→
ℱ
c
(
ℝ
) is called a fuzzy random variable in (
Ω
,
𝒜
,
Pr
); if α
∈
(
0,1
], the set-valued function ξ
¯
~
α
:
Ω
→
𝒦
c
(
ℝ
),
(2)
ξ
¯
~
α
(
ω
)
=
(
ξ
¯
~
(
ω
)
)
α
=
{
x
∣
x
∈
ℝ
,
μ
ξ
¯
~
(
ω
)
(
x
)
≥
α
}
,
∀
ω
∈
Ω
,
is ℱ measurable.Definition 2 (see [35]).
Ifξ
¯
~
1
,
ξ
¯
~
2
,
…
,
ξ
¯
~
n are fuzzy random variables defined in the probability space on (
Ω
,
𝒜
,
Pr
), then ξ
¯
~
=
(
ξ
¯
~
1
,
ξ
¯
~
2
,
…
,
ξ
¯
~
n
) is called fuzzy random vector.Lemma 3 (see [36]).
Letξ
¯
~
=
(
ξ
¯
~
1
,
ξ
¯
~
2
,
…
,
ξ
¯
~
n
) be a fuzzy random vector, and let f be a continuous function from ℝ
𝕞 to ℝ. Then f
(
ξ
¯
~
) is a fuzzy random variable.Definition 4 (see [14]).
In a given probability space(
Ω
,
𝒜
,
Pr
), if ω
∈
Ω,  α
∈
[
0,1
], the mapping ω
↦
(
ξ
¯
~
)
α
-
(
ω
) and ω
↦
(
ξ
¯
~
)
α
+
(
ω
) are integrable; then ξ
¯
~ is called the integrated bounded fuzzy random variable on the probability space (
Ω
,
𝒜
,
Pr
).Definition 5 (see [14]).
Letξ
¯
~ be an integrated bounded fuzzy random variable on the probability space (
Ω
,
𝒜
,
Pr
); the expected value E
(
ξ
¯
~
) of ξ
¯
~ is defined as the only fuzzy set in ℝ; for all α
∈
[
0,1
], it satisfies
(3)
E
(
ξ
¯
~
)
α
=
∫
Ω
‍
ξ
¯
~
α
d
p
=
{
∫
Ω
‍
f
(
ω
)
d
p
(
ω
)
:
f
∈
L
1
(
P
)
,
f
(
ω
)
∈
ξ
¯
~
α
(
ω
)
a
.
s
.
[
P
]
{
∫
Ω
‍
f
(
ω
)
d
p
(
ω
)
:
f
∈
L
1
(
P
)
,
}
,
where ∫
Ω
‍
ξ
¯
~
α
d
p is the Aumann integral of ξ
¯
~
α about P and L
1
(
P
) denote all of the integrable function f
:
Ω
→
ℝ about the probability measure P.Lemma 6 (see [37]).
Let(
Ω
,
𝒜
,
Pr
) be complete probability space; ξ
¯
~
:
Ω
→
ℱ
c
(
ℝ
) is an integrated bounded fuzzy random variable. Then for all α
∈
[
0,1
], the α-set of E
(
ξ
¯
~
) is the compact convex interval as follows:
(4)
E
(
ξ
¯
~
)
α
=
[
(
E
(
ξ
¯
~
)
)
α
-
,
(
E
(
ξ
¯
~
)
)
α
+
]
=
[
∫
Ω
‍
(
ξ
¯
~
(
ω
)
)
α
-
d
p
(
ω
)
,
(
ξ
¯
~
(
ω
)
)
α
+
d
p
(
ω
)
]
.Lemma 7.
Let(
Ω
,
𝒜
,
Pr
) be complete probability space; ξ
¯
~
1, ξ
¯
~
2 are integrated bounded fuzzy random variables on (
Ω
,
𝒜
,
Pr
), λ
,
γ
∈
ℝ, and then
(5)
E
(
λ
ξ
¯
~
1
+
γ
ξ
¯
~
2
)
=
λ
E
(
ξ
¯
~
1
)
+
γ
E
(
ξ
¯
~
2
)
.For the fuzzy random seismic damage scenario, according to [1], advanced structural analysis can lead to a probabilistic assessment of the structural damage for a given earthquake, in terms of a set of discrete probabilities associated with each of the five damage categories. Seismologists also have made predictions as to the probabilities of various earthquake occurrences. These two sets of probabilistic estimations from earthquake-structural engineers and seismologists can be combined to prepare the damage prediction [1]. For the convenience of discussion, seismic damage to a structure (i.e., LSCP transportation network) is usually classified into five categories ranging from no damage to complete collapse. However, a description of the perception result for seismic damage is a category which is vague. In this case, a vague perception of a crisp but unobservable random variable is used as in the following:
(6)
ξ
¯
~
=
(
a
i
L
,
a
i
C
,
a
i
R
)
with
probability
p
i
,
i
=
1
,
…
,
5
.Therefore, the seismic damage scenario can be viewed as a fuzzy random variable, which has a similar sense to the minor automobile collision damage outlined in [10]. See Figure 2 for a detailed description.Figure 2
Fuzzy random seismic damage scenario.An example can be used to explain how to use the fuzzy random variable and to describe the uncertainty in a seismic damage scenario. Suppose that there is a linka
∈
A in an LSCP transportation network. Seismic damage perception has five categories 1, 2, 3, 4, and 5 ranging from no damage to complete collapse and seismic damage randomly emerges with a certain probability. On the other hand, the description of the perception result is vague with values such as “about 1” and “about 3.” These denote the fuzzy sets and can be conveniently described using triangular fuzzy sets, as shown in Figure 3. Here it is assumed that the probabilities for the five categories are 0.1, 0.2, 0.3, 0.3, and 0.1, so the seismic damage scenario can be seen as a fuzzy random variable as in (7) and as shown in Figure 3. It should be mentioned that one damage scenario has different meanings for the different damage ranks (i.e., its membership is different for each different damage rank). Similar examples can be found in [10]. Consider
(7)
ξ
¯
~
a
=
{
(
0,1
,
2
)
with
probability
0.1
(
1,2
,
3
)
with
probability
0.2
(
2,3
,
4
)
with
probability
0.3
(
3,4
,
5
)
with
probability
0.3
(
4,5
,
6
)
with
probability
0.1
.Figure 3
Value of seismic damage scenario.It should be noted that the same category may have different possibilities for different links.


### 2.5. Model Formulation

Denote a transportation network asG
(
B
,
A
), where B is the set of nodes and A is the set of network links. The decision variable on the upper level is u
a
∈
{
0,1
,
2,3
,
4,5
}, which means that a decision has been done for link a to be retrofitted at rank u
a, a
∈
A. For each commodity k
∈
{
1
,
…
,
K
}, x
k
∈
R
+ is the commodity flow (i.e., the decision variable on the lower level), and c
a
b
∈
R
+ is the capacity of node b. Denote f
l
a as the total flow on link a (i.e., f
l
a
=
M
x, for all a
∈
A). To model the retrofit decision with seismic risk in this paper, the assumptions are as follows.(1)
The LSCP transportation network is composed of internal road systems and external road systems and has two types of links, permanent and temporary. In addition, two types of links are designated as critical or noncritical links.
(2)
The links under retrofit consideration are those which are either permanent or critical.
(3)
The retrofit activity process is the same for both permanent and temporary links.
(4)
In the retrofit, the variable environmental and reconstruction costs for the temporary links are considered of less importance than the permanent links.
(5)
The retrofit costs and the retrofit decision have a linear relationship which can be easily relaxed without changing the structure of the proposed model, as long as the data are available to support a more detailed analysis.
(6)
The variable environmental costs and the retrofit decision have a linear relationship which can be easily relaxed without changing the structure of the proposed model, as long as the data are available to support a more detailed analysis.
(7)
Origin-destination pairs (i.e., commodities) are determined in advance.
(8)
Traffic flow can be controlled to achieve system equilibrium [1].
(9)
The preretrofit link damage state is defined as the seismic damage scenario minus the retrofit decisionu
a developed from [1].
(10)
Reconstruction costs have a linear function with the postretrofit damage state which can be easily relaxed without changing the structure of the proposed model.


#### 2.5.1. Upper-Level Programming

The investor on the upper level makes a decision as to whether there should be a retrofit for each linka in the transportation network and what rank the retrofit should be. The decision needs to fully consider the link types (i.e., permanent and temporary, critical and noncritical).Objective Functions. One objective on the upper level is to minimize the retrofit costs, which include the environmental costs. In this paper, from a systems view, the retrofit costs are added directly to the objective function, which differs from [1]. The investor aims to minimize costs through their decision. Based on this assumption, the retrofit costs can be calculated using the sum of the variable and fixed costs for all links in the network. The retrofit costs for the temporary links are lower than the permanent links. In order to distinguish link types, 0-1 variables are introduced. m
a with 1 indicates a permanent link and is 0 otherwise; n
a with 1 indicates a critical link and is 0 otherwise. Therefore, the retrofit costs can be denoted as ∑
a
∈
A
(
m
a
∨
n
a
)
(
(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
+
(
c
f
i
t
+
m
a
c
f
i
p
)
). Here, ∨ is defined as max, namely, max
[
m
a
,
n
a
].Based on the ABC described above and the assumptions,∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a is the variable environment cost v and c
e
f
f denotes the fixed environmental cost f. Through an analysis of the activity processes, the activity definitions, and the cost allocations, c
e
j
c
=
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a is the variable environmental costs of activity cost center j and ∑
f
∈
F
p
e
i
f
f
c
e
f is the fixed environmental costs of output i. After determining the cost drivers for the activity cost centers and measuring the cost driver amounts, r
a
j
=
c
e
j
c
/
a
m
j is the cost driver rate for activity cost center j. The variable environmental costs for output i are ∑
j
∈
J
r
a
j
a
m
i
j. After this, the environmental costs can then be presented as ∑
i
∈
I
(
∑
j
∈
J
(
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
/
a
m
j
)
a
m
i
j
+
∑
f
∈
F
p
e
i
f
f
c
e
f
f
). Thus, the objective can be described as
(8)
C
(
u
)
=
∑
a
∈
A
‍
(
m
a
∨
n
a
)
×
(

(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
+
(
c
f
i
t
+
m
a
c
f
i
p
)

)
+
ρ
∑
i
∈
I
‍
(

∑
j
∈
J
‍

∑
v
∈
V
p
e
j
v
v

∑
a
∈
A

(
m
a
∨
n
a
)

(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
(
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
a
m
i
j
)
.
Here, ρ denotes the weight of the environmental costs and is determined by the investor.The maximization of the retrofit benefit is another upper-level objective. The decision result of the lower-level is denotedQ
(
x
,
ξ
¯
~
) and quantified as savings in reconstruction and travel delay costs. u is the vector of u
a
,
a
∈
A, and ξ
¯
~ is the vector of ξ
¯
~
a
,
a
∈
A. This can be described as
(9)
Q
(
x
,
ξ
¯
~
)
.Here, maximizing the retrofit benefit while minimizing reconstruction costs and travel time delay is denoted asQ
(
u
,
ξ
¯
~
), which will be described in detail in the objective for the lower level.Logical Constraints. To describe the discrete decision variables for practical sense, the constraints in the following are presented:
(10)
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
.The objective functions and constraints above make up the upper-level programming with lower-level programming as in the following:(11)
min
(
C
(
u
)
,
Q
(
x
,
ξ
¯
~
)
)
=
(
×
a
m
i
j
e
f
e
f
r
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
¯
~
)

)

∑
a
∈
A
‍
(
m
a
∨
n
a
)
(
+
(
c
f
i
t
+
m
a
c
f
i
p
)
)

(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
f
i
t
+
m
a
c
f
i
p
)
)
+
ρ
∑
i
∈
I
‍
(
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A

(
m
a
∨
n
a
)

(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
¯
~
)

)
s
.
t
.
{
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
,
lower-level
programming
.


#### 2.5.2. Lower-Level Programming

The administrator on the lower level decides on the commodity flowx
k. In transportation network literature, the flow between each origin-destination pair is often considered as one commodity. Different commodities represent travel between different origin-destination pairs. x
k is used to express the flow of k commodity. This optimal commodity flow decision seeks to achieve optimal retrofit benefits under a postretrofit state once an earthquake event has occurred and seismic damage sustained. First, it is necessary to introduce the postretrofit damage state before describing in detail the lower-level programming.Postretrofit Damage State. A fuzzy random vector Ξ
¯
~ is introduced to describe the damage to the link once an earthquake has occurred after the retrofit, which has been developed from [1]. Here, Ξ
¯
~ is the vector for Ξ
¯
~
a
a
∈
A. Assume that if a link is retrofitted at any rank, its damaged state (i.e., postretrofit damaged state in the earthquake) is denoted as a seismic damage scenario (i.e., preretrofit link damage state) minus the retrofit rank. Here, for demonstration, a negative postretrofit damaged state is not considered, so the negative state is treated as 0 indicating that the link will be intact. The relationship between the preretrofit link damage state ξ
¯
~
a, the retrofit decision u
a, and the postretrofit damage state Ξ
¯
~
a
(
ξ
¯
~
a
,
u
) is described as in the following:
(12)
Ξ
¯
~
a
(
ξ
¯
~
a
,
u
a
)
=
[
ξ
¯
~
a
-
u
a
]
+
,
∀
a
∈
A
.
For any scenario, the postretrofit damaged state of link a can describe the current damaged state of link a in an earthquake after a retrofit. Based on the above, the discussion for the lower-level programming is as follows.Objective Function. Retrofit benefits are the objective of the administrator. They are only quantified as savings in the minimization of reconstruction and travel delay costs [1]. To maximize benefits is to minimize costs. According to this assumption, the reconstruction costs can be presented as ∑
a
∈
A
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
). This is calculated using the sum of the variable and fixed costs for all the links in the network when links are damaged in an earthquake and need to be reconstructed. Travel delay costs are the total travel time of all the links in the network. The travel time of each link is the product of link travel time and link flow. The link travel time depends on the link flow. Their relationship is usually described using a nondecreasing function such as the bureau of public roads (BPR) function [1]. The BPR function is in the form of t
i
a
0
(
1
+
α
(
f
l
a
/
c
a
a
′
)
β
). Where t
i
a
0 and f
l
a are free flow travel time and flow for link a, respectively, c
a
a
′ is the “practical capacity” of link a and is set to be 90% of the design capacity. Thus, the travel delay costs of a can be denoted as t
i
a
0
(
1
+
α
(
f
l
a
/
c
a
a
′
)
β
)
f
l
a. The objective function is presented as shown below:
(13)
Q
(
x
,
ξ
¯
~
)
=
∑
a
∈
A
‍
(
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)


(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
iiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)
,
where γ is a weight coefficient converting the time to a monetary value, α, β are coefficients of the BPR function, and Ξ
¯
~
a is as [
ξ
¯
~
a
-
u
a
]
+.Node Capacity Constraint. Logistics in large scale postdisaster relief is very important [38]. Therefore, once an earthquake event occurs, a working transportation network for disaster relief and the LSCP are critically important, so the nodes in the network should be fully functioning. Therefore, the node capacity should be at capacity. The constraint is to keep transport in accordance with the flow and the capacity of node b as shown in the following:
(14)
W
x
=
c
a
b
,
∀
b
∈
B
,
where W represents the node-commodity adjacency matrix. x is the commodity flow vector for x
k
,
k
∈
K. c
a
b
∈
R
+ is the capacity of node b.Flow Equation Constraint. The total flow on each link a is the sum of all flows of all commodity k that contains a and is obtained using the link commodity adjacency matrix and the commodity flow vector x as in the following:
(15)
f
l
a
=
M
x
,
∀
a
∈
A
.Damaged Link Flow Constraint. This constraint restricts the link flow when a link is damaged by the earthquake as in (16). This constraint is applied to the postretrofit damaged state and the “practical capacity” of link a, which is set at 90% of the design capacity:
(16)
f
l
a
≤
(
1
-
Ξ
¯
~
a
5
)
c
a
a
′
,
∀
a
∈
A
,
where f
l
a is obtained in (15).Logical Constraints.In order to describe the nonnegative variables in the model, the constraints in the following are presented:
(17)
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.The objective function and constraints above compose the lower-level programming as in the following:(18)
Q
(
x
,
ξ
¯
~
)
:
=
min
∑
a
∈
A
(
∑
a
∈
A
(
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
‍

(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
‍
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
∑
a
∈
A
(
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
‍
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
)
s
.
t
.
{
W
x
=
c
a
b
,
∀
b
∈
B
,
f
l
a
=
M
x
,
∀
a
∈
A
,
f
l
a
≤
(
1
-
Ξ
¯
~
a
5
)
c
a
a
′
,
∀
a
∈
A
,
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.


#### 2.5.3. Fuzzy Random Multiobjective Bilevel Programming Model

The complete multiobjective bilevel programming model under a fuzzy random environment is formulated based on the previous discussion as in the following model:(19)
min
(
C
(
u
)
,
Q
(
x
,
ξ
¯
~
)
)
=
(
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
¯
~
)
)
(
x
,
ξ
¯
~
)
∑
a
∈
A
‍
(
m
a
∨
n
a
)
(
(
c
f
i
t
+
m
a
c
f
i
p
)


(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
f
i
t
+
m
a
c
f
i
p
)

)
+
ρ
∑
i
∈
I
(
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
¯
~
)
)
s
.
t
.
{
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
,
Q
(
x
,
ξ
¯
~
)
:
=
min
∑
a
∈
A
(
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
(
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
)


(
m
a
∨
n
a
)
(

(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)

)
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
(
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
)
iiiiiiiiiiiiiiiiiii
s
.
t
.
{
W
x
=
c
a
b
,
iiiii
∀
b
∈
B
,
f
l
a
=
M
x
,
iiiii
∀
a
∈
A
,
f
l
a
≤
(
1
-
Ξ
¯
~
a
5
)
c
a
a
′
,
iiiiiiiiiiiiiiiiiiii
∀
a
∈
A
,
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.


### 2.6. Transformation Approach for Fuzzy Random Variables

In this subsection, some basic knowledge for the fuzzy random variables is stated.Definition 8 (see [33]).
Given a domainU, if A
~ is a fuzzy set on U, then for any x
∈
U, see the following:
(20)
μ
A
~
:
U
⟶
[
0,1
]
,
x
⟶
μ
A
~
(
x
)
,
where μ
A
~ is called a membership function of x with respect to A
~ and μ
A
~ denoted the grade to each point in U with a real number in the interval [
0,1
] that represents the grade of membership of x in A
~. A
~ is called a fuzzy set and described as follows:
(21)
A
~
=
{
(
x
,
μ
A
~
(
x
)
)
∣
x
∈
U
}
.Definition 9 (see [33]).
Let there be a domainU. Let A
~ be a fuzzy set which is defined on U. If α is possibility level and 0
≤
α
≤
1, A
~
α consists of all elements whose degrees of membership in A
~ are greater than or equal to α as in the following:
(22)
A
~
α
=
{
x
∈
U
∣
μ
A
~
(
x
)
≥
α
}
;
then A
~
α is called the α-level set of fuzzy set A
~.Definition 10 (see [39]).
LetΘ be a nonempty set, and let P
(
Θ
) be the power set of Θ. For each A
⊆
P
(
Θ
), there is nonnegative number Pos
{
A
}, called its possibility, such that (1)
Pos
(
∅
)
=
0 and Pos
(
Θ
)
=
1,
(2)
Pos
(
⋃
k
‍
A
k
)
=
sup
k
Pos
(
A
k
) for any arbitrary collection {
A
k
} in P
(
Θ
).
The triple(
Θ
,
P
(
Θ
)
,
Pos
) is called a possibility space. The function Pos is referred to as a possibility measure.Definition 11 (see [40]).
A fuzzy variable is defined as a function from the possibility space(
Θ
,
P
(
Θ
)
,
Pos
) to the real number ℝ.Definition 12.
Letε be a discrete random variable defined on a probability space (
Ω
,
𝒜
,
Pr
) with the discrete distribution P
ε
(
x
)
=
P
{
x
=
x
n
},  n
=
1,2
,
…, and let θ be any given probability level and 0
≤
θ
≤
max
P
ε
(
x
). ε
θ consists of all elements whose value of P
ε
(
x
) for ε is greater than or equal to θ as the following:
(23)
ε
θ
=
{
x
∈
ℝ
∣
P
ε
(
x
)
≥
θ
}
;
then ε
θ is called the θ-level set of random variable ε.As stated in Section2.4, the definition proposed by [14] is used in this paper. Although there are many properties and transformation approaches for the fuzzy random variable, to conveniently convert programming with fuzzy random coefficients into crisp values, Xu and Liu [12] proposed a theorem which could transform fuzzy random variables into fuzzy variables similar to trapezoidal fuzzy numbers. In this paper, this theorem and proof are adjusted to a discrete random distribution with fluctuating lower, central, and upper parameters for the fuzzy properties and extended bounds of possibility for the fuzzy variable.Theorem 13.
Let(24)
ξ
¯
~
=
{
(
a
1
L
,
a
1
C
,
a
1
R
)
w
i
t
h
p
r
o
b
a
b
i
l
i
t
y
p
1
⋮
⋮
(
a
i
L
,
a
i
C
,
a
i
R
)
w
i
t
h
p
r
o
b
a
b
i
l
i
t
y
p
i
⋮
⋮
(
a
I
L
,
a
I
C
,
a
I
R
)
w
i
t
h
p
r
o
b
a
b
i
l
i
t
y
p
I
be a fuzzy random variable, which has discrete random distribution with fluctuating lower, central, and upper parameter for fuzzy property. The discrete distribution is P
ψ
(
x
). δ is any given probability level of random variable; η is any given possibility level of fuzzy variable; then the fuzzy random variable can be transformed into a (
δ
,
η
)-level trapezoidal fuzzy variable.Proof.
Let(25)
ξ
¯
~
=
{
(
a
1
L
,
a
1
C
,
a
1
R
)
with
probability
p
1
⋮
⋮
(
a
i
L
,
a
i
C
,
a
i
R
)
with
probability
p
i
⋮
⋮
(
a
I
L
,
a
I
C
,
a
I
R
)
with
probability
p
I
be a fuzzy random variable, which has discrete random distribution with fluctuating lower, central, and upper parameter for fuzzy property. The discrete distribution is P
ψ
(
x
). According to Definition 8, the δ-level sets (or δ-cuts) of the discrete random variable ψ can be denoted as follows:
(26)
ψ
δ
=
[
ψ
δ
L
,
ψ
δ
R
]
=
{
x
∈
ℝ
∣
P
ψ
(
x
)
≥
δ
}
.
Here, ψ
δ
L
=
min
{
x
∈
ℝ
∣
P
ψ
(
x
)
≥
δ
} and ψ
δ
R
=
max
{
x
∈
ℝ
∣
P
ψ
(
x
)
≥
δ
}. The parameter δ
∈
[
0
,
max
P
ψ
(
x
)
] here reflects the optimism degree for decision-maker. These intervals indicate where the range of the data lies at the probability level δ. Note that ψ
δ is crisp set.
LetX
=
{
x
ω
=
ψ
(
ω
)
∈
ℝ
∣
P
ψ
(
ψ
(
ω
)
)
≥
δ
,
ω
∈
Ω
}; it is not hard to prove that X
=
[
ψ
δ
L
,
ψ
δ
R
]
=
ψ
δ; namely, min
X
=
ψ
δ
L and max
X
=
ψ
δ
R. In other words, ψ
δ
L is the minimum value that ψ achieves with probability δ; ψ
δ
R is the maximum value that ψ achieves with probability δ. Therefore, the δ-level fuzzy random variable ξ
¯
~
δ can be defined as
(27)
ξ
¯
~
δ
=
{
ψ
δ
L
=
(
a
(
δ
,
L
)
L
,
a
(
δ
,
C
)
L
,
a
(
δ
,
R
)
L
)
with
probability
p
δ
L
⋮
⋮
ψ
δ
R
=
(
a
(
δ
,
L
)
R
,
a
(
δ
,
C
)
R
,
a
(
δ
,
R
)
R
)
with
probability
p
δ
R
.
It can also be denoted as follows:
(28)
ξ
¯
~
δ
=
{
ξ
~
δ
(
ω
)
=
(
a
(
δ
,
L
)
(
ω
)
,
a
(
δ
,
C
)
(
ω
)
,
a
(
δ
,
R
)
(
ω
)
)
with
probability
p
(
ω
)
∣
x
ω
∈
X
,
ω
∈
Ω
ξ
~
δ
(
ω
)
=
(
a
(
δ
,
L
)
(
ω
)
,
a
(
δ
,
C
)
(
ω
)
,
a
(
δ
,
R
)
(
ω
)
)
}
,
where ξ
~
δ
(
ω
) is a fuzzy variable. The variable ξ
¯
~
δ can be expressed in another form as ξ
¯
~
δ
=
⋃
ω
∈
Ω
‍
ξ
~
δ
(
ω
)
=
ξ
~
δ
(
Ω
); here ξ
~
δ
(
ω
)
(
ω
∈
Ω
) are fuzzy variables. So the fuzzy random variable ξ
¯
~ is transformed into a group of fuzzy variables ξ
~
δ
(
ω
)
(
ω
∈
Ω
), which is denoted as ξ
~
δ
(
Ω
). On the basis of the concept on fuzzy variable η-level sets (or η-cuts). The parameter 0
≤
η
≤
1 let
(29)
ξ
~
(
δ
,
η
)
(
ω
)
=
[
ξ
(
δ
,
η
)
L
(
ω
)
,
ξ
(
δ
,
η
)
R
(
ω
)
]
=
{
x
∈
U
∣
μ
ξ
~
δ
(
ω
)
(
x
)
≥
η
}
;
then the η-level sets (or η cuts) of ξ
~
δ
(
Ω
) are defined as follows:
(30)
ξ
~
(
δ
,
η
)
(
Ω
)
=
{
ξ
~
(
δ
,
η
)
(
ω
)
=
[
ξ
(
δ
,
η
)
L
(
ω
)
,
ξ
(
δ
,
η
)
R
(
ω
)
]
∣
ω
∈
Ω
}
;
here, ξ
(
δ
,
η
)
L
(
ω
)
=
inf
μ
ξ
~
δ
(
ω
)
-
1
(
η
), ξ
(
δ
,
η
)
R
(
ω
)
=
sup
μ
ξ
~
δ
(
ω
)
-
1
(
η
), ω
∈
Ω. Inspired by the fuzzy expected value of fuzzy random variable proposed by [10], it can be got as follows:
(31)
a
(
δ
,
L
)
=
∑
ω
‍
p
(
ω
)
a
(
δ
,
L
)
(
ω
)
,
a
(
δ
,
R
)
=
∑
ω
‍
p
(
ω
)
a
(
δ
,
R
)
(
ω
)
,
ξ
(
δ
,
η
)
L
=
∑
ω
‍
p
(
ω
)
ξ
(
δ
,
η
)
L
(
ω
)
,
ξ
(
δ
,
η
)
R
=
∑
ω
‍
p
(
ω
)
ξ
(
δ
,
η
)
R
(
ω
)
.
Consequently,ξ
¯
~ can be transformed into ξ
~
(
δ
,
η
) by the δ-cuts and η-cuts. See Figure 4.
Where0
≤
η
≤
1 and δ
∈
[
0
,
max
P
ψ
(
x
)
], let a
(
δ
,
L
)
=
[
s
]
L, a
(
δ
,
R
)
=
[
s
]
R, ξ
(
δ
,
η
)
L
=
s
_, and ξ
(
δ
,
η
)
R
=
s
¯; then the fuzzy random variable ξ
¯
~ can be transformed into the (
δ
,
η
)-level trapezoidal fuzzy variable ξ
~
(
δ
,
η
) by the following equation:
(32)
ξ
¯
~
⟶
ξ
~
(
δ
,
η
)
=
(
[
s
]
L
,
s
_
,
s
¯
,
[
s
]
R
)
.
The parametersδ and η both reflect optimism degree of the decision-maker. Thus, the fuzzy random variable ξ
¯
~ is transformed into a fuzzy variable which is a trapezoidal fuzzy number with the membership function μ
ξ
~
(
δ
,
η
)
(
x
). The value of μ
ξ
~
(
δ
,
η
)
(
x
) at x
∈
[
[
s
]
L
,
[
s
]
R
] is considered subjectively to be 1 as below:
(33)
μ
ξ
~
(
δ
,
η
)
(
x
)
=
{
1
,
if
s
_
≤
x
<
s
¯
,
x
-
[
s
]
L
s
_
-
[
m
]
L
,
if
[
s
]
L
≤
x
<
s
_
,
[
s
]
R
-
x
[
s
]
R
-
s
¯
,
if
s
¯
≤
x
<
[
s
]
R
,
0
,
if
x
<
[
s
]
L
,
x
>
[
s
]
R
.
Theorem13 is proved.Figure 4
The transformation process from fuzzy random variableξ
¯
~ to (
δ
,
η
)-level trapezoidal fuzzy variable ξ
¯
~
(
δ
,
η
).Through Theorem13, the fuzzy random seismic damage scenario, namely, ξ
¯
~, can be transformed into (
δ
,
η
)-level trapezoidal fuzzy variables ξ
~
(
δ
,
η
) and model (19) can be transformed into the following fuzzy multiobjective bilevel programming model:(34)
min
(
C
(
u
)
,
Q
(
x
,
ξ
~
(
δ
,
η
)
)
)
=
(
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
~
(
δ
,
η
)
)
)
∑
a
∈
A
‍
(
m
a
∨
n
a
)
(
(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
+
(
c
f
i
t
+
m
a
c
f
i
p
)
)
+
ρ
∑
i
∈
I
‍
(
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
~
(
δ
,
η
)
)
)
s
.
t
.
{
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
,
Q
(
x
,
ξ
~
(
δ
,
η
)
)
:
=
min
∑
a
∈
A
(
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)


(
m
a
∨
n
a
)
(
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)


(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
×
[
(
ξ
~
(
η
,
δ
)
)
a
-
u
a
]
+
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)
iiiiiiiiiiiiiiiiiiiiiiii
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)
iiiiiiiiii
s
.
t
.
{
W
x
=
c
a
b
,
iiiii
∀
b
∈
B
,
f
l
a
=
M
x
,
iiiii
∀
a
∈
A
,
f
l
a
≤
(
1
-
[
(
ξ
~
(
η
,
δ
)
)
a
-
u
a
]
+
5
)
c
a
a
′
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
∀
a
∈
A
,
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.


### 2.7. Approximation Decomposition of Fuzzy Variables

In model (34), ξ
¯
~ are coefficients, which when transformed to ξ
~
(
δ
,
η
) are fuzzy variables and so can be regarded as fuzzy numbers. Thus, an approximation decomposition method for fuzzy multiobjective linear bilevel programming model is introduced. This method is as in Zhang et al. [19] solution for fuzzy multiobjective bilevel programming, but with some further development done on the fuzzy multiobjective multifollower partial cooperative bilevel programming as outlined in [41].Definition 14 (see [19]).
A fuzzy numbera
~ is defined as a fuzzy set on ℝ, whose membership function μ
a
~ satisfies the following conditions. (1)
μ
a
~ is a mapping from ℝ to the closed interval [
0,1
].
(2)
It is normal; that is, there existsx
∈
ℝ such that μ
a
~
(
x
)
=
1.
(3)
For anyλ
∈
(
0,1
], a
λ
=
{
x
;
μ
a
~
(
x
)
≥
λ
} is a closed interval, denoted by [
a
λ
L
,
a
λ
R
].
Letℱ
(
ℝ
) be the set of all fuzzy numbers. By decomposition theorem of fuzzy sets [19], we have
(35)
a
~
=
⋃
λ
∈
[
0,1
]
‍
λ
[
a
λ
L
,
a
λ
R
]
,
for every a
~
∈
ℱ
(
ℝ
).From Theorems 17 and 18 in research of Zhang et al. [19], the optimal solution for the model can be determined by solving the equivalent crisp multiobjective bilevel programming model as shown below:
(36)
min
(
C
(
u
)
,
Q
(
x
,
ξ
λ
L
(
R
)
)
)
=
(
∑
a
∈
A
‍
(
m
a
∨
n
a
)
(
+
(
c
f
i
t
+
m
a
c
f
i
p
)
)

(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
f
i
t
+
m
a
c
f
i
p
)
)
+
ρ
∑
i
∈
I
‍
(
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
iiiiiiiiiiiiiiiiiiiiiiiiiiiiii
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiii
(
Q
(
x
,
ξ
λ
L
)
,
Q
(
x
,
ξ
λ
R
)
)
(
∑
a
∈
A
‍
(
m
a
∨
n
a
)
(
(
c
v
a
t
+
m
a
c
v
a
p
)
u
a


)
s
.
t
.
{
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
,
Q
(
x
,
ξ
λ
L
(
R
)
)
=
(
Q
(
x
,
ξ
λ
L
)
,
Q
(
x
,
ξ
λ
R
)
)
:
=
min
(
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)
∑
a
∈
A
(
(
f
l
a
c
a
a
′
)
β
(
m
a
∨
n
a
)
(
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)

(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
×
[
(
ξ
a
)
λ
L
-
u
a
]
+
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)
iiiiiiiiiiiiiiiiiiiiiiiiiii
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)
,
∑
a
∈
A
(
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)


(
m
a
∨
n
a
)
(
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)

(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
×
[
(
ξ
a
)
λ
R
-
u
a
]
+
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)
iiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)
)
iiiiiiiiiii
s
.
t
.
{
W
x
=
c
a
b
,
iiii
∀
b
∈
B
,
f
l
a
=
M
x
,
iiiii
∀
a
∈
A
,
f
l
a
≤
(
1
-
[
(
ξ
a
)
λ
L
-
u
a
]
+
5
)
c
a
a
′
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
∀
a
∈
A
,
f
l
a
≤
(
1
-
[
(
ξ
a
)
λ
R
-
u
a
]
+
5
)
c
a
a
′
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
∀
a
∈
A
,
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.After transforming the fuzzy random seismic damage scenarioξ
¯
~ into (
δ
,
η
)-level trapezoidal fuzzy variable ξ
~
(
δ
,
η
), an approximation progress of (
δ
,
η
)-level trapezoidal fuzzy variable ξ
~
(
δ
,
η
) is conducted until termination. During the approximation progress iterations, model (36) is solved within a series of λ valued by a decomposition of the interval [
0,1
] into equal subintervals.


## 2.1. LSCP Transportation Network

The LSCP transportation network is composed of an internal road system and an external road system and is always built based on the existing links around the site which are connected with the newly built links according to transportation need. There are two types of links (i.e., permanent links and temporary links) which vary considerably in terms of quality. In addition, according to the different functions, the links are divided into critical links and noncritical links. Critical links are those which have vital transport functions such as bridges and tunnels and which should be preferentially taken into account [1]. The retrofit decisions for the different link types vary. That is, the links being considered for the retrofit are either considered to be permanent or critical. The retrofit and reconstruction costs for the temporary links are lower than those for the permanent links. Further, the retrofit decision is specific with 0 (i.e., no retrofit) and there are several ranks according to the seismic damage scenarios.


## 2.2. Bilevel Decision Framework

In this paper, the seismic hazard retrofit decision for an LSCP transportation network involves two participants (i.e., the investor who pays for the retrofit and the administrator who controls the transportation). Therefore, these two participants are the decision-makers on two levels, both of whom successively make the retrofit decisions. The investor on upper level decides which retrofit rank should be taken for each link within the range and, therefore, the two objectives on this level are the retrofit costs including the environmental costs and the retrofit benefits. The administrator decides on the commodity flow (i.e., the transportation network flow once seismic damage has occurred) on the lower level according to the decision results of the upper level. On this level, the retrofit benefits are the primary objective. In this paper, the retrofit benefits are quantified as reconstruction and travel delay cost savings. The investor on the upper level affects the decisions of the administrator on lower level, but does not fully control them. The administrator makes their decision autonomously based on the scope of the decision of the upper level.


## 2.3. Environmental Costs

In recent years, more attention has been paid to environmental problems as these have begun to seriously affect both local communities and the economy. Thus, it is essential to consider the environmental costs in the LSCP. In addition, as environmental costs affect the overall project cost, it is necessary to effectively record and calculate environmental costs in the LSCP. Generally, many environmental costs are not usually tracked systematically or attributed to the related processes and outputs but are simply summed and added to total cost [25]. The fact that environmental costs are not fully recorded often leads to distorted calculations [25]. Activity based costing (ABC) is an effective method to record and calculate environmental costs [26]. Cooper provided a comprehensive discussion of ABC [27–30], following the pioneering work of [31, 32]. This method treats activities as accounting objects, and identifies and measures the amount of activities using cost drivers.The environmental costs for retrofitting LSCP transportation networks based on the ABC are as shown in Figure1. As can be seen, the environmental costs are fully recorded according to the environmental cost categories proposed in [25]. Note that some environmental cost categories are related to the work processes, and others are not. Then, all the outputs involved in retrofitting the LSCP transportation network are determined and an analysis of the activity processes and a definition of the activities are prepared to determine the environmental costs. Every activity corresponds to an activity cost center. Jasch [25] proposed that the costs should be more precisely allocated to cost centers. Therefore, environmental costs are either directly allocated to each activity cost center or systematically traced to the responsible environmental media. Of course, if those costs are attributed to outputs directly (i.e., they are not related to the work processes), it is not necessary to allocate them to activity cost centers. Then, the cost drivers for the activity cost centers are determined and the cost driver amounts measured to calculate the cost driver rates. Finally, the environmental costs of each output are determined.Figure 1
Obtain environmental costs for retrofitting transportation network in LSCP based on ABC.Therefore, by using a complete recording method, distortion in the environmental costs can be avoided and through this precise allocation it is easier to effectively manage these costs as it is possible to systematically trace them to the related processes and outputs.


## 2.4. Fuzzy Random Seismic Damage Scenario

To better understand the concepts for the fuzzy random seismic damage scenario, this subsection gives some basic knowledge for the definition and properties of the fuzzy random variables. After Zadeh [33] proposed the concept of fuzzy sets, many scholars have usually tied fuzziness to randomness as possible random outcomes have to be described using fuzzy sets. To describe this fuzziness and randomness, Kwakernaak [7] proposed the concept of fuzzy random variables in 1978. Kruse and Meyer [17] then worked on an expanded version of a similar model. In addition, Puri and Ralescu [14] and Klement et al. [15] also defined fuzzy random variables from other angles. In this paper, the fuzzy random variables are defined in the real number set. This makes the above definitions equivalent [34]. Here, the definition proposed by [14] is utilized.In the following,ℝ is denoted as the set of all real numbers, ℱ
c
(
ℝ
) is denoted as the set of all fuzzy variables, and 𝒦
c
(
ℝ
) is denoted as all of the nonempty bounded close intervals.Definition 1 (see [14]).
In a given probability space(
Ω
,
𝒜
,
Pr
), a mapping ξ
¯
~
:
Ω
→
ℱ
c
(
ℝ
) is called a fuzzy random variable in (
Ω
,
𝒜
,
Pr
); if α
∈
(
0,1
], the set-valued function ξ
¯
~
α
:
Ω
→
𝒦
c
(
ℝ
),
(2)
ξ
¯
~
α
(
ω
)
=
(
ξ
¯
~
(
ω
)
)
α
=
{
x
∣
x
∈
ℝ
,
μ
ξ
¯
~
(
ω
)
(
x
)
≥
α
}
,
∀
ω
∈
Ω
,
is ℱ measurable.Definition 2 (see [35]).
Ifξ
¯
~
1
,
ξ
¯
~
2
,
…
,
ξ
¯
~
n are fuzzy random variables defined in the probability space on (
Ω
,
𝒜
,
Pr
), then ξ
¯
~
=
(
ξ
¯
~
1
,
ξ
¯
~
2
,
…
,
ξ
¯
~
n
) is called fuzzy random vector.Lemma 3 (see [36]).
Letξ
¯
~
=
(
ξ
¯
~
1
,
ξ
¯
~
2
,
…
,
ξ
¯
~
n
) be a fuzzy random vector, and let f be a continuous function from ℝ
𝕞 to ℝ. Then f
(
ξ
¯
~
) is a fuzzy random variable.Definition 4 (see [14]).
In a given probability space(
Ω
,
𝒜
,
Pr
), if ω
∈
Ω,  α
∈
[
0,1
], the mapping ω
↦
(
ξ
¯
~
)
α
-
(
ω
) and ω
↦
(
ξ
¯
~
)
α
+
(
ω
) are integrable; then ξ
¯
~ is called the integrated bounded fuzzy random variable on the probability space (
Ω
,
𝒜
,
Pr
).Definition 5 (see [14]).
Letξ
¯
~ be an integrated bounded fuzzy random variable on the probability space (
Ω
,
𝒜
,
Pr
); the expected value E
(
ξ
¯
~
) of ξ
¯
~ is defined as the only fuzzy set in ℝ; for all α
∈
[
0,1
], it satisfies
(3)
E
(
ξ
¯
~
)
α
=
∫
Ω
‍
ξ
¯
~
α
d
p
=
{
∫
Ω
‍
f
(
ω
)
d
p
(
ω
)
:
f
∈
L
1
(
P
)
,
f
(
ω
)
∈
ξ
¯
~
α
(
ω
)
a
.
s
.
[
P
]
{
∫
Ω
‍
f
(
ω
)
d
p
(
ω
)
:
f
∈
L
1
(
P
)
,
}
,
where ∫
Ω
‍
ξ
¯
~
α
d
p is the Aumann integral of ξ
¯
~
α about P and L
1
(
P
) denote all of the integrable function f
:
Ω
→
ℝ about the probability measure P.Lemma 6 (see [37]).
Let(
Ω
,
𝒜
,
Pr
) be complete probability space; ξ
¯
~
:
Ω
→
ℱ
c
(
ℝ
) is an integrated bounded fuzzy random variable. Then for all α
∈
[
0,1
], the α-set of E
(
ξ
¯
~
) is the compact convex interval as follows:
(4)
E
(
ξ
¯
~
)
α
=
[
(
E
(
ξ
¯
~
)
)
α
-
,
(
E
(
ξ
¯
~
)
)
α
+
]
=
[
∫
Ω
‍
(
ξ
¯
~
(
ω
)
)
α
-
d
p
(
ω
)
,
(
ξ
¯
~
(
ω
)
)
α
+
d
p
(
ω
)
]
.Lemma 7.
Let(
Ω
,
𝒜
,
Pr
) be complete probability space; ξ
¯
~
1, ξ
¯
~
2 are integrated bounded fuzzy random variables on (
Ω
,
𝒜
,
Pr
), λ
,
γ
∈
ℝ, and then
(5)
E
(
λ
ξ
¯
~
1
+
γ
ξ
¯
~
2
)
=
λ
E
(
ξ
¯
~
1
)
+
γ
E
(
ξ
¯
~
2
)
.For the fuzzy random seismic damage scenario, according to [1], advanced structural analysis can lead to a probabilistic assessment of the structural damage for a given earthquake, in terms of a set of discrete probabilities associated with each of the five damage categories. Seismologists also have made predictions as to the probabilities of various earthquake occurrences. These two sets of probabilistic estimations from earthquake-structural engineers and seismologists can be combined to prepare the damage prediction [1]. For the convenience of discussion, seismic damage to a structure (i.e., LSCP transportation network) is usually classified into five categories ranging from no damage to complete collapse. However, a description of the perception result for seismic damage is a category which is vague. In this case, a vague perception of a crisp but unobservable random variable is used as in the following:
(6)
ξ
¯
~
=
(
a
i
L
,
a
i
C
,
a
i
R
)
with
probability
p
i
,
i
=
1
,
…
,
5
.Therefore, the seismic damage scenario can be viewed as a fuzzy random variable, which has a similar sense to the minor automobile collision damage outlined in [10]. See Figure 2 for a detailed description.Figure 2
Fuzzy random seismic damage scenario.An example can be used to explain how to use the fuzzy random variable and to describe the uncertainty in a seismic damage scenario. Suppose that there is a linka
∈
A in an LSCP transportation network. Seismic damage perception has five categories 1, 2, 3, 4, and 5 ranging from no damage to complete collapse and seismic damage randomly emerges with a certain probability. On the other hand, the description of the perception result is vague with values such as “about 1” and “about 3.” These denote the fuzzy sets and can be conveniently described using triangular fuzzy sets, as shown in Figure 3. Here it is assumed that the probabilities for the five categories are 0.1, 0.2, 0.3, 0.3, and 0.1, so the seismic damage scenario can be seen as a fuzzy random variable as in (7) and as shown in Figure 3. It should be mentioned that one damage scenario has different meanings for the different damage ranks (i.e., its membership is different for each different damage rank). Similar examples can be found in [10]. Consider
(7)
ξ
¯
~
a
=
{
(
0,1
,
2
)
with
probability
0.1
(
1,2
,
3
)
with
probability
0.2
(
2,3
,
4
)
with
probability
0.3
(
3,4
,
5
)
with
probability
0.3
(
4,5
,
6
)
with
probability
0.1
.Figure 3
Value of seismic damage scenario.It should be noted that the same category may have different possibilities for different links.


## 2.5. Model Formulation

Denote a transportation network asG
(
B
,
A
), where B is the set of nodes and A is the set of network links. The decision variable on the upper level is u
a
∈
{
0,1
,
2,3
,
4,5
}, which means that a decision has been done for link a to be retrofitted at rank u
a, a
∈
A. For each commodity k
∈
{
1
,
…
,
K
}, x
k
∈
R
+ is the commodity flow (i.e., the decision variable on the lower level), and c
a
b
∈
R
+ is the capacity of node b. Denote f
l
a as the total flow on link a (i.e., f
l
a
=
M
x, for all a
∈
A). To model the retrofit decision with seismic risk in this paper, the assumptions are as follows.(1)
The LSCP transportation network is composed of internal road systems and external road systems and has two types of links, permanent and temporary. In addition, two types of links are designated as critical or noncritical links.
(2)
The links under retrofit consideration are those which are either permanent or critical.
(3)
The retrofit activity process is the same for both permanent and temporary links.
(4)
In the retrofit, the variable environmental and reconstruction costs for the temporary links are considered of less importance than the permanent links.
(5)
The retrofit costs and the retrofit decision have a linear relationship which can be easily relaxed without changing the structure of the proposed model, as long as the data are available to support a more detailed analysis.
(6)
The variable environmental costs and the retrofit decision have a linear relationship which can be easily relaxed without changing the structure of the proposed model, as long as the data are available to support a more detailed analysis.
(7)
Origin-destination pairs (i.e., commodities) are determined in advance.
(8)
Traffic flow can be controlled to achieve system equilibrium [1].
(9)
The preretrofit link damage state is defined as the seismic damage scenario minus the retrofit decisionu
a developed from [1].
(10)
Reconstruction costs have a linear function with the postretrofit damage state which can be easily relaxed without changing the structure of the proposed model.


### 2.5.1. Upper-Level Programming

The investor on the upper level makes a decision as to whether there should be a retrofit for each linka in the transportation network and what rank the retrofit should be. The decision needs to fully consider the link types (i.e., permanent and temporary, critical and noncritical).Objective Functions. One objective on the upper level is to minimize the retrofit costs, which include the environmental costs. In this paper, from a systems view, the retrofit costs are added directly to the objective function, which differs from [1]. The investor aims to minimize costs through their decision. Based on this assumption, the retrofit costs can be calculated using the sum of the variable and fixed costs for all links in the network. The retrofit costs for the temporary links are lower than the permanent links. In order to distinguish link types, 0-1 variables are introduced. m
a with 1 indicates a permanent link and is 0 otherwise; n
a with 1 indicates a critical link and is 0 otherwise. Therefore, the retrofit costs can be denoted as ∑
a
∈
A
(
m
a
∨
n
a
)
(
(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
+
(
c
f
i
t
+
m
a
c
f
i
p
)
). Here, ∨ is defined as max, namely, max
[
m
a
,
n
a
].Based on the ABC described above and the assumptions,∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a is the variable environment cost v and c
e
f
f denotes the fixed environmental cost f. Through an analysis of the activity processes, the activity definitions, and the cost allocations, c
e
j
c
=
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a is the variable environmental costs of activity cost center j and ∑
f
∈
F
p
e
i
f
f
c
e
f is the fixed environmental costs of output i. After determining the cost drivers for the activity cost centers and measuring the cost driver amounts, r
a
j
=
c
e
j
c
/
a
m
j is the cost driver rate for activity cost center j. The variable environmental costs for output i are ∑
j
∈
J
r
a
j
a
m
i
j. After this, the environmental costs can then be presented as ∑
i
∈
I
(
∑
j
∈
J
(
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
/
a
m
j
)
a
m
i
j
+
∑
f
∈
F
p
e
i
f
f
c
e
f
f
). Thus, the objective can be described as
(8)
C
(
u
)
=
∑
a
∈
A
‍
(
m
a
∨
n
a
)
×
(

(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
+
(
c
f
i
t
+
m
a
c
f
i
p
)

)
+
ρ
∑
i
∈
I
‍
(

∑
j
∈
J
‍

∑
v
∈
V
p
e
j
v
v

∑
a
∈
A

(
m
a
∨
n
a
)

(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
(
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
a
m
i
j
)
.
Here, ρ denotes the weight of the environmental costs and is determined by the investor.The maximization of the retrofit benefit is another upper-level objective. The decision result of the lower-level is denotedQ
(
x
,
ξ
¯
~
) and quantified as savings in reconstruction and travel delay costs. u is the vector of u
a
,
a
∈
A, and ξ
¯
~ is the vector of ξ
¯
~
a
,
a
∈
A. This can be described as
(9)
Q
(
x
,
ξ
¯
~
)
.Here, maximizing the retrofit benefit while minimizing reconstruction costs and travel time delay is denoted asQ
(
u
,
ξ
¯
~
), which will be described in detail in the objective for the lower level.Logical Constraints. To describe the discrete decision variables for practical sense, the constraints in the following are presented:
(10)
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
.The objective functions and constraints above make up the upper-level programming with lower-level programming as in the following:(11)
min
(
C
(
u
)
,
Q
(
x
,
ξ
¯
~
)
)
=
(
×
a
m
i
j
e
f
e
f
r
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
¯
~
)

)

∑
a
∈
A
‍
(
m
a
∨
n
a
)
(
+
(
c
f
i
t
+
m
a
c
f
i
p
)
)

(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
f
i
t
+
m
a
c
f
i
p
)
)
+
ρ
∑
i
∈
I
‍
(
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A

(
m
a
∨
n
a
)

(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
¯
~
)

)
s
.
t
.
{
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
,
lower-level
programming
.


### 2.5.2. Lower-Level Programming

The administrator on the lower level decides on the commodity flowx
k. In transportation network literature, the flow between each origin-destination pair is often considered as one commodity. Different commodities represent travel between different origin-destination pairs. x
k is used to express the flow of k commodity. This optimal commodity flow decision seeks to achieve optimal retrofit benefits under a postretrofit state once an earthquake event has occurred and seismic damage sustained. First, it is necessary to introduce the postretrofit damage state before describing in detail the lower-level programming.Postretrofit Damage State. A fuzzy random vector Ξ
¯
~ is introduced to describe the damage to the link once an earthquake has occurred after the retrofit, which has been developed from [1]. Here, Ξ
¯
~ is the vector for Ξ
¯
~
a
a
∈
A. Assume that if a link is retrofitted at any rank, its damaged state (i.e., postretrofit damaged state in the earthquake) is denoted as a seismic damage scenario (i.e., preretrofit link damage state) minus the retrofit rank. Here, for demonstration, a negative postretrofit damaged state is not considered, so the negative state is treated as 0 indicating that the link will be intact. The relationship between the preretrofit link damage state ξ
¯
~
a, the retrofit decision u
a, and the postretrofit damage state Ξ
¯
~
a
(
ξ
¯
~
a
,
u
) is described as in the following:
(12)
Ξ
¯
~
a
(
ξ
¯
~
a
,
u
a
)
=
[
ξ
¯
~
a
-
u
a
]
+
,
∀
a
∈
A
.
For any scenario, the postretrofit damaged state of link a can describe the current damaged state of link a in an earthquake after a retrofit. Based on the above, the discussion for the lower-level programming is as follows.Objective Function. Retrofit benefits are the objective of the administrator. They are only quantified as savings in the minimization of reconstruction and travel delay costs [1]. To maximize benefits is to minimize costs. According to this assumption, the reconstruction costs can be presented as ∑
a
∈
A
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
). This is calculated using the sum of the variable and fixed costs for all the links in the network when links are damaged in an earthquake and need to be reconstructed. Travel delay costs are the total travel time of all the links in the network. The travel time of each link is the product of link travel time and link flow. The link travel time depends on the link flow. Their relationship is usually described using a nondecreasing function such as the bureau of public roads (BPR) function [1]. The BPR function is in the form of t
i
a
0
(
1
+
α
(
f
l
a
/
c
a
a
′
)
β
). Where t
i
a
0 and f
l
a are free flow travel time and flow for link a, respectively, c
a
a
′ is the “practical capacity” of link a and is set to be 90% of the design capacity. Thus, the travel delay costs of a can be denoted as t
i
a
0
(
1
+
α
(
f
l
a
/
c
a
a
′
)
β
)
f
l
a. The objective function is presented as shown below:
(13)
Q
(
x
,
ξ
¯
~
)
=
∑
a
∈
A
‍
(
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)


(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
iiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)
,
where γ is a weight coefficient converting the time to a monetary value, α, β are coefficients of the BPR function, and Ξ
¯
~
a is as [
ξ
¯
~
a
-
u
a
]
+.Node Capacity Constraint. Logistics in large scale postdisaster relief is very important [38]. Therefore, once an earthquake event occurs, a working transportation network for disaster relief and the LSCP are critically important, so the nodes in the network should be fully functioning. Therefore, the node capacity should be at capacity. The constraint is to keep transport in accordance with the flow and the capacity of node b as shown in the following:
(14)
W
x
=
c
a
b
,
∀
b
∈
B
,
where W represents the node-commodity adjacency matrix. x is the commodity flow vector for x
k
,
k
∈
K. c
a
b
∈
R
+ is the capacity of node b.Flow Equation Constraint. The total flow on each link a is the sum of all flows of all commodity k that contains a and is obtained using the link commodity adjacency matrix and the commodity flow vector x as in the following:
(15)
f
l
a
=
M
x
,
∀
a
∈
A
.Damaged Link Flow Constraint. This constraint restricts the link flow when a link is damaged by the earthquake as in (16). This constraint is applied to the postretrofit damaged state and the “practical capacity” of link a, which is set at 90% of the design capacity:
(16)
f
l
a
≤
(
1
-
Ξ
¯
~
a
5
)
c
a
a
′
,
∀
a
∈
A
,
where f
l
a is obtained in (15).Logical Constraints.In order to describe the nonnegative variables in the model, the constraints in the following are presented:
(17)
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.The objective function and constraints above compose the lower-level programming as in the following:(18)
Q
(
x
,
ξ
¯
~
)
:
=
min
∑
a
∈
A
(
∑
a
∈
A
(
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
‍

(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
‍
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
∑
a
∈
A
(
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
‍
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
)
s
.
t
.
{
W
x
=
c
a
b
,
∀
b
∈
B
,
f
l
a
=
M
x
,
∀
a
∈
A
,
f
l
a
≤
(
1
-
Ξ
¯
~
a
5
)
c
a
a
′
,
∀
a
∈
A
,
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.


### 2.5.3. Fuzzy Random Multiobjective Bilevel Programming Model

The complete multiobjective bilevel programming model under a fuzzy random environment is formulated based on the previous discussion as in the following model:(19)
min
(
C
(
u
)
,
Q
(
x
,
ξ
¯
~
)
)
=
(
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
¯
~
)
)
(
x
,
ξ
¯
~
)
∑
a
∈
A
‍
(
m
a
∨
n
a
)
(
(
c
f
i
t
+
m
a
c
f
i
p
)


(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
f
i
t
+
m
a
c
f
i
p
)

)
+
ρ
∑
i
∈
I
(
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
¯
~
)
)
s
.
t
.
{
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
,
Q
(
x
,
ξ
¯
~
)
:
=
min
∑
a
∈
A
(
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
(
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
)


(
m
a
∨
n
a
)
(

(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)

)
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
(
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
)
iiiiiiiiiiiiiiiiiii
s
.
t
.
{
W
x
=
c
a
b
,
iiiii
∀
b
∈
B
,
f
l
a
=
M
x
,
iiiii
∀
a
∈
A
,
f
l
a
≤
(
1
-
Ξ
¯
~
a
5
)
c
a
a
′
,
iiiiiiiiiiiiiiiiiiii
∀
a
∈
A
,
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.


## 2.5.1. Upper-Level Programming

The investor on the upper level makes a decision as to whether there should be a retrofit for each linka in the transportation network and what rank the retrofit should be. The decision needs to fully consider the link types (i.e., permanent and temporary, critical and noncritical).Objective Functions. One objective on the upper level is to minimize the retrofit costs, which include the environmental costs. In this paper, from a systems view, the retrofit costs are added directly to the objective function, which differs from [1]. The investor aims to minimize costs through their decision. Based on this assumption, the retrofit costs can be calculated using the sum of the variable and fixed costs for all links in the network. The retrofit costs for the temporary links are lower than the permanent links. In order to distinguish link types, 0-1 variables are introduced. m
a with 1 indicates a permanent link and is 0 otherwise; n
a with 1 indicates a critical link and is 0 otherwise. Therefore, the retrofit costs can be denoted as ∑
a
∈
A
(
m
a
∨
n
a
)
(
(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
+
(
c
f
i
t
+
m
a
c
f
i
p
)
). Here, ∨ is defined as max, namely, max
[
m
a
,
n
a
].Based on the ABC described above and the assumptions,∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a is the variable environment cost v and c
e
f
f denotes the fixed environmental cost f. Through an analysis of the activity processes, the activity definitions, and the cost allocations, c
e
j
c
=
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a is the variable environmental costs of activity cost center j and ∑
f
∈
F
p
e
i
f
f
c
e
f is the fixed environmental costs of output i. After determining the cost drivers for the activity cost centers and measuring the cost driver amounts, r
a
j
=
c
e
j
c
/
a
m
j is the cost driver rate for activity cost center j. The variable environmental costs for output i are ∑
j
∈
J
r
a
j
a
m
i
j. After this, the environmental costs can then be presented as ∑
i
∈
I
(
∑
j
∈
J
(
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
/
a
m
j
)
a
m
i
j
+
∑
f
∈
F
p
e
i
f
f
c
e
f
f
). Thus, the objective can be described as
(8)
C
(
u
)
=
∑
a
∈
A
‍
(
m
a
∨
n
a
)
×
(

(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
+
(
c
f
i
t
+
m
a
c
f
i
p
)

)
+
ρ
∑
i
∈
I
‍
(

∑
j
∈
J
‍

∑
v
∈
V
p
e
j
v
v

∑
a
∈
A

(
m
a
∨
n
a
)

(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
(
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
a
m
i
j
)
.
Here, ρ denotes the weight of the environmental costs and is determined by the investor.The maximization of the retrofit benefit is another upper-level objective. The decision result of the lower-level is denotedQ
(
x
,
ξ
¯
~
) and quantified as savings in reconstruction and travel delay costs. u is the vector of u
a
,
a
∈
A, and ξ
¯
~ is the vector of ξ
¯
~
a
,
a
∈
A. This can be described as
(9)
Q
(
x
,
ξ
¯
~
)
.Here, maximizing the retrofit benefit while minimizing reconstruction costs and travel time delay is denoted asQ
(
u
,
ξ
¯
~
), which will be described in detail in the objective for the lower level.Logical Constraints. To describe the discrete decision variables for practical sense, the constraints in the following are presented:
(10)
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
.The objective functions and constraints above make up the upper-level programming with lower-level programming as in the following:(11)
min
(
C
(
u
)
,
Q
(
x
,
ξ
¯
~
)
)
=
(
×
a
m
i
j
e
f
e
f
r
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
¯
~
)

)

∑
a
∈
A
‍
(
m
a
∨
n
a
)
(
+
(
c
f
i
t
+
m
a
c
f
i
p
)
)

(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
f
i
t
+
m
a
c
f
i
p
)
)
+
ρ
∑
i
∈
I
‍
(
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A

(
m
a
∨
n
a
)

(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
¯
~
)

)
s
.
t
.
{
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
,
lower-level
programming
.


## 2.5.2. Lower-Level Programming

The administrator on the lower level decides on the commodity flowx
k. In transportation network literature, the flow between each origin-destination pair is often considered as one commodity. Different commodities represent travel between different origin-destination pairs. x
k is used to express the flow of k commodity. This optimal commodity flow decision seeks to achieve optimal retrofit benefits under a postretrofit state once an earthquake event has occurred and seismic damage sustained. First, it is necessary to introduce the postretrofit damage state before describing in detail the lower-level programming.Postretrofit Damage State. A fuzzy random vector Ξ
¯
~ is introduced to describe the damage to the link once an earthquake has occurred after the retrofit, which has been developed from [1]. Here, Ξ
¯
~ is the vector for Ξ
¯
~
a
a
∈
A. Assume that if a link is retrofitted at any rank, its damaged state (i.e., postretrofit damaged state in the earthquake) is denoted as a seismic damage scenario (i.e., preretrofit link damage state) minus the retrofit rank. Here, for demonstration, a negative postretrofit damaged state is not considered, so the negative state is treated as 0 indicating that the link will be intact. The relationship between the preretrofit link damage state ξ
¯
~
a, the retrofit decision u
a, and the postretrofit damage state Ξ
¯
~
a
(
ξ
¯
~
a
,
u
) is described as in the following:
(12)
Ξ
¯
~
a
(
ξ
¯
~
a
,
u
a
)
=
[
ξ
¯
~
a
-
u
a
]
+
,
∀
a
∈
A
.
For any scenario, the postretrofit damaged state of link a can describe the current damaged state of link a in an earthquake after a retrofit. Based on the above, the discussion for the lower-level programming is as follows.Objective Function. Retrofit benefits are the objective of the administrator. They are only quantified as savings in the minimization of reconstruction and travel delay costs [1]. To maximize benefits is to minimize costs. According to this assumption, the reconstruction costs can be presented as ∑
a
∈
A
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
). This is calculated using the sum of the variable and fixed costs for all the links in the network when links are damaged in an earthquake and need to be reconstructed. Travel delay costs are the total travel time of all the links in the network. The travel time of each link is the product of link travel time and link flow. The link travel time depends on the link flow. Their relationship is usually described using a nondecreasing function such as the bureau of public roads (BPR) function [1]. The BPR function is in the form of t
i
a
0
(
1
+
α
(
f
l
a
/
c
a
a
′
)
β
). Where t
i
a
0 and f
l
a are free flow travel time and flow for link a, respectively, c
a
a
′ is the “practical capacity” of link a and is set to be 90% of the design capacity. Thus, the travel delay costs of a can be denoted as t
i
a
0
(
1
+
α
(
f
l
a
/
c
a
a
′
)
β
)
f
l
a. The objective function is presented as shown below:
(13)
Q
(
x
,
ξ
¯
~
)
=
∑
a
∈
A
‍
(
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)


(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
iiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)
,
where γ is a weight coefficient converting the time to a monetary value, α, β are coefficients of the BPR function, and Ξ
¯
~
a is as [
ξ
¯
~
a
-
u
a
]
+.Node Capacity Constraint. Logistics in large scale postdisaster relief is very important [38]. Therefore, once an earthquake event occurs, a working transportation network for disaster relief and the LSCP are critically important, so the nodes in the network should be fully functioning. Therefore, the node capacity should be at capacity. The constraint is to keep transport in accordance with the flow and the capacity of node b as shown in the following:
(14)
W
x
=
c
a
b
,
∀
b
∈
B
,
where W represents the node-commodity adjacency matrix. x is the commodity flow vector for x
k
,
k
∈
K. c
a
b
∈
R
+ is the capacity of node b.Flow Equation Constraint. The total flow on each link a is the sum of all flows of all commodity k that contains a and is obtained using the link commodity adjacency matrix and the commodity flow vector x as in the following:
(15)
f
l
a
=
M
x
,
∀
a
∈
A
.Damaged Link Flow Constraint. This constraint restricts the link flow when a link is damaged by the earthquake as in (16). This constraint is applied to the postretrofit damaged state and the “practical capacity” of link a, which is set at 90% of the design capacity:
(16)
f
l
a
≤
(
1
-
Ξ
¯
~
a
5
)
c
a
a
′
,
∀
a
∈
A
,
where f
l
a is obtained in (15).Logical Constraints.In order to describe the nonnegative variables in the model, the constraints in the following are presented:
(17)
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.The objective function and constraints above compose the lower-level programming as in the following:(18)
Q
(
x
,
ξ
¯
~
)
:
=
min
∑
a
∈
A
(
∑
a
∈
A
(
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
‍

(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
‍
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
∑
a
∈
A
(
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
‍
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
)
s
.
t
.
{
W
x
=
c
a
b
,
∀
b
∈
B
,
f
l
a
=
M
x
,
∀
a
∈
A
,
f
l
a
≤
(
1
-
Ξ
¯
~
a
5
)
c
a
a
′
,
∀
a
∈
A
,
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.


## 2.5.3. Fuzzy Random Multiobjective Bilevel Programming Model

The complete multiobjective bilevel programming model under a fuzzy random environment is formulated based on the previous discussion as in the following model:(19)
min
(
C
(
u
)
,
Q
(
x
,
ξ
¯
~
)
)
=
(
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
¯
~
)
)
(
x
,
ξ
¯
~
)
∑
a
∈
A
‍
(
m
a
∨
n
a
)
(
(
c
f
i
t
+
m
a
c
f
i
p
)


(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
f
i
t
+
m
a
c
f
i
p
)

)
+
ρ
∑
i
∈
I
(
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
¯
~
)
)
s
.
t
.
{
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
,
Q
(
x
,
ξ
¯
~
)
:
=
min
∑
a
∈
A
(
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
(
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
)


(
m
a
∨
n
a
)
(

(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)

)
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
(
(
m
a
∨
n
a
)
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
(
(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
Ξ
¯
~
a
)
iiiiiiiiiiiiiiiiiii
s
.
t
.
{
W
x
=
c
a
b
,
iiiii
∀
b
∈
B
,
f
l
a
=
M
x
,
iiiii
∀
a
∈
A
,
f
l
a
≤
(
1
-
Ξ
¯
~
a
5
)
c
a
a
′
,
iiiiiiiiiiiiiiiiiiii
∀
a
∈
A
,
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.


## 2.6. Transformation Approach for Fuzzy Random Variables

In this subsection, some basic knowledge for the fuzzy random variables is stated.Definition 8 (see [33]).
Given a domainU, if A
~ is a fuzzy set on U, then for any x
∈
U, see the following:
(20)
μ
A
~
:
U
⟶
[
0,1
]
,
x
⟶
μ
A
~
(
x
)
,
where μ
A
~ is called a membership function of x with respect to A
~ and μ
A
~ denoted the grade to each point in U with a real number in the interval [
0,1
] that represents the grade of membership of x in A
~. A
~ is called a fuzzy set and described as follows:
(21)
A
~
=
{
(
x
,
μ
A
~
(
x
)
)
∣
x
∈
U
}
.Definition 9 (see [33]).
Let there be a domainU. Let A
~ be a fuzzy set which is defined on U. If α is possibility level and 0
≤
α
≤
1, A
~
α consists of all elements whose degrees of membership in A
~ are greater than or equal to α as in the following:
(22)
A
~
α
=
{
x
∈
U
∣
μ
A
~
(
x
)
≥
α
}
;
then A
~
α is called the α-level set of fuzzy set A
~.Definition 10 (see [39]).
LetΘ be a nonempty set, and let P
(
Θ
) be the power set of Θ. For each A
⊆
P
(
Θ
), there is nonnegative number Pos
{
A
}, called its possibility, such that (1)
Pos
(
∅
)
=
0 and Pos
(
Θ
)
=
1,
(2)
Pos
(
⋃
k
‍
A
k
)
=
sup
k
Pos
(
A
k
) for any arbitrary collection {
A
k
} in P
(
Θ
).
The triple(
Θ
,
P
(
Θ
)
,
Pos
) is called a possibility space. The function Pos is referred to as a possibility measure.Definition 11 (see [40]).
A fuzzy variable is defined as a function from the possibility space(
Θ
,
P
(
Θ
)
,
Pos
) to the real number ℝ.Definition 12.
Letε be a discrete random variable defined on a probability space (
Ω
,
𝒜
,
Pr
) with the discrete distribution P
ε
(
x
)
=
P
{
x
=
x
n
},  n
=
1,2
,
…, and let θ be any given probability level and 0
≤
θ
≤
max
P
ε
(
x
). ε
θ consists of all elements whose value of P
ε
(
x
) for ε is greater than or equal to θ as the following:
(23)
ε
θ
=
{
x
∈
ℝ
∣
P
ε
(
x
)
≥
θ
}
;
then ε
θ is called the θ-level set of random variable ε.As stated in Section2.4, the definition proposed by [14] is used in this paper. Although there are many properties and transformation approaches for the fuzzy random variable, to conveniently convert programming with fuzzy random coefficients into crisp values, Xu and Liu [12] proposed a theorem which could transform fuzzy random variables into fuzzy variables similar to trapezoidal fuzzy numbers. In this paper, this theorem and proof are adjusted to a discrete random distribution with fluctuating lower, central, and upper parameters for the fuzzy properties and extended bounds of possibility for the fuzzy variable.Theorem 13.
Let(24)
ξ
¯
~
=
{
(
a
1
L
,
a
1
C
,
a
1
R
)
w
i
t
h
p
r
o
b
a
b
i
l
i
t
y
p
1
⋮
⋮
(
a
i
L
,
a
i
C
,
a
i
R
)
w
i
t
h
p
r
o
b
a
b
i
l
i
t
y
p
i
⋮
⋮
(
a
I
L
,
a
I
C
,
a
I
R
)
w
i
t
h
p
r
o
b
a
b
i
l
i
t
y
p
I
be a fuzzy random variable, which has discrete random distribution with fluctuating lower, central, and upper parameter for fuzzy property. The discrete distribution is P
ψ
(
x
). δ is any given probability level of random variable; η is any given possibility level of fuzzy variable; then the fuzzy random variable can be transformed into a (
δ
,
η
)-level trapezoidal fuzzy variable.Proof.
Let(25)
ξ
¯
~
=
{
(
a
1
L
,
a
1
C
,
a
1
R
)
with
probability
p
1
⋮
⋮
(
a
i
L
,
a
i
C
,
a
i
R
)
with
probability
p
i
⋮
⋮
(
a
I
L
,
a
I
C
,
a
I
R
)
with
probability
p
I
be a fuzzy random variable, which has discrete random distribution with fluctuating lower, central, and upper parameter for fuzzy property. The discrete distribution is P
ψ
(
x
). According to Definition 8, the δ-level sets (or δ-cuts) of the discrete random variable ψ can be denoted as follows:
(26)
ψ
δ
=
[
ψ
δ
L
,
ψ
δ
R
]
=
{
x
∈
ℝ
∣
P
ψ
(
x
)
≥
δ
}
.
Here, ψ
δ
L
=
min
{
x
∈
ℝ
∣
P
ψ
(
x
)
≥
δ
} and ψ
δ
R
=
max
{
x
∈
ℝ
∣
P
ψ
(
x
)
≥
δ
}. The parameter δ
∈
[
0
,
max
P
ψ
(
x
)
] here reflects the optimism degree for decision-maker. These intervals indicate where the range of the data lies at the probability level δ. Note that ψ
δ is crisp set.
LetX
=
{
x
ω
=
ψ
(
ω
)
∈
ℝ
∣
P
ψ
(
ψ
(
ω
)
)
≥
δ
,
ω
∈
Ω
}; it is not hard to prove that X
=
[
ψ
δ
L
,
ψ
δ
R
]
=
ψ
δ; namely, min
X
=
ψ
δ
L and max
X
=
ψ
δ
R. In other words, ψ
δ
L is the minimum value that ψ achieves with probability δ; ψ
δ
R is the maximum value that ψ achieves with probability δ. Therefore, the δ-level fuzzy random variable ξ
¯
~
δ can be defined as
(27)
ξ
¯
~
δ
=
{
ψ
δ
L
=
(
a
(
δ
,
L
)
L
,
a
(
δ
,
C
)
L
,
a
(
δ
,
R
)
L
)
with
probability
p
δ
L
⋮
⋮
ψ
δ
R
=
(
a
(
δ
,
L
)
R
,
a
(
δ
,
C
)
R
,
a
(
δ
,
R
)
R
)
with
probability
p
δ
R
.
It can also be denoted as follows:
(28)
ξ
¯
~
δ
=
{
ξ
~
δ
(
ω
)
=
(
a
(
δ
,
L
)
(
ω
)
,
a
(
δ
,
C
)
(
ω
)
,
a
(
δ
,
R
)
(
ω
)
)
with
probability
p
(
ω
)
∣
x
ω
∈
X
,
ω
∈
Ω
ξ
~
δ
(
ω
)
=
(
a
(
δ
,
L
)
(
ω
)
,
a
(
δ
,
C
)
(
ω
)
,
a
(
δ
,
R
)
(
ω
)
)
}
,
where ξ
~
δ
(
ω
) is a fuzzy variable. The variable ξ
¯
~
δ can be expressed in another form as ξ
¯
~
δ
=
⋃
ω
∈
Ω
‍
ξ
~
δ
(
ω
)
=
ξ
~
δ
(
Ω
); here ξ
~
δ
(
ω
)
(
ω
∈
Ω
) are fuzzy variables. So the fuzzy random variable ξ
¯
~ is transformed into a group of fuzzy variables ξ
~
δ
(
ω
)
(
ω
∈
Ω
), which is denoted as ξ
~
δ
(
Ω
). On the basis of the concept on fuzzy variable η-level sets (or η-cuts). The parameter 0
≤
η
≤
1 let
(29)
ξ
~
(
δ
,
η
)
(
ω
)
=
[
ξ
(
δ
,
η
)
L
(
ω
)
,
ξ
(
δ
,
η
)
R
(
ω
)
]
=
{
x
∈
U
∣
μ
ξ
~
δ
(
ω
)
(
x
)
≥
η
}
;
then the η-level sets (or η cuts) of ξ
~
δ
(
Ω
) are defined as follows:
(30)
ξ
~
(
δ
,
η
)
(
Ω
)
=
{
ξ
~
(
δ
,
η
)
(
ω
)
=
[
ξ
(
δ
,
η
)
L
(
ω
)
,
ξ
(
δ
,
η
)
R
(
ω
)
]
∣
ω
∈
Ω
}
;
here, ξ
(
δ
,
η
)
L
(
ω
)
=
inf
μ
ξ
~
δ
(
ω
)
-
1
(
η
), ξ
(
δ
,
η
)
R
(
ω
)
=
sup
μ
ξ
~
δ
(
ω
)
-
1
(
η
), ω
∈
Ω. Inspired by the fuzzy expected value of fuzzy random variable proposed by [10], it can be got as follows:
(31)
a
(
δ
,
L
)
=
∑
ω
‍
p
(
ω
)
a
(
δ
,
L
)
(
ω
)
,
a
(
δ
,
R
)
=
∑
ω
‍
p
(
ω
)
a
(
δ
,
R
)
(
ω
)
,
ξ
(
δ
,
η
)
L
=
∑
ω
‍
p
(
ω
)
ξ
(
δ
,
η
)
L
(
ω
)
,
ξ
(
δ
,
η
)
R
=
∑
ω
‍
p
(
ω
)
ξ
(
δ
,
η
)
R
(
ω
)
.
Consequently,ξ
¯
~ can be transformed into ξ
~
(
δ
,
η
) by the δ-cuts and η-cuts. See Figure 4.
Where0
≤
η
≤
1 and δ
∈
[
0
,
max
P
ψ
(
x
)
], let a
(
δ
,
L
)
=
[
s
]
L, a
(
δ
,
R
)
=
[
s
]
R, ξ
(
δ
,
η
)
L
=
s
_, and ξ
(
δ
,
η
)
R
=
s
¯; then the fuzzy random variable ξ
¯
~ can be transformed into the (
δ
,
η
)-level trapezoidal fuzzy variable ξ
~
(
δ
,
η
) by the following equation:
(32)
ξ
¯
~
⟶
ξ
~
(
δ
,
η
)
=
(
[
s
]
L
,
s
_
,
s
¯
,
[
s
]
R
)
.
The parametersδ and η both reflect optimism degree of the decision-maker. Thus, the fuzzy random variable ξ
¯
~ is transformed into a fuzzy variable which is a trapezoidal fuzzy number with the membership function μ
ξ
~
(
δ
,
η
)
(
x
). The value of μ
ξ
~
(
δ
,
η
)
(
x
) at x
∈
[
[
s
]
L
,
[
s
]
R
] is considered subjectively to be 1 as below:
(33)
μ
ξ
~
(
δ
,
η
)
(
x
)
=
{
1
,
if
s
_
≤
x
<
s
¯
,
x
-
[
s
]
L
s
_
-
[
m
]
L
,
if
[
s
]
L
≤
x
<
s
_
,
[
s
]
R
-
x
[
s
]
R
-
s
¯
,
if
s
¯
≤
x
<
[
s
]
R
,
0
,
if
x
<
[
s
]
L
,
x
>
[
s
]
R
.
Theorem13 is proved.Figure 4
The transformation process from fuzzy random variableξ
¯
~ to (
δ
,
η
)-level trapezoidal fuzzy variable ξ
¯
~
(
δ
,
η
).Through Theorem13, the fuzzy random seismic damage scenario, namely, ξ
¯
~, can be transformed into (
δ
,
η
)-level trapezoidal fuzzy variables ξ
~
(
δ
,
η
) and model (19) can be transformed into the following fuzzy multiobjective bilevel programming model:(34)
min
(
C
(
u
)
,
Q
(
x
,
ξ
~
(
δ
,
η
)
)
)
=
(
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
~
(
δ
,
η
)
)
)
∑
a
∈
A
‍
(
m
a
∨
n
a
)
(
(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
+
(
c
f
i
t
+
m
a
c
f
i
p
)
)
+
ρ
∑
i
∈
I
‍
(
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
Q
(
x
,
ξ
~
(
δ
,
η
)
)
)
s
.
t
.
{
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
,
Q
(
x
,
ξ
~
(
δ
,
η
)
)
:
=
min
∑
a
∈
A
(
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)


(
m
a
∨
n
a
)
(
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)


(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
×
[
(
ξ
~
(
η
,
δ
)
)
a
-
u
a
]
+
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)
iiiiiiiiiiiiiiiiiiiiiiii
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)
iiiiiiiiii
s
.
t
.
{
W
x
=
c
a
b
,
iiiii
∀
b
∈
B
,
f
l
a
=
M
x
,
iiiii
∀
a
∈
A
,
f
l
a
≤
(
1
-
[
(
ξ
~
(
η
,
δ
)
)
a
-
u
a
]
+
5
)
c
a
a
′
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
∀
a
∈
A
,
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.


## 2.7. Approximation Decomposition of Fuzzy Variables

In model (34), ξ
¯
~ are coefficients, which when transformed to ξ
~
(
δ
,
η
) are fuzzy variables and so can be regarded as fuzzy numbers. Thus, an approximation decomposition method for fuzzy multiobjective linear bilevel programming model is introduced. This method is as in Zhang et al. [19] solution for fuzzy multiobjective bilevel programming, but with some further development done on the fuzzy multiobjective multifollower partial cooperative bilevel programming as outlined in [41].Definition 14 (see [19]).
A fuzzy numbera
~ is defined as a fuzzy set on ℝ, whose membership function μ
a
~ satisfies the following conditions. (1)
μ
a
~ is a mapping from ℝ to the closed interval [
0,1
].
(2)
It is normal; that is, there existsx
∈
ℝ such that μ
a
~
(
x
)
=
1.
(3)
For anyλ
∈
(
0,1
], a
λ
=
{
x
;
μ
a
~
(
x
)
≥
λ
} is a closed interval, denoted by [
a
λ
L
,
a
λ
R
].
Letℱ
(
ℝ
) be the set of all fuzzy numbers. By decomposition theorem of fuzzy sets [19], we have
(35)
a
~
=
⋃
λ
∈
[
0,1
]
‍
λ
[
a
λ
L
,
a
λ
R
]
,
for every a
~
∈
ℱ
(
ℝ
).From Theorems 17 and 18 in research of Zhang et al. [19], the optimal solution for the model can be determined by solving the equivalent crisp multiobjective bilevel programming model as shown below:
(36)
min
(
C
(
u
)
,
Q
(
x
,
ξ
λ
L
(
R
)
)
)
=
(
∑
a
∈
A
‍
(
m
a
∨
n
a
)
(
+
(
c
f
i
t
+
m
a
c
f
i
p
)
)

(
c
v
a
t
+
m
a
c
v
a
p
)
u
a
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
f
i
t
+
m
a
c
f
i
p
)
)
+
ρ
∑
i
∈
I
‍
(
∑
j
∈
J
‍
∑
v
∈
V
p
e
j
v
v
∑
a
∈
A
(
m
a
∨
n
a
)
(
c
e
v
t
+
c
e
v
p
)
u
a
a
m
j
iiiiiiiiiiiiiiiiiiiiiiiiiiiiii
×
a
m
i
j
+
∑
f
∈
F
‍
p
e
i
f
f
c
e
f
f
)
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiii
(
Q
(
x
,
ξ
λ
L
)
,
Q
(
x
,
ξ
λ
R
)
)
(
∑
a
∈
A
‍
(
m
a
∨
n
a
)
(
(
c
v
a
t
+
m
a
c
v
a
p
)
u
a


)
s
.
t
.
{
u
a
∈
{
0,1
,
2,3
,
4,5
}
,
∀
a
∈
A
,
Q
(
x
,
ξ
λ
L
(
R
)
)
=
(
Q
(
x
,
ξ
λ
L
)
,
Q
(
x
,
ξ
λ
R
)
)
:
=
min
(
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)
∑
a
∈
A
(
(
f
l
a
c
a
a
′
)
β
(
m
a
∨
n
a
)
(
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)

(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
×
[
(
ξ
a
)
λ
L
-
u
a
]
+
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)
iiiiiiiiiiiiiiiiiiiiiiiiiii
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)
,
∑
a
∈
A
(
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)


(
m
a
∨
n
a
)
(
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)

(
c
r
v
a
t
+
m
a
c
r
v
a
p
)
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
×
[
(
ξ
a
)
λ
R
-
u
a
]
+
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
c
r
f
i
t
+
m
a
c
r
f
i
p
)
)
iiiiiiiiiiiiiiiiiiiiiiiiiiiii
+
γ
t
i
a
0
(
1
+
α
(
f
l
a
c
a
a
′
)
β
)
f
l
a
)
)
iiiiiiiiiii
s
.
t
.
{
W
x
=
c
a
b
,
iiii
∀
b
∈
B
,
f
l
a
=
M
x
,
iiiii
∀
a
∈
A
,
f
l
a
≤
(
1
-
[
(
ξ
a
)
λ
L
-
u
a
]
+
5
)
c
a
a
′
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
∀
a
∈
A
,
f
l
a
≤
(
1
-
[
(
ξ
a
)
λ
R
-
u
a
]
+
5
)
c
a
a
′
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
∀
a
∈
A
,
x
k
≥
0
,
∀
k
=
1
,
…
,
K
.After transforming the fuzzy random seismic damage scenarioξ
¯
~ into (
δ
,
η
)-level trapezoidal fuzzy variable ξ
~
(
δ
,
η
), an approximation progress of (
δ
,
η
)-level trapezoidal fuzzy variable ξ
~
(
δ
,
η
) is conducted until termination. During the approximation progress iterations, model (36) is solved within a series of λ valued by a decomposition of the interval [
0,1
] into equal subintervals.


## 3. An Approximation Decomposition-Based Multiobjective AGLNPSO

Bilevel programming problem is NP-hard, which loosely means that it cannot in general be solved with a polynomial time algorithm [42] and it is difficult to find numerical solutions [43]. Many methods have been proposed to solve these problems, such as the branch-and-bound methods [44, 45], the descent method [46], and the penalty function method [47]. In addition, heuristic algorithms [48] and evolutionary computation [49] have also been proposed to obtain a numerical optimal solution or numerical efficient solution. PSO has been adopted for dealing with multiobjective optimization problems and has been found to be very successful, in addition to heuristics [50]. Therefore, with these considerations, the PSO is the approach adopted in this study. An approximation decomposition-based multiobjective AGLNPSO made up of approximation decomposition [19], PAES [20], AGLNPSO [21], and MOPSO [24] is proposed to solve the problem. Of course, this proposed algorithm may not be the best; however, it can assist in obtaining an effective solution, which is demonstrated in the analysis of the case problem. In the future, in order to get even better solutions more effectively, alternative approaches and algorithms (e.g., other exact approaches, (meta)heuristics, evolutionary algorithms, etc.) will be explored for comparison.


### 3.1. Overall Procedure for the Proposed Algorithm

The procedure for the proposed algorithm is presented as follows; see Figure5.Figure 5
Overall procedure of the proposed method.Step 1.
Initialize approximation coefficientl
=
1, which is used to generate cut sets for the fuzzy numbers in model (36) and the error coefficient ϵ.Step 2.
Decompose interval[
0,1
] into 2
l
-
1 equal sub-intervals with (
2
l
-
1
+
1
) nodes λ
i  (
i
=
0
,
…
,
2
l
-
1
), which are arranged in the order 0
=
λ
0
<
λ
1
<
⋯
<
λ
2
l
-
1
=
1.Step 3.
Transform model (34) into a series of models for model (36) with l.Step 4.
Initialize the parameters: swarm_size, iteration_max, the range of velocity and position for the variables, the personal best position acceleration constant, the global best position acceleration constant, the local best position acceleration constant, the near neighbor best acceleration constant, and the inertia weight_max. Then, initialize the velocities and positions of the particle-represented solutions.Step 5.
Check the feasibility and decode the particles.Step 6.
Solve the lower-level programming with the feasible solutions of the upper level to determine the optimal objective value. Calculate the two objectives on the upper level to evaluate every particle.Step 7.
Calculate thep
b
e
s
t,  g
b
e
s
t, l
b
e
s
t, and n
b
e
s
t using the multiobjective method. Restore the Pareto optimal solutions (i.e., the (global) elite individuals), lower-level programming solutions, and objective values of upper level and lower level.Step 8.
Update the inertia weight for each iteration.Step 9.
Update the velocity and position of each particle.Step 10.
Check the multiobjective AGLNPSO termination. If the stopping criterion (i.e., iteration_max) is met, then end the multiobjective AGLNPSO procedure to obtain the optimal solution(
u
,
x
)
2
l and go to Step 10; otherwise, go to Step 5.Step 11.
Check the approximation termination. If a stabilization of the Pareto optimal solution is achieved, then the Pareto optimal solution for the complete multiobjective bilevel programming model under fuzzy random environments is obtained, and it terminates. Otherwise,l
=
l
+
1; go back to Step 3.Here, the set convergence is proposed in this paper to describe the stabilization of the Pareto optimal solution, which is defined asϖ and expressed as follows:
(37)
m
∈
M
,
n
∈
N
,
χ
=
0
.
If
traversal
M
for
any
m
there
is
n
=
m
,
n
∈
N
;
then
χ
=
χ
+
1
,
ϖ
=
χ
|
M
|
.
That is to say, if ϖ
≥
ϵ, the Pareto optimal solution set is stable and the approximation termination is achieved.The details for the multiobjective AGLNPSO are described as follows and the notations used are shown.s:
Particle index,s
=
1
,
…
,
S
τ:
Iteration index,τ
=
1
,
…
,
T
h:
Dimension index,h
=
1
,
…
,
H
u
r:
Uniform random number in the interval[
0,1
]
w
(
τ
):
Inertia weight in theτth iteration
w
max:
Maximum inertia weight value
w
min:
Minimum inertia weight value
ω
s
h
(
τ
):
Velocity of thesth particle at the hth dimension in the τth iteration
θ
s
h
(
τ
):
Position of thesth particle at thehth dimension in theτth  iteration
θ
s
h
0
(
τ
):
Position for temporary and noncritical link of thesth particle at the hth dimension in the τth iteration
ψ
s
h:
Personal best position of thesth particle at the hth dimension
ψ
g
h:
Global best position of thesth particle at thehth dimension
ψ
s
h
L:
Local best position of thesth particle at the hth dimension
ψ
s
h
N:
Near neighbor best position of thesth particle at the hth  dimension
c
p:
Personal best position acceleration constant
c
g:
Global best position acceleration constant
c
l:
Local best position acceleration constant
c
n:
Near neighbor best position acceleration constant
ω
max:
Maximum velocity value
ω
min:
Minimum velocity value
θ
max:
Maximum position value
θ
min:
Minimum position value
Θ
s:
Vector position of thesth particle [
θ
s
1
,
θ
s
2
,
…
,
θ
s
H
]
Ω
s:
Vector velocity of thesth particle [
ω
s
1
,
ω
s
2
,
…
,
ω
s
H
]
R
s:
Thesth set of solutions
c:
The current solution randomly selected one from the nondominated solutions
c
N:
New generated solution.


### 3.2. Solution Representation

In this paper, the particle-represented solution isA dimensions of retrofit rank u
a within [
0,1
,
2,3
,
4,5
] (i.e., a
∈
A) for all links in the LSCP transportation network.


### 3.3. Particle Swarm Initialization

InitializeS particles as a swarm; generate the sth particle with random position Θ
s in the range {
0,1
,
2,3
,
4,5
}. Randomly generate velocity for each particle in the range {
-
5
,
-
4
,
-
3
,
-
2
,
-
1,0
,
1,2
,
3,4
,
5
}. Set the iteration τ
=
1. Set swarm_size S, iteration_max T, personal best position acceleration constant c
p, global best position acceleration constant c
g, local best position acceleration constant c
l, near neighbor best position acceleration constant c
n, inertia weight_max w
max, and inertia weight_min w
min.


### 3.4. Feasibility Checking and Decoding Method

Since the links to be considered for retrofit are either permanent or critical, check and adjust the position of the temporary and noncritical links to 0. Then, the particle-represented solution can be directly decoded into a solution for the problem as shown in Figure6.Figure 6
Transformation from the particle-represented to the problem solution.


### 3.5. Particle Evaluation

Fors
=
1
,
…
,
S, set Θ
s
(
τ
) into the solution R
s that is u in the upper-level programming and put u into the lower-level programming to determine the optimal solution x and the optimal objective Q
(
x
). Calculate another objective for the upper-level programming C
(
u
).


### 3.6. Multiobjective Method

The multiobjective method is made up of the PAES procedure and the test procedure, and selection is introduced to calculatep
b
e
s
t, g
b
e
s
t, l
b
e
s
t, and n
b
e
s
t. This method uses a truncated archive to store the elite individuals (i.e., nondominated solutions), which is used to separate the objective function space into hypercubes, each of which has a score based on its density. Selection of the best is then based on a roulette wheel selection to select the best hypercube first and then uniformly choose a solution. Note that the initialized solution is regarded as the p
b
e
s
t and the nondominated solution of each particle at the 1th iteration. When the iteration updates, the updated solution and the nondominated solutions are used to calculate the p
b
e
s
t by the method. After the p
b
e
s
t has been confirmed at each iteration, the p
b
e
s
t nondominated solutions for all particles are considered with the g
b
e
s
t nondominated solutions (i.e., there is no g
b
e
s
t nondominated solution at initialization) to calculate the g
b
e
s
t by the method. Similar to the g
b
e
s
t, among all the p
b
e
s
t nondominated solutions from K neighbors of the s
th particle and l
b
e
s
t nondominated solutions, set the l
b
e
s
t is also set using this method. For each particle on each dimension, set ψ
s
h
N
=
ψ
o
h
N that maximizes (
Z
(
Θ
s
)
-
Z
(
Ψ
o
)
)
/
(
θ
s
h
-
ψ
o
h
) to get n
b
e
s
t, o
∈
S
∖
s. Here, the maximization process uses the multiobjective method for the calculation of the g
b
e
s
t and l
b
e
s
t above. The details for the PAES procedure, test procedure, and selection procedure are outlined similarly for the p
b
e
s
t, g
b
e
s
t, l
b
e
s
t, and n
b
e
s
t next and in Procedures 1 and 2, where c is the current solution randomly selected from the nondominated solutions. Note that c is randomly selected from the p
b
e
s
t nondominated solutions to calculate the g
b
e
s
t at the 1th iteration.Procedure 1:PAES.
generate a new solutionc
N
if (c dominates c
N)
discardc
N
else if (c
N dominates c)
replace c with c
N and add c
N to the archive
else if (c
N is dominated by any member of the archive)
discard c
N
else if (c
N dominates any member of the archive)
replace it with c
N and add c
N to the archive and all other members
dominated by c
N are discarded
else
apply test procedure to c, c
N and the archive to determine which
becomes the new current solution  and whether to add c
N to the archive
until a termination criterion has been reached, return to the beginningProcedure 2:Test.
if the archive is not full
addc
N to the archive
if (c
N is in a less crowded region of the archive than c)
accept c
N as the new current solution
else
maintain c as the current solution
else if (c
N is in a less crowded region of the archive than any other member on the archive)
addc
N to the archive, and remove a member of the archive from the most crowded region
if (c
N is in a less crowded region of the archive than c)
acceptc
N as the new current solution
else
maintainc as the current solution
else
don’t addc
N to the archiveSelectionStep 1.
Divide 10 by the number of particles in each hypercube to get its score.Step 2.
Apply roulette wheel selection to hypercube according to their scores and select a hypercube.Step 3.
Uniformly choose a member of that hypercube.Therefore, theg
b
e
s
t nondominated solutions at the T
th are the final solutions of the problem.


### 3.7. Inertia Weight Updating

Update the inertia weight for iterationτ using the following equations:
(38)
ω
¯
=
∑
s
=
1
S
∑
h
=
1
H
|
ω
s
h
|
S
·
H
,
ω
*
=
{
(
1
-
1.8
τ
T
)
ω
max
,
0
≤
τ
≤
T
2
,
(
0.2
-
0.2
τ
T
)
ω
max
,
T
2
≤
τ
≤
T
,
Δ
w
=
(
ω
*
-
ω
¯
)
ω
max
(
w
max
-
w
min
)
,
w
=
w
+
Δ
w
,
w
=
w
max
if
w
>
w
max
,
w
=
w
min
if
w
>
w
min
.


### 3.8. Velocity and Position Updating

Update the velocity and the position of eachs
th particle using the following equations:
(39)
ω
s
h
(
τ
+
1
)
=
w
(
τ
)
ω
s
h
(
τ
)
+
c
p
u
r
(
ψ
s
h
-
θ
s
h
(
τ
)
)
+
c
g
u
r
(
ψ
g
h
-
θ
s
h
(
τ
)
)
+
c
l
u
r
(
ψ
s
h
L
-
θ
s
h
(
τ
)
)
,
θ
s
h
(
τ
+
1
)
=
θ
s
h
(
τ
)
+
ω
s
h
(
τ
+
1
)
.
If
θ
s
h
(
τ
+
1
)
>
θ
max
,
then
set
θ
s
h
(
τ
+
1
)
=
θ
max
ω
s
h
(
τ
+
1
)
=
0
.
If
θ
s
h
(
τ
+
1
)
<
θ
min
,
then
set
θ
s
h
(
τ
+
1
)
=
θ
min
ω
s
h
(
τ
+
1
)
=
0
.


## 3.1. Overall Procedure for the Proposed Algorithm

The procedure for the proposed algorithm is presented as follows; see Figure5.Figure 5
Overall procedure of the proposed method.Step 1.
Initialize approximation coefficientl
=
1, which is used to generate cut sets for the fuzzy numbers in model (36) and the error coefficient ϵ.Step 2.
Decompose interval[
0,1
] into 2
l
-
1 equal sub-intervals with (
2
l
-
1
+
1
) nodes λ
i  (
i
=
0
,
…
,
2
l
-
1
), which are arranged in the order 0
=
λ
0
<
λ
1
<
⋯
<
λ
2
l
-
1
=
1.Step 3.
Transform model (34) into a series of models for model (36) with l.Step 4.
Initialize the parameters: swarm_size, iteration_max, the range of velocity and position for the variables, the personal best position acceleration constant, the global best position acceleration constant, the local best position acceleration constant, the near neighbor best acceleration constant, and the inertia weight_max. Then, initialize the velocities and positions of the particle-represented solutions.Step 5.
Check the feasibility and decode the particles.Step 6.
Solve the lower-level programming with the feasible solutions of the upper level to determine the optimal objective value. Calculate the two objectives on the upper level to evaluate every particle.Step 7.
Calculate thep
b
e
s
t,  g
b
e
s
t, l
b
e
s
t, and n
b
e
s
t using the multiobjective method. Restore the Pareto optimal solutions (i.e., the (global) elite individuals), lower-level programming solutions, and objective values of upper level and lower level.Step 8.
Update the inertia weight for each iteration.Step 9.
Update the velocity and position of each particle.Step 10.
Check the multiobjective AGLNPSO termination. If the stopping criterion (i.e., iteration_max) is met, then end the multiobjective AGLNPSO procedure to obtain the optimal solution(
u
,
x
)
2
l and go to Step 10; otherwise, go to Step 5.Step 11.
Check the approximation termination. If a stabilization of the Pareto optimal solution is achieved, then the Pareto optimal solution for the complete multiobjective bilevel programming model under fuzzy random environments is obtained, and it terminates. Otherwise,l
=
l
+
1; go back to Step 3.Here, the set convergence is proposed in this paper to describe the stabilization of the Pareto optimal solution, which is defined asϖ and expressed as follows:
(37)
m
∈
M
,
n
∈
N
,
χ
=
0
.
If
traversal
M
for
any
m
there
is
n
=
m
,
n
∈
N
;
then
χ
=
χ
+
1
,
ϖ
=
χ
|
M
|
.
That is to say, if ϖ
≥
ϵ, the Pareto optimal solution set is stable and the approximation termination is achieved.The details for the multiobjective AGLNPSO are described as follows and the notations used are shown.s:
Particle index,s
=
1
,
…
,
S
τ:
Iteration index,τ
=
1
,
…
,
T
h:
Dimension index,h
=
1
,
…
,
H
u
r:
Uniform random number in the interval[
0,1
]
w
(
τ
):
Inertia weight in theτth iteration
w
max:
Maximum inertia weight value
w
min:
Minimum inertia weight value
ω
s
h
(
τ
):
Velocity of thesth particle at the hth dimension in the τth iteration
θ
s
h
(
τ
):
Position of thesth particle at thehth dimension in theτth  iteration
θ
s
h
0
(
τ
):
Position for temporary and noncritical link of thesth particle at the hth dimension in the τth iteration
ψ
s
h:
Personal best position of thesth particle at the hth dimension
ψ
g
h:
Global best position of thesth particle at thehth dimension
ψ
s
h
L:
Local best position of thesth particle at the hth dimension
ψ
s
h
N:
Near neighbor best position of thesth particle at the hth  dimension
c
p:
Personal best position acceleration constant
c
g:
Global best position acceleration constant
c
l:
Local best position acceleration constant
c
n:
Near neighbor best position acceleration constant
ω
max:
Maximum velocity value
ω
min:
Minimum velocity value
θ
max:
Maximum position value
θ
min:
Minimum position value
Θ
s:
Vector position of thesth particle [
θ
s
1
,
θ
s
2
,
…
,
θ
s
H
]
Ω
s:
Vector velocity of thesth particle [
ω
s
1
,
ω
s
2
,
…
,
ω
s
H
]
R
s:
Thesth set of solutions
c:
The current solution randomly selected one from the nondominated solutions
c
N:
New generated solution.


## 3.2. Solution Representation

In this paper, the particle-represented solution isA dimensions of retrofit rank u
a within [
0,1
,
2,3
,
4,5
] (i.e., a
∈
A) for all links in the LSCP transportation network.


## 3.3. Particle Swarm Initialization

InitializeS particles as a swarm; generate the sth particle with random position Θ
s in the range {
0,1
,
2,3
,
4,5
}. Randomly generate velocity for each particle in the range {
-
5
,
-
4
,
-
3
,
-
2
,
-
1,0
,
1,2
,
3,4
,
5
}. Set the iteration τ
=
1. Set swarm_size S, iteration_max T, personal best position acceleration constant c
p, global best position acceleration constant c
g, local best position acceleration constant c
l, near neighbor best position acceleration constant c
n, inertia weight_max w
max, and inertia weight_min w
min.


## 3.4. Feasibility Checking and Decoding Method

Since the links to be considered for retrofit are either permanent or critical, check and adjust the position of the temporary and noncritical links to 0. Then, the particle-represented solution can be directly decoded into a solution for the problem as shown in Figure6.Figure 6
Transformation from the particle-represented to the problem solution.


## 3.5. Particle Evaluation

Fors
=
1
,
…
,
S, set Θ
s
(
τ
) into the solution R
s that is u in the upper-level programming and put u into the lower-level programming to determine the optimal solution x and the optimal objective Q
(
x
). Calculate another objective for the upper-level programming C
(
u
).


## 3.6. Multiobjective Method

The multiobjective method is made up of the PAES procedure and the test procedure, and selection is introduced to calculatep
b
e
s
t, g
b
e
s
t, l
b
e
s
t, and n
b
e
s
t. This method uses a truncated archive to store the elite individuals (i.e., nondominated solutions), which is used to separate the objective function space into hypercubes, each of which has a score based on its density. Selection of the best is then based on a roulette wheel selection to select the best hypercube first and then uniformly choose a solution. Note that the initialized solution is regarded as the p
b
e
s
t and the nondominated solution of each particle at the 1th iteration. When the iteration updates, the updated solution and the nondominated solutions are used to calculate the p
b
e
s
t by the method. After the p
b
e
s
t has been confirmed at each iteration, the p
b
e
s
t nondominated solutions for all particles are considered with the g
b
e
s
t nondominated solutions (i.e., there is no g
b
e
s
t nondominated solution at initialization) to calculate the g
b
e
s
t by the method. Similar to the g
b
e
s
t, among all the p
b
e
s
t nondominated solutions from K neighbors of the s
th particle and l
b
e
s
t nondominated solutions, set the l
b
e
s
t is also set using this method. For each particle on each dimension, set ψ
s
h
N
=
ψ
o
h
N that maximizes (
Z
(
Θ
s
)
-
Z
(
Ψ
o
)
)
/
(
θ
s
h
-
ψ
o
h
) to get n
b
e
s
t, o
∈
S
∖
s. Here, the maximization process uses the multiobjective method for the calculation of the g
b
e
s
t and l
b
e
s
t above. The details for the PAES procedure, test procedure, and selection procedure are outlined similarly for the p
b
e
s
t, g
b
e
s
t, l
b
e
s
t, and n
b
e
s
t next and in Procedures 1 and 2, where c is the current solution randomly selected from the nondominated solutions. Note that c is randomly selected from the p
b
e
s
t nondominated solutions to calculate the g
b
e
s
t at the 1th iteration.Procedure 1:PAES.
generate a new solutionc
N
if (c dominates c
N)
discardc
N
else if (c
N dominates c)
replace c with c
N and add c
N to the archive
else if (c
N is dominated by any member of the archive)
discard c
N
else if (c
N dominates any member of the archive)
replace it with c
N and add c
N to the archive and all other members
dominated by c
N are discarded
else
apply test procedure to c, c
N and the archive to determine which
becomes the new current solution  and whether to add c
N to the archive
until a termination criterion has been reached, return to the beginningProcedure 2:Test.
if the archive is not full
addc
N to the archive
if (c
N is in a less crowded region of the archive than c)
accept c
N as the new current solution
else
maintain c as the current solution
else if (c
N is in a less crowded region of the archive than any other member on the archive)
addc
N to the archive, and remove a member of the archive from the most crowded region
if (c
N is in a less crowded region of the archive than c)
acceptc
N as the new current solution
else
maintainc as the current solution
else
don’t addc
N to the archiveSelectionStep 1.
Divide 10 by the number of particles in each hypercube to get its score.Step 2.
Apply roulette wheel selection to hypercube according to their scores and select a hypercube.Step 3.
Uniformly choose a member of that hypercube.Therefore, theg
b
e
s
t nondominated solutions at the T
th are the final solutions of the problem.


## 3.7. Inertia Weight Updating

Update the inertia weight for iterationτ using the following equations:
(38)
ω
¯
=
∑
s
=
1
S
∑
h
=
1
H
|
ω
s
h
|
S
·
H
,
ω
*
=
{
(
1
-
1.8
τ
T
)
ω
max
,
0
≤
τ
≤
T
2
,
(
0.2
-
0.2
τ
T
)
ω
max
,
T
2
≤
τ
≤
T
,
Δ
w
=
(
ω
*
-
ω
¯
)
ω
max
(
w
max
-
w
min
)
,
w
=
w
+
Δ
w
,
w
=
w
max
if
w
>
w
max
,
w
=
w
min
if
w
>
w
min
.


## 3.8. Velocity and Position Updating

Update the velocity and the position of eachs
th particle using the following equations:
(39)
ω
s
h
(
τ
+
1
)
=
w
(
τ
)
ω
s
h
(
τ
)
+
c
p
u
r
(
ψ
s
h
-
θ
s
h
(
τ
)
)
+
c
g
u
r
(
ψ
g
h
-
θ
s
h
(
τ
)
)
+
c
l
u
r
(
ψ
s
h
L
-
θ
s
h
(
τ
)
)
,
θ
s
h
(
τ
+
1
)
=
θ
s
h
(
τ
)
+
ω
s
h
(
τ
+
1
)
.
If
θ
s
h
(
τ
+
1
)
>
θ
max
,
then
set
θ
s
h
(
τ
+
1
)
=
θ
max
ω
s
h
(
τ
+
1
)
=
0
.
If
θ
s
h
(
τ
+
1
)
<
θ
min
,
then
set
θ
s
h
(
τ
+
1
)
=
θ
min
ω
s
h
(
τ
+
1
)
=
0
.


## 4. A Case Study

In this section, computational experiments were carried out on a large-scale water conservancy and hydropower construction project. Through the illustrative example on the data set adopted from the case problem, the proposed approach is validated and the efficiency of the algorithm is tested.


### 4.1. Presentation of Case Problem

The XLD hydropower station LSCP is in XLD gorge section of the JS river located in LB county of SC province and YS county of YN province, an area which is earthquake prone. The Yingjiang, Wenchuan, and Panzhihua-Huili earthquakes all seriously affected the local area. Therefore, it is critical that the LSCP risks be controlled, especially in the transportation network. Therefore, the proposed approach is suitable for use on the transportation network at the XLD hydropower LSCP.The transportation network in the project has an internal road network and an external road network. The internal road network is composed of more than 20 major trunk roads and these roads form a solid network located on the left and right banks. There is a temporary traffic bridge upstream and a permanent traffic bridge downstream. The external road network is composed of several secondary roads used for automobiles, which begins at the project dam and terminates at the PED railway station. In order to apply the proposed approach more conveniently, adjacent roads of the same type have been combined and the concrete shapes of the roads have been ignored. A simplified transportation network illustration is shown in Figure7 which distinguishes the permanent and temporary, critical and noncritical nature of each road according to the practical location of the transportation network. The illustration has 24 nodes and 18 and 29 links. There are 12 commodities in total, which represent travel between different origin-destination pairs. The 16 nodes in these commodities have certain capacities (unit: number of vehicles (
n
)) and others have no capacity. Tables 1 and 2 give the detailed data.Table 1
Travel of the commodities.
Commodityk
Travel of commodity
1′
#24→ #23 → #22 → #21 → #20
2′
#20→ #21 → #22 → #23 →#24
3′
#20→ #19 → #7→ #5 → #2
4′
#2→ #5 → #7 → #19 → #20
5′
#6→ #7 → #19 → #20
6′
#20→ #19 → #7 → #6
7′
#20→ #18 → #17 → #16 → #13
8′
#13→ #16 → #17 → #18 →#20
9′
#20→ #18 → #17 → #16 → #14
10′
#14→ #16 → #17 → #18 →#20
11′
#20→ #18 → #17 → #16 → #15
12′
#15→ #16 → #17 → #18 → #20Table 2
Capacity of node.
Nodeb
#24
#23
#22
#21
#20
#19
#7
#5
#2
#6
#18
#17
#16
#13
#14
#15
Capacityc
a
b (n)
51
51
51
51
149
49
49
22
22
25
49
49
49
16
18
15Figure 7
Simplified transportation network illustration in XLD hydropower LSCP.For each link in the transportation network, there is a free flow travel timet
a
0 (unit: hour (
h
)), a “practical capacity” for the link c
a
′ (unit: number of vehicles (
n
)), which is set to be 90% of the design capacity [1], and a fuzzy random seismic damage scenario ξ
¯
~
a. The corresponding data for this case problem are stated in Table 3. It should be noted that the probabilities assigned to the fuzzy numbers for the fuzzy random numbers ξ
¯
~
a in Table 3 were obtained through a statistical analysis of the historical data in the area.Table 3
Free flow travel timet
i
a
0, practical capacity c
a
a
′, and fuzzy random seismic damage scenario ξ
¯
~
a of each link.
Linka
Corresponding nodesb
Free flow travel timet
i
a
0 (h)
Practical capacityc
a
a
′ (n)
Fuzzy random seismic damage scenarioξ
¯
~
a
1
#1, #2
0.10
72
ξ
¯
~
1
=
{
(
0,1
,
2
)
with  probability
13.7
%
(
1,2
,
3
)
with  probability
18.9
%
(
2,3
,
4
)
with  probability
25.8
%
(
3,4
,
5
)
with  probability
16.9
%
(
4,5
,
6
)
with  probability
24.7
%
2
#2, #3
0.08
75
ξ
¯
~
2
=
{
(
0,1
,
2
)
with  probability
7.8
%
(
1,2
,
3
)
with  probability
15.2
%
(
2,3
,
4
)
with  probability
23.4
%
(
3,4
,
5
)
with  probability
29.5
%
(
4,5
,
6
)
with  probability
24.1
%
34567891011
#1, #5#2, #5#2, #6#3, #4#4, #6#6, #7#5, #7#6, #8#7, #19
0.25
87
{
ξ
¯
~
3
ξ
¯
~
4
ξ
¯
~
5
ξ
¯
~
6
ξ
¯
~
7
ξ
¯
~
8
ξ
¯
~
9
ξ
¯
~
10
ξ
¯
~
11
=
{
(
0,1
,
2
)
with  probability
11.5
%
(
1,2
,
3
)
with  probability
8.9
%
(
2,3
,
4
)
with  probability
32.8
%
(
3,4
,
5
)
with  probability
27.4
%
(
4,5
,
6
)
with  probability
19.4
%
12
#19, #20
0.10
101
ξ
¯
~
12
=
{
(
0,1
,
2
)
with  probability
8.6
%
(
1,2
,
3
)
with  probability
20.3
%
(
2,3
,
4
)
with  probability
28.6
%
(
3,4
,
5
)
with  probability
21.5
%
(
4,5
,
6
)
with  probability
21.1
%
13
#10, #11
0.05
96
ξ
¯
~
13
=
{
(
0,1
,
2
)
with  probability
6.5
%
(
1,2
,
3
)
with  probability
17.2
%
(
2,3
,
4
)
with  probability
13.7
%
(
3,4
,
5
)
with  probability
27.2
%
(
4,5
,
6
)
with  probability
35.4
%
14151617181920212325
#1, #5#12, #14#14, #16#9, #14#11, #13#13, #16#10, #13#15, #16#16, #17#18, #20
0.30
90
{
ξ
¯
~
14
ξ
¯
~
15
ξ
¯
~
16
ξ
¯
~
17
ξ
¯
~
18
ξ
¯
~
19
ξ
¯
~
20
ξ
¯
~
21
ξ
¯
~
23
ξ
¯
~
25
=
{
(
0,1
,
2
)
with  probability
12.8
%
(
1,2
,
3
)
with  probability
20.3
%
(
2,3
,
4
)
with  probability
16.5
%
(
3,4
,
5
)
with  probability
31.2
%
(
4,5
,
6
)
with  probability
19.2
%
22
#8, #15
0.09
89
ξ
¯
~
22
=
{
(
0,1
,
2
)
with  probability
7.4
%
(
1,2
,
3
)
with  probability
19.4
%
(
2,3
,
4
)
with  probability
21.6
%
(
3,4
,
5
)
with  probability
23.8
%
(
4,5
,
6
)
with  probability
27.8
%
24
#17, #18
0.10
84
ξ
¯
~
24
=
{
(
0,1
,
2
)
with  probability
23.2
%
(
1,2
,
3
)
with  probability
9.1
%
(
2,3
,
4
)
with  probability
20.8
%
(
3,4
,
5
)
with  probability
21.6
%
(
4,5
,
6
)
with  probability
25.3
%
262729
#20, #21#21, #22#23, #24
0.40
126
{
ξ
¯
~
26
ξ
¯
~
27
ξ
¯
~
29
=
{
(
0,1
,
2
)
with  probability
16.6
%
(
1,2
,
3
)
with  probability
23.2
%
(
2,3
,
4
)
with  probability
7.9
%
(
3,4
,
5
)
with  probability
30.1
%
(
4,5
,
6
)
with  probability
22.2
%
28
#22, #23
0.15
96
ξ
¯
~
28
=
{
(
0,1
,
2
)
with  probability
14.7
%
(
1,2
,
3
)
with  probability
121.3
%
(
2,3
,
4
)
with  probability
10.8
%
(
3,4
,
5
)
with  probability
28.1
%
(
4,5
,
6
)
with  probability
25.1
%In the case problem, there are 2 outputs—the retrofit of permanent (i.e.,i
=
1) and temporary (i.e., i
=
2) links. The activity process and activities for the retrofit are the same for both types of links. Activities (i.e., j
=
1
,
…
,
10) are as (
1
) breaking pavement, (
2
) digging grooves, (
3
) laying pipes, (
4
) backfilling grooves, (
5
) strengthening earth-rock, (
6
) evening roadbeds, (
7
) digging gutters, (
8
) building kerbstones, (
9
) constructing bases, (
10
) constructing pavements. Every activity corresponds to an activity cost center. According to [25], 5 environmental media can be determined in the retrofit work: (
1
) air and climate, (
2
) waste water, (
3
) waste, (
4
) soil and ground water, (
5
) noise and vibration. Environmental cost categories are recorded as in [25].(1)
Waste and emission treatment: depreciation for related equipment; maintenance, operating materials and services; related personnel; fees, taxes, and charges; insurance for environmental liabilities.
(2)
Prevention and environmental management: external services for environmental management and personnel for general environmental management activities.
(3)
Material purchase value of nonproduct output: raw materials; auxiliary materials; operating materials; packaging; energy; and water.
(4)
Processing costs of nonproduct output: depreciation for machinery and labor hours.Since categories 1, 3, and 4 are variable costs (i.e.,v
=
1,2
,
3), category 2 is fixed costs (i.e., f
=
1); they are recoded as (unit:  ¥) [
c
e
1
p
,
c
e
2
p
,
c
e
3
p
]
=
[
11800,5764,2216
], [
c
e
1
t
,
c
e
2
t
,
c
e
3
t
]
=
[
73870,8665,3365
], and c
e
1
f
=
36538. Based on the descriptions above, the environmental costs of the retrofit are allocated as shown in Figure 8.Figure 8
Allocate environmental costs of retrofit work.Based on the practice of the case problem and Figure8, the percentage output for the fixed environmental costs 1 is [
p
e
1
f
,
p
e
2
f
]
=
[
90.8
%
,
9.2
%
]. To calculate the final environmental costs for each output, the corresponding data for the activity cost centers j
=
1
,
…
,
10 is used as in Table 4. The percentage of j in the variable environmental cost categories is denoted as [
p
e
j
1
v
,
p
e
j
2
v
,
p
e
j
3
v
]. a
m
j denotes the cost driver amount in j and [
a
m
1
j
,
a
m
2
j
] denotes the cost driver amount of outputs in j.Table 4
Corresponding data of activity cost centers.
Activity cost centerj
p
e
j
1
v
p
e
j
2
v
p
e
j
3
v
a
m
j
a
m
1
j
a
m
2
j
1
15.6%
0.0%
0.0%
426.30 m2
408.10 m2
18.20 m2
2
6.7%
0.0%
0.0%
9.74 m3
9.32 m3
0.42 m3
3
6.7%
25.0%
25.0%
60.90 m
58.30 m
2.60 m
4
6.7%
0.0%
0.0%
3.35 m3
3.21 m3
0.14 m3
5
5.4%
0.0%
0.0%
25.58 m3
24.49 m3
1.09 m3
6
15.6%
0.0%
0.0%
426.30 m2
408.10 m2
18.20 m2
7
6.7%
0.0%
0.0%
0.61 m3
0.58 m3
0.03 m3
8
5.4%
25.0%
25.0%
60.90 m
58.30 m
2.60 m
9
15.6%
25.0%
25.0%
426.30 m2
408.1 m2
18.20 m2
10
15.6%
25.0%
25.0%
426.30 m2
408.10 m2
18.20 m2In addition, the corresponding cost data for the retrofit and reconstruction are shown in Table5. The values for the other model parameters are as follows: δ
=
0.2, η
=
0.6, ρ
=
1, α
=
0.25, β
=
2, and γ
=
1.Table 5
Cost data of retrofit and reconstruction.
Item
Cost (¥)
Increased variable retrofit costs for permanent link by basic rank (i.e., rank 1)c
v
a
p
16732
Variable retrofit costs for temporary link by basic rank (i.e., rank 1)c
v
a
t
30528
Increased fixed retrofit costs for permanent linkc
f
i
p
14525
Fixed retrofit costs for temporary linkc
f
i
t
28637
Increased variable reconstruction cost for permanent link by basic rank (i.e., rank 1)c
r
v
a
p
107052
Variable reconstruction cost for temporary link by basic rank (i.e., rank 1)c
r
v
a
t
98063
Increased fixed reconstruction cost for permanent linkc
r
f
i
p
69894
Fixed reconstruction cost for temporary linkc
r
f
i
t
50183


### 4.2. Case Solution

The developed algorithm was adopted using MATLAB 7.0 on an Inter Core 2, 2.00 GHz clock pulse with 2048 MB memory. The algorithmic parameters for the case problem were set as follows: error coefficientϵ
=
0.9, swarm_size S
=
20, iteration_max T
=
100, inertia weight_max w
max
=
0.9, inertia weight_min w
min
=
0.1, personal best position acceleration constant c
p
=
0.5, global best position acceleration constant c
g
=
0.5, local best position acceleration constant c
l
=
0.2, and near neighbor best acceleration constant c
n
=
0.1.After 8 iterations of the approximation decomposition, the approximation termination was achieved within 36 minutes on an average of 10 runs, which is time acceptable. The optimal solutions are shown in Tables6 and 7. For all the Pareto optimal solutions on the upper level, the corresponding solutions on the lower level are the same. Table 6 shows a Pareto optimal solution set with 45 solutions on the upper level, where only 10 of the set are enumerated for convenience of expression. The investor is able to choose their preferred plan from the set. If they feel that the retrofit costs including environmental costs C are more important, they would choose the minimum costs plan and vice versa. Since there are fuzzy numbers in model (36), it cannot state crisp optimal objective values in final decision results. However, they are easy to be transformed into equivalent crisp forms by many fuzzy theories and it will not effect the decision.Table 6
Pareto optimal solution set on the upperlevel.
Solutionu
a
1*
2*
3*
4*
5*
6*
7*
8*
9*
10*
……
Linka
1
3
3
3
3
4
4
2
1
3
3
……
2
4
3
4
3
4
4
1
5
2
4
……
3
2
4
2
3
2
1
2
1
5
2
……
4
2
2
3
3
5
1
3
3
1
3
……
5
4
3
3
1
4
4
1
5
3
4
……
6
3
4
3
3
5
2
5
3
3
3
……
7
4
3
3
3
3
5
1
1
2
3
……
8
2
3
2
3
1
1
3
4
3
3
……
9
3
4
3
2
5
1
2
4
4
3
……
10
3
3
3
2
3
3
4
1
3
2
……
11
2
2
2
2
3
3
2
1
1
2
……
12
3
3
4
2
2
3
3
4
2
4
……
13
3
3
3
1
5
3
4
4
2
3
……
14
4
3
4
1
1
5
3
2
3
3
……
15
2
2
1
1
2
1
2
1
3
2
……
16
3
2
3
3
3
3
3
4
1
3
……
17
3
2
3
2
2
4
3
2
2
3
……
18
3
2
3
3
3
4
4
1
1
3
……
19
3
3
2
5
1
2
5
1
3
2
……
20
2
3
2
2
2
2
2
2
4
3
……
21
2
4
2
5
1
1
4
3
5
4
……
22
3
3
3
5
4
1
4
3
4
2
……
23
2
3
2
5
2
1
2
2
3
2
……
24
2
2
3
4
1
2
3
4
2
2
……
25
3
2
3
2
3
3
1
5
2
3
……
26
2
2
2
3
1
1
3
1
1
2
……
27
4
3
4
1
3
4
3
3
2
3
……
28
3
2
3
3
3
5
2
3
1
3
……
29
4
4
4
5
4
5
2
5
3
4
……Table 7
Optimal solutions on the lowerlevel.
Commodityk
1′
2′
3′
4′
5′
6′
7′
8′
9′
10′
11′
12′
Flow of commodityx
k (n/h)
25.50
25.50
12.00
12.00
12.50
12.50
8.00
8.00
9.00
9.00
7.50
7.50


### 4.3. Analytic Results of the Proposed Approach

(
1)  Worthiness of Modeling and Solutions. Fuzzy random programming approach explicitly considers the entire range of uncertain scenarios and thus is more conforming to reality for hedging better against uncertainty. Although it increases the complexity of modeling, the model is well transformed orienting computer implementation. In addition, the computing complexity of the model is not greater than that of the stochastic programming approach used in [1]. Therefore, the extra effort on modeling and solving fuzzy random bilevel programming is worthwhile.The multiobjective method is introduced to determine the Pareto optimal solution set for the upper level and provides more effective and nondominated alternatives for the decision-maker. Compared with the weight-sum method for multiobjective in [19], the solutions in this paper have more reference value for the decision-maker and reflect the users’ preference requirements. Therefore, it is more worthwhile.(
2)  Efficiency of Algorithm. This paper compares GA, an usually used metaheuristics algorithm, and the developed algorithm in this study. The merit of GA is its strong evolutionary process to find an optimal solution by the operation of crossover, mutation, and selection. However, the randomly generated initial generation at the algorithms’ beginning affects the solution quality because of the bad gene inherited from the parent generation. Moreover, the searching capability is reduced as GA does not rely on gradient or derivative information. In the developed algorithm, the particle-represented solutions closely connect the particles of PSO and the solutions to the problem. The hybrid particle-updating mechanism successfully enhances the searching capability. The developed algorithm is a useful tool for the solution to the problem. In contrast to the previous papers, such as [19, 21], both the bilevel and multiobjective environments are considered in this paper.For multiobjective optimization, the definition of quality is substantially more complex than for single-objective optimization problems. Figure9 shows the iterative results of Pareto optimal solutions in 8 iterations. Note that in each iteration of approximation decomposition, these fuzzy numbers in model (36) are decomposed into crisp forms. Therefore, the objective values can be evaluated. For further expression of the efficiency of the convergence, three metrics of performance are studied. Table 8 shows the metrics of performance for Pareto optimal sets proposed in [51] and the set convergence of Pareto optimal sets in 8 iterations after 10 runs.Table 8
Metrics of performance and set convergence for Pareto optimal sets.
Iteration
The average distance
The distribution
The extent
The set convergence
1
0.0812
0.5623
5.0088
—
2
0.0630
0.4615
5.7706
Iteration 1–2: 0.5610
3
0.1364
0.8600
6.8594
Iteration 2–3: 0.6765
4
0.1223
0.5250
6.3147
Iteration 3–4: 0.7632
5
0.1502
0.6820
6.7538
Iteration 4–5: 0.8750
6
0.2835
0.9600
6.6417
Iteration 5–6: 0.8571
7
0.0973
0.5289
5.4451
Iteration 6–7: 0.5712
8
0.1300
0.8635
6.6424
Iteration 7–8: 0.9143Iterative results of Pareto optimal solutions in 8 iterations.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)(
3)  Efficiency of Algorithm Parameters. Since the search space for the problem is so large and the computing process is time-consuming, it is necessary to choose reasonable algorithm parameters. Table 9 shows the comparison results of different swarm_size (i.e., S) and iteration_max (i.e., T) on average computing time and iterations after 10 runs. Looking into Table 9, when S
=
10, the program could not obtain results in more than 3600 s and more than 10 iterations with both T
=
100 and T
=
200. This is the same when S
=
50 with the lower iterations. When S
=
20 and T
=
200, the process is more time-consuming than the current algorithm parameters. Therefore, the current ones are considered suitable.Table 9
Comparison of different algorithm parameters.
S
=
10
S
=
20
S
=
50
T
=
100
T
=
200
T
=
100
T
=
200
T
=
100
T
=
200
Computing time (s)
≥3600
≥3600
2160
3480
≥3600
≥3600
Iterations
≥16
≥12
8
6
≥4
≥3


## 4.1. Presentation of Case Problem

The XLD hydropower station LSCP is in XLD gorge section of the JS river located in LB county of SC province and YS county of YN province, an area which is earthquake prone. The Yingjiang, Wenchuan, and Panzhihua-Huili earthquakes all seriously affected the local area. Therefore, it is critical that the LSCP risks be controlled, especially in the transportation network. Therefore, the proposed approach is suitable for use on the transportation network at the XLD hydropower LSCP.The transportation network in the project has an internal road network and an external road network. The internal road network is composed of more than 20 major trunk roads and these roads form a solid network located on the left and right banks. There is a temporary traffic bridge upstream and a permanent traffic bridge downstream. The external road network is composed of several secondary roads used for automobiles, which begins at the project dam and terminates at the PED railway station. In order to apply the proposed approach more conveniently, adjacent roads of the same type have been combined and the concrete shapes of the roads have been ignored. A simplified transportation network illustration is shown in Figure7 which distinguishes the permanent and temporary, critical and noncritical nature of each road according to the practical location of the transportation network. The illustration has 24 nodes and 18 and 29 links. There are 12 commodities in total, which represent travel between different origin-destination pairs. The 16 nodes in these commodities have certain capacities (unit: number of vehicles (
n
)) and others have no capacity. Tables 1 and 2 give the detailed data.Table 1
Travel of the commodities.
Commodityk
Travel of commodity
1′
#24→ #23 → #22 → #21 → #20
2′
#20→ #21 → #22 → #23 →#24
3′
#20→ #19 → #7→ #5 → #2
4′
#2→ #5 → #7 → #19 → #20
5′
#6→ #7 → #19 → #20
6′
#20→ #19 → #7 → #6
7′
#20→ #18 → #17 → #16 → #13
8′
#13→ #16 → #17 → #18 →#20
9′
#20→ #18 → #17 → #16 → #14
10′
#14→ #16 → #17 → #18 →#20
11′
#20→ #18 → #17 → #16 → #15
12′
#15→ #16 → #17 → #18 → #20Table 2
Capacity of node.
Nodeb
#24
#23
#22
#21
#20
#19
#7
#5
#2
#6
#18
#17
#16
#13
#14
#15
Capacityc
a
b (n)
51
51
51
51
149
49
49
22
22
25
49
49
49
16
18
15Figure 7
Simplified transportation network illustration in XLD hydropower LSCP.For each link in the transportation network, there is a free flow travel timet
a
0 (unit: hour (
h
)), a “practical capacity” for the link c
a
′ (unit: number of vehicles (
n
)), which is set to be 90% of the design capacity [1], and a fuzzy random seismic damage scenario ξ
¯
~
a. The corresponding data for this case problem are stated in Table 3. It should be noted that the probabilities assigned to the fuzzy numbers for the fuzzy random numbers ξ
¯
~
a in Table 3 were obtained through a statistical analysis of the historical data in the area.Table 3
Free flow travel timet
i
a
0, practical capacity c
a
a
′, and fuzzy random seismic damage scenario ξ
¯
~
a of each link.
Linka
Corresponding nodesb
Free flow travel timet
i
a
0 (h)
Practical capacityc
a
a
′ (n)
Fuzzy random seismic damage scenarioξ
¯
~
a
1
#1, #2
0.10
72
ξ
¯
~
1
=
{
(
0,1
,
2
)
with  probability
13.7
%
(
1,2
,
3
)
with  probability
18.9
%
(
2,3
,
4
)
with  probability
25.8
%
(
3,4
,
5
)
with  probability
16.9
%
(
4,5
,
6
)
with  probability
24.7
%
2
#2, #3
0.08
75
ξ
¯
~
2
=
{
(
0,1
,
2
)
with  probability
7.8
%
(
1,2
,
3
)
with  probability
15.2
%
(
2,3
,
4
)
with  probability
23.4
%
(
3,4
,
5
)
with  probability
29.5
%
(
4,5
,
6
)
with  probability
24.1
%
34567891011
#1, #5#2, #5#2, #6#3, #4#4, #6#6, #7#5, #7#6, #8#7, #19
0.25
87
{
ξ
¯
~
3
ξ
¯
~
4
ξ
¯
~
5
ξ
¯
~
6
ξ
¯
~
7
ξ
¯
~
8
ξ
¯
~
9
ξ
¯
~
10
ξ
¯
~
11
=
{
(
0,1
,
2
)
with  probability
11.5
%
(
1,2
,
3
)
with  probability
8.9
%
(
2,3
,
4
)
with  probability
32.8
%
(
3,4
,
5
)
with  probability
27.4
%
(
4,5
,
6
)
with  probability
19.4
%
12
#19, #20
0.10
101
ξ
¯
~
12
=
{
(
0,1
,
2
)
with  probability
8.6
%
(
1,2
,
3
)
with  probability
20.3
%
(
2,3
,
4
)
with  probability
28.6
%
(
3,4
,
5
)
with  probability
21.5
%
(
4,5
,
6
)
with  probability
21.1
%
13
#10, #11
0.05
96
ξ
¯
~
13
=
{
(
0,1
,
2
)
with  probability
6.5
%
(
1,2
,
3
)
with  probability
17.2
%
(
2,3
,
4
)
with  probability
13.7
%
(
3,4
,
5
)
with  probability
27.2
%
(
4,5
,
6
)
with  probability
35.4
%
14151617181920212325
#1, #5#12, #14#14, #16#9, #14#11, #13#13, #16#10, #13#15, #16#16, #17#18, #20
0.30
90
{
ξ
¯
~
14
ξ
¯
~
15
ξ
¯
~
16
ξ
¯
~
17
ξ
¯
~
18
ξ
¯
~
19
ξ
¯
~
20
ξ
¯
~
21
ξ
¯
~
23
ξ
¯
~
25
=
{
(
0,1
,
2
)
with  probability
12.8
%
(
1,2
,
3
)
with  probability
20.3
%
(
2,3
,
4
)
with  probability
16.5
%
(
3,4
,
5
)
with  probability
31.2
%
(
4,5
,
6
)
with  probability
19.2
%
22
#8, #15
0.09
89
ξ
¯
~
22
=
{
(
0,1
,
2
)
with  probability
7.4
%
(
1,2
,
3
)
with  probability
19.4
%
(
2,3
,
4
)
with  probability
21.6
%
(
3,4
,
5
)
with  probability
23.8
%
(
4,5
,
6
)
with  probability
27.8
%
24
#17, #18
0.10
84
ξ
¯
~
24
=
{
(
0,1
,
2
)
with  probability
23.2
%
(
1,2
,
3
)
with  probability
9.1
%
(
2,3
,
4
)
with  probability
20.8
%
(
3,4
,
5
)
with  probability
21.6
%
(
4,5
,
6
)
with  probability
25.3
%
262729
#20, #21#21, #22#23, #24
0.40
126
{
ξ
¯
~
26
ξ
¯
~
27
ξ
¯
~
29
=
{
(
0,1
,
2
)
with  probability
16.6
%
(
1,2
,
3
)
with  probability
23.2
%
(
2,3
,
4
)
with  probability
7.9
%
(
3,4
,
5
)
with  probability
30.1
%
(
4,5
,
6
)
with  probability
22.2
%
28
#22, #23
0.15
96
ξ
¯
~
28
=
{
(
0,1
,
2
)
with  probability
14.7
%
(
1,2
,
3
)
with  probability
121.3
%
(
2,3
,
4
)
with  probability
10.8
%
(
3,4
,
5
)
with  probability
28.1
%
(
4,5
,
6
)
with  probability
25.1
%In the case problem, there are 2 outputs—the retrofit of permanent (i.e.,i
=
1) and temporary (i.e., i
=
2) links. The activity process and activities for the retrofit are the same for both types of links. Activities (i.e., j
=
1
,
…
,
10) are as (
1
) breaking pavement, (
2
) digging grooves, (
3
) laying pipes, (
4
) backfilling grooves, (
5
) strengthening earth-rock, (
6
) evening roadbeds, (
7
) digging gutters, (
8
) building kerbstones, (
9
) constructing bases, (
10
) constructing pavements. Every activity corresponds to an activity cost center. According to [25], 5 environmental media can be determined in the retrofit work: (
1
) air and climate, (
2
) waste water, (
3
) waste, (
4
) soil and ground water, (
5
) noise and vibration. Environmental cost categories are recorded as in [25].(1)
Waste and emission treatment: depreciation for related equipment; maintenance, operating materials and services; related personnel; fees, taxes, and charges; insurance for environmental liabilities.
(2)
Prevention and environmental management: external services for environmental management and personnel for general environmental management activities.
(3)
Material purchase value of nonproduct output: raw materials; auxiliary materials; operating materials; packaging; energy; and water.
(4)
Processing costs of nonproduct output: depreciation for machinery and labor hours.Since categories 1, 3, and 4 are variable costs (i.e.,v
=
1,2
,
3), category 2 is fixed costs (i.e., f
=
1); they are recoded as (unit:  ¥) [
c
e
1
p
,
c
e
2
p
,
c
e
3
p
]
=
[
11800,5764,2216
], [
c
e
1
t
,
c
e
2
t
,
c
e
3
t
]
=
[
73870,8665,3365
], and c
e
1
f
=
36538. Based on the descriptions above, the environmental costs of the retrofit are allocated as shown in Figure 8.Figure 8
Allocate environmental costs of retrofit work.Based on the practice of the case problem and Figure8, the percentage output for the fixed environmental costs 1 is [
p
e
1
f
,
p
e
2
f
]
=
[
90.8
%
,
9.2
%
]. To calculate the final environmental costs for each output, the corresponding data for the activity cost centers j
=
1
,
…
,
10 is used as in Table 4. The percentage of j in the variable environmental cost categories is denoted as [
p
e
j
1
v
,
p
e
j
2
v
,
p
e
j
3
v
]. a
m
j denotes the cost driver amount in j and [
a
m
1
j
,
a
m
2
j
] denotes the cost driver amount of outputs in j.Table 4
Corresponding data of activity cost centers.
Activity cost centerj
p
e
j
1
v
p
e
j
2
v
p
e
j
3
v
a
m
j
a
m
1
j
a
m
2
j
1
15.6%
0.0%
0.0%
426.30 m2
408.10 m2
18.20 m2
2
6.7%
0.0%
0.0%
9.74 m3
9.32 m3
0.42 m3
3
6.7%
25.0%
25.0%
60.90 m
58.30 m
2.60 m
4
6.7%
0.0%
0.0%
3.35 m3
3.21 m3
0.14 m3
5
5.4%
0.0%
0.0%
25.58 m3
24.49 m3
1.09 m3
6
15.6%
0.0%
0.0%
426.30 m2
408.10 m2
18.20 m2
7
6.7%
0.0%
0.0%
0.61 m3
0.58 m3
0.03 m3
8
5.4%
25.0%
25.0%
60.90 m
58.30 m
2.60 m
9
15.6%
25.0%
25.0%
426.30 m2
408.1 m2
18.20 m2
10
15.6%
25.0%
25.0%
426.30 m2
408.10 m2
18.20 m2In addition, the corresponding cost data for the retrofit and reconstruction are shown in Table5. The values for the other model parameters are as follows: δ
=
0.2, η
=
0.6, ρ
=
1, α
=
0.25, β
=
2, and γ
=
1.Table 5
Cost data of retrofit and reconstruction.
Item
Cost (¥)
Increased variable retrofit costs for permanent link by basic rank (i.e., rank 1)c
v
a
p
16732
Variable retrofit costs for temporary link by basic rank (i.e., rank 1)c
v
a
t
30528
Increased fixed retrofit costs for permanent linkc
f
i
p
14525
Fixed retrofit costs for temporary linkc
f
i
t
28637
Increased variable reconstruction cost for permanent link by basic rank (i.e., rank 1)c
r
v
a
p
107052
Variable reconstruction cost for temporary link by basic rank (i.e., rank 1)c
r
v
a
t
98063
Increased fixed reconstruction cost for permanent linkc
r
f
i
p
69894
Fixed reconstruction cost for temporary linkc
r
f
i
t
50183


## 4.2. Case Solution

The developed algorithm was adopted using MATLAB 7.0 on an Inter Core 2, 2.00 GHz clock pulse with 2048 MB memory. The algorithmic parameters for the case problem were set as follows: error coefficientϵ
=
0.9, swarm_size S
=
20, iteration_max T
=
100, inertia weight_max w
max
=
0.9, inertia weight_min w
min
=
0.1, personal best position acceleration constant c
p
=
0.5, global best position acceleration constant c
g
=
0.5, local best position acceleration constant c
l
=
0.2, and near neighbor best acceleration constant c
n
=
0.1.After 8 iterations of the approximation decomposition, the approximation termination was achieved within 36 minutes on an average of 10 runs, which is time acceptable. The optimal solutions are shown in Tables6 and 7. For all the Pareto optimal solutions on the upper level, the corresponding solutions on the lower level are the same. Table 6 shows a Pareto optimal solution set with 45 solutions on the upper level, where only 10 of the set are enumerated for convenience of expression. The investor is able to choose their preferred plan from the set. If they feel that the retrofit costs including environmental costs C are more important, they would choose the minimum costs plan and vice versa. Since there are fuzzy numbers in model (36), it cannot state crisp optimal objective values in final decision results. However, they are easy to be transformed into equivalent crisp forms by many fuzzy theories and it will not effect the decision.Table 6
Pareto optimal solution set on the upperlevel.
Solutionu
a
1*
2*
3*
4*
5*
6*
7*
8*
9*
10*
……
Linka
1
3
3
3
3
4
4
2
1
3
3
……
2
4
3
4
3
4
4
1
5
2
4
……
3
2
4
2
3
2
1
2
1
5
2
……
4
2
2
3
3
5
1
3
3
1
3
……
5
4
3
3
1
4
4
1
5
3
4
……
6
3
4
3
3
5
2
5
3
3
3
……
7
4
3
3
3
3
5
1
1
2
3
……
8
2
3
2
3
1
1
3
4
3
3
……
9
3
4
3
2
5
1
2
4
4
3
……
10
3
3
3
2
3
3
4
1
3
2
……
11
2
2
2
2
3
3
2
1
1
2
……
12
3
3
4
2
2
3
3
4
2
4
……
13
3
3
3
1
5
3
4
4
2
3
……
14
4
3
4
1
1
5
3
2
3
3
……
15
2
2
1
1
2
1
2
1
3
2
……
16
3
2
3
3
3
3
3
4
1
3
……
17
3
2
3
2
2
4
3
2
2
3
……
18
3
2
3
3
3
4
4
1
1
3
……
19
3
3
2
5
1
2
5
1
3
2
……
20
2
3
2
2
2
2
2
2
4
3
……
21
2
4
2
5
1
1
4
3
5
4
……
22
3
3
3
5
4
1
4
3
4
2
……
23
2
3
2
5
2
1
2
2
3
2
……
24
2
2
3
4
1
2
3
4
2
2
……
25
3
2
3
2
3
3
1
5
2
3
……
26
2
2
2
3
1
1
3
1
1
2
……
27
4
3
4
1
3
4
3
3
2
3
……
28
3
2
3
3
3
5
2
3
1
3
……
29
4
4
4
5
4
5
2
5
3
4
……Table 7
Optimal solutions on the lowerlevel.
Commodityk
1′
2′
3′
4′
5′
6′
7′
8′
9′
10′
11′
12′
Flow of commodityx
k (n/h)
25.50
25.50
12.00
12.00
12.50
12.50
8.00
8.00
9.00
9.00
7.50
7.50


## 4.3. Analytic Results of the Proposed Approach

(
1)  Worthiness of Modeling and Solutions. Fuzzy random programming approach explicitly considers the entire range of uncertain scenarios and thus is more conforming to reality for hedging better against uncertainty. Although it increases the complexity of modeling, the model is well transformed orienting computer implementation. In addition, the computing complexity of the model is not greater than that of the stochastic programming approach used in [1]. Therefore, the extra effort on modeling and solving fuzzy random bilevel programming is worthwhile.The multiobjective method is introduced to determine the Pareto optimal solution set for the upper level and provides more effective and nondominated alternatives for the decision-maker. Compared with the weight-sum method for multiobjective in [19], the solutions in this paper have more reference value for the decision-maker and reflect the users’ preference requirements. Therefore, it is more worthwhile.(
2)  Efficiency of Algorithm. This paper compares GA, an usually used metaheuristics algorithm, and the developed algorithm in this study. The merit of GA is its strong evolutionary process to find an optimal solution by the operation of crossover, mutation, and selection. However, the randomly generated initial generation at the algorithms’ beginning affects the solution quality because of the bad gene inherited from the parent generation. Moreover, the searching capability is reduced as GA does not rely on gradient or derivative information. In the developed algorithm, the particle-represented solutions closely connect the particles of PSO and the solutions to the problem. The hybrid particle-updating mechanism successfully enhances the searching capability. The developed algorithm is a useful tool for the solution to the problem. In contrast to the previous papers, such as [19, 21], both the bilevel and multiobjective environments are considered in this paper.For multiobjective optimization, the definition of quality is substantially more complex than for single-objective optimization problems. Figure9 shows the iterative results of Pareto optimal solutions in 8 iterations. Note that in each iteration of approximation decomposition, these fuzzy numbers in model (36) are decomposed into crisp forms. Therefore, the objective values can be evaluated. For further expression of the efficiency of the convergence, three metrics of performance are studied. Table 8 shows the metrics of performance for Pareto optimal sets proposed in [51] and the set convergence of Pareto optimal sets in 8 iterations after 10 runs.Table 8
Metrics of performance and set convergence for Pareto optimal sets.
Iteration
The average distance
The distribution
The extent
The set convergence
1
0.0812
0.5623
5.0088
—
2
0.0630
0.4615
5.7706
Iteration 1–2: 0.5610
3
0.1364
0.8600
6.8594
Iteration 2–3: 0.6765
4
0.1223
0.5250
6.3147
Iteration 3–4: 0.7632
5
0.1502
0.6820
6.7538
Iteration 4–5: 0.8750
6
0.2835
0.9600
6.6417
Iteration 5–6: 0.8571
7
0.0973
0.5289
5.4451
Iteration 6–7: 0.5712
8
0.1300
0.8635
6.6424
Iteration 7–8: 0.9143Iterative results of Pareto optimal solutions in 8 iterations.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)(
3)  Efficiency of Algorithm Parameters. Since the search space for the problem is so large and the computing process is time-consuming, it is necessary to choose reasonable algorithm parameters. Table 9 shows the comparison results of different swarm_size (i.e., S) and iteration_max (i.e., T) on average computing time and iterations after 10 runs. Looking into Table 9, when S
=
10, the program could not obtain results in more than 3600 s and more than 10 iterations with both T
=
100 and T
=
200. This is the same when S
=
50 with the lower iterations. When S
=
20 and T
=
200, the process is more time-consuming than the current algorithm parameters. Therefore, the current ones are considered suitable.Table 9
Comparison of different algorithm parameters.
S
=
10
S
=
20
S
=
50
T
=
100
T
=
200
T
=
100
T
=
200
T
=
100
T
=
200
Computing time (s)
≥3600
≥3600
2160
3480
≥3600
≥3600
Iterations
≥16
≥12
8
6
≥4
≥3


## 5. Conclusions

This paper studies retrofit decision for transportation network of LSCP to hedge against seismic risk. On the consideration of the emerged challenges in the problem, using distinctions of various link types, bilevel decision, environmental costs, and fuzzy random seismic damage scenario, a fuzzy random multiobjective bilevel programming model is set up. A transforming approach is in use to obtain equivalent fuzzy bilevel programming model. Then, decomposition is utilized to these fuzzy variables by decomposition theorem of fuzzy number. An approximation decomposition-based multiobjective AGLNPSO is developed to solve the problem. A case study is presented as an illustrative example of this problem. The results validate the worthiness of modeling and solutions and test the efficiency of the algorithm and parameters.The contributions of this paper to literature are as the follows.(
1
) This study adopts the work of [1] to the filed of LSCP. Bilevel decision involves, with the investor and the administer, distinctions of various link types, and retrofit decision specified into several ranks according to the seismic damage scenario provides more reasonable and practical description of the problem. (
2
) Although there are many works on environmental cost, such as [25, 26], few papers consider it in transportation network in LSCP. Thus it enhances the focus points of management aims of the problem. (
3
) This paper uses fuzzy random seismic damage scenario to describe the hybrid uncertain situation. To the best of our knowledge, it has never been done before. (
4
) The approximation decomposition-based multiobjective AGLNPSO is developed as one of the useful tools to solve this problem. In contrast to the previous papers, such as [19, 21], both the bilevel and multiobjective environments are considered in this paper.There are three areas suggested for future research. First, more cost categories need to be investigated and the detailed relationships between the retrofit decision and the costs should be outlined to ensure the model is as practical as possible. Secondly, to determine better, more effective solutions with lower memory and computing time requirements, alternative approaches and algorithms (e.g., other exact approaches, (meta)heuristics, evolutionary algorithms, etc.) could be used to make comparisons. Finally, consideration of other behavior assumptions, such as the travelers learning or user equilibrium, may change the structure of the problem. Each of these areas is very important and equally worthy of our concern.


## Conflict of Interests

The authors declare no conflict of interests.



---

*Source: 505890-2014-01-12.xml*