A Risk-Averse Shelter Location and Evacuation Routing Assignment Problem in an Uncertain Environment

Disasters such as hurricanes, earthquakes and floods continue to have devastating socioeconomic impacts and endanger millions of lives. Shelters are safe zones that protect victims from possible damage, and evacuation routes are the paths from disaster zones toward shelter areas. To enable the timely evacuation of disaster zones, decisions regarding shelter location and routing assignment (i.e., traffic assignment) should be considered simultaneously. In this work, we propose a risk-averse stochastic programming model with a chance constraint that takes into account the uncertainty in the demand of disaster sites while minimizing the total evacuation time. The total evacuation time reflects the efficacy of emergency management from a system optimal (SO) perspective. A conditional value-at-risk (CVaR) is incorporated into the objective function to account for risk measures in the presence of uncertain post-disaster demand. We resolve the non-linear travel time function of traffic flow by employing a second-order cone programming (SOCP) approach and linearizing the non-linear chance constraints into a new mixed-integer linear programming (MILP) reformulation so that the problem can be directly solved by state-of-the-art optimization solvers. We illustrate the application of our model using two case studies. The first case study is used to demonstrate the difference between a risk-neutral model and our proposed model. An extensive computational study provides practical insight into the proposed modeling approach using another case study concerning the Black Saturday bushfire in Australia.


Introduction
Whether they are natural or man-made, disasters such as hurricanes, earthquakes and terrorist attacks inflict serious risk on humankind. In 2010, the Haitian government reported that that an estimated 316,000 people had died, 300,000 had been injured and 1,000,000 were homeless as a result of the Haiti earthquake. This caused a massive level of devastation and imposed enormous operational challenges on the humanitarian agencies and government. Scientific and reasonable shelter assignment and evacuation routing planning are two important safeguards to minimize the loss from disasters. Different types of disasters require different types of shelter assignments and evacuation planning. For example, earthquakes require emergency evacuation while typhoons allow early evacuation. By means of modern information technology, the relevant government department can analyze and forecast the path and intensity of a typhoon, which can give people sufficient time to evacuate. However, scientists cannot yet calculate when there will be a major earthquake. When an earthquake comes, many people need to evacuate immediately. To promote the efficiency of disasters rescue and reduce life and economic losses, shelter location is significant to the reasonable planning at the pre-disaster phase and evacuation plans are carried out relying upon the information gained (e.g., earthquake magnitude, volcanic eruptive intensity, fire dangerous grade, etc.) when disasters happen. Topics regarding shelter or evacuation, important aspects of disaster relief, have been frequently considered by researchers, and have appeared in emergency supply planning [1,2], shelter location problems [3][4][5], emergency evacuation planning problems  and shelter location and evacuation routing assignment problems [4,5,[8][9][10][11][12]16,27,28].
Shelter location and evacuation routing operations rely on the boundary between disaster preparedness and disaster response. In line with the framework proposed by [4], we assume that shelter location and evacuation routing assignments are disaster response operations.
A shelter is a safe facility that protects people from possible damaging effects of a disaster. It can provide people with necessary living conditions (e.g., food, water, medicine, etc.). People usually start leaving their own houses to seek shelters with a sense of panic when they either face or are experiencing dangerous circumstances. Under those circumstances, people may choose the same path to the shelters. This will reduce the traffic capacity of those roads and delay the time taken to reach the shelter. In addition, space within shelters is limited. Overcrowded shelters cause people who may arrive later to seek a new one. This will further aggravate the condition of people who were wounded by the disaster.
The existing evacuation planning and management models are mostly based on traffic assignment models [29], such as the user equilibrium (UE), the system optimal (SO) and the nearest allocation (NA) approaches. Recently, the Constrained System Optimal (CSO) model [4,5] was developed to provide a compromise between the SO and UE/NA approaches. The UE approach assumes that the evacuees act selfishly and have full information about the traffic conditions on every possible route in the evacuation network. Based on these assumptions, in the UE model, every evacuee can find the optimal route to their destination using an objective function of minimizing their travel time. Clearly, the perfect information assumption in the UE approach is strong and may not be valid in the case of a natural disaster, as the real traffic demand is unusual, and it is difficult for evacuees to guess the congestion level on a path. In the NA model, each evacuee uses the shortest path to reach the nearest open shelter sites, whereas this short-sighted assignment will lead to a poor system performance.
The main goal of evacuation agencies and the governments, however, is to avoid congestion and enable the evacuation of the disaster area in a timely manner. This can be achieved through effective planning strategies, the SO model, which distributes traffic over the road network, and shelter assignment. In the SO model, evacuees are assigned to designated safety and evacuation routes as soon as possible to minimize the total evacuation time. The SO solution may not favor the preferences of some evacuees to the benefit of the other evacuees. For example, some evacuees may be allocated to shelters much further away than they would like to go or to routes much longer than they would normally take.
In this article, we focus on emergency evacuation planning operations for a disaster response, which consists of shelter location decisions and evacuation routing assignments in an SO fashion. In this work, evacuation demand is an uncertain parameter, as shelter location decisions should be made before disasters happen. Contrary to the general expectation-based (i.e., risk-neutral) stochastic model, a risk-averse version that incorporates a risk measure called Conditional Value-at-Risk (CVaR) is established to capture the random parameters to minimize the total evacuation time and associated risk. Considering the utilization balance among shelters in uncertain contexts, on the other hand, chance constraints are introduced to ensure that the utilization level of shelters is not less than a threshold within a given level of reliability.
To summarize, the contribution of the present work is four-fold: (i) A two-stage stochastic programming shelter location and evacuation routing assignment (SLERA) model with a chance constraint is developed, in which the demand for evacuees is regarded as random. In the first stage, location decisions are made following some possible demand scenarios. The second-stage decision considers specific evacuation routing assignments concerning available shelters as the actual realization of demand. (ii) Considering the uncertainty of the demand, a risk measure (CVaR) is introduced into the stochastic model to balance the expectations related to the total evacuation time and inherent risk. (iii) To handle the non-linearity objective function of the traffic flow, a second-order cone programming (SOCP) approach is used. (iv) Extensive numerical experiments are carried out to provide practical insight into the proposed modeling approach.
The remainder of this work is organized as follows: Section 2 reviews some related work on emergency management, especially related to the disaster response, such as shelter locations and evacuation routing. Expectation-based and risk-averse mathematical models for the SLERA problem are given in Section 3. The solution approach is detailed in Section 4. We present the computational results of the case study in Section 5. Finally, Section 6 concludes this work and indicates some future research directions.

