A robust stochastic Casualty Collection Points location problem

In this paper, a Casualty Collection Points (CCPs) location problem is formulated as a two-stage robust stochastic optimization model in an uncertain environment. In this modelling approach, the network design decisions are integrated with the multi-period response operational decisions where the number of casualties with different levels of injuries coming from the affected areas is uncertain. Furthermore, the transportation capacity for the evacuation of casualties to CCPs and hospitals is also uncertain. To solve this complex problem, a robust sample average approximation method with the feasibility restoration technique is proposed, and its eﬃciency is examined through a statistical validation procedure. We then evaluate the proposed methodology in the backdrop of a hypothetical case of Bhopal gas tragedy (with the same hazard propagation proﬁle) at the present day. We also report the solution robustness and model robustness of 144 instances of the case-study to show the proﬁciency of our proposed solution approach. Results analysis reveals that our modelling approach enables the decision makers to design a humanitarian logistic network in which not only the proximity and accessibility to CCPs are improved, but also the number of lives lost is decreased. Moreover, it is shown that the proposed robust stochas- tic optimization approach converges rapidly and more eﬃciently. We hope that our methodology will encourage urban city planners to pre-identify CCP locations, and, in the event of a disaster, help them decide on the subset of these CCPs that could be rapidly mobilised for disaster response.


Introduction
Severe weather events and natural disasters have displaced approximately 32 million people globally in 2012 and numbers are projected to continue rising ( IPCC, 2014 ). According to the Centre for Research on the Epidemiology of Disasters , over the past ten years, natural disasters affected almost 1.7 billion people, including 0.7 million killed, and resulted in 1.4 trillion dollars in damages worldwide ( Guha-Sapir, Hoyois & Below, 2015 ). Similarly, manmade disasters have human, environmental and economic consequences. Examples of such disasters include stampede, nuclear or chemical plant explosion, emergencies resulting from incorrect handling/transportation of hazardous materials, water contamination and oil spill. Man-made disasters happen mainly due to accidents, negligence or incompetence. With the global increase in the number and severity of the disasters, researchers from differ-ent disciplines are increasingly paying attention to disaster management problems.
Alerts and early warning systems are among the tools available to city planners for dealing with emergencies. These inform the population of an impending disaster, e.g., tsunami warning system of the Japanese Meteorological Agency ( Tatehata, 1997 ) and COBRA alerts in the UK ( Thunhurst, Ritchie, Friend & Booker, 1992 ). Although these are useful for the advance warning, it is also essential to have, in place, existing strategies for humanitarian logistic network design that could be initiated after a disaster occurs. Analytical models may be developed to represent population centres with critical infrastructures like hospitals, power plants and transport networks. This will enable experimentation of humanitarian logistic operations and inform city planners if these are fit for purpose and where improvements can be made. In this paper, we have developed one such model. The model is motivated by a disaster which took place in Bhopal, India, in the year 1984. It is commonly referred to as the Bhopal gas tragedy and was caused due to a leak of toxic gas (methyl isocyanate) from a pesticide manufacturing plant. In OR literature, the case study of Bhopal disaster has been used once before to illustrate a methodology that can help identify root causes of disasters and facilitating allocation of resources to prevent their occurrence. In their work, Ishizaka and Labib (2014) propose a hybrid method consisting of problem structuring, visualisation, Analytic Hierarchy Process and mathematical programming, with the objective to calculate the optimal allocation of available funds in order to avoid a disaster.
For our model, we use the backdrop of the propagation of hazard that took place on the night of 2-3rd December 1984; we use population and other model-specific parameters from the latest available census data and other municipal reports for the city of Bhopal. We consider a hypothetical case of a gas leak taking place in Bhopal in today's date which follows the hazard propagation profile (e.g., wind direction) reported back in 1984. The number of people dead as a direct consequence of inhaling toxic gas is estimated to be between 3700 and 16,000. Considering the catastrophic loss of lives, our objective here is to design a humanitarian logistic network in which response planning and operations are taken into account for the evacuation of the entire population of the affected areas to facilities that provide temporary medical assessment and treatment (these are referred to as Casualty Collection Points or CCPs) and to the hospitals. In our model there are two uncertain parameters, namely, the number of casualties and the transportation capacity. The motivation for using these variables is based on the hazard profile that was associated with the Bhopal disaster. The direction of the wind determined the number of people that inhaled the toxic gas. If the wind movement was in the direction of build-up population centres (called as wards) then this would affect more people. Furthermore, the demographics associated with a ward could have a bearing on the severity associated with inhaling the gas. For example, inhalation of the gas had different sensitivities associated with children and elderly people compared to the rest of the population ( Bowonder, 1987 ). Our model, therefore, considered this uncertainty in the number and severity of casualties. The motivation for the second uncertain parameter (transport capacity) is based on the generally accepted fact that developing countries often have inadequate transportation and which is likely to affect emergency evacuation ( Bisarya & Puri, 2005;Bowonder, 1987 ). Bisarya and Puri (2005) recommend that the people living in the vicinity of hazardous plants should be made aware of the sources of transportation/ambulances for emergency evacuation. However, in a disaster of such magnitude, it is important to consider not only public transport but also private vehicles for the transportation of casualties (as happened in Bhopal). Ownership of private vehicles will usually depend on the socio-economic status of the people living in different wards. Further, public transport capacity will also be dictated by transport infrastructure available in different population centres. In order to account for these variations, our model includes transport capacity as an uncertain parameter.
In such uncertain environment, decision makers are to act without exact or complete information about number of casualties from the affected areas and the transportation capacity for moving casualties to CCPs and hospitals. These factors cannot be confidentially estimated due to the unpredictability of time, place and severeness of a disaster as well as the changing roadway infrastructure as a result of disaster impacts ( Bayram & Yaman, 2018 ). In the context considered here, the number of casualties with different levels of injuries coming from the affected areas over the planning horizon and the transportation capacity for moving casualties are uncertain parameters. The uncertainty about future realizations of these parameters are considered in the form of random sample of scenarios incorporated in the problem formulation.
The vast majority of studies in disaster and emergency management have focussed on the distribution of relief in the aftermath of disasters ( Anaya-Arenas, Renaud & Ruiz, 2014;Paul & Zhang, 2019 ). In this context, stock location, resource allocation, and commodity flow from predefined warehouse locations to affected areas have been the most impactful variables to optimize for the construction of relief distribution networks. Casualty management problems, such as the one presented in this paper, can similarly be construed in terms of CCP location, casualty medical treatments, and casualty flow from the affected areas to safer places and hospitals. In spite of the importance of casualty management in humanitarian logistics, relatively little attention has been paid to this subject ( Gupta, Starr, Farahani & Matinrad, 2016 ). Our work is, therefore, a contribution to this literature; specifically, we are concerned with the casualty management functions of disaster management that are caused by human error, such as industrial accidents, and which are implemented after a disaster strikes (response phase of disaster planning).
A disaster may result in numerous lives being lost. However, the severity of potential threats in the aftermath of a disaster can be mitigated by providing fast and essential aids through intermediary sites. As mentioned earlier, these sites with short-term missions and temporary locations are referred to as CCPs. An overall view of CCP establishment and operations is presented in Fig. 1 . In existing literature, several terms have been interchangeably used to denote these facilities, such as field treatment site ( Drezner, 2004 ) or alternative care facilities ( Caunhye, Li & Nie, 2015 ). However, for consistency, we have used CCP for Casualty Collection Point or facilities that are functionally similar to CCPs. CCP locations are identified before the disaster occurrence, i.e. during the preparedness phase, but selected after the disaster has occurred, i.e. in the response phase (see Fig. 1 ). After choosing the right location and establishing 1 the CCPs, the following operational and tactical decisions are to be made in the response phase at CCPs: (i) triage, (ii) casualty registration, (iii) casualty medical treatment, (iv) casualty evacuation, and finally (v) shutting down the site(s) ( Koehler, Foley & Jones, 1992 ).
Uncertainty affects strategic CCP location decisions, and which have a bearing on tactical and operational decisions. The network design decisions are strategic decisions that are made when forecasting uncertain parameters. Planning and operational decisions, on the other hand, are usually made when parameters are more obvious (e.g., the parts of the city that may be affected due to an unfolding weather-related event). It is arguable that including strategic decisions would improve the quality of casualty management and other operational decisions. In particular, optimizing the location and allocation decisions at the strategic level with the hierarchical integration of the periodical policy decisions lead us to a two-stage stochastic optimization model. With this motivation, we reflect on a robust stochastic optimization approach, which simultaneously optimizes the number of CCPs, location, allocation, and capacity decisions at the strategic level and scenario-based casualty triage, casualty registration, casualty holding, and casualty transportation decisions in a multi-period planning setting, while satisfying the system constraints enhancing the problem objective function.
The organization of the remainder of this paper is as follows. In Section 2 , we provide a literature review and highlight the main contributions of this paper. Section 3 represents a generic robust optimization modelling approach and a two-stage formulation for the problem context presented in the paper. Section 4 contains a robust stochastic optimization procedure as well as the validation procedure. In Section 5 , we study the application of the model to the case study; we provide experimental results for extensive realistic problem instances; we discuss these results and performance of the solution methodology. Section 6 is the concluding section and discusses future work.

