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

Authors: Yajie Liu; Bo Guo
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

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.


---

## Body

## 1. 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 helicopter-based 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. Section2 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 decomposition-based 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.


## 2. 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 two-stage 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 two-stage 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 mixed-integer 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 first-stage variables and 0-1 mixed-integer recourse variables [16–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-decomposition-based 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 location-allocation 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-decomposition-based 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 scenario-decomposition heuristic algorithm within an approximate lexicographic optimization frame are novel compared to existing work for emergency relief logistics planning.


## 3. Facility Location-Allocation Model

### 3.1. Notation Descriptions

Set and Indices
I: Set of candidate TFs; i
∈
I
J: Set of AAs; j
∈
J
K: Set of relief supplies; k
∈
K
K
R: Subset of relief supplies that can be classified as relief commodities, such as water, food, medicine and tent
K
S: Subset of relief supplies that can be classified as relief personnel, such as health personnel and relief workers
K
C: Subset of relief supplies that accommodates critical populations, that is, hospital beds
M: Set of helicopter modes; m
∈
M
M
G: Subset of general helicopter modes (utilized as delivering relief commodities, relief personnel, and critical population); m
∈
M
G
M
S: Subset of special helicopter modes (utilized as transporting relief personnel and critical population); m
∈
M
S
Ξ: Set of scenarios; ξ
∈
Ξ.Deterministic Parameters
U
i
k: Maximum mobilization levels for relief supply k at TF i (units)
h
-
k: Average heft of relief supply k  (k
∈
K
R, kilograms/units)
A
m: Total number of available helicopters of mode m
θ
-
m: Average capacity for relief supplies of helicopter mode m (kilograms/helicopter × trip)
η
-
m: Average capacity for relief personnel of helicopter mode m (persons/helicopter × trip)
ϑ
-
m: Average capacity for critical population of helicopter mode m (persons/helicopter × trip)
T
m: Average working time of helicopter mode m in the target responsive period (hours)
f
i
F: Fixed cost of establishing TF i ($)
f
m
H: Fixed cost of hiring a unit of helicopter mode m ($)
c
m
H: Operating cost for a unit of helicopter mode m ($/hour)
c
i
k: Mobilization cost for a unit of relief supply k at TF i ($/unit)
t
-
i
j
m: Average flight time between TF i and AA j by helicopter mode m (hours)
w
k: Weighting parameter for relief supply k in fill rate
ω
k
k
′
C: Weighting parameter on the penalty of discrepancy in fill rate for relief supplies k and k
′
L: A large number.Scenario-Specific Parameters, All under Scenario  ξ
d
j
k
ξ: Demand level for relief supply k  (
k
∈
K
R
∪
K
S
) at area j (units)
d
c
j
ξ: Amount of critical population at area j (persons)
p
ξ: Probability of scenario ξ occurrences.Scenario-Unrelated Decision Variables
x
i: Whether TF i is open (1 TF i open, 0 TF i closed)
s
i
k: Mobilization level of relief supply k at TF i (units)
y
i
m: Number of helicopter mode m that initially allocated to TF i (helicopters).Scenario-Specific Decision Variables (Units), All under Scenario  ξ
z
i
j
k
m
ξ: Amount of relief supply k  (
k
∈
K
R
∪
K
S
) delivered from TF i to AA j by helicopter mode m
z
c
j
i
m
ξ: Amount of critical population picked up by helicopter m from AA j to TF i
n
f
i
j
m
ξ: Number of trips from TF i to AA j by helicopter mode m
n
b
j
i
m
ξ: Number of trips from AA j to TF i by helicopter mode m
r
j
k
ξ: Fill rate of relief supply k at AA j
Δ
r
j
k
k
′
ξ: Fill rate difference between relief supplies k and k
′ at AA j.


### 3.2. 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:(1)
Max
⁡
f
1
=
∑
ξ
∈
Ξ
p
ξ
(
∑
k
∈
K
w
k
(
min
⁡
j
∈
J
{
r
j
k
ξ
}
)
)
-
∑
ξ
∈
Ξ
p
ξ
∑
j
∈
J
∑
k
∈
K
∑
k
′
∈
K
ω
k
k
′
C
Δ
r
j
k
k
′
ξ
,
(2)
Min
⁡
f
2
=
∑
i
∈
I
f
i
F
x
i
+
∑
i
∈
I
∑
k
∈
K
‍
c
i
k
s
i
k
+
∑
i
∈
I
∑
m
∈
M
f
m
H
y
i
m
+
∑
ξ
∈
Ξ
∑
i
∈
I
∑
j
∈
J
∑
m
∈
M
p
ξ
c
m
H
t
-
i
j
m
×
(
n
f
i
j
m
ξ
+
n
b
j
i
m
ξ
)
,
s.t. 								(3)
s
i
k
≤
L
x
i
,
∀
i
∈
I
,
k
∈
K
,
(4)
s
i
k
≤
U
i
k
,
∀
i
∈
I
,
k
∈
K
,
(5)
y
i
m
≤
A
m
x
i
,
∀
i
∈
I
,
m
∈
M
,
(6)
∑
i
∈
I
y
i
m
≤
A
m
,
∀
m
∈
M
,
(7)
∑
k
∈
K
R
(
h
-
k
z
i
j
k
m
ξ
)
θ
-
m
+
∑
k
′
∈
K
S
z
i
j
k
′
m
ξ
η
-
m
≤
n
f
i
j
m
ξ
,
∀
i
∈
I
,
j
∈
J
,
m
∈
M
G
,
ξ
∈
Ξ
,
(8)
∑
k
∈
K
S
z
i
j
k
m
ξ
≤
η
-
m
n
f
i
j
m
ξ
,
∀
i
∈
I
,
j
∈
J
,
m
∈
M
S
,
ξ
∈
Ξ
,
(9)
z
c
j
i
m
ξ
≤
ϑ
-
m
n
b
j
i
m
ξ
,
∀
i
∈
I
,
j
∈
J
,
m
∈
M
,
ξ
∈
Ξ
,
(10)
r
j
k
ξ
=
∑
i
∈
I
∑
m
∈
M
z
i
j
k
m
ξ
d
j
k
ξ
,
∀
j
∈
J
,
k
∈
K
R
∪
K
S
,
ξ
∈
Ξ
,
(11)
r
j
k
ξ
=
∑
i
∈
I
∑
m
∈
M
z
c
j
i
m
ξ
d
c
j
ξ
,
∀
j
∈
J
,
k
∈
K
C
,
ξ
∈
Ξ
,
(12)
Δ
r
j
k
k
′
ξ
≥
r
j
k
ξ
-
r
j
k
′
ξ
,
∀
j
∈
J
,
k
∈
K
,
k
′
∈
K
,
ξ
∈
Ξ
,
(13)
∑
i
∈
I
∑
m
∈
M
z
i
j
k
m
ξ
≤
d
j
k
ξ
,
∀
j
∈
J
,
k
∈
K
R
∪
K
S
,
ξ
∈
Ξ
,
(14)
∑
i
∈
I
∑
m
∈
M
z
c
j
i
m
ξ
≤
d
c
j
ξ
,
∀
j
∈
J
,
ξ
∈
Ξ
,
(15)
∑
i
∈
I
n
b
j
i
m
ξ
=
∑
i
∈
I
n
f
i
j
m
ξ
,
∀
j
∈
J
,
m
∈
M
,
ξ
∈
Ξ
,
(16)
∑
j
∈
J
n
f
i
j
m
ξ
≤
y
i
m
+
∑
j
∈
J
n
b
j
i
m
ξ
,
∀
i
∈
i
,
m
∈
M
,
ξ
∈
Ξ
,
(17)
∑
i
∈
I
∑
j
∈
J
t
-
i
j
m
(
n
f
i
j
m
ξ
+
n
b
j
i
m
ξ
)
≤
∑
i
∈
I
T
m
y
i
m
,
∀
m
∈
M
,
ξ
∈
Ξ
,
(18)
∑
j
∈
J
∑
m
∈
M
z
i
j
k
m
ξ
≤
s
i
k
,
∀
i
∈
I
,
k
∈
K
R
∪
K
S
,
ξ
∈
Ξ
,
(19)
∑
j
∈
J
∑
m
∈
M
z
c
j
i
m
ξ
≤
s
i
k
,
∀
i
∈
I
,
k
∈
K
C
,
ξ
∈
Ξ
,
(20)
x
i
∈
{
0,1
}
,
s
i
k
≥
0
,
y
i
m
≥
0
integer
,
∀
i
∈
I
,
k
∈
K
,
m
∈
M
,
(21)
z
i
j
k
m
ξ
≥
0
,
z
c
j
i
m
ξ
≥
0
,
n
f
i
j
m
ξ
≥
0
integer
,
n
b
j
i
m
ξ
≥
0
integer
,
0
≤
r
j
k
ξ
≤
1
,
Δ
r
j
k
k
′
ξ
≥
0
,
∀
i
∈
I
,
j
∈
J
,
k
∈
K
,
k
′
∈
K
,
m
∈
M
,
ξ
∈
Ξ
.
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 w
k
k
′
C is set larger if the synergistic requirement between k and k
′ is higher.Since this objective is explicitly nonlinear, the linear equivalent equations can be rewritten with the help of an auxiliary variableR
k
ξ as follows:
(22)
Max
⁡
f
1
=
∑
ξ
∈
Ξ
p
ξ
∑
k
∈
K
w
k
R
k
ξ
-
∑
ξ
∈
Ξ
p
ξ
∑
j
∈
J
∑
k
∈
K
‍
∑
k
′
∈
K
w
k
k
′
C
Δ
r
j
k
k
′
ξ
,
(23)
s
.
t
.
R
k
ξ
≤
r
j
k
ξ
,
∀
j
∈
J
,
k
∈
K
,
ξ
∈
Ξ
,
(24)
R
k
ξ
≥
0
,
∀
k
∈
K
,
ξ
∈
Ξ
.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 k and k
′ at AA j. 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 scenario-specific 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 scenario-specific decision variables, respectively.


