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.
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 Section 2, 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 Figure 1.
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 [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.
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
:index of origin node;:index of destination node;:set of carriers, is an index;:set of arcs in the transportation network, is an index;:set of paths in the transportation network, is an index;:set of origin-destination (OD) pairs, ;:set of arcs in transportation path ; :set of paths from origin node to destination node .
Certain Parameters
:maximal passing capacity of arc ;:weight for carrier in the transportation network;:volume capacity of carrier ;:transportation time constraint represents the time of transportation path in units in which the material demand between all OD pairs has to be transported;:a binary variable equal to 1 if and only if arc is a segment of transportation path for carrier ;:transportation demand of the material from origin node to destination node ;:direct cost of unit volume material using transportation ;ceiling operator rounding upward to integer.
Uncertain Parameters
:unit transportation cost of material flow on arc for carrier ;:free transportation time of material flow on arc for carrier ;:transportation time of material flow on arc for carrier .
Decision Variables
:volume of material flow on transportation path , which is the decision variable of the upper level;:volume of material flow transported by carrier through path , 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 to a collection of random variables such that for any Borel subset of the real line the induced function is a measurable function with respect to .
For each given Borel subset of the real line , the function is a random variable defined on the probability space .
From Definition 3.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 , and . Assume that is a function on as follows: where is a uniformly distributed random variable on , are normally distributed random variables with a mean 1 and a standard variance 0.5, and is a standard normally distributed random variable with a mean 0 and a standard variance 1, that is, , , and . From Definition 3.1, is clearly a birandom variable as shown in Figure 3.
Example 3.3. Assume that and are two random variables defined on , and for any , holds; then random variable is a birandom variable.
Example 3.4. Let be a random variable defined on the probability space satisfying , where is also a normally distributed random variable on with the mean and variance . 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: where the number of random variable of and is not less than one. The normal birandom variable is denoted by .
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 path 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 carrier through path can be denoted by . Since 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 ; then the objective function of the lower level model can be formulated as
Generally, a path can be represented by a sequence of adjacent arcs. A binary variable is introduced to determine whether an arc is a segment of path for carrier :
(2) Constraints of the Lower Level
For transportation time, each carrier requires the transport of material from the source to the destination on schedule . If not, a delay cost is applied. represents the total travel time of carrier on transportation path , in which is usually represented by a non-decreasing function (i.e., Bureau of Public Roads (BPR) function) [33] as follows:
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 does not exceed because of the birandom variable . In practical problem, the decision makers often provide an appropriate budget in advance, to ensure that the restriction is, to a certain extent, satisfied, that is to maximize the probability of the random event under a given confidence level, which can be written as follows:
Here the decision makers' aspiration level is indicated as , so we use a “” to ensure that the constraint holds at the predetermined confidence level. Additionally, based on probability theory, a further “” 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 arc cannot exceed the maximal capacity of the arc , which produces the following constraint:
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:
The path flow in path used by carrier should not be negative, such that
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 flow through transportation path . Thus, the decision variable on the upper level is .
(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
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 carrier in each path can be described as , and, since is a birandom variable, 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 , so the second objective function of the upper level can be described as follows:
(2) Constraints of the Upper Level
According to the basic assumptions in Section 3.1, it is stipulated that the demands between all OD pairs should be satisfied. Thus,
The following constraint ensures that the sum of the weights is equal to 1:
In order to describe the nonnegative variables, the constraints in (3.14) are presented:
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:
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 and 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 . Then the expected value of is defined by
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 ; then the expected value of birandom variable can be defined as follows: 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 . If the expected value of the random variable is finite for each , then is a random variable on .
Remark 4.4. It should be noted that the expected value operator , which appears on both sides of the previous definition of , is overloaded. In fact, symbol represents different meanings. That is to say, overloading allows us to use the same symbol for different expected value operators, because we can deduce the meaning from the type of argument.
Lemma 4.5. Let be a normal birandom variable, subject to a normal distribution where . Then the objective function in (3.15) is transformed into the following equivalent objective function:
Proof. From the assumption, for any , is an independent random variable, where . From Definition 4.2, since , and by Definition 4.1, function (4.5) is transformed as follows: and then This completes the proof.
Let be a normal birandom variable subject to a normal distribution where . Then similarly the objective function of the lower level is transformed into crisp equivalents:
Lemma 4.6. Assume that is a normal birandom variable, where is a normal distributed random variable characterized by ; then is also a random variable, where
Then the following constraint for the first constraint is derived in (3.15):
being equivalent to the following equation:
Proof. From the assumption, it is known that for any , is an independent random variable, so it follows that
is also a normally distributed random variable, where
Then
Since , then , where
for given confidence levels , so,
where .
The proof is completed.
Based on Lemmas 4.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: where , 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 , then the lower decision maker reacts 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: where is the velocity of th particle at the th dimension in the th iteration, is the inertia weight, is the position of th particle at the th dimension, and are random numbers in the range , and are personal and global best position acceleration constants, respectively, while, is personal best position of the th particle at the th dimension, and is the global best position at the th 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:: iteration index of the upper and lower level, and ;: particle index of the upper and lower level, and ;: index of transportation path, ;: index of carrier, ;: uniform distributed random number within [0,1];: inertia weight in the th iteration of the upper level;: inertia weight in the th iteration of the lower level;: particle selected randomly from the swarm at the th dimension in the th iteration;, : personal and global best position acceleration constant of the upper level;, : personal and global best position acceleration constant of the lower level;: passive congregation coefficient of the lower level;,: maximum and minimum position value of the upper level;: vector position of the th particle in the th iteration, ;: vector position of the th particle in the th iteration, ;: vector velocity of the th particle in the th iteration, ;: vector velocity of the th particle in the th iteration, ;: vector personal best position of the th particle in the th iteration, ;: vector personal best position of the th particle in the th iteration, ;: vector global best position in the th iteration, ;: the positions of the particles that represent nondominated vectors in the repository;: fitness value of ; and: fitness value of .
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 (Algorithms 1 and 2), in which is initially set equal to the initial position of particle .
| ||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||
The repository that stores the positions of the particles that represent nondominated vectors is denoted by ARC; then the velocity of each th particle of the upper level is updated using the following equations: where is a solution randomly selected from the repository in iteration , which can improve significantly the ability of global convergence by avoiding being trapped in a stagnant state in finite iterations [49, 50]. The index 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 th particle of the lower level:
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, , the output of PSOPC for lower level programming is the input of the MOPSO for upper level programming; and , 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 with for the upper level, , , and (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., ) and lower level (i.e., ) are set to be varying with iteration as follows:
where the iteration numbers are (i.e., for the upper level) and (i.e., for the lower level), and and (for ). 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 . For , generate the position of the particle with an integer random position (note that every particle in the upper level consists of dimensions in this study). In the lower level, set iteration . For , generate the position of the particle with an integer random position (note that every particle in the lower level consists of dimensions in this study).
Decode Particles into Solutions
Decode particles into solutions: for , decode to a solution as . For , decode to a solution as . Mapping between one potential solution for the upper and lower level of the MCNFP and particle representation is shown in Figure 4.
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; particle 's fitness value is a matrix, namely,
The fitness value of the particle in the lower level is
Figure 5 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.
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 of . 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 .
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 Figure 6). 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.
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 Table 1, and Table 2 shows the details of the paths in the whole road network.
The 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 Figure 8. For each link in the transportation network, there are a free flow birandom travel time and birandom transportation cost . 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, and .
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 Figure 9 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.
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 cost (i.e., , with ), for example. The stochastic nature of the expected value in the normal distribution 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 degenerated to a single randomness, in which the distribution of the transportation cost could be expressed as . Since the variance of the random variable was sufficiently small (i.e., ), the expectation of the random variable essentially reflected the most possible value over time. Thus, it was reasonable to select 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 Table 6, 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.
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 Figure 9 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.
The metrics of performance for the Pareto optimal set is shown in Table 7 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.
Table 8 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., and 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.
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.
Acknowledgments
This paper was supported by the Key Program of NSFC (Grant no. 70831005) and “985” Program of Sichuan University “Innovative Research Base for Economic Development and Management.” The authors would like to give their great appreciates to the editors and anonymous referees for their helpful and constructive comments and suggestions, which have helped to improve this paper.