Hindawi Publishing Corporation

Mathematical Problems in Engineering

Volume 2014, Article ID 939853, 17 pages
http://dx.doi.org/10.1155/2014/939853

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

Yajie Liu and Bo Guo

College of Information System and Management, National University of Defense Technology, Changsha, Hunan 410073, China

Correspondence should be addressed to Yajie Liu; [email protected]

Received 6 June 2013; Revised 2 September 2013; Accepted 9 September 2013; Published 16 January 2014

Academic Editor: Lu Zhen

Copyright © 2014 Y. Liu and B. Guo. This is an open access article distributed under the Creative Commons Attribution License,
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 proposed model: maximizing
the expected minimal fill rate of affected areas, where the mismatching distribution among correlated relief demands is penalized,
and minimizing 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.

1. Introduction in mountainous areas, where the transportation infrastruc-

ture, such as roads, bridges and railway lines, can suffer severe
Earthquakes are a special type of natural disaster that can damage after large earthquakes. The following three factors
result in huge catastrophes. The Great Sichuan Earthquake, must be taken into account. First, the demand for disaster
which occurred on May 12, 2008, in Wenchuan County, reliefs in scattered affected areas (AAs) can be uncertain
Sichuan Province, China, was one such case. It had a recorded early after earthquake [1]. This is because collecting actual
magnitude of 8.0 and imposed tremendous suffering on the and real-time information from AAs is relatively infeasible
local residents, causing over 69,000 deaths, 18,000 missing due to power disruption, damage to communication equip-
persons cases, 374,000 injuries, and huge loss of property. ment, and severe telecommunication traffic congestion in the
Although thousands of networked seismograph stations have disaster region. Second, as the transportation infrastructure
been installed around the world and although these stations is usually seriously damaged and repairing the damaged
use powerful computers, it is still difficult to predict when roadways quickly is very difficult within the mountainous
and where an earthquake will strike. Therefore, an effective regions, traditional transportation tools, such as trucks and
and necessary approach coping with an earthquake disaster is trains, are impractical. Instead, helicopters become the most
to plan and implement a disaster relief campaign that would useful transportation means within the disaster regions
reduce the damage right after the earthquake disaster takes during the initial postdisaster period (i.e., the first 72 hours
place. after a disaster, a golden relief period that is critical to
Emergency logistics plays a vital role in disaster relief provide relief since communities are not expected to stand
campaigns. However, the planning and implementation of on their own for much longer than that time). In fact, it
relief logistics are challenging, especially when help is needed has been reported that 99 helicopters were mobilized for
2 Mathematical Problems in Engineering

rescue missions after the Great Sichuan Earthquake [2]. The 2. Literature Survey
last factor is the optimization objectives of relief logistics.
Conflicting objectives, such as maximizing the fill rate of Relief logistics is essential to the timely and effective provision
AAs and minimizing costs, should be contemplated con- of resources to aid people affected by natural disasters.
currently. Although in this emergency situation minimizing To shorten the response time of disaster relief, a natural
costs should not be emphasized in planning the logistics, the approach is to establish facilities and stock prepositioning
excess accumulation of supplies can be effectively avoided by before the strikes of disasters. However, in making decisions
incorporating this objective. Otherwise, overaccumulation of for this strategic planning before the disaster actually takes
resources may delay or disrupt the activities of relief logistics. place, uncertainties, such as the degree of demands and trans-
In addition, because earthquakes are more unpredictable portation conditions, are the key challenges. The significance
than other types of natural disasters, such as floods and of these uncertainties has motivated a number of researchers
typhoons, it is usually uneconomical or even impossible to to address stochastic optimization in disaster relief planning
maintain enough capacity for relief operations before the using probabilistic scenarios representing disasters and their
occurrence of an earthquake. As a result, temporary facilities outcomes.
(TFs) must be established as distribution depots at the edges Considering the uncertainty of demands from potential
of disaster regions in order to promote these relief supplies disaster areas, Balcik and Beamon [3] developed a maximally
flowing into the affected areas in a controllable and fluid way. covered location model of prepositioning relief commodities
The purpose of this paper is to model the strategic to determine the number, locations and capacities of the
planning of relief logistics in response to a devastating relief distribution centers, in which the pre- and postdis-
earthquake disaster in a mountainous region. It addresses aster budget constraints were specially considered. Chang
specific issues such as the inherent uncertainty regarding et al. [4] presented a facility location-allocation model for
the degree of demand, conflicting objectives, and helicopter- stocking and distributing rescue resources taking possible
based distribution plans. The main contributions of this paper flood scenarios into account, and this model can be utilized
can be summarized as follows. as a decision-making tool in generating plans for flood
emergency logistics.
(i) Development of a novel optimization model to Beside the uncertainty on demand, some studies fur-
address the disaster relief logistics problem as a mul- ther considered other uncertainties, such as survival of
tiobjective, stochastic, mixed-integer, nonlinear pro- prepositioned stocks, transportation network conditions, and
gramming problem: this model combines temporary distribution costs. Ukkusuri and Yushimito [5] presented
facility locations, mobilization levels of relief supplies, a model to identify locations for prepositioned supplies,
and initial deployment of helicopters with uncertainty and its objective was to maximize the probability that
about demands. It also includes helicopter-based demand points could be reached from at least one supply
delivery plans for relief supplies and evacuation plans facility under conditions where facility locations and links
for critical population (referring to those in need of in transportation network may be destroyed. Rawls and
emergency medical evacuation) under each revealed Turnquist [6] developed a two-stage stochastic mixed-integer
scenario. The formulation not only takes into account program that determined the locations and quantities of
the expected total cost of relief logistics but also various types of emergency commodities, in which uncer-
considers the fill rate of AAs and the synergy among tainties of demand and transportation network availability
correlated relief demands and the fairness among AAs were considered. Rawls and Turnquist [7] also extended
in distribution process. this model with additional service quality constraints that
(ii) Proposal of a decomposition-based heuristic algo- ensured the probability of meeting all demands was at least
rithm based on the structure of the proposed model: a preset level, and the distance between a supply-demand
as a result, the stringent time requirements for solving pair was less than a specific limit. Mete and Zabinsky
the model can be satisfied in such an emergency [8] proposed a stochastic optimization model for disaster
environment of disaster response. preparedness and response with uncertainties on demand
(iii) Application of this model to a real-world disaster and cost, in order to assist in deciding the location and
relief chain, specially the “5.12” Wenchuan Earth- allocation of medical supplies to be used during emergencies
quake relief chain, which was dedicated to the relief in the Seattle area, which is known to be vulnerable to
supply and medical evacuation after the disaster. earthquakes. Salmerón and Apte [9] presented a two-stage
stochastic optimization model for planning the allocation
The rest of this paper is organized as follows. Section 2 of budget on acquiring and positioning relief assets; in this
provides a review of disaster relief optimization methods model the first-stage decisions included the expansion of
related to the present work. In Section 3, a multiobjective, resources such as warehouses, medical facilities, ramp spaces,
stochastic, mixed-integer, nonlinear programming disaster and shelters, whereas the second-stage decisions focused
relief model is described. Then, we develop a decomposition- on the transportation plans for commodities and critical
based heuristic algorithm suitable for solving the model in population; this model also contemplated the transportation
Section 4. This is followed by computational results of model network availability following a disaster under uncertainties
testing in Section 5. Conclusions and directions for future of demand and cost. Bozorgi-Amiri et al. [10] developed
research are provided in the last section. a multiobjective robust stochastic programming model for
Mathematical Problems in Engineering 3

disaster relief logistics, in which demands, supplies, and the stochastic model targets the strategic planning of postdisaster
costs of procurement and transportation were considered as relief logistics, there exists stringent requirement on time
uncertain parameters: it attempted to minimize the sum of for solving this model in such an emergency situation. As
the expected value of total costs and maximize the affected a result, we specially developed a scenario-decomposition-
areas’ satisfaction levels. Döyen et al. [11] developed a two- based heuristic algorithm to solve the proposed model, in
stage stochastic programming model for humanitarian relief order to satisfy the stringent requirement on making an
logistics problem where decisions are made for pre- and emergency decision in a short time.
postdisaster rescue centers, the amount of relief items to In conclusion, our contributions or differences compared
be stocked at the predisaster rescue centers, the amount of with related existing studies can be concluded as two aspects.
relief items flows at each echelon, and the shortage degree of The first is the novelty of the proposed biobjective stochastic
relief items with uncertainties on demand and transportation programming model itself. Most current stochastic emer-
time/cost. gency facility location and facility location-allocation models
Unlike these previous works, which gave extra con- were designed for predisaster preparedness and response
sideration to uncertainties other than demand, such as with single objective, such as minimizing expected cost [6,
supplies and transportation conditions, this paper merely 7, 11], maximizing the coverage of potential disaster areas
considers the uncertainty of demand. This is because our [3], or maximizing the probability on reaching the demand
work mainly focuses on the strategic planning of relief points [5], while our postdisaster stochastic facility location-
logistics that an earthquake disaster has struck, and under allocation model has two objectives on expected fill rate and
these conditions, the amounts of supplies available tend to cost. Although the model proposed by Bozorgi-Amiri et al.
be relatively definite. (However, there are also cases where [10] has similar objectives as ours, its two robust criteria, the
the availability of supplies is unpredictable even after an variance of total cost and satisfaction level among scenarios,
earthquake disaster, especially if it depends on the willingness and the solutions infeasibility due to parameter uncertainty,
of private international donors to financially support the were incorporated in the objectives of their model, whereas
relief operations. In most situations, the assumption that the the synergic issue among related demands was extraconsid-
availability of supplies is relatively definite after a disaster is ered in the first objective of our model. In addition, the coor-
reasonable.) In addition, because helicopters are utilized as dination between the mobilization levels of transportation
the preferred unique emergency vehicle transportation tools tools (i.e., helicopters) and the amount of relief supplies were
in this model and because the physical conditions of roads especially considered in our work for postdisaster relief logis-
have no significant effect on the ability of helicopters to travel, tics planning, and this coordination had not been contem-
the uncertainty of transportation conditions was not taken plated in the aforementioned literature because they mainly
into account here. focused on preparations to be made before disasters. The
Most existing stochastic models on prepositioning of second is about the solution approach to solve the proposed
emergency supplies for disaster response belong to two- biobjective model. An approximate lexicographic approach
stage mixed-integer stochastic programming problems, and was applied to our biobjective model, while the technique
decomposition algorithms can be adopt to solve them, such of Compromise Programming (CP) was adopted in Bozorgi-
as the L-shaped algorithm based on Bender decomposition Amiri et al. [10], in which identifying the value of the weight
in Birge and Louveaux [12] and the dual decomposition coefficient 𝜔 was not easy. Two scenario-decomposition-
algorithm based on scenario decomposition and Lagrangian based heuristic algorithms were further developed for the
relaxation in Carøe and Schultz [13]. Besides the above two transformed single objective stochastic programming
two classical decomposition algorithms, some other decom- problems, respectively, to satisfy the stringent requirement
position algorithms have also been developed for mixed- on time for solving the problem, and in contrast, there is no
integer stochastic programs with special structures: stochastic similar work on algorithm design in Bozorgi-Amiri et al. [10].
program with purely binary first-stage decision variables and Moreover, the problem decomposition technique utilized in
arbitrary second-stage decision variables [14], pure continu- this paper is scenario decomposition, whereas in Rawls and
ous or pure binary first-stage variables and 0-1 mixed-integer Turnquist [6] it is Benders decomposition and in Döyen
recourse variables [13], mixed-integer first-stage variables et al. [11] is by relaxing special constraints to obtain two
and pure integer second-stage variables [15], pure binary first- subproblems. In other words, the proposed solution methods
stage variables and 0-1 mixed-integer recourse variables [16– for our stochastic optimization problem based on scenario-
18], and 0-1 mixed-integer variables for both the first and decomposition heuristic algorithm within an approximate
the second-stage variables [19, 20]. However, as to the best lexicographic optimization frame are novel compared to
of our knowledge, only a few studies have been performed existing work for emergency relief logistics planning.
on the development of decomposition methods suitable for
humanitarian logistics problems [6, 11, 21]. Although the
computation complexity of most facility location models
3. Facility Location-Allocation Model
preparing for disasters is high when the sizes of optimization 3.1. Notation Descriptions
problems are large, the time required on solving such a model
is relatively insignificant because these models solely focus Set and Indices
on strategic planning of relief logistics for potential disasters
that have not yet stuck. In contrast, as our mixed-integer 𝐼 : Set of candidate TFs; 𝑖 ∈ 𝐼
4 Mathematical Problems in Engineering

