A Lexicographic Approach to Postdisaster Relief Logistics Planning Considering Fill Rates and Costs under Uncertainty

Predicting the occurrences of earthquakes is difficult, but because they often bring huge catastrophes, it is necessary to launch relief logistics campaigns soon after they occur. This paper proposes a stochastic optimization model for post-disaster relief logistics to guide the strategic planning with respect to the locations of temporary facilities, the mobilization levels of relief supplies, and the deployment of transportation assets with uncertainty on demands. In addition, delivery plans for relief supplies and evacuation plans for critical population have been developed for each scenario. Two objectives are featured in the proposedmodel: maximizing the expected minimal fill rate of affected areas, where the mismatching distribution among correlated relief demands is penalized, andminimizing the expected total cost. An approximate lexicographic approach is here used to transform the bi-objective stochastic programming model into a sequence of single objective stochastic programming models, and scenario-decomposition-based heuristic algorithms are furthermore developed to solve these transformed models. The feasibility of the proposed bi-objective stochastic model has been demonstrated empirically, and the effectiveness of the developed solution algorithms has also been evaluated and compared to that of commercial mixed-integer optimization software.


Introduction
Earthquakes are a special type of natural disaster that can result in huge catastrophes.The Great Sichuan Earthquake, which occurred on May 12, 2008, in Wenchuan County, Sichuan Province, China, was one such case.It had a recorded magnitude of 8.0 and imposed tremendous suffering on the local residents, causing over 69,000 deaths, 18,000 missing persons cases, 374,000 injuries, and huge loss of property.Although thousands of networked seismograph stations have been installed around the world and although these stations use powerful computers, it is still difficult to predict when and where an earthquake will strike.Therefore, an effective and necessary approach coping with an earthquake disaster is to plan and implement a disaster relief campaign that would reduce the damage right after the earthquake disaster takes place.
Emergency logistics plays a vital role in disaster relief campaigns.However, the planning and implementation of relief logistics are challenging, especially when help is needed in mountainous areas, where the transportation infrastructure, such as roads, bridges and railway lines, can suffer severe damage after large earthquakes.The following three factors must be taken into account.First, the demand for disaster reliefs in scattered affected areas (AAs) can be uncertain early after earthquake [1].This is because collecting actual and real-time information from AAs is relatively infeasible due to power disruption, damage to communication equipment, and severe telecommunication traffic congestion in the disaster region.Second, as the transportation infrastructure is usually seriously damaged and repairing the damaged roadways quickly is very difficult within the mountainous regions, traditional transportation tools, such as trucks and trains, are impractical.Instead, helicopters become the most useful transportation means within the disaster regions during the initial postdisaster period (i.e., the first 72 hours after a disaster, a golden relief period that is critical to provide relief since communities are not expected to stand on their own for much longer than that time).In fact, it has been reported that 99 helicopters were mobilized for rescue missions after the Great Sichuan Earthquake [2].The last factor is the optimization objectives of relief logistics.Conflicting objectives, such as maximizing the fill rate of AAs and minimizing costs, should be contemplated concurrently.Although in this emergency situation minimizing costs should not be emphasized in planning the logistics, the excess accumulation of supplies can be effectively avoided by incorporating this objective.Otherwise, overaccumulation of resources may delay or disrupt the activities of relief logistics.
In addition, because earthquakes are more unpredictable than other types of natural disasters, such as floods and typhoons, it is usually uneconomical or even impossible to maintain enough capacity for relief operations before the occurrence of an earthquake.As a result, temporary facilities (TFs) must be established as distribution depots at the edges of disaster regions in order to promote these relief supplies flowing into the affected areas in a controllable and fluid way.
The purpose of this paper is to model the strategic planning of relief logistics in response to a devastating earthquake disaster in a mountainous region.It addresses specific issues such as the inherent uncertainty regarding the degree of demand, conflicting objectives, and helicopterbased distribution plans.The main contributions of this paper can be summarized as follows.
(i) Development of a novel optimization model to address the disaster relief logistics problem as a multiobjective, stochastic, mixed-integer, nonlinear programming problem: this model combines temporary facility locations, mobilization levels of relief supplies, and initial deployment of helicopters with uncertainty about demands.It also includes helicopter-based delivery plans for relief supplies and evacuation plans for critical population (referring to those in need of emergency medical evacuation) under each revealed scenario.The formulation not only takes into account the expected total cost of relief logistics but also considers the fill rate of AAs and the synergy among correlated relief demands and the fairness among AAs in distribution process.
(ii) Proposal of a decomposition-based heuristic algorithm based on the structure of the proposed model: as a result, the stringent time requirements for solving the model can be satisfied in such an emergency environment of disaster response.
(iii) Application of this model to a real-world disaster relief chain, specially the "5.12"Wenchuan Earthquake relief chain, which was dedicated to the relief supply and medical evacuation after the disaster.
The rest of this paper is organized as follows.Section 2 provides a review of disaster relief optimization methods related to the present work.In Section 3, a multiobjective, stochastic, mixed-integer, nonlinear programming disaster relief model is described.Then, we develop a decompositionbased heuristic algorithm suitable for solving the model in Section 4. This is followed by computational results of model testing in Section 5. Conclusions and directions for future research are provided in the last section.