## 3.1. Notation Descriptions

Set and Indices
I: Set of candidate TFs; i
∈
I
J: Set of AAs; j
∈
J
K: Set of relief supplies; k
∈
K
K
R: Subset of relief supplies that can be classified as relief commodities, such as water, food, medicine and tent
K
S: Subset of relief supplies that can be classified as relief personnel, such as health personnel and relief workers
K
C: Subset of relief supplies that accommodates critical populations, that is, hospital beds
M: Set of helicopter modes; m
∈
M
M
G: Subset of general helicopter modes (utilized as delivering relief commodities, relief personnel, and critical population); m
∈
M
G
M
S: Subset of special helicopter modes (utilized as transporting relief personnel and critical population); m
∈
M
S
Ξ: Set of scenarios; ξ
∈
Ξ.Deterministic Parameters
U
i
k: Maximum mobilization levels for relief supply k at TF i (units)
h
-
k: Average heft of relief supply k  (k
∈
K
R, kilograms/units)
A
m: Total number of available helicopters of mode m
θ
-
m: Average capacity for relief supplies of helicopter mode m (kilograms/helicopter × trip)
η
-
m: Average capacity for relief personnel of helicopter mode m (persons/helicopter × trip)
ϑ
-
m: Average capacity for critical population of helicopter mode m (persons/helicopter × trip)
T
m: Average working time of helicopter mode m in the target responsive period (hours)
f
i
F: Fixed cost of establishing TF i ($)
f
m
H: Fixed cost of hiring a unit of helicopter mode m ($)
c
m
H: Operating cost for a unit of helicopter mode m ($/hour)
c
i
k: Mobilization cost for a unit of relief supply k at TF i ($/unit)
t
-
i
j
m: Average flight time between TF i and AA j by helicopter mode m (hours)
w
k: Weighting parameter for relief supply k in fill rate
ω
k
k
′
C: Weighting parameter on the penalty of discrepancy in fill rate for relief supplies k and k
′
L: A large number.Scenario-Specific Parameters, All under Scenario  ξ
d
j
k
ξ: Demand level for relief supply k  (
k
∈
K
R
∪
K
S
) at area j (units)
d
c
j
ξ: Amount of critical population at area j (persons)
p
ξ: Probability of scenario ξ occurrences.Scenario-Unrelated Decision Variables
x
i: Whether TF i is open (1 TF i open, 0 TF i closed)
s
i
k: Mobilization level of relief supply k at TF i (units)
y
i
m: Number of helicopter mode m that initially allocated to TF i (helicopters).Scenario-Specific Decision Variables (Units), All under Scenario  ξ
z
i
j
k
m
ξ: Amount of relief supply k  (
k
∈
K
R
∪
K
S
) delivered from TF i to AA j by helicopter mode m
z
c
j
i
m
ξ: Amount of critical population picked up by helicopter m from AA j to TF i
n
f
i
j
m
ξ: Number of trips from TF i to AA j by helicopter mode m
n
b
j
i
m
ξ: Number of trips from AA j to TF i by helicopter mode m
r
j
k
ξ: Fill rate of relief supply k at AA j
Δ
r
j
k
k
′
ξ: Fill rate difference between relief supplies k and k
′ at AA j.


## 3.2. 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:(1)
Max
⁡
f
1
=
∑
ξ
∈
Ξ
p
ξ
(
∑
k
∈
K
w
k
(
min
⁡
j
∈
J
{
r
j
k
ξ
}
)
)
-
∑
ξ
∈
Ξ
p
ξ
∑
j
∈
J
∑
k
∈
K
∑
k
′
∈
K
ω
k
k
′
C
Δ
r
j
k
k
′
ξ
,
(2)
Min
⁡
f
2
=
∑
i
∈
I
f
i
F
x
i
+
∑
i
∈
I
∑
k
∈
K
‍
c
i
k
s
i
k
+
∑
i
∈
I
∑
m
∈
M
f
m
H
y
i
m
+
∑
ξ
∈
Ξ
∑
i
∈
I
∑
j
∈
J
∑
m
∈
M
p
ξ
c
m
H
t
-
i
j
m
×
(
n
f
i
j
m
ξ
+
n
b
j
i
m
ξ
)
,
s.t. 								(3)
s
i
k
≤
L
x
i
,
∀
i
∈
I
,
k
∈
K
,
(4)
s
i
k
≤
U
i
k
,
∀
i
∈
I
,
k
∈
K
,
(5)
y
i
m
≤
A
m
x
i
,
∀
i
∈
I
,
m
∈
M
,
(6)
∑
i
∈
I
y
i
m
≤
A
m
,
∀
m
∈
M
,
(7)
∑
k
∈
K
R
(
h
-
k
z
i
j
k
m
ξ
)
θ
-
m
+
∑
k
′
∈
K
S
z
i
j
k
′
m
ξ
η
-
m
≤
n
f
i
j
m
ξ
,
∀
i
∈
I
,
j
∈
J
,
m
∈
M
G
,
ξ
∈
Ξ
,
(8)
∑
k
∈
K
S
z
i
j
k
m
ξ
≤
η
-
m
n
f
i
j
m
ξ
,
∀
i
∈
I
,
j
∈
J
,
m
∈
M
S
,
ξ
∈
Ξ
,
(9)
z
c
j
i
m
ξ
≤
ϑ
-
m
n
b
j
i
m
ξ
,
∀
i
∈
I
,
j
∈
J
,
m
∈
M
,
ξ
∈
Ξ
,
(10)
r
j
k
ξ
=
∑
i
∈
I
∑
m
∈
M
z
i
j
k
m
ξ
d
j
k
ξ
,
∀
j
∈
J
,
k
∈
K
R
∪
K
S
,
ξ
∈
Ξ
,
(11)
r
j
k
ξ
=
∑
i
∈
I
∑
m
∈
M
z
c
j
i
m
ξ
d
c
j
ξ
,
∀
j
∈
J
,
k
∈
K
C
,
ξ
∈
Ξ
,
(12)
Δ
r
j
k
k
′
ξ
≥
r
j
k
ξ
-
r
j
k
′
ξ
,
∀
j
∈
J
,
k
∈
K
,
k
′
∈
K
,
ξ
∈
Ξ
,
(13)
∑
i
∈
I
∑
m
∈
M
z
i
j
k
m
ξ
≤
d
j
k
ξ
,
∀
j
∈
J
,
k
∈
K
R
∪
K
S
,
ξ
∈
Ξ
,
(14)
∑
i
∈
I
∑
m
∈
M
z
c
j
i
m
ξ
≤
d
c
j
ξ
,
∀
j
∈
J
,
ξ
∈
Ξ
,
(15)
∑
i
∈
I
n
b
j
i
m
ξ
=
∑
i
∈
I
n
f
i
j
m
ξ
,
∀
j
∈
J
,
m
∈
M
,
ξ
∈
Ξ
,
(16)
∑
j
∈
J
n
f
i
j
m
ξ
≤
y
i
m
+
∑
j
∈
J
n
b
j
i
m
ξ
,
∀
i
∈
i
,
m
∈
M
,
ξ
∈
Ξ
,
(17)
∑
i
∈
I
∑
j
∈
J
t
-
i
j
m
(
n
f
i
j
m
ξ
+
n
b
j
i
m
ξ
)
≤
∑
i
∈
I
T
m
y
i
m
,
∀
m
∈
M
,
ξ
∈
Ξ
,
(18)
∑
j
∈
J
∑
m
∈
M
z
i
j
k
m
ξ
≤
s
i
k
,
∀
i
∈
I
,
k
∈
K
R
∪
K
S
,
ξ
∈
Ξ
,
(19)
∑
j
∈
J
∑
m
∈
M
z
c
j
i
m
ξ
≤
s
i
k
,
∀
i
∈
I
,
k
∈
K
C
,
ξ
∈
Ξ
,
(20)
x
i
∈
{
0,1
}
,
s
i
k
≥
0
,
y
i
m
≥
0
integer
,
∀
i
∈
I
,
k
∈
K
,
m
∈
M
,
(21)
z
i
j
k
m
ξ
≥
0
,
z
c
j
i
m
ξ
≥
0
,
n
f
i
j
m
ξ
≥
0
integer
,
n
b
j
i
m
ξ
≥
0
integer
,
0
≤
r
j
k
ξ
≤
1
,
Δ
r
j
k
k
′
ξ
≥
0
,
∀
i
∈
I
,
j
∈
J
,
k
∈
K
,
k
′
∈
K
,
m
∈
M
,
ξ
∈
Ξ
.
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 w
k
k
′
C is set larger if the synergistic requirement between k and k
′ is higher.Since this objective is explicitly nonlinear, the linear equivalent equations can be rewritten with the help of an auxiliary variableR
k
ξ as follows:
(22)
Max
⁡
f
1
=
∑
ξ
∈
Ξ
p
ξ
∑
k
∈
K
w
k
R
k
ξ
-
∑
ξ
∈
Ξ
p
ξ
∑
j
∈
J
∑
k
∈
K
‍
∑
k
′
∈
K
w
k
k
′
C
Δ
r
j
k
k
′
ξ
,
(23)
s
.
t
.
R
k
ξ
≤
r
j
k
ξ
,
∀
j
∈
J
,
k
∈
K
,
ξ
∈
Ξ
,
(24)
R
k
ξ
≥
0
,
∀
k
∈
K
,
ξ
∈
Ξ
.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 k and k
′ at AA j. 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 scenario-specific 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 scenario-specific decision variables, respectively.


