# A Nonlinear Multiobjective Bilevel Model for Minimum Cost Network Flow Problem in a Large-Scale Construction Project

Authors: Jiuping Xu; Yan Tu; Ziqiang Zeng
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

The aim of this study is to deal with a minimum cost network flow problem (MCNFP) in a large-scale construction project using a nonlinear multiobjective bilevel model with birandom variables. The main target of the upper level is to minimize both direct and transportation time costs. The target of the lower level is to minimize transportation costs. After an analysis of the birandom variables, an expectation multiobjective bilevel programming model with chance constraints is formulated to incorporate decision makers’ preferences. To solve the identified special conditions, an equivalent crisp model is proposed with an additional multiobjective bilevel particle swarm optimization (MOBLPSO) developed to solve the model. The Shuibuya Hydropower Project is used as a real-world example to verify the proposed approach. Results and analysis are presented to highlight the performances of the MOBLPSO, which is very effective and efficient compared to a genetic algorithm and a simulated annealing algorithm.


---

## Body

## 1. Introduction

Network flow optimization is a large part of combinatorial optimization. The minimum cost network flow problem (MCNFP) is made up of a wide category of problems [1, 2]. MCNFP plays a very important role in many real-world applications such as communications [3, 4], informatics [5], and transportation [6]. Other well-known problems like the shortest path problem and the assignment problem are considered to be special MCNFP cases [7].In recent decades, the MCNFP has been well researched with many models and algorithms being developed, for example, [8–13]. These studies, however, have not often taken carrier type selection and transportation time into account when looking at the transportation network. Yet, both cost and time control are important in construction projects, especially in large-scale construction projects where transportation costs are largely based on the rates charged by carriers, which have a significant influence on transportation time. In real conditions, there is increasing pressure to shorter transportation time to reduce or eliminate extra project expenses, with the early arrival of materials shortening the completion time of the construction project and improving construction efficiency. In these cases, it is necessary to include both the carrier type selection or transportation time in the transportation network analysis. In this paper, a multiobjective bilevel MCNFP is studied. On the upper level, the construction contractor determines the material flow of each transportation network path with the criteria being the minimization of both direct costs and total transportation time costs. On the lower level, the transportation manager controls each carrier's flow so that total transportation costs are minimized.The research presented previously has a common foundation in that they were all based on a deterministic transportation network. However, transportation systems are often complex, so decision makers inevitably encounter uncertain parameters when making a decision. Within the last two decades the use of multimodels with uncertain parameters in the study of network flow problems has been increasingly exploited. For example, Watling [14] studied a user equilibrium traffic network assignment problem with stochastic travel times and a late arrival penalty. Chen and Zhou [15] developed an α-reliable mean-excess traffic equilibrium model with stochastic travel times. Lin [16] constructed a revised stochastic flow network to model a realistic computer network in which each arc has a lead time and a stochastic capacity. Sumalee et al. [17] dealt with a reliable network design problem which was looked at uncertain demand and total travel time reliability. In actual analyses, randomness is considered one important source of uncertainty. Yet with MCNFP randomness is seen to be increasingly complex because of the often incomplete or uncertain information. To date, there has been little research which considers multilevel twofold uncertainty coefficients for MCNFP. Therefore this research concentrates on the problem under a birandom environment with the logic behind this choice of birandom variables illustrated in Section 2.The MCNFP proposed in this paper is a multiobjective bilevel programming problem, first introduced by Geoffrion and Hogan [18], and consequently developed by researchers such as Tarvainen and Haimes [19], Osman et al. [20], Zhang et al. [21], and Calvete and Galéb [22]. Multiobjective bilevel programming has been greatly improved in both the theoretical and practical areas. While these studies have significantly contributed to a variety of applications, to the best of our knowledge, there is still no known research considering the modeling for the MCNFP. With bilevel programming problems being intrinsically difficult, it is not surprising that most exact algorithms to date have focused on the simplest cases of bilevel programs, that is, problems with relatively easy to determine properties such as linear, quadratic, or convex objective and/or constraint functions [23]. Since the proposed bilevel MCNFP model is nonlinear, nonconvex, and nondifferentiable, it follows that the search for exact algorithms which are formally efficient is all but futile and it is necessary instead to search for effective heuristic algorithms to solve the MCNFP. Determining the global optimal solution is of great importance in MCNFP. Specifically, this paper deals with the multiple objectives by employing the concept of nondominated solutions instead of applying weighted sum scalarization. In this study, an effort is made to develop a multiobjective bilevel particle swarm optimization (MOBLPSO) to solve a real world MCNFP in the Shuibuya Hydropower Project.The remainder of this paper is structured as follows. In Section2, an introduction to the bilevel MCNFP is presented along with the motivation for employing birandom variables in the problem. An expectation multiobjective bilevel programming model with chance constraints under a birandom environment is established in Section 3, for which the equivalent crisp model is derived in Section 4. In addition, an MOBLPSO is illustrated in Section 5. In Section 6, an application to earth-rock work transportation in a large construction project is given in order to show the validity and efficiency of the proposed models and algorithms. Concluding remarks and further discussion are in Section 7.


## 2. Problem Statement

In construction projects, and especially in large-scale construction projects, the MCNFP is becoming increasingly important. Here we discuss a multiobjective bilevel MCNFP in a large-scale construction project. In order to establish the model, a description is given.


### 2.1. Bilevel Problem Description

The MCNFP discussed considers both the construction contractor and the transportation manager as the two participants. In a large-scale construction project, material often has supply (origin) and receipt (destination) nodes, with the construction contractor generally assigning a specialized Transportation Company. The bilevel model considers the construction contractor and the transportation manager in the specialized Transportation Company concurrently, gives priority to contractor benefit, and considers the influence of the contractor's decision making on the flow distribution of the transportation manager’s carriers. As both cost and time control are important in construction projects, their effectiveness needs to be considered. The construction contractor assigns the flow of material to each transportation path to minimize direct costs and transportation time costs, while the transportation manager aims to minimize transportation costs by making decisions about flow of material by each carrier through the transportation path based on the construction contractor’s decision making, which in turn influences contractor’s decision-making through adjustments to the flow of material by each carrier along the transportation path.Therefore, the MCNFP in this paper can be abstracted as a bilevel programming problem. To model the problem conveniently, the involved transportation network is considered a bipartite network represented by a graph with sets of nodes and arcs. In the network, a node represents the facilities in the network, for instance, a station or a yard, and an arc represents a line between two adjacent facilities. The model structure of the MCNFP is in Figure1.Figure 1
Model structure for MCNFP.


### 2.2. The Motivation for Considering Birandom Environment in the MCNFP

The birandom environment has been successfully studied and applied in many areas, such as flow shop scheduling problem [24], portfolio selection [25], and vendor selection [26]. These studies show the necessity of considering birandom environment in practical problems. There is a strong motivation for considering birandom environment for the MCNFP.In real conditions, the transportation plan is usually made before the occurrence of any transportation activity; thus the determined values of some parameters cannot be obtained in advance; so there is a strong need to consider uncertainty in transportation problems. An MCNFP in a large-scale construction project is considered in this paper, which may be subjected to twofold randomness with incomplete or uncertain information. For example, the transportation time in an arc of the project is not fixed because of the effect of the transportation environment which includes many uncertainties such as traffic accidents, traffic congestion, vehicle breakdowns, bad weather, natural disasters, and special events [27]. Therefore, it is often subject to a stochastic distribution. Generally, transportation time approximately follows a normal distribution expressed as 𝒩(μ̃,σ2) [28–30], whereas the normal distribution has to be truncated to avoid negative values. However, the expected value of transportation time μ̃ is also uncertain because it is determined by the carrier's speed, which in turn is influenced by such uncertainties as vehicle condition. It is possible to specify a realistic distribution (e.g., normal distribution) for the parameter μ̃ through statistical methods or related expertise and other knowledge. When the value of μ̃ is provided as a random variable which approximately follows a normal distribution, the pattern of overlapping randomness is said to be birandom as illustrated in Section 3.1. The flowchart of transportation time as a birandom variable is in Figure 2.Figure 2
Flowchart of transportation time as a birandom variable.The situation is similar with transportation costs. However, the mean values of transportation costs are considered as approximately following normal distributions due to the fluctuation of gasoline prices over time, which results in the birandomness of the transportation costs. From the previous description, the birandom variable is employed to take account of the hybrid uncertainty and obtain a more feasible network flow scheme.


## 2.1. Bilevel Problem Description

The MCNFP discussed considers both the construction contractor and the transportation manager as the two participants. In a large-scale construction project, material often has supply (origin) and receipt (destination) nodes, with the construction contractor generally assigning a specialized Transportation Company. The bilevel model considers the construction contractor and the transportation manager in the specialized Transportation Company concurrently, gives priority to contractor benefit, and considers the influence of the contractor's decision making on the flow distribution of the transportation manager’s carriers. As both cost and time control are important in construction projects, their effectiveness needs to be considered. The construction contractor assigns the flow of material to each transportation path to minimize direct costs and transportation time costs, while the transportation manager aims to minimize transportation costs by making decisions about flow of material by each carrier through the transportation path based on the construction contractor’s decision making, which in turn influences contractor’s decision-making through adjustments to the flow of material by each carrier along the transportation path.Therefore, the MCNFP in this paper can be abstracted as a bilevel programming problem. To model the problem conveniently, the involved transportation network is considered a bipartite network represented by a graph with sets of nodes and arcs. In the network, a node represents the facilities in the network, for instance, a station or a yard, and an arc represents a line between two adjacent facilities. The model structure of the MCNFP is in Figure1.Figure 1
Model structure for MCNFP.


## 2.2. The Motivation for Considering Birandom Environment in the MCNFP

The birandom environment has been successfully studied and applied in many areas, such as flow shop scheduling problem [24], portfolio selection [25], and vendor selection [26]. These studies show the necessity of considering birandom environment in practical problems. There is a strong motivation for considering birandom environment for the MCNFP.In real conditions, the transportation plan is usually made before the occurrence of any transportation activity; thus the determined values of some parameters cannot be obtained in advance; so there is a strong need to consider uncertainty in transportation problems. An MCNFP in a large-scale construction project is considered in this paper, which may be subjected to twofold randomness with incomplete or uncertain information. For example, the transportation time in an arc of the project is not fixed because of the effect of the transportation environment which includes many uncertainties such as traffic accidents, traffic congestion, vehicle breakdowns, bad weather, natural disasters, and special events [27]. Therefore, it is often subject to a stochastic distribution. Generally, transportation time approximately follows a normal distribution expressed as 𝒩(μ̃,σ2) [28–30], whereas the normal distribution has to be truncated to avoid negative values. However, the expected value of transportation time μ̃ is also uncertain because it is determined by the carrier's speed, which in turn is influenced by such uncertainties as vehicle condition. It is possible to specify a realistic distribution (e.g., normal distribution) for the parameter μ̃ through statistical methods or related expertise and other knowledge. When the value of μ̃ is provided as a random variable which approximately follows a normal distribution, the pattern of overlapping randomness is said to be birandom as illustrated in Section 3.1. The flowchart of transportation time as a birandom variable is in Figure 2.Figure 2
Flowchart of transportation time as a birandom variable.The situation is similar with transportation costs. However, the mean values of transportation costs are considered as approximately following normal distributions due to the fluctuation of gasoline prices over time, which results in the birandomness of the transportation costs. From the previous description, the birandom variable is employed to take account of the hybrid uncertainty and obtain a more feasible network flow scheme.


## 3. Modelling

In this section, the relevant assumptions and notations are first outlined, some basic knowledge about the birandom variable is introduced, and then the nonlinear multiobjective bilevel MCNFP model under a birandom environment is formulated.


### 3.1. Problem Assumptions

The mathematical model described has the following assumptions.(1)
The proposed transportation network is a single material transportation network composed of nodes and road sections.(2)
The capacities of the different arcs are satisfactorily independent, and the total flow of carriers on each arc cannot exceed its capacity.(3)
Flows on all transportation paths between OD pairs satisfy the feasible flow conservation [31].(4)
All transportation paths in the network are known.(5)
The transportation cost of every road section and transportation time are considered birandom variables, with the attributes determined from available statistics and historical data as well as forecast transportation environments. They are considered to be independent.(6)
The demand at every reception node must be met on schedule. The material has a given transportation duration. If the transportation time exceeds the given duration, then a delay cost is added.


### 3.2. Notations

The following mathematical notations are used to describe the MCNFP.Subscripts and Sets
o:
index of origin node;d:
index of destination node;Ψ:
set of carriers,k∈Ψ is an index;Φ:
set of arcs in the transportation network,i∈Φ is an index;Ω:
set of paths in the transportation network,j∈Ω is an index;E:
set of origin-destination (OD) pairs,(o,d)∈E;Aj:
set of arcs in transportation pathj;Pod:
set of paths from origin nodeo to destination node d.Certain Parameters
ri:
maximal passing capacity of arci;wk:
weight for carrierk in the transportation network;vk:
volume capacity of carrierk;Tj:
transportation time constraint represents the time of transportation pathj  in Tj units in which the material demand between all OD pairs has to be transported;γij:
a binary variable equal to 1 if and only if arci is a segment of transportation path j for carrier k;Qod:
transportation demand of the material from origin nodeo to destination node d;cj:
direct cost of unit volume material using transportationj;⌈⌉:
ceiling operator rounding upward to integer.Uncertain Parameters
e¯̃ik:
unit transportation cost of material flow on arcai for carrier k;t¯̃i0k:
free transportation time of material flow on arci for carrier k;t¯̃ik:
transportation time of material flow on arcai for carrier k.Decision Variables
xj:
volume of material flow on transportation pathj, which is the decision variable of the upper level;ykj:
volume of material flow transported by carrierk through path j, which is the decision variable of the lower level.


### 3.3. Birandom Variable