Literature review
This section presents a brief overview of research on casualty management and disaster response. ( Drezner, 2004 ) first introduced the CCP location problem in a discrete network and its application in disaster management in Orange County, California. Then, Drezner, Drezner and Salhi (2006) developed the problem to a multi-objective programming model to find appropriate locations for CCPs. Casualty transportation in cases of expected disasters and post-disaster, have been widely studied in the form of a transportation network design problem ( Ben-Tal, Do Chung, Mandala & Yao, 2011;Ozdamar, 2011;Shen, Pannala, Rai & Tsoi, 2008;Yao, Mandala & Do Chung, 2009 ). In this regard, An, Cui, Li and Ouyang (2013) and Kulshrestha, Lou and Yin (2014) developed a stochastic model that incorporates mass-transit casualty evacuation planning from pick-up locations. Goerigk and Grün (2014) , Najafi, Eshghi and Dullaert (2013) and Goerigk, Deghdak and T'Kindt (2015 ) studied the impact of multiple transportation modes including private vehicles, rapid transit, and mass-transit shuttle buses. Sacco et al. (2007) , Wilson, Hawe, Coates and Crouch (2013 ) and Kilic, Dincer and Gokce (2014 ) focused on processing operational decisions involving triage, transportation, and treatment for medical injuries over the planning period. He and Peeta (2014 ) and He, Zheng and Peeta (2015) underlined the impact of dynamic resource allocation on casualty transportation and evacuation.
In logistic network design, there exists temporal hierarchical structure between initial design considerations and the subsequent planning and operational decisions; this implies that these decisions are made under uncertainty ( Klibi, Lasalle, Martel & Ichoua, 2010;Shapiro, 2008 ). Klibi and Martel (2013 ) emphasized that in-dividual optimization of the logistical decisions may not guarantee an optimal solution for the whole operation. Amiri-Aref, Klibi and Babai (2018) showed that the integration of the design and planning decisions could improve the quality of solutions in network design when demand is uncertain. Due to unpredictability concerning the magnitude of a disaster, number and location of casualties, the availability of infrastructure, weather conditions, etc., providing a logistical response encounters a high level of difficulty and uncertainty ( Apte, 2010 ). Thus, for the construction of an optimization model, to enable the integration of design and planning decisions it is important to consider temporal hierarchical structures with uncertainty. Bayram, Tansel and Yaman, (2015) emphasised that disaster management models that do not take into account the uncertainties may lead towards inefficient logistical planning and operational decision making. According to Gupta et al. (2016) , who present the latest survey in this field, integrating decisions related to locating casualties and moving them to hospitals (or safer places) can save numerous lives and further research is required in this area.
In the existing literature, only a few authors have addressed the stochasticity in an integrated CCP network design problem with multi-period planning settings. Li, Nozick, Xu and Davidson (2012) developed a scenario-based bi-level programming model for the shelter location model with the evacuation consideration for a realistic case study of North Carolina and highlighted the impact of transportation when selecting the location decisions. Bayram and Yaman (2015) proposed a scenario-based two-stage stochastic shelter location model considering casualties (evacuees) allocation to the nearest facility to minimize the expected total evacuation time. Bayram and Yaman, (2017) provided the exact solution based on Benders-decomposition algorithm to the model formulated by Bayram and Yaman (2015) . They showed the importance of the inclusion of uncertainty in planning for evacuations. Despite the contribution of the abovementioned effort s on the interdependency of casualty transportation and shelter (CCP) location decisions in humanitarian logistics network design, the main shortcoming is the neglect of temporal hierarchy relationship between the strategic and planning decisions and the dynamicity of casualty arrivals as illustrated in Fig. 1 . Strategic decisions are adopted at the beginning of response phase in an uncertain environment where exact or complete information about the number of casualties are not available. Then, scenario-based multi-period decisions are made during the response phase in which we assume the horizon is composed of a set of discrete operational cycles. Note that a user makes periodical decisions on a timely basis (e.g., hourly, 8hourly, 12-hourly and daily). In fact, the consideration of casualty state transition from one operational cycle to the next and of hierarchical setting between decisions results to a multi-period twostage stochastic program model for the humanitarian logistics network design problem.
To the best of our knowledge, a two-stage stochastic modelling for the CCP location problem with uncertain number of casualties with different levels of injuries under a multi-period settings and uncertain transportation capacity is lacking in the literature. Given that this problem has the same NP-hardness property as a basic facility location problem, we have developed a heuristic robust method for solving the problem. This paper extends the literature related to the humanitarian logistic network design in the following three ways.
While several humanitarian logistical problems studied the response network design for providing medical supply from prepositioned warehouses or staging areas to the affected people through the points of distribution (POD), this paper focuses on a network design with casualty response planning from the affected areas to the evacuation points (EP) or safer places through the temporary CCPs. In the former context, relief items and supplies move towards affected areas, whereas the model presented in this paper relies on the flow of victims with life threating conditions from the affected areas to the EPs or hospitals through the intermediate CCPs in an uncertain environment. The additional contribution of this paper to the relevant articles in the literature (e.g., Apte, Heidtke & Salmerón, 2014;Yi & Ozdamar, 2007 ) is the explicit inclusion of the uncertainty inherent in CCP location-allocation decisions made at the design stage of the optimization model. In fact, the uncertainty is due to the time lag between the strategic design decisions in the first stage and the dynamic operational decisions in the second stage during the response phase. Strategic decisions on the number, location and allocation of CCPs are made through anticipating the plausible scenarios for the operational decisions in the second stage. Although several studies in the location-evacuation literature investigated the humanitarian logistic network design, they almost considered deterministic or meanvalue information. In this paper, we develop a two-stage stochastic programming modelling approach to cope adequately with the uncertainty inherent in disaster contexts, where the value of stochastic information is high. It has been shown that the inclusion of uncertainty at the strategic level improves the quality of the CCP design decisions ( Birge & Louveaux, 2011 ).
Second, the main aim of the problem considered here is to optimize CCP design decisions in view of the existence of the temporal hierarchy structure between the strategic and operational decisions over the planning period. The time setting between these decisions as well as the distinct time-horizon granularity are incorporated in the proposed model to capture the dynamic nature of lifesaving operations in the response phase. In this research, we deal with an integrated humanitarian logistic network problem in which strategic decisions are made in the first-stage model and operational decisions with anticipation of uncertain factors are made/revised during the multi-period planning horizon. This problem must not be confounded with problems which in fact optimizes the location and evacuation decisions simultaneously for achieving coordination, as pointed out in Sheu and Pan (2014); Yi and Ozdamar (2007) . It is important to note that in this modelling approach the objective is to use the anticipated decisions optimized for each operational cycle under all scenarios, so that more efficient and robust CCP design solutions are generated at the strategic level. From the practical point of view, strategic decisions include the number and location of CCPs to be opened, CCPs capacity allocation, allocation of affected areas to established CCPs, hospitals allocation to established CCPs and alternative CCP locations. These decisions, also known as design decisions, are made immediately following a disaster. Planning decisions such as casualty triage, casualty registration, casualty medical treatment and casualty transporting, then need to be made over the whole of the planning period. The number of casualties (with several levels of injuries) and the available transportation capacity are uncertain throughout the proposed network. Due to the hierarchical structure of strategic and planning decisions, finding an optimal solution for one activity is not usually sufficient for the whole of the response phase. Therefore, the focus of this paper is to present a model that reflects the hierarchical structure of the strategic and planning decisions in the presence of uncertain parameters in one unique problem in order to provide effective design solutions; Further, to formulate such a problem in the form of a two-stage robust stochastic optimization model .
Third, we proposed a robust stochastic optimization solution approach to cope with the infeasibility issues which may occur in the stochastic optimization problems. Our proposed approach returns robust solutions which are close to any given scenarios with minimum dispersion from the optimal values. This has been validated by the optimality gap analysis.
Our review of the literature on stochastic programming approaches specific to casualty management problems has shown that no existing model has taken into consideration the three features presented above. The main purpose of this study is to provide a specific representation of an integrated casualty management structure in an uncertain environment while the system constraints are met. To the best of our knowledge, the modelling of CCP logistic network design problem with the characteristics mentioned above has not been studied in the literature so far.

