Stochastic network design for disaster preparedness

Abstract

This article introduces a risk-averse stochastic modeling approach for a pre-disaster relief network design problem under uncertain demand and transportation capacities. The sizes and locations of the response facilities and the inventory levels of relief supplies at each facility are determined while guaranteeing a certain level of network reliability. A probabilistic constraint on the existence of a feasible flow is introduced to ensure that the demand for relief supplies across the network is satisfied with a specified high probability. Responsiveness is also accounted for by defining multiple regions in the network and introducing local probabilistic constraints on satisfying demand within each region. These local constraints ensure that each region is self-sufficient in terms of providing for its own needs with a large probability. In particular, the Gale–Hoffman inequalities are used to represent the conditions on the existence of a feasible network flow. The solution method rests on two pillars. A preprocessing algorithm is used to eliminate redundant Gale–Hoffman inequalities and then proposed models are formulated as computationally efficient mixed-integer linear programs by utilizing a method based on combinatorial patterns. Computational results for a case study and randomly generated problem instances demonstrate the effectiveness of the models and the solution method.

Keywords:

• Stochastic programming
• probabilistic constraint
• feasibility of second-stage
• feasible flow
• network reliability
• facility location
• disaster preparedness
• humanitarian logistics

Appendices

Appendix A: Elimination process

We briefly discuss the preprocessing method to eliminate the redundant Gale–Hoffman inequalities. The elimination procedure uses lower and upper bounds on the demand and capacity functions, which are assumed to be known. We denote the lower and upper bounds of the random capacity of arc (i, j) by v^l_ij and v^u_ij, respectively: P(v^l_ij ⩽ ξ^v_ij ⩽ v^u_ij) = 1, ∀(i, j) ∈ A. Recall that the random net demand at location i is obtained as ξ^d_i − r_i and M_l is the capacity of a facility of type l. Denoting the largest and the smallest possible values of the demand at node i by u^d_i and l^d_i, respectively, the random net demand is bounded from below by and bounded from above by u_i = u^d_i. The lower and upper bounds of the random total net demand associated with a subset B⊆I of nodes are

Example 1.

Consider the four-node network presented in Fig. A1

Fig. A1. Four-node network.

The lower and upper bounds of the random demand at each node and the random arc capacities are given in Table A1

Table A1. Upper and lower bounds of the random demand and arc capacities for the four-node network

.

Assuming ξ^v_ij = ξ^v_ji for all i < j and denoting the pre-positioned inventory levels by r_i, i = 1, …, 4, the corresponding stochastic Gale–Hoffman inequalities are given by (A1) (A2) (A3) (A4) (A5) (A6) (A7) (A8) (A9) (A10) (A11) (A12) (A13) (A14) (A15)

For this example, we can show that eight of 15 Gale–Hoffman inequalities are redundant.

The elimination procedure presented here consists of four stages; generating facets, elimination by upper bounds, elimination by lower bounds, and elimination by linear programming.

Stage 1. Generation of all the facets: We generate the facet inequalities using the recursive algorithm proposed by Wallace and Wets (1995). The method involves the introduction of a slack node with uncapacitated arcs connecting it to all other nodes. This allows us to formulate the so-called facet generating problem as a balanced supply and demand problem that is used to generate the facets of the feasible set of a network flow problem. In Stages 2, 3, and 4, we then withdraw the facets that become superfluous when we take into consideration the additional information about the lower and upper bounds of the demand and arc capacities. For the four-node network example, 12 inequalities ((A1(A2) )–(A6(A7) ), (A9(A10) )–(A13), and (A15(A16) )) are generated.

Stage 2. Elimination by upper bounds: (Prékopa and Boros, 1991): In contrast to Prékopa and Boros (1991), we model the arc capacities as random variables. We refer to the Gale–Hoffman inequality corresponding to the subset of nodes F⊆I as “inequality (F).” We eliminate an inequality (B) if . Indeed, since we have The three inequalities (A3(A4) ), (A5(A6) ), and (A10) are eliminated after executing Stage 2.

Stage 3. Elimination by lower bounds: If l(F) is finite, the inequality (F) can be rewritten as (A16) If G⊂F, then l(G) is also finite. Furthermore, if we have the inequality (A17) then (A16) implies that holds true. Indeed, the sum in d(F) − l(F) has non-negative terms and thus, The above result is the basis to eliminate the inequalities using the lower bounds on the net demand. Let H be the collection of the subsets of I, which have not been eliminated so far. For each F ∈ H, we identify all the subsets G⊆F for which (A17(A18) ) holds. Such subsets G ∈ H are eliminated and the set H representing the remaining inequalities are updated accordingly. The inequality (A15(A16) ) is eliminated after executing Stage 3.

Stage 4. Elimination by linear programming: Suppose that the set H is composed of the subsets of I which have not been eliminated so far. To check whether a subset F₀ ∈ H can be eliminated, we solve a linear programming problem. Let us introduce the auxiliary decision variables z_i, i ∈ I. Considering the random arc capacities, it is easy to see that the inequality (F₀) is a consequence of the other remaining inequalities (all F ∈ H, F ≠ F₀) only when the optimum value of the following linear problem is smaller than or equal to . In our example, the linear programming stage eliminates the inequality (A13(A14) ). Seven of the 15 inequalities remain after the elimination procedure.