𝐽: Set of AAs; 𝑗 ∈ 𝐽 Scenario-Specific Parameters, All under Scenario 𝜉

𝐾: Set of relief supplies; 𝑘 ∈ 𝐾 𝜉
𝑑𝑗𝑘 : Demand level for relief supply 𝑘 (𝑘 ∈ 𝐾𝑅 ∪ 𝐾𝑆 ) at
𝑅
𝐾 : Subset of relief supplies that can be classified as area 𝑗 (units)
relief commodities, such as water, food, medicine and
tent 𝑑𝑐𝑗𝜉 : Amount of critical population at area 𝑗 (persons)
𝐾𝑆 : Subset of relief supplies that can be classified as
𝑝𝜉 : Probability of scenario 𝜉 occurrences.
relief personnel, such as health personnel and relief
workers
𝐾𝐶: Subset of relief supplies that accommodates Scenario-Unrelated Decision Variables
critical populations, that is, hospital beds
𝑥𝑖 : Whether TF 𝑖 is open (1 TF 𝑖 open, 0 TF 𝑖 closed)
𝑀: Set of helicopter modes; 𝑚 ∈ 𝑀
𝑠𝑖𝑘 : Mobilization level of relief supply 𝑘 at TF 𝑖 (units)
𝑀𝐺: Subset of general helicopter modes (utilized as
delivering relief commodities, relief personnel, and 𝑦𝑖𝑚 : Number of helicopter mode 𝑚 that initially
critical population); 𝑚 ∈ 𝑀𝐺 allocated to TF 𝑖 (helicopters).

𝑀𝑆 : Subset of special helicopter modes (utilized as

transporting relief personnel and critical population); Scenario-Specific Decision Variables (Units), All under
𝑚 ∈ 𝑀𝑆 Scenario 𝜉
Ξ: Set of scenarios; 𝜉 ∈ Ξ. 𝜉
𝑧𝑖𝑗𝑘𝑚 : Amount of relief supply 𝑘 (𝑘 ∈ 𝐾𝑅 ∪ 𝐾𝑆 )
delivered from TF 𝑖 to AA 𝑗 by helicopter mode 𝑚
Deterministic Parameters 𝜉
𝑧𝑐𝑗𝑖𝑚 : Amount of critical population picked up by
𝑈𝑖𝑘 : Maximum mobilization levels for relief supply 𝑘 helicopter 𝑚 from AA 𝑗 to TF 𝑖
at TF 𝑖 (units) 𝜉
𝑅
𝑛𝑓𝑖𝑗𝑚 : Number of trips from TF 𝑖 to AA 𝑗 by helicopter
ℎ𝑘 : Average heft of relief supply 𝑘 (𝑘 ∈ 𝐾 , kilo- mode 𝑚
grams/units)
𝜉
𝑛𝑏𝑗𝑖𝑚 : Number of trips from AA 𝑗 to TF 𝑖 by helicopter
𝐴 𝑚 : Total number of available helicopters of mode 𝑚
mode 𝑚
𝜃𝑚 : Average capacity for relief supplies of helicopter 𝜉
𝑟𝑗𝑘 : Fill rate of relief supply 𝑘 at AA 𝑗
mode 𝑚 (kilograms/helicopter × trip)
𝜉
𝜂𝑚 : Average capacity for relief personnel of helicopter Δ𝑟𝑗𝑘𝑘󸀠 : Fill rate difference between relief supplies 𝑘
mode 𝑚 (persons/helicopter × trip) and 𝑘󸀠 at AA 𝑗.
𝜗𝑚 : Average capacity for critical population of heli-
copter mode 𝑚 (persons/helicopter × trip) 3.2. Mathematic Formulation. The model is constructed
𝑇𝑚 : Average working time of helicopter mode 𝑚 in the with two objectives based on the stochastic programming
target responsive period (hours) approach, in which uncertainty is represented by a set of finite
discrete scenarios:
𝑓𝑖𝐹 : Fixed cost of establishing TF 𝑖 ($)
𝑓𝑚𝐻: Fixed cost of hiring a unit of helicopter mode 𝑚 𝜉
Max 𝑓1 = ∑ 𝑝𝜉 ( ∑ 𝑤𝑘 (min {𝑟𝑗𝑘 }))
𝑗∈𝐽
($) 𝜉∈Ξ 𝑘∈𝐾
(1)
𝑐𝑚𝐻: Operating cost for a unit of helicopter mode 𝑚 𝜉 𝐶 𝜉
−∑𝑝 ∑ ∑ ∑ 𝜔𝑘𝑘󸀠 Δ𝑟
𝑗𝑘𝑘󸀠
,
($/hour) 𝑗∈𝐽 𝑘∈𝐾 𝑘󸀠 ∈𝐾
𝜉∈Ξ
𝑐𝑖𝑘 : Mobilization cost for a unit of relief supply 𝑘 at TF
𝑖 ($/unit) Min 𝑓2 = ∑ 𝑓𝑖𝐹 𝑥𝑖 + ∑ ∑ 𝑐𝑖𝑘 𝑠𝑖𝑘
𝑖∈𝐼 𝑖∈𝐼 𝑘∈𝐾
𝑡𝑖𝑗𝑚 : Average flight time between TF 𝑖 and AA 𝑗 by
helicopter mode 𝑚 (hours) + ∑ ∑ 𝑓𝑚𝐻𝑦𝑖𝑚
𝑖∈𝐼 𝑚∈𝑀
𝑤𝑘 : Weighting parameter for relief supply 𝑘 in fill rate (2)
𝐶
𝜔𝑘𝑘󸀠:Weighting parameter on the penalty of discrep- + ∑ ∑ ∑ ∑ 𝑝𝜉 𝑐𝑚𝐻𝑡𝑖𝑗𝑚
𝜉∈Ξ 𝑖∈𝐼 𝑗∈𝐽 𝑚∈𝑀
ancy in fill rate for relief supplies 𝑘 and 𝑘󸀠
𝜉 𝜉
𝐿: A large number. × (𝑛𝑓𝑖𝑗𝑚 + 𝑛𝑏𝑗𝑖𝑚 ),
Mathematical Problems in Engineering 5

𝜉
s.t. Δ𝑟𝑗𝑘𝑘󸀠 ≥ 0, ∀𝑖 ∈ 𝐼, 𝑗 ∈ 𝐽, 𝑘 ∈ 𝐾, 𝑘󸀠 ∈ 𝐾, 𝑚 ∈ 𝑀, 𝜉 ∈ Ξ.

𝑠𝑖𝑘 ≤ 𝐿𝑥𝑖 , ∀𝑖 ∈ 𝐼, 𝑘 ∈ 𝐾, (3) (21)