Modelling approach and problem formulation
In this section, we first present a generic robust stochastic optimization modelling approach and then apply this approach in the proposed CCP logistic network design problem where the uncertain number of casualties with different levels of injuries, and uncertain casualty transportation capacity, are described by a set of realizations or scenarios for their values.

Modelling approach
This problem is characterized by two decision variable sets: design variables and control variables. Let us assume x ∈ R n + 1 denotes the vector of design variables which need to be made hereand-now in the first-stage of decision-making problems with ndimensional integer space. The design variables have static nature during the planning horizon and are non-adjustable to the uncertain parameters. Let us further assume y ∈ R n + 2 denotes the vector of control variables in an n -dimensional nonnegative space that are subjected to adjustment once the actual data of the uncertain parameters reveals itself. This decision set is scenario-dependent and adjustable to the optimal value of the design variables. The control variables which are made in the second-stage of dynamic decisionmaking problems are so-called wait-and-see decisions. Considering the definition of the design and control variables, a general framework of a two-stage stochastic programming model with uncertain parameters is presented in the following, where C := { c, A , b } is a set of vectors of fixed coefficients of the first-stage decision-making problem and of free of noise in input data. Objective function (1) represents objective function of the first-stage decision-making problem and the expected optimal value of the second-stage decision-making problem, defined by Q ( x , ψ ) , and Eq.
(2) denotes the structural constraints with fixed parameters.
where, ψ := { q, B , C , e } defines a set of uncertain parameters, subjected to noisy input data, associated with the second-stage decision-making problem. Objective function (3) optimizes the control variables in the second-stage decision-making problem subject to noisy parameters. Eq. (4) denotes control constraints with uncertain parameters adjusted to the first stage variables. Let us denote the set of ψ as a random vector with corresponding probability space support and its particular realization. Suppose the expected value function E [ Q ( x , ψ ) ] with random vector ψ has finite support . That is to say ψ has a finite number of realizations or scenarios ψ (ω) Therefore, the expected value function is represented as follows: where, for each ω ∈ = { 1 , 2 , ..., | | } , Q ( x , ψ ( ω) ) denotes the optimal value of the deterministic-equivalent linear formulation of the second-stage decision making problem: If the set of constraints (7) has no feasible solution, the secondstage decision making problem is infeasible. Under this condition, there exists at least one scenario realization ω ∈ , for which q T (ω) y (ω) = + ∞ and so Q ( x , ψ ( ω) ) = + ∞ . On the other hand, this problem could be unbounded depending on the first-stage variables and scenario realizations and hence Q ( x , ψ ( ω) ) = −∞ .
The classical stochastic programming is likely to be infeasible especially when the distribution of uncertain parameter is unknown, or the uncertain parameter realizations do not follow a specific distribution ( Khor, Elkamel, Ponnambalam & Douglas, 2008 ). Due to the lack of information or imperfect data in disaster management, such as, location and time of disaster, it's severity in terms of number of casualties, and available transportation capacity subsequent to the disaster, the parameters are almost unpredictable or are forecasted with a wide range of variability. This, coupled with the need to execute a large number of scenarios, will most likely produce infeasible solutions to the stochastic programming ( Neyshabouri & Berg, 2017 ). To tackle possible infeasibility due to the presence of uncertain parameters and the risk attributed to the decision-maker, robust counterparts problem is proposed. Its purpose is to find an optimal solution that satisfies all constraints for any uncertainty realization while reducing the risk of dispersion of the objective function value. In the robust optimization literature, two performance metrics that have been widely applied are the concept of solution robustness and model robustness ( Mulvey, Vanderbei & Zenios, 1995 ). Our robust counterpart problem studies both solution robustness and the model robustness concepts simultaneously. Since the robust counterpart problem accounts for the second-stage decision-making problem with uncertain parameter realization, without loss of generality, we mainly focus on the formulations given in (6) and (7).
To achieve solution robustness, Mulvey and Ruszczy ński (1995 ) measured the dispersion of the objective values by minimizing the average of standard deviation (or absolute deviation) of the objective values over all scenarios. This metric guarantees that the second-stage solutions are close to any scenario realizations applied in the problem ( Ben-Tal, Goryashko, Guslitzer & Nemirovski, 2004 ). In order to avoid nonlinearity resulting from the standard deviation formulation, we instead utilize the absolute deviation one, denoted by (ω) in (8) for each scenario ω ∈ .
The model robustness, which focuses on infeasibility issue as a result of a violation of data-driven parameters, takes into account infeasibility penalty ρ in the objective function of the secondstage decision-making problem. In this modelling framework, the constraint violation is measured by an infeasibility variable vector z (ω) , where a positive value of z (ω) show the amount of infeasibility of the corresponding scenario ω ∈ in the model. It is clear that z (ω) = 0 if the model is feasible. The mean value of probable infeasibilities is then penalized in the objective function.
Considering both solution robustness and model robustness represented in (8) and (9) , respectively, the robust counterpart of the second-stage stochastic programming is given as follows.
s.t. constraints (8)-(9) , where the expected cost function is presented in the first term of (10) and solution robustness and model robustness are given in the second and third terms of (10) , respectively. Since the terms stated in (10) need to be unified, we use coefficients β 1 and β 2 to provide a compromised objective function, as denoted in (11) .
In the next subsection, an extended formulation of the robust two-stage stochastic programming model to design the casualty collection logistical network problem is presented.