## 4. 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:
(25)
(
P
1
)
:
f
1
*
=
Max
⁡
f
1
s.t
.
:
{
(
3
)
–
(
21
)
,
(
23
)
-
(
24
)
,
(26)
(
P
2
)
:
f
2
*
=
Min
⁡
f
2
s.t
.
:
{
(
3
)
–
(
21
)
,
(
23
)
-
(
24
)
,
f
1
≥
f
1
*
,
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 f
1
≥
f
1
* to (P2), a new model modified from (P2), which is named as (P3) here utilized to substitute (P2), is defined as follows:
(27)
(
P
3
)
:
f
2
*
=
Min
⁡
f
2
,
s.t
.
:
{
(
3
)
–
(
11
)
,
(
15
)
–
(
20
)
,
r
j
k
ξ
≥
(
1
-
α
)
r
~
j
k
ξ
,
∀
j
∈
J
,
k
∈
K
,
ξ
∈
Ξ
,
(
21
)
∖
{
Δ
r
j
k
k
′
ξ
}
,
where r
~
j
k
ξ is the solution of r
j
k
ξ 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 r
j
k
ξ
≥
(
1
-
α
)
r
~
j
k
ξ to (P3), and the value of f
1 can also be easily determined resorting to the solution of r
~
j
k
ξ 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.


### 4.1. Solution to (P1)

Since (P1) does not contain constraints related to expenditure, a scenario-decomposition-based heuristic approach was used to solve this problem. The details of this approach are as follows.Step 1.
Create a new problem named (P
1
′) by duplicating x
i
ξ of x
i, s
i
k
ξ of s
i
k, and y
i
m
ξ of y
i
m for all scenarios ξ
∈
Ξ on (P1).
As a result, all decision variables in (P
1
′) become scenario specific, meaning that (P
1
′) can be decomposed into independent subproblems by scenario.Step 2.
Decompose (P
1
′) into |
Ξ
| subproblems by scenario and solve each subproblem using a linear programming solver, denoting the value of x
i
ξ as x
^
i
ξ, the value of s
i
k
ξ as s
^
i
k
ξ, and the value of y
i
m
ξ as y
^
i
m
ξ, respectively. Let x
i
=
max
⁡
ξ
∈
Ξ
{
x
^
i
ξ
} and s
i
k
=
max
⁡
ξ
∈
Ξ
{
s
^
i
k
ξ
}.After these assignments, it can be easily deduced that the constraint sets related tox
i and s
i
k in (P
1
′) can be satisfied. However, to y
i
m, 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 y
i
m a feasible solution.Step 3.
For each helicopter modem, let y
~
i
m
=
∑
ξ
p
ξ
y
^
i
m
ξ  (
∀
i
∈
I
), y
^
i
m
=
⌊
y
~
i
m
⌋  (
∀
i
∈
I
), and Δ
y
i
m
=
y
~
i
m
-
y
^
i
m  (
∀
i
∈
I
), then sort Δ
y
i
m in nonincreasing order according to their values, and move the first elements with the number of min
⁡
{
A
m
,
⌈
∑
i
y
~
i
m
⌉
}
-
∑
i
y
^
i
m to a set (i.e., Q); to each element in Q, supposing its index on TF is i, let y
^
i
m
=
y
^
i
m
+
1. Finally, let y
i
m
=
y
^
i
m.Obviously, after executing this step, the constraint sets relating toy
i
m 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 r
j
k
ξ in (P1) with given values on all scenario-unrelated variables.Step 4.
Add the constraint setsx
i
=
max
⁡
ξ
∈
Ξ
{
x
^
i
ξ
}  (
∀
i
∈
I
), s
i
k
=
max
⁡
ξ
∈
Ξ
{
s
^
i
k
ξ
}  (
∀
i
∈
I
,
k
∈
K
), and y
i
m
=
y
^
i
m  (
∀
i
∈
I
,
m
∈
M
) to (P1), then decompose (P1) into subproblems by scenario, and solve each subproblem, denoting the value of r
j
k
ξ as r
~
j
k
ξ.


### 4.2. Solution to (P3)

#### 4.2.1. Model Transformation

To decompose (P3) effectively, a new model named (P
3
′) is defined by duplicating x
i
ξ of x
i, s
i
k
ξ of s
i
k, and y
i
m
ξ of y
i
m for all scenarios ξ
∈
Ξ on the model (P3); then the following three constraint sets are added to (P
3
′):
(28)
(
1
-
p
1
)
x
i
1
=
∑
ξ
∈
Ξ
∖
{
1
}
p
ξ
x
i
ξ
,
∀
i
∈
I
,
(29)
s
i
k
1
-
s
i
k
ξ
=
0
,
∀
i
∈
I
,
k
∈
K
,
ξ
∈
Ξ
∖
{
1
}
,
(30)
y
i
m
1
-
y
i
m
ξ
=
0
,
∀
i
∈
I
,
m
∈
M
,
ξ
∈
Ξ
∖
{
1
}
,
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 i under all scenarios are identical, that is, x
i
1
=
x
i
2
=
⋯
=
x
i
|
Ξ
|, and (29) and (30) state similar limitations on s
i
k
ξ and y
i
m
ξ, respectively. From (28)–(30), it can be easily concluded that the two problems, (P3) and (P
3
′), are equivalent.(P
3
′) can be decomposed into independent subproblems by scenario if constraint sets (28)–(30) are relaxed. Therefore, the Lagrangian relaxation of (P
3
′) with respect to the constraint sets (28)–(30) can be described as follows:
(31)
f
2
LR
(
u
,
w
,
v
)
=
min
⁡
{
f
2
+
∑
ξ
∈
Ξ
∖
{
1
}
(
p
ξ
∑
i
∈
I
‍
u
i
(
x
i
1
-
x
i
ξ
)
‍
+
∑
i
∈
I
‍
∑
k
∈
K
‍
w
i
k
ξ
(
s
i
k
1
-
s
i
k
ξ
)
+
∑
i
∈
I
‍
∑
m
∈
M
‍
v
i
m
ξ
(
y
i
m
1
-
y
i
m
ξ
)
)
∑
ξ
∈
Ξ
∖
{
1
}
(
p
ξ
∑
i
∈
I
‍
u
i
(
x
i
1
-
x
i
ξ
)
‍
}
s
.
t
.
:
all
the
constraints
of
(
P
3
′
)
∖
constraints
(
28
)
–
(
30
)
,
where f
2 is the objective function of (P3), and the vectors u, w, and v 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:
(32)
z
LD
=
max
⁡
u
,
w
,
v
f
2
LR
(
u
,
w
,
v
)
.


#### 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.Step 0 (initialization).
t
=
0 (iteration counter),
(33)
u
i
(
t
)
=
0
,
w
i
k
ξ
,
(
t
)
=
0
,
v
i
m
ξ
,
(
t
)
=
0
.
LetL
R
(
t
)
=
f
2
LR
(
u
(
t
)
,
w
(
t
)
,
v
(
t
)
) and U
B
(
t
)
= upper bound value of the objective of (P3). Initialize L
R
(
t
)
=
-
∞ and U
B
(
t
)
=
+
∞.Step 1 (solve each subproblem).
With givenu
(
t
), w
(
t
), and v
(
t
), decompose f
2
LR
(
u
,
w
,
v
) into subproblems by scenario, solve each subproblem using a linear programming solver, and add its objective value to L
R
(
t
+
1
).
Keepx
i
ξ
,
(
t
), s
i
k
ξ
,
(
t
), and y
i
m
ξ
,
(
t
)(the values of x
i
ξ, s
i
k
ξ, and y
i
m
ξ gotten in iteration t), respectively.
∀
i
∈
I, if (x
i
1
,
(
t
)
=
x
i
2
,
(
t
)
=
⋯
=
x
i
|
Ξ
|
,
(
t
)), then add constraint set x
i
ξ
=
x
i
1
,
(
t
)  (
∀
ξ
∈
Ξ
) to f
2
L
R
(
u
,
w
,
v
).
∀
i
∈
I and ∀
k
∈
K, if (s
i
k
1
,
(
t
)
=
s
i
k
2
,
(
t
)
=
⋯
=
s
i
k
|
Ξ
|
,
(
t
)), then add constraint set s
i
k
ξ
=
s
i
k
1
,
(
t
)  (
∀
ξ
∈
Ξ
) to f
2
LR
(
u
,
w
,
v
).
∀
i
∈
I and ∀
m
∈
M, if (y
i
m
1
,
(
t
)
=
y
i
m
2
,
(
t
)
=
⋯
=
y
i
m
|
Ξ
|
,
(
t
)), then add constraint set y
i
m
ξ
=
y
i
m
1
,
(
t
)  (
∀
ξ
∈
Ξ
) to f
2
LR
(
u
,
w
,
v
).Step 2 (compute subgradients).
(34)
subgrad
_
u
i
(
t
+
1
)
=
x
i
1
,
(
t
)
-
∑
ξ
∈
Ξ
‍
p
ξ
x
i
ξ
,
(
t
)
,
subgrad
_
w
i
k
ξ
,
(
t
+
1
)
=
s
i
k
1
,
(
t
)
-
s
i
k
ξ
,
(
t
)
,
subgrad
_
v
i
m
ξ
,
(
t
+
1
)
=
y
i
m
1
,
(
t
)
-
y
i
m
ξ
,
(
t
)
.Step 3 (update upper bound).
(35)
x
i
(
t
)
=
max
⁡
ξ
∈
Ξ
{
x
i
ξ
,
(
t
)
}
,
s
i
k
(
t
)
=
max
⁡
ξ
∈
Ξ
{
s
i
k
ξ
,
(
t
)
}
.
Compute the values ofy
i
m
(
t
)  (
∀
i
∈
I
,
m
∈
M
) following the adjustment process of Step 3 in Section 4.1 based on the values of y
i
m
ξ
,
(
t
)  (
∀
ξ
∈
Ξ
), and then solve (P3) where the scenario-unrelated variables, x
i, s
i
k, and y
i
m, are assigned to the values of x
i
(
t
), s
i
k
(
t
), and y
i
m
(
t
), respectively. Supposing the objective of (P3) is UB, let U
B
(
t
+
1
)
=
min
⁡
{
UB
,
U
B
(
t
)
}.
Record the values ofx
i, s
i
k, and y
i
m that have been utilized to generate the updated upper bound as x
^
i, s
^
i
k, and y
^
i
m, respectively.Step 4 (identify new step direction).
(36)
stepdirc
_
u
i
(
t
+
1
)
=
β
1
(
subgrad
_
u
i
(
t
+
1
)
)
+
β
2
(
subgrad
_
u
i
(
t
)
)
+
β
3
(
subgrad
_
u
i
(
t
-
1
)
)
,
stepdirc
_
w
i
k
ξ
,
(
t
+
1
)
=
β
1
(
subgrad
_
w
i
k
ξ
,
(
t
+
1
)
)
+
β
2
(
subgrad
_
w
i
k
ξ
,
(
t
)
)
+
β
3
(
subgrad
_
w
i
k
ξ
,
(
t
-
1
)
)
,
stepdirc
_
v
i
m
ξ
,
(
t
+
1
)
=
β
1
(
subgrad
_
v
i
m
ξ
,
(
t
+
1
)
)
+
β
2
(
subgrad
_
v
i
m
ξ
,
(
t
)
)
+
β
3
(
subgrad
_
v
i
m
ξ
,
(
t
-
1
)
)
,
where
β
1
=
0.6
,
β
2
=
0.3
,
β
3
=
0.1
.Step 5 (update multipliers).
(37)
sum
_
of
_
squares
=
∑
i
∈
I
(
subgrad
_
u
i
(
t
+
1
)
)
2
+
∑
i
∈
I
∑
k
∈
K
∑
ξ
∈
Ξ
(
subgrad
_
w
i
k
ξ
,
(
t
+
1
)
)
2
+
∑
i
∈
I
∑
m
∈
M
∑
ξ
∈
Ξ
(
subgrad
_
v
i
m
ξ
,
(
t
+
1
)
)
2
,
if (sum
_
of
_
squares
=
0)
goto Step6
else(38)
u
i
(
t
+
1
)
=
u
i
(
t
)
+
stepdirc
_
u
i
(
t
+
1
)
2.0
×
η
t
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
)
sum
_
of
_
squares
,
w
i
k
ξ
,
(
t
+
1
)
=
w
i
k
ξ
,
(
t
)
+
stepdirc
_
w
i
k
ξ
,
(
t
+
1
)
2.0
×
η
t
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
)
sum
_
of
_
squares
,
v
i
m
ξ
,
(
t
+
1
)
=
v
i
m
ξ
,
(
t
)
+
stepdirc
_
v
i
m
ξ
,
(
t
+
1
)
2.0
×
η
t
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
)
sum
_
of
_
squares
,
where η is a user-defined parameter and the value η
=
0.99 is used in the experiments.Step 6 (stopping criteria).
(39)
if
(
(
sum
_
of
_
squares
=
0
)
∥
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
U
B
(
t
+
1
)
≤
ε
)
∥
×
∥
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
U
B
(
t
+
1
)
≤
ε
)
∥
(
t
≥
N
2
)
)
stop;
else
t
=
t
+
1, goto Step 1
where the value ofN
2 is set to 200 and the value of ε is set to 0.001.
In Step0, 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 x
i, s
i
k, and y
i
m 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 (P
3
′) (which is equivalent to (P3)); that is, the values of x
i
ξ, s
i
k
ξ, and y
i
m
ξ 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.


