The earth is facing challenges to work for the survival of human life during domino effect disasters. The emergency resource storage locations should be selected considering the probability of domino effect disasters. The first purpose of this study is to select the storage locations where domino effect probability is less. And second, facility development cost and transportation costs and costs for unutilized capacity have been optimized.
The work is a multiobjective optimization problem and solved with weighted sum approach. At first, the probabilities of domino effect due to natural disasters are calculated based on the earthquake zones. Then with that result along with other necessary data, the location to set up storage facilities and the quantities of resources that need to be transported has been determined.
The work targeted a country, Bangladesh for example. The authors have noticed that Bangladesh is currently storing relief items at warehouse which is under the domino effect prone region. The authors are proposing to avoid this location and identified the optimized cost for setting up the facilities. In this work, the authors pointed out which location has high probability of domino effect and after avoiding this location whether cost can be optimized, and the result demonstrated that this decision can be economical.
Disaster response authorities should try to take necessary proactive steps during cascading disasters. The novelty of this work is determining the locations to select storage facilities if the authors consider the probability of the domino effect. Then a facility location optimization model has been developed to minimize the costs. This paper can support policymakers to assess the strategies for selecting the location of emergency resource facilities.
1. Introduction
Natural and man-made disasters like earthquakes, tornados, chemical explosions, etc., can create a domino effect in nature. Earthquake waves cause collapses of infrastructures which may lead to death and injury of human beings, property damages, etc. Short circuit and gas explosion due to the movement of ground can initiate fire hazard creating postearthquake threat. These can be termed as domino effect phenomena. If an earthquake leads to the domino effect, it is important to determine emergency resource storage facility locations to provide support to the affected people.
Facility location is to select a place for setting up facilities to conduct a firms’ operation. Determining the geographic site of facilities is a critical decision that should be made strategically. Humanitarian logistics include logistics support (warehousing supplies and organizing delivery) before and after the occurrence of any disaster (Yáñez-Sandivari et al., 2021). The activities can be observed in phases: preparedness, mitigation, response and recovery (Yáñez-Sandivari et al., 2021). The main goal of humanitarian logistics is to minimize life and property losses where demand is highly variable and inventory management is so much challenging (İlhan, 2012). These losses can be minimized by responding to sudden incidents promptly. Different studies indicated that, for having a successful emergency distribution management and reducing damage to human beings, identification of actual demand of relief items for the affected people, arranging required items and prompt distribution of emergency resources or reliefs are some of the critical activities in case of disaster response management (Shareef et al., 2019). The effects of pre- and postdisaster relief funding on the relief system’s performance had been analyzed by integrating facility location and inventory decisions (Balcik and Beamon, 2008).
If an accident in one area spreads to nearby areas, this can be characterized as a domino effect or chain of accidents (Khakzad et al., 2013). Here first accident is considered as a primary event that triggers other accidents in nearby areas, and the overall consequence is more severe than the primary event. This domino effect can also be observed in natural disasters when one disaster can lead to another disaster. At the time of the domino effect, the damage of materials, equipment, industrial systems, environment, etc., may lead to injury of people which is the concept for escalation.
The contribution of this paper is to consider the domino effect in the determination of locations to establish emergency storage facilities. The objective of this paper is to answer the following research questions:
How to understand the location where the domino effect probability is highest?
Which locations can be selected to minimize the facility setup cost, transportation cost and unutilized capacity cost?
The first research question is important to answer because we should avoid highest probability domino effect location to establish storage facility. This is because, we need to establish resource facility in any stable zone which is less affected by disasters. To answer the second research question, we have developed a multiobjective optimization problem with weighted sum approach. We separated facility setup cost, transportation cost and unutilized capacity cost to give importance on them.
In this paper, Bangladesh has taken for numerical model development, and Bangladesh with the surroundings are divided into grid areas. Then different areas are divided into different earthquake risk zones. The purpose of this part is to select the high-probability domino effect–prone regions which should be omitted from the selection of establishing resource storage facility.
The result from this paper will help to determine the locations to set up emergency resources storage facilities to satisfy the demand of affected people when an earthquake leads to disastrous incidents due to the domino effect. Characteristics of different zones have been considered to obtain the probabilities of the domino effect. Cascading failure in nature has severe impact on infrastructures, and as some infrastructures are interdependent (Laugé et al., 2013), the probability of failure becomes intensified (Toscano et al., 2022). In this work, we considered the interdependency between electricity, gas, hilly regions on domino effect in nature because these interdependencies will intensify the probability of domino effect.
The sequences of methodology are – at first, we calculated probabilities of the domino effect for different zones. Using this information, we considered which locations should not be selected for storage facility. Then we selected locations from the remaining ones to establish storage facility so that costs can be optimized. The multiobjective optimization model will optimize costs. Here, we considered three types of costs – facility development cost, transportation costs of resources and the cost of unutilized capacity. In multiobjective optimization, one objective is to minimize facility development cost and transportation costs. Other objective is to minimize unutilized capacity cost. We have applied the weighted sum method to solve the multiobjective optimization. The optimization result has been compared with two solvers – Gurobi and PuLP – to validate the results.
2. Background
The background of this paper will focus on previous research on domino effect probability, interdependency of critical infrastructure on cascading disasters and emergency resource location selection. The inclusion of these articles will help us to understand the previous studies on our topic.
2.1 Domino effect probability
Domino effect can be considered as a chain or series of incidents where an incident propagates to other incidents. This may also be called cascading effect of disasters which is the chain of interconnected failures. Emergency response becomes difficult for this chain effect. There are some models that addressed how the probability of domino effect can be determined in some scenarios. Damage probability or domino effect frequency estimation has been focused on the study of domino effects (Khakzad et al., 2013). That damage probability is also known as escalation probability and can be estimated using different models.
The probit model developed by Cozzani and Salzano (2004) is a popular model for determining that escalation probability because of its simplicity and flexibility. The contribution of this work was to develop escalation probability. And this probability was then used to estimate the probability of the domino effect. Bayesian network and coupling effects had been considered by Chen et al. (2018) to analyze the domino effect of pool fire in the storage areas of the petrochemical industry. Zhou et al. (2021) applied improved probit model for an oil storage farm to address domino effect induced by fire. Ramírez-Camacho et al. (2015) developed a mathematical model for estimating the probability of domino effect in parallel pipelines. Landucci et al. (2009) applied a simplified approach to estimate the escalation thresholds and escalation probabilities triggered by fire scenarios. All these works focused on how to determine domino effect probability.
The contribution of the work of Qie and Rong (2022) was to make emergency decisions regarding cascading disasters. They pointed out important disaster scenario that can worsen the cascading effect of disasters. Triyanti et al. (2022) pointed out a challenge that political will as well as local transformation of governance are necessary to adopt to systematic risk governance for managing cascading risk in disasters. According to them, one of the main challenges are to create local level governance system. If this local level governance system can be established, it will help actors and institutions work simultaneously to reduce the risk to cascading disasters.
2.2 Interdependency of critical infrastructure
We considered in this work the interdependency of critical infrastructure on domino effect, and the critical infrastructure includes telecommunication, transportation, energy supply, etc. (Toscano et al., 2022). The interdependency of critical infrastructure can be geographic, cyber, physical or logical interdependencies, and the risk associated with cascading disaster can be minimized by protecting the critical infrastructures (Toscano et al., 2022). In this regard, Rehak et al. (2018) pointed out that the improvement of the concerned critical infrastructures is a managerial character, and security measures can contribute to improve the situation. The challenges to address in cascading disasters are inadequate safety system and inappropriate emergency plans that damage the infrastructure (Ricci et al., 2022). However, to help policymakers to prepare for cascading disasters, Khalili-Damghani et al. (2022) contributed to simulate a bi-objective optimization model for optimal distribution centers and minimized relief inventory cost by using geographic information system. They conducted the case study at the capital of Iran, Iranian Red Crescent Society in Tehran. They faced problem because of unavailability of official database and unavailability of volcanic eruption data.
2.3 Facility location problem focusing emergency response
As we consider earthquake as primary incident leading to domino effect, this earthquake may accelerate accidental scenarios by releasing hazardous materials. Earthquake can also trigger natural–technological risk (Campedel et al., 2008). There are many models that focused on mixed-integer linear programming model on facility location. The proposed model of Syah et al. (2022) used multiobjective nonlinear mixed-integer programming model, and this can lower the probability of nonoperation of all facilities at several level. Their aim was to sum up the waiting time and the service that affected people can avail. On the other hand, Yang et al. (2013) proposed a model that maximizes the effectiveness of allocated resources of emergency rescue and minimizes the cost of resources. Minimizing cost for the changes in facility operation states and transferring demand from facilities to customers which is stochastic were represented with a model (Devin et al., 2015). Dynamic demand capacitated location problem for multiperiod, where facilities are allowed to be relocated in each period and are kept at a fixed location, had been considered by Torres-Soto and Üster (2010). Garrido and Aguirre (2017) proposed a mixed-integer model, which was a stochastic one, that reduced response time during the assistance of victims of the earthquake. Lu and Sun (2020) developed a multiobjective model that minimized penalty cost and the sum of allocation cost for allocating emergency resources in the rescue process of metro emergency, and they used a scenario-based analysis.
Su et al. (2016) proposed a constrained integer linear programming model and identified overlapping disaster response coalitions for concurrent incidents. The authors worked on the problem to allocate multiple emergency resources for multiple concurrent incidents in a parallel manner. They considered both response time and emergency resource cost for this model. To have the optimal allocation of available resources, an integrated simulation–optimization tool was proposed (Alsubaie et al., 2015). Their result will help to maximize the operational capacity of a hospital which is a critical infrastructure.
Some authors gave their research attention on selecting response center for disaster. Ghosh and Gosavi (2017) developed a semi-Markov model that identified how a system degrades due to an earthquake and how the responding center can be selected optimally. This model can be used in smart cities that will need the system to be responded automatically. However, Wang (2021) proposed a model from a cross-regional point of view by considering multiperiod emergency resource allocation. Determining the lowest allocation cost, shortest delivery time and allocating resources at maximum coverage rate were the targets of the model. Shavarani (2019) focused on developing model on minimizing response time after disaster and maximizing survival rate of the people by providing service. The author solved multilevel facility location problem using hybrid genetic algorithm. The novelty of this paper is it used drone to provide service.
From the above studies, most of the cascading disaster-related articles are trying to find out how the disaster can happen and how to mitigate them. The interdependency among the critical infrastructures is mostly studying the factors related to interdependency and how the interdependency will not intensify cascading failures. The state of the art of our paper is to integrate domino (cascading) effect probability with optimal emergency resource location selection by minimizing fixed cost, transportation cost and unutilized capacity cost.
3. Proposed model
3.1 Problem statement and assumptions
The work in this paper is related to humanitarian logistics when the domino effect is under consideration. For the domino effect, an earthquake is considered as a primary incident that may lead to secondary incidents – building collapse and fire. So, the probability of primary incident and probability of domino effect can be considered to optimize the storage of resources in facilities. This work is focused to select locations for setting up storage facilities considering the probability of domino effect and probability of a primary incident. At the same time, fixed cost, transportation cost and unutilized capacity cost are minimized using multiobjective optimization. Demand from population census is dependent on domino effect probability. Some assumptions and limitations are:
Demand is not dependent on earthquake magnitude.
The capacity of storage facilities is predefined.
To determine the domino effect probability, escalation probability is needed to calculate. Here, “probit function” has been used for calculating escalation probability. As there is no predefined value for probit coefficients, trial and error method has been applied to determine the value of the coefficients. The Appendix indicates the results.
Most of the input values has been assumed logically as there are no predefined values for them.
3.2 Model development
The step-by-step methodology of this model is as follows:
We have integrated probability of domino effect in the resource storage location selection. For this, we calculated probabilities of the domino effect based on different earthquake zones. As there are no preestablished values, we have assumed some values based on the risk zones of earthquake. To calculate the probability of domino effect, we need probability of primary event and probability of escalation. Here, the probability of primary event is the probability of earthquake at different risk levels. This has been written detailed in Section 3.2.1.
From the above information, we avoided high-probability domino effect zone for storage facility. We also focused on cost optimization and thus selected locations from the remaining ones to establish storage facility. We have emphasized on minimizing fixed cost, transportation cost and unutilized capacity cost. We applied the multiobjective optimization model to optimize costs. We have applied weighted sum approach to solve this. Section 3.2.2 describes this.
Section 4 contains a numerical example, and at last the optimization result has been compared with two solvers – Gurobi and PuLP.
3.2.1 Calculation of probability of domino effect
The probability of the domino effect can be calculated by multiplying the probability of the primary event and the escalation probability (Khakzad et al., 2013).
That is, Probability of domino effect = Probability of primary event × Escalation probability.
Earthquake is considered as the primary event, and the probability of primary event is given a symbol p(primary), whereas the escalation probability is given a symbol pE.
The escalation probability can be calculated using the probit model. The purpose of this probit model is to estimate probability where an observation would fall into a specific category. Probit analysis is a popular method for evaluating “dose–effect” relation in case of response to some unfavorable condition (Finney, 1971).
From Cozzani et al. (2005), we are defining escalation probability as the following:
where pE is the escalation probability (0 ≤ pE ≤ 1), and this is “cumulative expression for a normal Gaussian probability distribution function.” Here, m is defined as follows:
where factors is the independent variable, which is pointed out as “dose”; σ and μ are the variance and median of the Gaussian distribution. We have taken σ = 1 and μ = 0 for normal distribution. The Y variable in equation (1) is called “probit function” and is used to relate damage to the factors responsible for the domino effect after the primary event. The Y is calculated as below (Cozzani, et al., 2005):
Here, Y is the probit for damage, and K1 and K2 are probit coefficients.
We have set the factors as a numerical value which includes two things – number of buildings (per square kilometer) and absence or presence of some interdependent critical infrastructures. The selected infrastructures are gas line, industrial area, hilly area and coastal area. The presence of all these can intensify the probability of domino (cascading) effects. So, the following equation stands:
where N defines the number of buildings per square kilometer, and (Z × K) denotes a number that would include the impact of the absence or presence of gas line network, industries, a hilly area and/or a coastal area in the determination of Y. The N value will indicate whether the area is densely populated or not. (Z × K) includes the effect of these characteristics in the determination of domino probabilities. Z is a variable that depends on the characteristics of the areas. K is a constant and is assumed here as a factor to include the impact of the presence or absence of the abovementioned characteristics of the areas. For determining the values of Z and K, four levels have been assigned based on the characteristics of the areas. Level 0 means there is no interdependent infrastructures, whereas Level 4 means there is highest interdependency among infrastructures and highest probability of domino effects. The meanings of four levels are mentioned below:
Level 0: if there is no gas line network and no industry in the demand area and the area is not a hilly or a coastal area;
Level 1: if the demand area has a gas line network or some industries or the area is either a hilly or a coastal area;
Level 2: if the demand area is hilly and there is a gas line network; if the area has some industries and a gas line network; and if the area has large coastal areas;
Level 3: if the area has a high density of industries and has dense gas line networks; if the area has a high density of industries and there are large hilly areas; and if the area has a large coastal area and has some industries or a gas line network.
We did not find any predefined value for probit coefficients for this type of case. So, we used trial and error for determining probit coefficients K1 and K2 to obtain feasible value for escalation probability. Trial and error method is substituting different values of the variables for repeated attempts until the desired result is obtained. We randomly selected some values and see the results for escalation probability. In the Appendix, Table A1 and Table A2 show some trial and error. We have also done sensitivity analysis by analyzing the effect of some parameters while holding other parameters constant. We kept N, (Z × K) constant in Table A1, and in Table A2, the N, (Z × K) values vary and thus we observed the effect on values for probit coefficients.
Trial and error to determine the probit coefficients (keeping N, Z × K constant)
| Probability of primary event, pp | Escalation probability | |||||
|---|---|---|---|---|---|---|
| N (number of buildings per square kilometer) | Z × K | Probit function, Y | Escalation probability, pE | Probability of domino effect, pp × pE | ||
| 0.7 | 50 | 5 × 100 = 500 | K1 = −1, K2 = 5 | Y = −1 + 5(ln(50 + 500)) = 30.55 | 1 | 0.7 |
| 0.7 | 50 | 5 × 100 = 500 | K1 = −5, K2 = 1.5 | Y = −5 + 1.5(ln(50 + 500)) = 4.46 | 0.29 | 0.203 |
| 0.7 | 50 | 5 × 100 = 500 | K1 = −5, K2 = 0.5 | Y = −5 + 0.5(ln(50 + 500)) = −1.85 | 3.69 × 10−12 | 2.583 × 10−12 |
| 0.7 | 50 | 5 × 100 = 500 | K1 = −6, K2 = 0.5 | Y = −6 + 0.5(ln(50 + 500)) = −2.85 | 2.08 × 10−15 | 1.456 × 10−15 |
| 0.7 | 50 | 5 × 100 = 500 | K1 = 5, K2 = −1.5 | Y = 5 − 1.5(ln(50 + 500)) = −4.46 | 1.54 × 10−21 | 1.078 × 10−21 |
Trial and error to determine the probit coefficients (keeping K1, K2 as constant and by varying N, Z × K)
| Probability of primary event, pp | Escalation probability | |||||
|---|---|---|---|---|---|---|
| N | Z × K | Probit function, Y | Escalation probability, pE | Probability of domino effect, pp × pE | ||
| 0.7 | 100 | 5 × 100 = 500 | K1 = −1, K2 = 5 | Y = 30.98 | 1 | 0.7 |
| 6 × 100 = 600 | Y = 31.78 | 1 | 0.7 | |||
| 7 × 100 = 700 | Y = 32.42 | 1 | 0.7 | |||
| 950 | 5 × 100 = 500 | Y = 35.40 | 1 | 0.7 | ||
| 6 × 100 = 600 | Y = 35.73 | 1 | 0.7 | |||
| 7 × 100 = 700 | Y = 36.04 | 1 | 0.7 | |||
| 0.7 | 100 | 5 × 100 = 500 | K1 = −5, K2 = 1.5 | Y = 4.60 | 0.34 | 0.24 |
| 6 × 100 = 600 | Y = 4.83 | 0.43 | 0.30 | |||
| 7 × 100 = 700 | Y = 5.03 | 0.51 | 0.36 | |||
| 950 | 5 × 100 = 500 | Y = 5.92 | 0.82 | 0.57 | ||
| 6 × 100 = 600 | Y = 6.02 | 0.85 | 0.595 | |||
| 7 × 100 = 700 | Y = 6.11 | 0.87 | 0.61 | |||
| 0.7 | 100 | 5 × 100 = 500 | K1 = −5, K2 = 1 | Y = 1.40 | 0.00016 | 1.12 × 10−4 |
| 6 × 100 = 600 | Y = 1.55 | 0.00028 | 1.96 × 10−4 | |||
| 7 × 100 = 700 | Y = 1.68 | 0.00045 | 3.15 × 10−4 | |||
| 950 | 5 × 100 = 500 | Y = 2.28 | 0.0033 | 2.31 × 10−3 | ||
| 6 × 100 = 600 | Y = 2.34 | 0.0039 | 2.73 × 10−3 | |||
| 7 × 100 = 700 | Y = 2.41 | 0.0048 | 3.36 × 10−3 | |||
| 0.7 | 100 | 5 × 100 = 500 | K1 = −6, K2 = 0.5 | Y = −2.8 | 3.10 × 10−15 | 2.17 × 10−15 |
| 6 × 100 = 600 | Y = −2.72 | 5.82 × 10−15 | 4.07 × 10−15 | |||
| 7 × 100 = 700 | Y = −2.66 | 9.30 × 10−15 | 6.51 × 10−15 | |||
| 950 | 5 × 100 = 500 | Y = −2.36 | 9.20 × 10−14 | 6.44 × 10−14 | ||
| 6 × 100 = 600 | Y = −2.33 | 1.15 × 10−13 | 8.05 × 10−14 | |||
| 7 × 100 = 700 | Y = −2.30 | 1.44 × 10−13 | 1.008 × 10−13 | |||
It is observed from the trial and error that, the escalation probability pE remains 1 for all the following values of K1 = −1 and K2 = 5; K1 = −2 and K2 = 4; K1 = −3.5 and K2 = 4.5.
If K1 = −5 and K2 = 1.5, the escalation probability is 0.29, and the probability of domino effect is 0.203 which is feasible. Again, for K1 = −5 and K2 = 0.5, the probabilities become very small.
So, from the above calculation, we have taken probit coefficients as K1 = −5 and K2 = 1.5. Equation (4) becomes:
3.2.2 Objective function and constraints
This work is a mixed-integer linear optimization problem. This is a multiobjective optimization problem as we focused on minimizing fixed cost, transportation cost and unutilized capacity cost. The weighted sum approach has been applied to solve the multiobjective optimization.
Weightage w1 has been assigned for fixed cost and transportation cost, and weightage w2 = (1 − w1) has been assigned for unutilized capacity cost.
To formulate the model, we have used the article guidelines from Brown and Dell (2007).
Index use [∼Cardinality]:
i ∈ I candidate locations options [∼n]; and
j ∈ J demand area options [∼m].
Given data [Units]:
fixed_costi fixed cost for setting up storage facility at location i [Tk];
distancei,j distance from storage facility location i to demand area j [km];
unit_costi,j per unit distance cost for transporting per unit resource from location i to j [Tk/km/kg];
dj demand at demand area j [kg];
pj(primary) probability of primary incident at demand area j;
pj(domino) probability of domino effect at demand area j;
max_pj(domino) maximum probability of domino effect among all demand areas;
storage_costi storage cost for per unit unutilized capacity [Tk/kg]; and
capacityi capacity of the storage facility at location i [kg].
Decision variables [Units]:
Xi which location will be chosen to set up a storage facility [binary]; and
Fi,j total amount of resources transported from a storage facility i to demand area j [kg].
Formulation:
The objective function and constraints of the problem are as follows:
s.t.
The objective function, equation (6), minimizes fixed cost for facility development, transportation costs for transporting resources and unutilized capacity cost. This consists of two parts. In the first part, fixed cost and transportation costs have been considered. In the second part, a cost associated with unutilized capacity has been considered. Weightage w1 has been given to the first part and weightage (1 − w1) = w2 has been given to the second part. As it is a two-objective optimization problem, so while one weightage is w1, another one can be (1 − w1). The final objective is to minimize the combined costs.
Equation (7) describes demand for a demand area will be fulfilled from different storage facilities by transporting emergency resources. The transported amount will be equal to the summation of the probability of primary incident and probability of domino effect multiplied by demand of that area. So, demand is dependent on these probabilities.
Equation (8) means transported resource quantities from a storage facility will be less than or equal to the capacity of that facility. And if any facility will not be chosen, then the transported amount from that facility will be zero.
Decision variable domains are defined by equations (9)–(12). Equation (12) states that no storage facility would be set up in the location where the probability of domino effect is maximum if the location is one of the candidate locations. As high domino effect probability can also affect facilities, so, it would be better to skip that location.
4. A numerical example
A numerical example is presented here to have a light on the proposed model. For simplicity, a single resource problem has been considered.
We have chosen our country Bangladesh here for solving the model. Bangladesh and the surrounding are divided into nine grid areas. Grid center coordinates are taken as storage facility location coordinates. Both the storage facility location coordinate and the demand area location coordinate are considered same. Figure 1 shows the candidate facility locations which are denoted as 1–9.
The box area in Figure 1 indicates the grid area, and the number in the center of the grid (1–9) is the candidate facility locations and the demand area locations. The coordinates are same in the presented problem. But the developed model will also be applicable if these are not the same. We use these coordinates for obtaining distance of the locations which is distancei,j of Section 3.2.2. From this we will get transportation cost. We used Google Earth Pro Software for determining coordinates and QGIS software (version 3.4.5) for determining distance. Table 1 shows the coordinates of the candidate facility locations and demand areas indicated in Figure 1.
Coordinates for the candidate storage facility locations and demand areas
| Location, i/Demand area, j | Latitude | Longitude |
|---|---|---|
| 1 | 25.989 | 89.005 |
| 2 | 25.989 | 91.003 |
| 3 | 25.989 | 92.987 |
| 4 | 24.001 | 89.005 |
| 5 | 24.001 | 91.003 |
| 6 | 24.001 | 92.987 |
| 7 | 21.986 | 89.005 |
| 8 | 21.986 | 91.003 |
| 9 | 21.986 | 92.987 |
Table 2 indicates the distance between candidate facility location and demand areas in kilometers.
Distances (km) between storage facility location (i) and demand areas (j)
| i | j | Distance | i | j | Distance | i | j | Distance | i | j | Distance | i | j | Distance | i | j | Distance | i | j | Distance | i | j | Distance | i | j | Distance |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 1 | 0 | 2 | 1 | 201 | 3 | 1 | 399 | 4 | 1 | 221 | 5 | 1 | 299 | 6 | 1 | 459 | 7 | 1 | 444 | 8 | 1 | 488 | 9 | 1 | 601 |
| 1 | 2 | 201 | 2 | 2 | 0 | 3 | 2 | 199 | 4 | 2 | 299 | 5 | 2 | 221 | 6 | 2 | 298 | 7 | 2 | 488 | 8 | 2 | 444 | 9 | 2 | 488 |
| 1 | 3 | 399 | 2 | 3 | 199 | 3 | 3 | 0 | 4 | 3 | 459 | 5 | 3 | 298 | 6 | 3 | 221 | 7 | 3 | 601 | 8 | 3 | 488 | 9 | 3 | 444 |
| 1 | 4 | 221 | 2 | 4 | 299 | 3 | 4 | 459 | 4 | 4 | 0 | 5 | 4 | 204 | 6 | 4 | 406 | 7 | 4 | 224 | 8 | 4 | 303 | 9 | 4 | 466 |
| 1 | 5 | 299 | 2 | 5 | 221 | 3 | 5 | 298 | 4 | 5 | 204 | 5 | 5 | 0 | 6 | 5 | 202 | 7 | 5 | 303 | 8 | 5 | 224 | 9 | 5 | 302 |
| 1 | 6 | 459 | 2 | 6 | 298 | 3 | 6 | 221 | 4 | 6 | 406 | 5 | 6 | 202 | 6 | 6 | 0 | 7 | 6 | 466 | 8 | 6 | 302 | 9 | 6 | 224 |
| 1 | 7 | 444 | 2 | 7 | 488 | 3 | 7 | 601 | 4 | 7 | 224 | 5 | 7 | 303 | 6 | 7 | 466 | 7 | 7 | 0 | 8 | 7 | 207 | 9 | 7 | 412 |
| 1 | 8 | 488 | 2 | 8 | 444 | 3 | 8 | 488 | 4 | 8 | 303 | 5 | 8 | 224 | 6 | 8 | 302 | 7 | 8 | 207 | 8 | 8 | 0 | 9 | 8 | 205 |
| 1 | 9 | 601 | 2 | 9 | 488 | 3 | 9 | 444 | 4 | 9 | 466 | 5 | 9 | 302 | 6 | 9 | 224 | 7 | 9 | 412 | 8 | 9 | 205 | 9 | 9 | 0 |
Per unit transportation cost has been assumed as 250 Tk for per unit (kg) relief item to travel per kilometer distance. So, unit_costi,j = 250 (Tk/km/kg) for transporting resources from location i to j. Our model will work for any logical unit_costi,j value. Emergency resources will be assumed as necessary items given in a box in the unit of kg. As unit–distance transportation cost is multiplied with distance, where the distance is zero for the same facility location and demand area, the total transportation cost will be zero, which is not actually logical in real scenario. So, we have taken 50 km distance to travel at the same coordinate for calculating transportation cost at these locations.
The fixed cost of setting up and operating a facility at each candidate location is assumed to be different for different locations, as in some locations, it is cheaper to establish storage facilities, whereas some locations may experience high cost to set up facilities. Table 3 shows the data.
Fixed cost of each candidate location
| Location, i | Fixed cost (×105) (Tk) |
|---|---|
| 1 | 800 |
| 2 | 700 |
| 3 | 600 |
| 4 | 500 |
| 5 | 400 |
| 6 | 500 |
| 7 | 700 |
| 8 | 900 |
| 9 | 800 |
For calculation convenience, values are converted to the unit of 105. Location 8 shows the highest fixed cost for setting up facility because it is situated near the coastal area where extra arrangements are needed for safety.
Storage cost for per unit unutilized capacity is taken as 2 (Tk/kg) for all locations.
Capacity has been taken using random number generation. Table 4 describes this.
Random capacity ranges for the storage facilities
| Capacity (×105) (kg) | Random number between 200 and 400 |
| Random number between 300 and 600 |
Now we have to estimate the probability of primary event (p(primary)) and the probability of domino effect (p(domino)). We considered earthquake as a primary incident that may lead to cascading disasters like building collapse or fire. As there is no established value for the probability of occurring earthquake, we have chosen the probability according to the earthquake risk zones.
Earthquake-prone regions are expressed by seismic zone. Frequent earthquake-faced regions have high seismicity than less-frequent earthquake regions. Ali and Choudhury divided Bangladesh into three seismic zones – Zone I, Zone II and Zone III (Earthquake – Banglapedia [WWW Document], 2021). The most active zone for earthquakes is Zone I which is the high-risk zone, Zone II is a moderate risk zone and Zone III is a seismically quiet zone. This seismic zonation has been done based on the behavior of tectonic blocks of Bangladesh rather than the ground conditions. From the grid areas of Figure 1, we selected different earthquake risk (high, medium and low) zones. Areas 1, 2, 3, 6 and 9 are assumed to be in a high-risk zone; Areas 4 and 5 are in a medium-risk zone; and Areas 7 and 8 are in a low-risk zone.
As there are no previous data or any specific data for calculation of domino effect, we have assumed these data based on the earthquake risk zones. To obtain probability of domino effect, we need primary event probability and escalation probability. To have escalation probability, pE, we must determine “probit function.” We have to determine values for equations (1) and (5). We first obtain values for equation (5), and we assigned levels (mentioned in Section 3.2.1) to different areas.
The values of Z and K (symbol meaning of Z and K mentioned in Section 3.2.1) have been taken in a way so that the escalation probability can have a feasible value where the probability value remains between 0 and 1, and it varies with the change of levels. So, values for Z are taken randomly for different levels:
Z = 0 for Level 0;
Z = 5 for Level 1;
Z = 6 for Level 2; and
Z = 7 for Level 3.
The constant K is assumed to have a value of 100. We have also tried with other values, but other values did not give feasible result for probability, so we have chosen this value for K.
Table 5 shows the values taken here for equation (5). Table 5 is showing the characteristics of different selected areas and assigning a level for determining (Z × K) values. With the help of this value, we will get escalation probability and probability of domino effect.
Assigning level and value of (Z × K) and the probability of domino effect
| Location, i/Demand area, j | Characteristics of the area | Level | (Z × K) | Probability of primary event (pj(primary)) | N (number of buildings per square kilometer) | Probit function (Y = −5 + 1.5 ln(N + (Z × K))) | Escalation probability | Probability of domino effect (pj(domino) = pj(primary) × pE) |
|---|---|---|---|---|---|---|---|---|
| Area 1 | No gas line network, no industrial zone, not a hilly or a coastal area | Level 0 | 0 × 100 = 000 | 0.7 | 100 | 1.91 | 0.001 | 0.0007 |
| Area 2 | Presence of gas line network | Level 1 | 5 × 100 = 500 | 0.7 | 50 | 4.46 | 0.29 | 0.203 |
| Area 3 | Presence of hilly area and gas line network | Level 2 | 6 × 100 = 600 | 0.7 | 70 | 4.76 | 0.41 | 0.287 |
| Area 4 | Presence of some industries and gas line network | Level 2 | 6 × 100 = 600 | 0.6 | 500 | 5.50 | 0.69 | 0.414 |
| Area 5 | Presence of high density of industries and dense gas line network | Level 3 | 7 × 100 = 700 | 0.6 | 950 | 6.11 | 0.87 | 0.522 |
| Area 6 | Presence of hilly area and gas line network | Level 2 | 6 × 100 = 600 | 0.7 | 60 | 4.74 | 0.40 | 0.28 |
| Area 7 | Presence of coastal area | Level 1 | 5 × 100 = 500 | 0.5 | 50 | 4.46 | 0.29 | 0.145 |
| Area 8 | Presence of large coastal area | Level 2 | 6 × 100 = 600 | 0.5 | 60 | 4.74 | 0.40 | 0.2 |
| Area 9 | Presence of hilly area | Level 1 | 5 × 100 = 500 | 0.7 | 70 | 4.52 | 0.32 | 0.224 |
For escalation probability, we also need N which is the number of buildings per square kilometer. We have chosen high N value for densely populated areas and low N values for low densely populated areas.
The probability of occurring primary event (p(primary)) has been assigned three values (70%, 60% and 50%) according to the risk of earthquake of different zones. As there is no specific data for the probability of earthquake, we have chosen these values randomly based on the intensity of risk of occurring earthquake; this is how p(primary) has been defined. The place where there is high risk of occurring earthquake, it has given highest percentage (randomly 70%) and the area for low risk of earthquake given lowest percentage (randomly 50%). Probability of primary event 70% means that the probability of occurring earthquake at that zone is assumed (as there is no predefined value) to be 70%. Table 5 also shows the probability of primary event and the probability of domino effect.
From the table we can see that maximum probability of domino effect (max_p(domino)) among all demand areas is Area 5, which is a densely populated area, and many industries and dense gas line network are also present here.
We need demand to understand the real-life scenario; so, Table 6 describes the estimated demand.
Population and estimated demands of the demand areas
| Demand area, j | Districts | Population | Demand, dj (×105) (kg) | Final estimated demandj (×105) (kg) (values rounded up) |
|---|---|---|---|---|
| Area 1 | Panchagarh, Thakurgaon, Lalmonirhat, Kurigram, Nilphamari, Dinajpur, Rangpur, Joypurhat, Gaibandha | 16,701,526 | 167.02 | 118 |
| Area 2 | Sherpur, Sunamganj | 3,826,293 | 38.26 | 35 |
| Area 3 | Sylhet | 3,434,188 | 34.34 | 34 |
| Area 4 | Naogaon, Chapainawabganj, Rajshahi, Natore, Sirajganj, Bogra, Pabna, Jamalpur, Tangail, Kustia, Meherpur, Chuadanga, Rajbari, Jhenaidah, Faridpur, Magura, Narail, Jessore, Gopalganj, Manikganj | 38,904,059 | 389.04 | 395 |
| Area 5 | Netrokona, Mymensingh, Gazipur, Narsingdi, Munshiganj, Kishoreganj, Habiganj, Dhaka, Narayanganj, Brahamanbaria, Cumilla, Shariatpur, Chandpur, Khagrachhari, Feni, Madaripur | 49,424,400 | 494.24 | 555 |
| Area 6 | Maulvibazar | 1,919,062 | 19.19 | 19 |
| Area 7 | Satkhira, Khulna, Pirojpur, Patuakhali, Barguna, Bagerhat | 9,322,468 | 93.22 | 61 |
| Area 8 | Barisal, Bhola, Jhalkathi, Laxmipur, Noakhali, Chattogram, Cox’s Bazar | 19,527,387 | 195.27 | 137 |
| Area 9 | Rangamati, Bandarban | 984,314 | 9.84 | 10 |
For estimating demand, dj, we collected information of population of Bangladesh from the population census of 2011 (Bangladesh: Districts and Cities – Population Statistics, Maps, Charts, Weather and Web Information [WWW Document], 2023). Demand dj is taken as 1 kg light weighted emergency items for one person for one month. The final demand is dependent on the probability of occurrence of the primary incident and the probability of domino effect.
So, final estimated demand for any demand area j will be:
These data have been provided in Table 7.
Gurobi and PuLP solver results for (w1 = 0.5 and (1 − w1) = 0.5)
| Weightage | Capacity range (×105) (kg) | Solver name | Objective function (×105) (Tk) | Selected facility |
|---|---|---|---|---|
| w1 = 0.5 and (1 − w1) = 0.5 | Random between 300 and 600 | Gurobi | 1,135.73 | Area 2, Area 4, Area 6 |
| PuLP | 1,147.73 | Area 2, Area 4, Area 6 | ||
| Random between 200 and 400 | Gurobi | 1,478.27 | Area 2, Area 3, Area 4, Area 6 | |
| PuLP | 1,564.26 | Area 2, Area 3, Area 4, Area 6 |
5. Results and discussion
Currently, Bangladesh purchase the relief items (non-food, food and cash) centrally and then store in three central warehouses (here Areas 5, 7 and 8 in Figure 1) and then send to the affected areas (Shareef et al., 2019). Though till now no huge domino effect disasters happened here, but this country is at the high risk of cascading disasters due to the changed climatic conditions and interdependent infrastructures. So, we need to focus on this sector. Area 5, which is currently using as a warehouse, is at the highest risk of domino effect according to the calculation. We are proposing to avoid this location as domino effect can also affect the storage facility. It would be better if we could evacuate the population from the domino effect affected area but in a short time this could not be possible. So, we are focusing on delivering relief items to the affected people from the storage facilities. And for this, we should establish storage facility in a location that is not affected by domino effect.
The numerical problem is solved using the Python programming language. The integrated development environment used for this is PyCharm (version 2020.3.1). To verify the results, Gurobi and PuLP solvers have been used. PuLP is an open-source solver, whereas Gurobi is a commercial solver which is faster than PuLP. PuLP has default solver CBC (COIN-OR Branch and Cut), and Gurobi has different solvers for different models which has no limit for the decision variable numbers. We obtained results for different combinations of weightages (giving weights for two objective functions) and generated Pareto front. If we change the weights, the results also change. We analyzed two different capacity ranges; random between 300 (kg) to 600 (kg) and 200 (kg) to 400 (kg) (all are in the unit of 105).
The shape of Pareto front is concave here. This means we cannot both optimize fixed cost and transportation cost with unutilized capacity cost. There is a trade-off here.
Figure 2(a) and (b) shows the Pareto front for capacity range 300 (kg) to 600 (kg) (×105) for two solvers. The solvers give almost the same results. The difference in results is very much negligible, which is the validation of the model.
In Figure 2(a) and (b), F1 is objective function for fixed cost and transportation cost and F2 is objective function for unutilized capacity cost. The sloped linear line is the Pareto front, and if we move along the points, we can either minimize fixed cost and transportation cost at the expense of unutilized capacity cost and vice versa. We cannot improve both at the same time. The negative slope also denotes it.
For this model, from the slope of Pareto curve, we can determine Pareto optimal solutions.
When we give equal weightage for two objectives, weight w1 = 0.5 for objective one which is minimizing fixed cost and transportation cost, and weight (1 − w1) = 0.5 for objective two which is minimizing unutilized capacity cost, we received the following results of Table 7.
Figure 3(a) and (b) indicates the proposed locations for two capacity ranges.
(Circled) Proposed storage facility locations for different capacity ranges (×105) (kg)
(Circled) Proposed storage facility locations for different capacity ranges (×105) (kg)
We can see that, if we decrease capacity range for each facility, we need to establish more facilities. We want to avoid the location which has highest probability of domino effect, which is Area 5 here. This Area 5 is not selected to establish storage facility. And both solvers give almost same results.
Table 8 shows the transported amount (×105) kg for capacity range taken random between 300 and 600 (×105) (kg).
Quantity transported (×105) (kg) to demand area j from a storage facility at location i (Gurobi solver and PuLP solver)
Quantity transported (×105) (kg) to demand area j from a storage facility at location i (Gurobi solver and PuLP solver)
In Table 8, the boxed areas are selected locations and the amount of transfer resources.
6. Conclusion and recommendation
In case of domino effect due to disasters, this model would give an idea of where to establish emergency storage facilities that will minimize the costs. In conclusion, we can answer the research questions as follows:
To understand the location where the probability of domino effect is highest, we have focused on estimating the probability of domino effect from the earthquake-risk zones. So, this model will help to preestablish the storage facilities considering the domino effect. This research has a positive impact on society that will affect quality of life by appropriate delivery of emergency resources. The work is thus bridging the gaps between the past research and practice.
The work also focused on the optimization model on minimizing fixed and transportation cost and unutilized capacity cost. Fixed and transportation cost are objective function one and unutilized capacity cost is objective function two. The cost is considered separately to understand the situation. Weighted sum method for multiobjective optimization has been applied, and it has been noticed that we have to either minimize fixed and transportation cost or unutilized capacity cost. We cannot minimize both at the same time.
The true benefit of the model is, the decision-makers of disaster management can use this in practice because the model is economical. The location where the domino effect probability is highest has been omitted from the selection as a high probability of domino effect indicates the high chance of damage if any facility is set up. The past researchers did not consider this phenomenon till now. We also vary the capacity range of the storage facilities. When the capacity decreases, we need to establish more facility to fill up the demand, so, the total cost increases. Thus, the decision-maker must be careful about selecting the capacity range.
Finally, Pareto optimal front has been determined for both Gurobi solver and PuLP solver. Both solvers give almost the same result, which validate the model.
While determining the domino effect probability, we needed to calculate escalation probability. Here, probit coefficients were taken after trial and error (which is one of the limitations) and to the best of our knowledge, there are no available optimal values of these coefficients in case of an earthquake. During trial and error, there were some other values for the coefficients, which were not considered (another limitation). So, there may be some possible directions to which this research can be extended. Further study can be done to search for the realistic values of the probit coefficients. More studies can be done to choose the optimal capacity range by adding more random capacity ranges to the calculation. Thus, the model can be more appropriate to the decision-makers.
Conflict of interest: The authors declare no potential conflict of interest.
