# Hazmats Transportation Network Design Model with Emergency Response under Complex Fuzzy Environment

Authors: Jiuping Xu; Jun Gang; Xiao Lei
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

A bilevel optimization model for a hazardous materials transportation network design is presented which considers an emergency response teams location problem. On the upper level, the authority designs the transportation network to minimize total transportation risk. On the lower level, the carriers first choose their routes so that the total transportation cost is minimized. Then, the emergency response department locates their emergency service units so as to maximize the total weighted arc length covered. In contrast to prior studies, the uncertainty associated with transportation risk has been explicitly considered in the objective function of our mathematical model. Specifically, our research uses a complex fuzzy variable to model transportation risk. An improved artificial bee colony algorithm with priority-based encoding is also applied to search for the optimal solution to the bilevel model. Finally, the efficiency of the proposed model and algorithm is evaluated using a practical case and various computing attributes.


---

## Body

## 1. Introduction

The transportation of hazardous materials (hazmats) has always been an important issue because of the risk associated with an accidental release of hazardous materials during transportation [1]. Globally, more than 4 billion tons of hazmat are transported annually [2] with more than 2 billion tons of hazmat cargo, 6773 companies, and over 116.4 thousand vehicles in China alone [3]. Due to the potential magnitude of accidents, and the risks associated with incidents involving hazmat shipments, the public is very sensitive to the dangers. As a consequence, these risks have attracted considerable attention from governments, the public, and scholars [4–6]. Today, it has been realized that hazmat transportation risks can be significantly reduced through the design of better transportation networks.Hazmat transportation network design has been increasingly studied since 2004. Kara and Verter [7] first considered a hazmat transportation network design problem using a bilevel integer programming model. Erkut and Gzara [8] generalized Kara and Verter’s model to an undirected network case. Bianco et al. [9] proposed a bilevel network flow model which aimed to minimize total risk and to guarantee risk equity. Most of the hazmat network design models in previous researches considered the governments’ and carriers’ points of view in trying to mitigate the total risk. However, they ignored the effect of the emergency response services. When a hazmat incident occurs, effective emergency response is crucial for mitigating the undesirable incident consequences [10]. That is, hazmats transportation risk may be drastically reduced by designing a transportation network which secures a timely and unobstructed provision of emergency response services [11]. Therefore, the goal of network design is not only to encourage carriers to choose routes with a lower risk, but also to ensure that the chosen road is close to emergency response service sites.For hazmat transportation network design, risk assessment is the basis of the decision making. There are two types of methods for risk assessment. One is called quantitative risk analysis method [12] such as frequency analysis [13], logical diagram-based techniques [14], and modelling of impact area [1, 15]. They can work well when data is sufficient to model the risk. However, in many cases, there are no enough available historical data. As a result, fuzzy approaches [16–18] have been used to model the risk. In these fuzzy approaches, fuzzy theory was often used to model only a single parameter such as accident probability or accident consequence which tended to lead to a lack of a comprehensive and uniform description for risk uncertainty. However, transportation risk is a complex uncertain parameter, which has much imperfect information. For example, both of the accident consequence and the accident probability may be fuzzy; as a result, the risk sometimes could be stated as “it is about in the interval [a,b] dollars/mile.” This is a complex fuzzy problem which contains two fuzzy factors. To deal with such complex fuzzy problem, a number of complex fuzzy theories have been proposed such as type-2 fuzzy sets [19], level-2 fuzzy sets [20], bifuzzy variable [21], and Fu-Fu variable [22]. In this paper, a Fu-Fu variable is used to model the risk.By comparison to previous researches (see Table1), the effect of the emergency response service is first considered in this paper. Thus, here, a new model for hazmat transportation network design is presented which considers an emergency response teams location problem. In this problem, there are three considered actors: the authority, the carriers, and the emergency response department. The government authority hopes to reduce the hazmats transportation risk by designing an effective transportation network. However, it must also consider the decisions of the carriers and the emergency response department which are closely involved with the transportation risk. Hence, this is considered as a bilevel problem. On the upper level, the authority designs the transportation network. On the lower level, the carriers and emergency response department select the routes and locates the emergency service teams, respectively. In contrast to prior studies, the uncertainty of transportation risk is also explicitly considered in the mathematical model. The remainder of this paper is organized as in Figure 1.Table 1
Comparison among these network design models for transportation of hazardous materials.
Author and year
Model
Decision maker
Objective
Main contribution
Kara and Verter, 2004 [7]
Bi-level integernetwork design model
Leader: local authorityFollower: carriers
Minimum of total riskMinimum of total cost
The transportation network design is first considered as a bi-level model.The bi-level model is generalized to the undirected case.
Erkut and Gzara, 2008 [8]
Bi-level undirectednetwork design model
Bianco et al., 2009 [9]
Bi-level networkflow model
Leader: metalocal authorityFollower: regional authority
Risk equityMinimum of total risk
Multiple layers of authorities and risk equity are first considered.
Our paper
Bi-level uncertain network design model
Leader: local authorityFollower 1: emergency dept. Follower 2: carriers
Minimum of total riskMaximum of coverage Minimum of total cost
Emergency response and uncertainty of risk are first considered in network design model.Figure 1
The structure of the remainder of this paper.


## 2. Problem Statement

Every day around the world large quantities of hazmats are carried by truck. Potential dangers from the hazmats transportation always concern the general public. Therefore, governmental authorities always seek to take measures to reduce and mitigate the risks associated with hazmat transportation. One tool available to authorities is to design a transportation network that restricts hazmat transport to nominated routes [23]. At the same time, the emergency response service network is also important in mitigating total transportation risk. Hence, a hazmat transportation network design problem which considers the emergency response service is discussed in this paper.


### 2.1. Bilevel Problem Description

In the problem presented here, three actors are considered: government authorities, transportation carriers, and emergency response departments. The government authority is often the Traffic and Transport Department. Their main concern is to control the hazmat transportation risk that may be caused to the general population and the environment. This risk is often dependent on both the choice of transportation routes and the available emergency response network. After the authority has designed a transportation network, transportation carriers select transportation routes on this network. They tend to select minimum cost routes rather than those of minimum risk. This choice complicates the hazmat network design problem, since the government authority has to consider the global problem by taking into account all shipments that may pass through its jurisdictional area (Bianco et al. [9]). Emergency response departments are often a fire-fighting department, a first-aid department, a police station, or their union. Though they are all public service departments, there are no leader-follower relationships between the transportation department and the emergency response department. Hence, the emergency response department is independent of the transportation department and in charge of the design of the emergency response network. The response network can be designed by locating accident service teams to specified sites near the highest risk links. This can be solved via a maximal arc-covering location model [10].From the above, the three actors can be seen to have a special relationship. The government authority hopes to mitigate the transportation risk by designing a new network. They know that the risk is dependent on the carriers’ selected routes and the designed emergency response network. However, they do not have the right to impose specific routes on individual carriers or to demand specific service teams locations on emergency response departments and only have the authority to close certain roads to hazmat vehicles. Therefore, the considered problem can be abstracted as a bilevel programming problem as in Figure2. The decision maker on the upper level is the government authority, while the lower level decision makers are the hazmat carriers and the emergency response departments.Figure 2
The abstracted structure of the proposed bilevel hazmat transportation network design problem.


### 2.2. Modelling of Uncertain Transportation Risk

Hazardous materials transportation is an important issue in industrialized societies because of the inherent dangers. What differentiates hazmat shipments from shipments of other materials is the risk associated with an accidental release of these materials during transportation.


#### 2.2.1. Complex Fuzzy Transportation Risk

Most people would agree that risk has to do with the probability and the consequence of an undesirable event. Although some authors define risk as only one of these terms (i.e., probability or consequence), it is more common to define risk as the product of both the probability and consequence of an undesirable event [24]. That is, the transportation risk on link (i,j) is often represented as
(1)rij=pijCij,
where pij is the accident probability on link (i,j) and Cij is the accident consequences.Based on this equation, Erkut and Verter [1] proposed a risk model that takes into account the dependency on the impedances of preceding road segments. It fact, the above equation supposes the existence of the probabilities and consequences of an accident occurring on a route section. Usually, incident probabilities can be estimated using historical frequencies analysis and the logical diagrams methods. Accident consequences are often modelled as a function of the impact area and population, property, and environmental assets within the impact area. Erkut et al. [12] summarized the frequency analysis methods and consequence model in detail. The common characteristic of these methods is the demand for historical/characteristic data. However, it is often difficult to calculate this risk because the precise accident probability and the consequence of an accident are not known as a result of insufficient information. An accident probability is often determined from evidence recorded in past or experimental data. Unfortunately, too many factors impact the probability of an accident such as the volume of traffic, the air exchange rate, the type of hazardous material, and the drivers’ skill, so it is difficult to determine the precise accident probability in any given road section from experiments or past records. Therefore, most often, accident probability is based on subjective management judgment and it is vaguely expressed. For example, it can be said that “the accident probability at X place is about 4×10-5 per mile.” With this information a fuzzy set can be used to describe the uncertain accident probability. The fuzzy probability results in a fuzzy risk, which can be described as “it is about 4×10-5×C dollars/mile.”Let the fuzzy accident probabilityp~ij be given by a domain Up and p~ij={(x,μp~ij(x))∣x∈Up}. Then the fuzzy risk can be given as r~ij,μr~ij(pijCij)=μp~ij(pij), p∈Up, where μp~ij and μr~ij are the membership functions of fuzzy accident probability p~ij and fuzzy transportation risk rijk, respectively.On the other hand, one accident may result in a variety of impacts which need to be factored into the accident consequences such as number of fatalities, size of economic losses, damage to road network, and the effect on the population. Most existing research only looks at the effect on the population. In many cases, such as in construction project transportation, all consequences which affect project cost must be considered. However, for many consequences it is difficult to give a precise evaluation. Manager prefers to give a statement such as “the cost of the consequence of an accident is between 1 and 2 million dollars, and the most likely value is 1.5 million dollars.” This is an example of imprecise information in the decision-making process, which can be translated into a triangular fuzzy number[1,1.5,2] million dollars. Boulmakoul [17] presented a fuzzy approach which uses a fuzzy set to model uncertain consequences. If we consider accident probability and consequences at the same time, then the risk can be described as “it is possible that the transportation risk is [4,6,8] dollars/mile.”Assume we considerN types of fuzzy consequence C~ijn with membership function μC~ijn, n=1,2,…,N. Then the risk for a given pij∈Up will result in fuzzy values:
(2)r~ij=∑n=1N‍C~ijnpij,pij∈Up
with membership function
(3)μr~ij(∑n=1NC~ijnpij)=μp~ij(pij),pij∈Up.
This is a complex fuzzy variable, namely, a fuzzy variable with fuzzy parameters, also called Fu-Fu variable. Therefore, this risk can be modelled as a Fu-Fu variable which means that it has fuzzy values and there are corresponding membership degrees of the risk taking these fuzzy values. Actually, this type of complex fuzzy variables has been applied to some important fields such as database modelling [25], inventory management [26, 27], and vendor selection [28]. By comparison to these quantitative risk analysis methods involved in Erkut et al. [12], the complex fuzzy variables do need less empirical or historical data.


#### 2.2.2. Data Fuzzification Method