Problem formulation
The context of the study is based on the 1984 Bhopal gas tragedy and our methodology for problem formulation is inspired by the guidelines provided in the technical report of Haynes and Freeman (1989 ). One of the key recommendations of this report is the importance of designing an efficient logistical network in cases of disasters with mass causality. With this motivation, a robust two-stage stochastic programming model is formulated to de-velop a logistics network design problem for the Casualty Collection Points in the event of a disaster.
The first-stage objective function follows the recommendations of the report ( Haynes & Freeman, 1989 ), where the CCP locations should be close enough to both affected areas and hospitals so as to facilitate the efficient movement of casualties (here, the casualties are people who have been affected by the disaster). This objective function considers a fixed cost that can be assigned to each potential location for establishing a CCP and relative cost associated with distance to travel from the affected areas to the established CCPs and/or to the hospitals. According to the technical report, when injury severity is minor and the hospitals are in the vicinity of the affected areas, casualties can travel directly to established CCPs or hospitals without assistance (so-called self-evacuees ). On the other hand, casualties with the need for intermediate or immediate medical care are directed to established CCPs by either emergency, mass-transit or even private vehicles (so-called emergency-evacuees ). The emergency-evacuees go through four stages-triage, registration, treatment and evacuation. Incoming casualties to the CCPs are diagnosed for severity of their injuries (triage) and are registered subsequently. Temporary hospitalization and first aid medical services are then provided to the casualties (treatment). They are then transferred to the hospitals or other health and care facilities for further treatment (evacuation). The distance to travel for both self-evacuees (moving from affected areas to the established CCPs or hospitals) and emergency-evacuees (first travelling from affected areas to the CCPs, and then from the CCPs to the hospitals) is formulated in the form of a travelled distance cost minimization function . The second-stage objective function minimizes the expected operating costs and the penalty cost due to lives lost. The expected operating cost consists of the cost incurred by periodical decisions for casualty holding and transportation as well as the penalty cost due to system inefficiency and lack of adequate resources. In this study, the first-stage objective function incorporates the design decisions, while the second-stage objective considers the planning decisions incorporated into the design decisions. A typical CCP logistical network is graphically illustrated in Fig. 2 .

The first stage model
The first-stage of casualty collection logistical network design problem focuses only on the humanitarian objective optimization with the aim of locating the CCPs where the cost of travelled distance to affected areas and hospitals are minimal. Let us introduce the set of affected areas by I = { 1 , 2 , ..., |I| } , where the casualties come from. This set can also be viewed as demand points in the business context. The set of potential locations for establishing CCPs are denoted by J = { 1 , 2 , ..., |J | } and the set of existing hospitals to serve CCPs are indicated by K = { 1 , 2 , ..., |K| } .
The set of first-stage design decisions is composed of (1) determining the number of required CCPs to meet demand, (2) selecting the location of CCPs among potential locations where each CCP is characterized by its capacity, and (3) allocating the affected areas as well as hospitals to every established CCP. The input parameters, according to the first-stage requirements, contain the distance matrices in our designated network. Let us denote the distance from affected area i ∈ I to potential CCP location j ∈ J by D IJ = [ d i j ] |I|×|J | , the distance from potential CCP location j ∈ J to hospital k ∈ K by D J K = [ d jk ] |J |×|K| , and the distance from affected area i ∈ I to hospital k ∈ K by D IK = [ d ik ] |I|×|K| . The binary decision variables used in the first stage model are also represented in the following: X j = 1 if potential location j is selected as a CCP, and 0 otherwise, Y i j = 1 if affected area i is allocated to potential location j, and 0 otherwise, V jk = 1 if operating CCP j is allocated to hospital k , and 0 otherwise, U ik = 1 if affected area i is allocated to hospital k , and 0 otherwise.
These strategic decisions are made considering the uncertain parameters for the whole planning period. The uncertain parameter that is used most often in network design is demand value, which corresponds to the flow of casualties in the humanitarian context. In addition to this, we also incorporated uncertain transportation capacity into the model to achieve more realistic and more reliable results. The set of all possible flow of casualty scenarios and available transportation capacity scenarios are denoted by ϒ and , respectively. At the second stage, while realizing the possible scenarios υ ∈ ϒ and γ ∈ , the response decisions, including (i) triage, (ii) registration, (iii) treatment, and (iv) evacuation, are adopted over the planning period. Let Q ( x , ω ) be the solution of planning and operating decisions at the second-stage depending on the scenario ω = ( υ, γ ) ∈ = ϒ × , where represents a set of all combinations of scenarios υ ∈ ϒ and γ ∈ . Let assume π (ω) is the probability of each scenario occurrence, where π (ω) = π (υ ) .π (γ ) . Thus, we can introduce π (ω) q ( x , ω ) as the expected value of the objective function of the second stage, where α represents the provisional cost of travelling per distance unit. The objective function of the first stage model, presented in (12) , contains the expected cost of the second-stage problem after uncertainty realization and the cost related to the travelled distance. Constraints (13) ensure that each affected area can be allocated to each established CCP. Constraints (14) show that only the established CCPs are allowed to be allocated to the existing hospitals. Constraints (15) represent that each affected area must be allocated to either established CCPs or existing hospitals. Constraints (16) guarantee that each operating CCP must be allocated to at least one of the existing hospitals. The binary decision variables are given in (17) .

The second stage model
Once CCP location identification is done, the casualty response operations, including (i) triage, (ii) registration, (iii) treatment, and (iv) evacuation, commence with scenario realizations over the planning period T = { 1 , 2 , ..., |T | } . However, the purpose of the second stage stochastic formulation is to generate a robust design solution by involving different scenarios at the planning and operational level. To address the uncertainty of the number of casualties, a set of possible scenarios are generated and are then used in the model. Let ξ ilt (υ ) be the number of casualties identified with the injury severity level l ∈ L = { 1 , 2 , . . . , |L| } at the affected area i ∈ I on day t ∈ T under scenario υ ∈ ϒ. Although a wide range of injury severity level can be used, we divide the set of injury severity L into the following three subsets, L mi , L in , and L im , indicating minor injury severity, intermediate injury severity, and immediate injury severity, respectively. Note that the injury severity subsets are independent pairwise and that  (ω) as the flow of casualties characterized by all injury severity levels l ∈ { L mi ∪ L in ∪ L im } , from affected area i to CCP j in period t under scenario ω and F dir iklt (ω) as the flow casualties with minor injury severity level, i.e. l ∈ L mi , from affected area i directly to hospital k in period t under scenario ω. We also denote the outflows of casualties with injury severity level l ∈ { L in ∪ L im } from CCP j to hospital k in period t under scenario ω by nonnegative con- In the second stage, the set of all possible scenarios ω = ( υ, γ ) ∈ = ϒ × associated with the flow of casualties υ ∈ ϒ and the available transportation capacity γ ∈ are randomly generated from the historical data outside the optimization procedure.
For each scenario ω ∈ , the objective function (18) of the secondstage is the expected value of the total response planning and operational costs, involving casualty transportation cost from/to CCPs ( 18.1 ), casualty holding cost ( 18.2 ), and mortality cost ( 18.3 ) as follows. This objective function is subject to the system constraints (19) to (35) , as described in afterwards.
where, for each ω ∈ ,  (21) guarantee that the flow of casualties in the network is considered where the allocation in the network is certified.
where M is a positive large number.   (25) indicates the maximum treatment capacity of hospitals, indicated by c k , k ∈ K, for providing the required medical services to casualties coming directly from the affected areas and casualties transporting from the established CCPs.
Moreover, we consider the situation wherein the available transportation capacity at CCP j to cover inflows and outflows of casualties at injury severity level l ∈ L is uncertain due to failures, traffic congestion, accident, etc., in the roadways. This uncertain parameter is denoted by ζ jl (γ ) , where γ ∈ is the set of scenarios for the available transportation capacity. Constraint (26) indicates that the inflows and outflows of casualties, i.e. F in i jlt (ω) and F out jklt (ω) , respectively, at CCP j for injury severity level l ∈ { L in ∪ L im } cannot exceed the available transportation capacity under scenario γ . (27) takes into account the current uncertain flow of casualties under scenario υ ∈ ϒ with injury severity level l transporting from the affected areas to the established CCPs and the hospitals.