The birandom variable, proposed by Peng and Liu [32], can be used to explain the proposed problem. Some basic knowledge about birandom variable is as follows.Definition 3.1 (see [32]).
A birandom variableξ is a mapping from a probability space (Ω,𝒜,Pr) to a collection S of random variables such that for any Borel subset B of the real line ℜ the induced function Pr{ξ(ω)∈B} is a measurable function with respect to ξ.
For each given Borel subsetB of the real line ℜ, the function Pr{ξ(ω)∈B} is a random variable defined on the probability space (Ω,𝒜,Pr).From Definition3.1, a birandom variable is a mapping from a probability space to a collection of random variables. Roughly speaking, a birandom can be seen as a random variable with random parameter(s). Here we give three examples of birandom variable.Example 3.2.
LetΩ={ω1,ω2,…,ωm}, and Pr{ω1}+Pr{ω2}+⋯+Pr{ωm-1}+Pr{ωm}=1. Assume that ξ is a function on (Ω,𝒜,Pr) as follows:
(3.1)ξ(ω)={ξ̃1,ifω=ω1,ξ̃2,ifω=ω2,⋮ξ̃m-1,ifω=ωm-1,ξ̃m,ifω=ωm,
where ξ̃1 is a uniformly distributed random variable on [0,1], ξ̃2,…,ξ̃m-1 are normally distributed random variables with a mean 1 and a standard variance 0.5, and ξ̃m is a standard normally distributed random variable with a mean 0 and a standard variance 1, that is, ξ̃1~U[0,1], ξ̃2~𝒩[1,0.5],…,ξ̃m-1~𝒩[1,0.5], and ξ̃m~𝒩[0,1]. From Definition 3.1, ξ is clearly a birandom variable as shown in Figure 3.Figure 3
Representation of a Birandom variable in Example3.2.Example 3.3.
Assume thata¯ and b¯ are two random variables defined on (Ω′,𝒜′,Pr′), and for any ω*∈Ω′, b¯(ω*)≥a¯(ω*) holds; then random variable ξ~𝒩[a¯(ω*),b¯(ω*)] is a birandom variable.Example 3.4.
Letξ be a random variable defined on the probability space (Ω,𝒜,Pr) satisfying ξ∈𝒩(μ̃,σ2), where μ̃ is also a normally distributed random variable on (Ω′,𝒜′,Pr′) with the mean μ and variance σ*2. Then ξ is a birandom variable.Definition 3.5.
A birandom variableξ is said to be normal, if for each ω, ξ(ω) is a random variable with the following probability density function:
(3.2)ϕ(x)=1σ(ω)2πexp⁡(-(x-μ(ω))22σ(ω)2),
where the number of random variable of μ(ω) and σ(ω) is not less than one. The normal birandom variable is denoted by 𝒩(μ(ω)andσ(ω)).


### 3.4. Bilevel Model Formulation

In this section, the expectation multiobjective bilevel programming model with chance constraints under a birandom environment based on the philosophy is proposed: decisions are selected by optimizing the expected values of the objective functions subject to chance constraints with some predetermined confidence levels given by the actual decision makers.


#### 3.4.1. Lower Level Model for the Bilevel MCNFP

The problem posed on the lower level is how to make decisions on material flow by each type of carrier on transportation pathykj while satisfying all capacity constraints, with the main objective being to minimize expected total transportation cost.(1) Objective Functions of the Lower Level
The total transportation cost of the material is calculated by taking the sum of the carriers of each arc's transportation cost and the number of carriers needed to transport the material across the network. In real conditions, it is desirable that each service carrier be fully loaded, so the numbers of carrierk through path j can be denoted by ⌈ykj/vk⌉. Since e¯̃ik is considered a birandom variable, the total transportation cost of material is considered under a birandom environment. Generally, it is difficult to completely minimize total transportation costs because of the birandom variables. Because decision-makers expect minimal cost, the expected value of the total transportation cost is the objective of the lower level. Denote the expected total transportation cost of material as C(ykj); then the objective function of the lower level model can be formulated as
(3.3)min⁡C(ykj)=E[∑k∈Ψ∑j∈Ω∑i∈Φγije¯̃ik⌈ykjvk⌉].
Generally, a path can be represented by a sequence of adjacent arcs. A binary variableγij is introduced to determine whether an arc i is a segment of path j for carrier k:
(3.4)γij={1,ifi∈Aj,j∈Ω,0,otherwise.(2) Constraints of the Lower Level
For transportation time, each carrier requires the transport of material from the source to the destination on scheduleTj. If not, a delay cost is applied. ∑i∈Φγijt¯̃ik represents the total travel time of carrier k on transportation path j, in which t¯̃ik is usually represented by a non-decreasing function (i.e., Bureau of Public Roads (BPR) function) [33] as follows:
(3.5)t¯̃ik=t¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β],i∈Φ,k∈Ψ,
where α and β are user-defined parameters and, in this problem, are set to 0.15 and 2.0, respectively.
Technically, it is not possible to strictly ensure that the random event∑i∈Φγijt¯̃ik does not exceed Tj because of the birandom variable t¯̃i0k. In practical problem, the decision makers often provide an appropriate budget Tj in advance, to ensure that the restriction is, to a certain extent, satisfied, that is to maximize the probability of the random event Pr{∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β]≤Tj}under a given confidence level, which can be written as follows:
(3.6)Pr{ω∣Pr{∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]≤Tj}≥θ}≥δ,j∈Ω,k∈Ψ.
Here the decision makers' aspiration level is indicated as θ, so we use a “Pr” to ensure that the constraint holds at the predetermined confidence level. Additionally, based on probability theory, a further “Pr” is needed to describe the random elements presented in Section 2, which guarantee the establishment of a certain confidence level δ, resembling the P-model (probability maximization model) presented in [34].
The transportation flow may exceed some arcs' capacity because of uncertainties such as the condition of the construction project road. Such conditions may require the manager to select another path. Thus, the total amount of capacity on arci cannot exceed the maximal capacity of the arc i, which produces the following constraint:(3.7)∑j∈Ω∑k∈Ψγij⌈ykjvk⌉≤ri,i∈Φ.
Actually, it is difficult to ensure that each service carrier is fully loaded; so the sum of all flows transported by all kinds of carriers through each path cannot be less than the material flow assigned to it; thus the following constraint is obtained:(3.8)∑k∈Ψykj≥xj,j∈Ω.
The path flow in pathj used by carrier k should not be negative, such that
(3.9)ykj≥0,j∈Ω,k∈Ψ.


#### 3.4.2. Upper Level Model for the Bilevel MCNFP

The problem the construction contractor on the upper level faces is how to assign material flow among the transportation paths across the complete transportation network, that is, how to decide the material flowxj through transportation path j. Thus, the decision variable on the upper level is xj.(1) Objective Functions of the Upper Level
For large-scale construction projects, cost and time control are both important, so minimizing total direct costs and total transportation time costs is the two objectives of the upper level model. The two objectives of the upper level can be described as follows.
Firstly, the upper level decision maker attempts to minimize the direct costs of the complete network by assigning the flow of the material to each transportation path to achieve a system optimized flow pattern; thus, the total direct cost is the sum of all transportation costs from different transportation paths. The first objective function of the upper level model can be formulated as follows(3.10)C(xj,ykj)=∑(o,d)∈E∑j∈Podcjxj.
In real conditions, there is increasing pressure to shorter transportation time to reduce or eliminate extra project expenses, with the early arrival of materials shortening the completion time of the construction project and improving construction efficiency. Thus, the total transportation time for carrierk in each path can be described as ∑i∈Φγijt¯̃ik, and, since t¯̃ik is a birandom variable, ∑i∈Φγijt¯̃ik can be regarded as a special birandom variable. Similarly, the expected value of the total transportation time cost is one of the objectives on the upper level. Different carriers are given different weights, and t¯̃ik=t¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β], so the second objective function of the upper level can be described as follows:
(3.11)T(xj,ykj)=E[∑k∈Ψwk∑j∈Ω∑i∈Φγijt¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]].(2) Constraints of the Upper Level
According to the basic assumptions in Section3.1, it is stipulated that the demands between all OD pairs should be satisfied. Thus,
(3.12)∑j∈Podxj=Qod,(o,d)∈E.
The following constraint ensures that the sum of the weights is equal to 1:(3.13)∑k∈Ψwk=1.
In order to describe the nonnegative variables, the constraints in (3.14) are presented:
(3.14)xj≥0,j∈Ω.


#### 3.4.3. Global Model for the Bilevel MCNFP

Based on the previous discussion, by integrating (3.3)–(3.14), the following global model for the nonlinear multiobjective bilevel programming with birandom variables is formulated for the MCNFP in a large-scale construction project:
(3.15)min⁡C(xj,ykj)=∑(o,d)∈E∑j∈Podcjxj,min⁡T(xj,ykj)=E[∑k∈Ψwk∑j∈Ω∑i∈Φγijt¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]],s.t.{∑j∈Podxj=Qod,(o,d)∈E,∑k∈Ψwk=1,xj≥0,j∈Ω,min⁡C(ykj)=E[∑k∈Ψ∑j∈Ω∑i∈Φγije¯̃ik⌈ykjvk⌉],s.t.{Pr{ω∣Pr{∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]≤Tj}≥θ}≥δ,j∈Ω,k∈Ψ,∑j∈Ω∑k∈Ψγij⌈ykjvk⌉≤ri,i∈Φ,γij={1,ifi∈Aj,j∈Ω,0,otherwise,∑k∈Ψykj≥xj,j∈Ω,ykj≥0,j∈Ω,k∈Ψ.


## 3.1. Problem Assumptions

The mathematical model described has the following assumptions.(1)
The proposed transportation network is a single material transportation network composed of nodes and road sections.(2)
The capacities of the different arcs are satisfactorily independent, and the total flow of carriers on each arc cannot exceed its capacity.(3)
Flows on all transportation paths between OD pairs satisfy the feasible flow conservation [31].(4)
All transportation paths in the network are known.(5)
The transportation cost of every road section and transportation time are considered birandom variables, with the attributes determined from available statistics and historical data as well as forecast transportation environments. They are considered to be independent.(6)
The demand at every reception node must be met on schedule. The material has a given transportation duration. If the transportation time exceeds the given duration, then a delay cost is added.


## 3.2. Notations

The following mathematical notations are used to describe the MCNFP.Subscripts and Sets
o:
index of origin node;d:
index of destination node;Ψ:
set of carriers,k∈Ψ is an index;Φ:
set of arcs in the transportation network,i∈Φ is an index;Ω:
set of paths in the transportation network,j∈Ω is an index;E:
set of origin-destination (OD) pairs,(o,d)∈E;Aj:
set of arcs in transportation pathj;Pod:
set of paths from origin nodeo to destination node d.Certain Parameters
ri:
maximal passing capacity of arci;wk:
weight for carrierk in the transportation network;vk:
volume capacity of carrierk;Tj:
transportation time constraint represents the time of transportation pathj  in Tj units in which the material demand between all OD pairs has to be transported;γij:
a binary variable equal to 1 if and only if arci is a segment of transportation path j for carrier k;Qod:
transportation demand of the material from origin nodeo to destination node d;cj:
direct cost of unit volume material using transportationj;⌈⌉:
ceiling operator rounding upward to integer.Uncertain Parameters
e¯̃ik:
unit transportation cost of material flow on arcai for carrier k;t¯̃i0k:
free transportation time of material flow on arci for carrier k;t¯̃ik:
transportation time of material flow on arcai for carrier k.Decision Variables
xj:
volume of material flow on transportation pathj, which is the decision variable of the upper level;ykj:
volume of material flow transported by carrierk through path j, which is the decision variable of the lower level.


## 3.3. Birandom Variable

The birandom variable, proposed by Peng and Liu [32], can be used to explain the proposed problem. Some basic knowledge about birandom variable is as follows.Definition 3.1 (see [32]).
A birandom variableξ is a mapping from a probability space (Ω,𝒜,Pr) to a collection S of random variables such that for any Borel subset B of the real line ℜ the induced function Pr{ξ(ω)∈B} is a measurable function with respect to ξ.
For each given Borel subsetB of the real line ℜ, the function Pr{ξ(ω)∈B} is a random variable defined on the probability space (Ω,𝒜,Pr).From Definition3.1, a birandom variable is a mapping from a probability space to a collection of random variables. Roughly speaking, a birandom can be seen as a random variable with random parameter(s). Here we give three examples of birandom variable.Example 3.2.
LetΩ={ω1,ω2,…,ωm}, and Pr{ω1}+Pr{ω2}+⋯+Pr{ωm-1}+Pr{ωm}=1. Assume that ξ is a function on (Ω,𝒜,Pr) as follows:
(3.1)ξ(ω)={ξ̃1,ifω=ω1,ξ̃2,ifω=ω2,⋮ξ̃m-1,ifω=ωm-1,ξ̃m,ifω=ωm,
where ξ̃1 is a uniformly distributed random variable on [0,1], ξ̃2,…,ξ̃m-1 are normally distributed random variables with a mean 1 and a standard variance 0.5, and ξ̃m is a standard normally distributed random variable with a mean 0 and a standard variance 1, that is, ξ̃1~U[0,1], ξ̃2~𝒩[1,0.5],…,ξ̃m-1~𝒩[1,0.5], and ξ̃m~𝒩[0,1]. From Definition 3.1, ξ is clearly a birandom variable as shown in Figure 3.Figure 3
Representation of a Birandom variable in Example3.2.Example 3.3.
Assume thata¯ and b¯ are two random variables defined on (Ω′,𝒜′,Pr′), and for any ω*∈Ω′, b¯(ω*)≥a¯(ω*) holds; then random variable ξ~𝒩[a¯(ω*),b¯(ω*)] is a birandom variable.Example 3.4.
Letξ be a random variable defined on the probability space (Ω,𝒜,Pr) satisfying ξ∈𝒩(μ̃,σ2), where μ̃ is also a normally distributed random variable on (Ω′,𝒜′,Pr′) with the mean μ and variance σ*2. Then ξ is a birandom variable.Definition 3.5.
A birandom variableξ is said to be normal, if for each ω, ξ(ω) is a random variable with the following probability density function:
(3.2)ϕ(x)=1σ(ω)2πexp⁡(-(x-μ(ω))22σ(ω)2),
where the number of random variable of μ(ω) and σ(ω) is not less than one. The normal birandom variable is denoted by 𝒩(μ(ω)andσ(ω)).


## 3.4. Bilevel Model Formulation

In this section, the expectation multiobjective bilevel programming model with chance constraints under a birandom environment based on the philosophy is proposed: decisions are selected by optimizing the expected values of the objective functions subject to chance constraints with some predetermined confidence levels given by the actual decision makers.


### 3.4.1. Lower Level Model for the Bilevel MCNFP

The problem posed on the lower level is how to make decisions on material flow by each type of carrier on transportation pathykj while satisfying all capacity constraints, with the main objective being to minimize expected total transportation cost.(1) Objective Functions of the Lower Level
The total transportation cost of the material is calculated by taking the sum of the carriers of each arc's transportation cost and the number of carriers needed to transport the material across the network. In real conditions, it is desirable that each service carrier be fully loaded, so the numbers of carrierk through path j can be denoted by ⌈ykj/vk⌉. Since e¯̃ik is considered a birandom variable, the total transportation cost of material is considered under a birandom environment. Generally, it is difficult to completely minimize total transportation costs because of the birandom variables. Because decision-makers expect minimal cost, the expected value of the total transportation cost is the objective of the lower level. Denote the expected total transportation cost of material as C(ykj); then the objective function of the lower level model can be formulated as
(3.3)min⁡C(ykj)=E[∑k∈Ψ∑j∈Ω∑i∈Φγije¯̃ik⌈ykjvk⌉].
Generally, a path can be represented by a sequence of adjacent arcs. A binary variableγij is introduced to determine whether an arc i is a segment of path j for carrier k:
(3.4)γij={1,ifi∈Aj,j∈Ω,0,otherwise.(2) Constraints of the Lower Level
For transportation time, each carrier requires the transport of material from the source to the destination on scheduleTj. If not, a delay cost is applied. ∑i∈Φγijt¯̃ik represents the total travel time of carrier k on transportation path j, in which t¯̃ik is usually represented by a non-decreasing function (i.e., Bureau of Public Roads (BPR) function) [33] as follows:
(3.5)t¯̃ik=t¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β],i∈Φ,k∈Ψ,
where α and β are user-defined parameters and, in this problem, are set to 0.15 and 2.0, respectively.
Technically, it is not possible to strictly ensure that the random event∑i∈Φγijt¯̃ik does not exceed Tj because of the birandom variable t¯̃i0k. In practical problem, the decision makers often provide an appropriate budget Tj in advance, to ensure that the restriction is, to a certain extent, satisfied, that is to maximize the probability of the random event Pr{∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β]≤Tj}under a given confidence level, which can be written as follows:
(3.6)Pr{ω∣Pr{∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]≤Tj}≥θ}≥δ,j∈Ω,k∈Ψ.
Here the decision makers' aspiration level is indicated as θ, so we use a “Pr” to ensure that the constraint holds at the predetermined confidence level. Additionally, based on probability theory, a further “Pr” is needed to describe the random elements presented in Section 2, which guarantee the establishment of a certain confidence level δ, resembling the P-model (probability maximization model) presented in [34].
The transportation flow may exceed some arcs' capacity because of uncertainties such as the condition of the construction project road. Such conditions may require the manager to select another path. Thus, the total amount of capacity on arci cannot exceed the maximal capacity of the arc i, which produces the following constraint:(3.7)∑j∈Ω∑k∈Ψγij⌈ykjvk⌉≤ri,i∈Φ.
Actually, it is difficult to ensure that each service carrier is fully loaded; so the sum of all flows transported by all kinds of carriers through each path cannot be less than the material flow assigned to it; thus the following constraint is obtained:(3.8)∑k∈Ψykj≥xj,j∈Ω.
The path flow in pathj used by carrier k should not be negative, such that
(3.9)ykj≥0,j∈Ω,k∈Ψ.


### 3.4.2. Upper Level Model for the Bilevel MCNFP

The problem the construction contractor on the upper level faces is how to assign material flow among the transportation paths across the complete transportation network, that is, how to decide the material flowxj through transportation path j. Thus, the decision variable on the upper level is xj.(1) Objective Functions of the Upper Level
For large-scale construction projects, cost and time control are both important, so minimizing total direct costs and total transportation time costs is the two objectives of the upper level model. The two objectives of the upper level can be described as follows.
Firstly, the upper level decision maker attempts to minimize the direct costs of the complete network by assigning the flow of the material to each transportation path to achieve a system optimized flow pattern; thus, the total direct cost is the sum of all transportation costs from different transportation paths. The first objective function of the upper level model can be formulated as follows(3.10)C(xj,ykj)=∑(o,d)∈E∑j∈Podcjxj.
In real conditions, there is increasing pressure to shorter transportation time to reduce or eliminate extra project expenses, with the early arrival of materials shortening the completion time of the construction project and improving construction efficiency. Thus, the total transportation time for carrierk in each path can be described as ∑i∈Φγijt¯̃ik, and, since t¯̃ik is a birandom variable, ∑i∈Φγijt¯̃ik can be regarded as a special birandom variable. Similarly, the expected value of the total transportation time cost is one of the objectives on the upper level. Different carriers are given different weights, and t¯̃ik=t¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β], so the second objective function of the upper level can be described as follows:
(3.11)T(xj,ykj)=E[∑k∈Ψwk∑j∈Ω∑i∈Φγijt¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]].(2) Constraints of the Upper Level
According to the basic assumptions in Section3.1, it is stipulated that the demands between all OD pairs should be satisfied. Thus,
(3.12)∑j∈Podxj=Qod,(o,d)∈E.
The following constraint ensures that the sum of the weights is equal to 1:(3.13)∑k∈Ψwk=1.
In order to describe the nonnegative variables, the constraints in (3.14) are presented:
(3.14)xj≥0,j∈Ω.