Literature Review
In this section, we provide a brief introduction to the existing literature, including evacuation planning, shelter location, and an integrated approach of both. Table 1 gives an overview of these works.

Emergency Evacuation Planning Problem
In contrast to conventional emergency evacuation plans that often assign evacuees to fixed routes or shelters, Han et al. [6] proposed a static model that decides on the shelter and route assignment simultaneously by not constraining evacuees to pre-specified shelter sites to ensure the quality of plans while roads are disturbed. In [18], optimal evacuation destination, traffic assignment, and evacuation departure schedule decisions were integrated into a unified optimal traffic flow optimization model.
To represent the time-varying characteristics of the evacuation traffic conditions, dynamic network flow models were employed in the evacuation planning problem [17,20,30]. Specifically, Lim et al. [20] presented a capacity-constrained network flow optimization approach for finding evacuation paths, flows and schedules to maximize the total number of evacuees for short notice evacuation planning against hurricanes. In [17], a conflict-based path-generation approach was proposed for the evacuation planning problem by decomposing the problem into a master and a sub-problem. The sub-problem generates new evacuation paths for each evacuate area, while the master problem optimizes the flow of evacuees and produces an evacuation plan. In [15], the demand management strategies of aggregate-level staging and routing were studied in an attempt to manage large regional evacuations against disaster events by reducing or eliminating congestion in the evacuation network. In [22], the dynamic resource allocation problem was addressed for transportation evacuation planning on large scale networks. In a system optimal dynamic traffic assignment problem, Tuydes-Yaman and Ziliaskopoulos [23] focused on modeling demand management strategies regarding to the optimal departure time, optimal destination choice and optimal evacuation zone scheduling (also known as staggered evacuation) under a given fixed evacuation time assumption.
Most existing evacuation planning problems (i.e., traffic assignment problems) assume that the input parameters, such as evacuation demand and capacity, are deterministic. Notably, some works have accounted for uncertainties in parameters under stochastic settings such as stochastic programming [7], robust optimization [19,24] and chance constraint programming [30,31].

Shelter Site Location Problem
Recently, Mostajabdaveh et al. [14] studied a stochastic scenario-based problem involving the determination of a set of shelter locations in the preparedness stage for natural disasters, which takes into consideration the equity quantified by Gini's Mean Absolute Difference (called Gini index) under uncertainty. For large-scale problems, a genetic algorithm to solve the sub-problem by mixed-integer programming was proposed.
Motivated by the realization that shelters are shared facilities between the source (relief supply) and demand sides, designing a strategic preparedness network consisting of the evacuation source, shelters and distribution centers [1] naturally integrates the evacuation and relief distribution problem at the strategic and tactical levels simultaneously. Their formulation can make strategic and tactical decisions at the same time before a foreseen disaster, such as a hurricane reaches the area. A benders decomposition approach is used because of the complexity of the integrated model.