#### 4.2.3. 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.(i)
Add the constraint setsx
i
=
x
^
i  (
∀
i
∈
I
), s
i
k
=
s
^
i
k  (
∀
i
∈
I
,
k
∈
K
), and y
i
m
=
y
^
i
m  (
∀
i
∈
I
,
m
∈
M
) to (P3), in which x
^
i, s
^
i
k, and y
^
i
m are the recorded values of x
i, s
i
k, and y
i
m in generating the UB in Step 3 of the subgradient algorithm. After adding these three constraint sets, (P3) merely contains scenario-specific decision variables. In this way, (P3) can further be decomposed into independent subproblems by scenario;
(ii)
Decompose (P3) into|
Ξ
| subproblems by scenario, and then solve each subproblem using a linear programming solver. Supposing the values of n
f
i
j
m
ξ and n
b
j
i
m
ξ are n
f
^
i
j
m
ξ and n
b
^
j
i
m
ξ, respectively, the objective value of (P3) equals ∑
i
∈
I
f
i
F
x
^
i
+
∑
i
∈
I
∑
k
∈
K
c
i
k
s
^
i
k
+
∑
i
∈
I
∑
m
∈
M
f
m
H
y
^
i
m
ξ
+
∑
ξ
∈
Ξ
∑
i
∈
I
∑
j
∈
J
∑
m
∈
M
p
ξ
c
m
H
t
-
i
j
m
(
n
f
^
i
j
m
ξ
+
n
b
^
j
i
m
ξ
).


## 4.1. Solution to (P1)

