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

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 tomaximize 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.


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][5][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][17][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 [, ] 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 Table 1), 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.

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.

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 Figure 2. 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.

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.

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 (, ) is often represented as where   is the accident probability on link (, ) and   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 ×  dollars/mile." Let the fuzzy accident probability p be given by a domain   and p = {(,  p ()) |  ∈   }.Then the fuzzy risk can be given as r ,  r (    ) =  pj (  ),  ∈   , where  p and  r are the membership functions of fuzzy accident probability p and fuzzy transportation risk   , 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] with membership function 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.

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 [ min ,  max ] 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.
Main fuzzification steps for accident probability are as follows.
Step 2. Determine the weight of these impacts using analytic hierarchy process method,  = [ 1 ,  2 , . . .,   ].Step 3. Calculate the fuzzy relation matrix Ψ = (  ) × ,   =   /  , where   is the number of experts invited to nominate their evaluation to a specific term for impact factor  and   is the nominated number of linguistic terms  for factor .
Step 5. Describe these linguistic terms with fuzzy numbers H = [ h 1 , h 2 , . . ., h  ], and calculate the fuzzy results p =  ∘ H, where h  is the fuzzy number associated with comment term V  , 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 probability p is modelled using a discrete fuzzy number and the consequence C is modelled using a triangular fuzzy number C  = [   ,    ,    ], then the risk can be described as (4)

Modelling
The hazmat network design problem is a graph theoretical problem defined on a directed graph  = (, ), where  is the set of vertices and  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  commodities between their respective origins and destinations.Each commodity corresponds to an OD pair.Let ((), ()) be the OD pair of commodity ,  ∈ {1 ⋅ ⋅ ⋅ } and let   be the corresponding number of shipments.The parameters   and   refer to the risk and cost associated with a unit flow of commodity  on arc (, ), respectively.Each link is assumed to be a network segment.

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.
Example 3. ξ = (, γ, ) with γ ∼ (  ,   ,   ) is called Fu-Fu variable, see Figure 4, if the outer-layer and inner-layer membership functions are as follows: where γ is the center of ξ, which is a triangular fuzzy variable, and  and  are the smallest possible value and the largest possible value of ξ.   ,   , and   are the smallest possible value, the most promising value, and the largest possible value of γ, respectively.Definition 4 (see [22]).The expected value of a Fu-Fu variable is defined by 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) [()] and [()] are independent fuzzy variables, then for any real numbers  and , one has Lemma 6.If transportation risk r is a Fu-Fu variable characterized as follows (see Figure 5): where [  where the weights   ,  = 1, 2, . . .,  are given by Proof.For any  = 1, 2, . . ., , r ] is a triangular fuzzy variable.It follows from the definition of expected value of fuzzy variable that we have From Definition 4 and Theorem 5, Here, r  () ( = 1, 2, 3) is a discrete fuzzy variable whose membership function is given by From the definition of expected value of fuzzy variable, we have where And then, This completes the proof.

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 are   .
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 cleanup 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  to a demand link (, ) can be described using the shortest distance from the midpoint of link (, ) to node , that is,    .Meanwhile, this distance cannot exceed the maximum service distance  max of the emergency response teams.Therefore, if link (, ) is covered by node , that is, ℎ   = 1, then the reduced risk can be described as a function of the distance    and  max .For convenience, it is assumed that the reduced risk and the travel distance meet a linear relation such as Here, ( max −    )ℎ   describes the percentage of reduced risk when link (, ) is covered by emergency response team located in node . is a coefficient related to the category of undesirable accident consequences and the power of the emergency response teams. max 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  max 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 In the risk described previously, r 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:  In the expected total risk, only r is a Fu-Fu variable.Based on Definition 4 and Theorem 5, the expected total risk can be transformed into following equation: For this objective, the variables   and ℎ  are solved in the lower level model.Hence, the network design model can be modelled as follows: where   and ℎ  are solved in the lower-level model. (21)

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 ∑  =1 ∑ (,)∈         .For the flow problem, two constraints must be met.One is the flow equilibrium constraint which ensures the flow of commodity  from its origin to its destination: 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: The objective and constraints constitute the routes choice model of the carriers as follows: 1,  =  () , ∀ ∈ ,  = 1, 2, . . .,  0, otherwise   ≤   ,   ∈ {0, 1} , ∀ (, ) ∈ ,  = 1, 2, . . ., . (24)

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 firstaid 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 (, ) describes the gain when link (, ) is covered by emergency response teams, that is, the reduced risk ∑ ∈ ∑  =1     ( max −    ) r   .The expected value operation is used to deal with the Fu-Fu risk r ; where ℎ   takes value 1 if (, ) is covered by emergency response team located at note  and 0 otherwise.Then the objective can be modelled as follows: At the same time, two constraints must also be considered.First, all located teams must be equal to its sum.Let   describe whether a team is located in node , and  max describes the total number of located emergency response teams; then the constraint can be stated as follows: Moreover, a link (, ) is covered by node  only if a team is located at node ; then this constraint can be stated as follows: 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: Finally, the maximum travel distance of an emergency response team to any link traversed by hazardous materials within its jurisdiction should not exceed the threshold  max .It can be described as the following constraint: From the aforementioned, the maximal arc-covering location problem can be modelled as max 3.5.Bilevel Programming Model.As outlined in Section 2, 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: 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 ] is modelled as a discrete triangular Fu-Fu variable (see ( 9)), and its expected value can be calculated by (10).

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 NPhard 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.