Literature Survey
Relief logistics is essential to the timely and effective provision of resources to aid people affected by natural disasters.To shorten the response time of disaster relief, a natural approach is to establish facilities and stock prepositioning before the strikes of disasters.However, in making decisions for this strategic planning before the disaster actually takes place, uncertainties, such as the degree of demands and transportation conditions, are the key challenges.The significance of these uncertainties has motivated a number of researchers to address stochastic optimization in disaster relief planning using probabilistic scenarios representing disasters and their outcomes.
Considering the uncertainty of demands from potential disaster areas, Balcik and Beamon [3] developed a maximally covered location model of prepositioning relief commodities to determine the number, locations and capacities of the relief distribution centers, in which the pre-and postdisaster budget constraints were specially considered.Chang et al. [4] presented a facility location-allocation model for stocking and distributing rescue resources taking possible flood scenarios into account, and this model can be utilized as a decision-making tool in generating plans for flood emergency logistics.
Beside the uncertainty on demand, some studies further considered other uncertainties, such as survival of prepositioned stocks, transportation network conditions, and distribution costs.Ukkusuri and Yushimito [5] presented a model to identify locations for prepositioned supplies, and its objective was to maximize the probability that demand points could be reached from at least one supply facility under conditions where facility locations and links in transportation network may be destroyed.Rawls and Turnquist [6] developed a two-stage stochastic mixed-integer program that determined the locations and quantities of various types of emergency commodities, in which uncertainties of demand and transportation network availability were considered.Rawls and Turnquist [7] also extended this model with additional service quality constraints that ensured the probability of meeting all demands was at least a preset level, and the distance between a supply-demand pair was less than a specific limit.Mete and Zabinsky [8] proposed a stochastic optimization model for disaster preparedness and response with uncertainties on demand and cost, in order to assist in deciding the location and allocation of medical supplies to be used during emergencies in the Seattle area, which is known to be vulnerable to earthquakes.Salmerón and Apte [9] presented a two-stage stochastic optimization model for planning the allocation of budget on acquiring and positioning relief assets; in this model the first-stage decisions included the expansion of resources such as warehouses, medical facilities, ramp spaces, and shelters, whereas the second-stage decisions focused on the transportation plans for commodities and critical population; this model also contemplated the transportation network availability following a disaster under uncertainties of demand and cost.Bozorgi-Amiri et al. [10] developed a multiobjective robust stochastic programming model for disaster relief logistics, in which demands, supplies, and the costs of procurement and transportation were considered as uncertain parameters: it attempted to minimize the sum of the expected value of total costs and maximize the affected areas' satisfaction levels.Döyen et al. [11] developed a twostage stochastic programming model for humanitarian relief logistics problem where decisions are made for pre-and postdisaster rescue centers, the amount of relief items to be stocked at the predisaster rescue centers, the amount of relief items flows at each echelon, and the shortage degree of relief items with uncertainties on demand and transportation time/cost.Unlike these previous works, which gave extra consideration to uncertainties other than demand, such as supplies and transportation conditions, this paper merely considers the uncertainty of demand.This is because our work mainly focuses on the strategic planning of relief logistics that an earthquake disaster has struck, and under these conditions, the amounts of supplies available tend to be relatively definite.(However, there are also cases where the availability of supplies is unpredictable even after an earthquake disaster, especially if it depends on the willingness of private international donors to financially support the relief operations.In most situations, the assumption that the availability of supplies is relatively definite after a disaster is reasonable.)In addition, because helicopters are utilized as the preferred unique emergency vehicle transportation tools in this model and because the physical conditions of roads have no significant effect on the ability of helicopters to travel, the uncertainty of transportation conditions was not taken into account here.
Most existing stochastic models on prepositioning of emergency supplies for disaster response belong to twostage mixed-integer stochastic programming problems, and decomposition algorithms can be adopt to solve them, such as the L-shaped algorithm based on Bender decomposition in Birge and Louveaux [12] and the dual decomposition algorithm based on scenario decomposition and Lagrangian relaxation in Carøe and Schultz [13].Besides the above two classical decomposition algorithms, some other decomposition algorithms have also been developed for mixedinteger stochastic programs with special structures: stochastic program with purely binary first-stage decision variables and arbitrary second-stage decision variables [14], pure continuous or pure binary first-stage variables and 0-1 mixed-integer recourse variables [13], mixed-integer first-stage variables and pure integer second-stage variables [15], pure binary firststage variables and 0-1 mixed-integer recourse variables [16][17][18], and 0-1 mixed-integer variables for both the first and the second-stage variables [19,20].However, as to the best of our knowledge, only a few studies have been performed on the development of decomposition methods suitable for humanitarian logistics problems [6,11,21].Although the computation complexity of most facility location models preparing for disasters is high when the sizes of optimization problems are large, the time required on solving such a model is relatively insignificant because these models solely focus on strategic planning of relief logistics for potential disasters that have not yet stuck.In contrast, as our mixed-integer stochastic model targets the strategic planning of postdisaster relief logistics, there exists stringent requirement on time for solving this model in such an emergency situation.As a result, we specially developed a scenario-decompositionbased heuristic algorithm to solve the proposed model, in order to satisfy the stringent requirement on making an emergency decision in a short time.
In conclusion, our contributions or differences compared with related existing studies can be concluded as two aspects.The first is the novelty of the proposed biobjective stochastic programming model itself.Most current stochastic emergency facility location and facility location-allocation models were designed for predisaster preparedness and response with single objective, such as minimizing expected cost [6,7,11], maximizing the coverage of potential disaster areas [3], or maximizing the probability on reaching the demand points [5], while our postdisaster stochastic facility locationallocation model has two objectives on expected fill rate and cost.Although the model proposed by Bozorgi-Amiri et al. [10] has similar objectives as ours, its two robust criteria, the variance of total cost and satisfaction level among scenarios, and the solutions infeasibility due to parameter uncertainty, were incorporated in the objectives of their model, whereas the synergic issue among related demands was extraconsidered in the first objective of our model.In addition, the coordination between the mobilization levels of transportation tools (i.e., helicopters) and the amount of relief supplies were especially considered in our work for postdisaster relief logistics planning, and this coordination had not been contemplated in the aforementioned literature because they mainly focused on preparations to be made before disasters.The second is about the solution approach to solve the proposed biobjective model.An approximate lexicographic approach was applied to our biobjective model, while the technique of Compromise Programming (CP) was adopted in Bozorgi-Amiri et al. [10], in which identifying the value of the weight coefficient  was not easy.Two scenario-decompositionbased heuristic algorithms were further developed for the two transformed single objective stochastic programming problems, respectively, to satisfy the stringent requirement on time for solving the problem, and in contrast, there is no similar work on algorithm design in Bozorgi-Amiri et al. [10].Moreover, the problem decomposition technique utilized in this paper is scenario decomposition, whereas in Rawls and Turnquist [6] it is Benders decomposition and in Döyen et al. [11] is by relaxing special constraints to obtain two subproblems.In other words, the proposed solution methods for our stochastic optimization problem based on scenariodecomposition heuristic algorithm within an approximate lexicographic optimization frame are novel compared to existing work for emergency relief logistics planning.