Since (P1) does not contain constraints related to expenditure, a scenario-decomposition-based heuristic approach was used to solve this problem. The details of this approach are as follows.Step 1.
Create a new problem named (P
1
′) by duplicating x
i
ξ of x
i, s
i
k
ξ of s
i
k, and y
i
m
ξ of y
i
m for all scenarios ξ
∈
Ξ on (P1).
As a result, all decision variables in (P
1
′) become scenario specific, meaning that (P
1
′) can be decomposed into independent subproblems by scenario.Step 2.
Decompose (P
1
′) into |
Ξ
| subproblems by scenario and solve each subproblem using a linear programming solver, denoting the value of x
i
ξ as x
^
i
ξ, the value of s
i
k
ξ as s
^
i
k
ξ, and the value of y
i
m
ξ as y
^
i
m
ξ, respectively. Let x
i
=
max
⁡
ξ
∈
Ξ
{
x
^
i
ξ
} and s
i
k
=
max
⁡
ξ
∈
Ξ
{
s
^
i
k
ξ
}.After these assignments, it can be easily deduced that the constraint sets related tox
i and s
i
k in (P
1
′) can be satisfied. However, to y
i
m, 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 y
i
m a feasible solution.Step 3.
For each helicopter modem, let y
~
i
m
=
∑
ξ
p
ξ
y
^
i
m
ξ  (
∀
i
∈
I
), y
^
i
m
=
⌊
y
~
i
m
⌋  (
∀
i
∈
I
), and Δ
y
i
m
=
y
~
i
m
-
y
^
i
m  (
∀
i
∈
I
), then sort Δ
y
i
m in nonincreasing order according to their values, and move the first elements with the number of min
⁡
{
A
m
,
⌈
∑
i
y
~
i
m
⌉
}
-
∑
i
y
^
i
m to a set (i.e., Q); to each element in Q, supposing its index on TF is i, let y
^
i
m
=
y
^
i
m
+
1. Finally, let y
i
m
=
y
^
i
m.Obviously, after executing this step, the constraint sets relating toy
i
m 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 r
j
k
ξ in (P1) with given values on all scenario-unrelated variables.Step 4.
Add the constraint setsx
i
=
max
⁡
ξ
∈
Ξ
{
x
^
i
ξ
}  (
∀
i
∈
I
), s
i
k
=
max
⁡
ξ
∈
Ξ
{
s
^
i
k
ξ
}  (
∀
i
∈
I
,
k
∈
K
), and y
i
m
=
y
^
i
m  (
∀
i
∈
I
,
m
∈
M
) to (P1), then decompose (P1) into subproblems by scenario, and solve each subproblem, denoting the value of r
j
k
ξ as r
~
j
k
ξ.


## 4.2. Solution to (P3)

### 4.2.1. Model Transformation

To decompose (P3) effectively, a new model named (P
3
′) is defined by duplicating x
i
ξ of x
i, s
i
k
ξ of s
i
k, and y
i
m
ξ of y
i
m for all scenarios ξ
∈
Ξ on the model (P3); then the following three constraint sets are added to (P
3
′):
(28)
(
1
-
p
1
)
x
i
1
=
∑
ξ
∈
Ξ
∖
{
1
}
p
ξ
x
i
ξ
,
∀
i
∈
I
,
(29)
s
i
k
1
-
s
i
k
ξ
=
0
,
∀
i
∈
I
,
k
∈
K
,
ξ
∈
Ξ
∖
{
1
}
,
(30)
y
i
m
1
-
y
i
m
ξ
=
0
,
∀
i
∈
I
,
m
∈
M
,
ξ
∈
Ξ
∖
{
1
}
,
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 i under all scenarios are identical, that is, x
i
1
=
x
i
2
=
⋯
=
x
i
|
Ξ
|, and (29) and (30) state similar limitations on s
i
k
ξ and y
i
m
ξ, respectively. From (28)–(30), it can be easily concluded that the two problems, (P3) and (P
3
′), are equivalent.(P
3
′) can be decomposed into independent subproblems by scenario if constraint sets (28)–(30) are relaxed. Therefore, the Lagrangian relaxation of (P
3
′) with respect to the constraint sets (28)–(30) can be described as follows:
(31)
f
2
LR
(
u
,
w
,
v
)
=
min
⁡
{
f
2
+
∑
ξ
∈
Ξ
∖
{
1
}
(
p
ξ
∑
i
∈
I
‍
u
i
(
x
i
1
-
x
i
ξ
)
‍
+
∑
i
∈
I
‍
∑
k
∈
K
‍
w
i
k
ξ
(
s
i
k
1
-
s
i
k
ξ
)
+
∑
i
∈
I
‍
∑
m
∈
M
‍
v
i
m
ξ
(
y
i
m
1
-
y
i
m
ξ
)
)
∑
ξ
∈
Ξ
∖
{
1
}
(
p
ξ
∑
i
∈
I
‍
u
i
(
x
i
1
-
x
i
ξ
)
‍
}
s
.
t
.
:
all
the
constraints
of
(
P
3
′
)
∖
constraints
(
28
)
–
(
30
)
,
where f
2 is the objective function of (P3), and the vectors u, w, and v 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:
(32)
z
LD
=
max
⁡
u
,
w
,
v
f
2
LR
(
u
,
w
,
v
)
.


### 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.Step 0 (initialization).
t
=
0 (iteration counter),
(33)
u
i
(
t
)
=
0
,
w
i
k
ξ
,
(
t
)
=
0
,
v
i
m
ξ
,
(
t
)
=
0
.
LetL
R
(
t
)
=
f
2
LR
(
u
(
t
)
,
w
(
t
)
,
v
(
t
)
) and U
B
(
t
)
= upper bound value of the objective of (P3). Initialize L
R
(
t
)
=
-
∞ and U
B
(
t
)
=
+
∞.Step 1 (solve each subproblem).
With givenu
(
t
), w
(
t
), and v
(
t
), decompose f
2
LR
(
u
,
w
,
v
) into subproblems by scenario, solve each subproblem using a linear programming solver, and add its objective value to L
R
(
t
+
1
).
Keepx
i
ξ
,
(
t
), s
i
k
ξ
,
(
t
), and y
i
m
ξ
,
(
t
)(the values of x
i
ξ, s
i
k
ξ, and y
i
m
ξ gotten in iteration t), respectively.
∀
i
∈
I, if (x
i
1
,
(
t
)
=
x
i
2
,
(
t
)
=
⋯
=
x
i
|
Ξ
|
,
(
t
)), then add constraint set x
i
ξ
=
x
i
1
,
(
t
)  (
∀
ξ
∈
Ξ
) to f
2
L
R
(
u
,
w
,
v
).
∀
i
∈
I and ∀
k
∈
K, if (s
i
k
1
,
(
t
)
=
s
i
k
2
,
(
t
)
=
⋯
=
s
i
k
|
Ξ
|
,
(
t
)), then add constraint set s
i
k
ξ
=
s
i
k
1
,
(
t
)  (
∀
ξ
∈
Ξ
) to f
2
LR
(
u
,
w
,
v
).
∀
i
∈
I and ∀
m
∈
M, if (y
i
m
1
,
(
t
)
=
y
i
m
2
,
(
t
)
=
⋯
=
y
i
m
|
Ξ
|
,
(
t
)), then add constraint set y
i
m
ξ
=
y
i
m
1
,
(
t
)  (
∀
ξ
∈
Ξ
) to f
2
LR
(
u
,
w
,
v
).Step 2 (compute subgradients).
(34)
subgrad
_
u
i
(
t
+
1
)
=
x
i
1
,
(
t
)
-
∑
ξ
∈
Ξ
‍
p
ξ
x
i
ξ
,
(
t
)
,
subgrad
_
w
i
k
ξ
,
(
t
+
1
)
=
s
i
k
1
,
(
t
)
-
s
i
k
ξ
,
(
t
)
,
subgrad
_
v
i
m
ξ
,
(
t
+
1
)
=
y
i
m
1
,
(
t
)
-
y
i
m
ξ
,
(
t
)
.Step 3 (update upper bound).
(35)
x
i
(
t
)
=
max
⁡
ξ
∈
Ξ
{
x
i
ξ
,
(
t
)
}
,
s
i
k
(
t
)
=
max
⁡
ξ
∈
Ξ
{
s
i
k
ξ
,
(
t
)
}
.
Compute the values ofy
i
m
(
t
)  (
∀
i
∈
I
,
m
∈
M
) following the adjustment process of Step 3 in Section 4.1 based on the values of y
i
m
ξ
,
(
t
)  (
∀
ξ
∈
Ξ
), and then solve (P3) where the scenario-unrelated variables, x
i, s
i
k, and y
i
m, are assigned to the values of x
i
(
t
), s
i
k
(
t
), and y
i
m
(
t
), respectively. Supposing the objective of (P3) is UB, let U
B
(
t
+
1
)
=
min
⁡
{
UB
,
U
B
(
t
)
}.
Record the values ofx
i, s
i
k, and y
i
m that have been utilized to generate the updated upper bound as x
^
i, s
^
i
k, and y
^
i
m, respectively.Step 4 (identify new step direction).
(36)
stepdirc
_
u
i
(
t
+
1
)
=
β
1
(
subgrad
_
u
i
(
t
+
1
)
)
+
β
2
(
subgrad
_
u
i
(
t
)
)
+
β
3
(
subgrad
_
u
i
(
t
-
1
)
)
,
stepdirc
_
w
i
k
ξ
,
(
t
+
1
)
=
β
1
(
subgrad
_
w
i
k
ξ
,
(
t
+
1
)
)
+
β
2
(
subgrad
_
w
i
k
ξ
,
(
t
)
)
+
β
3
(
subgrad
_
w
i
k
ξ
,
(
t
-
1
)
)
,
stepdirc
_
v
i
m
ξ
,
(
t
+
1
)
=
β
1
(
subgrad
_
v
i
m
ξ
,
(
t
+
1
)
)
+
β
2
(
subgrad
_
v
i
m
ξ
,
(
t
)
)
+
β
3
(
subgrad
_
v
i
m
ξ
,
(
t
-
1
)
)
,
where
β
1
=
0.6
,
β
2
=
0.3
,
β
3
=
0.1
.Step 5 (update multipliers).
(37)
sum
_
of
_
squares
=
∑
i
∈
I
(
subgrad
_
u
i
(
t
+
1
)
)
2
+
∑
i
∈
I
∑
k
∈
K
∑
ξ
∈
Ξ
(
subgrad
_
w
i
k
ξ
,
(
t
+
1
)
)
2
+
∑
i
∈
I
∑
m
∈
M
∑
ξ
∈
Ξ
(
subgrad
_
v
i
m
ξ
,
(
t
+
1
)
)
2
,
if (sum
_
of
_
squares
=
0)
goto Step6
else(38)
u
i
(
t
+
1
)
=
u
i
(
t
)
+
stepdirc
_
u
i
(
t
+
1
)
2.0
×
η
t
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
)
sum
_
of
_
squares
,
w
i
k
ξ
,
(
t
+
1
)
=
w
i
k
ξ
,
(
t
)
+
stepdirc
_
w
i
k
ξ
,
(
t
+
1
)
2.0
×
η
t
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
)
sum
_
of
_
squares
,
v
i
m
ξ
,
(
t
+
1
)
=
v
i
m
ξ
,
(
t
)
+
stepdirc
_
v
i
m
ξ
,
(
t
+
1
)
2.0
×
η
t
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
)
sum
_
of
_
squares
,
where η is a user-defined parameter and the value η
=
0.99 is used in the experiments.Step 6 (stopping criteria).
(39)
if
(
(
sum
_
of
_
squares
=
0
)
∥
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
U
B
(
t
+
1
)
≤
ε
)
∥
×
∥
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
U
B
(
t
+
1
)
≤
ε
)
∥
(
t
≥
N
2
)
)
stop;
else
t
=
t
+
1, goto Step 1
where the value ofN
2 is set to 200 and the value of ε is set to 0.001.
In Step0, 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 x
i, s
i
k, and y
i
m 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 (P
3
′) (which is equivalent to (P3)); that is, the values of x
i
ξ, s
i
k
ξ, and y
i
m
ξ 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.