Shelter Location and Evacuation Routing Assignment Problem
Within the emergency management context, the shelter location and evacuation routing problems have been studied by numerous optimization researchers, either in a separate or integrated fashion. Bi-level programming, in which shelter location decisions and route assignments are given at the upper and low levels respectively, is one of the most common approaches to address the above two problems separately in the early period [9,10,16,27,28].
Specifically, Kongsomsaksakul et al. [27] studied an optimal shelter location problem for flood evacuation planning with the assumption that the planning authority can control the traffic in certain parts of the evacuation network comprising the evacuation routes of evacuees. Bi-level programming was deployed to model this problem, where the upper lever minimized the total evacuation time and the evacuees' decisions (e.g., shelter assignment) at the lower level were formulated by a combined distribution and assignment model in line with the shelter location information of the upper level. Similarly, Ng and Waller [16] developed a hybrid bi-level model solved using a simulated annealing algorithm to balance the optimal system and individual components, where the shelter assignment is made first in a system optimal manner and then evacuees are free to choose the optimal route to reach their assigned shelter. Furthermore, bi-level programming was extended to a stochastic context by [10,28] in which a scenario-based shelter location model was proposed to identify a set of robust shelter locations during hurricane events. In [10], user equilibrium-based traffic dynamics were incorporated into the optimization of shelter locations under uncertainty to identify a set of robust shelter locations. The objective of the upper level was to minimize the weighted sum of the expected unmet shelter demand and the expected total travel time of network users. By confining the uncertain demand to an uncertainty set, Kulshrestha et al. [9] presented a robust optimization approach for both shelter location and capacity decisions under demand uncertainty and their proposed model was formulated as a mathematical program with complementarity constraints and was solved using a cutting-plane algorithm.
In order to balance different performance criteria or take into account diverse interest stakeholders, multi-objective models are always formulated [8,11,12,32,33]. A multi-objective shelter location and evacuation routing for an urban evacuation planning problem was proposed by [32], in which primary and secondary routes for evacuees are identified simultaneously. These optimization approaches are embedded into a Geographic Information System-based decision support system so planners can access them via the Internet. In [8], a three-stage multi-objective programming model was used to design a centralized emergency supply network for emergency logistics operations. A shelter network for evacuation was designed in the first-stage phase with a minimized total travel distance.
In contrast to the separate optimization method, the integrated shelter location and evacuation routing problems have received recent attention [4,5,11,12,33,34]. In the realm of emergency evacuation planning under urban area circumstances, multi-criteria optimization models were developed in [11,33]. In particular, Coutinho-Rodrigues et al. [33] introduced a multi-objective location-routing model based on the work proposed by [32] to identify the locations of shelters and evacuation paths for urban evacuation planning. The factors considered in this model are six-fold and include the total travel distance, evacuee risk, total evacuation time, and so on. This model was tested for a simulated fire situation and solved by an off-the-shelf optimization solver. In [11], a comprehensive optimization model based on the multiple commodities flow problem using public and private transport was presented to address the location and routing problems simultaneously. A genetic algorithm based on nondominated sorting genetic algorithm-II (NSGA-II) was designed to solve the multi-objective mixed-integer program, in which the authors modeled the dynamic aspects of an evacuation process and took into account the total evacuation time, risk of the evacuees, and number of shelters to be opened simultaneously. In [4,5], a non-linear mixed-integer model was presented for vehicle-evacuation toward shelter destinations, in which a constrained system optimization (CSO) methodology was introduced to guarantee a level of tolerance so that evacuees are willing to accept the travel routes, even though the assigned routes to shelters are not the shortest ones. The objective of the problem was to minimize the total evacuation time of the entire system, which was modeled through a non-linear function of traffic flow. To handle the non-convex effects of the non-linear objective function, a second-order cone programming (SOCP) method was deployed by transferring the non-linear part of the objective function to a constraint.
Furthermore, Bayram and Yaman [4] extended the CSO shelter location and routing assignment problem based on [5] to a scenario-based stochastic case where the evacuation demand and the potential alterations to the network infrastructure are uncertain before a disaster occurred. The benders decomposition algorithm with diverse cut strategies and the cutting plane algorithm was developed to deal with large-scale problems. In [34], the Turkish Red Crescent (TRC) approach, which is used to address the temporary shelter location and self-evacuation problem, was improved by developing a mathematical model to select the best shelter area locations from a set of criteria. The aim of this model is to determine the best shelter location sites and match evacuation zones to shelter sites while taking into account the utilization of shelters. Focusing on a specific category of support evacuation of late evacuees, in [12], a bi-objective integer programming model was developed to support shelter location and routing decisions during short-notice evacuations. For this model, the support of public authorities is required, and the model is solved by an -constrained algorithm with two conflicting objectives: maximizing the number of evacuees and minimizing the utilization of resources.