Mathematic Formulation.
The model is constructed with two objectives based on the stochastic programming approach, in which uncertainty is represented by a set of finite discrete scenarios: s.t.
The two objectives are represented by ( 1) and ( 2), respectively.They are explained as follows.
Objective 1.The first objective maximizes the expected minimal fill rate of affected areas while penalizing the mismatching distribution among correlated relief demands.The first term in ( 1) is meant to maximize the AAs' expected weighted fill rate by maximizing the sum of the minimum fill rate at AAs, and this max-min method is employed to avoid unfairness in distribution among AAs.The second term in (1) is the expected penalty due to mismatching distribution among correlated relief supplies, and the weight coefficient     is set larger if the synergistic requirement between  and   is higher.
Since this objective is explicitly nonlinear, the linear equivalent equations can be rewritten with the help of an auxiliary variable    as follows: Objective 2. The second objective function minimizes the sum of four types of costs: the cost of establishing facilities, the cost of procuring relief supplies, the cost of recruiting helicopters, and the expected cost of transportation.Equations ( 3)-( 6) define all the constraints of the scenario-unrelated decision variables.Equation (3) ensures that the mobilization level of any relief supply at any closed TF should be zero.Equation (4) bounds the maximum mobilization levels of relief supplies at any open TF.Equation (5) states that helicopters can only be assigned to open TFs.Equation (6) guarantees that the number of helicopters of each mode allocated to all TFs cannot exceed the total number of helicopters available.Equations ( 7)-( 15) define all the constraints related to scenario-specific decision variables.Equation ( 7) depicts the total capacity of general helicopter modes as a linear function of relief commodities and relief personnel.Equation (8) ensures that the number of relief personnel carried by special modes of helicopters does not exceed the total capacity of these modes of helicopters, while (9) makes the similar constraint on the amount of critical population evacuated from AAs. Equation (10) defines the fill rate, for a pair of AA and relief commodity or relief personnel, as a fraction of demand that is satisfied.Equation (11) states the fill rate with respect to evacuating critical population from AAs. Equation (12) shows the difference in fill rates between commodity  and   at AA .Equation (13) indicates that the amount of any relief material allocated to an AA should not exceed the demand level of this AA, whereas (14) depicts similar constraint with respect to the critical population.Equation (15) states the balance of flow on any helicopter mode at AAs. Equations ( 16)-( 19) are designed to connect scenariospecific decisions to scenario-unrelated decisions.Equation (16) guarantees that the inbound and outbound helicopter flows at any AA are balanced.Equation (17) ensures that the total flight duration of any helicopter mode does not exceed its available working time, allowing a helicopter return to any TF different from its departure TF.Equation (18) states that the relief commodities and personnel delivered from an open TF may not exceed its mobilization level at this TF, and (19) bounds the amount of critical population evacuated to any TF.
Finally, ( 20) and ( 21) define the appropriate domains for the scenario-unrelated decision variables and scenariospecific decision variables, respectively.