Casualties management constraints.
Casualty management operations emphasize the necessary functions including ( i ) registration, ( ii ) temporary hospitalization, and ( iii ) evacuation to hospitals or safer places, in the humanitarian logistics ( Lejeune & Margot, 2018 ). Considering this sequence of operations explained in the context of the problem, we define the scenario-based decision variables accordingly. Let C reg jlt (ω) indicate the number of casualties registered with injury severity level l at CCP j in period t under scenario ω. Constraint (28) guarantees that this latter does not exceed the inflows of casualties from the affected areas to a CCP.
For the medical treatment, the available capacity dedicated to each injury severity level of a CCP should be taken into account. This matter is represented in constraint (29) . Let us recall that C hos jlt (ω) represents the number of casualties for temporary hospitalization . It states that the number of casualties receiving temporary hospitalization services cannot be more than the dedicated capacity divisions at a CCP. Note that C hos jlt (ω) refers to the cumulative hospitalized individuals that corresponds to constraint (32) .

C hos
The medical services are immediately provided to the registered casualties diagnosed with injury severity level l. Depending on the severity of injuries l, the length of the hospitalization period, during which the casualties have to be kept and treated at CCPs, is denoted by τ l . After completing the hospitalization period τ l , these casualties become ready-to-evacuate to the corresponding hospitals. Constraint (30) reflects on the evacuation operations.
where C e v a jlt (ω) denotes the number of ready-to-evacuate casualties with injury severity level l at CCP j in period t under scenario ω. Constraint (31) certifies that the number of casualties transported from a CCP to the allocated hospitals cannot exceed the number of ready-to-evacuates . Note that only casualties with injury levels of intermediate and immediate have to be evacuated to hospitals, since they require further medical treatments.
Constraints (32) verifies the equilibrium casualty state transition in the consecutive periods in which the number of hospitalized casualties from the previous period plus the number of registered casualties of the current period is equal to the number of ready-to-evacuate casualties and the hospitalized casualties of the current period.

C hos
The most impactful output of humanitarian logistic network design is to save lives and reduce human suffering. This critical output is measured in our model by the following variables, M R jlt (ω) indicating the number of lives lost with injury severity level l ∈ L due to facility capacity limitation at CCP j in period t under scenario ω and M T jlt (ω) indicating the number of lives lost with injury severity level l ∈ { L in ∪ L im } due to transportation capacity limitation passing through CCP j in period t under scenario ω. Constraint (33) states that when casualty inflows are more than the CCP capacity to register, lives lost due to facility capacity limitation occurs. Similarly, Constraint (34) states that when the number of ready-to-evacuate exceeds the casualty outflows, lives lost due to the limitation in transportation capacity occurs.
The nonnegative continuous decision variables are given in (35) .

Solution approach
The solution approach proposed in this section is partly inspired from the sample average approximation (SAA) technique ( Shapiro, 2008 ), which is based on an approximation of the stochastic model by an equivalent deterministic mixed-integer programming (MIP) model. The methodology incorporates the SAA method, the robust counterpart problem and the feasibility restoration technique to solve the stochastic CCP network design problem with uncertain parameters.

Sample average approximation method
The scenario-based two-stage stochastic programming model represented above is a complex large-scale optimization problem, as a large number of scenarios is involved for uncertain parameters realization. To solve the two-stage stochastic CCP network design problem represented above, we are inspired by the SAA technique ( Shapiro, 2008 ), which is based on an approximation of the stochastic model by an equivalent deterministic mixedinteger programming (MIP) model. The SAA model incorporates the equivalent deterministic mixed-integer program of the secondstage decision-making problem into the first-stage decision-making problem. The SAA method has mainly been used to find nearoptimal solutions for two-stage stochastic problems ( Amiri-Aref et al., 2018; Klibi & Martel, 2013;Schütz, Tomasgard & Ahmed, 2009 ).
Since two sets of uncertain parameters are concerned in this paper, i.e. the number of casualties and available transportation capacity, two sets of scenario generation should be realized in this model. By generating N 1 independent number of casualty scenarios given as { υ 1 , υ 2 , . . . , υ N 1 } = ϒ N 1 ⊂ ϒ, and N 2 independent available transportation capacity sce- Given the original two-stage stochastic model (12) -(35) , the SAA program is constructed in the following: where, the first three terms in (36) denote the first-stage objective function and the last term denotes the expected objective function of the second-stage problem. The SAA method is performed when a feasible solution exists and the problem has a finite objective value ( Shapiro, 2008 ). However, the uncertain parameters in humanitarian logistics may not have an identical distribution or a known distribution parameter. In such a situation, the SAA method is prone to return infeasible solutions by violating some of the constraints in at least one scenario. To tackle this challenge, we provide a robust counterpart problem for the represented SAA method, involving robust solution and robust model, proposed by  , in the following subsection.

Robust SAA method
A robust solution is characterized by its proximity to the optimal solution of a stochastic programming model. We incorporate solution robustness by the inclusion of the mean absolute deviation of the second-stage solutions, indicated by (ω) , over the number of scenarios in the SAA model, as follows: Let us recall that q ( x , ω ) is the second-stage decision-making problem. As discussed earlier, expression (37) has to be minimized to achieve solution robustness. Therefore, it is included in the objective function of the SAA model. As it contains the absolute function which makes the SAA model nonlinear, we apply a linearization approach to guarantee the convexity of the solution space. ( 37 ) is included in the minimization objective function, we can substitute it by the following expressions:

Proposition 1. As the expression
where, when minimizing expression ( 38 ). In this case, . For more information regarding this linearization method, refer to Yu and Li (20 0 0 ).
A robust model is regarded as a model that returns solutions which are feasible for any given scenario realizations. Due to the variability of the uncertain parameters, a stochastic programming model might be infeasible for some scenario realizations. One of the most probable reasons for infeasibility in a stochastic programming model is the variability of scenario realizations, which corresponds to the inflows of casualties ( Birge & Louveaux, 2011 ). This issue, which is coupled with the limited available physical capacity of each potential node for establishing a CCP, corresponds to constraint (23) . In fact, this constraint verifies the additional CCP nodes are required for accommodating the inflows, as the existing areas of potential CCP nodes are not sufficient. To overcome this issue, we apply a model robustness approach, in which an infeasibility variable z j (ω) is taken into account in the system constraints, as represented in (9) . The infeasibility variable z j (ω) shows the amount of infeasibility of each scenario ω ∈ in the model. It is clear that z j (ω) = 0 if the model is feasible. Otherwise, it returns a positive value. However, a huge penalty number ρ is assigned to the infeasibility variable z j (ω) in the objective function of the model to avoid being infeasible for all scenarios. We then modify the constraint (23) , which refers to the j-th CCP capacity limitation, by adding the infeasibility variable z j (ω) , as follows: Considering that, the SAA model with the robust optimization techniques, namely solution robustness and model robustness, is represented in the following: If there exists at least one infeasibility variable with a nonzero value for any scenarios, the results are meaningless and inapplicable. That is to say, the model does not guarantee that obtaining solutions satisfy the system constraints for all scenario realizations and do not converge to the optimal solution. This failure can be partly due to the inappropriate set of location and allocation decisions or inadequate capacity acquisitions in the network structure in the first-stage decision-making problem. In other words, not all choices of design decisions x ∈ R n + 1 give rise to feasible solutions. To achieve feasible solutions when the infeasibility variables return nonzero values, we apply a feasibility restoration technique on the non-algebraic constraints, i.e. design decisions, which is discussed in the following section.

Feasibility restoration technique
As casualty flow in humanitarian logistics is unpredictable, in case of failure in the robust SAA model, the feasibility restoration technique proposed in this paper allows us to reconsider the CCP logistic network structure and adopt appropriate design decisions accordingly. It is clear that the operational decisions at the second stage will improve as a result of improvement in the design decisions.
The feasibility restoration technique is inspired from the work of Abramson and Randall (1999) , which has been further developed in Casey and Sen (2005) , and applied in Huang and Mehrotra (2016), Kim and Wright (2016) and Lee, Liu, Mehrotra and Bie (2015) . This technique is characterized by detecting infeasibility and incorporating auxiliary design decision variables in the twostage program to tackle the issue while considering all scenario realizations. The key feature of this technique is to expand the network configuration so that feasible solution is enhanced and can reproduce more efficient objective value.
Proposition 2. Let us recall that x = ( X j , Y i j , V jk , U ik ) ∈ R n + 1 denotes the vector of the first-stage binary decision variables, where j ∈ J = { 1 , 2 , . . . , |J | } represents potential locations to establish CCPs.
To redesign the network structure, we need a modified set of potential locations. We introduce 1 is a vector of binary variables where j s ∈ J s = { 1 , 2 , . . . , | J s | } represents the modified set of potential locations.
Note that J and J r are independent sets and that J ∩ J r = ∅ . Also, note that J ∪ J r = J s and that | J s | = |J | + | J r | .
As a result of Proposition 2 , the new pooling of design decision variables gives rise to the evolution of the control variables accordingly. We introduce the evolving control variables by y s (ω) ∈ R n + 2 as a vector of non-negative variables and subsequently the evolving absolute deviation function and infeasibility variable, s (ω) and z s (ω) , respectively. According to the Proposition 2 , we then reconstruct the robust SAA programming, as represented in (44) -(48) .
As the extended (44) -(48) are partly similar to those already described, we avoid repeating the description of the above model.
Increasing the size of feasibility restoration variable set to | J r | , where r r , allows the model to choose the most appropriate locations among the available nodes, although it increases the problem complexity.
A general computational framework for the robust stochastic optimization under uncertainty is outlined in Fig. 3 .

Validation analysis
In this section, we discuss a validation analysis which is based on optimality gap estimation between the objective value at a solution found by the proposed algorithm and the optimal value of the true problem. The optimality gap estimation is a way to evaluate the quality of stochastic solutions in two-stage programming where the true objective value is finite and the second-stage solution is feasible for almost every realization of the random data.
We suppose x T and y T denote the true optimal solutions of the first-stage and the second-stage problem and f ( x T , y T ) is the true optimal objective value. According to Shapiro (2008) , since finding the value of f ( x T , y T ) is almost impossible, as enormously large number of scenarios are required, statistical lower and upper bounds for the true optimal objective value using the valid inequality can qualify the solution procedure. The statistical lower bound is estimated by averaging the solutions of the algorithm in M independent times based on N generated scenarios and a valid statistical upper bound for the true optimal objective value is given by sampling . This latter can be done through solving the secondstage problem using a large enough sample of scenarios N N, where the solution of the first-stage problem is given as input.

Averaging procedure
Let x m N and y m N , m = 1 , . . . , M, denote the optimal solution vector of the two-stage stochastic problem found by the algorithm with scenario sample size N in the m th replication of sample generation, and f ( x m N , y m N ) be the optimal objective value corresponding solution values. We then provide average and standard-deviation estimators for the true objective values. An unbiased estimator of the statistical lower bound of the expected true objective value, Validation procedure: Step 1.
Calculate Optimality gap Calculate the statistical optimality gap percentage given in (55). If the gap is acceptable, stop; otherwise increase N and/or M and return to step 1 Output Statistically valid bounds on the true objective value (with confidence at least 1 − 2 ). as follows: Considering M independent scenario generations, the standard deviation is estimated in the following: Using the average and standard deviation estimators for M replications of samples of size N , an approximate ( 1 − α) × 100% confidence lower bound of the true objective value, denoted by L M N , is given as follows: where t α,M−1 represents the α-critical value of the t -distribution with M − 1 degrees of freedom.

Sampling procedure
The statistical upper bound of the expected true optimal objective value can be estimated by sampling procedure. Let x be the best optimal solution vector of the first-stage problem found by the algorithm with a scenario sample size N among M replications. We then solve the problem with x as an input and generate sample scenarios { ω 1 , ω 2 , . . . , ω N } ∈ N ⊂ , where N N, which are independent to samples used in computing x . It is clear that when x is given as an input, the problem can be decomposed into N deterministic problems. We denote the optimal objective value based on a sample size N by ˆ f ( x , y * N ) and the optimal objective value solved one at a time by ˆ f ω ( x , y * ω ) , where ω ∈ N ⊂ . Note that y * N and y * ω represent the solution of the second-stage problem when N sample scenarios are involved and the scenario-wise solution of the second-stage problem, respectively. One can calculate the standard deviation of ˆ An approximate ( 1 − α) × 100% confidence upper bound of the true objective value, denoted by U N , is then given as where n α represents the standard normal critical value with ( 1 − α) × 100% confidence level. Therefore, an approximate ( 1 − α) × 100% confidence interval for the expected true objective value is represented in the form of ( L M N, 1 −α , U N , 1 −α ) , using Eqs. (51) The validation procedure discussed above is then summarized in Fig. 4 .

Computational study
The robust stochastic optimization modelling approach described in Section 3 is implemented through a computational study using data scenarios modelled on the Bhopal gas tragedy that occurred in India over three decades ago. More specifically, we consider a hypothetical case of a gas leak in Bhopal in today's date and which follows the hazard propagation profile (e.g., wind direction, affected wards) reported back in 1984. The underlying data for the study, which includes the population of specific wards (population areas/catchments), available transportation in the city, existing infrastructure (including schools and hospitals), open spaces, and other model-specific parameters, was obtained through census data and from local municipal reports. We conducted one field trip to get access to some of this information. The data thus obtained was used to estimate the required parameters, which were then used to model the scenarios for the computational study. In this section, we also discuss the efficiency of our proposed modelling approach and present the solution sensitivity analysis to provide further insights to humanitarian logistics planners and practitioners.