Problem Description and Model Formulation
The aim of our emergency evacuation planning problem is to decide the location of S shelters among the potential sites at the tactical level and the routing assignment of evacuees from stricken regions to safe destinations operationally so that the disaster areas are evacuated as quickly as possible. In this work, we focus on evacuation with private vehicles (i.e., car-based evacuation). Figure 1 illustrates the evacuation network. The objective of our problem is to minimize the total evacuation time on all road segments. Before introducing the mathematical model of the SLERA problem, we first detail how to model the travel time. It is generally known that there is a positive and monotonically increasing relationship between travel time and traffic flow, as a large traffic volume will lead to congestion and increase the travel time on roads. Getting a reasonable estimate of the entire evacuation time is the main reason for modeling traffic time on road segments in accordance with corresponding traffic flows. In the literature, there are different functional forms that model the relationship between travel time and traffic flow on a road segment. In particular, the U.S. Bureau of Public Roads (BPR) function has been extensively adopted in evacuation models to estimate the traffic times, for example, in [4,5]. Thus, the BPR function is applied in this work as well. In BPR function, the travel time t is not constant, as the level of congestion has an impact on the travel speed of vehicles, especially during disaster evacuation. According to [5], travel time t(x) is a monotonically increasing function of the traffic flow x. The relationship between travel time and traffic flow is formulated as where t(x) is the travel time when the traffic volume on the road is x, c denotes the maximum flow capacity, t 0 is the free-flow travel time at zero volume. The parameters γ ≥ 0 and β ≥ 0 are the tuning coefficients, which determine the shape of the function. The exact values of the coefficients are defined in accordance with the road characteristics, and they are generally taken as 0.15 and 1-5, respectively. When higher flows occur on the link, the coefficient β determines the threshold at which the BPR function rises significantly [35]. Considering that the evacuated tool is the same type of vehicle, a reasonable smaller value of β is preferred. So, we use β = 2 in this work.

Notation
The notation used throughout the paper is summarized as follows. Sets:

Variables:
y j : = 1, if a shelter is located/opened at alternative shelter site j ∈ J; = 0, otherwise.
x r ij (ω): the number of individuals that evacuate from demand node i ∈ I to shelter j ∈ J via route r ∈ R in scenario ω ∈ Ω. f l (ω): the traffic flow on road l ∈ L in scenario ω ∈ Ω.

Expectation-Based Model
Based on the above definition and travel time estimation Equation (1), firstly, the mathematical formulation of the expectation-based problem is written as follows: The objective Function (2) minimizes the expected total evacuation time spent by the evacuees in the network. Constraint (3) limits the number of open shelters to a pre-specified number S. Constraint (4) ensures that evacuees are not assigned to a non-open shelter. Constraint (5) guarantees that the capacity of each opened shelter is not exceeded. The minimum utilization rate of each opened shelter j is constrained by chance Constraint (6), which ensures that the utilization rate is not less than a threshold value θ with a specified high probability of 1 − . Constraint (7) ensures that every evacuee from every origin i is assigned to a shelter and a route leading to that shelter. The set of Constraint (8) computes the traffic flow on every road segment l ∈ L. Finally, the domains of the variables are restricted by Constraints (9)- (11).
In M1, the objective criterion is expectation-based (i.e., risk-neutral), and the objective function value entirely depends on the scenario set Ω, which captures the randomness under a finite sample distribution space assumption. The expectation-motivated objective function is reasonable and reliable only if the number of realizations of random demand is large enough. In practice, however, it is difficult to take a large number of scenarios into consideration, due to the restriction of solution techniques. Unavoidably, a case exists where the solutions obtained according to the expectation-based model do not perform well in some extreme demand cases, which does not occur in scenario set Ω. To determine the risk of solutions performing badly in exceptional cases, the introduction of a risk measure method to assess the expectation-based objective value is preferred.

Risk Measure by CVaR
In this subsection, we briefly introduce the concept of CVaR and present a linear programming formulation of CVaR under the assumption that the distribution space of random variables is finite. We recommend readers to refer to [36][37][38] for more details.
As a modeling method, it should be noted that the expectation-based objective Function (2) only accounts for random scenarios in a finite set Ω. The performance of solutions obtained from given scenarios will not be guaranteed under certain realizations of random parameters. To ensure an efficient emergency response operation in which optimal decisions present relatively good stability even in worst-scenarios, especially in the context of disaster management, authors have applied indicators such as the Mean-Variance [39] and CVaR [38] as risk measures. In particular, CVaR has gradually been adopted in pre-disaster relief network design problems [40], emergency logistics planning problems during disasters [2] and disaster location and allocation problems [41]. In these works, risk measure was adopted to avoid a case where the solutions provided by the risk-neutral model might perform poorly in extreme scenarios (e.g., a very high level of unsatisfied demand under certain realizations of random data).
We firstly introduce the concept of the mean-risk model in relation to general two-stage stochastic programming model considering risk aversion. Then, we present the risk measure definition and the risk-averse model, and finally, we conclude this subsection with a linear representation of the risk measure.
(Ω, F , P) denotes an abstract probability space, in which Ω is the sample space, F is a σ-algebra on Ω, and P is a probability measure on Ω. In the general risk-neutral two-stage stochastic model, (i.e., expectation-based), the sample space Ω is taken as a finite probability distribution where Ω = {ω 1 , ..., ω N }. Then, the common expression form of the expectation-based stochastic linear programming problem is written as min where f (x, ω) = c T x + Q(x, ξ(ω)) is the total cost function of the first-stage problem and Q (x, ξ(ω)) is the optimal value of the second-stage problem (also known as the recourse function), which is formulated as In the above-mentioned problem, "x" is the first-stage decision variable vector (also referred to as a here-and-now decision), such as the location decision variables in this work, which is made before the realization of the actual random parameter ξ(ω) for the elementary event ω ∈ Ω; "y" is the second-stage variable vector, like the routing assignment variables, that occurs after a realization of ξ(ω) has been identified. Indeed, this risk-neutral approach that optimizes the cost c T x of the first-stage decision plus the expected cost of the recourse function does not take the risk element into account. The mean-risk model is one of the classical approaches to incorporate the risk measure into a general two-stage stochastic programming model while simultaneously minimizing the risk function. The risk-neutral two-stage framework can be extended to the risk-averse mean-risk version formulated below: where λ is a non-negative trade-off coefficient representing the exchange rate between the mean value and risk, called the risk coefficient or risk level. Specially, the objective Function (14) is risk-neutral when For the random variable f (x, ω), CVaR α denotes the conditional value-at-risk of f (x, ω) at a given confidence level α. The mean-risk model minimizes both the expected cost E(·) and risk measure function CVaR α (·). The relevant definitions of the risk measure CVaR are presented as follows: Definition 1. Let F Z (·) represent the cumulative distribution function of a random variable Z. According to [36], the value-at-risk (VaR) at a given confidence level α, denoted by VaR α (Z), is defined as where α ∈ (0, 1].
The conditional value-at-risk at a given confidence level α is related to the expectation of Z under the condition that it exceeds the α-quantile threshold. Its definition is formulated as CVaR can capture a wide range of risk preferences by altering another risk coefficient α, in particular, the risk-neutral case (α = 0) and the pessimistic worst-case for a sufficiently large value of α → 1. Z is a random variable that obeys a discretion distribution with finite realizations {z 1 , ..., z N } and corresponding probabilities {p 1 , ..., p N }. An alternative definition of CVaR α (Z) , which gives a linearization description of CVaR, is restated below [36].