Solution Procedure
To the proposed biobjective stochastic programming problem, as the importance of objective ( 1) is known to be higher than that of objective (2) according to the performance metrics for humanitarian relief chain in Beamon and Balcik [22], the lexicographic optimization technique is applied here to transform the biobjective problem into single objective problems in a hierarchical order.As a result, the original model is transformed into two sequential models, (P1) and (P2), as follows: where (P2) is defined following the standard procedure of lexicographic optimization.However, as we can merely get a unique solution on the Pareto hyper surface (see Prats et al. [23] and the references therein) after solving (P1) and (P2) and the synergic requirement among related demands cannot be definitely satisfied by adding the constraint  1 ≥  * 1 to (P2), a new model modified from (P2), which is named as (P3) here utilized to substitute (P2), is defined as follows: where r  is the solution of    after solving (P1) and  is an adjustment coefficient to make a tradeoff between fill rate and cost.In addition, the fairness and synergy issues featured in the objective (1) can be inherited by adding the constraint    ≥ (1 − )r   to (P3), and the value of  1 can also be easily determined resorting to the solution of r  after solving (P3).
From (25) and ( 27) it can be seen that both (P1) and (P3) belong to mixed-integer stochastic programming models.For relatively small problem instances, it is reasonable to solve these two models directly using a commercial MIP solver like CPLEX.However, as the size of the problem instance increases, solving such size of a problem directly will become computationally unattractive, especially considering the stringent requirement on solving time in such an emergency situation.To make the model much more usable, two scenario-decomposition-based heuristic algorithms that can be used to solve (P1) and (P3), respectively, are proposed as follows.