𝑠𝑖𝑘 ≤ 𝑈𝑖𝑘 , ∀𝑖 ∈ 𝐼, 𝑘 ∈ 𝐾, (4) The two objectives are represented by (1) and (2), respectively.
They are explained as follows.
𝑦𝑖𝑚 ≤ 𝐴 𝑚 𝑥𝑖 , ∀𝑖 ∈ 𝐼, 𝑚 ∈ 𝑀, (5)
Objective 1. The first objective maximizes the expected mini-
∑ 𝑦𝑖𝑚 ≤ 𝐴 𝑚 , ∀𝑚 ∈ 𝑀, (6)
mal fill rate of affected areas while penalizing the mismatch-
𝑖∈𝐼 ing 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
𝜃𝑚 𝜂𝑚 (7) 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
𝐶
𝜉
∑ 𝑧𝑖𝑗𝑘𝑚 𝜉
≤ 𝜂𝑚 𝑛𝑓𝑖𝑗𝑚 , ∀𝑖 ∈ 𝐼, 𝑗 ∈ 𝐽, 𝑚 ∈ 𝑀𝑆 , 𝜉 ∈ Ξ, (8) 𝑤𝑘𝑘 󸀠 is set larger if the synergistic requirement between 𝑘 and
󸀠
𝑘∈𝐾𝑆 𝑘 is higher.
Since this objective is explicitly nonlinear, the linear
𝜉 𝜉
𝑧𝑐𝑗𝑖𝑚 ≤ 𝜗𝑚 𝑛𝑏𝑗𝑖𝑚 , ∀𝑖 ∈ 𝐼, 𝑗 ∈ 𝐽, 𝑚 ∈ 𝑀, 𝜉 ∈ Ξ, (9) equivalent equations can be rewritten with the help of an
auxiliary variable 𝑅𝑘𝜉 as follows:
𝜉
𝜉
𝑧𝑖𝑗𝑘𝑚
𝑟𝑗𝑘 =∑ ∑ 𝜉
, ∀𝑗 ∈ 𝐽, 𝑘 ∈ 𝐾𝑅 ∪ 𝐾𝑆 , 𝜉 ∈ Ξ, (10) Max 𝑓1 = ∑ 𝑝𝜉 ∑ 𝑤𝑘 𝑅𝑘𝜉
𝑖∈𝐼 𝑚∈𝑀 𝑑𝑗𝑘 𝜉∈Ξ 𝑘∈𝐾
(22)
𝜉 𝜉
𝜉
𝑧𝑐𝑗𝑖𝑚 − ∑ 𝑝𝜉 ∑ ∑ ∑ 𝑤𝑘𝑘
𝐶
󸀠 Δ𝑟 ,
𝑟𝑗𝑘 =∑ ∑ , ∀𝑗 ∈ 𝐽, 𝑘 ∈ 𝐾𝐶, 𝜉 ∈ Ξ, (11) 𝜉∈Ξ
𝑗𝑘𝑘󸀠
𝑗∈𝐽 𝑘∈𝐾 𝑘󸀠 ∈𝐾
𝑖∈𝐼 𝑚∈𝑀 𝑑𝑐𝑗𝜉
𝜉 𝜉 𝜉 s.t. 𝑅𝑘𝜉 ≤ 𝑟𝑗𝑘
𝜉
, ∀𝑗 ∈ 𝐽, 𝑘 ∈ 𝐾, 𝜉 ∈ Ξ, (23)
Δ𝑟𝑗𝑘𝑘󸀠 ≥ 𝑟𝑗𝑘 − 𝑟𝑗𝑘󸀠 , ∀𝑗 ∈ 𝐽, 𝑘 ∈ 𝐾, 𝑘󸀠 ∈ 𝐾, 𝜉 ∈ Ξ, (12)
𝜉 𝜉
𝑅𝑘𝜉 ≥ 0, ∀𝑘 ∈ 𝐾, 𝜉 ∈ Ξ. (24)
∑ ∑ 𝑧𝑖𝑗𝑘𝑚 ≤ 𝑑𝑗𝑘 , ∀𝑗 ∈ 𝐽, 𝑘 ∈ 𝐾𝑅 ∪ 𝐾𝑆 , 𝜉 ∈ Ξ, (13)
𝑖∈𝐼 𝑚∈𝑀 Objective 2. The second objective function minimizes the sum
𝜉
of four types of costs: the cost of establishing facilities, the cost
∑ ∑ 𝑧𝑐𝑗𝑖𝑚 ≤ 𝑑𝑐𝑗𝜉 , ∀𝑗 ∈ 𝐽, 𝜉 ∈ Ξ, (14) of procuring relief supplies, the cost of recruiting helicopters,
𝑖∈𝐼 𝑚∈𝑀 and the expected cost of transportation.
𝜉 𝜉 Equations (3)–(6) define all the constraints of the
∑ 𝑛𝑏𝑗𝑖𝑚 = ∑ 𝑛𝑓𝑖𝑗𝑚 , ∀𝑗 ∈ 𝐽, 𝑚 ∈ 𝑀, 𝜉 ∈ Ξ, (15) 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
∑ 𝑛𝑓𝑖𝑗𝑚 ≤ 𝑦𝑖𝑚 + ∑ 𝑛𝑏𝑗𝑖𝑚 , ∀𝑖 ∈ 𝑖, 𝑚 ∈ 𝑀, 𝜉 ∈ Ξ, (16)
𝑗∈𝐽 𝑗∈𝐽 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
(17) 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
∑ ∑ 𝑧𝑖𝑗𝑘𝑚 ≤ 𝑠𝑖𝑘 , ∀𝑖 ∈ 𝐼, 𝑘 ∈ 𝐾𝑅 ∪ 𝐾𝑆 , 𝜉 ∈ Ξ, (18) 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
(19) these modes of helicopters, while (9) makes the similar
𝑗∈𝐽 𝑚∈𝑀
constraint on the amount of critical population evacuated
𝑥𝑖 ∈ {0, 1} , 𝑠𝑖𝑘 ≥ 0, 𝑦𝑖𝑚 ≥ 0 integer, from AAs. Equation (10) defines the fill rate, for a pair of
(20) 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
𝑧𝑖𝑗𝑘𝑚 ≥ 0, 𝑧𝑐𝑗𝑖𝑚 ≥ 0, 𝑛𝑓𝑖𝑗𝑚 ≥ 0 integer, (12) shows the difference in fill rates between commodity
𝑘 and 𝑘󸀠 at AA 𝑗. Equation (13) indicates that the amount
𝜉 𝜉
𝑛𝑏𝑗𝑖𝑚 ≥ 0 integer, 0 ≤ 𝑟𝑗𝑘 ≤ 1, of any relief material allocated to an AA should not exceed
6 Mathematical Problems in Engineering