Definition 2.
The conditional value-at-risk of random variable Z at the confidence level α is given as where [Z − η] + := max{0, Z − η} and α ∈ (0, 1]. By means of linearizing the term ([Z − η] + ) in Equation (17), the conditional value-at-risk of random variable Z at the confidence level α can be rewritten as subject to

Risk-Averse Model
Given the definition of risk measure, a mathematical formulation of the stochastic version considering risk measures can be represented as follows. This is based on model M1 and Equations (14) and (18) subject to constraints (3)-(11) In this risk-averse model, smaller values of the risk measure part are preferred, as this indicates the corresponding risk value of the solution's performance loss in exceptional demand cases. λ and α are two important coefficients that are used to tune the risk preferences for decision makers or the uncertainty level of the decision environment. However, we cannot directly solve this risk-averse stochastic model, because the objection function contains a cubic item f β+1 l (ω), which is a nonlinear objective function. Heuristic techniques or the piecewise linear approximation method are frequently used to handle related nonlinear objective functions occurring in mixed-integer programming (MIP) models. Recently, motivated by the advances in second-order cone programming (SOCP), the nonlinear MIP models were reformulated as corresponding second-order conic mixed-integer programs and then directly solved by off-the-shelf solvers (e.g., CPLEX, Gurobi). In the following section, we describe the SOCP solution approach to problems with the nonlinear objective function.

Solution Approach
In this section, we first give a brief review of second-order cone programming and provide the mixed-integer SOCP formulation of our model. Then, a sample average approximation of the uncertain problem as well as the sample based mixed-integer linear programming (MILP) reformulation of the non-convex chance constraint are presented.

Second-Order Cone Programming Approach
Conic optimization refers to the optimization of a linear function over conic inequalities. With the large number of practical applications and efficient algorithms (e.g., interior algorithms), SOCP has become a state-of-the-art technique used in commercial solvers such as CPLEX. With the availability of commercial solvers, conic integer models have recently started to be used in airline scheduling problem applications [42] and emergency response planning for disasters [4,5].
By employing SOCP, the nonlinearity in the objective function is transferred to the constraint set in the form of second-order quadratic constraints. One of the common uses of conic quadratic inequalities is to represent a hyperbolic inequality, as in the form of Equation (25) (please refer to [43] for more details): Then, each hyperbolic inequality can easily be transformed into a second-order cone inequality, as follows: Here, || · || denotes the Euclidean norm. In our work, we define an auxiliary variable π l for each l ∈ L, and the expectation part of the objective function becomes The road characteristic coefficient β is generally adopted in a range [1,5], which is explained in Equation (1). In this work, β = 2 is adopted to represent f 3 l (ω) ≤ π l (ω) with hyperbolic inequalities of the form where m l , π l , µ l are auxiliary variables that are used to define hyperbolic inequalities. These hyperbolic inequalities are represented by their respective quadratic cone constraints: Despite the problem still being a nonlinear programming problem in which the nonlinearity element is transferred to the constraint set, it can be solved efficiently when represented as a SOCP. Note that SOCP problems are solved with the barrier algorithm in CPLEX. The resulting conic problem with SOCP constraints is given below: s.t. constraints (3)-(11)