Solution to (P1).
Since (P1) does not contain constraints related to expenditure, a scenario-decompositionbased heuristic approach was used to solve this problem.The details of this approach are as follows.
Step 1. Create a new problem named (P1  ) by duplicating    of   ,    of   , and    of   for all scenarios  ∈ Ξ on (P1).As a result, all decision variables in (P1  ) become scenario specific, meaning that (P1  ) can be decomposed into independent subproblems by scenario.
Step 2. Decompose (P1  ) into |Ξ| subproblems by scenario and solve each subproblem using a linear programming solver, denoting the value of    as x  , the value of    as ŝ  , and the value of    as ŷ  , respectively.Let   = max ∈Ξ { x  } and   = max ∈Ξ {ŝ   }.
After these assignments, it can be easily deduced that the constraint sets related to   and   in (P1  ) can be satisfied.However, to   , we cannot directly assign it a value in this way, otherwise some constraints, such as (6), will be violated.The next step is to give   a feasible solution.
Obviously, after executing this step, the constraint sets relating to   can be satisfied.By now, as the solutions of all scenario-unrelated variables in (P1) have been identified, (P1) can be decomposed into independent subproblems by scenario.The next step is to obtain the solution of    in (P1) with given values on all scenario-unrelated variables. Step where the number 1 in ( 28)-( 30) represents the first scenario.Equation ( 28) ensures that the decisions on whether establishing a facility at candidate location  under all scenarios are identical, that is,  1  =  2  = ⋅ ⋅ ⋅ =  |Ξ|  , and ( 29) and ( 30) state similar limitations on    and    , respectively.From ( 28)-( 30), it can be easily concluded that the two problems, (P3) and (P3  ), are equivalent.
(P3  ) can be decomposed into independent subproblems by scenario if constraint sets (28)- (30) are relaxed.Therefore, the Lagrangian relaxation of (P3  ) with respect to the constraint sets ( 28)-(30) can be described as follows: where  2 is the objective function of (P3), and the vectors u, w, and k in (31) are the Lagrangian multipliers associated with relaxed constraint sets (28), (29), and (30), respectively.The Lagrangian dual problem can be stated as follows: 4.2.2.Subgradient Algorithm.The lagrangian dual (32) is a convex nonsmooth program that can be solved by subgradient algorithm and the result can be taken as a lower bound (LB) on the optimal objective value of the original problem, (P3).A major advantage of this process is that it allows the problem to be split into separate subproblems for each scenario.The details of the subgradient algorithm are given below, in which the settings of all parameters employed were obtained using preliminary experimentation.
Record the values of   ,   , and   that have been utilized to generate the updated upper bound as x , ŝ , and ŷ , respectively.
where  is a user-defined parameter and the value  = 0.99 is used in the experiments.
if ((sum of squares = 0)          ( stop; else  =  + 1, goto Step 1 where the value of  2 is set to 200 and the value of  is set to 0.001. In Step 0, the Lagrangian multipliers are initialized with 0, and the lower and upper bounds of (P3) are initialized with −∞ and +∞, respectively.With the values of the Lagrangian multipliers assigned in Step 0, each subproblem is solved optimally using a linear programming solver, and the values of   ,   , and   are duly fixed in Step 1 if their solutions under all scenarios are identical at this intermediate iteration, and in this way the computation complexity can be effectively decreased in latter iterations.Step 2 computes new subgradient values and Step 3 updates the upper bound used in Steps 5 and 6.Here the upper bound is updated by using a heuristic according to the intermediate result to generate a primal feasible solution at intermediate iterations, and by updating UB the convergence speed can be accelerated [24].After updating the upper bound, a new step direction is found in Step 4 by taking the convex combinations of the last three subgradient values, and the typical oscillating behavior of the subgradient algorithm can be reduced adopting this combination operation [25,26].
The subgradient algorithm does not always generate primal feasible solutions to (P3  ) (which is equivalent to (P3)); that is, the values of    ,    , and    obtained from the subgradient algorithm may violate the constraint sets ( 28), (29), and (30), respectively.For this reason, a final heuristic is needed to restore feasibility starting from the dual feasible solution provided by the subgradient algorithm.

Feasibility Algorithm.
A primal feasible solution to (P3) can be obtained in a greedy manner starting from the values that generate the UB by the subgradient algorithm.The details of the feasibility heuristic are as follows.respectively, the objective value of (P3) equals

Test Case.
A case study focused on the disaster relief logistics of the Great Wenchuan Earthquake is here used to illustrate the utility of the temporary facility locationallocation model and verify the proposed heuristic decomposition algorithms.This great earthquake disaster occurred at the Long Menshan fault, a thrust structure located along the border of the Eurasian Plate and Indo-Australian Plate, and this fault consists of three parallel subfaults, that is, the front fault (Guanxian-Jiangyou fault), the central fault (Yingxiu-Beichuan fault), and the back fault (Wenxian-Maoxian fault), along the north-east direction.The epicenter was located near Wenchuan County at a depth of about 14 km.According to the above information and the discussion results with subject matter experts of Earthquake Agency of P. R. China, 15 disaster areas in the disaster region were chosen as AAs and 8 cities or towns, in which the inbound transportation condition by land remained relatively good was selected as candidate TFs around this region.The map of AAs and TFs is illustrated in Figure 1 where each AA is represented by a unique alphabet from "A" to "O" and each candidate TF is here marked by a unique number from "1" to "8, " and all the identified AAs were found to be aligned nearly in three rows corresponding to the three subfaults.8 types of demand, including water, food, medicine, doctors, nurses, rescuers, tents, and critical population, were considered in this case study.Table 1 gives the demand levels of AAs constructed from the loss assessment results of NDRC of P. R.China [27] and Yuan [28], in which the number of requisitions in depots H, K, L, and M is relatively higher than those of other AAs because these AAs are affected more seriously.Water, food, medicine, and tents were set in units of 15, 50, 60, and 45 kilograms, respectively.The data set in Table 1 here is represented as the first scenario in this case study, and other 4 scenarios, which represented more serious damage than scenario 1, were generated by multiplying a random coefficient between 1.0 and 1.5 with each element of the data set of scenario 1, respectively.Similarly, another 4 scenarios which represented lighter damage were produced using random coefficients between 0.5 and 1.0.The occurrence probabilities of these 9 scenarios were hypothetically set to 0.20, 0.15, 0.15, 0.15, 0.15, 0.05, 0.05, 0.05, and 0.05, respectively.
The parameters of candidate TFs are given in Tables 2 and 3.The maximum capacity of each relief supply at each candidate TF mainly depends on its original inventory level and its inbound transportation capacity.For example, candidate TF1 is a large urban city that not only has larger initial and potential capacities with respect to warehouse and health care facilities but also has more accessible inbound transportation capabilities, and thus its maximum capacities on all categories are larger than TF2 to TF5.In contrast, because the urban areas of TF6 to TF8 are relatively small and inbound capacities of these TFs are also relatively low, the maximum capacities of these TFs are set lower than the capacities of other TFs.Similarly, the differences on mobilization cost are also set among the 8 candidate TFs, and these differences mainly reflect the variations in inbound transportation costs.The fixed establishment costs of TF1-TF8 were set as 25,30,30,30,45,45,45, and 45 thousand US dollars, respectively, indicating that the cost of establishing TFs with good inbound transportation conditions (i.e., TF1, TF2, and TF3) was lower than these with limited inbound transportation conditions (i.e., TF6, TF7, and TF8).The setting of fixed cost of establishing a TF was not associated with its total mobilization levels of relief supplies as indicated in previous works (i.e., Özdamar et al. [29], Rawls and Turnquist [6], Salmerón and Apte [9], and Bozorgi-Amiri et al. [10]).This is because these open TFs in our model were used solely as temporary distribution depots for postdisaster relief and not as supply facilities for natural disaster asset prepositioning as in previous works.Table 4 shows a summary of the data used for helicopters.To the incorporated 5 modes of helicopters, Z-8S and M-170S were classified as special mode of helicopters and the remaining 3 were classified as general helicopters.The number of hours that each helicopter was available here included available flight time, not time required for maintenance, refueling, shift changes, or land operations.The average flight time between each TF and each AA was calculated based on the flight speed of each helicopter mode and the distance between each TF and each AA (i.e., Figure 1).
The adjustment coefficient on tradeoff, , was initialized as 0.01, and the weighting parameters of all relief supplies were set as follows: the total criticality weight of the three subsets of relief supplies (i.e., relief commodities, relief personnel, and hospital beds) was first set to 0.65, 0.25, and 0.1, respectively; then, to each type of relief supplies within the relief commodities subset and the relief personnel subset, its weighting parameter was set to 0.65 ∑ ∈ ℎ    / ∑ ∈  ∑ ∈ ℎ    (∀ ∈   ) and 0.25 ∑ ∈   / ∑ ∈  ∑ ∈   (∀ ∈   ), respectively.The penalty of discrepancy in fill rates for water and food and for medicine, doctors, and nurses was initially set as 0.001, denoting that there exists synergic requirement on distribution of food and water as well as on distribution of medicine, doctors, and nurses.The penalty of discrepancy among other relief supplies was set as zero.These values are useful for illustration purposes, but they may not match values that might be estimated by ultimate users of the model.

Computational Results
. The proposed decomposition algorithm was implemented using CPLEX 9.0 Callable Library embedded in C++ on a 2.4 GHz desktop computer with 4 GB of RAM.The parameter of EPGAP in CPLEX, which determines the relative tolerance on the MIP gap between the best integer objective and the objective of the best node remaining, was revised from default setting, 1.0-4, to 1.0-3 in the present case study in order to increase the convergence speed of computation.5  and 6 summarize the solutions of the scenario-unrelated decision variables.Table 5 shows that 6 TFs had opened and distributed widely around the disaster region.Among the 6 open TFs, only TF 8 was not arranged to receive and transship all types of relief supplies.Comparing with the maximum capacities for relief supplies depicted in Table 1, the As shown in Table 6, all modes of helicopters were allocated to TFs 1, 6, and 7.However, no helicopter was allocated to TFs 3, 4, and 8, and the total number of helicopters allocated to open TFs, such as zero helicopter allocated to TFs 3, 4, and 8, did not match the overall level of mobilized relief supplies of these TFs.This is because Table 6 only represents the initial deployment results of helicopters, and in the succedent distribution process, each helicopter was permitted to return to a depot other than the one from which it had departed.According to Tables 4 and 6, all available helicopters had been initially deployed to open TFs, indicating that helicopters may be a type of deficient resource to further improve the expected fill rate of AAs.The fill rates of all demand types at AAs under scenario 1 are shown in Table 7.As shown, to any type of demand, the fill rates from all the 15 AAs are equivalent.That is, the distribution of relief among AAs was fair in fill rate, indicating that the max-min method employed in (1) reaches its original goal.The fill rates among relief demands with any nonzero setting (    = 0.001) on the penalty of discrepancy (such as, the choosing relief supplies between food and water or among medicine, doctors, and nurses) are also equal.For example, the fill rates of water and food at any AA were equal to 51.6%, and the demands for medicine, doctors, and nurses from all the AAs also had the same fill rate value, 48.3%.In contrast, the fill rates among those relief demands that were not associated with any penalty of discrepancy, such as rescuers and tents, were different.This means the synergic requirement with respect to fill rate among correlated relief demands can be completely satisfied at a lower setting on the penalty coefficient of     , that is, 0.001.Table 8 depicts the fill rates on the 8 categories of demands and the integrated fill rate (which incorporates the criticality weight of each category of demands) under 3 representative scenarios: the levels of demand of scenario 1 are summarized in Table 1; and those of scenario 2 and scenario 6 individually represent a scenario with more severe damage than scenario 1 and a scenario with lighter damage, respectively.Regarding the 3 representative scenarios, the individual fill rate on any type of demands as well as the integrated fill rate under scenario 6 was the highest while the corresponding value under scenario 2 was the lowest.This is because the overall level of demand was highest under scenario 2 and lowest under scenario 6.As the integrated fill rates under the 3 representative scenarios were not high (i.e., 56.8%, 44.6%, and 77.0%under scenario 1, 2, and 6, resp.) but the mobilization levels of most relief supplies, at the same time, had not reach their maximum capacities at the open TFs (Table 2), this inconsistency was probably caused by the insufficiency of available helicopters, which has been analyzed before.Indeed, after adding 5 helicopters on each mode, the integrated fill rates of the 3 representative scenarios increased from 56.8%, 44.6%, and 77.0% to 87.1%, 72.0%, and 99%, respectively.This indicates that the proposed model can help disaster relief decision makers identify potential resource bottlenecks, such as the maximum capacities of relief supplies at open TFs, the overall distribution capacities of helicopters.In summary, to any type of demands whose fill rate must be improved, if the corresponding mobilization levels of all open TFs reach their maximum capacities, then the total capacities of open TFs should first be considered to be expanded, or more TFs should be established.Otherwise, more helicopters should be recruited.

Effects of the Penalty Coefficient 𝜔 𝐶
and the Adjustment Coefficient .As the synergic requirements among correlated relief supplies (i.e., medicine, doctors, and nurses) could be achieved after setting the penalty coefficient     = 0.001, to further investigate the influence of this penalty coefficient, we modified this coefficient's value while keeping the settings of other parameters unchanged.Figure 2 shows the results  Table 6: Initial deployment results of helicopter modes at TFs. with the modification of     , in which the correlation assumption among relief supplies was kept unchanged as aforementioned.The indices of the -axis in Figures 2(a) and 2(b), total discrepancy in fill rate and the expected fill rate on 3 representative scenarios, were calculated through respectively.    was modified by multiplying with 10 − , and the range of  is depicted by the -axis of Figures 2(a) and 2(b).According to Figure 2(a), there is no discrepancy in fill rate among correlated relief supplies when     is set as 0.001.As the value of penalty coefficient decreases further, the discrepancy among correlated relief supplies increases.When     falls below 1×10 −5 , the total discrepancy ceases to increase.However, as depicted in Figure 2(b), the system's expected fill rate does not change markedly as     varies.These phenomena can provide valuable information on setting the value of this penalty coefficient of discrepancy.
In order to arrive at an appropriate solution that would allow the decision maker to assess the tradeoff between the objective (1) and the objective (2), the aforementioned problem was solved several times while varying the adjustment coefficient, , and the results are illustrated in Figure 3.According to Figure 3(a), decreasing the value of  causes the values of the objective (1) and the objective (2) to increase synchronously.In other words, as the overall fill rate of  AAs increases, the procuring cost of relief supplies and the distribution cost will increase, and more TFs are needed to be established (the number of open TFs is depicted in Figure 3(b)).

Solution Quality and Evaluation of Computational
Speed.To evaluate the performance of the proposed decomposition-based heuristic algorithm, the problem described in Section 5.1 was solved using an ILOG CPLEX engine with the same settings (i.e., EPGAP).Tables 9  and 10 summarize the solutions of scenario-unrelated decision variables, and Table 11 depicts the fill rates under 3 representative scenarios.As shown in Table 9, the location results of open TFs are equivalent by utilizing the heuristic algorithm and ILOG CPLEX.Tables 5 and 9 show the presence of minor differences in mobilization levels of relief supplies between the two approaches.In contrast, there exist obvious differences with respect to the initial deployment result of helicopters at open TFs according to Tables 6  and 10, but, as shown in Tables 8 and 11, the expected fill rates under the 3 representative scenarios achieved by the heuristic algorithm are very similar to the corresponding fill rates computed through CPLEX.For example, the difference  is approximately 0.1% under each of the 3 scenarios.This is because the individual fill rates had been obtained by solving (P1) using the method described in Section 4.1, and in the procedure of that method, only the scenariounrelated variable   (∀ ∈ ,  ∈ ) was assigned with an approximate value.However,   denotes only the initial deployment results of all modes of helicopters at open TFs, and during the detailed distribution process helicopters were permitted to return to any open TF after finishing its current forward delivery mission.As a result, our proposed heuristic algorithm was able to produce approximate fill rates compared with these using CPLEX.The expected value of total cost computed using our proposed algorithm is 1.80% higher than the value gotten by CPLEX (31.302 million dollars versus 30.749 million dollars), but the computation time is 15.8% of that required by CPLEX (1696 seconds versus 10600 seconds) on the same computer.
To further investigate the performance of the proposed algorithm on the expected value of total cost and on computation time, the aforementioned problem was solved several times by varying the number of helicopters of each mode with the proposed algorithm and the CPLEX solver, respectively.The results are depicted in Figure 4, in which the numerical value of the -axis denotes the number of added helicopters on each mode.As shown in Figure 4(a), the expected value of total cost increased as the number of helicopters increased, and the total cost computed by the heuristic algorithm was slightly higher (from 1.28% to 2.59%) than the cost produced by the CPLEX solver.However, the computation time of the heuristic algorithm was much shorter than the time required by the CPLEX solver (the ratios between the two ranged from 16.2% to 25.7%), denoting that our heuristic algorithm exhibited a marked advantage in computation time than a commercial MIP solver.

Conclusions and Future Work
For the difficulty in predicting earthquakes, the mobilization of emergency resources is essential after disasters, especially when an earthquake disaster, for example, the Great Wenchuan Earthquake that occurred on May 12, 2008, in China, was devastating.We developed a disaster relief optimization model whose solution provides a strategic plan for locations of TFs, mobilization levels of relief supplies, initial deployment of helicopters at open TFs, and the distribution plans under each revealed scenario where the demands of AAs are presented by probabilistic scenarios.The model is a stochastic mixed-integer programming problem featuring two objectives: maximizing the expected minimal fill rate of affected areas where the mismatching distribution among correlated relief demands is penalized and minimizing the expected value of total cost.A lexicographic multiobjective optimization technique was used to transform the proposed model into two sequential single objective stochastic models.For relatively small data instances, the transformed models can be solved using a commercial mixed-integer solver such as CPLEX, but this approach does not scale well to large problems.As a result, scenario-decomposition-based heuristic algorithms had also been proposed to solve the transformed models, respectively, and this approach scales well, to allow large problems to be solved.
A numerical case, which is based on scenarios of the Great Wenchuan Earthquake of China, illustrates both the formulation and the character of the problem solutions in a practical context.The analysis of the results from the numerical case shows, for example, that the proposed model can help government officials to identify the bottlenecks hindering the improvement on the expected overall fill rate, such as the mismatching between the transportation capacity and the qualities of relief supplies to be distributed.The varied settings of the penalty coefficient,     , have obvious impact on the fill rate of correlated relief demands; when     reaches 0.001, the discrepancy in fill rate among correlated relief demands almost disappears.In addition, the adjustment coefficient, , can have good effect on making tradeoffs between the expected fill rate and the expected value of total cost.Across a test problem, the proposed algorithm consistently produces solutions whose expected overall fill rates are almost the same as the optimal values and whose total costs remained within 2.59% higher than the costs of the optimal, with computation times that are between 16.2% and 25.7% of the time required using commercial MIP software.These experiments indicate that the proposed algorithm can be utilized as a large-scale disaster relief logistics planning tool.
The character of the discrete distribution in probability for uncertain demand also implies a desirable direction for further enhancing the applicability of the model.That is, as the forecasting results on uncertain demands describing as scenarios may be inaccurate early after a disaster, the postdisaster conditions could become more and more clear over time.Meanwhile, considering that the mobilization campaign itself requires a period of time to be completely realized, we have opportunities to modify the original Mathematical Problems in Engineering mobilization plan according to updated information and the actual progress of the relief logistics.This offers an important avenue for enhancement of the current work by establishing dynamic facility location and allocation models.
Another desirable direction for future work is to further enhance the efficiency of the proposed heuristic algorithm.For example, after decomposing the original problem into subproblems by scenario, it may be worthy to specially design heuristic algorithms to solve each subproblem, such as separating the forward distribution problem as commodity flow subproblem and vehicle flow subproblem as in Haghani and Oh [30] with given values on all scenario-unrelated decision variables or directly proposing heuristic algorithms to solve each facility location and allocation subproblem within each scenario, in order to further reduce the time required for computation.

Figure 1 :
Figure 1: Map of candidate TFs and AAs.

Figure 2 :
Figure 2: Total discrepancy and expected fill rate with varied     .

Figure 3 :
Figure 3: Expected value of total cost and the number of open TFs for different values of .

Figure 4 :
Figure 4: Performance comparing results between the heuristic algorithm and the CPLEX solver.
Scenario-Specific Parameters, All under Scenario    : Demand level for relief supply  ( ∈   ∪   ) at area  (units) : Maximum mobilization levels for relief supply  at TF  (units) ℎ  : Average heft of relief supply  ( ∈   , kilograms/units)   : Total number of available helicopters of mode    : Average capacity for relief supplies of helicopter mode  (kilograms/helicopter × trip)   : Average capacity for relief personnel of helicopter mode  (persons/helicopter × trip)   : Average capacity for critical population of helicopter mode  (persons/helicopter × trip)   : Average working time of helicopter mode  in the target responsive period (hours)    : Fixed cost of establishing TF  ($)    : Fixed cost of hiring a unit of helicopter mode  ($)    : Operating cost for a unit of helicopter mode  ($/hour)   : Mobilization cost for a unit of relief supply  at TF  ($/unit)   : Average flight time between TF  and AA  by helicopter mode  (hours)   : Weighting parameter for relief supply  in fill rate     : Weighting parameter on the penalty of discrepancy in fill rate for relief supplies  and   : A large number. : Amount of critical population at area  (persons)   : Probability of scenario  occurrences.Scenario-Unrelated Decision Variables   : Whether TF  is open (1 TF  open, 0 TF  closed)   : Mobilization level of relief supply  at TF  (units)   : Number of helicopter mode  that initially allocated to TF  (helicopters).Scenario-Specific Decision Variables (Units), All under Scenario     : Amount of relief supply  ( ∈   ∪   ) delivered from TF  to AA  by helicopter mode     : Amount of critical population picked up by helicopter  from AA  to TF     : Number of trips from TF  to AA  by helicopter mode     : Number of trips from AA  to TF  by helicopter mode     : Fill rate of relief supply  at AA  Δ    : Fill rate difference between relief supplies  and   at AA .

Table 2 :
Maximum capacities of relief supplies at TFs (units).

Table 4 :
Summary data of helicopter modes.

Table 5 :
Mobilization levels of relief supplies at TFs (units).

Table 7 :
Fill rates of all demands at AAs under scenario 1 (percent).

Table 8 :
Fill rates of all demands under three representative scenarios (percent).