Context for study
This section briefly describes the Bhopal Tragedy in India, often known as the worst industrial accident in the world and provides a computational investigation into the humanitarian logistics network design for establishing CCPs in the affected areas. On December 3, 1984, a highly toxic cloud of methyl isocyanate (MIC) leaked from a pesticide plant in Bhopal, the capital city of the state of Madhya Pradesh, the second largest state in India. The leak was the consequence of a large volume of water entering one of the methyl isocyanate storage tanks around 9:30 pm the day before. This triggered off a chemical reaction resulting in a tremendous increase of temperature and pressure in the tank and consequently led to an explosion. More than thirty years have passed since the gas explosion, but the Bhopal saga is far from over. During our trips to the plant site and conversations with the volunteers at the NGO clinics as well as the local slum dwellers, we were told that of the 80 0,0 0 0 people living in Bhopal at that time, no one knows exactly how many people were affected that night (Veron and Nanda, 2014) .
The geographical scope of our study focuses on the affected areas in the city of Bhopal. According to the technical report of the Indian Council of Medical Research ( ICMR, 1985 ) on the Bhopal disaster, |I| = 33 wards have been identified as affected areas with more than 70 0,0 0 0 population ( Pradesh, 2011 ) (each ward is shown as an orange icon in Fig. 5 ). Within this area, a set of predesignated locations have been selected as the candidate points to establish CCPs. These CCP points are usually sites that can accommodate a large number of casualties ( Drezner, 2004 ), for exam-ple, college and university campuses, high schools with a football field, mosques, malls and large parks. We identified a total of 65 CCP candidate points, including existing buildings and open-spaces (shown as blue icons in Fig. 5 ). The capacity of each potential CCP location to provide medical services to casualties is estimated by its total available area divided by the space required to treat per person. We considered the latter equal to μ = 7 m ² per person as reported in the statistical report ( Moore, Levit & Elixhauser, 2014 ). Moreover, the network includes |K| = 9 hospitals and medical care centres as safe places to evacuate the casualties for further treatment (hospitals are shown as a white cross in a purple circle). Union Carbide plant, i.e. the disaster point, is shown using a yellow icon. For more details about the case study, refer to Appendix A.
We assume that the available transportation capacity by means of ambulances for immediate severity injury level is 500 people per trip, 2 which was far below the required capacity to move mass casualty in the disaster we considered. Therefore, we considered the public/private vehicles (including mini buses, standard buses, and private cars) into the transportation capacity to move mass casualty with minor and intermediate severity injury level to CCPs and hospitals. Using both public and private modes of transport, the available transportation capacity reached more than 20 0,0 0 0 people. To generate random number of casualties, we utilized the simulation procedure provided in Singh and Ghosh (1987) . Analysing the data, we observed that the coefficient of variation of generated number of casualties of 7 wards out of 33 was about 80%, while it was above 120% for 26 wards out of 33, which represents a considerable uncertainty inherent in the generated number of casualties.
Note that t r H and t r L refers to high and low transportation capacity scenarios, respectively. Therefore, the combination of various values of compromising coefficients β 0 , β 1 , and β 2 , and two usable transportation capacity percentages t r H and t r L yields 144 problem instances. For each problem instance, we generated, based on the population of each ward and the number of casualties reported in Singh and Ghosh (1987) , N 1 = 5 independent number of casualty scenarios in |L| = 3 level of injury severity for each ward and N 2 = 3 independent available transportation capacity scenarios for each CCP, over a planning period of |T | = 7 days. In other words, for each problem instance, N 1 . N 2 . |T | . |L| = 315 sample scenarios are generated to represent the number of casualties for each ward.

Numerical results and discussion
The instances described in Section 4.1 are solved after scenario generations on a 64-bit operating system server with a 2.7 gigahertz CPU on Intel(R) processor and 72 gigabytes of RAM. The proposed robust stochastic optimization approach, shown in Fig. 3 , is performed using the optimization solver GAMS with a MIP Relative Tolerance of 0.005 within a 5-hour computation time. The detailed numerical results, including the solution value and computational time, related to the 144 instances are represented in Tables B1-B8 of Appendix B.
In order to measure the efficiency of the proposed logistic network design and related operations, we applied important metrics related to disaster management. We present the results in the following sections.

Locational decisions
We compare CCP location decisions found by SAA method and the proposed robust stochastic optimization with feasibility restoration variables with respect to coefficients of β 0 = { 0 . 1 , 1 , 10 , 100 } . This comparison is illustrated in Figs. 6 (a) and (b) showing two levels of available transportation capacity after a disaster strikes. In these figures, the coefficients of β 0 can be considered as the risk aversion attitude of a decision maker (DM), where 0.1 attributes to a risk incentive DM and 100 relates to a risk aversive DM and is represented on the x-axis. The average number of opened CCPs over the number of involved instances is represented on the y -axis. Results illustrated in Fig. 6 (a) show that, when on average 90% of casualty transportation capacity is available, the SAA method opens 39 locations for establishing CCPs, on average, and is not sensitive to the risk aversion attitude of a DM. Compared to this, our proposed methodology suggests opening up to 43, on average, locations for establishing CCPs and is fairly relative to the risk aversion attitude of a DM. CCP location decisions are more of the essence when the available casualty transportation capacity decreases to 80%, on average. Our findings show that the output of the SAA method remains unchanged even when casualty transportation capacity is reduced by 10%. However, by using our pro-posed algorithm a significant increase in the number of CCPs is observed, which contributes to 47 locations for establishing CCPs in the case of risk averse DM ( β 0 = 100 ) -refer to Fig. 6 (b). In other words, the results reveal that the more conservative a DM is, the more the number of CCPs that will need to be operationalised.
Further, using our proposed algorithm, as the coefficients of β 0 increases, the number of existing building selected for establishing CCPs decreased and instead more potential locations are chosen from open-space spots as locations to set-up CCPs. The information on buildings and open spaces was based on data from our Bhopal case study.
From Fig. 6 (a) and (b), it can be concluded that our proposed methodology, which is based on robust stochastic programming model, enables a DM to cope with the infeasibility issues due to the dispersion of scenarios and generates more efficient solutions which are feasible for any scenarios. Furthermore, being more risk averse in an uncertain decision-making environment results in opening more CCPs among the existing buildings and open-spaces and therefore being closer to the affected areas. This fact emphasizes the necessity of providing fast and efficient medical services to the casualties from the shortest possible distance. In the following section, the role of accessibility to the services in CCPs and its impact on the number of lives lost is explored.