Introduction to ABC Algorithm.
The artificial bee colony (ABC) algorithm, proposed by Karaboga in 2005 for realparameter 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][38][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.

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 lowerlevel 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 Figure 7 and the main procedure is presented as follows.
Main procedure for the proposed ABC algorithm is as follows.
Step 1. Initialize the number of food sources SN, the maximum cycle number MCN, and the control parameter for the abandoned food source limit.Let iteration  = 1.For

Begin
Step 1: initialize the food source and parameters Step 2: encode the positions of food sources into solutions  Step 2. Encode each food source position   into the solution using priority-based encoding method: the set of upper-level decision variables {  }, that is, the designed transportation.
Step 3. Evaluate each food source site   by solving the lowerlevel 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 position   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   , the employed bee memorizes it as the new position   .
(4.2) An onlooker bee evaluates the fitness value information fit  of food source  taken from all employed bees and modifies the food source site with a probability   related to its fitness value.If the position of the modified food source is better than   , replace it.
(4.3) Determine an abandoned solution for the scout.If the position of food source  is not modified after the employed bees and onlooker bees phases, its trials counter cou  is increased by 1; otherwise, the counter is reset to 0. Check the counter cou  : if cou  > limits, replace   with a new randomly produced position and reset cou  to 0.
Step 5.If the stopping criterion is met, that is,  > MCN, stop.Otherwise, let  =  + 1 and go to Step 2.

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 Node priority PSO [40] and GA [41].Here, it is proposed to encode the paths into a food source.

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 to .Hence, the position of each food source   can be represented as a matrix with  rows and  columns.That is,   = [ 11 ,  12 ,. ..,  1 ;  21 ,  22 ,. ..,  2 ;. ..;  1 ,  2 ,. ..,   ], where each row   = [ 1 ,  2 , . . .,   ] is a priority sequence of 1 to .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.

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 node () 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 () 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 ( = 1, 2, 3), three paths can be easily determined using the given encoding method: ( In the equation,   and ℎ  are dependent on the lowerlevel 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,   is determined.After this, the "bintprog" function in Matlab is used to solve the emergency response teams location model.Step 1: choose food source  randomly Step 2: choose modified position  randomly:  = 4 Step 3: update the priority of chosen position

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.

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 a -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  rows and  columns (i.e.,  random sequences from 1 to ), 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 , 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   , a food source position   ( ̸ = ) is first randomly chosen.Then for every row , a position  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 , that is,   =   , the position whose value is equal to   must be found and its value replaced by   .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.

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 probability   related to its fitness value.The probability value is calculated using the following expression: produces a position modification where fit  is the fitness value of the solution  evaluated by the employed bee.
After that, a parameter value is randomly generated between 0 and 1.If it is smaller than   , 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   and another food source  randomly.Then it will replace the path with the similar path from another food source, that is,   =   .An example is given in Figure 10.

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 counter cou  is used.If the position of food source  is not modified after the employed bees and onlooker bees phases, its trials counter cou  is increased by 1; otherwise, the counter is reset to 0. Check the counter cou  : if cou  > limits, replace   with a new randomly produced position and reset cou  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.

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.,640,000 m 3 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 Figure 11.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 Mathematical Problems in Engineering 15  and the application of the proposed methodologies does not affect the result of the case problem.
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 Table 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 Table 3.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

Qingjiang River
The roads used to transport hazardous explosives The sites selected to locate emergency response teams discrete triangular Fu-Fu number, and the detailed data is shown in Table 4.
In 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 set   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 sources SN = 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.

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.Figure 14 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.

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 variables   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 2 53 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 2 23 + 2 25 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 Table 7.
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.

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.

Figure 1 :
Figure 1: The structure of the remainder of this paper.

Figure 2 :
Figure 2: The abstracted structure of the proposed bilevel hazmat transportation network design problem.

Figure 3 :
Figure 3: The fuzzification of accident probability linguistic terms.

Figure 5 :
Figure 5: An illustration of Fu-Fu fuzzy risk.

Figure 6 :
Figure 6: Estimation of the value of .

Step 3 :
calculate the fitness values and memorize the position of the best food source The stopping criterion is met?End Step 4.1: search the neighborhood by employed bees and conduct a position modification Step 4.3: produce new food sources randomly by scouts Yes Find abandoned food sources?No No Yes Step 4: update the food sources Step 4.2: select food source sites by onlookers with a probability

Figure 7 :
Figure 7: The flow chart of the proposed algorithm.

Figure 8 :
Figure 8: An example of food source representation in a network with 10 nodes.

Figure 9 :
Figure 9: An example of neighbourhood search of employed bees.

Figure 10 :
Figure 10: An example of the neighbourhood search for onlooker bees.

Figure 11 :
Figure 11: The floor plan of the Shuibuya project.

Figure 12 :
Figure 12: An abstracted network of the Shuibuya project.

Figure 13 :
Figure 13: The optimal network for the explosives transportation.

Figure 14 :
Figure 14: Difference between the two proposed networks.

Figure 15 :
Figure 15: The average convergence curves of the proposed ABC, PSO, and GA.

Table 2 :
Detailed data on the eight transportation commodities.

Table 3 :
Data on the road network of Shuibuya hydropower project.
Analysis to .In Section 3, it is assumed that the risk of link (, ) can be reduced by ( max −

Table 4 :
Transportation risk on each road section.

Table 7 :
Parameter values used in the experiments.MCN Pop size Limit Weight  best > best  crossover  mutation

Table 8 :
Performance of the proposed ABC, PSO, GA, and ES based on 50 experiments.
FBS: Frequency to find best solution; Memory: required memory space to represent a solution.