Often there is little historical data to describe the accident probability and consequences. Hence, a complex fuzzy variable is proposed to model the risk. Here, how to obtain the complex fuzzy transportation risk from insufficient data using a fuzzification method is introduced. The essence of fuzzification is to find an approximate membership function to describe the fuzzy number [29].Transportation risk is made up of accident probability and accident consequences. In order to determine the membership function of the two types of parameters, a fuzzy evaluation method is proposed. The fuzzy evaluation has been used in many areas such as performance management [30] and students’ evaluation [31]. With accident probability, for example, it is difficult to assign a determined value to each link because there are too many influencing factors. However, it is easy to evaluate these links for a certain impact factor using some fuzzy linguistic term such as low or high. So an evaluation term set is first given; then some experts are invited to give fuzzy evaluation for each impact factor and generate a fuzzy evaluation matrix. And then, the weight for these impact factors is calculated using analytic hierarchy process method. By a fuzzy operation, the evaluation matrix and the weight can be integrated into a set of membership grades of the probability. Finally, each linguistic term for accident probability is modelled as a fuzzy number estimated using historical data, and the final fuzzy result are calculated by the fuzzy product of the menbership grades of the probability and the fuzzy numbers associated with these comment terms. For example, an interval [pmin,pmax] can be determined from a historical frequency analysis. Based on this interval, five different subintervals can be defined to describe five linguistic terms {very  low,  low,  medium,  high,  very  high}. Each interval can be modelled as a fuzzy number such as triangle or a discrete fuzzy number. An example is given in Figure 3. In conclusion, the main fuzzification steps for accident probability are shown as follows.Figure 3
The fuzzification of accident probability linguistic terms.Main fuzzification steps for accident probability are as follows.Step 1.
Selectm impacts as a set U={u1,u2,…,um}. Define a comment linguistic term set for these impacts of accident probability: V={verylow,low,medium,high,veryhigh}.Step 2.
Determine the weight of these impacts using analytic hierarchy process method,W=[w1,w2,…,wm].Step 3.
Calculate the fuzzy relation matrixΨ=(ψij)m×n, ψij=nij/ni, where ni is the number of experts invited to nominate their evaluation to a specific term for impact factor i and nij is the nominated number of linguistic terms j for factor i.Step 4.
Calculate the linguistic evaluation results:B=W∘Ψ, where “∘” represents the fuzzy operator.Step 5.
Describe these linguistic terms with fuzzy numbersHP~=[hp~1,hp~2,…,hp~J], and calculate the fuzzy results p~=B∘HP~, where hp~j is the fuzzy number associated with comment term vj, which can be derived from historical data on accident probabilities.In comparison to accident probability, accident consequences are more difficult to estimate. Usually, accident consequences are composed of injuries and fatalities, property damage, traffic incident delays, and environmental damage. Some consequences can be estimated easily according to available observation data. Taking fatalities as an example, it is common to assume that fatality consequences are proportional to the size of the population in the neighborhood of the considered road link. Hence the consequences cost can be obtained directly through the product of population, mortality rate, and unit cost. The proportion and unit cost can be estimated from historical information. Other consequences, such as environmental damage, are difficult to estimate from observation data. For this type of consequences, the same method as accident probability is used to evaluate the consequence for each link. Then, the sum of all these consequences is the final accident consequence.If the accident probabilityp~ is modelled using a discrete fuzzy number and the consequence C~n is modelled using a triangular fuzzy number C~ijn=[Cln,Cmn,Crn], then the risk can be described as
(4)r~~=p~[Cl,Cm,Cr]={[Clp1,Cmp1,Crp1],withmembershipgradeofμp~(p1),[Clp2,Cmp2,Crp2],withmembershipgradeofμp~(p2),……[Clpn,Cmpn,Crpn],withmembershipgradeofμp~(pn).


## 2.1. Bilevel Problem Description

In the problem presented here, three actors are considered: government authorities, transportation carriers, and emergency response departments. The government authority is often the Traffic and Transport Department. Their main concern is to control the hazmat transportation risk that may be caused to the general population and the environment. This risk is often dependent on both the choice of transportation routes and the available emergency response network. After the authority has designed a transportation network, transportation carriers select transportation routes on this network. They tend to select minimum cost routes rather than those of minimum risk. This choice complicates the hazmat network design problem, since the government authority has to consider the global problem by taking into account all shipments that may pass through its jurisdictional area (Bianco et al. [9]). Emergency response departments are often a fire-fighting department, a first-aid department, a police station, or their union. Though they are all public service departments, there are no leader-follower relationships between the transportation department and the emergency response department. Hence, the emergency response department is independent of the transportation department and in charge of the design of the emergency response network. The response network can be designed by locating accident service teams to specified sites near the highest risk links. This can be solved via a maximal arc-covering location model [10].From the above, the three actors can be seen to have a special relationship. The government authority hopes to mitigate the transportation risk by designing a new network. They know that the risk is dependent on the carriers’ selected routes and the designed emergency response network. However, they do not have the right to impose specific routes on individual carriers or to demand specific service teams locations on emergency response departments and only have the authority to close certain roads to hazmat vehicles. Therefore, the considered problem can be abstracted as a bilevel programming problem as in Figure2. The decision maker on the upper level is the government authority, while the lower level decision makers are the hazmat carriers and the emergency response departments.Figure 2
The abstracted structure of the proposed bilevel hazmat transportation network design problem.


## 2.2. Modelling of Uncertain Transportation Risk

Hazardous materials transportation is an important issue in industrialized societies because of the inherent dangers. What differentiates hazmat shipments from shipments of other materials is the risk associated with an accidental release of these materials during transportation.


### 2.2.1. Complex Fuzzy Transportation Risk

Most people would agree that risk has to do with the probability and the consequence of an undesirable event. Although some authors define risk as only one of these terms (i.e., probability or consequence), it is more common to define risk as the product of both the probability and consequence of an undesirable event [24]. That is, the transportation risk on link (i,j) is often represented as
(1)rij=pijCij,
where pij is the accident probability on link (i,j) and Cij is the accident consequences.Based on this equation, Erkut and Verter [1] proposed a risk model that takes into account the dependency on the impedances of preceding road segments. It fact, the above equation supposes the existence of the probabilities and consequences of an accident occurring on a route section. Usually, incident probabilities can be estimated using historical frequencies analysis and the logical diagrams methods. Accident consequences are often modelled as a function of the impact area and population, property, and environmental assets within the impact area. Erkut et al. [12] summarized the frequency analysis methods and consequence model in detail. The common characteristic of these methods is the demand for historical/characteristic data. However, it is often difficult to calculate this risk because the precise accident probability and the consequence of an accident are not known as a result of insufficient information. An accident probability is often determined from evidence recorded in past or experimental data. Unfortunately, too many factors impact the probability of an accident such as the volume of traffic, the air exchange rate, the type of hazardous material, and the drivers’ skill, so it is difficult to determine the precise accident probability in any given road section from experiments or past records. Therefore, most often, accident probability is based on subjective management judgment and it is vaguely expressed. For example, it can be said that “the accident probability at X place is about 4×10-5 per mile.” With this information a fuzzy set can be used to describe the uncertain accident probability. The fuzzy probability results in a fuzzy risk, which can be described as “it is about 4×10-5×C dollars/mile.”Let the fuzzy accident probabilityp~ij be given by a domain Up and p~ij={(x,μp~ij(x))∣x∈Up}. Then the fuzzy risk can be given as r~ij,μr~ij(pijCij)=μp~ij(pij), p∈Up, where μp~ij and μr~ij are the membership functions of fuzzy accident probability p~ij and fuzzy transportation risk rijk, respectively.On the other hand, one accident may result in a variety of impacts which need to be factored into the accident consequences such as number of fatalities, size of economic losses, damage to road network, and the effect on the population. Most existing research only looks at the effect on the population. In many cases, such as in construction project transportation, all consequences which affect project cost must be considered. However, for many consequences it is difficult to give a precise evaluation. Manager prefers to give a statement such as “the cost of the consequence of an accident is between 1 and 2 million dollars, and the most likely value is 1.5 million dollars.” This is an example of imprecise information in the decision-making process, which can be translated into a triangular fuzzy number[1,1.5,2] million dollars. Boulmakoul [17] presented a fuzzy approach which uses a fuzzy set to model uncertain consequences. If we consider accident probability and consequences at the same time, then the risk can be described as “it is possible that the transportation risk is [4,6,8] dollars/mile.”Assume we considerN types of fuzzy consequence C~ijn with membership function μC~ijn, n=1,2,…,N. Then the risk for a given pij∈Up will result in fuzzy values:
(2)r~ij=∑n=1N‍C~ijnpij,pij∈Up
with membership function
(3)μr~ij(∑n=1NC~ijnpij)=μp~ij(pij),pij∈Up.
This is a complex fuzzy variable, namely, a fuzzy variable with fuzzy parameters, also called Fu-Fu variable. Therefore, this risk can be modelled as a Fu-Fu variable which means that it has fuzzy values and there are corresponding membership degrees of the risk taking these fuzzy values. Actually, this type of complex fuzzy variables has been applied to some important fields such as database modelling [25], inventory management [26, 27], and vendor selection [28]. By comparison to these quantitative risk analysis methods involved in Erkut et al. [12], the complex fuzzy variables do need less empirical or historical data.


### 2.2.2. Data Fuzzification Method

Often there is little historical data to describe the accident probability and consequences. Hence, a complex fuzzy variable is proposed to model the risk. Here, how to obtain the complex fuzzy transportation risk from insufficient data using a fuzzification method is introduced. The essence of fuzzification is to find an approximate membership function to describe the fuzzy number [29].Transportation risk is made up of accident probability and accident consequences. In order to determine the membership function of the two types of parameters, a fuzzy evaluation method is proposed. The fuzzy evaluation has been used in many areas such as performance management [30] and students’ evaluation [31]. With accident probability, for example, it is difficult to assign a determined value to each link because there are too many influencing factors. However, it is easy to evaluate these links for a certain impact factor using some fuzzy linguistic term such as low or high. So an evaluation term set is first given; then some experts are invited to give fuzzy evaluation for each impact factor and generate a fuzzy evaluation matrix. And then, the weight for these impact factors is calculated using analytic hierarchy process method. By a fuzzy operation, the evaluation matrix and the weight can be integrated into a set of membership grades of the probability. Finally, each linguistic term for accident probability is modelled as a fuzzy number estimated using historical data, and the final fuzzy result are calculated by the fuzzy product of the menbership grades of the probability and the fuzzy numbers associated with these comment terms. For example, an interval [pmin,pmax] can be determined from a historical frequency analysis. Based on this interval, five different subintervals can be defined to describe five linguistic terms {very  low,  low,  medium,  high,  very  high}. Each interval can be modelled as a fuzzy number such as triangle or a discrete fuzzy number. An example is given in Figure 3. In conclusion, the main fuzzification steps for accident probability are shown as follows.Figure 3
The fuzzification of accident probability linguistic terms.Main fuzzification steps for accident probability are as follows.Step 1.
Selectm impacts as a set U={u1,u2,…,um}. Define a comment linguistic term set for these impacts of accident probability: V={verylow,low,medium,high,veryhigh}.Step 2.
Determine the weight of these impacts using analytic hierarchy process method,W=[w1,w2,…,wm].Step 3.
Calculate the fuzzy relation matrixΨ=(ψij)m×n, ψij=nij/ni, where ni is the number of experts invited to nominate their evaluation to a specific term for impact factor i and nij is the nominated number of linguistic terms j for factor i.Step 4.
Calculate the linguistic evaluation results:B=W∘Ψ, where “∘” represents the fuzzy operator.Step 5.
Describe these linguistic terms with fuzzy numbersHP~=[hp~1,hp~2,…,hp~J], and calculate the fuzzy results p~=B∘HP~, where hp~j is the fuzzy number associated with comment term vj, which can be derived from historical data on accident probabilities.In comparison to accident probability, accident consequences are more difficult to estimate. Usually, accident consequences are composed of injuries and fatalities, property damage, traffic incident delays, and environmental damage. Some consequences can be estimated easily according to available observation data. Taking fatalities as an example, it is common to assume that fatality consequences are proportional to the size of the population in the neighborhood of the considered road link. Hence the consequences cost can be obtained directly through the product of population, mortality rate, and unit cost. The proportion and unit cost can be estimated from historical information. Other consequences, such as environmental damage, are difficult to estimate from observation data. For this type of consequences, the same method as accident probability is used to evaluate the consequence for each link. Then, the sum of all these consequences is the final accident consequence.If the accident probabilityp~ is modelled using a discrete fuzzy number and the consequence C~n is modelled using a triangular fuzzy number C~ijn=[Cln,Cmn,Crn], then the risk can be described as
(4)r~~=p~[Cl,Cm,Cr]={[Clp1,Cmp1,Crp1],withmembershipgradeofμp~(p1),[Clp2,Cmp2,Crp2],withmembershipgradeofμp~(p2),……[Clpn,Cmpn,Crpn],withmembershipgradeofμp~(pn).


## 2.2.1. Complex Fuzzy Transportation Risk

Most people would agree that risk has to do with the probability and the consequence of an undesirable event. Although some authors define risk as only one of these terms (i.e., probability or consequence), it is more common to define risk as the product of both the probability and consequence of an undesirable event [24]. That is, the transportation risk on link (i,j) is often represented as
(1)rij=pijCij,
where pij is the accident probability on link (i,j) and Cij is the accident consequences.Based on this equation, Erkut and Verter [1] proposed a risk model that takes into account the dependency on the impedances of preceding road segments. It fact, the above equation supposes the existence of the probabilities and consequences of an accident occurring on a route section. Usually, incident probabilities can be estimated using historical frequencies analysis and the logical diagrams methods. Accident consequences are often modelled as a function of the impact area and population, property, and environmental assets within the impact area. Erkut et al. [12] summarized the frequency analysis methods and consequence model in detail. The common characteristic of these methods is the demand for historical/characteristic data. However, it is often difficult to calculate this risk because the precise accident probability and the consequence of an accident are not known as a result of insufficient information. An accident probability is often determined from evidence recorded in past or experimental data. Unfortunately, too many factors impact the probability of an accident such as the volume of traffic, the air exchange rate, the type of hazardous material, and the drivers’ skill, so it is difficult to determine the precise accident probability in any given road section from experiments or past records. Therefore, most often, accident probability is based on subjective management judgment and it is vaguely expressed. For example, it can be said that “the accident probability at X place is about 4×10-5 per mile.” With this information a fuzzy set can be used to describe the uncertain accident probability. The fuzzy probability results in a fuzzy risk, which can be described as “it is about 4×10-5×C dollars/mile.”Let the fuzzy accident probabilityp~ij be given by a domain Up and p~ij={(x,μp~ij(x))∣x∈Up}. Then the fuzzy risk can be given as r~ij,μr~ij(pijCij)=μp~ij(pij), p∈Up, where μp~ij and μr~ij are the membership functions of fuzzy accident probability p~ij and fuzzy transportation risk rijk, respectively.On the other hand, one accident may result in a variety of impacts which need to be factored into the accident consequences such as number of fatalities, size of economic losses, damage to road network, and the effect on the population. Most existing research only looks at the effect on the population. In many cases, such as in construction project transportation, all consequences which affect project cost must be considered. However, for many consequences it is difficult to give a precise evaluation. Manager prefers to give a statement such as “the cost of the consequence of an accident is between 1 and 2 million dollars, and the most likely value is 1.5 million dollars.” This is an example of imprecise information in the decision-making process, which can be translated into a triangular fuzzy number[1,1.5,2] million dollars. Boulmakoul [17] presented a fuzzy approach which uses a fuzzy set to model uncertain consequences. If we consider accident probability and consequences at the same time, then the risk can be described as “it is possible that the transportation risk is [4,6,8] dollars/mile.”Assume we considerN types of fuzzy consequence C~ijn with membership function μC~ijn, n=1,2,…,N. Then the risk for a given pij∈Up will result in fuzzy values:
(2)r~ij=∑n=1N‍C~ijnpij,pij∈Up
with membership function
(3)μr~ij(∑n=1NC~ijnpij)=μp~ij(pij),pij∈Up.
This is a complex fuzzy variable, namely, a fuzzy variable with fuzzy parameters, also called Fu-Fu variable. Therefore, this risk can be modelled as a Fu-Fu variable which means that it has fuzzy values and there are corresponding membership degrees of the risk taking these fuzzy values. Actually, this type of complex fuzzy variables has been applied to some important fields such as database modelling [25], inventory management [26, 27], and vendor selection [28]. By comparison to these quantitative risk analysis methods involved in Erkut et al. [12], the complex fuzzy variables do need less empirical or historical data.


## 2.2.2. Data Fuzzification Method

Often there is little historical data to describe the accident probability and consequences. Hence, a complex fuzzy variable is proposed to model the risk. Here, how to obtain the complex fuzzy transportation risk from insufficient data using a fuzzification method is introduced. The essence of fuzzification is to find an approximate membership function to describe the fuzzy number [29].Transportation risk is made up of accident probability and accident consequences. In order to determine the membership function of the two types of parameters, a fuzzy evaluation method is proposed. The fuzzy evaluation has been used in many areas such as performance management [30] and students’ evaluation [31]. With accident probability, for example, it is difficult to assign a determined value to each link because there are too many influencing factors. However, it is easy to evaluate these links for a certain impact factor using some fuzzy linguistic term such as low or high. So an evaluation term set is first given; then some experts are invited to give fuzzy evaluation for each impact factor and generate a fuzzy evaluation matrix. And then, the weight for these impact factors is calculated using analytic hierarchy process method. By a fuzzy operation, the evaluation matrix and the weight can be integrated into a set of membership grades of the probability. Finally, each linguistic term for accident probability is modelled as a fuzzy number estimated using historical data, and the final fuzzy result are calculated by the fuzzy product of the menbership grades of the probability and the fuzzy numbers associated with these comment terms. For example, an interval [pmin,pmax] can be determined from a historical frequency analysis. Based on this interval, five different subintervals can be defined to describe five linguistic terms {very  low,  low,  medium,  high,  very  high}. Each interval can be modelled as a fuzzy number such as triangle or a discrete fuzzy number. An example is given in Figure 3. In conclusion, the main fuzzification steps for accident probability are shown as follows.Figure 3
The fuzzification of accident probability linguistic terms.Main fuzzification steps for accident probability are as follows.Step 1.
Selectm impacts as a set U={u1,u2,…,um}. Define a comment linguistic term set for these impacts of accident probability: V={verylow,low,medium,high,veryhigh}.Step 2.
Determine the weight of these impacts using analytic hierarchy process method,W=[w1,w2,…,wm].Step 3.
Calculate the fuzzy relation matrixΨ=(ψij)m×n, ψij=nij/ni, where ni is the number of experts invited to nominate their evaluation to a specific term for impact factor i and nij is the nominated number of linguistic terms j for factor i.Step 4.
Calculate the linguistic evaluation results:B=W∘Ψ, where “∘” represents the fuzzy operator.Step 5.
Describe these linguistic terms with fuzzy numbersHP~=[hp~1,hp~2,…,hp~J], and calculate the fuzzy results p~=B∘HP~, where hp~j is the fuzzy number associated with comment term vj, which can be derived from historical data on accident probabilities.In comparison to accident probability, accident consequences are more difficult to estimate. Usually, accident consequences are composed of injuries and fatalities, property damage, traffic incident delays, and environmental damage. Some consequences can be estimated easily according to available observation data. Taking fatalities as an example, it is common to assume that fatality consequences are proportional to the size of the population in the neighborhood of the considered road link. Hence the consequences cost can be obtained directly through the product of population, mortality rate, and unit cost. The proportion and unit cost can be estimated from historical information. Other consequences, such as environmental damage, are difficult to estimate from observation data. For this type of consequences, the same method as accident probability is used to evaluate the consequence for each link. Then, the sum of all these consequences is the final accident consequence.If the accident probabilityp~ is modelled using a discrete fuzzy number and the consequence C~n is modelled using a triangular fuzzy number C~ijn=[Cln,Cmn,Crn], then the risk can be described as
(4)r~~=p~[Cl,Cm,Cr]={[Clp1,Cmp1,Crp1],withmembershipgradeofμp~(p1),[Clp2,Cmp2,Crp2],withmembershipgradeofμp~(p2),……[Clpn,Cmpn,Crpn],withmembershipgradeofμp~(pn).


## 3. Modelling

The hazmat network design problem is a graph theoretical problem defined on a directed graphG=(V,A), where V is the set of vertices and A is the set of arcs on the graph. A vertex corresponds to a road intersection, and an arc corresponds to a road segment on the network. The network design problem finds a network to transport K commodities between their respective origins and destinations. Each commodity corresponds to an OD pair. Let (s(k),t(k)) be the OD pair of commodity k,k∈{1⋯K} and let dk be the corresponding number of shipments. The parameters rijk and cijk refer to the risk and cost associated with a unit flow of commodity k on arc (i,j), respectively. Each link is assumed to be a network segment.


### 3.1. Fu-Fu Variable

Considering the lack of historical data used to describe the accident probability and consequences, the transportation risks are modelled as Fu-Fu variables. Some basic knowledge about Fu-Fu variable is introduced as follows.Definition 1 (see [22]).
A Fu-Fu variableξ is a fuzzy variable with fuzzy parameters.Example 2.
Letξ~1,ξ~2,…,ξ~n be triangular fuzzy number and let μ1,μ2,…,μn be real numbers in [0,1] such that μ1∨μ2∨⋯∨μn=1. Then
(5)ξ={ξ~1withmembershipμ1ξ~2withmembershipμ2…ξ~nwithmembershipμn
is a Fu-Fu variable.Example 3.
ξ
~
~
=
(
L
,
γ
~
,
R
) with γ~~(γL,γM,γR) is called Fu-Fu variable, see Figure 4, if the outer-layer and inner-layer membership functions are as follows:
(6)μξ~~(x)={(x-L)γ~-L,ifL≤x≤γ~(R-x)R-γ~,ifγ~≤x≤R0,otherwise,μγ~(x′)={(x-γL)γM-γL,ifγL≤x≤γM(γR-x)γR-γM,ifγM≤x≤γR0,otherwise,
where γ~ is the center of ξ~~, which is a triangular fuzzy variable, and L and R are the smallest possible value and the largest possible value of ξ~~. γL, γM, and γR are the smallest possible value, the most promising value, and the largest possible value of γ~, respectively.Figure 4
A triangular Fu-Fu variable.Definition 4 (see [22]).
The expected value of a Fu-Fu variable is defined by(7)E[ξ]=∫0+∞Cr{θ∈Θ∣E[ξ(θ)]≥r}dr-∫-∞0Cr{θ∈Θ∣E[ξ(θ)]≤r}dr
provided that at least one of the two integrals is finite.Theorem 5 (see [22]).
Assume thatξ and η are Fu-Fu variables with finite expected values. If (i) for each θ∈Θ, the fuzzy variables ξ(θ) and η(θ) are independent, and (ii) E[ξ(θ)] and E[η(θ)] are independent fuzzy variables, then for any real numbers a and b, one has
(8)E[aξ+bη]=aE[ξ]+bE[η].Lemma 6.
If transportation riskr~~ijk is a Fu-Fu variable characterized as follows (see Figure 5):
(9)r~~ijk={[rijk11,rijk12,rijk13],withmembershipgrade  ofμijk1[rijk21,rijk22,rijk23],withmembershipgrade  ofμijk2……[rijkn1,rijkn2,rijkn3],withmembershipgrade  ofμijkn,
where [rijkm1,rijkm2,rijkm3] is a triangular fuzzy number, μijk1 is the degree of membership associated with r~~ijk taking fuzzy value [rijk11,rijk12,rijk13], and rijk11, rijk12, and rijk13 are the smallest possible value, the most promising value, and the largest possible value of the triangular fuzzy number, respectively, then the expected value of r~~ijk is
(10)E[r~~ijk]=14∑t=1n‍(rijkt1+2rijkt2+rijkt3)wijkt,
where the weights wijkt,t=1,2,…,n are given by
(11)wijk1=12(μijk1+max1≤m≤nμijkm-max1<m≤nμijkm),wijkl=12(max1≤m≤lμijkm-max1≤m<lμijkm+maxl≤m≤nμijkm-maxl<m≤nμijkm),2≤l<n,wijkn=12(max1≤m≤nμijkm-max1≤m<nμijkm+μijkn).Figure 5
An illustration of Fu-Fu fuzzy risk.Proof.
For anym=1,2,…,n, r~~ijk(xm)=[rijkm1,rijkm2,rijkm3] is a triangular fuzzy variable. It follows from the definition of expected value of fuzzy variable that we have
(12)E[r~~ijk(xm)]=14(rijkm1+2rijkm2+rijkm3).
From Definition 4 and Theorem 5,
(13)E[r~~ijk]=E[E[r~~ijk(x)]]=E[14(r~ijk1(x)+2r~ijk2(x)+r~ijk3(x))]=14(E[r~ijk1(x)]+2E[r~ijk2(x)]+E[r~ijk3(x)]).
Here, r~ijks(x)(s=1,2,3) is a discrete fuzzy variable whose membership function is given by
(14)μr~ijks(x)(θ)={μijk1,ifθ=rijk1sμijk2,ifθ=rijk2s……μijkn,ifθ=rijkns.
From the definition of expected value of fuzzy variable, we have
(15)E[r~ijks]=∑t=1n‍(rijktswijkt),s=1,2,3,
where
(16)wijk1=12(μijk1+max1≤m≤nμijkm-max1<m≤nμijkm),wijkl=12(max1≤m≤lμijkm-max1≤m<lμijkm+maxl≤m≤nμijkm-maxl<m≤nμijkm),2≤l<n,wijkn=12(max1≤m≤nμijkm-max1≤m<nμijkm+μijkn).
And then,
(17)E[r~~ijk]=14(E[r~ijk1(x)]+2E[r~ijk2(x)]+E[r~ijk3(x)])=14∑t=1n‍(rijkt1+2rijkt2+rijkt3)wijkt.
This completes the proof.


### 3.2. A Network Design Model

The problem the government authority on the upper level faces is how to design a transportation network which minimizes the total risk of the shipments; in other words they need to decide whether it is to be used to transport the hazardous materials to each road section. With this in mind, the decision variables for the upper level areyij.In most transport planning models, the objective is to move products from the origins to the destinations at minimal cost. However, for hazmat shipments, a cost-minimizing objective is usually not suitable. The risk associated with hazmats means that the problems are more complicated (and more interesting) than many other transport problems. The total risk of the network is the sum of the risks of each link. The risk of each link is dependent on the number of shipments, the unit risk. It also relies on whether it is covered by emergency response teams. An emergency response team is often made up of various emergency response facilities and staffs. Many tasks such as fire fighting, ambulance and police services, and hazmat containment and clean-up involve the emergency response teams. These activities can have a positive effect on most accident consequences. Hence, effective emergency response is crucial to contain the impact on the smallest possible area and mitigate undesirable consequences when a hazmat incident occurs [10, 11]. All of these produce an effect on link risk.Hamouda et al. [32] developed a risk assessment model which considers emergency response. They assumed that risk is reduced if a demand node/link can be responded to by an emergency team. Moreover, the reduced risk also depends on the type of material transported and the travel distance from the accident site to the response team location. It is reasonable to assume that the response teams always drive to the accident sites along a shortest path, and let the midpoint be the concentrated point of a road link. Then the travel distance from a response node q to a demand link (i,j) can be described using the shortest distance from the midpoint of link (i,j) to node q, that is, Dijq. Meanwhile, this distance cannot exceed the maximum service distance Dmax of the emergency response teams. Therefore, if link (i,j) is covered by node q, that is, hijq=1, then the reduced risk can be described as a function of the distance Dijq and Dmax. For convenience, it is assumed that the reduced risk and the travel distance meet a linear relation such as α(Dmax-Dijq)hijqrijk.Here,α(Dmax-Dijq)hijq describes the percentage of reduced risk when link (i,j) is covered by emergency response team located in node q. α is a coefficient related to the category of undesirable accident consequences and the power of the emergency response teams. αDmax reflects the maximum power of the emergency response teams to service a hazardous accident, that is, how much the accident consequences can be reduced when an accident is serviced by a very timely emergency team. The maximum service distance Dmax can be obtained from the experience of the emergency response teams. The parameter α can be estimated by an analysis of the category of undesirable accident consequences and the power of the emergency response teams as outlined in Figure 6. First, the service ability of the emergency teams should be evaluated. Then, the undesirable consequences are classified and their weights are determined based on the analysis of past accidents. Next, for each category of accident consequences, estimate a possible range of the decrease if an accident is serviced by a very timely emergency response team. The median values of these ranges are taken as the most possible values. Finally, the weight sum of these values is determined as the final value of parameter α. In fact, it is very difficult to get an accurate value of α because it is also dependent on the managers’ attitude. However, a sensitivity analysis on α can assist in the managers’ decision making. Thus, the total uncertain transportation risk can be described as
(18)∑k=1K‍∑(i,j)∈A‍∑q∈H‍dkr~~ijklijxijk(1-α(Dmax-Dijq)hijq).Figure 6
Estimation of the value ofα.In the risk described previously,r~~ijk is a Fu-Fu variable, so the total risk is also a Fu-Fu variable. However, it is difficult to make a decision when it involves uncertain information, so it is necessary to transform the Fu-Fu risk to a determinate one. In this case, the authorities tend to design a network with minimal expected risk. That is, the Fu-Fu risk can be transformed to a determinate one by an expected value operation; the expected total transportation risk can be described as follows:
(19)Risk=E[∑k=1K‍∑(i,j)∈A‍∑q∈H‍dkr~~ijklijxijkRisk=E[∑k=1K‍∑(i,j)∈A‍∑q∈H‍×(1-α(Dmax-Dijq)hijq)∑k=1K‍∑(i,j)∈A‍∑q∈H‍].In the expected total risk, onlyr~~ijk is a Fu-Fu variable. Based on Definition 4 and Theorem 5, the expected total risk can be transformed into following equation:
(20)Risk=∑k=1K‍∑(i,j)∈A‍∑q∈H‍dkE[r~~ijk]lijxijk×(1-α(Dmax-Dijq)hijq).For this objective, the variablesxijk and hij are solved in the lower level model. Hence, the network design model can be modelled as follows:
(21)minyijRisk=∑k=1K‍∑(i,j)∈A‍∑q∈H‍dkE[r~~ijk]lijxijkminyijRisk=∑k=1K‍∑(i,j)∈A‍∑q∈H‍×(1-α(Dmax-Dijq)hijq)s.t.{yij=yji,yij∈{0,1}wherexijkandhijaresolvedinthelower-levelmodel.


### 3.3. A Network Flow Model

The main concern of a government authority is to minimize the total transportation risk by designing a new transportation network. However, risk is also affected by the carriers’ transportation routes and the location of the emergency response teams. The hazmat carriers choose their routes from the upper-level designed network. Thus, this is a multicommodity network flow assignment problem. To these carriers, the minimization of their total transportation cost is a primary objective which conflicts with that of the government authority as cost is not their primary concern. Hence, the flow model should be considered as a separate model, rather than being merged into the network design model. Therefore, the authority faces a bilevel decision problem. That is, while designing transportation network, the authority must also consider the actual use of the hazmat network by the carriers and emergency response teams. Erkut and Gzara [8] also carried out a comparison analysis to explain why it is better to consider the hazmat network design problem as a bilevel model.In the model, the total cost can be described as∑k=1K∑(i,j)∈Adkcijklijxijk. For the flow problem, two constraints must be met. One is the flow equilibrium constraint which ensures the flow of commodity k from its origin to its destination:
(22)∑(i,j)∈A‍xijk-∑(i,j)∈A‍xjik={-1,j=s(k)1,j=t(k),∀j∈V,k=1,2,…,K0,otherwise.The others are the road capacity constraint, which ensures that only the routes selected by the government can be used by the carriers, and the logic constraint, which ensures that the variables only take a binary value:(23)xijk≤yij,xijk∈{0,1},∀(i,j)∈A,k=1,2,…,K.The objective and constraints constitute the routes choice model of the carriers as follows:(24)minxijkCost=∑k=1K‍∑(i,j)∈A‍dkcijklijxijks.t.{∑(i,j)∈A‍xijk-∑(i,j)∈A‍xjik={-1,j=s(k)1,j=t(k),∀j∈V,k=1,2,…,K0,otherwisexijk≤yij,xijk∈{0,1},∀(i,j)∈A,k=1,2,…,K.


### 3.4. A Maximal Arc-Covering Location Model

Usually, a transportation network is designed by a local traffic and transportation department. However, the location of emergency response teams is often chosen by the emergency response departments such as the fire department, the first-aid department, the police office, or their union. The location of these emergency response teams plays an important role in mitigating transportation risk, so the location is considered to be a decision on the lower level when designing the transportation network. The emergency response departments decide the teams location based on the road network design and the carriers’ chosen routes. They hope to cover all the road links, so maximizing the total weighted arc length covered is the objective. This is a maximal covering location problem for hazardous materials transportation. The maximal covering location model was originally developed by Church and ReVelle [33] and was used to locate hazmat response teams in [10, 32].Here, the weight of link(i,j) describes the gain when link (i,j) is covered by emergency response teams, that is, the reduced risk ∑q∈H∑k=1Kdklijα(Dmax-Dijq)r~~ijkxijk. The expected value operation is used to deal with the Fu-Fu risk r~~ijk; where hijq takes value 1 if (i,j) is covered by emergency response team located at note q and 0 otherwise. Then the objective can be modelled as follows:
(25)Cover=∑(i,j)∈A‍∑q∈H‍∑k=1K‍dklijα(Dmax-Dijq)E[r~~ijk]xijkhijq.At the same time, two constraints must also be considered. First, all located teams must be equal to its sum. Letzq describe whether a team is located in node q, and fmax describes the total number of located emergency response teams; then the constraint can be stated as follows:
(26)∑q∈H‍zq=fmax.Moreover, a link(i,j) is covered by node q only if a team is located at node q; then this constraint can be stated as follows:
(27)hijq≤zq.Meanwhile, in most cases, when an accident occurs, only one nearest emergency response team is arranged to service. Therefore, it is assumed that each link is only serviced by a single response team:(28)∑q∈H‍hijq≤1.Finally, the maximum travel distance of an emergency response team to any link traversed by hazardous materials within its jurisdiction should not exceed the thresholdDmax. It can be described as the following constraint:
(29)∑q∈H‍Dijqhijq≤Dmax.From the aforementioned, the maximal arc-covering location problem can be modelled as(30)maxzq,hijqCover=∑(i,j)∈A‍∑q∈H‍∑k=1K‍dklijαmaxzq,hijqCover=∑(i,j)∈A‍∑q∈H‍∑k=1K‍y×(Dmax-Dijq)E[r~~ijk]xijkhijqs.t.{∑q∈H‍zq=fmaxhijq≤zq,∀(i,j)∈A,q∈H∑q∈H‍hijq≤1,∀(i,j)∈A∑q∈H‍Dijqhijq≤Dmax,∀(i,j)∈Azq,hijq∈{0,1}.


### 3.5. Bilevel Programming Model

As outlined in Section2, the considered problem should be modelled as a bilevel programming model. The decision maker on the upper level is the authority who hopes to mitigate the hazmats transportation risk by designing a new road network. They know that the risk is dependent on the carriers’ selected routes and the designed emergency response network which they do not have the right to determine. Fortunately, they also know that the carriers and the emergency response department will make decisions based on the designed network. Therefore, the authority can consider their decision making and influence the lower-level model. On the lower level, the carriers first choose their routes so that total transportation cost is minimized. Then, the emergency response department locates their emergency service teams with the objective of maximizing the total weighted arc length covered. Hence, the complete bilevel programming model can be established based on the previous discussion as follows:
(31)minyijRisk=∑k=1K‍∑(i,j)∈A‍∑q∈H‍dkE[r~~ijk]lijxijkminyijRisk=∑k=1K‍∑(i,j)∈A‍∑q∈H‍×(1-α(Dmax-Dijq)hijq)s.t.{yij=yji,yij∈{0,1}maxzq,hijqCover=∑(i,j)∈A∑q∈H∑k=1Kdklijαmaxzq,hijqCover=∑∑∑∑∑×(Dmax-Dijq)E[r~~ijk]maxzq,hijqCover=∑∑∑∑∑×xijkhijqs.t.{∑q∈H‍zq=fmaxhijq≤zq,∀(i,j)∈A,q∈H∑q∈H‍hijq≤1,∀(i,j)∈A∑q∈H‍Dijqhijq≤Dmax,∀(i,j)∈Azq,hijq∈{0,1}minxijkCost=∑k=1K‍∑(i,j)∈A‍dkcijklijxijks.t.{∑(i,j)∈A‍xijk-∑(i,j)∈A‍xjik=-1,∑(i,j)∈A‍xijk-xjikj=s(k)∑(i,j)∈A‍xijk-∑(i,j)∈A‍xjik=1,∑(i,j)∈A‍xijk-∑(i,j)∈A‍j=t(k)∑(i,j)∈A‍xijk-∑(i,j)∈A‍xjik=0,∑(i,j)∈A‍xijk-j≠s(k),t(k)xijk≤yijxijk∈{0,1}.In this model, the Fu-Fu risks have been dealt with using expected value operation defined in (7). For some especial types of Fu-Fu variables, the expected risk can be calculated directly. As for others, they can be calculated by some simulation methods as introduced in [22]. In this paper, the risk [r~~ijk] is modelled as a discrete triangular Fu-Fu variable (see (9)), and its expected value can be calculated by (10).


## 3.1. Fu-Fu Variable

Considering the lack of historical data used to describe the accident probability and consequences, the transportation risks are modelled as Fu-Fu variables. Some basic knowledge about Fu-Fu variable is introduced as follows.Definition 1 (see [22]).
A Fu-Fu variableξ is a fuzzy variable with fuzzy parameters.Example 2.
Letξ~1,ξ~2,…,ξ~n be triangular fuzzy number and let μ1,μ2,…,μn be real numbers in [0,1] such that μ1∨μ2∨⋯∨μn=1. Then
(5)ξ={ξ~1withmembershipμ1ξ~2withmembershipμ2…ξ~nwithmembershipμn
is a Fu-Fu variable.Example 3.
ξ
~
~
=
(
L
,
γ
~
,
R
) with γ~~(γL,γM,γR) is called Fu-Fu variable, see Figure 4, if the outer-layer and inner-layer membership functions are as follows:
(6)μξ~~(x)={(x-L)γ~-L,ifL≤x≤γ~(R-x)R-γ~,ifγ~≤x≤R0,otherwise,μγ~(x′)={(x-γL)γM-γL,ifγL≤x≤γM(γR-x)γR-γM,ifγM≤x≤γR0,otherwise,
where γ~ is the center of ξ~~, which is a triangular fuzzy variable, and L and R are the smallest possible value and the largest possible value of ξ~~. γL, γM, and γR are the smallest possible value, the most promising value, and the largest possible value of γ~, respectively.Figure 4
A triangular Fu-Fu variable.Definition 4 (see [22]).
The expected value of a Fu-Fu variable is defined by(7)E[ξ]=∫0+∞Cr{θ∈Θ∣E[ξ(θ)]≥r}dr-∫-∞0Cr{θ∈Θ∣E[ξ(θ)]≤r}dr
provided that at least one of the two integrals is finite.Theorem 5 (see [22]).
Assume thatξ and η are Fu-Fu variables with finite expected values. If (i) for each θ∈Θ, the fuzzy variables ξ(θ) and η(θ) are independent, and (ii) E[ξ(θ)] and E[η(θ)] are independent fuzzy variables, then for any real numbers a and b, one has
(8)E[aξ+bη]=aE[ξ]+bE[η].Lemma 6.
If transportation riskr~~ijk is a Fu-Fu variable characterized as follows (see Figure 5):
(9)r~~ijk={[rijk11,rijk12,rijk13],withmembershipgrade  ofμijk1[rijk21,rijk22,rijk23],withmembershipgrade  ofμijk2……[rijkn1,rijkn2,rijkn3],withmembershipgrade  ofμijkn,
where [rijkm1,rijkm2,rijkm3] is a triangular fuzzy number, μijk1 is the degree of membership associated with r~~ijk taking fuzzy value [rijk11,rijk12,rijk13], and rijk11, rijk12, and rijk13 are the smallest possible value, the most promising value, and the largest possible value of the triangular fuzzy number, respectively, then the expected value of r~~ijk is
(10)E[r~~ijk]=14∑t=1n‍(rijkt1+2rijkt2+rijkt3)wijkt,
where the weights wijkt,t=1,2,…,n are given by
(11)wijk1=12(μijk1+max1≤m≤nμijkm-max1<m≤nμijkm),wijkl=12(max1≤m≤lμijkm-max1≤m<lμijkm+maxl≤m≤nμijkm-maxl<m≤nμijkm),2≤l<n,wijkn=12(max1≤m≤nμijkm-max1≤m<nμijkm+μijkn).Figure 5
An illustration of Fu-Fu fuzzy risk.Proof.
For anym=1,2,…,n, r~~ijk(xm)=[rijkm1,rijkm2,rijkm3] is a triangular fuzzy variable. It follows from the definition of expected value of fuzzy variable that we have
(12)E[r~~ijk(xm)]=14(rijkm1+2rijkm2+rijkm3).
From Definition 4 and Theorem 5,
(13)E[r~~ijk]=E[E[r~~ijk(x)]]=E[14(r~ijk1(x)+2r~ijk2(x)+r~ijk3(x))]=14(E[r~ijk1(x)]+2E[r~ijk2(x)]+E[r~ijk3(x)]).
Here, r~ijks(x)(s=1,2,3) is a discrete fuzzy variable whose membership function is given by
(14)μr~ijks(x)(θ)={μijk1,ifθ=rijk1sμijk2,ifθ=rijk2s……μijkn,ifθ=rijkns.
From the definition of expected value of fuzzy variable, we have
(15)E[r~ijks]=∑t=1n‍(rijktswijkt),s=1,2,3,
where
(16)wijk1=12(μijk1+max1≤m≤nμijkm-max1<m≤nμijkm),wijkl=12(max1≤m≤lμijkm-max1≤m<lμijkm+maxl≤m≤nμijkm-maxl<m≤nμijkm),2≤l<n,wijkn=12(max1≤m≤nμijkm-max1≤m<nμijkm+μijkn).
And then,
(17)E[r~~ijk]=14(E[r~ijk1(x)]+2E[r~ijk2(x)]+E[r~ijk3(x)])=14∑t=1n‍(rijkt1+2rijkt2+rijkt3)wijkt.
This completes the proof.


## 3.2. A Network Design Model

The problem the government authority on the upper level faces is how to design a transportation network which minimizes the total risk of the shipments; in other words they need to decide whether it is to be used to transport the hazardous materials to each road section. With this in mind, the decision variables for the upper level areyij.In most transport planning models, the objective is to move products from the origins to the destinations at minimal cost. However, for hazmat shipments, a cost-minimizing objective is usually not suitable. The risk associated with hazmats means that the problems are more complicated (and more interesting) than many other transport problems. The total risk of the network is the sum of the risks of each link. The risk of each link is dependent on the number of shipments, the unit risk. It also relies on whether it is covered by emergency response teams. An emergency response team is often made up of various emergency response facilities and staffs. Many tasks such as fire fighting, ambulance and police services, and hazmat containment and clean-up involve the emergency response teams. These activities can have a positive effect on most accident consequences. Hence, effective emergency response is crucial to contain the impact on the smallest possible area and mitigate undesirable consequences when a hazmat incident occurs [10, 11]. All of these produce an effect on link risk.Hamouda et al. [32] developed a risk assessment model which considers emergency response. They assumed that risk is reduced if a demand node/link can be responded to by an emergency team. Moreover, the reduced risk also depends on the type of material transported and the travel distance from the accident site to the response team location. It is reasonable to assume that the response teams always drive to the accident sites along a shortest path, and let the midpoint be the concentrated point of a road link. Then the travel distance from a response node q to a demand link (i,j) can be described using the shortest distance from the midpoint of link (i,j) to node q, that is, Dijq. Meanwhile, this distance cannot exceed the maximum service distance Dmax of the emergency response teams. Therefore, if link (i,j) is covered by node q, that is, hijq=1, then the reduced risk can be described as a function of the distance Dijq and Dmax. For convenience, it is assumed that the reduced risk and the travel distance meet a linear relation such as α(Dmax-Dijq)hijqrijk.Here,α(Dmax-Dijq)hijq describes the percentage of reduced risk when link (i,j) is covered by emergency response team located in node q. α is a coefficient related to the category of undesirable accident consequences and the power of the emergency response teams. αDmax reflects the maximum power of the emergency response teams to service a hazardous accident, that is, how much the accident consequences can be reduced when an accident is serviced by a very timely emergency team. The maximum service distance Dmax can be obtained from the experience of the emergency response teams. The parameter α can be estimated by an analysis of the category of undesirable accident consequences and the power of the emergency response teams as outlined in Figure 6. First, the service ability of the emergency teams should be evaluated. Then, the undesirable consequences are classified and their weights are determined based on the analysis of past accidents. Next, for each category of accident consequences, estimate a possible range of the decrease if an accident is serviced by a very timely emergency response team. The median values of these ranges are taken as the most possible values. Finally, the weight sum of these values is determined as the final value of parameter α. In fact, it is very difficult to get an accurate value of α because it is also dependent on the managers’ attitude. However, a sensitivity analysis on α can assist in the managers’ decision making. Thus, the total uncertain transportation risk can be described as
(18)∑k=1K‍∑(i,j)∈A‍∑q∈H‍dkr~~ijklijxijk(1-α(Dmax-Dijq)hijq).Figure 6
Estimation of the value ofα.In the risk described previously,r~~ijk is a Fu-Fu variable, so the total risk is also a Fu-Fu variable. However, it is difficult to make a decision when it involves uncertain information, so it is necessary to transform the Fu-Fu risk to a determinate one. In this case, the authorities tend to design a network with minimal expected risk. That is, the Fu-Fu risk can be transformed to a determinate one by an expected value operation; the expected total transportation risk can be described as follows:
(19)Risk=E[∑k=1K‍∑(i,j)∈A‍∑q∈H‍dkr~~ijklijxijkRisk=E[∑k=1K‍∑(i,j)∈A‍∑q∈H‍×(1-α(Dmax-Dijq)hijq)∑k=1K‍∑(i,j)∈A‍∑q∈H‍].In the expected total risk, onlyr~~ijk is a Fu-Fu variable. Based on Definition 4 and Theorem 5, the expected total risk can be transformed into following equation:
(20)Risk=∑k=1K‍∑(i,j)∈A‍∑q∈H‍dkE[r~~ijk]lijxijk×(1-α(Dmax-Dijq)hijq).For this objective, the variablesxijk and hij are solved in the lower level model. Hence, the network design model can be modelled as follows:
(21)minyijRisk=∑k=1K‍∑(i,j)∈A‍∑q∈H‍dkE[r~~ijk]lijxijkminyijRisk=∑k=1K‍∑(i,j)∈A‍∑q∈H‍×(1-α(Dmax-Dijq)hijq)s.t.{yij=yji,yij∈{0,1}wherexijkandhijaresolvedinthelower-levelmodel.


## 3.3. A Network Flow Model

The main concern of a government authority is to minimize the total transportation risk by designing a new transportation network. However, risk is also affected by the carriers’ transportation routes and the location of the emergency response teams. The hazmat carriers choose their routes from the upper-level designed network. Thus, this is a multicommodity network flow assignment problem. To these carriers, the minimization of their total transportation cost is a primary objective which conflicts with that of the government authority as cost is not their primary concern. Hence, the flow model should be considered as a separate model, rather than being merged into the network design model. Therefore, the authority faces a bilevel decision problem. That is, while designing transportation network, the authority must also consider the actual use of the hazmat network by the carriers and emergency response teams. Erkut and Gzara [8] also carried out a comparison analysis to explain why it is better to consider the hazmat network design problem as a bilevel model.In the model, the total cost can be described as∑k=1K∑(i,j)∈Adkcijklijxijk. For the flow problem, two constraints must be met. One is the flow equilibrium constraint which ensures the flow of commodity k from its origin to its destination:
(22)∑(i,j)∈A‍xijk-∑(i,j)∈A‍xjik={-1,j=s(k)1,j=t(k),∀j∈V,k=1,2,…,K0,otherwise.The others are the road capacity constraint, which ensures that only the routes selected by the government can be used by the carriers, and the logic constraint, which ensures that the variables only take a binary value:(23)xijk≤yij,xijk∈{0,1},∀(i,j)∈A,k=1,2,…,K.The objective and constraints constitute the routes choice model of the carriers as follows:(24)minxijkCost=∑k=1K‍∑(i,j)∈A‍dkcijklijxijks.t.{∑(i,j)∈A‍xijk-∑(i,j)∈A‍xjik={-1,j=s(k)1,j=t(k),∀j∈V,k=1,2,…,K0,otherwisexijk≤yij,xijk∈{0,1},∀(i,j)∈A,k=1,2,…,K.


## 3.4. A Maximal Arc-Covering Location Model

Usually, a transportation network is designed by a local traffic and transportation department. However, the location of emergency response teams is often chosen by the emergency response departments such as the fire department, the first-aid department, the police office, or their union. The location of these emergency response teams plays an important role in mitigating transportation risk, so the location is considered to be a decision on the lower level when designing the transportation network. The emergency response departments decide the teams location based on the road network design and the carriers’ chosen routes. They hope to cover all the road links, so maximizing the total weighted arc length covered is the objective. This is a maximal covering location problem for hazardous materials transportation. The maximal covering location model was originally developed by Church and ReVelle [33] and was used to locate hazmat response teams in [10, 32].Here, the weight of link(i,j) describes the gain when link (i,j) is covered by emergency response teams, that is, the reduced risk ∑q∈H∑k=1Kdklijα(Dmax-Dijq)r~~ijkxijk. The expected value operation is used to deal with the Fu-Fu risk r~~ijk; where hijq takes value 1 if (i,j) is covered by emergency response team located at note q and 0 otherwise. Then the objective can be modelled as follows:
(25)Cover=∑(i,j)∈A‍∑q∈H‍∑k=1K‍dklijα(Dmax-Dijq)E[r~~ijk]xijkhijq.At the same time, two constraints must also be considered. First, all located teams must be equal to its sum. Letzq describe whether a team is located in node q, and fmax describes the total number of located emergency response teams; then the constraint can be stated as follows:
(26)∑q∈H‍zq=fmax.Moreover, a link(i,j) is covered by node q only if a team is located at node q; then this constraint can be stated as follows:
(27)hijq≤zq.Meanwhile, in most cases, when an accident occurs, only one nearest emergency response team is arranged to service. Therefore, it is assumed that each link is only serviced by a single response team:(28)∑q∈H‍hijq≤1.Finally, the maximum travel distance of an emergency response team to any link traversed by hazardous materials within its jurisdiction should not exceed the thresholdDmax. It can be described as the following constraint:
(29)∑q∈H‍Dijqhijq≤Dmax.From the aforementioned, the maximal arc-covering location problem can be modelled as(30)maxzq,hijqCover=∑(i,j)∈A‍∑q∈H‍∑k=1K‍dklijαmaxzq,hijqCover=∑(i,j)∈A‍∑q∈H‍∑k=1K‍y×(Dmax-Dijq)E[r~~ijk]xijkhijqs.t.{∑q∈H‍zq=fmaxhijq≤zq,∀(i,j)∈A,q∈H∑q∈H‍hijq≤1,∀(i,j)∈A∑q∈H‍Dijqhijq≤Dmax,∀(i,j)∈Azq,hijq∈{0,1}.


## 3.5. Bilevel Programming Model

As outlined in Section2, the considered problem should be modelled as a bilevel programming model. The decision maker on the upper level is the authority who hopes to mitigate the hazmats transportation risk by designing a new road network. They know that the risk is dependent on the carriers’ selected routes and the designed emergency response network which they do not have the right to determine. Fortunately, they also know that the carriers and the emergency response department will make decisions based on the designed network. Therefore, the authority can consider their decision making and influence the lower-level model. On the lower level, the carriers first choose their routes so that total transportation cost is minimized. Then, the emergency response department locates their emergency service teams with the objective of maximizing the total weighted arc length covered. Hence, the complete bilevel programming model can be established based on the previous discussion as follows:
(31)minyijRisk=∑k=1K‍∑(i,j)∈A‍∑q∈H‍dkE[r~~ijk]lijxijkminyijRisk=∑k=1K‍∑(i,j)∈A‍∑q∈H‍×(1-α(Dmax-Dijq)hijq)s.t.{yij=yji,yij∈{0,1}maxzq,hijqCover=∑(i,j)∈A∑q∈H∑k=1Kdklijαmaxzq,hijqCover=∑∑∑∑∑×(Dmax-Dijq)E[r~~ijk]maxzq,hijqCover=∑∑∑∑∑×xijkhijqs.t.{∑q∈H‍zq=fmaxhijq≤zq,∀(i,j)∈A,q∈H∑q∈H‍hijq≤1,∀(i,j)∈A∑q∈H‍Dijqhijq≤Dmax,∀(i,j)∈Azq,hijq∈{0,1}minxijkCost=∑k=1K‍∑(i,j)∈A‍dkcijklijxijks.t.{∑(i,j)∈A‍xijk-∑(i,j)∈A‍xjik=-1,∑(i,j)∈A‍xijk-xjikj=s(k)∑(i,j)∈A‍xijk-∑(i,j)∈A‍xjik=1,∑(i,j)∈A‍xijk-∑(i,j)∈A‍j=t(k)∑(i,j)∈A‍xijk-∑(i,j)∈A‍xjik=0,∑(i,j)∈A‍xijk-j≠s(k),t(k)xijk≤yijxijk∈{0,1}.In this model, the Fu-Fu risks have been dealt with using expected value operation defined in (7). For some especial types of Fu-Fu variables, the expected risk can be calculated directly. As for others, they can be calculated by some simulation methods as introduced in [22]. In this paper, the risk [r~~ijk] is modelled as a discrete triangular Fu-Fu variable (see (9)), and its expected value can be calculated by (10).


## 4. An Improved ABC Algorithm with Priority-Based Encoding

The proposed model is a bilevel programming model which is considered an NP-hard [34] problem and is a strong NP-hard problem [35]. It is often difficult to obtain an analytical optimal solution for such problems and the most commonly used methods are to obtain a numerically optimal solution or a numerically efficient solution using an approximation algorithm. In addition, transportation networks contain numerous nodes and links, which greatly increase the computing complexity. In such conditions, many heuristic algorithms have been developed to solve the bilevel programming problems. In this paper, an improved artificial bee colony algorithm with priority-based encoding is proposed.


### 4.1. Introduction to ABC Algorithm

The artificial bee colony (ABC) algorithm, proposed by Karaboga in 2005 for real-parameter optimization, is a recently introduced optimization algorithm which simulates the foraging behaviour of a bee colony [36]. In the ABC algorithm, the position of a food source represents a possible solution to the optimization problem and the nectar amount of a food source corresponds to the quality (fitness) of the associated solution. The colony of artificial bees is made up of three groups of bees: employed bees, onlookers, and scouts. Half of the colony consists of employed artificial bees and the others are composed of the onlookers. Each food source has one employed bee. After an initial population (food source positions), each employed bee starts to exploit the discovered source and then returns to the hive with the nectar to onlooker bees. Onlooker bees wait in the hive and decide on a food source to exploit based on the information shared by the employed bees. The exploitation of employed bees and onlooker bees represents the modification of solutions. In order to find better solution, if a source is exhausted (i.e., the position of the food source has not been modified through a predetermined number of cycles), the employed bee will become a scout to search for a new source. Like this, the population is subjected to a repeat of the cycles of the search processes of the employed, onlooker, and scout bees, respectively, until the termination criteria are met.The ABC algorithm has been successfully applied to such areas as scheduling, clustering, and engineering design. Results have showed that the performance of the ABC algorithm is better than or similar to other optimization algorithms although it uses less control parameters and it can be efficiently used for solving multivariable and multidimensional optimization problems [37–39]. The problem considered in this paper is a multivariable problem, so the ABC algorithm was chosen. To express the problem more effectively and solve the model more quickly, a priority-based encoding method is also proposed.


### 4.2. Overall Procedure of the Proposed Algorithm

For the proposed algorithm, the ABC algorithm is primarily used to solve the upper-level programming. First food source positions are randomly generated and then encoded into the upper-level road networks. Next, after solving the lower-level programming using existing methods, these sources are evaluated. Then the neighbourhood is searched and the sources are updated by employed bees, onlookers, and scouts, respectively. The flow chart of the proposed algorithm is illustrated in Figure7 and the main procedure is presented as follows.Figure 7
The flow chart of the proposed algorithm.Main procedure for the proposed ABC algorithm is as follows.Step 1.
Initialize the number of food sourcesSN, the maximum cycle number MCN, and the control parameter for the abandoned food source limit. Let iteration t=1. For i=1,…,SN, generate the position of food source Fi and let the holding trials counter coui=0.Step 2.
Encode each food source positionFi into the solution using priority-based encoding method: the set of upper-level decision variables {yij}, that is, the designed transportation.Step 3.
Evaluate each food source siteFi by solving the lower-level flow programming and teams location problem with the encoded initialization result and calculating the fitness value of each food source. Memorize the position of the best food source.Step 4.
Update the food sources using employed bees, onlookers, and scouts, respectively.(4.1)
An employed bee produces a modification on the food source positionFi depending on local information and finds a neighbouring food source and then evaluates its quality using the procedures in Step 2 and Step 3. If it represents better than Fi, the employed bee memorizes it as the new position Fi.(4.2)
An onlooker bee evaluates the fitness value informationfiti of food source i taken from all employed bees and modifies the food source site with a probability pi related to its fitness value. If the position of the modified food source is better than Fi, replace it.(4.3)
Determine an abandoned solution for the scout. If the position of food sourcei is not modified after the employed bees and onlooker bees phases, its trials counter coui is increased by 1; otherwise, the counter is reset to 0. Check the counter coui: if coui>limits, replace Fi with a new randomly produced position and reset coui to 0.Step 5.
If the stopping criterion is met, that is,t>MCN, stop. Otherwise, let t=t+1 and go to Step 2.


### 4.3. Food Source Representation and Encoding Method

The problem here is composed of three nested decision problems. The outer decision problem is the decision of the government (leader). The leader decides a feasible transportation network so that each commodity can be transported to the destination node from its source node. Therefore, the feasible solution can be given by the set of feasible paths, and each path is associated with a transportation commodity. In these heuristic algorithms for the path/network problem, encoding methods are the focus. Priority-based encoding method is one of the most popular methods which has been used with PSO [40] and GA [41]. Here, it is proposed to encode the paths into a food source.


#### 4.3.1. Matrix Representation of Food Source

In the problem, a feasible network is made up of several transportation paths. Each path can be first stated as a priority sequence of 1 toN. Hence, the position of each food source Fi can be represented as a matrix with K rows and N columns. That is,Fi=[Fi11, Fi12,…, Fi1N; Fi21, Fi22,…, Fi2N;…; FiK1, FiK2,…, FiKN], where each row Fik=[Fik1,Fik2,…,FikN] is a priority sequence of 1 to N. At the beginning of the algorithm, all food sources are generated randomly. Figure 8 illustrates an example of food source representation in a network with 10 nodes.Figure 8
An example of food source representation in a network with 10 nodes.


#### 4.3.2. Encoding Food Source into Transportation Network

The food sources can be encoded into a feasible network using a priority-based encoding method. For the path construction of each commodity, the source nodes(k) is taken as the first node. Then there are usually several nodes available (connected with the given node) at each step; the one with priority is chosen into the path, and the priority of the chosen node is reset to 0. This process is iterated until the destination node t(k) is found. Take the network in Figure 8 as an example and assume that there are three considered commodities: 1–10, 2–9, and 3–8. Then according to the food source codes (k=1,2,3), three paths can be easily determined using the given encoding method: 1→4→5→6→9→7→10,  2→4→7→9,  3→5→7→6→8. Further, the feasible solution can be found: y14=y45=y56=y69=y97=y7,10=y24=y47=y79=y35=y57=y76=y68=1. Let s(k), t(k) be the original node and destination node of commodity k, let Pk record the encoded path of commodity k, let η record the chosen node, and let Xk(j) record the priority of node j for commodity k; then the detailed encoding procedure is as follows in Algorithm 1.Algorithm 1:Procedure: priority-based encoding.
Begin: put in food sourceFi
fork=1 to K
letXk=Fik, n=0, Pk←{s(k)}, η=s(k), Xk(η)=0
while (n<N) do
if (η=t(k))
break;
else
J={j∣j∈V and (η,j)∈A}, γ=argmaxj∈J{Xij}, Pk←{Pk,γ};
n=n+1, Xk(j)=0, η=γ;
end if
end while
end for
End: put out the complete pathPk


### 4.4. Evaluation of Food Sources

In this study, the objective function of the upper programming is directly used as the fitness function of each food source. At the first step, each food source is encoded into the solution variables{yij}. Then the following equation is used to calculate the fitness value:
(32)fit=∑k=1K‍∑(i,j)∈A‍dkEMe(SF(r~~ijk))xijk(1-αhij).In the equation,xijk and hij are dependent on the lower-level flow programming and emergency response teams location model. In fact, the upper-level solution has defined a designed network. The Dijkstra algorithm is used to find the shortest path for each transportation commodity; that is, xijk is determined. After this, the “bintprog” function in Matlab is used to solve the emergency response teams location model.


### 4.5. Modification of Food Source Position

In ABC algorithm, the modification of food source position (i.e., search for new solution) consists of three steps. First, an employed bee produces a modification of the position of the food source (solution) in her memory depending on local information. Next, an onlooker bee evaluates the nectar information taken from all employed bees and chooses a food source with a probability related to its nectar amount. Finally, the food source whose nectar is abandoned by the bees is replaced with a new food source by the scout.


#### 4.5.1. Employed Bees Phase

In the ABC algorithm, each employed bee is associated with only one food source site. After being sent to a food source site, each employed bee will produce a modification of the position of the food source depending on local information and find a neighboring food source. Generally, the production of a new food source position is based on a comparison of food source positions. For each food source, a position is randomly chosen and modified according to a comparison with another food source [37].In the original ABC algorithm, each food source was aD-dimensional vector, and the modification is only executed in one dimension of the vector. The considered problem is a discrete optimization problem. Different from the original ABC algorithm, the food sources in the proposed algorithm are represented by a matrix with K rows and N columns (i.e., K random sequences from 1 to N), and each row is a representation of a commodity path. In order to produce a search with a bigger range, each path (i.e., each row of the matrix) should execute a modification; otherwise it is very difficult to find a better food source because of the tiny search space. Moreover, since each row is a sequence from 1 to N, then each position alteration will lead to another corresponding position modification. Therefore, the method described in previous research cannot be applied to this problem. In this study, in order to produce a modification to food source position Fi, a food source position Fj(j≠i) is first randomly chosen. Then for every row k, a position n is also randomly determined. In succession, when the priority of the chosen position is replaced by one of the same positions of the food source j, that is, Fikn=Fjkn, the position whose value is equal to Fikn must be found and its value replaced by Fikn. Finally, the new source position is memorized and the old one forgotten if its fitness value is smaller. An example with 10 nodes is illustrated in Figure 9.Figure 9
An example of neighbourhood search of employed bees.


#### 4.5.2. Onlooker Bees Phase

After all employed bees complete the search process, they share the food source nectar information (solutions) and their position with the onlooker bees on the dance area. Each onlooker bee evaluates the fitness values taken from all employed bees and modifies the food source with a probabilitypi related to its fitness value. The probability value is calculated using the following expression: produces a position modification
(33)pi=fiti∑n=1SNfitn,
where fiti is the fitness value of the solution i evaluated by the employed bee.After that, a parameter value is randomly generated between 0 and 1. If it is smaller thanpi, then the onlooker bee produces a position modification to the food source; otherwise, it chooses the food source from the employed bee. In the employed bees phase, the position modification is only executed to a node of one path. This is a tiny modification and leads to a slow convergence speed. Hence, a greater modification, that is, a path modification, is proposed in this phase. The onlooker first chooses a path Fik and another food source j randomly. Then it will replace the path with the similar path from another food source, that is, Fik=Fjk. An example is given in Figure 10.Figure 10
An example of the neighbourhood search for onlooker bees.


#### 4.5.3. Scout Bees Phase

In each cycle, after all employed bees and onlooker bees complete their searches, the algorithm checks to see if there is any exhausted food source to be abandoned. In order to decide if a source is to be abandoned, a trial countercoui is used. If the position of food source i is not modified after the employed bees and onlooker bees phases, its trials counter coui is increased by 1; otherwise, the counter is reset to 0. Check the counter coui: if coui>limits, replace Fi with a new randomly produced position and reset coui to 0. Here, limits are a predetermined control parameter. Since the scouts can accidentally discover rich, entirely unknown food sources, they plays the role of fast discovery of the group of feasible solutions.


## 4.1. Introduction to ABC Algorithm

The artificial bee colony (ABC) algorithm, proposed by Karaboga in 2005 for real-parameter optimization, is a recently introduced optimization algorithm which simulates the foraging behaviour of a bee colony [36]. In the ABC algorithm, the position of a food source represents a possible solution to the optimization problem and the nectar amount of a food source corresponds to the quality (fitness) of the associated solution. The colony of artificial bees is made up of three groups of bees: employed bees, onlookers, and scouts. Half of the colony consists of employed artificial bees and the others are composed of the onlookers. Each food source has one employed bee. After an initial population (food source positions), each employed bee starts to exploit the discovered source and then returns to the hive with the nectar to onlooker bees. Onlooker bees wait in the hive and decide on a food source to exploit based on the information shared by the employed bees. The exploitation of employed bees and onlooker bees represents the modification of solutions. In order to find better solution, if a source is exhausted (i.e., the position of the food source has not been modified through a predetermined number of cycles), the employed bee will become a scout to search for a new source. Like this, the population is subjected to a repeat of the cycles of the search processes of the employed, onlooker, and scout bees, respectively, until the termination criteria are met.The ABC algorithm has been successfully applied to such areas as scheduling, clustering, and engineering design. Results have showed that the performance of the ABC algorithm is better than or similar to other optimization algorithms although it uses less control parameters and it can be efficiently used for solving multivariable and multidimensional optimization problems [37–39]. The problem considered in this paper is a multivariable problem, so the ABC algorithm was chosen. To express the problem more effectively and solve the model more quickly, a priority-based encoding method is also proposed.


## 4.2. Overall Procedure of the Proposed Algorithm

For the proposed algorithm, the ABC algorithm is primarily used to solve the upper-level programming. First food source positions are randomly generated and then encoded into the upper-level road networks. Next, after solving the lower-level programming using existing methods, these sources are evaluated. Then the neighbourhood is searched and the sources are updated by employed bees, onlookers, and scouts, respectively. The flow chart of the proposed algorithm is illustrated in Figure7 and the main procedure is presented as follows.Figure 7
The flow chart of the proposed algorithm.Main procedure for the proposed ABC algorithm is as follows.Step 1.
Initialize the number of food sourcesSN, the maximum cycle number MCN, and the control parameter for the abandoned food source limit. Let iteration t=1. For i=1,…,SN, generate the position of food source Fi and let the holding trials counter coui=0.Step 2.
Encode each food source positionFi into the solution using priority-based encoding method: the set of upper-level decision variables {yij}, that is, the designed transportation.Step 3.
Evaluate each food source siteFi by solving the lower-level flow programming and teams location problem with the encoded initialization result and calculating the fitness value of each food source. Memorize the position of the best food source.Step 4.
Update the food sources using employed bees, onlookers, and scouts, respectively.(4.1)
An employed bee produces a modification on the food source positionFi depending on local information and finds a neighbouring food source and then evaluates its quality using the procedures in Step 2 and Step 3. If it represents better than Fi, the employed bee memorizes it as the new position Fi.(4.2)
An onlooker bee evaluates the fitness value informationfiti of food source i taken from all employed bees and modifies the food source site with a probability pi related to its fitness value. If the position of the modified food source is better than Fi, replace it.(4.3)
Determine an abandoned solution for the scout. If the position of food sourcei is not modified after the employed bees and onlooker bees phases, its trials counter coui is increased by 1; otherwise, the counter is reset to 0. Check the counter coui: if coui>limits, replace Fi with a new randomly produced position and reset coui to 0.Step 5.
If the stopping criterion is met, that is,t>MCN, stop. Otherwise, let t=t+1 and go to Step 2.


## 4.3. Food Source Representation and Encoding Method

The problem here is composed of three nested decision problems. The outer decision problem is the decision of the government (leader). The leader decides a feasible transportation network so that each commodity can be transported to the destination node from its source node. Therefore, the feasible solution can be given by the set of feasible paths, and each path is associated with a transportation commodity. In these heuristic algorithms for the path/network problem, encoding methods are the focus. Priority-based encoding method is one of the most popular methods which has been used with PSO [40] and GA [41]. Here, it is proposed to encode the paths into a food source.


### 4.3.1. Matrix Representation of Food Source

In the problem, a feasible network is made up of several transportation paths. Each path can be first stated as a priority sequence of 1 toN. Hence, the position of each food source Fi can be represented as a matrix with K rows and N columns. That is,Fi=[Fi11, Fi12,…, Fi1N; Fi21, Fi22,…, Fi2N;…; FiK1, FiK2,…, FiKN], where each row Fik=[Fik1,Fik2,…,FikN] is a priority sequence of 1 to N. At the beginning of the algorithm, all food sources are generated randomly. Figure 8 illustrates an example of food source representation in a network with 10 nodes.Figure 8
An example of food source representation in a network with 10 nodes.


### 4.3.2. Encoding Food Source into Transportation Network

The food sources can be encoded into a feasible network using a priority-based encoding method. For the path construction of each commodity, the source nodes(k) is taken as the first node. Then there are usually several nodes available (connected with the given node) at each step; the one with priority is chosen into the path, and the priority of the chosen node is reset to 0. This process is iterated until the destination node t(k) is found. Take the network in Figure 8 as an example and assume that there are three considered commodities: 1–10, 2–9, and 3–8. Then according to the food source codes (k=1,2,3), three paths can be easily determined using the given encoding method: 1→4→5→6→9→7→10,  2→4→7→9,  3→5→7→6→8. Further, the feasible solution can be found: y14=y45=y56=y69=y97=y7,10=y24=y47=y79=y35=y57=y76=y68=1. Let s(k), t(k) be the original node and destination node of commodity k, let Pk record the encoded path of commodity k, let η record the chosen node, and let Xk(j) record the priority of node j for commodity k; then the detailed encoding procedure is as follows in Algorithm 1.Algorithm 1:Procedure: priority-based encoding.
Begin: put in food sourceFi
fork=1 to K
letXk=Fik, n=0, Pk←{s(k)}, η=s(k), Xk(η)=0
while (n<N) do
if (η=t(k))
break;
else
J={j∣j∈V and (η,j)∈A}, γ=argmaxj∈J{Xij}, Pk←{Pk,γ};
n=n+1, Xk(j)=0, η=γ;
end if
end while
end for
End: put out the complete pathPk


## 4.3.1. Matrix Representation of Food Source

In the problem, a feasible network is made up of several transportation paths. Each path can be first stated as a priority sequence of 1 toN. Hence, the position of each food source Fi can be represented as a matrix with K rows and N columns. That is,Fi=[Fi11, Fi12,…, Fi1N; Fi21, Fi22,…, Fi2N;…; FiK1, FiK2,…, FiKN], where each row Fik=[Fik1,Fik2,…,FikN] is a priority sequence of 1 to N. At the beginning of the algorithm, all food sources are generated randomly. Figure 8 illustrates an example of food source representation in a network with 10 nodes.Figure 8
An example of food source representation in a network with 10 nodes.


## 4.3.2. Encoding Food Source into Transportation Network

The food sources can be encoded into a feasible network using a priority-based encoding method. For the path construction of each commodity, the source nodes(k) is taken as the first node. Then there are usually several nodes available (connected with the given node) at each step; the one with priority is chosen into the path, and the priority of the chosen node is reset to 0. This process is iterated until the destination node t(k) is found. Take the network in Figure 8 as an example and assume that there are three considered commodities: 1–10, 2–9, and 3–8. Then according to the food source codes (k=1,2,3), three paths can be easily determined using the given encoding method: 1→4→5→6→9→7→10,  2→4→7→9,  3→5→7→6→8. Further, the feasible solution can be found: y14=y45=y56=y69=y97=y7,10=y24=y47=y79=y35=y57=y76=y68=1. Let s(k), t(k) be the original node and destination node of commodity k, let Pk record the encoded path of commodity k, let η record the chosen node, and let Xk(j) record the priority of node j for commodity k; then the detailed encoding procedure is as follows in Algorithm 1.Algorithm 1:Procedure: priority-based encoding.
Begin: put in food sourceFi
fork=1 to K
letXk=Fik, n=0, Pk←{s(k)}, η=s(k), Xk(η)=0
while (n<N) do
if (η=t(k))
break;
else
J={j∣j∈V and (η,j)∈A}, γ=argmaxj∈J{Xij}, Pk←{Pk,γ};
n=n+1, Xk(j)=0, η=γ;
end if
end while
end for
End: put out the complete pathPk


## 4.4. Evaluation of Food Sources

In this study, the objective function of the upper programming is directly used as the fitness function of each food source. At the first step, each food source is encoded into the solution variables{yij}. Then the following equation is used to calculate the fitness value:
(32)fit=∑k=1K‍∑(i,j)∈A‍dkEMe(SF(r~~ijk))xijk(1-αhij).In the equation,xijk and hij are dependent on the lower-level flow programming and emergency response teams location model. In fact, the upper-level solution has defined a designed network. The Dijkstra algorithm is used to find the shortest path for each transportation commodity; that is, xijk is determined. After this, the “bintprog” function in Matlab is used to solve the emergency response teams location model.


## 4.5. Modification of Food Source Position

In ABC algorithm, the modification of food source position (i.e., search for new solution) consists of three steps. First, an employed bee produces a modification of the position of the food source (solution) in her memory depending on local information. Next, an onlooker bee evaluates the nectar information taken from all employed bees and chooses a food source with a probability related to its nectar amount. Finally, the food source whose nectar is abandoned by the bees is replaced with a new food source by the scout.


### 4.5.1. Employed Bees Phase

In the ABC algorithm, each employed bee is associated with only one food source site. After being sent to a food source site, each employed bee will produce a modification of the position of the food source depending on local information and find a neighboring food source. Generally, the production of a new food source position is based on a comparison of food source positions. For each food source, a position is randomly chosen and modified according to a comparison with another food source [37].In the original ABC algorithm, each food source was aD-dimensional vector, and the modification is only executed in one dimension of the vector. The considered problem is a discrete optimization problem. Different from the original ABC algorithm, the food sources in the proposed algorithm are represented by a matrix with K rows and N columns (i.e., K random sequences from 1 to N), and each row is a representation of a commodity path. In order to produce a search with a bigger range, each path (i.e., each row of the matrix) should execute a modification; otherwise it is very difficult to find a better food source because of the tiny search space. Moreover, since each row is a sequence from 1 to N, then each position alteration will lead to another corresponding position modification. Therefore, the method described in previous research cannot be applied to this problem. In this study, in order to produce a modification to food source position Fi, a food source position Fj(j≠i) is first randomly chosen. Then for every row k, a position n is also randomly determined. In succession, when the priority of the chosen position is replaced by one of the same positions of the food source j, that is, Fikn=Fjkn, the position whose value is equal to Fikn must be found and its value replaced by Fikn. Finally, the new source position is memorized and the old one forgotten if its fitness value is smaller. An example with 10 nodes is illustrated in Figure 9.Figure 9
An example of neighbourhood search of employed bees.


### 4.5.2. Onlooker Bees Phase

After all employed bees complete the search process, they share the food source nectar information (solutions) and their position with the onlooker bees on the dance area. Each onlooker bee evaluates the fitness values taken from all employed bees and modifies the food source with a probabilitypi related to its fitness value. The probability value is calculated using the following expression: produces a position modification
(33)pi=fiti∑n=1SNfitn,
where fiti is the fitness value of the solution i evaluated by the employed bee.After that, a parameter value is randomly generated between 0 and 1. If it is smaller thanpi, then the onlooker bee produces a position modification to the food source; otherwise, it chooses the food source from the employed bee. In the employed bees phase, the position modification is only executed to a node of one path. This is a tiny modification and leads to a slow convergence speed. Hence, a greater modification, that is, a path modification, is proposed in this phase. The onlooker first chooses a path Fik and another food source j randomly. Then it will replace the path with the similar path from another food source, that is, Fik=Fjk. An example is given in Figure 10.Figure 10
An example of the neighbourhood search for onlooker bees.


### 4.5.3. Scout Bees Phase

In each cycle, after all employed bees and onlooker bees complete their searches, the algorithm checks to see if there is any exhausted food source to be abandoned. In order to decide if a source is to be abandoned, a trial countercoui is used. If the position of food source i is not modified after the employed bees and onlooker bees phases, its trials counter coui is increased by 1; otherwise, the counter is reset to 0. Check the counter coui: if coui>limits, replace Fi with a new randomly produced position and reset coui to 0. Here, limits are a predetermined control parameter. Since the scouts can accidentally discover rich, entirely unknown food sources, they plays the role of fast discovery of the group of feasible solutions.


## 4.5.1. Employed Bees Phase

In the ABC algorithm, each employed bee is associated with only one food source site. After being sent to a food source site, each employed bee will produce a modification of the position of the food source depending on local information and find a neighboring food source. Generally, the production of a new food source position is based on a comparison of food source positions. For each food source, a position is randomly chosen and modified according to a comparison with another food source [37].In the original ABC algorithm, each food source was aD-dimensional vector, and the modification is only executed in one dimension of the vector. The considered problem is a discrete optimization problem. Different from the original ABC algorithm, the food sources in the proposed algorithm are represented by a matrix with K rows and N columns (i.e., K random sequences from 1 to N), and each row is a representation of a commodity path. In order to produce a search with a bigger range, each path (i.e., each row of the matrix) should execute a modification; otherwise it is very difficult to find a better food source because of the tiny search space. Moreover, since each row is a sequence from 1 to N, then each position alteration will lead to another corresponding position modification. Therefore, the method described in previous research cannot be applied to this problem. In this study, in order to produce a modification to food source position Fi, a food source position Fj(j≠i) is first randomly chosen. Then for every row k, a position n is also randomly determined. In succession, when the priority of the chosen position is replaced by one of the same positions of the food source j, that is, Fikn=Fjkn, the position whose value is equal to Fikn must be found and its value replaced by Fikn. Finally, the new source position is memorized and the old one forgotten if its fitness value is smaller. An example with 10 nodes is illustrated in Figure 9.Figure 9
An example of neighbourhood search of employed bees.


## 4.5.2. Onlooker Bees Phase

After all employed bees complete the search process, they share the food source nectar information (solutions) and their position with the onlooker bees on the dance area. Each onlooker bee evaluates the fitness values taken from all employed bees and modifies the food source with a probabilitypi related to its fitness value. The probability value is calculated using the following expression: produces a position modification
(33)pi=fiti∑n=1SNfitn,
where fiti is the fitness value of the solution i evaluated by the employed bee.After that, a parameter value is randomly generated between 0 and 1. If it is smaller thanpi, then the onlooker bee produces a position modification to the food source; otherwise, it chooses the food source from the employed bee. In the employed bees phase, the position modification is only executed to a node of one path. This is a tiny modification and leads to a slow convergence speed. Hence, a greater modification, that is, a path modification, is proposed in this phase. The onlooker first chooses a path Fik and another food source j randomly. Then it will replace the path with the similar path from another food source, that is, Fik=Fjk. An example is given in Figure 10.Figure 10
An example of the neighbourhood search for onlooker bees.


## 4.5.3. Scout Bees Phase

In each cycle, after all employed bees and onlooker bees complete their searches, the algorithm checks to see if there is any exhausted food source to be abandoned. In order to decide if a source is to be abandoned, a trial countercoui is used. If the position of food source i is not modified after the employed bees and onlooker bees phases, its trials counter coui is increased by 1; otherwise, the counter is reset to 0. Check the counter coui: if coui>limits, replace Fi with a new randomly produced position and reset coui to 0. Here, limits are a predetermined control parameter. Since the scouts can accidentally discover rich, entirely unknown food sources, they plays the role of fast discovery of the group of feasible solutions.


## 5. Case Study: HTND for Shuibuya Project

In this section, computational experiments that were carried out on a real application are presented. Through an illustrative example on the data set adopted from a case study, the proposed method is validated and the efficiency of the algorithm is tested. The data for material requisition, transportation cost, transportation risk, road network, and others involved in the case are from the Shuibuya Hydropower Station large-scale construction project. The case is introduced to demonstrate the potential real world applications of the proposed methods.


### 5.1. Presentation of Case Problem

The Shuibuya hydropower project is located in Badong County in the middle reaches of the Qingjiang River. It is the first cascaded project on the Qingjiang mainstream. At 233 m tall and containing 15,640,000 m3 of material, it is the tallest concrete face rockfill dam in the world. This is a large scale project which has a complex transportation network for construction materials within the construction yard. In the network, large quantities of hazardous materials such as explosives are also shipped every day. However, there is no specified network for hazmat transportation. As a result, the project managers and public in the yard are concerned about the potential risk. Therefore, it is necessary to design a hazmat transportation network in the construction yard. In this case, the project manager is the decision maker on the upper level who hopes to mitigate the hazmat transportation risk by designing a new road network. Meanwhile, he also considers the decision making of the carriers and the emergency response department. Here, the emergency response department is a union of the fire-fighting department, the first-aid department, and the police station.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 has 17 preexisting trunk roads and 9 temporary roads located on the left and right banks forming a solid network of cycle traffic. There is a Cross River bridge connecting the left and right banks. The actual construction floor plan is as in Figure11. In order to apply the proposed methods more conveniently, adjacent roads of the same type have been combined and the road shapes have been ignored. An abstracted transportation network illustration is in Figure 12. The illustration has 37 nodes and 53 links in total and the application of the proposed methodologies does not affect the result of the case problem.Figure 11
The floor plan of the Shuibuya project.Figure 12
An abstracted network of the Shuibuya project.In the network, there are two hazardous materials warehouses (node no. 2 and no. 23), eight main demand nodes (node no. 1, no. 7, no. 9, no. 16, no. 18, no. 25, no. 32, and no. 36), and eight transportation commodities. Detailed data on the transportation commodities is shown in Table2.Table 2
Detailed data on the eight transportation commodities.
Commodityk
1
2
3
4
5
6
7
8
Sources(k)
No. 2
No. 2
No. 2
No. 2
No. 2
No. 23
No. 23
No. 23
Destinationt(k)
No. 1
No. 7
No. 9
No. 16
No. 18
No. 25
No. 32
No. 36
d
k (kg)
600
420
1500
2000
800
420
2000
480


### 5.2. Data Collection and Computing Results

In this case, all data on the transportation network, hazardous material (explosive), and emergency response were obtained from the Hubei Qingjiang Shuibuya Project Construction Company. The transportation data is shown in Table3. However, there was no ready-made data on transportation risk. In this paper, transportation risk is described as a Fu-Fu variable composed of two fuzzy factors: fuzzy accident possibility and fuzzy accident consequences. In a construction project, an accident may lead to various losses such as fatalities, damage to road networks and construction facilities, project duration delay, and social impacts. In this paper, fatalities, project duration delay, and damage to road network and construction facilities are considered accident consequences. The data fuzzification method outlined in Section 2 is used to obtain the accident probability and consequences of fuzzy data. Based on this method, and considering the convenience of use and interpretation, the probabilities are modelled as discrete fuzzy numbers and the consequences are transformed into costs and modelled as triangle fuzzy numbers. Hence, the unit transportation risk for each road section is modelled as a discrete triangular Fu-Fu number, and the detailed data is shown in Table 4.Table 3
Data on the road network of Shuibuya hydropower project.
Links
Nodes
Distant(miles)
Cost(CNY)
1
#1, #4
0.95
5.2
2
#2, #3
0.76
3.8
3
#2, #4
0.86
4.8
4
#3, #6
1.98
9.5
5
#4, #5
0.98
5.0
6
#4, #7
0.75
3.0
7
#5, #6
0.27
1.2
8
#6, #9
0.87
4.0
9
#7, #8
1.02
5.0
10
#8, #9
0.95
5.0
11
#9, #10
0.53
2.5
12
#10, #16
0.82
3.6
13
#7, #11
1.78
9.6
14
#8, #13
0.85
3.5
15
#9, #13
0.82
4.6
16
#11, #12
0.22
1.0
17
#12, #13
0.67
3.5
18
#13, #14
0.67
3.5
19
#10, #14
0.87
4.6
20
#14, #15
1.33
6.2
21
#15, #16
0.33
2.1
22
#11, #17
0.22
2.0
23
#17, #18
1.90
10.8
24
#12, #18
0.33
2.1
25
#14, #19
2.92
10.3
26
#18, #19
1.58
6.8
27
#17, #20
5.24
29.7
28
#18, #20
4.11
26.2
29
#19, #33
0.86
2.5
30
#21, #22
0.90
4.0
31
#22, #23
0.88
2.5
32
#21, #24
1.69
7.5
33
#21, #26
0.55
2.5
34
#22, #26
0.60
3.0
35
#25, #26
0.55
2.9
36
#24, #25
0.46
3.6
37
#25, #29
0.27
1.9
38
#29, #30
0.68
3.2
39
#24, #28
0.22
1.9
40
#28, #29
0.53
3.8
41
#28, #30
0.68
3.2
42
#30, #31
1.05
5.4
43
#31, #32
0.27
0.9
44
#23, #27
0.70
3.9
45
#27, #32
2.00
12.9
46
#25, #27
0.76
3.8
47
#32, #33
1.32
7.6
48
#32, #34
2.00
14.8
49
#34, #36
1.56
8.6
50
#33, #35
0.36
2.1
51
#35, #36
2.22
14.6
52
#34, #35
0.81
5.7
53
#36, #37
0.86
5.8Table 4
Transportation risk on each road section.
Links
Transportation risk (CNY/kg)
Links
Transportation risk (CNY/kg)
Links
Transportation risk
(
6,12,18
)
,
0.2
(
6,9
,
12
)
,
0.3
(
12,15,18
)
,
0.3
(
8,16,24
)
,
0.7
(
8,12,16
)
,
0.7
(
16,20,24
)
,
0.7
1
(
10,20,30
)
,
1.0
2, 3
(
10,15,20
)
,
1.0
4
(
20,25,30
)
,
1.0
(
12,24,36
)
,
0.6
(
12,18,24
)
,
0.6
(
24,30,36
)
,
0.6
(
14,28,42
),0.3
(
14,21,28
)
,
0.4
(
28,35,42
)
,
0.4
(
24,32,40
)
,
0.2
(
48,60,72
)
,
0.3
(
8,12,16
)
,
0.3
(
30,40,50
)
,
0.7
(
56,70,84
)
,
0.7
(
10,15,20
)
,
0.7
5, 6
(
36,48,60
)
,
1.0
7, 8
(
64,80,96
)
,
1.0
9
(
12,18,24
)
,
1.0
(
42,56,70
)
,
0.6
(
72,90,108
)
,
0.6
(
14,21,28
)
,
0.6
(
48,64,80
)
,
0.3
(
80,100,120
)
,
0.4
(
16,24,32
)
,
0.4
(
8,12,16
)
,
0.3
(
80,90,100
)
,
0.3
(
8,10,12
)
,
0.3
(
10,15,20
)
,
0.7
(
96,108,120
)
,
0.7
(
12,15,18
)
,
0.7
10
(
12,18,24
)
,
1.0
11, 12
(
112,126,140
)
,
1.0
13
(
16,20,24
)
,
1.0
(
14,21,28
)
,
0.6
(
128,144,160
)
,
0.6
(
20,25,30
)
,
0.6
(
16,24,32
)
,
0.4
(
144,162,180
)
,
0.4
(
24,30,36
)
,
0.4
(
8,10,12
)
,
0.3
(
8,10,12
)
,
0.3
(
12,15,18
)
,
0.3
(
12,15,18
)
,
0.7
(
12,15,18
)
,
0.7
(
16,20,24
)
,
0.7
14
(
16,20,24
)
,
1.0
15
(
16,20,24
)
,
1.0
16, 17
(
20,25,30
)
,
1.0
(
20,25,30
)
,
0.6
(
20,25,30
)
,
0.6
(
24,30,36
)
,
0.6
(
24,30,36
)
,
0.4
(
24,30,36
)
,
0.4
(
28,35,42
)
,
0.4
(
12,15,18
)
,
0.3
(
24,32,40
)
,
0.2
(
12,15,18
)
,
0.3
(
16,20,24
)
,
0.7
(
30,40,50
)
,
0.7
(
16,20,24
)
,
0.7
18, 19
(
20,25,30
)
,
1.0
20, 21
(
36,48,60
)
,
1.0
22, 23
(
20,25,30
)
,
1.0
(
24,30,36
)
,
0.6
(
42,56,70
)
,
0.6
(
24,30,36
)
,
0.6
(
28,35,42
)
,
0.4
(
48,64,80
)
,
0.3
(
28,35,42
)
,
0.4
(
12,15,18
)
,
0.3
(
18,21,24
)
,
0.3
(
12,15,18
)
,
0.3
(
16,20,24
)
,
0.7
(
24,28,32
)
,
0.7
(
16,20,24
)
,
0.7
24, 25
(
20,25,30
)
,
1.0
26
(
30,35,40
)
,
1.0
27, 28
(
20,25,30
)
,
1.0
(
24,30,36
)
,
0.6
(
36,42,48
)
,
0.6
(
24,30,36
)
,
0.6
(
28,35,42
)
,
0.4
(
42,49,56
)
,
0.4
(
28,35,42
)
,
0.4
(
6,12,18
)
,
0.2
(
12,15,18
)
,
0.3
(
6,9
,
12
)
,
0.3
(
8,16,24
)
,
0.7
(
16,20,24
)
,
0.7
(
8,12,16
)
,
0.7
29
(
10,20,30
)
,
1.0
30, 33
(
20,25,30
)
,
1.0
31
(
10,15,20
)
,
1.0
(
12,24,36
)
,
0.6
(
24,30,36
)
,
0.6
(
12,18,24
)
,
0.6
(
14,28,42
)
,
0.3
(
28,35,42
)
,
0.4
(
14,21,28
)
,
0.4
(
6,9
,
12
)
,
0.3
(
6,12,18
)
,
0.2
(
6,12,18
)
,
0.2
(
8,12,16
)
,
0.7
(
8,16,24
)
,
0.7
(
8,16,24
)
,
0.7
32
(
10,15,20
)
,
1.0
34, 35
(
10,20,30
)
,
1.0
36, 40
(
10,20,30
)
,
1.0
(
12,18,24
)
,
0.6
(
12,24,36
)
,
0.6
(
12,24,36
)
,
0.6
(
14,21,28
)
,
0.4
(
14,28,42
)
,
0.3
(
14,28,42
)
,
0.3
(
12,15,18
)
,
0.3
(
6,9
,
12
)
,
0.3
(
6,12,18
)
,
0.2
(
16,20,24
)
,
0.7
(
8,12,16
)
,
0.7
(
8,16,24
)
,
0.7
37, 38
(
20,25,30
)
,
1.0
39, 41
(
10,15,20
)
,
1.0
42, 43
(
10,20,30
)
,
1.0
(
24,30,36
)
,
0.6
(
12,18,24
)
,
0.6
(
12,24,36
)
,
0.6
(
28,35,42
)
,
0.4
(
14,21,28
)
,
0.4
(
14,28,42
)
,
0.3
(
12,15,18
)
,
0.3
(
6,9
,
12
)
,
0.3
(
6,12,18
)
,
0.2
(
16,20,24
)
,
0.7
(
8,12,16
)
,
0.7
(
8,16,24
)
,
0.7
44–46
(
20,25,30
)
,
1.0
47
(
10,15,20
)
,
1.0
48, 49
(
10,20,30
)
,
1.0
(
24,30,36
)
,
0.6
(
12,18,24
)
,
0.6
(
12,24,36
)
,
0.6
(
28,35,42
)
,
0.4
(
14,21,28
)
,
0.4
(
14,28,42
)
,
0.3
(
12,15,18
)
,
0.3
(
6,9
,
12
)
,
0.3
(
16,20,24
)
,
0.7
(
8,12,16
)
,
0.7
50, 51
(
20,25,30
)
,
1.0
52, 53
(
10,15,20
)
,
1.0
(
24,30,36
)
,
0.6
(
12,18,24
)
,
0.6
(
28,35,42
)
,
0.4
(
14,21,28
)
,
0.4In addition, the emergency response department plans to locate two emergency response teams for potential hazmat transportation accidents while six nodes are optional for this location (i.e., no. 7, no. 9, no. 13, no. 23, no. 25, and no. 32). The service range of the emergency response teams is 2 kilometers. Based on this coverage distance, all covered links by each optional location are predetermined. Hence, the setNij is also determined. In order to achieve the maximal coverage, if a link is covered partially, then it is divided into two links (one is entirely covered, and another is not covered) by adding a new node.Using this data, after running the computer program for the improved ABC algorithm using MATLAB 2007, the solutions for the case problem and the efficiency of the proposed algorithm were obtained.The algorithm and model parameters for the case problem were set as follows: the number of the food sourcesSN=40, the value of limit=20, the maximum cycle number MCN=200, and the mitigated risk coefficient for emergency response α=0.3. The computer running environment was an Intel Core 2 Duo 2.26 GHz clock pulse with 2048 MB memory. After 10.46 minutes on average, the optimal solutions for the bilevel programming were determined. The optimal transportation network is shown in Figure 13.Figure 13
The optimal network for the explosives transportation.


### 5.3. Model Analysis

#### 5.3.1. Two Different Networks

As discussed, the authority needs to consider the emergency response network when the hazardous materials transportation network is designed, because an emergency response network has an impact on the transportation network.In this study, a comparison is given for a network considering emergency response and one without regard to emergency response. Figure14 illustrates the two different networks. From this, it can be seen that there are some distinct differences between the two networks. Links No. 7, No. 8, No. 15, No. 16, and No. 17 would be selected if emergency response were considered though their unit risks are higher than other links, while they would not be included in the network without considering the emergency response. In addition, there is also a lower risk (806.44 thousand CNY and 1049.06 thousand CNY in the two networks, resp.) if the emergency response is considered when designing the transportation network. Therefore, the total transportation risk is reduced by 23.13% due to the existence of emergency response teams.Figure 14
Difference between the two proposed networks.


#### 5.3.2. Sensitivity Analysis toα

In Section3, it is assumed that the risk of link (i,j) can be reduced by α(Dmax-Dijq)hijqrijk if it is covered by an emergency response team located at node q. It is known that the value of parameter α is a key factor impacting the results. In fact, the decision maker is able to adjust the parameter to obtain different solutions. A sensitivity analysis is represented in Table 5. From Table 5, it can be seen that total risk decreases as the value of α increases. The designed networks are also related to α. With an adjustment in α, three different optimal networks were found. Each optimal network is related to a value range of α. It shows that the optimal network is not sensitive to α and a range of α is enough for the proposed network designing problem. Hence, the authority only needs to determine the range of α when designing the transportation network.Table 5
Sensitivity analysis of parameterα.
Value ofα
Total risk
Location
Road links for hazmat transportation
α
=
0.00
1049.06
—
{1, 2, 3, 5, 6, 9, 10, 13, 14, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.05
1016.37
#7, #32
{1, 2, 3, 5, 6, 9, 10, 13, 14, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.10
983.67
#7, #32
{1, 2, 3, 5, 6, 9, 10, 13, 14, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.15
949.17
#7, #32
{1, 2, 3, 5, 6, 9, 10, 13, 14, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.20
901.59
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.25
854.02
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.30
806.44
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.35
758.87
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.40
711.29
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.45
663.19
#9, #25
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 37, 38, 42, 43, 48, 49}
α
=
0.50
612.34
#9, #25
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 37, 38, 42, 43, 48, 49}


#### 5.3.3. Uncertainty Analysis

In addition, uncertainty is also an important consideration in this study. Fu-Fu variables are used to model the transportation risk because of the lack of precise data such as the accident probabilities and the accident consequences on specific roads, so the proposed model contains a complex fuzzy factor. Besides the model considering a Fu-Fu factor, there are two other kinds considered, one without uncertainty and one that includes one fuzzy factor. Truly determinate data does not exist in project practice, so the managers make decisions from historical statistical data and their own experience. This factor makes it very difficult to select decision-making data with most of it obtained by ignoring the uncertainty and using an average value or data chosen at random. If only one fuzzy factor is considered such as the fuzzy accident consequences, then other uncertain information needs to be ignored, such as the accident probabilities. Then the accident probabilities have to be randomly chosen from the given interval. All these will lead to an imprecise solution. Table6 shows a comparison result of the three model types. In the comparison, the model considering a complex fuzzy factor finds different solutions by adjusting the optimistic and pessimistic index λ when the others find different solutions by choosing different transportation risk randomly. It is clear that the proposed model considering the Fu-Fu factors has a much better performance than the others, not only in the average value of the results, but also in their stability.Table 6
Comparisons among the three model types based on the risk-minimizing objective.
Type
Complex fuzzy type
Fuzzy type
Determined type
Best result
736.68
720.22
698.36
Worst result
854.35
902.84
932.07
Average result
806.04
834.52
842.36


### 5.4. Algorithm Evaluation

In this paper, an improved artificial bee colony algorithm with priority-based encoding has been proposed. In order to test the efficiency of the algorithm, a comparison with other solution methods was conducted.The proposed problem is modeled as a linear bilevel programming model in this paper. The most common solution strategy for linear bilevel problems is to transform the bilevel model into a single one by the use of its Karush-Kuhn-Tucker (KKT) conditions. However, it is difficult when there are binary variables in the inner models. This also results in that it cannot be solved by common commercial solvers. In theory, it is possible to solve this problem by conducting an enumeration search (ES) over the design variablesyij on the upper level and solving a series of linear programs on the lower level for each selection of design variables. However, this search would be conducted 253 times because of the 57 binary design variables. This is impractical for large-scale problems. Hence, it is not appropriate to solve the problem by KKT conditions and enumeration search in most cases. Fortunately, it is found that the transportation network in this case can be divided into two subnetworks along the river because all these hazmats transportation would not go through the river. Moreover, some road links such as (1,4), (2,3), (17,20), (18,20), and (36,37)  can be determined by observing the network. Therefore, the network can be divided into two networks with 25 and 23 links, respectively. The problem in the case study can be solved by conducting 223+225 searches. In experiments, it spends 186 minutes in running the search program. Although the solution is optimal, the too long time is not acceptable.In order to show the efficiency of the proposed ABC algorithm, a comparative experiment among the proposed ABC algorithm, a particle swarm optimization (PSO), a genetic algorithm (GA), and the enumeration search method was also conducted. In the experiment, the same encoding method (i.e., priority-based encoding) was used for the three algorithms. For the GA, the partially mapped crossover and local search-based mutation method are used. For the PSO, a hybrid particle-updating mechanism is used in the PSO. The other parameters for the algorithms are shown in Table7. The experiments for the algorithms were carried out over 50 times. Figure 15 shows the convergence histories for the three heuristic algorithms. The detailed performance of the four algorithms is stated in Table 8. The results indicate that the ABC needs less time than the GA and has a more stable tendency than the PSO. Moreover, the accuracy of both the ABC and the GA is very high; even the optimal solution can be found with a percentage more than 14%. Therefore the performance of the ABC is on par with that of both the PSO and the GA. The proposed algorithm is also useful and efficient for solving the proposed case.Table 7
Parameter values used in the experiments.
MCN
Pop size
Limit
Weight
c
p
best
>cgbest
r
crossover
r
mutation
ABC
200
40
20
—
—
—
—
—
PSO
200
40
—
1
0.5
0.8
—
—
GA
200
40
—
—
—
—
0.8
0.05Table 8
Performance of the proposed ABC, PSO, GA, and ES based on 50 experiments.
Performance
ABCalgorithm
PSOalgorithm
GAalgorithm
Enumerationsearch
Accuracy
95.47%
87.99%
96.45%
100%
FBS
14%
5%
16%
100%
CPU usage
5.23
4.15
6.84
83.62
Memory
424
424
424
53
FBS: Frequency to find best solution; Memory: required memory space to represent a solution.Figure 15
The average convergence curves of the proposed ABC, PSO, and GA.


## 5.1. Presentation of Case Problem

The Shuibuya hydropower project is located in Badong County in the middle reaches of the Qingjiang River. It is the first cascaded project on the Qingjiang mainstream. At 233 m tall and containing 15,640,000 m3 of material, it is the tallest concrete face rockfill dam in the world. This is a large scale project which has a complex transportation network for construction materials within the construction yard. In the network, large quantities of hazardous materials such as explosives are also shipped every day. However, there is no specified network for hazmat transportation. As a result, the project managers and public in the yard are concerned about the potential risk. Therefore, it is necessary to design a hazmat transportation network in the construction yard. In this case, the project manager is the decision maker on the upper level who hopes to mitigate the hazmat transportation risk by designing a new road network. Meanwhile, he also considers the decision making of the carriers and the emergency response department. Here, the emergency response department is a union of the fire-fighting department, the first-aid department, and the police station.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 has 17 preexisting trunk roads and 9 temporary roads located on the left and right banks forming a solid network of cycle traffic. There is a Cross River bridge connecting the left and right banks. The actual construction floor plan is as in Figure11. In order to apply the proposed methods more conveniently, adjacent roads of the same type have been combined and the road shapes have been ignored. An abstracted transportation network illustration is in Figure 12. The illustration has 37 nodes and 53 links in total and the application of the proposed methodologies does not affect the result of the case problem.Figure 11
The floor plan of the Shuibuya project.Figure 12
An abstracted network of the Shuibuya project.In the network, there are two hazardous materials warehouses (node no. 2 and no. 23), eight main demand nodes (node no. 1, no. 7, no. 9, no. 16, no. 18, no. 25, no. 32, and no. 36), and eight transportation commodities. Detailed data on the transportation commodities is shown in Table2.Table 2
Detailed data on the eight transportation commodities.
Commodityk
1
2
3
4
5
6
7
8
Sources(k)
No. 2
No. 2
No. 2
No. 2
No. 2
No. 23
No. 23
No. 23
Destinationt(k)
No. 1
No. 7
No. 9
No. 16
No. 18
No. 25
No. 32
No. 36
d
k (kg)
600
420
1500
2000
800
420
2000
480


## 5.2. Data Collection and Computing Results

In this case, all data on the transportation network, hazardous material (explosive), and emergency response were obtained from the Hubei Qingjiang Shuibuya Project Construction Company. The transportation data is shown in Table3. However, there was no ready-made data on transportation risk. In this paper, transportation risk is described as a Fu-Fu variable composed of two fuzzy factors: fuzzy accident possibility and fuzzy accident consequences. In a construction project, an accident may lead to various losses such as fatalities, damage to road networks and construction facilities, project duration delay, and social impacts. In this paper, fatalities, project duration delay, and damage to road network and construction facilities are considered accident consequences. The data fuzzification method outlined in Section 2 is used to obtain the accident probability and consequences of fuzzy data. Based on this method, and considering the convenience of use and interpretation, the probabilities are modelled as discrete fuzzy numbers and the consequences are transformed into costs and modelled as triangle fuzzy numbers. Hence, the unit transportation risk for each road section is modelled as a discrete triangular Fu-Fu number, and the detailed data is shown in Table 4.Table 3
Data on the road network of Shuibuya hydropower project.
Links
Nodes
Distant(miles)
Cost(CNY)
1
#1, #4
0.95
5.2
2
#2, #3
0.76
3.8
3
#2, #4
0.86
4.8
4
#3, #6
1.98
9.5
5
#4, #5
0.98
5.0
6
#4, #7
0.75
3.0
7
#5, #6
0.27
1.2
8
#6, #9
0.87
4.0
9
#7, #8
1.02
5.0
10
#8, #9
0.95
5.0
11
#9, #10
0.53
2.5
12
#10, #16
0.82
3.6
13
#7, #11
1.78
9.6
14
#8, #13
0.85
3.5
15
#9, #13
0.82
4.6
16
#11, #12
0.22
1.0
17
#12, #13
0.67
3.5
18
#13, #14
0.67
3.5
19
#10, #14
0.87
4.6
20
#14, #15
1.33
6.2
21
#15, #16
0.33
2.1
22
#11, #17
0.22
2.0
23
#17, #18
1.90
10.8
24
#12, #18
0.33
2.1
25
#14, #19
2.92
10.3
26
#18, #19
1.58
6.8
27
#17, #20
5.24
29.7
28
#18, #20
4.11
26.2
29
#19, #33
0.86
2.5
30
#21, #22
0.90
4.0
31
#22, #23
0.88
2.5
32
#21, #24
1.69
7.5
33
#21, #26
0.55
2.5
34
#22, #26
0.60
3.0
35
#25, #26
0.55
2.9
36
#24, #25
0.46
3.6
37
#25, #29
0.27
1.9
38
#29, #30
0.68
3.2
39
#24, #28
0.22
1.9
40
#28, #29
0.53
3.8
41
#28, #30
0.68
3.2
42
#30, #31
1.05
5.4
43
#31, #32
0.27
0.9
44
#23, #27
0.70
3.9
45
#27, #32
2.00
12.9
46
#25, #27
0.76
3.8
47
#32, #33
1.32
7.6
48
#32, #34
2.00
14.8
49
#34, #36
1.56
8.6
50
#33, #35
0.36
2.1
51
#35, #36
2.22
14.6
52
#34, #35
0.81
5.7
53
#36, #37
0.86
5.8Table 4
Transportation risk on each road section.
Links
Transportation risk (CNY/kg)
Links
Transportation risk (CNY/kg)
Links
Transportation risk
(
6,12,18
)
,
0.2
(
6,9
,
12
)
,
0.3
(
12,15,18
)
,
0.3
(
8,16,24
)
,
0.7
(
8,12,16
)
,
0.7
(
16,20,24
)
,
0.7
1
(
10,20,30
)
,
1.0
2, 3
(
10,15,20
)
,
1.0
4
(
20,25,30
)
,
1.0
(
12,24,36
)
,
0.6
(
12,18,24
)
,
0.6
(
24,30,36
)
,
0.6
(
14,28,42
),0.3
(
14,21,28
)
,
0.4
(
28,35,42
)
,
0.4
(
24,32,40
)
,
0.2
(
48,60,72
)
,
0.3
(
8,12,16
)
,
0.3
(
30,40,50
)
,
0.7
(
56,70,84
)
,
0.7
(
10,15,20
)
,
0.7
5, 6
(
36,48,60
)
,
1.0
7, 8
(
64,80,96
)
,
1.0
9
(
12,18,24
)
,
1.0
(
42,56,70
)
,
0.6
(
72,90,108
)
,
0.6
(
14,21,28
)
,
0.6
(
48,64,80
)
,
0.3
(
80,100,120
)
,
0.4
(
16,24,32
)
,
0.4
(
8,12,16
)
,
0.3
(
80,90,100
)
,
0.3
(
8,10,12
)
,
0.3
(
10,15,20
)
,
0.7
(
96,108,120
)
,
0.7
(
12,15,18
)
,
0.7
10
(
12,18,24
)
,
1.0
11, 12
(
112,126,140
)
,
1.0
13
(
16,20,24
)
,
1.0
(
14,21,28
)
,
0.6
(
128,144,160
)
,
0.6
(
20,25,30
)
,
0.6
(
16,24,32
)
,
0.4
(
144,162,180
)
,
0.4
(
24,30,36
)
,
0.4
(
8,10,12
)
,
0.3
(
8,10,12
)
,
0.3
(
12,15,18
)
,
0.3
(
12,15,18
)
,
0.7
(
12,15,18
)
,
0.7
(
16,20,24
)
,
0.7
14
(
16,20,24
)
,
1.0
15
(
16,20,24
)
,
1.0
16, 17
(
20,25,30
)
,
1.0
(
20,25,30
)
,
0.6
(
20,25,30
)
,
0.6
(
24,30,36
)
,
0.6
(
24,30,36
)
,
0.4
(
24,30,36
)
,
0.4
(
28,35,42
)
,
0.4
(
12,15,18
)
,
0.3
(
24,32,40
)
,
0.2
(
12,15,18
)
,
0.3
(
16,20,24
)
,
0.7
(
30,40,50
)
,
0.7
(
16,20,24
)
,
0.7
18, 19
(
20,25,30
)
,
1.0
20, 21
(
36,48,60
)
,
1.0
22, 23
(
20,25,30
)
,
1.0
(
24,30,36
)
,
0.6
(
42,56,70
)
,
0.6
(
24,30,36
)
,
0.6
(
28,35,42
)
,
0.4
(
48,64,80
)
,
0.3
(
28,35,42
)
,
0.4
(
12,15,18
)
,
0.3
(
18,21,24
)
,
0.3
(
12,15,18
)
,
0.3
(
16,20,24
)
,
0.7
(
24,28,32
)
,
0.7
(
16,20,24
)
,
0.7
24, 25
(
20,25,30
)
,
1.0
26
(
30,35,40
)
,
1.0
27, 28
(
20,25,30
)
,
1.0
(
24,30,36
)
,
0.6
(
36,42,48
)
,
0.6
(
24,30,36
)
,
0.6
(
28,35,42
)
,
0.4
(
42,49,56
)
,
0.4
(
28,35,42
)
,
0.4
(
6,12,18
)
,
0.2
(
12,15,18
)
,
0.3
(
6,9
,
12
)
,
0.3
(
8,16,24
)
,
0.7
(
16,20,24
)
,
0.7
(
8,12,16
)
,
0.7
29
(
10,20,30
)
,
1.0
30, 33
(
20,25,30
)
,
1.0
31
(
10,15,20
)
,
1.0
(
12,24,36
)
,
0.6
(
24,30,36
)
,
0.6
(
12,18,24
)
,
0.6
(
14,28,42
)
,
0.3
(
28,35,42
)
,
0.4
(
14,21,28
)
,
0.4
(
6,9
,
12
)
,
0.3
(
6,12,18
)
,
0.2
(
6,12,18
)
,
0.2
(
8,12,16
)
,
0.7
(
8,16,24
)
,
0.7
(
8,16,24
)
,
0.7
32
(
10,15,20
)
,
1.0
34, 35
(
10,20,30
)
,
1.0
36, 40
(
10,20,30
)
,
1.0
(
12,18,24
)
,
0.6
(
12,24,36
)
,
0.6
(
12,24,36
)
,
0.6
(
14,21,28
)
,
0.4
(
14,28,42
)
,
0.3
(
14,28,42
)
,
0.3
(
12,15,18
)
,
0.3
(
6,9
,
12
)
,
0.3
(
6,12,18
)
,
0.2
(
16,20,24
)
,
0.7
(
8,12,16
)
,
0.7
(
8,16,24
)
,
0.7
37, 38
(
20,25,30
)
,
1.0
39, 41
(
10,15,20
)
,
1.0
42, 43
(
10,20,30
)
,
1.0
(
24,30,36
)
,
0.6
(
12,18,24
)
,
0.6
(
12,24,36
)
,
0.6
(
28,35,42
)
,
0.4
(
14,21,28
)
,
0.4
(
14,28,42
)
,
0.3
(
12,15,18
)
,
0.3
(
6,9
,
12
)
,
0.3
(
6,12,18
)
,
0.2
(
16,20,24
)
,
0.7
(
8,12,16
)
,
0.7
(
8,16,24
)
,
0.7
44–46
(
20,25,30
)
,
1.0
47
(
10,15,20
)
,
1.0
48, 49
(
10,20,30
)
,
1.0
(
24,30,36
)
,
0.6
(
12,18,24
)
,
0.6
(
12,24,36
)
,
0.6
(
28,35,42
)
,
0.4
(
14,21,28
)
,
0.4
(
14,28,42
)
,
0.3
(
12,15,18
)
,
0.3
(
6,9
,
12
)
,
0.3
(
16,20,24
)
,
0.7
(
8,12,16
)
,
0.7
50, 51
(
20,25,30
)
,
1.0
52, 53
(
10,15,20
)
,
1.0
(
24,30,36
)
,
0.6
(
12,18,24
)
,
0.6
(
28,35,42
)
,
0.4
(
14,21,28
)
,
0.4In addition, the emergency response department plans to locate two emergency response teams for potential hazmat transportation accidents while six nodes are optional for this location (i.e., no. 7, no. 9, no. 13, no. 23, no. 25, and no. 32). The service range of the emergency response teams is 2 kilometers. Based on this coverage distance, all covered links by each optional location are predetermined. Hence, the setNij is also determined. In order to achieve the maximal coverage, if a link is covered partially, then it is divided into two links (one is entirely covered, and another is not covered) by adding a new node.Using this data, after running the computer program for the improved ABC algorithm using MATLAB 2007, the solutions for the case problem and the efficiency of the proposed algorithm were obtained.The algorithm and model parameters for the case problem were set as follows: the number of the food sourcesSN=40, the value of limit=20, the maximum cycle number MCN=200, and the mitigated risk coefficient for emergency response α=0.3. The computer running environment was an Intel Core 2 Duo 2.26 GHz clock pulse with 2048 MB memory. After 10.46 minutes on average, the optimal solutions for the bilevel programming were determined. The optimal transportation network is shown in Figure 13.Figure 13
The optimal network for the explosives transportation.


## 5.3. Model Analysis

### 5.3.1. Two Different Networks

As discussed, the authority needs to consider the emergency response network when the hazardous materials transportation network is designed, because an emergency response network has an impact on the transportation network.In this study, a comparison is given for a network considering emergency response and one without regard to emergency response. Figure14 illustrates the two different networks. From this, it can be seen that there are some distinct differences between the two networks. Links No. 7, No. 8, No. 15, No. 16, and No. 17 would be selected if emergency response were considered though their unit risks are higher than other links, while they would not be included in the network without considering the emergency response. In addition, there is also a lower risk (806.44 thousand CNY and 1049.06 thousand CNY in the two networks, resp.) if the emergency response is considered when designing the transportation network. Therefore, the total transportation risk is reduced by 23.13% due to the existence of emergency response teams.Figure 14
Difference between the two proposed networks.


### 5.3.2. Sensitivity Analysis toα

In Section3, it is assumed that the risk of link (i,j) can be reduced by α(Dmax-Dijq)hijqrijk if it is covered by an emergency response team located at node q. It is known that the value of parameter α is a key factor impacting the results. In fact, the decision maker is able to adjust the parameter to obtain different solutions. A sensitivity analysis is represented in Table 5. From Table 5, it can be seen that total risk decreases as the value of α increases. The designed networks are also related to α. With an adjustment in α, three different optimal networks were found. Each optimal network is related to a value range of α. It shows that the optimal network is not sensitive to α and a range of α is enough for the proposed network designing problem. Hence, the authority only needs to determine the range of α when designing the transportation network.Table 5
Sensitivity analysis of parameterα.
Value ofα
Total risk
Location
Road links for hazmat transportation
α
=
0.00
1049.06
—
{1, 2, 3, 5, 6, 9, 10, 13, 14, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.05
1016.37
#7, #32
{1, 2, 3, 5, 6, 9, 10, 13, 14, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.10
983.67
#7, #32
{1, 2, 3, 5, 6, 9, 10, 13, 14, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.15
949.17
#7, #32
{1, 2, 3, 5, 6, 9, 10, 13, 14, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.20
901.59
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.25
854.02
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.30
806.44
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.35
758.87
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.40
711.29
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.45
663.19
#9, #25
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 37, 38, 42, 43, 48, 49}
α
=
0.50
612.34
#9, #25
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 37, 38, 42, 43, 48, 49}


### 5.3.3. Uncertainty Analysis

In addition, uncertainty is also an important consideration in this study. Fu-Fu variables are used to model the transportation risk because of the lack of precise data such as the accident probabilities and the accident consequences on specific roads, so the proposed model contains a complex fuzzy factor. Besides the model considering a Fu-Fu factor, there are two other kinds considered, one without uncertainty and one that includes one fuzzy factor. Truly determinate data does not exist in project practice, so the managers make decisions from historical statistical data and their own experience. This factor makes it very difficult to select decision-making data with most of it obtained by ignoring the uncertainty and using an average value or data chosen at random. If only one fuzzy factor is considered such as the fuzzy accident consequences, then other uncertain information needs to be ignored, such as the accident probabilities. Then the accident probabilities have to be randomly chosen from the given interval. All these will lead to an imprecise solution. Table6 shows a comparison result of the three model types. In the comparison, the model considering a complex fuzzy factor finds different solutions by adjusting the optimistic and pessimistic index λ when the others find different solutions by choosing different transportation risk randomly. It is clear that the proposed model considering the Fu-Fu factors has a much better performance than the others, not only in the average value of the results, but also in their stability.Table 6
Comparisons among the three model types based on the risk-minimizing objective.
Type
Complex fuzzy type
Fuzzy type
Determined type
Best result
736.68
720.22
698.36
Worst result
854.35
902.84
932.07
Average result
806.04
834.52
842.36


## 5.3.1. Two Different Networks

As discussed, the authority needs to consider the emergency response network when the hazardous materials transportation network is designed, because an emergency response network has an impact on the transportation network.In this study, a comparison is given for a network considering emergency response and one without regard to emergency response. Figure14 illustrates the two different networks. From this, it can be seen that there are some distinct differences between the two networks. Links No. 7, No. 8, No. 15, No. 16, and No. 17 would be selected if emergency response were considered though their unit risks are higher than other links, while they would not be included in the network without considering the emergency response. In addition, there is also a lower risk (806.44 thousand CNY and 1049.06 thousand CNY in the two networks, resp.) if the emergency response is considered when designing the transportation network. Therefore, the total transportation risk is reduced by 23.13% due to the existence of emergency response teams.Figure 14
Difference between the two proposed networks.


## 5.3.2. Sensitivity Analysis toα

In Section3, it is assumed that the risk of link (i,j) can be reduced by α(Dmax-Dijq)hijqrijk if it is covered by an emergency response team located at node q. It is known that the value of parameter α is a key factor impacting the results. In fact, the decision maker is able to adjust the parameter to obtain different solutions. A sensitivity analysis is represented in Table 5. From Table 5, it can be seen that total risk decreases as the value of α increases. The designed networks are also related to α. With an adjustment in α, three different optimal networks were found. Each optimal network is related to a value range of α. It shows that the optimal network is not sensitive to α and a range of α is enough for the proposed network designing problem. Hence, the authority only needs to determine the range of α when designing the transportation network.Table 5
Sensitivity analysis of parameterα.
Value ofα
Total risk
Location
Road links for hazmat transportation
α
=
0.00
1049.06
—
{1, 2, 3, 5, 6, 9, 10, 13, 14, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.05
1016.37
#7, #32
{1, 2, 3, 5, 6, 9, 10, 13, 14, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.10
983.67
#7, #32
{1, 2, 3, 5, 6, 9, 10, 13, 14, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.15
949.17
#7, #32
{1, 2, 3, 5, 6, 9, 10, 13, 14, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.20
901.59
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.25
854.02
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.30
806.44
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.35
758.87
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.40
711.29
#9, #32
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 44, 45, 48, 49}
α
=
0.45
663.19
#9, #25
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 37, 38, 42, 43, 48, 49}
α
=
0.50
612.34
#9, #25
{1, 2, 3, 5, 6, 7, 8, 15, 16, 17, 18, 20, 21, 23, 31, 34, 35, 37, 38, 42, 43, 48, 49}


## 5.3.3. Uncertainty Analysis

In addition, uncertainty is also an important consideration in this study. Fu-Fu variables are used to model the transportation risk because of the lack of precise data such as the accident probabilities and the accident consequences on specific roads, so the proposed model contains a complex fuzzy factor. Besides the model considering a Fu-Fu factor, there are two other kinds considered, one without uncertainty and one that includes one fuzzy factor. Truly determinate data does not exist in project practice, so the managers make decisions from historical statistical data and their own experience. This factor makes it very difficult to select decision-making data with most of it obtained by ignoring the uncertainty and using an average value or data chosen at random. If only one fuzzy factor is considered such as the fuzzy accident consequences, then other uncertain information needs to be ignored, such as the accident probabilities. Then the accident probabilities have to be randomly chosen from the given interval. All these will lead to an imprecise solution. Table6 shows a comparison result of the three model types. In the comparison, the model considering a complex fuzzy factor finds different solutions by adjusting the optimistic and pessimistic index λ when the others find different solutions by choosing different transportation risk randomly. It is clear that the proposed model considering the Fu-Fu factors has a much better performance than the others, not only in the average value of the results, but also in their stability.Table 6
Comparisons among the three model types based on the risk-minimizing objective.
Type
Complex fuzzy type
Fuzzy type
Determined type
Best result
736.68
720.22
698.36
Worst result
854.35
902.84
932.07
Average result
806.04
834.52
842.36


## 5.4. Algorithm Evaluation

In this paper, an improved artificial bee colony algorithm with priority-based encoding has been proposed. In order to test the efficiency of the algorithm, a comparison with other solution methods was conducted.The proposed problem is modeled as a linear bilevel programming model in this paper. The most common solution strategy for linear bilevel problems is to transform the bilevel model into a single one by the use of its Karush-Kuhn-Tucker (KKT) conditions. However, it is difficult when there are binary variables in the inner models. This also results in that it cannot be solved by common commercial solvers. In theory, it is possible to solve this problem by conducting an enumeration search (ES) over the design variablesyij on the upper level and solving a series of linear programs on the lower level for each selection of design variables. However, this search would be conducted 253 times because of the 57 binary design variables. This is impractical for large-scale problems. Hence, it is not appropriate to solve the problem by KKT conditions and enumeration search in most cases. Fortunately, it is found that the transportation network in this case can be divided into two subnetworks along the river because all these hazmats transportation would not go through the river. Moreover, some road links such as (1,4), (2,3), (17,20), (18,20), and (36,37)  can be determined by observing the network. Therefore, the network can be divided into two networks with 25 and 23 links, respectively. The problem in the case study can be solved by conducting 223+225 searches. In experiments, it spends 186 minutes in running the search program. Although the solution is optimal, the too long time is not acceptable.In order to show the efficiency of the proposed ABC algorithm, a comparative experiment among the proposed ABC algorithm, a particle swarm optimization (PSO), a genetic algorithm (GA), and the enumeration search method was also conducted. In the experiment, the same encoding method (i.e., priority-based encoding) was used for the three algorithms. For the GA, the partially mapped crossover and local search-based mutation method are used. For the PSO, a hybrid particle-updating mechanism is used in the PSO. The other parameters for the algorithms are shown in Table7. The experiments for the algorithms were carried out over 50 times. Figure 15 shows the convergence histories for the three heuristic algorithms. The detailed performance of the four algorithms is stated in Table 8. The results indicate that the ABC needs less time than the GA and has a more stable tendency than the PSO. Moreover, the accuracy of both the ABC and the GA is very high; even the optimal solution can be found with a percentage more than 14%. Therefore the performance of the ABC is on par with that of both the PSO and the GA. The proposed algorithm is also useful and efficient for solving the proposed case.Table 7
Parameter values used in the experiments.
MCN
Pop size
Limit
Weight
c
p
best
>cgbest
r
crossover
r
mutation
ABC
200
40
20
—
—
—
—
—
PSO
200
40
—
1
0.5
0.8
—
—
GA
200
40
—
—
—
—
0.8
0.05Table 8
Performance of the proposed ABC, PSO, GA, and ES based on 50 experiments.
Performance
ABCalgorithm
PSOalgorithm
GAalgorithm
Enumerationsearch
Accuracy
95.47%
87.99%
96.45%
100%
FBS
14%
5%
16%
100%
CPU usage
5.23
4.15
6.84
83.62
Memory
424
424
424
53
FBS: Frequency to find best solution; Memory: required memory space to represent a solution.Figure 15
The average convergence curves of the proposed ABC, PSO, and GA.


## 6. Conclusion

In this paper, a bilevel optimization model for an integrated hazardous materials transportation and emergency response network design was proposed. In the model, three decisions were considered. On the upper level, the authority (the leader) designs the transportation network with the criterion of minimizing total transportation risk. On the lower level, the carriers and the emergency response department (the followers) make their decision based on the leader’s decision. The carriers first choose their routes so that total transportation cost is minimized. Then, the emergency response department locates their emergency service teams so as to maximize the total weighted arc length covered. In contrast to prior studies, the uncertainties associated with transportation risk were explicitly considered in the objective function of the mathematical model. Specifically, this research uses a Fu-Fu variable to model the transportation risk and an expected value operation is proposed to transform the uncertain risk to a determinate one. Then, an improved artificial bee colony algorithm was applied to search for the optimal solution of the bilevel model. Finally, the efficiency of the proposed model and algorithm was evaluated using a practical case and various computing attributes. Two comparisons for the model were conducted: one looking at three uncertainty types and the other between the networks taking into consideration the emergency response and the other without this consideration. The results show that it is significant to consider a network design problem with emergency response in a complex fuzzy environment. The efficiency of the proposed algorithm was also evaluated by comparing it with the GA, the PSO, and the enumeration search method.The area for future research has four aspects. Firstly, investigate the methods to deal with other types of uncertain risks for hazmat transportation such as the uncertainty risk including fuzzy factor and random factor at the same time Secondly, apply the proposed model and algorithm to more complex road network such as an urban traffic network. Thirdly, develop more efficient heuristic methods to solve such bilevel problems. Fourthly, consider other cases for controlling hazmat transportation risk, for example, only to close the roads in some time range. Each of these areas is very important and equally worthy of attention.



---

*Source: 517372-2013-05-07.xml*