### 4.2.3. 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.(i)
Add the constraint setsx
i
=
x
^
i  (
∀
i
∈
I
), s
i
k
=
s
^
i
k  (
∀
i
∈
I
,
k
∈
K
), and y
i
m
=
y
^
i
m  (
∀
i
∈
I
,
m
∈
M
) to (P3), in which x
^
i, s
^
i
k, and y
^
i
m are the recorded values of x
i, s
i
k, and y
i
m in generating the UB in Step 3 of the subgradient algorithm. After adding these three constraint sets, (P3) merely contains scenario-specific decision variables. In this way, (P3) can further be decomposed into independent subproblems by scenario;
(ii)
Decompose (P3) into|
Ξ
| subproblems by scenario, and then solve each subproblem using a linear programming solver. Supposing the values of n
f
i
j
m
ξ and n
b
j
i
m
ξ are n
f
^
i
j
m
ξ and n
b
^
j
i
m
ξ, respectively, the objective value of (P3) equals ∑
i
∈
I
f
i
F
x
^
i
+
∑
i
∈
I
∑
k
∈
K
c
i
k
s
^
i
k
+
∑
i
∈
I
∑
m
∈
M
f
m
H
y
^
i
m
ξ
+
∑
ξ
∈
Ξ
∑
i
∈
I
∑
j
∈
J
∑
m
∈
M
p
ξ
c
m
H
t
-
i
j
m
(
n
f
^
i
j
m
ξ
+
n
b
^
j
i
m
ξ
).


## 4.2.1. Model Transformation