Sample Average Approximate Algorithm and MILP Reformulation of Chance Constraint
The sample average approximate (SAA) algorithm is one of the predominantly used solution methods to tackle random characteristics in stochastic programming. The basic idea of this approach is to replace the original distribution of a random parameter d with a sample Ω of size |Ω| = N, that is, d(ω), ω ∈ {1, 2, ..., N}, where N is much smaller than the real number of scenarios of d. In SAA, Monte Carlo simulation is used to generate scenarios to represent the demand for evacuees. It is assumed that the demand in a disaster zone is a multivariate uniform distribution U(µ, ∑), where µ and ∑ are the mean and covariance matrix of the random parameter, respectively. |Ω| scenarios are generated after the Monte Carlo simulation has been run, and each scenario has the same probability 1 |Ω| . After the |Ω| scenarios are generated, the expected objective Function (34) is estimated by the sample average Function (48) On the other hand, in line with [44], the chance Constraint (6) can be estimated by an indicator function 1 |Ω| ∑ |Ω| n=1 I (0,∞) ∑ i∈I ∑ r∈R ij x r ij (ω) − θb j y j ≤ , which requires a certain percentage of the samples to satisfy the chance constraint, where I (0,∞) (·) is an indicator function denoted as The corresponding problem based on the sample average approximate algorithm (refer to as the SAAP) is formulated as follows: Although the chance constraint is converted to a form involving the indicator Function (50) in the SAAP, it is still difficult to be solved directly, as it is not convex. To solve the sample-based chance constraint directly, an equivalent mixed-integer linear programming formulation is deployed to substitute Equation (50) into the SAAP. For a given sample of size |Ω|, auxiliary binary decision variables z(ω), ω ∈ Ω are introduced. For each scenario ω, z(ω) = 0 if the chance constraint is feasible or z(ω) = 1, otherwise. Additionally, a big constant number "M" is used to model the case when the chance constraint is violated. Therefore, the function of the chance constraint is realized by limiting the sum of binary variable z(ω), ω ∈ Ω. Finally, the equivalent expression of the chance Constraint (50) is

Computational Results
In this section, we describe computational experiments that were carried out for two case studies. In the first part of the computational study, we investigate the comparison results with the risk-neutral model to provide insight into our proposed model considering the risk measure. In the second part, we show the effectiveness of the modeling approach and conduct a sensitivity analysis of the risk parameters. In this work, we considered the evacuation pattern as recommended for the movement of victims toward a safe place with the support of vehicles, which should be planned systematically. We coded the proposed model using MATLAB 2018b and solved it using CPLEX solver 12.8 on a DELL PC with a CPU capacity of 4.00 GHz and 32G RAM.

Instance Data
We solved the risk-averse model with different test networks. The test problems we used were from two sources. The first instance, called Sioux Falls, was inherited from [5], which comes from an online library called the Transportation Network Test Problems. Sioux Falls were downloaded from [45], and this instance was used to provide results comparisons and contrasts with the risk-neutral model in the literature [4]. We focused on the second case that was originally developed by [12] to determine the location and evacuation of a bushfire in Australia. In this case study, the region was represented by a network consisting of 11 nodes and 90 links. The detailed information related to instance cases are presented in Appendix A, and the associated codes can be downloaded from the following link: https://pan.baidu.com/s/1QrG2auqDfCbKepnh8EIWOQ.