Appendix B: Illustrative example for the combinatorial pattern-based reformulation method

Consider the following probabilistically constrained problem (see Lejeune, 2012a): (A18) where the random vector has 10 equally likely realizations represented by . The corresponding bivariate probability distribution is presented in Table A2.

Table A2. Probability distribution

The sufficient-equivalent set of the cut points is obtained as (A19) which includes three cut points associated with each component (n₁ = n₂ = 3 and n = 6). The central part of Table A3 displays the binarization of the probability distribution of with the set C^e. The right-hand side of Table A3 displays the truth table of the pdBf obtained with the set C^e for p = 0.7.

Table A3. Realizations, binary images, and truth table of pdBf

Table A4. Location-allocation decisions with model SP2_IP.

Table A5. Location-allocation decisions with model SP3_IP - p′ = 0.60.

Table A6. Location-allocation decisions with model SP3_IP - p′ = 0.70.

Table A7. Location-allocation policies obtained by model RTM for various values of the unit shortage Cost.

The MIP reformulation (see formulation SP2_IP) of the probabilistically constrained optimization problem (A18(A19) ) obtained with the combinatorial pattern reformulation method is as follows: (A20) The optimal solution is obtained as (u*, x*) = (0, 0, 1, 1, 0, 0, 1, 0) with the objective value of one.

Appendix C: Scenario generation

In this section, we present a scenario generation method that takes into account the dependency structures inherent in disaster relief networks. As in Rawls and Turnquist (2010), we consider five severity levels for hurricanes. We refer to levels 1 and 5 as the least and the most severe levels, respectively. We postulate that the magnitude of the commodity demand at a node and the damage inflicted to the capacity of an arc depend on the proximity to the epicenter of the disaster. Furthermore, we assume that the epicenter is at a coastal location.

To generate a scenario, the first step is to randomly select a coastal node as the epicenter and determine the severity level of the hurricane. In the second step, we define the nodes and arcs that are affected by the disaster and its impact on them. To that end, we determine three areas A₁, A₂, and A₃ that are defined as concentric areas with radius R₁, R₂, and R₃ centered at the epicenter A of the disaster (see Fig. A2). Next, an intensity level is assigned to each concentric area. The innermost area is most severely impacted by the disaster. Nodes within R₁ (resp., within R₂ but not within R₁; within R₃ but not within R₂) are impacted with the intensity assigned to A₁ (resp., A₂; A₃). Nodes that are not within distance R₃ of the epicenter are not affected. Considering the dependency structure based on the proximity to the epicenter, the intensity levels of three areas under five severity levels of hurricanes are specified. The demand at a node is a function of the population size associated with the node, the severity of the disaster, and a random coefficient based on the intensity level of the area to which the node belongs. This random coefficient is generated from a uniform distribution with lower and upper limits defined according to the corresponding intensity level. To determine how the disaster affects the capacity of an arc (i, j), we verify whether (i, j) crosses the areas A₁, A₂, and A₃. If the link (i, j) does not cross A₁, A₂, and A₃, then its capacity is not affected. Otherwise, we calculate the segment c¹_ij (resp. c²_ij, c³_ij) of link (i, j) that spans over A₁ (resp., A₂, A₃) and define . Thus, is the disaster-stricken area that is “most crossed” by link (i, j). We use the intensity level of to define the proportion by which the nominal capacity of link (i, j) is reduced. For example, for the capacity of the link BC (i.e., connecting B and C) in Fig. A2, the segment BD is longer than DC. Thus, the largest proportion of link BC is in area A₂ and BC is affected by the intensity level of A₂. We calculate the length of segments using the cosine formula for triangles.

Fig. A2. Illustration of the region for the scenario generation method.

Appendix D: Detailed numerical results

Additional information

Notes on contributors

Xing Hong

Xing Hong is a Senior Analyst of Continuous Improvement and Enterprise Optimization at United Airlines. His current focus is providing mathematical modeling and analysis services to support critical financial, operational, and strategic planning decisions. He received his Ph.D. in Systems Engineering from George Washington University.

Miguel A. Lejeune

Miguel Lejeune is a tenured Associate Professor of Decision Sciences at the George Washington University (GWU). He is the recipient of the CAREER/Young Investigator Research Grant from the Army Research Office and of the IBM Smarter Planet Faculty Innovation Award. Prior to joining GWU, he was a Visiting Assistant Professor in Operations Research at Carnegie Mellon University. He obtained his Ph.D. at Rutgers University. His research interests include stochastic programming, large-scale optimization, financial risk, and supply chain and disaster management.

Nilay Noyan

Nilay Noyan is an Associate Professor in the Manufacturing Systems and Industrial Engineering Program at Sabancı University, Turkey. She is a recipient of the Young Scientist (BAGEP) Research Award of Bilim Akademisi–the Science Academy, Turkey, and the Career Award of the Scientific and Technological Research Council of Turkey. She received her Ph.D. degree in Operations Research from Rutgers University in 2006. Her research interests include optimization, stochastic programming, risk modeling, and stochastic optimization applications with an emphasis on sustainable urban transportation, airline revenue management, and disaster relief network design.