To decompose (P3) effectively, a new model named (P
3
′) is defined by duplicating x
i
ξ of x
i, s
i
k
ξ of s
i
k, and y
i
m
ξ of y
i
m for all scenarios ξ
∈
Ξ on the model (P3); then the following three constraint sets are added to (P
3
′):
(28)
(
1
-
p
1
)
x
i
1
=
∑
ξ
∈
Ξ
∖
{
1
}
p
ξ
x
i
ξ
,
∀
i
∈
I
,
(29)
s
i
k
1
-
s
i
k
ξ
=
0
,
∀
i
∈
I
,
k
∈
K
,
ξ
∈
Ξ
∖
{
1
}
,
(30)
y
i
m
1
-
y
i
m
ξ
=
0
,
∀
i
∈
I
,
m
∈
M
,
ξ
∈
Ξ
∖
{
1
}
,
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 i under all scenarios are identical, that is, x
i
1
=
x
i
2
=
⋯
=
x
i
|
Ξ
|, and (29) and (30) state similar limitations on s
i
k
ξ and y
i
m
ξ, respectively. From (28)–(30), it can be easily concluded that the two problems, (P3) and (P
3
′), are equivalent.(P
3
′) can be decomposed into independent subproblems by scenario if constraint sets (28)–(30) are relaxed. Therefore, the Lagrangian relaxation of (P
3
′) with respect to the constraint sets (28)–(30) can be described as follows:
(31)
f
2
LR
(
u
,
w
,
v
)
=
min
⁡
{
f
2
+
∑
ξ
∈
Ξ
∖
{
1
}
(
p
ξ
∑
i
∈
I
‍
u
i
(
x
i
1
-
x
i
ξ
)
‍
+
∑
i
∈
I
‍
∑
k
∈
K
‍
w
i
k
ξ
(
s
i
k
1
-
s
i
k
ξ
)
+
∑
i
∈
I
‍
∑
m
∈
M
‍
v
i
m
ξ
(
y
i
m
1
-
y
i
m
ξ
)
)
∑
ξ
∈
Ξ
∖
{
1
}
(
p
ξ
∑
i
∈
I
‍
u
i
(
x
i
1
-
x
i
ξ
)
‍
}
s
.
t
.
:
all
the
constraints
of
(
P
3
′
)
∖
constraints
(
28
)
–
(
30
)
,
where f
2 is the objective function of (P3), and the vectors u, w, and v 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:
(32)
z
LD
=
max
⁡
u
,
w
,
v
f
2
LR
(
u
,
w
,
v
)
.


## 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.Step 0 (initialization).
t
=
0 (iteration counter),
(33)
u
i
(
t
)
=
0
,
w
i
k
ξ
,
(
t
)
=
0
,
v
i
m
ξ
,
(
t
)
=
0
.
LetL
R
(
t
)
=
f
2
LR
(
u
(
t
)
,
w
(
t
)
,
v
(
t
)
) and U
B
(
t
)
= upper bound value of the objective of (P3). Initialize L
R
(
t
)
=
-
∞ and U
B
(
t
)
=
+
∞.Step 1 (solve each subproblem).
With givenu
(
t
), w
(
t
), and v
(
t
), decompose f
2
LR
(
u
,
w
,
v
) into subproblems by scenario, solve each subproblem using a linear programming solver, and add its objective value to L
R
(
t
+
1
).
Keepx
i
ξ
,
(
t
), s
i
k
ξ
,
(
t
), and y
i
m
ξ
,
(
t
)(the values of x
i
ξ, s
i
k
ξ, and y
i
m
ξ gotten in iteration t), respectively.
∀
i
∈
I, if (x
i
1
,
(
t
)
=
x
i
2
,
(
t
)
=
⋯
=
x
i
|
Ξ
|
,
(
t
)), then add constraint set x
i
ξ
=
x
i
1
,
(
t
)  (
∀
ξ
∈
Ξ
) to f
2
L
R
(
u
,
w
,
v
).
∀
i
∈
I and ∀
k
∈
K, if (s
i
k
1
,
(
t
)
=
s
i
k
2
,
(
t
)
=
⋯
=
s
i
k
|
Ξ
|
,
(
t
)), then add constraint set s
i
k
ξ
=
s
i
k
1
,
(
t
)  (
∀
ξ
∈
Ξ
) to f
2
LR
(
u
,
w
,
v
).
∀
i
∈
I and ∀
m
∈
M, if (y
i
m
1
,
(
t
)
=
y
i
m
2
,
(
t
)
=
⋯
=
y
i
m
|
Ξ
|
,
(
t
)), then add constraint set y
i
m
ξ
=
y
i
m
1
,
(
t
)  (
∀
ξ
∈
Ξ
) to f
2
LR
(
u
,
w
,
v
).Step 2 (compute subgradients).
(34)
subgrad
_
u
i
(
t
+
1
)
=
x
i
1
,
(
t
)
-
∑
ξ
∈
Ξ
‍
p
ξ
x
i
ξ
,
(
t
)
,
subgrad
_
w
i
k
ξ
,
(
t
+
1
)
=
s
i
k
1
,
(
t
)
-
s
i
k
ξ
,
(
t
)
,
subgrad
_
v
i
m
ξ
,
(
t
+
1
)
=
y
i
m
1
,
(
t
)
-
y
i
m
ξ
,
(
t
)
.Step 3 (update upper bound).
(35)
x
i
(
t
)
=
max
⁡
ξ
∈
Ξ
{
x
i
ξ
,
(
t
)
}
,
s
i
k
(
t
)
=
max
⁡
ξ
∈
Ξ
{
s
i
k
ξ
,
(
t
)
}
.
Compute the values ofy
i
m
(
t
)  (
∀
i
∈
I
,
m
∈
M
) following the adjustment process of Step 3 in Section 4.1 based on the values of y
i
m
ξ
,
(
t
)  (
∀
ξ
∈
Ξ
), and then solve (P3) where the scenario-unrelated variables, x
i, s
i
k, and y
i
m, are assigned to the values of x
i
(
t
), s
i
k
(
t
), and y
i
m
(
t
), respectively. Supposing the objective of (P3) is UB, let U
B
(
t
+
1
)
=
min
⁡
{
UB
,
U
B
(
t
)
}.
Record the values ofx
i, s
i
k, and y
i
m that have been utilized to generate the updated upper bound as x
^
i, s
^
i
k, and y
^
i
m, respectively.Step 4 (identify new step direction).
(36)
stepdirc
_
u
i
(
t
+
1
)
=
β
1
(
subgrad
_
u
i
(
t
+
1
)
)
+
β
2
(
subgrad
_
u
i
(
t
)
)
+
β
3
(
subgrad
_
u
i
(
t
-
1
)
)
,
stepdirc
_
w
i
k
ξ
,
(
t
+
1
)
=
β
1
(
subgrad
_
w
i
k
ξ
,
(
t
+
1
)
)
+
β
2
(
subgrad
_
w
i
k
ξ
,
(
t
)
)
+
β
3
(
subgrad
_
w
i
k
ξ
,
(
t
-
1
)
)
,
stepdirc
_
v
i
m
ξ
,
(
t
+
1
)
=
β
1
(
subgrad
_
v
i
m
ξ
,
(
t
+
1
)
)
+
β
2
(
subgrad
_
v
i
m
ξ
,
(
t
)
)
+
β
3
(
subgrad
_
v
i
m
ξ
,
(
t
-
1
)
)
,
where
β
1
=
0.6
,
β
2
=
0.3
,
β
3
=
0.1
.Step 5 (update multipliers).
(37)
sum
_
of
_
squares
=
∑
i
∈
I
(
subgrad
_
u
i
(
t
+
1
)
)
2
+
∑
i
∈
I
∑
k
∈
K
∑
ξ
∈
Ξ
(
subgrad
_
w
i
k
ξ
,
(
t
+
1
)
)
2
+
∑
i
∈
I
∑
m
∈
M
∑
ξ
∈
Ξ
(
subgrad
_
v
i
m
ξ
,
(
t
+
1
)
)
2
,
if (sum
_
of
_
squares
=
0)
goto Step6
else(38)
u
i
(
t
+
1
)
=
u
i
(
t
)
+
stepdirc
_
u
i
(
t
+
1
)
2.0
×
η
t
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
)
sum
_
of
_
squares
,
w
i
k
ξ
,
(
t
+
1
)
=
w
i
k
ξ
,
(
t
)
+
stepdirc
_
w
i
k
ξ
,
(
t
+
1
)
2.0
×
η
t
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
)
sum
_
of
_
squares
,
v
i
m
ξ
,
(
t
+
1
)
=
v
i
m
ξ
,
(
t
)
+
stepdirc
_
v
i
m
ξ
,
(
t
+
1
)
2.0
×
η
t
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
)
sum
_
of
_
squares
,
where η is a user-defined parameter and the value η
=
0.99 is used in the experiments.Step 6 (stopping criteria).
(39)
if
(
(
sum
_
of
_
squares
=
0
)
∥
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
U
B
(
t
+
1
)
≤
ε
)
∥
×
∥
(
U
B
(
t
+
1
)
-
L
R
(
t
+
1
)
U
B
(
t
+
1
)
≤
ε
)
∥
(
t
≥
N
2
)
)
stop;
else
t
=
t
+
1, goto Step 1
where the value ofN
2 is set to 200 and the value of ε is set to 0.001.
In Step0, 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 x
i, s
i
k, and y
i
m 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 (P
3
′) (which is equivalent to (P3)); that is, the values of x
i
ξ, s
i
k
ξ, and y
i
m
ξ 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.


## 4.2.3. 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.(i)
Add the constraint setsx
i
=
x
^
i  (
∀
i
∈
I
), s
i
k
=
s
^
i
k  (
∀
i
∈
I
,
k
∈
K
), and y
i
m
=
y
^
i
m  (
∀
i
∈
I
,
m
∈
M
) to (P3), in which x
^
i, s
^
i
k, and y
^
i
m are the recorded values of x
i, s
i
k, and y
i
m in generating the UB in Step 3 of the subgradient algorithm. After adding these three constraint sets, (P3) merely contains scenario-specific decision variables. In this way, (P3) can further be decomposed into independent subproblems by scenario;
(ii)
Decompose (P3) into|
Ξ
| subproblems by scenario, and then solve each subproblem using a linear programming solver. Supposing the values of n
f
i
j
m
ξ and n
b
j
i
m
ξ are n
f
^
i
j
m
ξ and n
b
^
j
i
m
ξ, respectively, the objective value of (P3) equals ∑
i
∈
I
f
i
F
x
^
i
+
∑
i
∈
I
∑
k
∈
K
c
i
k
s
^
i
k
+
∑
i
∈
I
∑
m
∈
M
f
m
H
y
^
i
m
ξ
+
∑
ξ
∈
Ξ
∑
i
∈
I
∑
j
∈
J
∑
m
∈
M
p
ξ
c
m
H
t
-
i
j
m
(
n
f
^
i
j
m
ξ
+
n
b
^
j
i
m
ξ
).


## 5. Case Study

### 5.1. 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 location-allocation 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 Figure1 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.Figure 1
Map of candidate TFs and AAs.8 types of demand, including water, food, medicine, doctors, nurses, rescuers, tents, and critical population, were considered in this case study. Table1 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.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
318The parameters of candidate TFs are given in Tables2 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 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
1054Table 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.60Table4 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).Table 4
Summary data of helicopter modes.
Parameter
Z-8G
Z-8S
S-70G
M-170G
M-170S
Number of helicopters available
20
15
16
15
12
Available hours
18
18
21
18
18
Capacity for commodities (tons)
4.0
0
3.5
4.5
0
Capacity for personnel (persons)
18
20
16
28
32
Capacity for critical population (persons)
8
15
6
9
16
Hiring cost (103 $/unit)
12
12
10
15
15
Operating cost (103 $/unit)
0.5
0.5
0.4
0.6
0.6The 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
∑
j
∈
J
h
-
k
d
j
k
/
∑
k
∈
K
R
∑
j
∈
J
h
-
k
d
j
k  (
∀
k
∈
K
R
) and 0.25
∑
j
∈
J
d
j
k
/
∑
k
∈
K
R
∑
j
∈
J
d
j
k  (
∀
k
∈
K
S
), 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.


### 5.2. 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
E
-
4, to 1.0
E
-
3 in the present case study in order to increase the convergence speed of computation.


#### 5.2.1. Explanation and Analysis of the Solution