### 3.4.3. Global Model for the Bilevel MCNFP

Based on the previous discussion, by integrating (3.3)–(3.14), the following global model for the nonlinear multiobjective bilevel programming with birandom variables is formulated for the MCNFP in a large-scale construction project:
(3.15)min⁡C(xj,ykj)=∑(o,d)∈E∑j∈Podcjxj,min⁡T(xj,ykj)=E[∑k∈Ψwk∑j∈Ω∑i∈Φγijt¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]],s.t.{∑j∈Podxj=Qod,(o,d)∈E,∑k∈Ψwk=1,xj≥0,j∈Ω,min⁡C(ykj)=E[∑k∈Ψ∑j∈Ω∑i∈Φγije¯̃ik⌈ykjvk⌉],s.t.{Pr{ω∣Pr{∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]≤Tj}≥θ}≥δ,j∈Ω,k∈Ψ,∑j∈Ω∑k∈Ψγij⌈ykjvk⌉≤ri,i∈Φ,γij={1,ifi∈Aj,j∈Ω,0,otherwise,∑k∈Ψykj≥xj,j∈Ω,ykj≥0,j∈Ω,k∈Ψ.


## 3.4.1. Lower Level Model for the Bilevel MCNFP

The problem posed on the lower level is how to make decisions on material flow by each type of carrier on transportation pathykj while satisfying all capacity constraints, with the main objective being to minimize expected total transportation cost.(1) Objective Functions of the Lower Level
The total transportation cost of the material is calculated by taking the sum of the carriers of each arc's transportation cost and the number of carriers needed to transport the material across the network. In real conditions, it is desirable that each service carrier be fully loaded, so the numbers of carrierk through path j can be denoted by ⌈ykj/vk⌉. Since e¯̃ik is considered a birandom variable, the total transportation cost of material is considered under a birandom environment. Generally, it is difficult to completely minimize total transportation costs because of the birandom variables. Because decision-makers expect minimal cost, the expected value of the total transportation cost is the objective of the lower level. Denote the expected total transportation cost of material as C(ykj); then the objective function of the lower level model can be formulated as
(3.3)min⁡C(ykj)=E[∑k∈Ψ∑j∈Ω∑i∈Φγije¯̃ik⌈ykjvk⌉].
Generally, a path can be represented by a sequence of adjacent arcs. A binary variableγij is introduced to determine whether an arc i is a segment of path j for carrier k:
(3.4)γij={1,ifi∈Aj,j∈Ω,0,otherwise.(2) Constraints of the Lower Level
For transportation time, each carrier requires the transport of material from the source to the destination on scheduleTj. If not, a delay cost is applied. ∑i∈Φγijt¯̃ik represents the total travel time of carrier k on transportation path j, in which t¯̃ik is usually represented by a non-decreasing function (i.e., Bureau of Public Roads (BPR) function) [33] as follows:
(3.5)t¯̃ik=t¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β],i∈Φ,k∈Ψ,
where α and β are user-defined parameters and, in this problem, are set to 0.15 and 2.0, respectively.
Technically, it is not possible to strictly ensure that the random event∑i∈Φγijt¯̃ik does not exceed Tj because of the birandom variable t¯̃i0k. In practical problem, the decision makers often provide an appropriate budget Tj in advance, to ensure that the restriction is, to a certain extent, satisfied, that is to maximize the probability of the random event Pr{∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β]≤Tj}under a given confidence level, which can be written as follows:
(3.6)Pr{ω∣Pr{∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]≤Tj}≥θ}≥δ,j∈Ω,k∈Ψ.
Here the decision makers' aspiration level is indicated as θ, so we use a “Pr” to ensure that the constraint holds at the predetermined confidence level. Additionally, based on probability theory, a further “Pr” is needed to describe the random elements presented in Section 2, which guarantee the establishment of a certain confidence level δ, resembling the P-model (probability maximization model) presented in [34].
The transportation flow may exceed some arcs' capacity because of uncertainties such as the condition of the construction project road. Such conditions may require the manager to select another path. Thus, the total amount of capacity on arci cannot exceed the maximal capacity of the arc i, which produces the following constraint:(3.7)∑j∈Ω∑k∈Ψγij⌈ykjvk⌉≤ri,i∈Φ.
Actually, it is difficult to ensure that each service carrier is fully loaded; so the sum of all flows transported by all kinds of carriers through each path cannot be less than the material flow assigned to it; thus the following constraint is obtained:(3.8)∑k∈Ψykj≥xj,j∈Ω.
The path flow in pathj used by carrier k should not be negative, such that
(3.9)ykj≥0,j∈Ω,k∈Ψ.


## 3.4.2. Upper Level Model for the Bilevel MCNFP

The problem the construction contractor on the upper level faces is how to assign material flow among the transportation paths across the complete transportation network, that is, how to decide the material flowxj through transportation path j. Thus, the decision variable on the upper level is xj.(1) Objective Functions of the Upper Level
For large-scale construction projects, cost and time control are both important, so minimizing total direct costs and total transportation time costs is the two objectives of the upper level model. The two objectives of the upper level can be described as follows.
Firstly, the upper level decision maker attempts to minimize the direct costs of the complete network by assigning the flow of the material to each transportation path to achieve a system optimized flow pattern; thus, the total direct cost is the sum of all transportation costs from different transportation paths. The first objective function of the upper level model can be formulated as follows(3.10)C(xj,ykj)=∑(o,d)∈E∑j∈Podcjxj.
In real conditions, there is increasing pressure to shorter transportation time to reduce or eliminate extra project expenses, with the early arrival of materials shortening the completion time of the construction project and improving construction efficiency. Thus, the total transportation time for carrierk in each path can be described as ∑i∈Φγijt¯̃ik, and, since t¯̃ik is a birandom variable, ∑i∈Φγijt¯̃ik can be regarded as a special birandom variable. Similarly, the expected value of the total transportation time cost is one of the objectives on the upper level. Different carriers are given different weights, and t¯̃ik=t¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β], so the second objective function of the upper level can be described as follows:
(3.11)T(xj,ykj)=E[∑k∈Ψwk∑j∈Ω∑i∈Φγijt¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]].(2) Constraints of the Upper Level
According to the basic assumptions in Section3.1, it is stipulated that the demands between all OD pairs should be satisfied. Thus,
(3.12)∑j∈Podxj=Qod,(o,d)∈E.
The following constraint ensures that the sum of the weights is equal to 1:(3.13)∑k∈Ψwk=1.
In order to describe the nonnegative variables, the constraints in (3.14) are presented:
(3.14)xj≥0,j∈Ω.


## 3.4.3. Global Model for the Bilevel MCNFP