Network structure decisions
In order to measure the quality of a complex emergency network design, we introduce the proximity metric which is defined as the total distance travelled in the network to the number of links associated with all pair nodes, i.e. from the affected areas to the established CCPs, also known as the average path length. Proximity is an important metric in humanitarian logistics and has been extensively used in this context ( Muggy & Stamm, 2017 ). Let us indicate the solution value of allocation decision variables by ˆ Y i j . The proximity metric is then formulated as i which represents the average path length to reach a CCP. Results illustrated in Fig. 7 reveal that our proposed robust optimization method designs a network in which the average path length (shown as a dash-line in Fig. 7 ) has improved in comparison to the SAA method. This can be confirmed by the results of increasing in the number of CCPs that are opened, as illustrated in Fig. 6 . We then investigate the number of lives loss, also known as mortality in this work, to see whether it is influenced by the average path length improvement. As shown in Fig. 7 , on average, the mortality rate experienced a significant reduction from 438 individuals to 294 individuals due to the decrease in the average path length. In general, Fig. 7 suggests that a small improvement in proximity to CCPs can result in a significant decrease in the number of lives saved. This analysis also addresses the Equity , also known as fairness, which tackles the discoordination of operational decisions for providing appropriate emergency services to casualties. When it comes to relief contexts, this metric measures the unsatisfied demand associated with each operational decision over the planning period ( Marsh & Schilling, 1994 ), which refers to mortality in our case. Moreover, the average path length which denotes the rapidity is also widely considered as the equity metric ( Anaya-Arenas, Ruiz & Renaud, 2013 ). It can be interpreted from Fig. 7 that overall, the equity metric has been improved by the modelling approach we proposed in this work.

Robust performance metrics
As discussed earlier, robust optimization approach enables DMs to generate solutions while reducing the risk of dispersion and ensuring the solution concentration in an uncertain environment. In this work, we measure the dispersion of the objective func-   tion values over all given scenarios, found by the proposed model, as a metric to evaluate the solution robustness. This metric gives rise to the standard deviation of the objective values which represents the closeness between them. Results illustrate that the dispersion of objective values in the uncertain environment increases, as the transportation capacity contributes to 10% reduction, on average. Results also suggest that the standard deviation of our proposed optimization approach is slightly larger than the stochastic programming method; this can be due to network expansion and the resultant distribution of entities throughout the optimized network. Overall, the dispersion of the objective value in both cases, i.e. SAA method and the proposed stochastic robust optimization method, are negligible (less than 10 -8 ).
Another important factor that is used in robust optimization approaches, also known as model robustness, is to generate solutions values which satisfy all system constraints for any given scenarios. Due to the uncertainty inherent in mass casualty flow management, it is very likely to observe the infeasible solutions. We also evaluate the infeasibility produced in the model for the 144 instances when using the stochastic modelling approach and compare with the corresponding values when applying our proposed solution algorithm by the usable transportation capacities. It has been observed that the stochastic programming (SAA method) approach results in solutions where positive values of infeasibility exist, on average 3.5907445 ×10 5 , whereas the proposed approach ended up with zero infeasibility values. It is also found out that the infeasibility values corresponding to the SAA method increase as transportation capacity tends to decrease. For detailed information, refer to Tables B1-B8 of Appendix B.
In relation to robust optimization, the overall performance and reliability of solutions are measured by calculating the coefficient of variation, i.e. standard deviation-to-mean ratio, through all the scenarios ( Birge, 1982 ). We then calculate the coefficient of variation corresponding to decision variables used in the model over 144 instances and represent the minimum, mean, and maximum value of the coefficient of variation of each variable over all instances (see Table 1 ). Results show that the coefficient of variation of all operational decisions are considerably low, such that, for the majority of them it is less than 1%. However, the coefficient of variation value corresponding to the strategic decision of capacity allocation is 2.15% on average and which is not too large. In general, it shows that the solution values have low variability and are quite reliable.

Validation metrics
In this section, the validation procedure, represented in Fig. 4 , is used to examine the accuracy of the solutions found by the proposed robust optimization solution method. All instances are tested and their associated statistical optimality gap values are computed according to the validation procedure. A lower bound solution with 95% confidence level is computed using the averaging procedure with replication size M = 4 and scenario size N = 15 . Then, using the best solution found from the average, the sampling procedure is applied with sample evaluation scenario size N = 150 to generate an upper bound with 95% confidence level. We then calculate the statistical optimality gap percentage for each instance. The results are reported in Tables B1-B8, in Appendix B. To provide a clear view of the optimality gap percentage over the instances and its relationship with DM risk aversion attitude, we represent the average of optimality gap percentage over the instances for each corresponding value of β 0 = { 0 . 1 , 1 , 10 , 100 } in Fig. 8 . As can be observed, the optimality gap has a decreasing trend as the weight corresponding to DM risk aversion attitude increases. It is due to the fact that instances with higher weight of DM risk aversion attitude have the objective function with low variability and therefore with less optimality gap. It can be concluded that the more conservative the DM is, the less the optimality gap that exists.
We finally compare the convergence rate of the proposed solution methodology to that of SAA method by reporting the dual objective value and the best integer bound found by the solver in each iteration corresponding to an instance in Fig. 9 . Results represent that the proposed algorithm converges to optimal solutions after about 10 0,0 0 0 iterations while the corresponding number re- lated to the SAA method is over 20 0,0 0 0, which shows the fast convergence rate of the proposed algorithm. This is due to the fact that the feasibility restoration technique is able to facilitate the proposed stochastic robust optimization approach to perform more efficiently and rapidly.

Conclusion
In this paper, a two-stage stochastic programming model has been formulated for the Casualty Collection Point network design problem that is based on the 1994 Bhopal gas tragedy. The number of causalities and the available transportation capacity were the uncertain parameters of this problem; they were generated using an existing simulation model from literature and resulted in a high variability of number of casualty scenario realization. To tackle this issue, we have proposed a stochastic robust optimization approach with the feasibility restoration technique, inspired by the SAA method, and an extensive computational experiment has been conducted for this case problem. The performance of the solution approach has been tested by the validation procedure commonly used in stochastic programming.
The experimental results reveal some practical and managerial insights confirming the importance of CCP logistical network design and operational response decisions in an uncertain environment. The findings show that the network configuration obtained by our proposed methodology has a significant difference with the SAA method. More specifically, the proposed approach opens more CCPs and is more sensitive to the transportation capacity; this can be contrasted with the SAA method where no significant sensitivity has been observed. We notice that a conservative decision maker (DM), with risk aversion attitude, tends to open more CCPs in an uncertain decision-making environment.
The proximity metric has been quantified as the average path length to a CCP in the network structure for all instances. It has been observed that the network configuration by our methodology enables a DM to improve the proximity metric in a CCP logistical network design. We notify that a small improvement in the proximity metric can result in a significant increase in the number of lives saved. Results also show that reduction in transportation capacity in stochastic programming can lead to increasing the dispersion of the solutions, however, our stochastic robust optimiza-tion approach is able to achieve solution and model robustness approaching the optimal solutions. We realize that the optimality gap in the stochastic programming can be improved by taking risk aversion attitude which results in less variability of the objective values.
Our future research will investigate a hybrid simulationoptimization approach for casualty evacuation based on CCP network structures that has been identified in this work. The inclusion of the medical supply flow from the multiple available hospitals to the established CCPs for the purpose of casualty treatment can be another direction to develop this problem towards a more realistic context ( Haynes & Freeman, 1989 ). In this regard, simulation approaches like Discrete-event Simulation (DES) could be used for modelling of healthcare supply chains ( Mustafee, Taylor, Katsaliaki & Brailsford, 2009 ). Yet area of interest is the use of qualitative system dynamic at the tactical level as an alternative to the scenario generation in the optimization model to overcome the complexity of the problem ( Powell, Mustafee, Chen & Hammond, 2016 ). An extension to the robust minmax regret stochastic programming model can be another interesting research topic to consider in the humanitarian logistics network problem ( Feizollahi & Averbakh, 2013 ). As the casualty accessibility to CCPs plays an important role in humanitarian logistics, a maximal accessibility network design can be further extended ( Aboolian, Berman & Verter, 2015 ). These are all future directions to the work presented in this paper.

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.ejor.2019.06.018 .