Tables5 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 mobilization levels of most relief supplies did not reach their maximum capacities at any open TF, indicating that the total capacities of the 6 open TFs had not been completely utilized. 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.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
730Table 6
Initial deployment results of helicopter modes at TFs.
TFs
Z-8G
Z-8S
S-70G
M-170G
M-170S
1
11
—
10
9
—
2
—
—
—
—
—
3
—
—
—
—
—
4
—
—
—
—
—
5
—
—
—
—
—
6
3
—
—
1
11
7
6
15
6
5
1
8
—
—
—
—
—
Total
20
15
16
15
12The fill rates of all demand types at AAs under scenario 1 are shown in Table7. 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 (ω
k
k
′
C
=
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 ω
k
k
′
C, that is, 0.001.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.9Table8 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.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


#### 5.2.2. Effects of the Penalty Coefficientω
k
k
′
C 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ω
k
k
′
C
=
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 with the modification of ω
k
k
′
C, in which the correlation assumption among relief supplies was kept unchanged as aforementioned. The indices of the y-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 ∑
ξ
∈
Ξ
∑
j
∈
J
∑
k
∈
K
∑
k
′
∈
K
p
ξ
Δ
r
j
k
k
′
ξ
/
|
J
| and ∑
ξ
∈
Ξ
∑
k
∈
K
p
ξ
w
k
R
k
ξ, respectively. ω
k
k
′
C was modified by multiplying with 1
0
-
β, and the range of β is depicted by the x-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 ω
k
k
′
C is set as 0.001. As the value of penalty coefficient decreases further, the discrepancy among correlated relief supplies increases. When ω
k
k
′
C falls below 1
×
1
0
-
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 ω
k
k
′
C varies. These phenomena can provide valuable information on setting the value of this penalty coefficient of discrepancy.Total discrepancy and expected fill rate with variedw
k
k
′
C.
(a)
(b)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)).Expected value of total cost and the number of open TFs for different values ofα.
(a)
(b)


#### 5.2.3. Solution Quality and Evaluation of Computational Speed

To evaluate the performance of the proposed decomposition-based heuristic algorithm, the problem described in Section5.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 scenario-unrelated variable y
i
m  (
∀
i
∈
I
,
m
∈
M
) was assigned with an approximate value. However, y
i
m 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.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
943Table 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
12Table 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.1To 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 Figure4, in which the numerical value of the x-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.Performance comparing results between the heuristic algorithm and the CPLEX solver.
(a)
(b)


## 5.1. 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 location-allocation 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 Figure1 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.Figure 1
Map of candidate TFs and AAs.8 types of demand, including water, food, medicine, doctors, nurses, rescuers, tents, and critical population, were considered in this case study. Table1 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.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
318The parameters of candidate TFs are given in Tables2 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 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
1054Table 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.60Table4 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).Table 4
Summary data of helicopter modes.
Parameter
Z-8G
Z-8S
S-70G
M-170G
M-170S
Number of helicopters available
20
15
16
15
12
Available hours
18
18
21
18
18
Capacity for commodities (tons)
4.0
0
3.5
4.5
0
Capacity for personnel (persons)
18
20
16
28
32
Capacity for critical population (persons)
8
15
6
9
16
Hiring cost (103 $/unit)
12
12
10
15
15
Operating cost (103 $/unit)
0.5
0.5
0.4
0.6
0.6The 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
∑
j
∈
J
h
-
k
d
j
k
/
∑
k
∈
K
R
∑
j
∈
J
h
-
k
d
j
k  (
∀
k
∈
K
R
) and 0.25
∑
j
∈
J
d
j
k
/
∑
k
∈
K
R
∑
j
∈
J
d
j
k  (
∀
k
∈
K
S
), 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.


## 5.2. 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
E
-
4, to 1.0
E
-
3 in the present case study in order to increase the convergence speed of computation.


### 5.2.1. Explanation and Analysis of the Solution

Tables5 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 mobilization levels of most relief supplies did not reach their maximum capacities at any open TF, indicating that the total capacities of the 6 open TFs had not been completely utilized. 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.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
730Table 6
Initial deployment results of helicopter modes at TFs.
TFs
Z-8G
Z-8S
S-70G
M-170G
M-170S
1
11
—
10
9
—
2
—
—
—
—
—
3
—
—
—
—
—
4
—
—
—
—
—
5
—
—
—
—
—
6
3
—
—
1
11
7
6
15
6
5
1
8
—
—
—
—
—
Total
20
15
16
15
12The fill rates of all demand types at AAs under scenario 1 are shown in Table7. 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 (ω
k
k
′
C
=
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 ω
k
k
′
C, that is, 0.001.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.9Table8 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.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


### 5.2.2. Effects of the Penalty Coefficientω
k
k
′
C 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ω
k
k
′
C
=
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 with the modification of ω
k
k
′
C, in which the correlation assumption among relief supplies was kept unchanged as aforementioned. The indices of the y-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 ∑
ξ
∈
Ξ
∑
j
∈
J
∑
k
∈
K
∑
k
′
∈
K
p
ξ
Δ
r
j
k
k
′
ξ
/
|
J
| and ∑
ξ
∈
Ξ
∑
k
∈
K
p
ξ
w
k
R
k
ξ, respectively. ω
k
k
′
C was modified by multiplying with 1
0
-
β, and the range of β is depicted by the x-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 ω
k
k
′
C is set as 0.001. As the value of penalty coefficient decreases further, the discrepancy among correlated relief supplies increases. When ω
k
k
′
C falls below 1
×
1
0
-
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 ω
k
k
′
C varies. These phenomena can provide valuable information on setting the value of this penalty coefficient of discrepancy.Total discrepancy and expected fill rate with variedw
k
k
′
C.
(a)
(b)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)).Expected value of total cost and the number of open TFs for different values ofα.
(a)
(b)


### 5.2.3. Solution Quality and Evaluation of Computational Speed

To evaluate the performance of the proposed decomposition-based heuristic algorithm, the problem described in Section5.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 scenario-unrelated variable y
i
m  (
∀
i
∈
I
,
m
∈
M
) was assigned with an approximate value. However, y
i
m 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.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
943Table 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
12Table 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.1To 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 Figure4, in which the numerical value of the x-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.Performance comparing results between the heuristic algorithm and the CPLEX solver.
(a)
(b)


## 5.2.1. Explanation and Analysis of the Solution

Tables5 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 mobilization levels of most relief supplies did not reach their maximum capacities at any open TF, indicating that the total capacities of the 6 open TFs had not been completely utilized. 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.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
730Table 6
Initial deployment results of helicopter modes at TFs.
TFs
Z-8G
Z-8S
S-70G
M-170G
M-170S
1
11
—
10
9
—
2
—
—
—
—
—
3
—
—
—
—
—
4
—
—
—
—
—
5
—
—
—
—
—
6
3
—
—
1
11
7
6
15
6
5
1
8
—
—
—
—
—
Total
20
15
16
15
12The fill rates of all demand types at AAs under scenario 1 are shown in Table7. 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 (ω
k
k
′
C
=
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 ω
k
k
′
C, that is, 0.001.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.9Table8 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.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


## 5.2.2. Effects of the Penalty Coefficientω
k
k
′
C 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ω
k
k
′
C
=
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 with the modification of ω
k
k
′
C, in which the correlation assumption among relief supplies was kept unchanged as aforementioned. The indices of the y-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 ∑
ξ
∈
Ξ
∑
j
∈
J
∑
k
∈
K
∑
k
′
∈
K
p
ξ
Δ
r
j
k
k
′
ξ
/
|
J
| and ∑
ξ
∈
Ξ
∑
k
∈
K
p
ξ
w
k
R
k
ξ, respectively. ω
k
k
′
C was modified by multiplying with 1
0
-
β, and the range of β is depicted by the x-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 ω
k
k
′
C is set as 0.001. As the value of penalty coefficient decreases further, the discrepancy among correlated relief supplies increases. When ω
k
k
′
C falls below 1
×
1
0
-
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 ω
k
k
′
C varies. These phenomena can provide valuable information on setting the value of this penalty coefficient of discrepancy.Total discrepancy and expected fill rate with variedw
k
k
′
C.
(a)
(b)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)).Expected value of total cost and the number of open TFs for different values ofα.
(a)
(b)


## 5.2.3. Solution Quality and Evaluation of Computational Speed

To evaluate the performance of the proposed decomposition-based heuristic algorithm, the problem described in Section5.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 scenario-unrelated variable y
i
m  (
∀
i
∈
I
,
m
∈
M
) was assigned with an approximate value. However, y
i
m 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.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
943Table 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
12Table 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.1To 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 Figure4, in which the numerical value of the x-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.Performance comparing results between the heuristic algorithm and the CPLEX solver.
(a)
(b)


## 6. 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,ω
k
k
′
C, have obvious impact on the fill rate of correlated relief demands; when ω
k
k
′
C 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 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.



---

*Source: 939853-2014-01-16.xml*