Based on the previous discussion, by integrating (3.3)–(3.14), the following global model for the nonlinear multiobjective bilevel programming with birandom variables is formulated for the MCNFP in a large-scale construction project:
(3.15)min⁡C(xj,ykj)=∑(o,d)∈E∑j∈Podcjxj,min⁡T(xj,ykj)=E[∑k∈Ψwk∑j∈Ω∑i∈Φγijt¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]],s.t.{∑j∈Podxj=Qod,(o,d)∈E,∑k∈Ψwk=1,xj≥0,j∈Ω,min⁡C(ykj)=E[∑k∈Ψ∑j∈Ω∑i∈Φγije¯̃ik⌈ykjvk⌉],s.t.{Pr{ω∣Pr{∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]≤Tj}≥θ}≥δ,j∈Ω,k∈Ψ,∑j∈Ω∑k∈Ψγij⌈ykjvk⌉≤ri,i∈Φ,γij={1,ifi∈Aj,j∈Ω,0,otherwise,∑k∈Ψykj≥xj,j∈Ω,ykj≥0,j∈Ω,k∈Ψ.


## 4. Equivalent Crisp Model

Problem (3.15) is an expectation multiobjective bilevel programming model with chance constraints which has clear nonlinear objectives and constraints. However, this is also difficult to solve with complete confidence because the measures E{·} and Pr{·} are difficult to obtain. In the following section, the normal birandom variables are focused on and the equivalent crisp model of (3.15) is presented.Definition 4.1 (see [35]).
Letξ be a random variable on the probability space (Ω,𝒜,Pr). Then the expected value of ξ is defined by
(4.1)E[ξ]=∫0+∞Pr{ξ≥r}dr-∫-∞0Pr{ξ≤r}dr.
Note that the terms expected value, expectation, and mean value can be used interchangeably.Definition 4.2 (see [36]).
Letξ be a birandom variable on the probability space (Ω,𝒜,Pr); then the expected value of birandom variable ξ can be defined as follows:
(4.2)E[ξ]=∫0+∞Pr{ω∈Ω∣E[ξ(ω)]≥t}dt-∫-∞0Pr{ω∈Ω∣E[ξ(ω)]≤t}dt,
provided that at least one of the aforementioned two integrals is finite.Lemma 4.3 (see [32]).
Letξ be a birandom variable on the probability space (Ω,𝒜,Pr). If the expected value E[ξ(ω)] of the random variable ξ(ω) is finite for each ω, then E[ξ(ω)] is a random variable on (Ω,𝒜,Pr).Remark 4.4.
It should be noted that the expected value operatorE, which appears on both sides of the previous definition of E[ξ], is overloaded. In fact, symbol E represents different meanings. That is to say, overloading allows us to use the same symbol E for different expected value operators, because we can deduce the meaning from the type of argument.Lemma 4.5.
Lett¯̃i0k(ω) be a normal birandom variable, subject to a normal distribution 𝒩(μ̃ti0k(ω),σ1ti0k2(ω)) where μ̃ti0k~𝒩(μti0k,σ2ti0k2). Then the objective function in (3.15)
(4.3)T(xj,ykj)=E[∑k∈Ψwk∑j∈Ω∑i∈Φγijt¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]]
is transformed into the following equivalent objective function:
(4.4)T(xj,ykj)=∑k∈Ψwk∑j∈Ω∑i∈Φγijμti0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β].Proof.
From the assumption, for anyω∈Ω, t¯̃i0k(ω)~𝒩(μ̃ti0k(ω),σ1ti0k2(ω)) is an independent random variable, where μ̃ti0k~𝒩(μti0k,σ2ti0k2). From Definition 4.2,
(4.5)E[t¯̃i0k]=∫0+∞Pr{ω∈Ω∣E[t¯̃i0k(ω)]≥t}dt-∫-∞0Pr{ω∈Ω∣E[t¯̃i0k(ω)]≤t}dt,
since t¯̃1i0k(ω)~𝒩(μ̃ti0k(ω),σti0k2(ω)), and by Definition 4.1, function (4.5) is transformed as follows:
(4.6)E[t¯̃i0k]=∫0+∞Pr{μ̃ti0k(ω)≥t}dt-∫-∞0Pr{μ̃ti0k(ω)≤t}dt=μti0k,
and then
(4.7)T(xj,ykj)=E[∑k∈Ψwk∑j∈Ω∑i∈Φγijt¯̃i0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]]=T(xj,ykj)=∑k∈Ψwk∑j∈Ω∑i∈ΦγijE[t¯̃i0k][1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]=∑k∈Ψwk∑j∈Ω∑i∈Φγijμti0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β].
This completes the proof.Lete¯̃ik(ω) be a normal birandom variable subject to a normal distribution 𝒩(μ̃eik(ω),σ1eik2(ω)) where μ̃eik~𝒩(μeik,σ2eik2). Then similarly the objective function of the lower level is transformed into crisp equivalents:(4.8)C(ykj)=∑k∈Ψ∑j∈Ω∑i∈Φγijμeik⌈ykjvk⌉.Lemma 4.6.
Assume thatt¯̃i0k(ω)~𝒩(μ̃ti0k(ω),σ1ti0k2(ω)) is a normal birandom variable, where μ̃ti0k is a normal distributed random variable characterized by μ̃ti0k~𝒩(μti0k,σ2ti0k2); then ∑i∈Φγijμ̃ti0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β]-Tj∈𝒩(μr,σr2) is also a random variable, where
(4.9)μr=∑i∈Φγijμti0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]-Tj,σr=∑i∈Φ(γij)2σ2ti0k2[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]2.
Then the following constraint for the first constraint is derived in (3.15):
(4.10)Pr{ω∣Pr{∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]≤Tj}≥θ}≥δ,
being equivalent to the following equation:
(4.11)Φ-1(θ)∑i∈Φ(γij)2σ1ti0k2(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]2+μr+σrΦ-1(δ)≤0.Proof.
From the assumption, it is known that for anyω∈Ω, t¯̃i0k(ω)~𝒩(μ̃ti0k(ω),σ1ti0k2(ω)) is an independent random variable, so it follows that
(4.12)∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]-Tj~N(μ̃r(ω),σ̃r2(ω))
is also a normally distributed random variable, where
(4.13)μ̃r(ω)=∑i∈Φγijμ̃ti0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]-Tj,σ̃r2(ω)=∑i∈Φ(γij)2σ1ti0k2(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]2.
Then
(4.14)Pr{∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]≤Tj}≥θ⟺Pr{(∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β]-Tj)-μ̃r(ω)σr(ω)≤-μ̃r(ω)σ̃r(ω)}≥θ⟺Φ(-μ̃r(ω)σ̃r(ω))≥θ⟺μ̃r(ω)≤-Φ-1(θ1r1)σ̃r(ω)⟺∑i∈Φγijμ̃ti0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]-Tj≤-Φ-1(θ)∑i∈Φ(γij)2σ1ti0k2(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]2.
Sinceμ̃ti0k~𝒩(μti0k,σ2ti0k2), then ∑i∈Φγijμ̃ti0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β]-Tj∈𝒩(μr,σr2), where
(4.15)μr=∑i∈Φγijμti0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]-Tj,σr=∑i∈Φ(γij)2σ2ti0k2[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]2,
for given confidence levels θ,δ∈[0,1], so,
(4.16)Pr{ω∣Pr{∑i∈Φγijt¯̃i0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]≤Tj}≥θ}≥δ⟺Pr{ω∣∑i∈Φγijμ̃ti0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]-Tj≤M}≥δ⟺Pr{ω∣∑i∈Φγijμ̃ti0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β]-Tj-μrσr≤M-μrσr}≤δ⟺M≥μr+σrΦ-1(δ)⟺Φ-1(θ)∑i∈Φ(γij)2σ1ti0k2(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]2+μr+σrΦ-1(δ)≤0,
where M=-Φ-1(θ)∑i∈Φ(γij)2σ1ti0k2(ω)∑[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β]2.
The proof is completed.Based on Lemmas4.3 and 4.5, it is determined that the expectation multiobjective bilevel programming model with chance constraints (3.15) is equivalent to the following crisp nonlinear multiobjective bilevel programming problem:
(4.17)min⁡C(xj,ykj)=∑(o,d)∈E∑j∈Podcjxj,min⁡T(xj,ykj)=∑k∈Ψwk∑j∈Ω∑i∈Φγijμti0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β],s.t.{∑j∈Podxj=Qod,(o,d)∈E,∑k∈Ψwk=1,xj≥0,j∈Ω,min⁡C(ykj)=∑k∈Ψ∑j∈Ω∑i∈Φγijμeik⌈ykjvk⌉,s.t{Φ-1(θ)∑i∈Aj(γij)2σ1ti0k2(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]2+μr+σrΦ-1(δ)≤0,j∈Ω,k∈Ψ,∑j∈Ω∑k∈Ψγij⌈ykjvk⌉≤ri,i∈Φ,γij={1,ifi∈Aj,i∈Ω,0,otherwise,∑k∈Ψykj≥xj,j∈Ω,ykj≥0,j∈Ω,k∈Ψ,
where μr=∑i∈Φγijμti0k(ω)[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β]-Tj,σr=∑i∈Φ(γij)2σ2ti0k2[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉/ri)β]2 and  θ and δ are predicted confidence levels specified by the DM.When (4.17) is used to make decisions, the upper level decision maker first makes a decision xj, then the lower decision maker reacts ykj according to the cost, and the upper level decision maker makes a proper adjustment based on the lower level feedback, finally making the upper level objective optimal. Thus, the upper programming and the lower programming influence and restrict each other in bilevel programming.


## 5. Solution Approach

To obtain an analytical optimal solution for a bilevel programming problem (BLPP) is difficult, yet there is theoretical evidence supporting these observations since BLPP is in fact NP-hard even in its linear form [37]. Moreover, this problem is nonlinear and nondifferentiable, and the MCNFP in a large-scale construction project has various nodes and links. On the other hand, the nondifferentiable piecewise objective functions and constraints presented in the MCNFP bring computational difficulties. Thus, the possibility of finding a solution to the complexity is increased, and it is difficult to solve with any exact algorithm. Therefore, here, an MOBLPSO is proposed to solve the MCNFP in a large-scale construction project. Because particle swarm optimization (PSO) is computationally tractable compared with other heuristic algorithms, it is easy to implement, and does not require any gradient information for an objective function but the value.PSO is a population-based self-adaptive search stochastic optimization technique proposed by Kennedy and Eberhart [38], which was inspired by the social behaviour of animals such as fish schooling and birds flocking. Similar to other population-based algorithms, such as evolutionary algorithms, PSO can solve a variety of difficult optimization problems but has shown a faster convergence rate than other evolutionary algorithms on some problems [39]. PSO is influenced little by the continuity of the objective function; it just uses the primary math operators and receives good results in static, noisy, and continuously changing environments [40]. Another advantage of PSO is that it has very few parameters to adjust, which makes it particularly easy to implement. Since PSO can be implemented easily and effectively, it has been applied in solving real-world optimization problems in recent years, such as [41–45]. In PSO, the following formulas [38] are applied to update the position and velocity of each particle:
(5.1)vdl(τ1+1)=w(τ)vdl(τ1)+cpr1[pdl,best(τ)-pdl(τ)]+cgr2[gdl,best(τ)-pdl(τ)],pdl(τ1+1)=pdl(τ)+vdl(τ+1),
where vdl is the velocity of lth particle at the dth dimension in the τth iteration, w(τ) is the inertia weight, pdl(τ) is the position of lth particle at the dth dimension, r1 and r2 are random numbers in the range [0,1], cp and cg are personal and global best position acceleration constants, respectively, while, pdl,best(τ) is personal best position of the lth particle at the dth dimension, and gdl,best(τ) is the global best position at the dth dimension.As mentioned before, it is very difficult to solve the bilevel model, especially when the model is nonlinear. The contribution of this paper is that a universal effective algorithm for solving the bilevel model is put forward, which is based on the hierarchical iteration. The key idea of the algorithm is that the optimum of the bilevel model can be approached step by step through repeatedly iterative calculations between the upper and lower models. The main body of the proposed approach is a type of PSO—multiobjective PSO (MOPSO) with Pareto-Archived Evolution Strategy (PAES)—which is designed only to cope with the upper level programming problem based on the follower's optimal response. Another type of improved PSO—PSO with passive congregation (PSOPC)—is embedded to deal with the lower level programming problem and obtain the optimal response of the follower for each given decision variable of the upper level programming. The follower's optimal reaction is then returned to upper level programming problem as the implementation base for the MOPSO. The algorithm is called the MOBLPSO, the notations of which for the MCNFP are listed as follows:τ1,τ2: iteration index of the upper and lower level, τ1=1,2,…,T1 and τ2=1,2,…,T2;l1,l2: particle index of the upper and lower level, l1=1,2,…,L1 and l2=1,2,…,L2;j: index of transportation path, j=1,2,…,J;k: index of carrier, k=1,2,…,K;r1,r2,r3,r4,r5: uniform distributed random number within [0,1];w1(τ1): inertia weight in the τ1th iteration of the upper level;w2(τ2): inertia weight in the τ2th iteration of the lower level;Rkjl2(τ2): particle selected randomly from the swarm at the kjth dimension in the τ2th iteration;c1p, c1g: personal and global best position acceleration constant of the upper level;c2p, c2g: personal and global best position acceleration constant of the lower level;cpc: passive congregation coefficient of the lower level;p1max,p1min: maximum and minimum position value of the upper level;Pl1(τ1): vector position of the l1th particle in the τ1th iteration, Pl1(τ1)=[p1l1(τ1),p2l1(τ1),…,pJl1(τ1)];Pl2(τ2): vector position of the l2th particle in the τ2th iteration, Pl2(τ2)=[p1l2(τ2),p2l2(τ2),…,pKJl1(τ2)];Vl1(τ1): vector velocity of the l1th particle in the τ1th iteration, Vl1(τ1)=[v1l1(τ1),v2l1(τ1),…,vJl1(τ1)];Vl2(τ2): vector velocity of the l2th particle in the τ2th iteration, Vl2(τ2)=[v1l2(τ2),v2l2(τ2),…,vKJl2(τ2)];Pl1,best(τ1): vector personal best position of the l1th particle in the τ1th iteration, Pl1,best(τ1)=[p1l1,best(τ1),p2l1,best(τ1),…,pJl1,best(τ1)];Pl2,best(τ1): vector personal best position of the l2th particle in the τ2th iteration, Pl2,best(τ2)=[p1l2,best(τ2),p2l2,best(τ2),…,pJl2,best(τ2)];Gbest(τ2): vector global best position in the τ2th iteration, Gbest(τ2)=[g1best(τ2),g2best(τ2),…,gKJbest(τ2)];ARC: the positions of the particles that represent nondominated vectors in the repository;F(Pl1(τ1)): fitness value of Pl1(τ1); andF(Pl2(τ2)): fitness value of Pl2(τ2).


### 5.1. Multiobjective Methods for the Upper Level Programming

Researchers are also seeing PSO as a very strong competitor to other algorithms in solving multiobjective optimal problems [46] and it has been proved to be especially suitable for multiobjective optimization [47]. The method applied here to deal with the upper level problem is a multiobjective method which combines a MOPSO with PAES. The PAES [48] is one of the Pareto-based approaches to update the best position. The methods use a truncated archive to store the elite individuals. This approach uses leader selection techniques based on Pareto dominance. The basic idea is to select as leaders to the particles that are nondominated with respect to the swarm.The details for the PAES procedure, test procedure, and selection procedure are stated hereinafter (Algorithms1 and 2), in which Pl1,best(τ1) is initially set equal to the initial position of particle l1.Algorithm 1
Procedure The updating of the best position in the subsequent iterationsgenerate an initial random solutionpjl1,best(τ1) and add it to the archive;updatepjl1,best(τ1) to produce pjl1(τ1+1) and evaluate pjl1(τ1+1);If(pjl1,best(τ1)dominatespjl1(τ1+1))  then discard pjl1(τ1+1);else if(pjl1(τ1+1)dominatespjl1,best(τ1))  then replace pjl1,best(τ1) with pjl1(τ1+1), and add pjl1(τ1+1)to the archive pjl1,best(τ1+1)=pjl1(τ1+1);else if(pjl1(τ1+1)isdominatedbyanymemberofthearchive) then discard pjl(τ1+1);else apply test (pjl1,best(τ1),pjl1(τ1+1),archive) to determine which becomes the new currentsolution and whether to addpjl1(τ1+1) to the archive;until a termination criterion has been reached, return to line 2.Algorithm 2
Procedure Testif the archive is not full;addpjl1(τ1+1) to the archive;if(pjl1(τ1+1)isinalesscrowdedregionofthearchivethanpjl1,best(τ1));acceptpjl1,best(τ1+1)=pjl1(τ1+1);else maintain pjl1,best(τ1+1)=pjl1,best(τ1);elseIf(pjl1(τ1+1) is in a less crowded region of the archive than pjl1,best(τ1)forsomememberpjl1,best(τ1)on the archive);addpjl1(τ1+1) to the archive, and remove a member of the archive from the mostcrowded region;if(pjl1(τ1+1)isinalesscrowdedregionofthearchivethanpjl1,best(τ1));acceptpjl1,best(τ1+1)=pjl1(τ1+1);else maintain pjl1,best(τ1+1)=pjl1,best(τ1);elseif(pjl1(τ1+1)isinalesscrowdedregionofthearchivethanpjl1,best(τ1));acceptpjl1,best(τ1+1)=pjl1(τ1+1);else maintain pjl1,best(τ1+1)=pjl1,best(τ1).The repository that stores the positions of the particles that represent nondominated vectors is denoted by ARC; then the velocity of eachl1th particle of the upper level is updated using the following equations:
(5.2)vjl1(τ1+1)=w1(τ1)vjl1(τ1)+c1pr1[pjl1,best(τ1)-pjl1(τ1)]+c1gr2[ARCh(τ1)-pjl1(τ1)],
where ARCh(τ1) is a solution randomly selected from the repository in iteration τ1, which can improve significantly the ability of global convergence by avoiding being trapped in a stagnant state in finite iterations [49, 50]. The index h is selected in the following way: the hypercubes containing more than one particle are assigned a fitness equal to the result of dividing any number into the number of particles they contain. This aims to decrease the fitness of those hypercubes that contain more particles which is seen as a form of fitness sharing [51]. Then, a roulette-wheel selection is applied using these fitness values to select the hypercube from which the corresponding particle is taken. Once the hypercube has been selected, a particle is selected randomly within the hypercube.


### 5.2. PSOPC for the Lower Level Programming

From the previous description, because there are many constraints in the lower level, the standard PSO can easily fall into premature convergence, so a PSOPC based on the standard PSO is adopted to solve the lower level problem. He et al. [52] proved that PSOPC can avoid the premature convergence problem, in the running of the PSO algorithm by adding a passive congregation coefficient to the standard PSO, as this helps the algorithm move out of the local optimal solution improving the convergence of the algorithm, and thus improving the global search ability. The following equation is adopted to update the velocity of each l2th particle of the lower level:(5.3)vkjl2(τ2+1)=w2(τ2)vkjl2(τ2)+c2pr3[pkjl2,best(τ2)-pkjl2(τ2)]+c2gr4[gkjbest(τ2)-pkjl2(τ2)]+cpcr5(Rkjl2(τ2)-pkjl2(τ2)).


### 5.3. Framework of the Proposed MOBLPSO for the MCNFP

In our proposed algorithm, solving multiobjective bilevel programming is transformed to solve the upper level and lower level programming problem interactively while determining, respectively, the decision variable of the upper level or the lower level. To be mentionable, the PSOPC for the lower level is nested in the MOPSO for the upper level, and the MOPSO is the main body of the MOBLPSO. The solution information is exchanged between the two types of PSO, and the output of one algorithm is the input of another algorithm, namely,y, the output of PSOPC for lower level programming is the input of the MOPSO for upper level programming; and x, the output of the MOPSO for the upper level programming is the input of PSOPC for the lower level programming. These form a hierarchical and sequential framework.Parameter Selection
In order to guarantee the convergence of MOBLPSO, the parameters are selected on the basis of empirical results that are carried out to observe the behaviour of the algorithm in different parameter settings. By comparing several sets of parameters, including population size, iteration number, acceleration coefficients, and inertia weight, the empirical results have shown that the constant acceleration coefficients withc1p=c1g=1.5 for the upper level, c2p=1, c2g=2, and cpc=1.5 (i.e., passive congregation coefficient [52]) for the lower level, and the adaptive inertia weights [53] provide good convergent behaviour in this study, which is in accordance with the results provided by Eberhart and Shi [54]. The adaptive inertia weights for the upper level (i.e., e=1) and lower level (i.e., e=2) are set to be varying with iteration as follows:
(5.4)we(τe)=we(Te)-τe-Te1-Te[we(1)-we(Te)],
where the iteration numbers are T1=500 (i.e., for the upper level) and T2=100 (i.e., for the lower level), and we(1)=0.9 and we(Te)=0.1 (for e=1,2). Since the probability of becoming trapped in the stagnant state can be reduced dramatically by using a large number of particles [49], the population sizes are set to be 300 for the upper level and 100 for the lower level.Initialize
In the upper level, set iterationτ1=1. For l1=1,2,…,L1, generate the position of the particle l1 with an integer random position (note that every particle in the upper level consists of j dimensions in this study). In the lower level, set iteration τ2=1. For l2=1,2,…,L2, generate the position of the particle l2 with an integer random position (note that every particle in the lower level consists of k×j dimensions in this study).Decode Particles into Solutions
Decode particles into solutions: forl1=1,2,…,L1, decode pl1(τ1) to a solution as pjl1(τ1)=xj(τ1). For l2=1,2,…,L2, decode pl2(τ2) to a solution as pkjl2(τ2)=ykj(τ2). Mapping between one potential solution for the upper and lower level of the MCNFP and particle representation is shown in Figure 4.Figure 4
Decoding method and mapping between PSO particles and solutions of two levels.Check the Feasibility
All particles of the upper level satisfy the constraints of the upper level. All particles of the lower level satisfy the constraints of the lower level. Then, the particles are feasible.Fitness Value
The fitness value used to evaluate the particle is the value of objective function in each level. There are two objectives in the upper level; particlePl1(τ1)'s fitness value F(Pl1(τ1)) is a 1×2 matrix, namely,
(5.5)F(Pl1(τ1))=[∑(o,d)∈E∑j∈Podcjxj,∑k∈Ψwk∑j∈Ω∑i∈Φγijμti0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]].
The fitness value of the particle in the lower level is
(5.6)F(Pl2(τ2))=∑k∈Ψ∑j∈Ω∑i∈Φγijμeik⌈ykjvk⌉.
Figure5 shows the schematic procedure for the MOBLPSO to generate solutions for the proposed multiobjective bilevel model. Such a repeated interaction and cooperation between two types of PSO reflects and simulates the decision process of multiobjective bilevel programming and is able to solve multiobjective bilevel programming sequentially.Figure 5
Schematic diagram of MOBLPSO for MCNFP.


## 5.1. Multiobjective Methods for the Upper Level Programming

Researchers are also seeing PSO as a very strong competitor to other algorithms in solving multiobjective optimal problems [46] and it has been proved to be especially suitable for multiobjective optimization [47]. The method applied here to deal with the upper level problem is a multiobjective method which combines a MOPSO with PAES. The PAES [48] is one of the Pareto-based approaches to update the best position. The methods use a truncated archive to store the elite individuals. This approach uses leader selection techniques based on Pareto dominance. The basic idea is to select as leaders to the particles that are nondominated with respect to the swarm.The details for the PAES procedure, test procedure, and selection procedure are stated hereinafter (Algorithms1 and 2), in which Pl1,best(τ1) is initially set equal to the initial position of particle l1.Algorithm 1
Procedure The updating of the best position in the subsequent iterationsgenerate an initial random solutionpjl1,best(τ1) and add it to the archive;updatepjl1,best(τ1) to produce pjl1(τ1+1) and evaluate pjl1(τ1+1);If(pjl1,best(τ1)dominatespjl1(τ1+1))  then discard pjl1(τ1+1);else if(pjl1(τ1+1)dominatespjl1,best(τ1))  then replace pjl1,best(τ1) with pjl1(τ1+1), and add pjl1(τ1+1)to the archive pjl1,best(τ1+1)=pjl1(τ1+1);else if(pjl1(τ1+1)isdominatedbyanymemberofthearchive) then discard pjl(τ1+1);else apply test (pjl1,best(τ1),pjl1(τ1+1),archive) to determine which becomes the new currentsolution and whether to addpjl1(τ1+1) to the archive;until a termination criterion has been reached, return to line 2.Algorithm 2
Procedure Testif the archive is not full;addpjl1(τ1+1) to the archive;if(pjl1(τ1+1)isinalesscrowdedregionofthearchivethanpjl1,best(τ1));acceptpjl1,best(τ1+1)=pjl1(τ1+1);else maintain pjl1,best(τ1+1)=pjl1,best(τ1);elseIf(pjl1(τ1+1) is in a less crowded region of the archive than pjl1,best(τ1)forsomememberpjl1,best(τ1)on the archive);addpjl1(τ1+1) to the archive, and remove a member of the archive from the mostcrowded region;if(pjl1(τ1+1)isinalesscrowdedregionofthearchivethanpjl1,best(τ1));acceptpjl1,best(τ1+1)=pjl1(τ1+1);else maintain pjl1,best(τ1+1)=pjl1,best(τ1);elseif(pjl1(τ1+1)isinalesscrowdedregionofthearchivethanpjl1,best(τ1));acceptpjl1,best(τ1+1)=pjl1(τ1+1);else maintain pjl1,best(τ1+1)=pjl1,best(τ1).The repository that stores the positions of the particles that represent nondominated vectors is denoted by ARC; then the velocity of eachl1th particle of the upper level is updated using the following equations:
(5.2)vjl1(τ1+1)=w1(τ1)vjl1(τ1)+c1pr1[pjl1,best(τ1)-pjl1(τ1)]+c1gr2[ARCh(τ1)-pjl1(τ1)],
where ARCh(τ1) is a solution randomly selected from the repository in iteration τ1, which can improve significantly the ability of global convergence by avoiding being trapped in a stagnant state in finite iterations [49, 50]. The index h is selected in the following way: the hypercubes containing more than one particle are assigned a fitness equal to the result of dividing any number into the number of particles they contain. This aims to decrease the fitness of those hypercubes that contain more particles which is seen as a form of fitness sharing [51]. Then, a roulette-wheel selection is applied using these fitness values to select the hypercube from which the corresponding particle is taken. Once the hypercube has been selected, a particle is selected randomly within the hypercube.


## 5.2. PSOPC for the Lower Level Programming

From the previous description, because there are many constraints in the lower level, the standard PSO can easily fall into premature convergence, so a PSOPC based on the standard PSO is adopted to solve the lower level problem. He et al. [52] proved that PSOPC can avoid the premature convergence problem, in the running of the PSO algorithm by adding a passive congregation coefficient to the standard PSO, as this helps the algorithm move out of the local optimal solution improving the convergence of the algorithm, and thus improving the global search ability. The following equation is adopted to update the velocity of each l2th particle of the lower level:(5.3)vkjl2(τ2+1)=w2(τ2)vkjl2(τ2)+c2pr3[pkjl2,best(τ2)-pkjl2(τ2)]+c2gr4[gkjbest(τ2)-pkjl2(τ2)]+cpcr5(Rkjl2(τ2)-pkjl2(τ2)).


## 5.3. Framework of the Proposed MOBLPSO for the MCNFP

In our proposed algorithm, solving multiobjective bilevel programming is transformed to solve the upper level and lower level programming problem interactively while determining, respectively, the decision variable of the upper level or the lower level. To be mentionable, the PSOPC for the lower level is nested in the MOPSO for the upper level, and the MOPSO is the main body of the MOBLPSO. The solution information is exchanged between the two types of PSO, and the output of one algorithm is the input of another algorithm, namely,y, the output of PSOPC for lower level programming is the input of the MOPSO for upper level programming; and x, the output of the MOPSO for the upper level programming is the input of PSOPC for the lower level programming. These form a hierarchical and sequential framework.Parameter Selection
In order to guarantee the convergence of MOBLPSO, the parameters are selected on the basis of empirical results that are carried out to observe the behaviour of the algorithm in different parameter settings. By comparing several sets of parameters, including population size, iteration number, acceleration coefficients, and inertia weight, the empirical results have shown that the constant acceleration coefficients withc1p=c1g=1.5 for the upper level, c2p=1, c2g=2, and cpc=1.5 (i.e., passive congregation coefficient [52]) for the lower level, and the adaptive inertia weights [53] provide good convergent behaviour in this study, which is in accordance with the results provided by Eberhart and Shi [54]. The adaptive inertia weights for the upper level (i.e., e=1) and lower level (i.e., e=2) are set to be varying with iteration as follows:
(5.4)we(τe)=we(Te)-τe-Te1-Te[we(1)-we(Te)],
where the iteration numbers are T1=500 (i.e., for the upper level) and T2=100 (i.e., for the lower level), and we(1)=0.9 and we(Te)=0.1 (for e=1,2). Since the probability of becoming trapped in the stagnant state can be reduced dramatically by using a large number of particles [49], the population sizes are set to be 300 for the upper level and 100 for the lower level.Initialize
In the upper level, set iterationτ1=1. For l1=1,2,…,L1, generate the position of the particle l1 with an integer random position (note that every particle in the upper level consists of j dimensions in this study). In the lower level, set iteration τ2=1. For l2=1,2,…,L2, generate the position of the particle l2 with an integer random position (note that every particle in the lower level consists of k×j dimensions in this study).Decode Particles into Solutions
Decode particles into solutions: forl1=1,2,…,L1, decode pl1(τ1) to a solution as pjl1(τ1)=xj(τ1). For l2=1,2,…,L2, decode pl2(τ2) to a solution as pkjl2(τ2)=ykj(τ2). Mapping between one potential solution for the upper and lower level of the MCNFP and particle representation is shown in Figure 4.Figure 4
Decoding method and mapping between PSO particles and solutions of two levels.Check the Feasibility
All particles of the upper level satisfy the constraints of the upper level. All particles of the lower level satisfy the constraints of the lower level. Then, the particles are feasible.Fitness Value
The fitness value used to evaluate the particle is the value of objective function in each level. There are two objectives in the upper level; particlePl1(τ1)'s fitness value F(Pl1(τ1)) is a 1×2 matrix, namely,
(5.5)F(Pl1(τ1))=[∑(o,d)∈E∑j∈Podcjxj,∑k∈Ψwk∑j∈Ω∑i∈Φγijμti0k[1+α(∑j∈Ω∑k∈Ψγij⌈ykj/vk⌉ri)β]].
The fitness value of the particle in the lower level is
(5.6)F(Pl2(τ2))=∑k∈Ψ∑j∈Ω∑i∈Φγijμeik⌈ykjvk⌉.
Figure5 shows the schematic procedure for the MOBLPSO to generate solutions for the proposed multiobjective bilevel model. Such a repeated interaction and cooperation between two types of PSO reflects and simulates the decision process of multiobjective bilevel programming and is able to solve multiobjective bilevel programming sequentially.Figure 5
Schematic diagram of MOBLPSO for MCNFP.


## 6. Practical Application

### 6.1. Project Description

In this section, an earth-rock work transportation project in a large-scale water conservancy and hydropower construction project is taken as an example for our optimization method. The Shuibuya Hydropower Project was conducted in Badong County, which is located in the middle reaches of Qingjiang River in Sichuan province, China. The project is the first cascaded project on the Qingjiang main stream and the third important project following Geheyan and Gaobazhou in China. Once completed, it will provide a major power source to meet the peak load demand in the Central China Power Grid. The installed capacity and annual output of Shuibuya Power Plant are 1,600 MW and 3.92 GWh, respectively. The project has a powerful regulating ability with a normal pool level of 400 m and reservoir capacity of4.58×109m3. The project consists of a concrete-faced rock fill dam (CFRD), underground power house, a chute spillway on the left bank, and the sluice tunnel on the right bank. The dam is 233 m high and is the tallest of its kind in the world at present with a total volume of 15.64×106m3.


### 6.2. Data Collection

All detailed data for the Shuibuya Hydropower Project were obtained from the Hubei Qingjiang Shuibuya Project Construction Company. In a large-scale construction project, especially in a large water conservancy and hydropower construction project, earth-rock work is usually a primary material, and earth-rock work transportation occurs every day in excavation projects, borrow areas, filling projects, dumpling sites, and stockpile areas as it is turned over and needs to be replaced frequently (see Figure6). In the Shuibuya Hydropower Project, there are 4 excavation projects, 2 borrow areas, 3 stockpile areas, and 2 dumpling sites, The location and detailed information of borrow areas, dumpling sites, and stockpile areas, of the Shuibuya Hydropower Project is illustrated in Figure 7.Figure 6
Earth-rock work transportation relations in large-scale construction project.Figure 7
Layout of borrow areas, dumpling sites and stockpile areas of Shuibuya Hydropower Project.Three types of dump trucks (carriers) in the construction project are considered, which transport earth-rock work along different paths connecting the OD pairs, with the destination nodes having the practical demand of timeliness. All necessary data for every kind of carrier were calculated as shown in Table1, and Table 2 shows the details of the paths in the whole road network.Table 1
Information of carriers in Shuibuya Hydropower Project.
CarrierKind indexkType (heaped capacity; maximum payload)WeightwkDump truck1K30N-8.4 (16m3; 26 t)0.32Terex TA28 (17m3; 28 t)0.33Perlini DP366 (20m3; 36 t)0.4Table 2
Related data about information of paths in earth-rock work transportation network.
OD pairsDemandQod (m3)Transportation pathjDirect cost (CNY/m3)Constraint time (h)Excavation project 1,2→dumpling site 19501⇔#5→#7→#91.801.45Excavation project 1,2→filling project 1,226502⇔ #5→#3→#1→#62.502.153⇔#5→#1→#61.701.354⇔ #5→#7→#8→#62.402.205⇔ #5→#7→#61.901.40Excavation project 1,2→stockpile area 117006⇔ #5→#71.201.25Excavation project 3,4→filling project 3,423007⇔#11→#10→#14→#162.302.258⇔ #11→#13→#161.901.459⇔ #11→#10→#13→#162.702.3010⇔ #11→#10→#14→#15→#163.003.50Excavation project 3,4→dumpling site 2165011⇔ #11→#13→#121.601.5012⇔ #11→#10→#13→#122.332.50Excavation project 3,4→stockpile area 2148013⇔ #11→#10→#132.601.7514⇔ #11→#131.101.25Excavation project 3,4→stockpile area 3138015⇔ #11→#10→#14→#152.102.33Borrow area 1→filling project 1,2170016⇔ #2→#4→#61.501.50Borrow area 2→filling project 3,4125017⇔ #14→#15→#161.801.5518⇔ #14→#161.201.20Stockpile area 1→filling project 1,2115019⇔ #7→#8→#61.901.45Stockpile area 2→filling project 3,4130020⇔ #13→#161.401.17Stockpile area 3→filling project 1,2126021⇔ #15→#7→#8→#62.402.1722⇔ #15→#7→#62.201.60The transportation network in the project includes an internal road network and an external road network. In this case, only the internal road network is considered. The internal road network is composed of 20 trunk roads located on the left and right banks, with 11 on the left and 9 on the right, with a cross-river bridge connecting the left and right banks. Therefore, 21 links are considered in this paper. In order to apply the proposed methods more conveniently, adjacent roads of the same type have been combined and road shapes have been ignored. An abstracted transportation network is illustrated in Figure8. For each link in the transportation network, there are a free flow birandom travel time t¯̃i0k and birandom transportation cost e¯̃ik. The corresponding data of which are stated in Table 3. In order to collect the data of transportation time and costs, investigations and surveys were made to obtain historical data from the financial department and experienced engineers of construction team in Hubei Qingjiang Shuibuya Project Construction Company. Since the transportation time and costs for each arc of the path change over time, the data are classified into different parts based on different periods. Both transportation time and costs are assumed to approximately follow normal distributions for each period, and the two parameters (expected value and variance) for the normal distributions are estimated by using maximum likelihood estimation, which justified by a chi-square goodness-of-fit test. By comparing the normal distributions for the same transportation time and costs in different periods, it can be found that the expected values of the aforementioned normal distributions also approximately follow a similar random distribution pattern, which has also been justified using a chi-square goodness-of-fit test. It should be noted that since the variance fluctuations are quite insignificant, the median values of the variances in different periods are selected as the variance for the previous normal distribution. The predicted confidence levels given by the decision maker are, respectively, θ=0.9 and δ=0.85.Table 3
Birandom free flow travel timet¯̃i0k and birandom transportation cost e¯̃ik.
ArciCorresponding nodesArc capacityri (n/h)Transportation costeik (CNY/unit)Free transportation timeti0k (h)Carrier 1Carrier 2Carrier 3Carrier 1Carrier 2Carrier 31#1, #3105𝒩(μ, 0.64), withμ~𝒩 (5.2, 0.10)𝒩(μ, 1.00), withμ~𝒩 (6.0, 0.09)𝒩(μ, 0.01), withμ~𝒩 (6.8, 1.00)𝒩(μ, 0.25), withμ~𝒩 (0.2, 0.01)𝒩(μ, 0.81),withμ~𝒩 (0.24, 0.09)𝒩(μ, 0.16), withμ~𝒩(0.28, 0.04)2#1, #5110𝒩(μ, 1.00), withμ~𝒩(3.2, 0.64)𝒩(μ, 0.49), withμ~𝒩 (3.5, 0.16)𝒩(μ, 0.81), withμ~𝒩 (3.8, 0.25)𝒩(μ, 0.36), withμ~𝒩 (0.15, 0.04)𝒩(μ, 0.25), withμ~𝒩 (0.16, 0.04)𝒩(μ, 0.09), withμ~𝒩 (0.18, 0.01)3#1, #6135𝒩(μ, 4.00), withμ~𝒩 (7.2, 1.00)𝒩(μ, 0.81), withμ~𝒩 (7.5, 0.64)𝒩(μ, 0.64), withμ~𝒩(7.8, 0.49)𝒩(μ, 0.25), withμ~𝒩 (0.30, 0.16)𝒩(μ, 0.36), withμ~𝒩 (0.35, 0.09)𝒩(μ, 0.16), withμ~𝒩 (0.38, 0.01)4#2, #4115𝒩(μ, 1.00), withμ~𝒩 (3.7, 0.49)𝒩(μ, 4.00), withμ~𝒩 (4.0, 0.81)𝒩(μ, 0.01), withμ~𝒩 (4.3, 0.01)𝒩(μ, 0.81), withμ~𝒩 (0.14, 0.16)𝒩(μ, 0.36), withμ~𝒩 (0.18, 0.09)𝒩(μ, 0.64), withμ~𝒩 (0.19, 0.01)5#3, #5108𝒩(μ, 0.36), withμ~𝒩 (2.2, 0.01)𝒩(μ, 0.81), withμ~𝒩 (2.4, 0.09)𝒩(μ, 0.64), withμ~𝒩 (2.6, 0.16)𝒩(μ, 0.01), withμ~𝒩 (0.14, 0.01)𝒩(μ, 0.25), withμ~𝒩 (0.18, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.22, 0.01)6#4, #6112𝒩(μ, 0.81), withμ~𝒩 (3.2, 0.16)𝒩(μ, 1.00), withμ~𝒩 (4.0, 0.09)𝒩(μ, 0.49), withμ~𝒩 (4.8, 0.09)𝒩(μ, 0.64), withμ~𝒩 (0.12, 0.16)𝒩(μ, 0.09), withμ~𝒩 (0.15, 0.01)𝒩(μ, 0.25), withμ~𝒩 (0.18, 0.01)7#5, #7180𝒩(μ, 1.00), withμ~𝒩(5.0, 0.01)𝒩(μ, 0.64), withμ~𝒩 (5.2, 0.25)𝒩(μ, 0.36), withμ~𝒩 (5.4, 0.16)𝒩(μ, 0.01), withμ~𝒩 (0.25, 0.01)𝒩(μ, 0.04), withμ~𝒩 (0.28, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.32, 0.01)8#6, #7145𝒩(μ, 0.49), withμ~𝒩 (4.5, 0.01)𝒩(μ, 0.36), withμ~𝒩 (4.2, 0.16)𝒩(μ, 0.49), withμ~𝒩 (3.4, 0.16)𝒩(μ, 0.09), withμ~𝒩 (0.16, 0.01)𝒩(μ, 0.04), withμ~𝒩 (0.19, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.23, 0.01)9#6, #8160𝒩(μ, 0.16), withμ~𝒩 (3.6, 0.25)𝒩(μ, 0.81), withμ~𝒩 (4.0, 1.00)𝒩(μ, 1.00), withμ~𝒩 (4.4, 0.09)𝒩μ, 0.04), withμ~𝒩 (0.13, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.15, 0.04)𝒩(μ, 0.09), withμ~𝒩 (0.17, 0.01)10#7, #8165𝒩(μ, 1.00), withμ~𝒩 (4.0, 0.81)𝒩(μ, 1.00), withμ~𝒩 (4.4, 0.49)𝒩(μ, 0.64), withμ~𝒩 (4.8, 0.25)𝒩(μ, 0.01), withμ~𝒩 (0.16, 0.01)𝒩(μ, 0.04), withμ~𝒩 (0.19, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.23, 0.04)11#7, #9100𝒩(μ, 4.00), withμ~𝒩 (7.2, 0.09)𝒩(μ, 0.16), withμ~𝒩 (7.6, 0.81)𝒩(μ, 0.25), withμ~𝒩 (8.0, 0.16)𝒩(μ, 0.01), withμ~𝒩 (0.36, 0.01)𝒩(μ, 0.04), withμ~𝒩 (0.39,0.01)𝒩(μ, 0.09), withμ~𝒩 (0.42, 0.01)12#7, #15130𝒩(μ, 0.81), withμ~𝒩 (8.0, 0.81)𝒩(μ, 0.64), withμ~𝒩 (8.5, 0.50)𝒩(μ, 1.00), withμ~𝒩 (8.7, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.33, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.36, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.38, 0.01)13#10, #11225𝒩(μ, 4.00), withμ~𝒩 (5.0, 0.81)𝒩(μ, 0.64), withμ~𝒩 (6.0, 0.16)𝒩(μ, 0.25), withμ~𝒩 (6.5, 0.09)𝒩(μ, 0.01), withμ~𝒩 (0.30, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.34, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.36, 0.01)14#10, #13155𝒩(μ, 0.81), withμ~𝒩 (5.5, 0.36)𝒩(μ, 1.00), withμ~𝒩 (5.6, 0.49)𝒩(μ, 0.04), withμ~𝒩 (5.8, 0.16)𝒩(μ, 0.01), withμ~𝒩 (0.22, 0.09)𝒩(μ, 0.04), withμ~𝒩 (0.26, 0.16)𝒩(μ, 0.04), withμ~𝒩 (0.30, 0.09)15#10, #14150𝒩(μ, 0.64), withμ~𝒩 (4.5, 0.49)𝒩(μ, 1.00), withμ~𝒩 (4.6, 0.25)𝒩(μ, 0.09), withμ~𝒩 (4.8, 0.09)𝒩(μ, 0.04), withμ~𝒩 (0.22, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.25, 0.04)𝒩(μ, 0.01), withμ~𝒩 (0.27, 0.01)16#11, #13170𝒩(μ, 0.64), withμ~𝒩 (3.6, 0.16)𝒩(μ, 0.09), withμ~𝒩 (4.0, 0.01)𝒩(μ, 0.16), withμ~𝒩 (4.4, 0.25)𝒩(μ, 0.04), withμ~𝒩 (0.16, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.18, 0.04)𝒩(μ, 0.01), withμ~𝒩 (0.20, 0.01)17#12, #13125𝒩(μ, 0.25), withμ~𝒩 (4.2, 0.16)𝒩(μ, 0.49), withμ~𝒩 (4.3, 0.09)𝒩(μ, 1.00), withμ~𝒩(4.4, 0.16)𝒩(μ, 0.01), withμ~𝒩 (0.22, 0.01)𝒩(μ, 0.04), withμ~𝒩 (0.24, 0.09)𝒩(μ, 0.09), withμ~𝒩 (0.26, 0.01)18#13, #16145𝒩(μ, 1.00), withμ~𝒩 (4.5, 0.16)𝒩(μ, 0.04), withμ~𝒩 (4.7, 0.81)𝒩(μ, 0.49), withμ~𝒩 (4.9, 1.00)𝒩(μ, 0.01), withμ~𝒩 (0.21, 0.64)𝒩(μ, 1.00), withμ~𝒩 (0.23, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.25, 0.16)19#14, #15140𝒩(μ, 0.25), withμ~𝒩 (5.6, 0.49)𝒩(μ, 0.09), withμ~𝒩 (6.0, 0.01)𝒩(μ, 0.04), withμ~𝒩 (6.5, 0.04)𝒩(μ, 0.01), withμ~𝒩 (0.30, 0.01)𝒩(μ, 0.16), withμ~𝒩 (0.33, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.35, 0.04)20#14, #16135𝒩(μ, 0.36), withμ~𝒩 (4.6, 0.25)𝒩(μ, 0.04), withμ~𝒩 (4.8, 0.16)𝒩(μ, 0.04), withμ~𝒩 (5.0, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.18, 0.04)𝒩(μ, 0.16), withμ~𝒩 (0.22, 0.09)𝒩(μ, 0.25), withμ~𝒩 (0.25, 0.49)21#15, #16140𝒩(μ, 1.00), withμ~𝒩 (5.8, 0.09)𝒩(μ, 1.00), withμ~𝒩 (6.2, 0.16)𝒩(μ, 0.81), withμ~𝒩 (6.5, 0.04)𝒩(μ, 0.01), withμ~𝒩 (0.30, 0.09)𝒩(μ, 0.16), withμ~𝒩 (0.32, 0.25)𝒩(μ, 0.04), withμ~𝒩 (0.35, 0.01)Figure 8
Illustrations of road network of Shuibuya Hydropower Project.


### 6.3. Computational Results

To verify the practicality and efficiency of the MCNFP optimization model under a birandom environment presented previously, the proposed MOBLPSO was implemented to determine the flow assignment amongst the transportation paths and amongst the carriers over a certain period using actual data from the Hubei Qingjiang Shuibuya Project Construction Company. After running the proposed MOBLPSO using MATLAB 7.0, the computational results were obtained and the efficiency of the proposed algorithm was proven.The computer running environment was an intel core 2 Duo 2.26 GHz clock pulse with 2048 MB memory. The problem was solved using the proposed algorithm with satisfactory solutions within 21 minutes on average, and the optimal solutions for lower level programming and the Pareto optimal solution set for upper level programming were worked out.The red dots in Figure9 show the Pareto optimal solutions, while the blue dots show the best position of particles in this iteration. The decision maker can choose a plan from these Pareto optimal solutions depending on their preference. For example, if the decision makers feel that the total direct cost objective is more important, they may sacrifice transportation time for more economical scheme. Thus they choose the absolute left of the Pareto optimal solution in the minimum cost network flow plan in Table 4. On the contrary, if decision makers feel that the transportation time objective is more important, they may choose minimum total transportation time cost and sacrifice more costs for the MCNFP, choosing the lowest of the Pareto optimal solutions. The minimum transportation time plan is in Table 5.Table 4
The most cost-effective plan for MCNFP in Shuibuya case.
Transportation pathj1234567891011Material volumexj950.000543.7631672.488256.317177.4321700.000766.6461411.95857.96663.4301279.675Material volumey1j426.432309.723576.20671.02093.874505.161371.885449.94336.22446.715520.477Material volumey2j413.345107.514587.05192.52144.360620.743267.585521.37842.17248.536364.225Material volumey3j369.912159.467509.23052.741127.025574.096285.792440.63734.15448.542394.974Transportation pathj12131415161718192021122Material volumexj370.325161.3571318.6431380.0001700.000346.384903.6161150.0001300.000747.349512.651Material volumey1j144.94743.373487.119398.457650.263158.874427.931381.462382.515363.82763.158Material volumey2j195.90142.686421.532476.339598.634111.618200.388430.262446.545264.17280.990Material volumey3j155.66246.672409.992505.204451.103171.539325.558338.276470.940373.284484.426Table 5
The most time-saving plan for MCNFP in Shuibuya case.
Transportation pathj1234567891011Material volumexj950.0001013.517758.507820.24757.7351700.000948.7621056.185253.39241.661647.444Material volumey1j265.783354.500308.790329.11545.111573.575492.092422.106157.70721.322363.277Material volumey2j427.952324.088400.510476.40339.031523.953354.168405.05448.49128.744222.921Material volumey3j293.735334.929390.005432.16538.571602.472398.627315.83156.69639.224199.717Transportation pathj1213141516171819202122Material volumexj1002.556200.4351279.5651380.0001700.0001057.049192.9511150.0001300.000500.306759.694Material volumey1j489.996100.614337.177417.547552.807274.888127.423337.068377.977274.467369.143Material volumey2j402.297109.078535.589454.315539.657493.531131.399451.603525.488255.822299.626Material volumey3j265.07431.314406.799508.137607.536439.66047.452361.329396.53565.138331.525Pareto optimal solutions of upper level programming for Shuibuya Hydropower Project.


### 6.4. Model Comparison to an MCNFP with Single Randomness

To highlight the advantages of our mathematical model (4.17), additional computational work was done using the proposed MOBLPSO to solve a similar MCNFP under a different uncertain environment, that is, the random environment. To carry out comparisons under a similar circumstance, analyses were conducted based on results from running the test problem 10 times. A detailed analysis follows.In order to guarantee a fair comparison between the MCNFP model with birandomness (denoted by MCNFP-birm) and a model just considering single randomness (denoted by MCNFP-rm), the random distribution for each related uncertain parameter in the MCNFP-rm was selected in the following way. Take the transportation coste¯̃11 (i.e., e¯̃11~𝒩(μ,0.64), with μ~𝒩(5.2,0.10)), for example. The stochastic nature of the expected value μ in the normal distribution 𝒩(μ,0.64) was ignored by using its expectation 5.2 as a representation while the variance 0.64 was retained. Thus, for the MCNFP-rm, the birandomness of the transportation cost e¯̃11 degenerated to a single randomness, in which the distribution of the transportation cost could be expressed as 𝒩(5.2,0.64). Since the variance of the random variable μ was sufficiently small (i.e., 0.1⩽1), the expectation of the random variable μ essentially reflected the most possible value over time. Thus, it was reasonable to select 𝒩(5.2,0.64) as the normal distribution for the transportation cost in the MCNFP-rm to compare with that in the MCNFP-birm. The transformation of the other related uncertain parameters followed a similar pattern. Thus the model for the MCNFP-rm was formulated and solved using the proposed MOBLPSO 10 times.As shown in Table6, the best, the worst, and the average results for the random type are higher than their counterparts for the MCNFP-birm. It is worth noting that the gaps between the best and the worst and between the best and the average solutions for MCNFP-rm are wider than the gaps of their counterparts for the birandom type. This shows that randomness creates a much larger solution space when uncertainty is introduced. Fortunately, the widened solution space with a further stochastic nature in the MCNFP-birm provides better solutions and they were successfully located by MOBLPSO, evidenced by the narrower gaps between the best and the worst and between the best and the average solutions found for the MCNFP-birm. This suggests that MOBLPSO is an effective and relatively efficient approach for solving the MCNFP under a birandom environment.Table 6
Comparisons of two types of model.
The direct cost objective (×103)The transportation time cost objective (×103)The transportation cost objective (×103)Birandom typeRandom typeBirandom typeRandom typeBirandom typeRandom typeBest result33.66233.83212.800312.799718.03618.067Average result34.54635.02212.802712.804518.23918.285Worst result35.42936.21112.805112.809218.44118.503


### 6.5. Algorithm Evaluation

Since the PSOPC for the lower level is nested in the MOPSO for the upper level, and the MOPSO is the main body of the proposed MOBLPSO, the evaluation of the MOPSO was mainly paid attention to. In the MOPSO, a multiobjective method is introduced to derive the Pareto optimal solution set for the upper level programming. This provides effective and nondominated alternate schemes for the construction contractor. Compared to the weight-sum method dealing with multiobjectives in [21], the solutions here are confirmed to be more practical.Comparing different optimization techniques experimentally always involves the notion of performance. In the case of multiobjective optimization, the definition of quality is substantially more complex than for single-objective optimization problems. There are many metrics of performance to measure the distance of the resulting nondominated set to the Paretooptimal front, the distribution of the solution found, and the extent of the obtained nondominated front [55].To gain further insight into the performance of the multiobjective method in the proposed algorithm, the procedure was run different times and the results are summarized in Figure9 and Table 7. As is shown in Figure 9, the amount and the distribution of Pareto optimal solutions in each iteration are satisfactory. For further expression of the efficiency of the convergence, three metrics of performance proposed by Zitzler et al. [55] were introduced: (1) the average distances of the resulting nondominated set to the Pareto-optimal front: the value of which decreases with an increase in the iterations, meaning that the program results come in toward the Pareto optimal front, which expresses the convergence of the algorithm; (2) the distribution in combination with the number of non-dominated solutions found: the higher the value, the better the distribution for an appropriate neighbourhood parameter; (3) the extent of the obtained nondominated fronts: it uses a maximum extent in each dimension to estimate the range to which the fronts spread out, which in this paper equals the distance of the two outer solutions.Table 7
Algorithm evaluation by metrics of performance for Pareto optimal sets.
IterationThe average distance metricThe distribution metricThe extent metric10.07810.147448.62271000.07380.346058.14232000.08780.403161.17223000.05260.696769.29704000.05540.808264.77735000.04350.893868.8648The metrics of performance for the Pareto optimal set is shown in Table7 to provide a satisfactory result for the efficiency of the convergence. Although there are some fluctuations in the three metrics in the 500 iterations, these do not affect the final result.To asses the efficiency and effectiveness of the MOBLPSO for the proposed MCNFP, the MOBLPSO results for the MCNFP in the Shuibuya Hydropower Project are compared with two other state-of-the-art heuristic algorithms, that is, a genetic algorithm for a multiobjective bilevel model (denoted by MOBLGA) [56] and a simulated annealing algorithm for a multiobjective bilevel model (denoted by MOBLSA) [57].In order to carry out the comparisons under a similar circumstance, the parameter selections for the MOBLGA and MOBLSA refer to those of the MOBLPSO, and nondominated alternate schemes are also employed for both. To measure the quality of the results obtained by the three algorithms, a weight sum method was introduced to detrmine one minimal weight sum for the objectives from the nondominated solutions. Thus the comparison could be implemented based on the unique measured criterion (i.e., the minimal weight sum of the objectives). To ensure the conformity validity of the multiobjectives, the division of the dimensions and a unifying of the order of magnitude need to be performed before the weight-sum procedure.Table8 shows the comparison results, that is, the minimal weight sum value of the two objectives, and the average computation times, obtained using the preceding approaches for the different combination of weights (i.e., ω1 and ω2 represent the weights of the two objectives, resp.). It is demonstrated that the MOBLPSO for the MCNFP can perform optimizing better than the MOBLGA, since MOBLGA may lead to a local search and need more computation time. On the other hand, the MOBLSA could get similar results to those from MOBLPSO, but computation was much slower than the MOBLPSO.Table 8
Computation time and memory used by MOBLPSO, MOBLGA and MOBLSA.
TypeCombination of weightsMinimal weight sum value of the two objectivesAverage computation time (min)ω1ω2MOBLPSOMOBLGAMOBLSAMOBLPSOMOBLGAMOBLSA1ω1=0.1ω2=0.91152.0271160.9891152.53120.85625.89934.4582ω1=0.2ω2=0.81024.2721032.4191024.27320.92626.03233.9213ω1=0.3ω2=0.7896.238903.848897.21421.05225.67934.0684ω1=0.4ω2=0.6769.472775.277769.55521.15825.32935.5685ω1=0.5ω2=0.5641.786646.706641.89720.96326.16935.0126ω1=0.6ω2=0.4514.148518.136514.21620.67826.96935.1427ω1=0.7ω2=0.3386.486388.529386.52121.20326.13534.7918ω1=0.8ω2=0.2258.792260.836258.79220.73425.47535.2919ω1=0.9ω2=0.1131.081132.143131.08121.13426.18934.815


## 6.1. Project Description

In this section, an earth-rock work transportation project in a large-scale water conservancy and hydropower construction project is taken as an example for our optimization method. The Shuibuya Hydropower Project was conducted in Badong County, which is located in the middle reaches of Qingjiang River in Sichuan province, China. The project is the first cascaded project on the Qingjiang main stream and the third important project following Geheyan and Gaobazhou in China. Once completed, it will provide a major power source to meet the peak load demand in the Central China Power Grid. The installed capacity and annual output of Shuibuya Power Plant are 1,600 MW and 3.92 GWh, respectively. The project has a powerful regulating ability with a normal pool level of 400 m and reservoir capacity of4.58×109m3. The project consists of a concrete-faced rock fill dam (CFRD), underground power house, a chute spillway on the left bank, and the sluice tunnel on the right bank. The dam is 233 m high and is the tallest of its kind in the world at present with a total volume of 15.64×106m3.


## 6.2. Data Collection

All detailed data for the Shuibuya Hydropower Project were obtained from the Hubei Qingjiang Shuibuya Project Construction Company. In a large-scale construction project, especially in a large water conservancy and hydropower construction project, earth-rock work is usually a primary material, and earth-rock work transportation occurs every day in excavation projects, borrow areas, filling projects, dumpling sites, and stockpile areas as it is turned over and needs to be replaced frequently (see Figure6). In the Shuibuya Hydropower Project, there are 4 excavation projects, 2 borrow areas, 3 stockpile areas, and 2 dumpling sites, The location and detailed information of borrow areas, dumpling sites, and stockpile areas, of the Shuibuya Hydropower Project is illustrated in Figure 7.Figure 6
Earth-rock work transportation relations in large-scale construction project.Figure 7
Layout of borrow areas, dumpling sites and stockpile areas of Shuibuya Hydropower Project.Three types of dump trucks (carriers) in the construction project are considered, which transport earth-rock work along different paths connecting the OD pairs, with the destination nodes having the practical demand of timeliness. All necessary data for every kind of carrier were calculated as shown in Table1, and Table 2 shows the details of the paths in the whole road network.Table 1
Information of carriers in Shuibuya Hydropower Project.
CarrierKind indexkType (heaped capacity; maximum payload)WeightwkDump truck1K30N-8.4 (16m3; 26 t)0.32Terex TA28 (17m3; 28 t)0.33Perlini DP366 (20m3; 36 t)0.4Table 2
Related data about information of paths in earth-rock work transportation network.
OD pairsDemandQod (m3)Transportation pathjDirect cost (CNY/m3)Constraint time (h)Excavation project 1,2→dumpling site 19501⇔#5→#7→#91.801.45Excavation project 1,2→filling project 1,226502⇔ #5→#3→#1→#62.502.153⇔#5→#1→#61.701.354⇔ #5→#7→#8→#62.402.205⇔ #5→#7→#61.901.40Excavation project 1,2→stockpile area 117006⇔ #5→#71.201.25Excavation project 3,4→filling project 3,423007⇔#11→#10→#14→#162.302.258⇔ #11→#13→#161.901.459⇔ #11→#10→#13→#162.702.3010⇔ #11→#10→#14→#15→#163.003.50Excavation project 3,4→dumpling site 2165011⇔ #11→#13→#121.601.5012⇔ #11→#10→#13→#122.332.50Excavation project 3,4→stockpile area 2148013⇔ #11→#10→#132.601.7514⇔ #11→#131.101.25Excavation project 3,4→stockpile area 3138015⇔ #11→#10→#14→#152.102.33Borrow area 1→filling project 1,2170016⇔ #2→#4→#61.501.50Borrow area 2→filling project 3,4125017⇔ #14→#15→#161.801.5518⇔ #14→#161.201.20Stockpile area 1→filling project 1,2115019⇔ #7→#8→#61.901.45Stockpile area 2→filling project 3,4130020⇔ #13→#161.401.17Stockpile area 3→filling project 1,2126021⇔ #15→#7→#8→#62.402.1722⇔ #15→#7→#62.201.60The transportation network in the project includes an internal road network and an external road network. In this case, only the internal road network is considered. The internal road network is composed of 20 trunk roads located on the left and right banks, with 11 on the left and 9 on the right, with a cross-river bridge connecting the left and right banks. Therefore, 21 links are considered in this paper. In order to apply the proposed methods more conveniently, adjacent roads of the same type have been combined and road shapes have been ignored. An abstracted transportation network is illustrated in Figure8. For each link in the transportation network, there are a free flow birandom travel time t¯̃i0k and birandom transportation cost e¯̃ik. The corresponding data of which are stated in Table 3. In order to collect the data of transportation time and costs, investigations and surveys were made to obtain historical data from the financial department and experienced engineers of construction team in Hubei Qingjiang Shuibuya Project Construction Company. Since the transportation time and costs for each arc of the path change over time, the data are classified into different parts based on different periods. Both transportation time and costs are assumed to approximately follow normal distributions for each period, and the two parameters (expected value and variance) for the normal distributions are estimated by using maximum likelihood estimation, which justified by a chi-square goodness-of-fit test. By comparing the normal distributions for the same transportation time and costs in different periods, it can be found that the expected values of the aforementioned normal distributions also approximately follow a similar random distribution pattern, which has also been justified using a chi-square goodness-of-fit test. It should be noted that since the variance fluctuations are quite insignificant, the median values of the variances in different periods are selected as the variance for the previous normal distribution. The predicted confidence levels given by the decision maker are, respectively, θ=0.9 and δ=0.85.Table 3
Birandom free flow travel timet¯̃i0k and birandom transportation cost e¯̃ik.
ArciCorresponding nodesArc capacityri (n/h)Transportation costeik (CNY/unit)Free transportation timeti0k (h)Carrier 1Carrier 2Carrier 3Carrier 1Carrier 2Carrier 31#1, #3105𝒩(μ, 0.64), withμ~𝒩 (5.2, 0.10)𝒩(μ, 1.00), withμ~𝒩 (6.0, 0.09)𝒩(μ, 0.01), withμ~𝒩 (6.8, 1.00)𝒩(μ, 0.25), withμ~𝒩 (0.2, 0.01)𝒩(μ, 0.81),withμ~𝒩 (0.24, 0.09)𝒩(μ, 0.16), withμ~𝒩(0.28, 0.04)2#1, #5110𝒩(μ, 1.00), withμ~𝒩(3.2, 0.64)𝒩(μ, 0.49), withμ~𝒩 (3.5, 0.16)𝒩(μ, 0.81), withμ~𝒩 (3.8, 0.25)𝒩(μ, 0.36), withμ~𝒩 (0.15, 0.04)𝒩(μ, 0.25), withμ~𝒩 (0.16, 0.04)𝒩(μ, 0.09), withμ~𝒩 (0.18, 0.01)3#1, #6135𝒩(μ, 4.00), withμ~𝒩 (7.2, 1.00)𝒩(μ, 0.81), withμ~𝒩 (7.5, 0.64)𝒩(μ, 0.64), withμ~𝒩(7.8, 0.49)𝒩(μ, 0.25), withμ~𝒩 (0.30, 0.16)𝒩(μ, 0.36), withμ~𝒩 (0.35, 0.09)𝒩(μ, 0.16), withμ~𝒩 (0.38, 0.01)4#2, #4115𝒩(μ, 1.00), withμ~𝒩 (3.7, 0.49)𝒩(μ, 4.00), withμ~𝒩 (4.0, 0.81)𝒩(μ, 0.01), withμ~𝒩 (4.3, 0.01)𝒩(μ, 0.81), withμ~𝒩 (0.14, 0.16)𝒩(μ, 0.36), withμ~𝒩 (0.18, 0.09)𝒩(μ, 0.64), withμ~𝒩 (0.19, 0.01)5#3, #5108𝒩(μ, 0.36), withμ~𝒩 (2.2, 0.01)𝒩(μ, 0.81), withμ~𝒩 (2.4, 0.09)𝒩(μ, 0.64), withμ~𝒩 (2.6, 0.16)𝒩(μ, 0.01), withμ~𝒩 (0.14, 0.01)𝒩(μ, 0.25), withμ~𝒩 (0.18, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.22, 0.01)6#4, #6112𝒩(μ, 0.81), withμ~𝒩 (3.2, 0.16)𝒩(μ, 1.00), withμ~𝒩 (4.0, 0.09)𝒩(μ, 0.49), withμ~𝒩 (4.8, 0.09)𝒩(μ, 0.64), withμ~𝒩 (0.12, 0.16)𝒩(μ, 0.09), withμ~𝒩 (0.15, 0.01)𝒩(μ, 0.25), withμ~𝒩 (0.18, 0.01)7#5, #7180𝒩(μ, 1.00), withμ~𝒩(5.0, 0.01)𝒩(μ, 0.64), withμ~𝒩 (5.2, 0.25)𝒩(μ, 0.36), withμ~𝒩 (5.4, 0.16)𝒩(μ, 0.01), withμ~𝒩 (0.25, 0.01)𝒩(μ, 0.04), withμ~𝒩 (0.28, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.32, 0.01)8#6, #7145𝒩(μ, 0.49), withμ~𝒩 (4.5, 0.01)𝒩(μ, 0.36), withμ~𝒩 (4.2, 0.16)𝒩(μ, 0.49), withμ~𝒩 (3.4, 0.16)𝒩(μ, 0.09), withμ~𝒩 (0.16, 0.01)𝒩(μ, 0.04), withμ~𝒩 (0.19, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.23, 0.01)9#6, #8160𝒩(μ, 0.16), withμ~𝒩 (3.6, 0.25)𝒩(μ, 0.81), withμ~𝒩 (4.0, 1.00)𝒩(μ, 1.00), withμ~𝒩 (4.4, 0.09)𝒩μ, 0.04), withμ~𝒩 (0.13, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.15, 0.04)𝒩(μ, 0.09), withμ~𝒩 (0.17, 0.01)10#7, #8165𝒩(μ, 1.00), withμ~𝒩 (4.0, 0.81)𝒩(μ, 1.00), withμ~𝒩 (4.4, 0.49)𝒩(μ, 0.64), withμ~𝒩 (4.8, 0.25)𝒩(μ, 0.01), withμ~𝒩 (0.16, 0.01)𝒩(μ, 0.04), withμ~𝒩 (0.19, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.23, 0.04)11#7, #9100𝒩(μ, 4.00), withμ~𝒩 (7.2, 0.09)𝒩(μ, 0.16), withμ~𝒩 (7.6, 0.81)𝒩(μ, 0.25), withμ~𝒩 (8.0, 0.16)𝒩(μ, 0.01), withμ~𝒩 (0.36, 0.01)𝒩(μ, 0.04), withμ~𝒩 (0.39,0.01)𝒩(μ, 0.09), withμ~𝒩 (0.42, 0.01)12#7, #15130𝒩(μ, 0.81), withμ~𝒩 (8.0, 0.81)𝒩(μ, 0.64), withμ~𝒩 (8.5, 0.50)𝒩(μ, 1.00), withμ~𝒩 (8.7, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.33, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.36, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.38, 0.01)13#10, #11225𝒩(μ, 4.00), withμ~𝒩 (5.0, 0.81)𝒩(μ, 0.64), withμ~𝒩 (6.0, 0.16)𝒩(μ, 0.25), withμ~𝒩 (6.5, 0.09)𝒩(μ, 0.01), withμ~𝒩 (0.30, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.34, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.36, 0.01)14#10, #13155𝒩(μ, 0.81), withμ~𝒩 (5.5, 0.36)𝒩(μ, 1.00), withμ~𝒩 (5.6, 0.49)𝒩(μ, 0.04), withμ~𝒩 (5.8, 0.16)𝒩(μ, 0.01), withμ~𝒩 (0.22, 0.09)𝒩(μ, 0.04), withμ~𝒩 (0.26, 0.16)𝒩(μ, 0.04), withμ~𝒩 (0.30, 0.09)15#10, #14150𝒩(μ, 0.64), withμ~𝒩 (4.5, 0.49)𝒩(μ, 1.00), withμ~𝒩 (4.6, 0.25)𝒩(μ, 0.09), withμ~𝒩 (4.8, 0.09)𝒩(μ, 0.04), withμ~𝒩 (0.22, 0.01)𝒩(μ, 0.01), withμ~𝒩 (0.25, 0.04)𝒩(μ, 0.01), withμ~𝒩 (0.27, 0.01)16#11, #13170𝒩(μ, 0.64), withμ~𝒩 (3.6, 0.16)𝒩(μ, 0.09), withμ~𝒩 (4.0, 0.01)𝒩(μ, 0.16), withμ~𝒩 (4.4, 0.25)𝒩(μ, 0.04), withμ~𝒩 (0.16, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.18, 0.04)𝒩(μ, 0.01), withμ~𝒩 (0.20, 0.01)17#12, #13125𝒩(μ, 0.25), withμ~𝒩 (4.2, 0.16)𝒩(μ, 0.49), withμ~𝒩 (4.3, 0.09)𝒩(μ, 1.00), withμ~𝒩(4.4, 0.16)𝒩(μ, 0.01), withμ~𝒩 (0.22, 0.01)𝒩(μ, 0.04), withμ~𝒩 (0.24, 0.09)𝒩(μ, 0.09), withμ~𝒩 (0.26, 0.01)18#13, #16145𝒩(μ, 1.00), withμ~𝒩 (4.5, 0.16)𝒩(μ, 0.04), withμ~𝒩 (4.7, 0.81)𝒩(μ, 0.49), withμ~𝒩 (4.9, 1.00)𝒩(μ, 0.01), withμ~𝒩 (0.21, 0.64)𝒩(μ, 1.00), withμ~𝒩 (0.23, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.25, 0.16)19#14, #15140𝒩(μ, 0.25), withμ~𝒩 (5.6, 0.49)𝒩(μ, 0.09), withμ~𝒩 (6.0, 0.01)𝒩(μ, 0.04), withμ~𝒩 (6.5, 0.04)𝒩(μ, 0.01), withμ~𝒩 (0.30, 0.01)𝒩(μ, 0.16), withμ~𝒩 (0.33, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.35, 0.04)20#14, #16135𝒩(μ, 0.36), withμ~𝒩 (4.6, 0.25)𝒩(μ, 0.04), withμ~𝒩 (4.8, 0.16)𝒩(μ, 0.04), withμ~𝒩 (5.0, 0.01)𝒩(μ, 0.09), withμ~𝒩 (0.18, 0.04)𝒩(μ, 0.16), withμ~𝒩 (0.22, 0.09)𝒩(μ, 0.25), withμ~𝒩 (0.25, 0.49)21#15, #16140𝒩(μ, 1.00), withμ~𝒩 (5.8, 0.09)𝒩(μ, 1.00), withμ~𝒩 (6.2, 0.16)𝒩(μ, 0.81), withμ~𝒩 (6.5, 0.04)𝒩(μ, 0.01), withμ~𝒩 (0.30, 0.09)𝒩(μ, 0.16), withμ~𝒩 (0.32, 0.25)𝒩(μ, 0.04), withμ~𝒩 (0.35, 0.01)Figure 8
Illustrations of road network of Shuibuya Hydropower Project.


## 6.3. Computational Results

To verify the practicality and efficiency of the MCNFP optimization model under a birandom environment presented previously, the proposed MOBLPSO was implemented to determine the flow assignment amongst the transportation paths and amongst the carriers over a certain period using actual data from the Hubei Qingjiang Shuibuya Project Construction Company. After running the proposed MOBLPSO using MATLAB 7.0, the computational results were obtained and the efficiency of the proposed algorithm was proven.The computer running environment was an intel core 2 Duo 2.26 GHz clock pulse with 2048 MB memory. The problem was solved using the proposed algorithm with satisfactory solutions within 21 minutes on average, and the optimal solutions for lower level programming and the Pareto optimal solution set for upper level programming were worked out.The red dots in Figure9 show the Pareto optimal solutions, while the blue dots show the best position of particles in this iteration. The decision maker can choose a plan from these Pareto optimal solutions depending on their preference. For example, if the decision makers feel that the total direct cost objective is more important, they may sacrifice transportation time for more economical scheme. Thus they choose the absolute left of the Pareto optimal solution in the minimum cost network flow plan in Table 4. On the contrary, if decision makers feel that the transportation time objective is more important, they may choose minimum total transportation time cost and sacrifice more costs for the MCNFP, choosing the lowest of the Pareto optimal solutions. The minimum transportation time plan is in Table 5.Table 4
The most cost-effective plan for MCNFP in Shuibuya case.
Transportation pathj1234567891011Material volumexj950.000543.7631672.488256.317177.4321700.000766.6461411.95857.96663.4301279.675Material volumey1j426.432309.723576.20671.02093.874505.161371.885449.94336.22446.715520.477Material volumey2j413.345107.514587.05192.52144.360620.743267.585521.37842.17248.536364.225Material volumey3j369.912159.467509.23052.741127.025574.096285.792440.63734.15448.542394.974Transportation pathj12131415161718192021122Material volumexj370.325161.3571318.6431380.0001700.000346.384903.6161150.0001300.000747.349512.651Material volumey1j144.94743.373487.119398.457650.263158.874427.931381.462382.515363.82763.158Material volumey2j195.90142.686421.532476.339598.634111.618200.388430.262446.545264.17280.990Material volumey3j155.66246.672409.992505.204451.103171.539325.558338.276470.940373.284484.426Table 5
The most time-saving plan for MCNFP in Shuibuya case.
Transportation pathj1234567891011Material volumexj950.0001013.517758.507820.24757.7351700.000948.7621056.185253.39241.661647.444Material volumey1j265.783354.500308.790329.11545.111573.575492.092422.106157.70721.322363.277Material volumey2j427.952324.088400.510476.40339.031523.953354.168405.05448.49128.744222.921Material volumey3j293.735334.929390.005432.16538.571602.472398.627315.83156.69639.224199.717Transportation pathj1213141516171819202122Material volumexj1002.556200.4351279.5651380.0001700.0001057.049192.9511150.0001300.000500.306759.694Material volumey1j489.996100.614337.177417.547552.807274.888127.423337.068377.977274.467369.143Material volumey2j402.297109.078535.589454.315539.657493.531131.399451.603525.488255.822299.626Material volumey3j265.07431.314406.799508.137607.536439.66047.452361.329396.53565.138331.525Pareto optimal solutions of upper level programming for Shuibuya Hydropower Project.


## 6.4. Model Comparison to an MCNFP with Single Randomness

To highlight the advantages of our mathematical model (4.17), additional computational work was done using the proposed MOBLPSO to solve a similar MCNFP under a different uncertain environment, that is, the random environment. To carry out comparisons under a similar circumstance, analyses were conducted based on results from running the test problem 10 times. A detailed analysis follows.In order to guarantee a fair comparison between the MCNFP model with birandomness (denoted by MCNFP-birm) and a model just considering single randomness (denoted by MCNFP-rm), the random distribution for each related uncertain parameter in the MCNFP-rm was selected in the following way. Take the transportation coste¯̃11 (i.e., e¯̃11~𝒩(μ,0.64), with μ~𝒩(5.2,0.10)), for example. The stochastic nature of the expected value μ in the normal distribution 𝒩(μ,0.64) was ignored by using its expectation 5.2 as a representation while the variance 0.64 was retained. Thus, for the MCNFP-rm, the birandomness of the transportation cost e¯̃11 degenerated to a single randomness, in which the distribution of the transportation cost could be expressed as 𝒩(5.2,0.64). Since the variance of the random variable μ was sufficiently small (i.e., 0.1⩽1), the expectation of the random variable μ essentially reflected the most possible value over time. Thus, it was reasonable to select 𝒩(5.2,0.64) as the normal distribution for the transportation cost in the MCNFP-rm to compare with that in the MCNFP-birm. The transformation of the other related uncertain parameters followed a similar pattern. Thus the model for the MCNFP-rm was formulated and solved using the proposed MOBLPSO 10 times.As shown in Table6, the best, the worst, and the average results for the random type are higher than their counterparts for the MCNFP-birm. It is worth noting that the gaps between the best and the worst and between the best and the average solutions for MCNFP-rm are wider than the gaps of their counterparts for the birandom type. This shows that randomness creates a much larger solution space when uncertainty is introduced. Fortunately, the widened solution space with a further stochastic nature in the MCNFP-birm provides better solutions and they were successfully located by MOBLPSO, evidenced by the narrower gaps between the best and the worst and between the best and the average solutions found for the MCNFP-birm. This suggests that MOBLPSO is an effective and relatively efficient approach for solving the MCNFP under a birandom environment.Table 6
Comparisons of two types of model.
The direct cost objective (×103)The transportation time cost objective (×103)The transportation cost objective (×103)Birandom typeRandom typeBirandom typeRandom typeBirandom typeRandom typeBest result33.66233.83212.800312.799718.03618.067Average result34.54635.02212.802712.804518.23918.285Worst result35.42936.21112.805112.809218.44118.503


## 6.5. Algorithm Evaluation

Since the PSOPC for the lower level is nested in the MOPSO for the upper level, and the MOPSO is the main body of the proposed MOBLPSO, the evaluation of the MOPSO was mainly paid attention to. In the MOPSO, a multiobjective method is introduced to derive the Pareto optimal solution set for the upper level programming. This provides effective and nondominated alternate schemes for the construction contractor. Compared to the weight-sum method dealing with multiobjectives in [21], the solutions here are confirmed to be more practical.Comparing different optimization techniques experimentally always involves the notion of performance. In the case of multiobjective optimization, the definition of quality is substantially more complex than for single-objective optimization problems. There are many metrics of performance to measure the distance of the resulting nondominated set to the Paretooptimal front, the distribution of the solution found, and the extent of the obtained nondominated front [55].To gain further insight into the performance of the multiobjective method in the proposed algorithm, the procedure was run different times and the results are summarized in Figure9 and Table 7. As is shown in Figure 9, the amount and the distribution of Pareto optimal solutions in each iteration are satisfactory. For further expression of the efficiency of the convergence, three metrics of performance proposed by Zitzler et al. [55] were introduced: (1) the average distances of the resulting nondominated set to the Pareto-optimal front: the value of which decreases with an increase in the iterations, meaning that the program results come in toward the Pareto optimal front, which expresses the convergence of the algorithm; (2) the distribution in combination with the number of non-dominated solutions found: the higher the value, the better the distribution for an appropriate neighbourhood parameter; (3) the extent of the obtained nondominated fronts: it uses a maximum extent in each dimension to estimate the range to which the fronts spread out, which in this paper equals the distance of the two outer solutions.Table 7
Algorithm evaluation by metrics of performance for Pareto optimal sets.
IterationThe average distance metricThe distribution metricThe extent metric10.07810.147448.62271000.07380.346058.14232000.08780.403161.17223000.05260.696769.29704000.05540.808264.77735000.04350.893868.8648The metrics of performance for the Pareto optimal set is shown in Table7 to provide a satisfactory result for the efficiency of the convergence. Although there are some fluctuations in the three metrics in the 500 iterations, these do not affect the final result.To asses the efficiency and effectiveness of the MOBLPSO for the proposed MCNFP, the MOBLPSO results for the MCNFP in the Shuibuya Hydropower Project are compared with two other state-of-the-art heuristic algorithms, that is, a genetic algorithm for a multiobjective bilevel model (denoted by MOBLGA) [56] and a simulated annealing algorithm for a multiobjective bilevel model (denoted by MOBLSA) [57].In order to carry out the comparisons under a similar circumstance, the parameter selections for the MOBLGA and MOBLSA refer to those of the MOBLPSO, and nondominated alternate schemes are also employed for both. To measure the quality of the results obtained by the three algorithms, a weight sum method was introduced to detrmine one minimal weight sum for the objectives from the nondominated solutions. Thus the comparison could be implemented based on the unique measured criterion (i.e., the minimal weight sum of the objectives). To ensure the conformity validity of the multiobjectives, the division of the dimensions and a unifying of the order of magnitude need to be performed before the weight-sum procedure.Table8 shows the comparison results, that is, the minimal weight sum value of the two objectives, and the average computation times, obtained using the preceding approaches for the different combination of weights (i.e., ω1 and ω2 represent the weights of the two objectives, resp.). It is demonstrated that the MOBLPSO for the MCNFP can perform optimizing better than the MOBLGA, since MOBLGA may lead to a local search and need more computation time. On the other hand, the MOBLSA could get similar results to those from MOBLPSO, but computation was much slower than the MOBLPSO.Table 8
Computation time and memory used by MOBLPSO, MOBLGA and MOBLSA.
TypeCombination of weightsMinimal weight sum value of the two objectivesAverage computation time (min)ω1ω2MOBLPSOMOBLGAMOBLSAMOBLPSOMOBLGAMOBLSA1ω1=0.1ω2=0.91152.0271160.9891152.53120.85625.89934.4582ω1=0.2ω2=0.81024.2721032.4191024.27320.92626.03233.9213ω1=0.3ω2=0.7896.238903.848897.21421.05225.67934.0684ω1=0.4ω2=0.6769.472775.277769.55521.15825.32935.5685ω1=0.5ω2=0.5641.786646.706641.89720.96326.16935.0126ω1=0.6ω2=0.4514.148518.136514.21620.67826.96935.1427ω1=0.7ω2=0.3386.486388.529386.52121.20326.13534.7918ω1=0.8ω2=0.2258.792260.836258.79220.73425.47535.2919ω1=0.9ω2=0.1131.081132.143131.08121.13426.18934.815


## 7. Conclusions and Future Research

In this paper, a multiobjective bilevel programming model under birandom environment for an MCNFP in a large-scale construction project was formulated. The contributions of this paper to the literature are as follows. Firstly, the multiobjective bilevel model for the minimum cost network flow problem in a large-scale construction project focused on here was found to provide a more reasonable expression of the proposed problem, where the upper level aims at optimizing the material flow assignment along the transportation paths and the lower level decides on the flow of each carrier transports on the paths. Secondly, because of the complicated realistic decision systems, this study employs birandom variables to characterize the hybrid uncertain environment. The application of birandom variables makes the proposed programming model more suitable for describing a vague and uncertain situation in the real world. Further, the birandom uncertainty model was converted into an expectation multiobjective bilevel programming model with chance constraints. Thirdly, in order to solve the NP-hard multiobjective bilevel problem, a very effective and relatively efficient algorithm (i.e., MOBLPSO) was developed by employing both a MOPSO and a PSOPC. Finally, the Shuibuya Hydropower Project was used here as a practical application example. The MOBLPSO results for the preceding project example were compared with MOBLGA and MOBLSA methods, which demonstrated the validity of the proposed mathematical model and the effectiveness of the proposed MOBLPSO method in handling complex problems.Further research is necessary to identify further properties to develop a more effective method for solving other practical problems: (1) the formulation of an MCNFP for manifold materials rather than only one type of material transportation network in large-scale construction projects, (2) the investigation of other new approaches such as an automated design methodology and dependent chance programming to handle the birandom variables more reasonably and effectively, (3) the development of more efficient solution methods to solve multiobjective bilevel programming problems. Each of these areas is very important and equally worthy of attention. It should be mentioned that there are several commercial solvers that can efficiently solve large-scale nonlinear problems such as MINOS, CONOPT and SNOPT. However, when solving bilevel programming with nonlinear and non-differentiable piecewise objective functions and constraints like the MCNFP discussed in this paper, these solvers may face difficulties to deal with the nondifferentiability and nonconvexity by employing the exact techniques such as enumeration method, Karush-Kuhn-Tucker method, and penalty function approach. The future research may seek to address this issue with alternative exact techniques.



---

*Source: 463976-2012-06-18.xml*