Comparison of the Results between Risk-Averse and Risk-Neutral Models
Sioux Falls comprises 24 nodes and 76 arches. The details of the network are presented in Appendix A, Table A1. In line with [5], we set nine potential shelter sites (i.e., destination nodes) at nodes 2, 6-8, and 16-20, and the other nodes were sites affected by disaster (i.e., origin/demand nodes). So, there were 135 |O-D| pairs in total. For each |O-D| pair, there were many possible evacuation routes. It was difficult to solve this problem as all available routes were taken into account. Therefore, we chose K shortest paths for each |O-D| pair by the Deletion Algorithm [46], which is based on the classical Dijkstra Algorithm. We determined the demand for each origin node randomly between 200 and 600, and we assumed that the maximum traffic flow rate (capacity) of the road segment was 60 vehicles per hour with a free flow speed of 80 km/h in an uninterrupted traffic flow. We considered three types of facility capacity (1000, 2000, and 3000), and the detailed capacity of potential shelter sites was [3000 2000 1000 1000 3000 3000 2000 2000 1000].
In the risk aversion setting, the solution performance under some worst demand scenarios is emphasized so that the obtained solutions are more conservative than those obtained with the risk-neutral setting in the literature [4]. To illustrate this difference, we compared the proposed risk-averse model with the risk-neutral (i.e., expectation based) approach presented in [4]. The experiment setting was S ∈ {3, 4, 5, 6} and three types of risk coefficients were used to reflect different risk preference levels λ ∈ {0, 0.1, 0.5, 0.9}. Note that the risk-averse model is transferred into the risk-neutral one when λ = 0. Firstly, we took |Ω| = 10 as the scenario's sample size to generate the shelter location decisions for both the risk-averse and risk-neutral models. The corresponding results of the shelter usage rate is reported in Tables 2-4. Then, to measure the performance of the obtained solutions in a bigger sample space, we applied the Monte Carlo simulation to generate a bigger sample to test the performance of solutions over a wide range of situations. These simulation results are presented in Table 5. Intuitively, a sound solution not only performs better in a small sample space, but is also guaranteed under general circumstances. The utilization rate of opened shelters is summarized in Tables 2 and 3. The column "#shelter" lists the particular opened shelters, and the symbol "-" indicates that the corresponding utilization rate is 0 or the alternative shelter site is not chosen as opened shelter. As shown in Tables 2 and 3, we can observe that the shelter site with a bigger capacity is preferred when the risk coefficient λ increases such that it exceeds a threshold value. For example, in S = 5, as the risk parameter λ increases to 0.1, compared with the location decision of the risk-neutral model, shelter #2, which has a capacity of 3000, is more preferred over shelter #6, which has a capacity of 2000. Moreover, in the risk-neutral model, when λ takes a value of 0.5, shelter site #16, which has a capacity of 3000, is chosen as an open shelter to substitute shelter #17 and λ = 0.1. Turning to the objective values presented in Table 4, the total evacuation time of the risk-averse model is greater than that of the risk-neutral one, and it increases without decreasing as the risk coefficient λ increases. In addition, the time spent on evacuation decreases as the number of open shelters increases. However, based on the effect of S on the total evacuation time, what can be estimated is that establishing one more shelter adds little to the evacuation of the entire system when the available number of safety facilities is greater than 6. On the other hand, as shown in the right part of Table 4, the CVaR of total evacuation time decreases as λ increases. This means that the risk that the obtained solutions will perform poorly under other unexpected scenarios is relatively reduced. In light of the location decisions of λ = 0.5 and λ = 0.9 being identical, we carried out a simulation of cases of "risk-neutral" and "risk-averse" where λ ∈ {0.1, 0.5}. Table 5 presents the simulation results for the total evacuation time under sample sizes of |Ω| ∈ {50, 100, 200}. From Table 5, we can observe that the risk-averse model performed slightly better than the risk-neutral one, as the total evacuation time generated by the simulation of the risk-averse approach was less than that of the risk-neutral approach. This rightly correlates with the results of CVaR for the total evacuation time (TET) shown in Table 4, whereby the solution quality is guaranteed in more general scenarios.

Effectiveness and Sensitivity Analysis
In this section, we reference the instance of a bushfire event in Australia that was presented in [12] using their network structure. The Shire of Murrindindi, which experienced a series of bushfires on 7 February 2009, Black Saturday, is located approximately 100 km northeast of Melbourne in Victoria, Australia. It covers an area of 3889 km 2 and includes the towns of Alexandra, Buxton, Eildon, Flowerdale, Kinglake, Marysville, Molesworth, Strath Creek, Taggerty, Yarck, and Yea. The study area consists of six main townships which suffer from fire risk: Narbethong i 1 , Marysville i 2 , Taggerty i 3 , Buxton i 4 , Cambarville i 5 , and Rubicon i 6 . In [12], the demand of evacuees was regarded as deterministic, as given in Table 6. Conversely, we assumed that the actual demands in terms of disaster zones were random and independent. We used the deterministic demand reported in Table 6 as the expectation d of uncertain demand. Table 6. The case study demand and capacity data.

Index i Evacuation Nodes Demand (d i ) Index j Shelters Capacity (b j )
1 Narbethong  240  1  Yea  1500  2  Marysville  260  2  Alexandra  500  3  Taggerty  170  3  Thornton  500  4  Buxton  130  4  Eildon  1000  5  Cambarville  110  5  Yarra Glen  1000  6 Rubicon 190 Five shelters in adjacent townships (destinations) were nominated by Country Fire Authority (CFA): Goughs Bay Fire Station in Alexandra j 1 , a recreation reserve oval in Thornton j 2 , a basketball court in Eildon j 3 , a skate park in Yea j 4 , and a racecourse track in Yarra Glen j 5 . Figure A2 in Appendix A gives the location distribution of these disaster regions and safe destinations. Each shelter has a maximum capacity for evacuees; these capacity values are reported in Table 6. The details of the road segments and road network information are provided in Appendix A, Tables A1 and A2. We set the capacity for each road segment randomly with a maximum traffic flow rate of 30 vehicles per minute and a free flow speed of 100 km/h. To balance the utilization rates among opened shelters, we set θ = 0.2, and the reliability level was = 0.1. Note that despite a much smaller value of than 0.1 being desirable, we chose = 0.1 as a result of the limited solution ability in this work. Specifically, the reliable level depends on the sample size of |Ω| in the sample-based chance Constraint (50), so the larger |Ω| is, the smaller can be, for example, = 0.05 was feasible since the sample size was |Ω| > 20.