the demand level of this AA, whereas (14) depicts similar cost. In addition, the fairness and synergy issues featured in
constraint with respect to the critical population. Equation the objective (1) can be inherited by adding the constraint
𝜉 𝜉
(15) states the balance of flow on any helicopter mode at AAs. 𝑟𝑗𝑘 ≥ (1 − 𝛼)̃𝑟𝑗𝑘 to (P3), and the value of 𝑓1 can also be
Equations (16)–(19) are designed to connect scenario- 𝜉
specific decisions to scenario-unrelated decisions. Equation easily determined resorting to the solution of 𝑟̃𝑗𝑘 after solving
(16) guarantees that the inbound and outbound helicopter (P3).
flows at any AA are balanced. Equation (17) ensures that the From (25) and (27) it can be seen that both (P1) and (P3)
total flight duration of any helicopter mode does not exceed belong to mixed-integer stochastic programming models.
its available working time, allowing a helicopter return to any For relatively small problem instances, it is reasonable to
TF different from its departure TF. Equation (18) states that solve these two models directly using a commercial MIP
the relief commodities and personnel delivered from an open solver like CPLEX. However, as the size of the problem
TF may not exceed its mobilization level at this TF, and (19) instance increases, solving such size of a problem directly will
bounds the amount of critical population evacuated to any become computationally unattractive, especially considering
TF. the stringent requirement on solving time in such an emer-
Finally, (20) and (21) define the appropriate domains gency situation. To make the model much more usable, two
for the scenario-unrelated decision variables and scenario- scenario-decomposition-based heuristic algorithms that can
specific decision variables, respectively. be used to solve (P1) and (P3), respectively, are proposed as
follows.
4. Solution Procedure
4.1. Solution to (P1). Since (P1) does not contain con-
To the proposed biobjective stochastic programming prob- straints related to expenditure, a scenario-decomposition-
lem, as the importance of objective (1) is known to be based heuristic approach was used to solve this problem. The
higher than that of objective (2) according to the performance details of this approach are as follows.
metrics for humanitarian relief chain in Beamon and Balcik
[22], the lexicographic optimization technique is applied here Step 1. Create a new problem named (P1󸀠 ) by duplicating 𝑥𝑖𝜉
to transform the biobjective problem into single objective 𝜉 𝜉
problems in a hierarchical order. As a result, the original of 𝑥𝑖 , 𝑠𝑖𝑘 of 𝑠𝑖𝑘 , and 𝑦𝑖𝑚 of 𝑦𝑖𝑚 for all scenarios 𝜉 ∈ Ξ on (P1).
model is transformed into two sequential models, (P1) and As a result, all decision variables in (P1󸀠 ) become scenario
(P2), as follows: specific, meaning that (P1󸀠 ) can be decomposed into inde-
pendent subproblems by scenario.
(P1) : 𝑓1∗ = Max 𝑓1
(25) Step 2. Decompose (P1󸀠 ) into |Ξ| subproblems by scenario
(3) – (21) , and solve each subproblem using a linear programming
s.t. : {
(23) - (24) ,
solver, denoting the value of 𝑥𝑖𝜉 as 𝑥̂𝑖𝜉 , the value of 𝑠𝑖𝑘
𝜉 𝜉
as 𝑠̂𝑖𝑘 ,
𝜉 𝜉 𝜉
(P2) : 𝑓2∗ = Min 𝑓2 and the value of 𝑦𝑖𝑚 as 𝑦̂𝑖𝑚 , respectively. Let 𝑥𝑖 = max𝜉∈Ξ {𝑥̂𝑖 }
𝜉
and 𝑠𝑖𝑘 = max𝜉∈Ξ {̂𝑠𝑖𝑘 }.
{(3) – (21) , (26)
s.t. : {(23) - (24) , After these assignments, it can be easily deduced that the
∗
{𝑓1 ≥ 𝑓1 , constraint sets related to 𝑥𝑖 and 𝑠𝑖𝑘 in (P1󸀠 ) can be satisfied.
However, to 𝑦𝑖𝑚 , we cannot directly assign it a value in this
where (P2) is defined following the standard procedure of
way, otherwise some constraints, such as (6), will be violated.
lexicographic optimization. However, as we can merely get
The next step is to give 𝑦𝑖𝑚 a feasible solution.
a unique solution on the Pareto hyper surface (see Prats et
al. [23] and the references therein) after solving (P1) and (P2) 𝜉
Step 3. For each helicopter mode 𝑚, let 𝑦̃𝑖𝑚 = ∑𝜉 𝑝𝜉 𝑦̂𝑖𝑚 (∀𝑖 ∈
and the synergic requirement among related demands cannot
𝐼), 𝑦̂𝑖𝑚 = ⌊𝑦̃𝑖𝑚 ⌋ (∀𝑖 ∈ 𝐼), and Δ𝑦𝑖𝑚 = 𝑦̃𝑖𝑚 − 𝑦̂𝑖𝑚 (∀𝑖 ∈ 𝐼),
be definitely satisfied by adding the constraint 𝑓1 ≥ 𝑓1∗ to
then sort Δ𝑦𝑖𝑚 in nonincreasing order according to their
(P2), a new model modified from (P2), which is named as
values, and move the first elements with the number of
(P3) here utilized to substitute (P2), is defined as follows:
min{𝐴 𝑚 , ⌈∑𝑖 𝑦̃𝑖𝑚 ⌉} − ∑𝑖 𝑦̂𝑖𝑚 to a set (i.e., 𝑄); to each element
(P3) : 𝑓2∗ = Min 𝑓2 , in 𝑄, supposing its index on TF is 𝑖, let 𝑦̂𝑖𝑚 = 𝑦̂𝑖𝑚 + 1. Finally,
let 𝑦𝑖𝑚 = 𝑦̂𝑖𝑚 .
(3) – (11) ,
{
{
{
{(15) – (20) , (27) Obviously, after executing this step, the constraint sets
s.t. : {𝑟𝜉 ≥ (1 − 𝛼) 𝑟̃𝜉 , ∀𝑗 ∈ 𝐽, 𝑘 ∈ 𝐾, 𝜉 ∈ Ξ, relating to 𝑦𝑖𝑚 can be satisfied. By now, as the solutions
{
{ 𝑗𝑘 𝑗𝑘 of all scenario-unrelated variables in (P1) have been iden-
{ 𝜉
{(21) \ {Δ𝑟𝑗𝑘𝑘󸀠 } , tified, (P1) can be decomposed into independent subprob-
lems by scenario. The next step is to obtain the solution
𝜉 𝜉 𝜉
where 𝑟̃𝑗𝑘 is the solution of 𝑟𝑗𝑘 after solving (P1) and 𝛼 is an of 𝑟𝑗𝑘 in (P1) with given values on all scenario-unrelated
adjustment coefficient to make a tradeoff between fill rate and variables.
Mathematical Problems in Engineering 7

Step 4. Add the constraint sets 𝑥𝑖 = max𝜉∈Ξ {𝑥̂𝑖𝜉 } (∀𝑖 ∈ 𝐼), bound (LB) on the optimal objective value of the original
𝜉 problem, (P3). A major advantage of this process is that it
𝑠𝑖𝑘 = max𝜉∈Ξ {̂𝑠𝑖𝑘 } (∀𝑖 ∈ 𝐼, 𝑘 ∈ 𝐾), and 𝑦𝑖𝑚 = 𝑦̂𝑖𝑚 (∀𝑖 ∈
𝐼, 𝑚 ∈ 𝑀) to (P1), then decompose (P1) into subproblems allows the problem to be split into separate subproblems for
by scenario, and solve each subproblem, denoting the value each scenario. The details of the subgradient algorithm are
𝜉 𝜉 given below, in which the settings of all parameters employed
of 𝑟𝑗𝑘 as 𝑟̃𝑗𝑘 .
were obtained using preliminary experimentation.

4.2. Solution to (P3) Step 0 (initialization). 𝑡 = 0 (iteration counter),

4.2.1. Model Transformation. To decompose (P3) effectively, a 𝜉,(𝑡) 𝜉,(𝑡)

new model named (P3󸀠 ) is defined by duplicating 𝑥𝑖𝜉 of 𝑥𝑖 , 𝑠𝑖𝑘𝜉 𝑢𝑖(𝑡) = 0, 𝑤𝑖𝑘 = 0, V𝑖𝑚 = 0. (33)
𝜉
of 𝑠𝑖𝑘 , and 𝑦𝑖𝑚 of 𝑦𝑖𝑚 for all scenarios 𝜉 ∈ Ξ on the model (P3);
then the following three constraint sets are added to (P3󸀠 ): Let LR(𝑡) = 𝑓2LR (u(𝑡) , w(𝑡) , k(𝑡) ) and UB(𝑡) = upper bound
value of the objective of (P3). Initialize LR(𝑡) = −∞ and
(1 − 𝑝1 ) 𝑥𝑖1 = ∑ 𝑝𝜉 𝑥𝑖𝜉 , ∀𝑖 ∈ 𝐼, (28) UB(𝑡) = +∞.
𝜉∈Ξ\{1}
Step 1 (solve each subproblem). With given u(𝑡) , w(𝑡) , and k(𝑡) ,
1 𝜉
𝑠𝑖𝑘 − 𝑠𝑖𝑘 = 0, ∀𝑖 ∈ 𝐼, 𝑘 ∈ 𝐾, 𝜉 ∈ Ξ \ {1} , (29) decompose 𝑓2LR (u, w, k) into subproblems by scenario, solve
1 𝜉
each subproblem using a linear programming solver, and add
𝑦𝑖𝑚 − 𝑦𝑖𝑚 = 0, ∀𝑖 ∈ 𝐼, 𝑚 ∈ 𝑀, 𝜉 ∈ Ξ \ {1} , (30)
its objective value to LR(𝑡+1) .
where the number 1 in (28)–(30) represents the first scenario. Keep 𝑥𝑖𝜉,(𝑡) , 𝑠𝑖𝑘
𝜉,(𝑡)
, and 𝑦𝑖𝑚 𝜉,(𝑡)
(the values of 𝑥𝑖𝜉 , 𝑠𝑖𝑘 𝜉
, and 𝑦𝑖𝑚 𝜉

Equation (28) ensures that the decisions on whether estab- gotten in iteration 𝑡), respectively.
lishing a facility at candidate location 𝑖 under all scenarios are ∀𝑖 ∈ 𝐼, if (𝑥𝑖1,(𝑡) = 𝑥𝑖2,(𝑡) = ⋅ ⋅ ⋅ = 𝑥𝑖|Ξ|,(𝑡) ), then add
identical, that is, 𝑥𝑖1 = 𝑥𝑖2 = ⋅ ⋅ ⋅ = 𝑥𝑖|Ξ| , and (29) and (30) state constraint set 𝑥𝑖𝜉 = 𝑥𝑖1,(𝑡) (∀𝜉 ∈ Ξ) to 𝑓2𝐿𝑅 (u, w, k).
𝜉 𝜉 1,(𝑡) 2,(𝑡) |Ξ|,(𝑡)
similar limitations on 𝑠𝑖𝑘 and 𝑦𝑖𝑚 , respectively. From (28)– ∀𝑖 ∈ 𝐼 and ∀𝑘 ∈ 𝐾, if (𝑠𝑖𝑘 = 𝑠𝑖𝑘 = ⋅ ⋅ ⋅ = 𝑠𝑖𝑘 ), then
(30), it can be easily concluded that the two problems, (P3) 𝜉 1,(𝑡)
add constraint set 𝑠𝑖𝑘 = 𝑠𝑖𝑘 (∀𝜉 ∈ Ξ) to 𝑓2LR (u, w, k).
and (P3󸀠 ), are equivalent. 1,(𝑡) 2,(𝑡) |Ξ|,(𝑡)
∀𝑖 ∈ 𝐼 and ∀𝑚 ∈ 𝑀, if (𝑦𝑖𝑚 = 𝑦𝑖𝑚 = ⋅ ⋅ ⋅ = 𝑦𝑖𝑚 ), then
(P3󸀠 ) can be decomposed into independent subproblems 𝜉 1,(𝑡)
by scenario if constraint sets (28)–(30) are relaxed. Therefore, add constraint set 𝑦𝑖𝑚 = 𝑦𝑖𝑚 (∀𝜉 ∈ Ξ) to 𝑓2 (u, w, k). LR

the Lagrangian relaxation of (P3󸀠 ) with respect to the con-

straint sets (28)–(30) can be described as follows: Step 2 (compute subgradients).

𝑓2LR (u, w, k) subgrad 𝑢𝑖(𝑡+1) = 𝑥𝑖1,(𝑡) − ∑ 𝑝𝜉 𝑥𝑖𝜉,(𝑡) ,

𝜉∈Ξ
{
= min {𝑓2 + ∑ (𝑝𝜉 ∑ 𝑢𝑖 (𝑥𝑖1 − 𝑥𝑖𝜉 ) 𝜉,(𝑡+1) 1,(𝑡) 𝜉,(𝑡) (34)
𝜉∈Ξ\{1} 𝑖∈𝐼
subgrad 𝑤𝑖𝑘 = 𝑠𝑖𝑘 − 𝑠𝑖𝑘 ,
{
𝜉 1 𝜉 𝜉,(𝑡+1) 1,(𝑡) 𝜉,(𝑡)
+ ∑ ∑ 𝑤𝑖𝑘 (𝑠𝑖𝑘 − 𝑠𝑖𝑘 ) subgrad V𝑖𝑚 = 𝑦𝑖𝑚 − 𝑦𝑖𝑚 .
𝑖∈𝐼 𝑘∈𝐾

Step 3 (update upper bound).

𝜉 1 𝜉 }
+ ∑ ∑ V𝑖𝑚 (𝑦𝑖𝑚 − 𝑦𝑖𝑚 )) }
𝑖∈𝐼 𝑚∈𝑀
} 𝑥𝑖(𝑡) = max {𝑥𝑖𝜉,(𝑡) } ,
𝜉∈Ξ
s.t.: all the constraints of (P3󸀠 ) \ constraints (28) – (30) , (𝑡) 𝜉,(𝑡)
(35)
(31) 𝑠𝑖𝑘 = max {𝑠𝑖𝑘 }.
𝜉∈Ξ

where 𝑓2 is the objective function of (P3), and the vectors u,

(𝑡)
w, and k in (31) are the Lagrangian multipliers associated with Compute the values of 𝑦𝑖𝑚 (∀𝑖 ∈ 𝐼, 𝑚 ∈ 𝑀) following
relaxed constraint sets (28), (29), and (30), respectively. The the adjustment process of Step 3 in Section 4.1 based on the
𝜉,(𝑡)
Lagrangian dual problem can be stated as follows: values of 𝑦𝑖𝑚 (∀𝜉 ∈ Ξ), and then solve (P3) where the
scenario-unrelated variables, 𝑥𝑖 , 𝑠𝑖𝑘 , and 𝑦𝑖𝑚 , are assigned to
𝑧LD = max 𝑓LR (u, w, k) .
u,w,k 2
(32) the values of 𝑥𝑖(𝑡) , 𝑠𝑖𝑘
(𝑡) (𝑡)
, and 𝑦𝑖𝑚 , respectively. Supposing the
objective of (P3) is UB, let UB(𝑡+1) = min{UB, UB(𝑡) }.
4.2.2. Subgradient Algorithm. The lagrangian dual (32) is Record the values of 𝑥𝑖 , 𝑠𝑖𝑘 , and 𝑦𝑖𝑚 that have been utilized
a convex nonsmooth program that can be solved by sub- to generate the updated upper bound as 𝑥̂𝑖 , 𝑠̂𝑖𝑘 , and 𝑦̂𝑖𝑚 ,
gradient algorithm and the result can be taken as a lower respectively.
8 Mathematical Problems in Engineering

Step 4 (identify new step direction). Step 6 (stopping criteria).

󵄩󵄩 󵄩󵄩
󵄩 UB(𝑡+1) − LR(𝑡+1) 󵄩󵄩
stepdirc 𝑢𝑖(𝑡+1) = 𝛽1 (subgrad 𝑢𝑖(𝑡+1) ) if ((sum of squares = 0) 󵄩󵄩󵄩󵄩( ≤ 𝜀) 󵄩󵄩
󵄩󵄩
󵄩󵄩 UB(𝑡+1) 󵄩
(39)
+ 𝛽2 (subgrad 𝑢𝑖(𝑡) )
× (𝑡 ≥ 𝑁2 ))
+ 𝛽3 (subgrad 𝑢𝑖(𝑡−1) ) ,
𝜉,(𝑡+1) 𝜉,(𝑡+1) stop;
stepdirc 𝑤𝑖𝑘 = 𝛽1 (subgrad 𝑤𝑖𝑘 ) else
𝜉,(𝑡) 𝑡 = 𝑡 + 1, goto Step 1
+ 𝛽2 (subgrad 𝑤𝑖𝑘 ) where the value of 𝑁2 is set to 200 and the value of 𝜀 is set to
(36) 0.001.
𝜉,(𝑡−1)
+ 𝛽3 (subgrad 𝑤𝑖𝑘 ), In Step 0, the Lagrangian multipliers are initialized with
𝜉,(𝑡+1) 𝜉,(𝑡+1)
0, and the lower and upper bounds of (P3) are initialized
stepdirc V𝑖𝑚 = 𝛽1 (subgrad V𝑖𝑚 ) with −∞ and +∞, respectively. With the values of the
𝜉,(𝑡)
Lagrangian multipliers assigned in Step 0, each subproblem
+ 𝛽2 (subgrad V𝑖𝑚 ) is solved optimally using a linear programming solver, and
the values of 𝑥𝑖 , 𝑠𝑖𝑘 , and 𝑦𝑖𝑚 are duly fixed in Step 1 if their
𝜉,(𝑡−1)
+ 𝛽3 (subgrad V𝑖𝑚 ), solutions under all scenarios are identical at this intermediate
iteration, and in this way the computation complexity can
where 𝛽1 = 0.6, 𝛽2 = 0.3, 𝛽3 = 0.1. be effectively decreased in latter iterations. Step 2 computes
new subgradient values and Step 3 updates the upper bound
Step 5 (update multipliers). 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,
sum of squares and by updating UB the convergence speed can be accelerated
2 [24]. After updating the upper bound, a new step direction is
= ∑ (subgrad 𝑢𝑖(𝑡+1) ) found in Step 4 by taking the convex combinations of the last
𝑖∈𝐼 three subgradient values, and the typical oscillating behavior
𝜉,(𝑡+1) 2 (37) of the subgradient algorithm can be reduced adopting this
+ ∑ ∑ ∑ (subgrad 𝑤𝑖𝑘 ) combination operation [25, 26].
𝑖∈𝐼 𝑘∈𝐾 𝜉∈Ξ The subgradient algorithm does not always generate
𝜉,(𝑡+1) 2 primal feasible solutions to (P3󸀠 ) (which is equivalent to
+ ∑ ∑ ∑ (subgrad V𝑖𝑚 ), (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
if (sum of squares = 0) is needed to restore feasibility starting from the dual feasible
goto Step 6 solution provided by the subgradient algorithm.
else
4.2.3. Feasibility Algorithm. A primal feasible solution to (P3)
2.0 × 𝜂𝑡 (UB(𝑡+1) − LR(𝑡+1) ) can be obtained in a greedy manner starting from the values
𝑢𝑖(𝑡+1) = 𝑢𝑖(𝑡) + stepdirc 𝑢𝑖 (𝑡+1)
, that generate the UB by the subgradient algorithm. The details
sum of squares
of the feasibility heuristic are as follows.
𝜉,(𝑡+1) 𝜉,(𝑡)
𝑤𝑖𝑘 = 𝑤𝑖𝑘 (i) Add the constraint sets 𝑥𝑖 = 𝑥̂𝑖 (∀𝑖 ∈ 𝐼), 𝑠𝑖𝑘 = 𝑠̂𝑖𝑘 (∀𝑖 ∈
𝑡 (𝑡+1) (𝑡+1) 𝐼, 𝑘 ∈ 𝐾), and 𝑦𝑖𝑚 = 𝑦̂𝑖𝑚 (∀𝑖 ∈ 𝐼, 𝑚 ∈ 𝑀) to
𝜉,(𝑡+1)
2.0 × 𝜂 (UB − LR ) (P3), in which 𝑥̂𝑖 , 𝑠̂𝑖𝑘 , and 𝑦̂𝑖𝑚 are the recorded values
+ stepdirc 𝑤𝑖𝑘 ,
sum of squares of 𝑥𝑖 , 𝑠𝑖𝑘 , and 𝑦𝑖𝑚 in generating the UB in Step 3 of
the subgradient algorithm. After adding these three
𝜉,(𝑡+1) 𝜉,(𝑡)
V𝑖𝑚 = V𝑖𝑚 constraint sets, (P3) merely contains scenario-specific
decision variables. In this way, (P3) can further
𝜉,(𝑡+1)
2.0 × 𝜂𝑡 (UB(𝑡+1) − LR(𝑡+1) ) be decomposed into independent subproblems by
+ stepdirc V𝑖𝑚 ,
sum of squares scenario;
(38) (ii) Decompose (P3) into |Ξ| subproblems by
scenario, and then solve each subproblem using
where 𝜂 is a user-defined parameter and the value 𝜂 = 0.99 is a linear programming solver. Supposing the
used in the experiments.
𝜉
values of 𝑛𝑓𝑖𝑗𝑚 𝜉
and 𝑛𝑏𝑗𝑖𝑚 are 𝑛𝑓̂𝑖𝑗𝑚
𝜉
and 𝑛̂𝑏𝑗𝑖𝑚
𝜉
,
Mathematical Problems in Engineering 9

J
6
H
I
K F
7 5
G
L O 4
E
M
N D
C 3
8 A
B
2

Figure 1: Map of candidate TFs and AAs.

respectively, the objective value of (P3) equals 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
5. Case Study the data set of scenario 1, respectively. Similarly, another 4
5.1. Test Case. A case study focused on the disaster relief scenarios which represented lighter damage were produced
logistics of the Great Wenchuan Earthquake is here used using random coefficients between 0.5 and 1.0. The occur-
to illustrate the utility of the temporary facility location- rence probabilities of these 9 scenarios were hypothetically
allocation model and verify the proposed heuristic decompo- set to 0.20, 0.15, 0.15, 0.15, 0.15, 0.05, 0.05, 0.05, and 0.05,
sition algorithms. This great earthquake disaster occurred at respectively.
the Long Menshan fault, a thrust structure located along the The parameters of candidate TFs are given in Tables
border of the Eurasian Plate and Indo-Australian Plate, and 2 and 3. The maximum capacity of each relief supply at
this fault consists of three parallel subfaults, that is, the front each candidate TF mainly depends on its original inventory
fault (Guanxian-Jiangyou fault), the central fault (Yingxiu- level and its inbound transportation capacity. For example,
Beichuan fault), and the back fault (Wenxian-Maoxian fault), candidate TF1 is a large urban city that not only has larger
along the north-east direction. The epicenter was located initial and potential capacities with respect to warehouse and
near Wenchuan County at a depth of about 14 km. According health care facilities but also has more accessible inbound
to the above information and the discussion results with transportation capabilities, and thus its maximum capacities
subject matter experts of Earthquake Agency of P. R. China, on all categories are larger than TF2 to TF5. In contrast,
15 disaster areas in the disaster region were chosen as AAs because the urban areas of TF6 to TF8 are relatively small
and 8 cities or towns, in which the inbound transportation and inbound capacities of these TFs are also relatively low,
condition by land remained relatively good was selected as the maximum capacities of these TFs are set lower than
candidate TFs around this region. The map of AAs and TFs the capacities of other TFs. Similarly, the differences on
is illustrated in Figure 1 where each AA is represented by a mobilization cost are also set among the 8 candidate TFs,
unique alphabet from “A” to “O” and each candidate TF is and these differences mainly reflect the variations in inbound
here marked by a unique number from “1” to “8,” and all the transportation costs. The fixed establishment costs of TF1–
identified AAs were found to be aligned nearly in three rows TF8 were set as 25, 30, 30, 30, 45, 45, 45, and 45 thousand US
corresponding to the three subfaults. dollars, respectively, indicating that the cost of establishing
8 types of demand, including water, food, medicine, TFs with good inbound transportation conditions (i.e., TF1,
doctors, nurses, rescuers, tents, and critical population, were TF2, and TF3) was lower than these with limited inbound
considered in this case study. Table 1 gives the demand levels transportation conditions (i.e., TF6, TF7, and TF8). The
of AAs constructed from the loss assessment results of NDRC setting of fixed cost of establishing a TF was not associated
of P. R. China [27] and Yuan [28], in which the number of with its total mobilization levels of relief supplies as indicated
requisitions in depots H, K, L, and M is relatively higher in previous works (i.e., Özdamar et al. [29], Rawls and
than those of other AAs because these AAs are affected more Turnquist [6], Salmerón and Apte [9], and Bozorgi-Amiri et
10 Mathematical Problems in Engineering

Table 1: Demand levels of AAs within scenario 1 (units).

AAs Water Food Medicine Doctors Nurses Rescuers Tents Critical population
A 1660 700 360 68 98 106 488 240
B 1826 770 396 75 108 117 537 264
C 2490 1050 540 102 147 159 732 360
D 2241 945 486 92 132 143 659 324
E 1378 581 299 56 81 88 405 199
F 2216 935 481 91 131 142 651 320
G 2407 1015 522 99 142 154 708 348
H 9960 4200 2160 408 588 636 2928 1440
I 2822 1190 612 116 167 283 1303 641
J 4980 2100 1080 204 294 318 1464 720
K 10375 4375 2250 425 613 663 3050 1500
L 12948 5460 2808 530 764 827 3806 1872
M 6225 2625 1350 255 368 398 1830 900
N 2075 875 450 85 123 133 610 300
O 2200 928 477 90 130 140 647 318

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

TFs Water Food Medicine Doctors Nurses Rescuers Tents Hospital beds
1 39500 19700 8400 1950 3260 2750 9750 6200
2 25675 12805 5460 1268 2119 1788 6338 4030
3 33575 16745 7140 1658 2771 2338 8288 5270
4 32390 16154 6888 1599 2673 2255 7995 5084
5 27650 13790 5880 1365 2282 1925 6825 4340
6 9875 4925 2100 488 815 688 2438 1550
7 8295 4137 1764 410 685 578 2048 1302
8 6715 3349 1428 332 554 468 1658 1054

al. [10]). This is because these open TFs in our model were and for medicine, doctors, and nurses was initially set as
used solely as temporary distribution depots for postdisaster 0.001, denoting that there exists synergic requirement on
relief and not as supply facilities for natural disaster asset distribution of food and water as well as on distribution of
prepositioning as in previous works. medicine, doctors, and nurses. The penalty of discrepancy
Table 4 shows a summary of the data used for helicopters. among other relief supplies was set as zero. These values
To the incorporated 5 modes of helicopters, Z-8S and M- are useful for illustration purposes, but they may not match
170S were classified as special mode of helicopters and the values that might be estimated by ultimate users of the model.
remaining 3 were classified as general helicopters. The num-
ber of hours that each helicopter was available here included 5.2. Computational Results. The proposed decomposition
available flight time, not time required for maintenance, algorithm was implemented using CPLEX 9.0 Callable
refueling, shift changes, or land operations. The average flight Library embedded in C++ on a 2.4 GHz desktop computer
time between each TF and each AA was calculated based with 4 GB of RAM. The parameter of EPGAP in CPLEX,
on the flight speed of each helicopter mode and the distance which determines the relative tolerance on the MIP gap
between each TF and each AA (i.e., Figure 1). between the best integer objective and the objective of the
The adjustment coefficient on tradeoff, 𝛼, was initial- best node remaining, was revised from default setting, 1.0𝐸-4,
ized as 0.01, and the weighting parameters of all relief to 1.0𝐸-3 in the present case study in order to increase the
supplies were set as follows: the total criticality weight of convergence speed of computation.
the three subsets of relief supplies (i.e., relief commodi-
ties, relief personnel, and hospital beds) was first set to
5.2.1. Explanation and Analysis of the Solution. Tables 5
0.65, 0.25, and 0.1, respectively; then, to each type of
and 6 summarize the solutions of the scenario-unrelated
relief supplies within the relief commodities subset and
decision variables. Table 5 shows that 6 TFs had opened and
the relief personnel subset, its weighting parameter was
distributed widely around the disaster region. Among the
set to 0.65 ∑𝑗∈𝐽 ℎ𝑘 𝑑𝑗𝑘 / ∑𝑘∈𝐾𝑅 ∑𝑗∈𝐽 ℎ𝑘 𝑑𝑗𝑘 (∀𝑘 ∈ 𝐾𝑅 ) and 6 open TFs, only TF 8 was not arranged to receive and
0.25 ∑𝑗∈𝐽 𝑑𝑗𝑘 / ∑𝑘∈𝐾𝑅 ∑𝑗∈𝐽 𝑑𝑗𝑘 (∀𝑘 ∈ 𝐾𝑆 ), respectively. The transship all types of relief supplies. Comparing with the
penalty of discrepancy in fill rates for water and food maximum capacities for relief supplies depicted in Table 1, the
Mathematical Problems in Engineering 11

Table 3: Mobilization cost for relief supplies at the TFs (103 $/unit).

TFs Water Food Medicine Doctors Nurses Rescuers Tents Hospital beds
1 0.006 0.22 0.80 1.70 0.96 1.12 0.20 0.50
2 0.007 0.24 0.88 1.87 1.06 1.23 0.22 0.55
3 0.007 0.24 0.88 1.87 1.06 1.23 0.22 0.55
4 0.007 0.24 0.88 1.87 1.06 1.23 0.22 0.55
5 0.007 0.24 0.88 1.87 1.06 1.23 0.22 0.55
6 0.008 0.26 0.96 2.04 1.15 1.34 0.24 0.60
7 0.008 0.26 0.96 2.04 1.15 1.34 0.24 0.60
8 0.008 0.26 0.96 2.04 1.15 1.34 0.24 0.60

Table 4: Summary data of helicopter modes. were not associated with any penalty of discrepancy, such as
rescuers and tents, were different. This means the synergic
Parameter Z-8G Z-8S S-70G M-170G M-170S requirement with respect to fill rate among correlated relief
Number of helicopters demands can be completely satisfied at a lower setting on the
20 15 16 15 12
available 𝐶
penalty coefficient of 𝜔𝑘𝑘 󸀠 , that is, 0.001.
Available hours 18 18 21 18 18 Table 8 depicts the fill rates on the 8 categories of demands
Capacity for commodities and the integrated fill rate (which incorporates the criticality
4.0 0 3.5 4.5 0
(tons) weight of each category of demands) under 3 representative
Capacity for personnel scenarios: the levels of demand of scenario 1 are summarized
18 20 16 28 32
(persons) in Table 1; and those of scenario 2 and scenario 6 individually
Capacity for critical represent a scenario with more severe damage than scenario
8 15 6 9 16
population (persons) 1 and a scenario with lighter damage, respectively. Regarding
3
Hiring cost (10 $/unit) 12 12 10 15 15 the 3 representative scenarios, the individual fill rate on any
3
Operating cost (10 $/unit) 0.5 0.5 0.4 0.6 0.6 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
mobilization levels of most relief supplies did not reach their level of demand was highest under scenario 2 and lowest
maximum capacities at any open TF, indicating that the total under scenario 6. As the integrated fill rates under the 3
capacities of the 6 open TFs had not been completely utilized. representative scenarios were not high (i.e., 56.8%, 44.6%,
As shown in Table 6, all modes of helicopters were allocated and 77.0% under scenario 1, 2, and 6, resp.) but the mobi-
to TFs 1, 6, and 7. However, no helicopter was allocated to TFs lization levels of most relief supplies, at the same time, had
3, 4, and 8, and the total number of helicopters allocated to not reach their maximum capacities at the open TFs (Table 2),
open TFs, such as zero helicopter allocated to TFs 3, 4, and this inconsistency was probably caused by the insufficiency of
8, did not match the overall level of mobilized relief supplies available helicopters, which has been analyzed before. Indeed,
of these TFs. This is because Table 6 only represents the after adding 5 helicopters on each mode, the integrated fill
initial deployment results of helicopters, and in the succedent rates of the 3 representative scenarios increased from 56.8%,
distribution process, each helicopter was permitted to return 44.6%, and 77.0% to 87.1%, 72.0%, and 99%, respectively. This
to a depot other than the one from which it had departed. indicates that the proposed model can help disaster relief
According to Tables 4 and 6, all available helicopters had been decision makers identify potential resource bottlenecks, such
initially deployed to open TFs, indicating that helicopters as the maximum capacities of relief supplies at open TFs, the
may be a type of deficient resource to further improve the overall distribution capacities of helicopters. In summary, to
expected fill rate of AAs. any type of demands whose fill rate must be improved, if
The fill rates of all demand types at AAs under scenario the corresponding mobilization levels of all open TFs reach
1 are shown in Table 7. As shown, to any type of demand, their maximum capacities, then the total capacities of open
the fill rates from all the 15 AAs are equivalent. That is, the TFs should first be considered to be expanded, or more TFs
distribution of relief among AAs was fair in fill rate, indicating should be established. Otherwise, more helicopters should be
that the max-min method employed in (1) reaches its original recruited.
goal. The fill rates among relief demands with any nonzero
𝐶 𝐶
setting (𝜔𝑘𝑘 󸀠 = 0.001) on the penalty of discrepancy (such 5.2.2. Effects of the Penalty Coefficient 𝜔𝑘𝑘󸀠 and the Adjustment

as, the choosing relief supplies between food and water or Coefficient 𝛼. As the synergic requirements among correlated
among medicine, doctors, and nurses) are also equal. For relief supplies (i.e., medicine, doctors, and nurses) could be
𝐶
example, the fill rates of water and food at any AA were equal achieved after setting the penalty coefficient 𝜔𝑘𝑘 󸀠 = 0.001,

to 51.6%, and the demands for medicine, doctors, and nurses to further investigate the influence of this penalty coefficient,
from all the AAs also had the same fill rate value, 48.3%. we modified this coefficient’s value while keeping the settings
In contrast, the fill rates among those relief demands that of other parameters unchanged. Figure 2 shows the results
12 Mathematical Problems in Engineering

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

TFs Selected Water Food Medicine Doctors Nurses Rescuers Tents Hospital beds
1 True 2223 1947 1789 293 214 344 821 648
2 False — — — — — — — —
3 True 4331 1677 915 155 278 367 1107 1208
4 True 8609 3918 2754 617 736 1370 3379 2548
5 False — — — — — — — —
6 True 9875 2944 1822 279 534 688 2393 1550
7 True 8295 3863 1728 396 544 548 2048 1302
8 True 2018 — — — 88 255 1009 730

2 100

80

Expected fill rate on 3 scenarios (%)

1.5
Total discrepancy in fill rate

60

40

0.5
20

0 0
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
𝛽 𝛽
(a) (b)

𝐶
Figure 2: Total discrepancy and expected fill rate with varied 𝑤𝑘𝑘󸀠.

𝐶 −𝛽
Table 6: Initial deployment results of helicopter modes at TFs. respectively. 𝜔𝑘𝑘 󸀠 was modified by multiplying with 10 ,
and the range of 𝛽 is depicted by the 𝑥-axis of Figures 2(a)
TFs Z-8G Z-8S S-70G M-170G M-170S and 2(b). According to Figure 2(a), there is no discrepancy
1 11 — 10 9 — 𝐶
in fill rate among correlated relief supplies when 𝜔𝑘𝑘 󸀠 is
2 — — — — — set as 0.001. As the value of penalty coefficient decreases
3 — — — — — further, the discrepancy among correlated relief supplies
𝐶 −5
4 — — — — — increases. When 𝜔𝑘𝑘 󸀠 falls below 1×10 , the total discrepancy
5 — — — — — ceases to increase. However, as depicted in Figure 2(b), the
𝐶
6 3 — — 1 11 system’s expected fill rate does not change markedly as 𝜔𝑘𝑘 󸀠

7 6 15 6 5 1 varies. These phenomena can provide valuable information

8 — — — — — on setting the value of this penalty coefficient of discrepancy.
Total 20 15 16 15 12
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
with the modification of 𝜔𝑘𝑘 󸀠 , in which the correlation

assumption among relief supplies was kept unchanged as problem was solved several times while varying the adjust-
aforementioned. The indices of the 𝑦-axis in Figures 2(a) ment coefficient, 𝛼, and the results are illustrated in Figure 3.
and 2(b), total discrepancy in fill rate and the expected fill According to Figure 3(a), decreasing the value of 𝛼 causes the
rate on 3 representative scenarios, were calculated through values of the objective (1) and the objective (2) to increase
𝜉 𝜉
∑𝜉∈Ξ ∑𝑗∈𝐽 ∑𝑘∈𝐾 ∑𝑘󸀠 ∈𝐾 𝑝𝜉 Δ𝑟𝑗𝑘𝑘 𝜉
󸀠 /|𝐽| and ∑𝜉∈Ξ ∑𝑘∈𝐾 𝑝 𝑤𝑘 𝑅𝑘 ,
synchronously. In other words, as the overall fill rate of
Mathematical Problems in Engineering 13

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

AAs Water Food Medicine Doctors Nurses Rescuers Tents Critical population
A 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
B 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
C 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
D 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
E 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
F 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
G 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
H 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
I 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
J 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
K 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
L 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
M 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
N 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9
O 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9

×104
3.25 9

3.2 𝛼=0 8

3.15 7
𝛼 = 0.01
The value of objective (2)

Number of opened TFs

3.1 6
𝛼 = 0.02
3.05 5
𝛼 = 0.03
3 4
𝛼 = 0.04
2.95 3
𝛼 = 0.05
2.9 2
𝛼 = 0.06
2.85 1
𝛼 = 0.07
2.8 0
0.54 0.55 0.56 0.57 0.58 0.59 0 0.02 0.04 0.06 0.08
The value of objective (1) 𝛼
(a) (b)

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

AAs increases, the procuring cost of relief supplies and the decision variables, and Table 11 depicts the fill rates under 3
distribution cost will increase, and more TFs are needed representative scenarios. As shown in Table 9, the location
to be established (the number of open TFs is depicted in results of open TFs are equivalent by utilizing the heuristic
Figure 3(b)). 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
5.2.3. Solution Quality and Evaluation of Computational obvious differences with respect to the initial deployment
Speed. To evaluate the performance of the proposed result of helicopters at open TFs according to Tables 6
decomposition-based heuristic algorithm, the problem and 10, but, as shown in Tables 8 and 11, the expected fill
described in Section 5.1 was solved using an ILOG CPLEX rates under the 3 representative scenarios achieved by the
engine with the same settings (i.e., EPGAP). Tables 9 heuristic algorithm are very similar to the corresponding fill
and 10 summarize the solutions of scenario-unrelated rates computed through CPLEX. For example, the difference
14 Mathematical Problems in Engineering

Expected value of total cost (million dollars) 50 12000

10000
40

8000

Solve time (s)

30
6000

20
4000

10 2000

0 0
0 1 2 3 4 5 0 1 2 3 4 5
Number of helicopters added of each mode Number of helicopters added of each mode

CPLEX CPLEX
Heuristic algorithm Heuristic algorithm
(a) (b)

Figure 4: Performance comparing results between the heuristic algorithm and the CPLEX solver.

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

Scenario Water Food Medicine Doctors Nurses Rescuers Tents Critical population Integrated
1 51.6 51.6 48.3 48.3 48.3 78.5 53.4 81.9 56.8
2 29.2 29.2 50.9 50.9 50.9 68.4 42.2 61.4 44.6
6 73.1 73.1 77.8 77.8 77.8 99.0 53.0 98.8 77.0

is approximately 0.1% under each of the 3 scenarios. This and the total cost computed by the heuristic algorithm was
is because the individual fill rates had been obtained by slightly higher (from 1.28% to 2.59%) than the cost produced
solving (P1) using the method described in Section 4.1, by the CPLEX solver. However, the computation time of the
and in the procedure of that method, only the scenario- heuristic algorithm was much shorter than the time required
unrelated variable 𝑦𝑖𝑚 (∀𝑖 ∈ 𝐼, 𝑚 ∈ 𝑀) was assigned with by the CPLEX solver (the ratios between the two ranged
an approximate value. However, 𝑦𝑖𝑚 denotes only the initial from 16.2% to 25.7%), denoting that our heuristic algorithm
deployment results of all modes of helicopters at open TFs, exhibited a marked advantage in computation time than a
and during the detailed distribution process helicopters commercial MIP solver.
were permitted to return to any open TF after finishing its
current forward delivery mission. As a result, our proposed 6. Conclusions and Future Work
heuristic algorithm was able to produce approximate fill
rates compared with these using CPLEX. The expected value For the difficulty in predicting earthquakes, the mobilization
of total cost computed using our proposed algorithm is of emergency resources is essential after disasters, espe-
1.80% higher than the value gotten by CPLEX (31.302 million cially when an earthquake disaster, for example, the Great
dollars versus 30.749 million dollars), but the computation Wenchuan Earthquake that occurred on May 12, 2008, in
time is 15.8% of that required by CPLEX (1696 seconds versus China, was devastating. We developed a disaster relief opti-
10600 seconds) on the same computer. mization model whose solution provides a strategic plan for
To further investigate the performance of the proposed locations of TFs, mobilization levels of relief supplies, initial
algorithm on the expected value of total cost and on compu- deployment of helicopters at open TFs, and the distribution
tation time, the aforementioned problem was solved several plans under each revealed scenario where the demands of
times by varying the number of helicopters of each mode with AAs are presented by probabilistic scenarios. The model is
the proposed algorithm and the CPLEX solver, respectively. a stochastic mixed-integer programming problem featuring
The results are depicted in Figure 4, in which the numerical two objectives: maximizing the expected minimal fill rate
value of the 𝑥-axis denotes the number of added helicopters of affected areas where the mismatching distribution among
on each mode. As shown in Figure 4(a), the expected value of correlated relief demands is penalized and minimizing the
total cost increased as the number of helicopters increased, expected value of total cost. A lexicographic multiobjective
Mathematical Problems in Engineering 15

Table 9: Mobilization levels of relief supplies at TFs using CPLEX (units).

TFs Selected Water Food Medicine Doctors Nurses Rescuers Tents Hospital beds
1 True 2235 1331 1022 159 141 248 535 525
2 False — — — — — — — —
3 True 4676 1689 1114 156 324 432 1403 1083
4 True 9905 3942 2510 600 770 1320 3397 2628
5 False — — — — — — — —
6 True 9875 3362 2100 324 537 688 2438 1550
7 True 8295 4104 1764 410 635 578 2048 1302
8 True 561 — 547 101 — 327 994 943

Table 10: Initial deployment results of helicopter modes at TFs using CPLEX.

TFs Z-8G Z-8S S-70G M-170G M-170S

1 2 — 6 1 —
2 — — — — —
3 — — — — —
4 — — — — —
5 — — — — —
6 8 2 8 1 12
7 10 13 2 13 —
8 — — — — —
Total 20 15 16 15 12

Table 11: Fill rates of all demands under the 3 representative scenarios using CPLEX (percent).

Scenario Water Food Medicine Doctor Nurse Rescuer Tent Critical population Integrated
1 52.8 52.8 45.7 45.7 45.7 82.8 53.4 81.8 56.9
2 27.6 27.6 51.5 51.5 51.5 72.3 42.2 60.8 44.6
6 73.5 73.5 77.8 77.8 77.8 99.9 52.7 98.7 77.1

optimization technique was used to transform the proposed relief demands almost disappears. In addition, the adjustment
model into two sequential single objective stochastic models. coefficient, 𝛼, can have good effect on making tradeoffs
For relatively small data instances, the transformed models between the expected fill rate and the expected value of
can be solved using a commercial mixed-integer solver total cost. Across a test problem, the proposed algorithm
such as CPLEX, but this approach does not scale well to consistently produces solutions whose expected overall fill
large problems. As a result, scenario-decomposition-based rates are almost the same as the optimal values and whose
heuristic algorithms had also been proposed to solve the total costs remained within 2.59% higher than the costs of the
transformed models, respectively, and this approach scales optimal, with computation times that are between 16.2% and
well, to allow large problems to be solved. 25.7% of the time required using commercial MIP software.
A numerical case, which is based on scenarios of the These experiments indicate that the proposed algorithm can
Great Wenchuan Earthquake of China, illustrates both the be utilized as a large-scale disaster relief logistics planning
formulation and the character of the problem solutions in tool.
a practical context. The analysis of the results from the The character of the discrete distribution in probability
numerical case shows, for example, that the proposed model for uncertain demand also implies a desirable direction for
can help government officials to identify the bottlenecks further enhancing the applicability of the model. That is,
hindering the improvement on the expected overall fill rate, as the forecasting results on uncertain demands describing
such as the mismatching between the transportation capacity as scenarios may be inaccurate early after a disaster, the
and the qualities of relief supplies to be distributed. The postdisaster conditions could become more and more clear
𝐶
varied settings of the penalty coefficient, 𝜔𝑘𝑘 󸀠 , have obvious over time. Meanwhile, considering that the mobilization
𝐶
impact on the fill rate of correlated relief demands; when 𝜔𝑘𝑘 󸀠 campaign itself requires a period of time to be completely
reaches 0.001, the discrepancy in fill rate among correlated realized, we have opportunities to modify the original
16 Mathematical Problems in Engineering

mobilization plan according to updated information and the [8] H. O. Mete and Z. B. Zabinsky, “Stochastic optimization of med-
actual progress of the relief logistics. This offers an important ical supply location and distribution in disaster management,”
avenue for enhancement of the current work by establishing International Journal of Production Economics, vol. 126, no. 1, pp.
dynamic facility location and allocation models. 76–84, 2010.
Another desirable direction for future work is to further [9] J. Salmerón and A. Apte, “Stochastic optimization for natural
enhance the efficiency of the proposed heuristic algorithm. disaster asset prepositioning,” Production and Operations Man-
For example, after decomposing the original problem into agement, vol. 19, no. 5, pp. 561–574, 2010.
subproblems by scenario, it may be worthy to specially [10] A. Bozorgi-Amiri, M. S. Jabalameli, and S. M. J. Mirzapour Al-
design heuristic algorithms to solve each subproblem, such e-Hashem, “A multi-objective robust stochastic programming
as separating the forward distribution problem as commodity model for disaster relief logistics under uncertainty,” OR Spec-
flow subproblem and vehicle flow subproblem as in Haghani trum, 2011.
and Oh [30] with given values on all scenario-unrelated [11] A. Döyen, N. Aras, and G. Barbarosoğlu, “A two-echelon
decision variables or directly proposing heuristic algorithms stochastic facility location model for humanitarian relief logis-
to solve each facility location and allocation subproblem tics,” Optimization Letters, vol. 6, no. 6, pp. 1123–1145, 2012.
within each scenario, in order to further reduce the time [12] J. R. Birge and F. Louveaux, Introduction to Stochastic Program-
required for computation. ming, Springer, New York, NY, USA, 1997.
[13] C. C. Carøe and R. Schultz, “Dual decomposition in stochastic
Conflict of Interests integer programming,” Operations Research Letters, vol. 24, no.
1-2, pp. 37–45, 1999.
The authors declare that there is no conflict of interests [14] G. Laporte and F. V. Louveaux, “The integer 𝐿-shaped method
regarding the publication of this paper. for stochastic integer programs with complete recourse,” Oper-
ations Research Letters, vol. 13, no. 3, pp. 133–142, 1993.
Acknowledgments [15] S. Ahmed, M. Tawarmalani, and N. V. Sahinidis, “A finite
branch-and-bound algorithm for two-stage stochastic integer
This work was partly supported by the National Natural programs,” Mathematical Programming A, vol. 100, no. 2, pp.
Science Foundation of China under Grants 71371181 and 355–377, 2004.
91024006 and by the China Postdoctoral Science Foundation [16] H. D. Sherali and B. M. P. Fraticelli, “A modification of
under Grant 2012M521918. The authors thank Professor Chi- Benders’ decomposition algorithm for discrete subproblems:
Guhn Lee at the University of Toronto and Dr. Jian-Mai an approach for stochastic programs with integer recourse,”
Shi at the National University of Defense Technology of Journal of Global Optimization, vol. 22, no. 1–4, pp. 319–342,
2002.
China for their assistance in this research. They are also
grateful to the Guest Editor, Professor Lu Zhen, and to the [17] L. Ntaimo and S. Sen, “The million-variable “march” for
anonymous referees whose insightful comments led to a stochastic combinatorial optimization,” Journal of Global Opti-
significant improvement in this paper. mization, vol. 32, no. 3, pp. 385–400, 2005.
[18] S. Sen and J. L. Higle, “The C3 theorem and a D2 algorithm for
large scale stochastic mixed-integer programming: set convexi-
References fication,” Mathematical Programming A, vol. 104, no. 1, pp. 1–20,
2005.
[1] J.-B. Sheu, “Challenges of emergency logistics management,”
Transportation Research E, vol. 43, no. 6, pp. 655–659, 2007. [19] H. D. Sherali and X. Zhu, “On solving discrete two-stage
stochastic programs having mixed-integer first- and second-
[2] S. Qing, “Ten major sciences and technologies for disaster stage variables,” Mathematical Programming B, vol. 108, no. 2-3,
rescue work in earthquake,” Science & Technology Review, vol. pp. 597–616, 2006.
26, pp. 19–24, 2008.
[20] S. Sen and H. D. Sherali, “Decomposition with branch-and-cut
[3] B. Balcik and B. M. Beamon, “Facility location in humanitarian approaches for two-stage stochastic mixed-integer program-
relief,” International Journal of Logistics Research and Applica- ming,” Mathematical Programming A, vol. 106, no. 2, pp. 203–
tions, vol. 11, no. 2, pp. 101–121, 2008. 223, 2006.
[4] M.-S. Chang, Y.-L. Tseng, and J.-W. Chen, “A scenario planning [21] L. Li, M. Jin, and L. Zhang, “Sheltering network planning and
approach for the flood emergency logistics preparation problem management with a case in the Gulf Coast region,” International
under uncertainty,” Transportation Research E, vol. 43, no. 6, pp. Journal of Production Economics, vol. 131, no. 2, pp. 431–440,
737–754, 2007. 2011.
[5] S. V. Ukkusuri and W. F. Yushimito, “Location routing approach [22] B. M. Beamon and B. Balcik, “Performance measurement in
for the humanitarian prepositioning problem,” Transportation humanitarian relief chains,” International Journal of Public
Research Record, vol. 2089, pp. 18–25, 2008. Sector Management, vol. 21, no. 1, pp. 4–25, 2008.
[6] C. G. Rawls and M. A. Turnquist, “Pre-positioning of emer- [23] X. Prats, V. Puig, J. Quevedo, and F. Nejjari, “Multi-objective
gency supplies for disaster response,” Transportation Research optimisation for aircraft departure trajectories minimising
B, vol. 44, no. 4, pp. 521–534, 2010. noise annoyance,” Transportation Research C, vol. 18, no. 6, pp.
[7] C. G. Rawls and M. A. Turnquist, “Pre-positioning planning 975–989, 2010.
for emergency response with service quality constraints,” OR [24] L. A. Wolsey, Integer Programming, John Wiley & Sons, New
Spectrum, vol. 33, no. 3, pp. 481–498, 2011. York, NY, USA, 1998.
Mathematical Problems in Engineering 17

[25] A. Kokott and A. Lobel, “Lagrangian relaxations and subgradi-

ent methods for multiple-depot vehicle scheduling problems,”
Working Paper, Konrad-Zuse-Zentrum fur Information stech-
nik Berlin, 1996.
[26] F. Fumero, “A modified subgradient algorithm for Lagrangean
relaxation,” Computers & Operations Research, vol. 28, no. 1, pp.
33–52, 2001.
[27] National Disaster Reduction Commission of the P. R. China,
Integrated Analysis and Assessment on Wenchuan Earthquake
Disaster, Science Press, Beijing, China, 2008.
[28] Y. Yuan, “Impact of intensity and loss assessment following
the great Wenchuan Earthquake,” Earthquake Engineering and
Engineering Vibration, vol. 7, no. 3, pp. 247–254, 2008.
[29] L. Özdamar, E. Ekinci, and B. Küçükyazici, “Emergency
logistics planning in natural disasters,” Annals of Operations
Research, vol. 129, pp. 217–245, 2004.
[30] A. Haghani and S.-C. Oh, “Formulation and solution of a multi-
commodity, multi-modal network flow model for disaster relief
operations,” Transportation Research A, vol. 30, no. 3, pp. 231–
250, 1996.
Advances in Advances in Journal of Journal of
Operations Research
Hindawi Publishing Corporation
Decision Sciences
Hindawi Publishing Corporation
Applied Mathematics
Hindawi Publishing Corporation
Algebra
Hindawi Publishing Corporation
Probability and Statistics
Hindawi Publishing Corporation
http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014

The Scientific International Journal of

World Journal
Hindawi Publishing Corporation
Differential Equations
Hindawi Publishing Corporation
http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014

Submit your manuscripts at

http://www.hindawi.com

International Journal of Advances in

Combinatorics
Hindawi Publishing Corporation
Mathematical Physics
Hindawi Publishing Corporation
http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014

Journal of Journal of Mathematical Problems Abstract and Discrete Dynamics in

Complex Analysis
Hindawi Publishing Corporation
Mathematics
Hindawi Publishing Corporation
in Engineering
Hindawi Publishing Corporation
Applied Analysis
Hindawi Publishing Corporation
Nature and Society
Hindawi Publishing Corporation
http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014

International
Journal of Journal of
Mathematics and
Mathematical
Discrete Mathematics
Sciences

Journal of International Journal of Journal of

Hindawi Publishing Corporation Hindawi Publishing Corporation Volume 2014

Function Spaces
Hindawi Publishing Corporation
Stochastic Analysis
Hindawi Publishing Corporation
Optimization
Hindawi Publishing Corporation
http://www.hindawi.com Volume 2014 http://www.hindawi.com http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014