Effectiveness Analysis
To illustrate the validity of our proposed model, it is necessary to compare the output of computation results with real evacuation outcome. However, it is difficult to obtain the real evacuation schemes used in the bushfire disaster event, as the features of evacuation activity were real-time and spontaneous. In general emergency evacuation events, the nearest allocation (NA) fashion is mostly used to assign shelters and corresponding paths to evacuees, as the information on path lengths or free-flow travel time is more accessible to evacuees compared with actual travel times in a disaster [4]. On the other hand, in an emergency or chaotic environment, it is easy for evacuation personnel to self-select the shortest route based on the distance/travel time to reach a safe destination as quick as possible. Therefore, we used the calculation results generated by the NA approach to approximate what happened in real evacuation case and compared the computational results between our model and the NA method to test the validity of our model.
We assumed that demand distribution of origin nodes was uniform with the means given in Table 6, and the distribution interval was [(1 − e)d, (1 + e)d] where 0 < e < 1. Taking e ∈ {0.1, 0.3, 0.5} to capture different demand fluctuations, we randomly generated 4 scenarios for each distribution pattern with different numbers of opened shelters S ∈ {2, 3, 4, 5}. The set of total evacuation times is summarized in Table 7. We can make several observations from the table: In order to gain insight into the performance difference, Figures 2 and 3 presents the traffic flow distribution in the evacuation process. Compared with Figures 2 and 3, we can see that the traffic flow in the evacuation paths generated by our model is evenly distributed, and the upper bound of flow in each road segment is less than 250. However, in the traffic flow distribution of NA, road segments 2, 6, 10, 12, and 17 are the main converging roads to safe destinations and road segments 27-42 are idle. Certainly, compared with Figure 2, the probability and severity of traffic congestion in the real evacuation case on road segments 2, 6, 10, 12, and 17 is large, as the traffic flow exceeds the road capacity. The unbalanced distribution of traffic flow directly results in poor system efficiency during emergency evacuations.

The Impact of Risk Coefficients on CVaR and the Total Evacuation Time
In the risk-averse shelter location and routing assignment problem, the parameter λ and the confidence level α are two key coefficients related to shelter location decisions and emergency evacuation planning. For example, a higher λ indicates that more emphasis is placed on the risk part and the corresponding decisions tend to be more conservative to avoid cases of extreme demand. In this part of the computational study, we analyzed the results based on the bushfire evacuation problem with |Ω| = 10 and the risk parameters λ ∈ {0.1, 0.  Tables 8 and 9 report the impacts of risk parameters on the expected total evacuation time and CVaR. Figures 4 and 5 give graphic illustrations of the impacts of λ on the trade-off between the expected total evacuation time (TET) and the associated CVaR under different available shelter cases. Figures 6 and 7 show the impacts of the confidence level α on the risk and objective function value. Figure 8 represents an efficient frontier of risk and expectation.
In the risk-averse model, higher values of λ and α indicate a more conservative decision preference and suggest that more attention should be paid to some extreme realizations of uncertain demand to avoid loss of performance in the solutions. Thus, increasing risk parameter λ brings about more weight to the risk measure value (i.e., CVaR) of the total evacuation time at the expense of increasing the expected TET. As shown in Figures 4 and 5 (as well as Tables 8 and 9), the CVaR of TET decreases and gradually converges to a fixed value as λ increases, while the expectation of TET increases. These increments, shown in Figure 5 and Table 8, become more pronounced when the number of opened shelters S is small, such as S ≤ 3. In addition, we can also observe that the graphs of the expected TET and corresponding CVaR of larger S values are always below those of the smaller S values, which implies that improvements in CVaR, and thus the level of compromise from the expectation, increase with S. Besides, it is obvious that increasing the number of opened shelters will not promote a significant decrease in the total evacuation time once S is greater than 3. In contrast, it is resource-wasteful.
Turning to Figure 6, we can observe that a higher α value leads to a greater CVaR, which illustrates that CVaR is a non-decreasing function of α. Similarly, the same relationship is valid for the objective function value graphs in Figure 7. This consequence can be explained by the definition of CVaR in which a smaller confidence level α indicates a worse realization. According to both figures, we can say that the influence of the risk coefficient λ becomes more significant as larger α values are given. Figure 6 and Figure 7 further confirm that the contribution to improving the emergency evacuation effect of simply increasing the number of shelters is small when S exceeds 4.
Furthermore, to facilitate the decision makers who are seeking to reach the right balance between expectation of TET and accompanying risk, Figure 8 lists solutions to the expected TET and corresponding risk with different S values when α = 0.95. As shown in Figure 8, along with the left tails of curves (starting from smaller λ values), the CVaR of TET decreases significantly at the cost of increasing the expected TET. On the other hand, it is obvious that the right tails of these curves tend to be flat. This means that higher compromises from the expected TET are required, especially in the case where the risk coefficient λ is close to 1, to further decrease the risk.

Conclusions
In this work, we studied a shelter location and routing assignment problem under random demand circumstances and introduced a risk-averse stochastic programming optimization model. Our aim was to minimize the total evacuation time using a system optimal style. To handle the non-linearity in the objective value arising from using a traffic assignment model, second-order cone programming was deployed to transfer the non-linearity from the objective function to constraints. In addition, we performed a series of numerical analyses to validate our proposed model and the impacts of key features were discussed for decision making purposes. This emergency management strategy involving risk measure is valuable as it could support plan making on both tactical and operational levels.