Two-Stage Stochastic Programming Model for Improving Transportation Network Resilience of Relief Supplies in Sequential Hazards Scenario

. Cascading failure of road transport networks caused by complex sequential natural hazards adversely a ﬀ ects the use of pre-positioning relief supplies. Therefore, it is vital to improve the transportation network resilience of relief supplies by using an optimization model. The purpose of this paper is to pre-position and distribute relief supplies in uncertain scenarios of sequential hazards. A two-stage stochastic programming model to maximize the total resilience is proposed to provide an optimal plan against the uncertain impact of sequential natural hazards. The combined impact of the Jiuzhaigou 7.0 magnitude earthquake and its associated landslide is a prototype disaster scenario for the implementation of the method proposed in this paper. The model is solved by a neighborhood search-based genetic algorithm (NS-based GA), which has both the global search capability of a genetic algorithm and the local search capability of a large-scale neighborhood search algorithm, can improve the so-lution ﬁnding capability. A case study focusing on ﬁnding the optimal solution for the pre-position and distribution of relief supplies in the sequential hazard of Jiuzhaigou earthquake is conducted to illustrate the validity of the proposed model


Introduction
The occurrence of natural hazards often brings massive losses to society. Most catastrophic hazards are caused by sequential hazards that occur sequentially instead of singly, resulting in more severe consequences [1]. The Great East Japan Earthquake that happened in the Tohoku region of Japan on 11th March 2011 is an example of a sequential disaster. The earthquake has a huge impact on Japan. Besides, it triggered a tsunami that devastated northeastern Japan and triggered a nuclear meltdown at the Fukushima Daiichi nuclear power plant [2]. More than 15,000 people were dead [2], around 130,000 residential buildings were demolished or washed away, and 260,000 houses and 59,000 non-residential buildings were partially destroyed in this disaster [3]. Hurricane Ida, which made landfall in coastal Louisiana on 29th August 2021, is another example of a sequential hazard. Hurricane Ida caused life-threatening

Literature Review
The word resilience originally comes from the Latin word "resiliere", which means "rebound" [13]. It refers to the ability of an entity or system to withstand, adapt to, and recover from damage. Different definitions of resilience can be obtained from different angles. The National Academy of Sciences [14] defined resilience as "the ability to prepare for, absorb, recover from, and adapt to disturbances. "The National Infrastructure Advisory Council (NIAC) [1] defined the resilience of an infrastructure system as the ability to anticipate, absorb, adapt to, and recover quickly from destructive events such as natural disasters. The cyber resilience discussed in this paper mainly refers to the ability of the network to maintain or restore its original function after suffering from sequential disasters and cascading failures.
At present, resilience is often used in the study of infrastructure networks such as power grids [15] and road networks [15], but it is rarely applied in the transportation network of emergency supplies [16]. In previous studies, there are generally two strategies to improve network resilience: (1) improving the physical performance of the network, such as improving road construction standards [15] and hardening of power grid lines [17]; and (2) improving the flexibility of the network, such as placing distributed energy sources in the grid in advance [11], and storing relief supplies in the network in advance [18]. Among the above countermeasures, the best way to improve the resilience of emergency supply transportation network is to pre-place relief supplies in the network and distribute them appropriately after disasters. So, a two-stage stochastic programming model is needed to find the best solution. Due to the diversity of network types and functional applications, there is currently no consensus on the specific evaluation indicators of network resilience. Kong [19] uses the difference between expected system performance and actual system performance as a measure of system resilience. Based on this indicator, this paper takes the difference between the expected material requirements of the system and the material requirements met by the system as the system resilience.
Because the model proposed in this paper has binary and integer types of decision variables and has multiple constraints, it can be regarded as a NP-hard problem to find an algorithm to solve. Since the solution of NP-hard problems is relatively complicated, its solving algorithm has always been a hot spot and difficult point in research. At present, most heuristic algorithms, especially improved heuristic algorithms, are used to solve NP-hard problems.
The main heuristic algorithms include genetic algorithm [20], ant colony algorithm [21], improved genetic algorithm [22], etc. Since the genetic algorithm is particularly prominent in system optimization, it is selected as the solution algorithm for the proposed model. For the problem of weak local search ability [20], we choose to combine the large-scale neighborhood search algorithm with strong local search ability to improve the reliability of the solution [23].

Problem Statement
According to the classification of earthquake damage to lifeline engineering, we can divide road ground damage caused by earthquakes into five levels. Furthermore, low-level damage will not affect the transportation of relief supplies, moderately damaged roads can be used after a short period of repair, and severe damage will make the road to be fixed for a long time before available [24]. No matter what degree of road damage occurs, it will adversely affect the transportation of relief materials, make it impossible for the affected points to get the required materials in time. Besides, it will impact the resilience of the disaster relief logistics network. In this paper, a two-stage stochastic programming model is proposed to optimize the above problems.
Due to the high uncertainty and complexity of earthquakes, the demand for various locations and the extent of road damage is uncertain during the planning stage. We use a set of discrete scenarios simulated with the Jiuzhaigou earthquake as a prototype to represent the uncertainty. The composition of the scenario includes the need for relief supplies of each node, the probability of occurrence, road damage due to earthquakes, road damage due to landslides. The anti-earthquake characteristics of the warehouse make it nearly impossible to damage in an earthquake.
After the earthquake, pre-placed relief supplies are transported through a traffic network to satisfy the demand. Each warehouse serves only one node at a time. Due to the damage of roads caused by sequential disasters and the interdependence in a traffic network, the cascading failure of other roads will occur during the transportation of materials. Furthermore, the occurrence of landslides that caused by earthquakes may again cause the roads to be damaged, nodes to generate demand, and network resilience to be down. The aim of this paper is to determine the best strategy for maximising the resilience of the disaster relief logistics network. In doing so, we determine the location of warehouses, the reserve amount of relief supplies, and the relief supply allocation decisions during and after sequential disasters. Meanwhile, we also consider the uncertainty of relief supplies and road damage in sequential disasters. Due to the long distance between nodes, it is difficult for vehicles to travel back and forth in a short period of time. Therefore, we choose a longer time period for analysis when calculating cascading failures. In this paper, we have only used two points in time, before the landslide occurred and at the end of the relief transport, for cascading failure analysis.

Two-stage Stochastic Programming Model
The problem described above is represented in this paper as a two-stage stochastic programming model. The first stage focuses on making decisions about the location, capacity and quantity of goods to be stored in the warehouse. The second stage is to make decisions on the route and quantity of relief supplies to be transported. The symbolic representation of the model is as follows.
The sets in the model are defined as follows: I Set of potential locations where relief supplies can be stored, indexed by i; J Set of demand locations that have needs for commodities, indexed by j; S Set of possible sequential disaster scenarios, indexed by s; K Set of storehouse size categories in potential locations, indexed by k; T Set of the overall time span when the sequential hazards effect lasts, indexed by t.
The parameters in the model are defined as follows: B Available investment for setting up storehouse and purchasing initial supplies; C Capacity of each car; W k Capacity of k category storehouse; C k f ixed Fixed cost of building a k category storehouse; C transp One unit distance transportation cost for per unit relief supplies; C p The cost of purchasing one unit of relief supplies; P s The occurrence probability of scenario s; D i j The distance from supply location i to demand location j; d s jt Target demand for relief supplies of location j at time t in scenario s; d s jt Actual demand that has been satisfied of location j at time t in scenario s; L ts i j The load that transported from supply location i to demand location j at time t in scenario s; e ts i j The finite capacity of traffic route from supply location i to demand location j at time t in scenario s; R The distance limitation of a car. The decision variables in the model are defined as follows: x ts i j The quantity of supplies shipped from supply location i to demand location j at time t in scenario s; y k i A binary variable equal to 1 if a k category storehouse is built at supply location i, and 0 otherwise; z ts i j A binary variable equal to 1 if transport from supply location i to demand location j at time t in scenario s, and 0 otherwise; (1) Subject to: The objective function 1 maximizes the expected total resilience over all scenarios. The total resilience represented as the ratio of satisfied demand to target demand. Constraint 2 assures that the cost consisting of setting up storehouse and purchasing initial supplies does not exceed the total budget. Constraint 3 the satisfied demand in time t equal to the quantity shipped from the warehouse to. Constraint 4 is to limit the number of warehouses in one location. Constraint 5 states that the load in a side must equal to the initial load plus transport load. Constraint 6 and Constraint 7 enforce the quality shipped from the warehouse and pre-positioned not exceeds the capacity of warehouse respectively. Constraint 8 makes that the quantity shipped out of the warehouse does not exceed the quantity in existing storage. Constraint 9 is the distance constraints. Constraint 10 is to make sure that the satisfied demand will not exceed the real demand. Constraint 11 and Constraint 12 are restrictions on decision variables.

NS-based GA
An NS-based GA is used to solve for the proposed two-stage stochastic programming model. The algorithm can improve the local search capability using a large-scale neighborhood search algorithm while retaining the global search capability of the genetic algorithm. Only two time points in the study, before the occurrence of the landslide and at the end of transport, may experience cascading failures. Therefore, only these two time points need to be considered in the solution to see if cascading failures occur. The steps of this solution method are as follows: Step 1. Coding of decision variables for each scenario for the time period after the earthquake and before the landslide occurred.
Step 2. Compute initial solutions for the coded chromosomes.
Step 3. Select, crossover and mutate the initial solutions.
Step 4. Apply the large-scale neighborhood search algorithm to the evolved chromosomes.
Step 5. Repeat steps 3 and 4 above until a solution satisfying the conditions is obtained.
Step 6. Determine if any roads have cascading failures, and if so, set the capacity of the road to 0.
Step 7. Code the decision variables for each scenario for the time period between the occurrence of the landslide and the end of transport.
Step 8. 7. Repeat step Step 2. to Step 6. above. The type of decision variables is the main factor that determines the encoding method of chromosomes in the genetic algorithm. The decision variables in the model of this paper include x ts i j , y k i and h i . The coding method will be chosen with respect to the meaning and type of each decision variable.
Since x ts i j is an integer variable and represents the quantity of relief supplies delivered, it is coded as an integer, and the coded value is the quantity of relief supplies delivered. y k i is coded as a binary variable since it represents the type of warehouse built and is of the binary type. A code of 1 for y k i means that a warehouse of type k is built at point i. A code of 0 for y k i means that a warehouse of type k is not built at point i. h i represents the quantity of relief supplies stored in the warehouse and is an integer variable, so it is also coded as an integer. The final structure of chromosome is shown in the figure 1.
The fitness function determines the degree of superiority or inferiority of each chromosome. Generally speaking, the objective function can be directly used as the fitness function [25]. In this paper, we use the calculation of the objective function to perform the adaptation degree calculation. For the chromosomes encoded in the time period from the occurrence of the earthquake to the occurrence of the landslide, the elasticity of this period is calculated as f s 1 . For chromosomes encoded in the time period from the occurrence of the landslide to the  (13) In addition, the determination of some important parameters in the genetic algorithm is also very important. Vital parameters including the initial population size, crossover probability, variation probability and the number of iterations. The initial population size should not be too large or too small, so it is set at 100 in this paper. Crossover probability and variation probability will have an impact on the search ability of the algorithm. It is generally larger for crossover probability and smaller for variation probability, so they are set at 0.9 and 0.1 in this paper. The number of iterations should also not be too large or too small, so it is set at 200 in this paper.
Constraints are set for the generation of some of these decision variables in order to improve the quality of the initial solution as well as to better satisfy the constraints in the model. This allows the generated initial solution to naturally satisfy some of the constraints. The constraints include the following four conditions. One is that only a maximum of one type of warehouse is built at each supply location. Secondly, the pre-placement of supplies in the warehouse does not exceed the capacity of the warehouse. Thirdly, relief supplies will only be delivered to locations where there is a need. Fourthly, the supplies delivered by the warehouse will not exceed the pre-stored supplies in the warehouse. The next step is to determine if the remaining constraints are satisfied until an initial solution is obtained that satisfies all the constraints.
The tournament method is used in the selection step of the genetic algorithm species. Two chromosomes are selected from the population at a time for comparison of fitness values, and the individual with the greater fitness value is selected to enter the next generation of optimization. To ensure that the solution after this step is still a feasible solution, we use a two-slice crossover to exchange the segments of the two parents where t and s are the same. The swapped segments are from x t s i 1 to x t s i j . Also, to make the crossed solutions also satisfy the constraints, the decision variables y k i and h i need to be swapped at the corresponding k and i. The part of the crossover of the decision variable x ts i j is shown in figure 2. In this paper, to ensure that the solution after this step is still a feasible solution, we use a mutation with two segment exchange. The two fragments from the decision variable are exchanged with those from x t s i 1 to x t s i j . Also, to make the solutions after mutation also satisfy the constraints, The y k i and h i decision variables need to be transformed accordingly. The part of the mutation of the decision variable x ts i j is shown in figure 3. For the destruction session in the large-scale neighborhood search algorithm, we use the approach of randomly deleting one x ts i j decision variable fragment and moving all the behind x ts i j forward. The repair link is a way to insert the deleted fragment into the last part of the x ts i j decision variables. Where the deleted segment is an arbitrary selection of t and s identical, from x t s i 1 to x t s i j . In addition, the same is needed for y k i to move in correspondence with the h i decision variable. The part of the large-scale neighborhood search algorithm of the decision variable x ts i j is shown in figure 4.

Data Preparation
In this paper, the 7.0 magnitude earthquake data of Jiuzhaigou is applied to generate sequential disaster scenarios to verify the effectiveness of the two-stage stochastic programming model. The magnitude of earthquakes is usually measured by intensity, and the intensity distribution of the area near Jiuzhaigou in this earthquake is shown in the figure 5 (from http://news.ceic.ac.cn/). In figure 5 we use green closed curves represent isoseismic lines. The four green closed curves from inside to outside in the figure represent the intensities of IX, VIII, VII and VI, respectively. Red circles in figure 5 represent major cities and counties in the region that were affected by the earthquake. Colored curves with double borders in figure 5 represent the main roads in the area [26]. In this paper, the 14 affected areas are used as both demand locations and candidate locations for setup warehouses. The location names and indexes of that areas are shown in the table 1. Study areas can be thought of as a directed network diagram. In this network, nodes can reach each other. The distance data between two directly connected nodes comes from Baidu Map. Detailed distance data information is provided in table 2 in kilometers. When relief supplies are transported between two nodes, the shortest route between the two nodes is used. The vehicle needs to go back and forth between two nodes in a single transport task. According to Chinese road regulations, the speed of vehicles is limited to 30 to 120 kilometers per hour. The average speed of the vehicle is set to 30 kilometers per hour, and the maximum distance for the transportation of emergency supplies is 180 km. In this paper, for the convenience of calculations, only the case of water as relief supplies is discussed. When calculating, the units of water are unified into cubic meters. In (3,5) 120 18 (8,9) 387 27 (13,14) 142 the transportation of relief supplies, trucks with a capacity of 50 cubic meters are used. The average person needs to drink 0.002 cubic meters of water a day (available from http://dg.cnsoc.org/article/04/wDCyy7cWSJCN6pwKHOo5Dw.html). The price of purchasing a unit of water is $118 (available from https://www.shuishang wang.com/jiage/pinpai/). In this paper, there are three capacities of warehouses to choose from, with a capacity of 70, 90 and 110 cubic meters, and the costs are $1924, $3702 and $5923, respectively [18]. According to the demand for relief materials and the cost of building the warehouse, the funds available for the construction of the warehouse and the purchase of relief materials are set to $80,000.
According to the highway vulnerability function and related data, the probability that each road cannot be used under the Jiuzhaigou earthquake can be obtained [27]. In addition, 11 possible road damage situations and their probability of occurrence are calculated. In this paper, each situation is used as a scenario. The probability of earthquake-triggered landslides occurring is determined based on the incidence of Jiuzhaigou earthquake-triggered landslides in previous studies [28]. In these scenarios, earthquake-induced road failures occur at the very beginning of the study, and landslides-induced road failures are assumed to occur in the middle of the study. Table 3 lists the characteristics of the 11 scenarios, including probability of occurrence, road damage due to the earthquake, road damage due to landslides, and the total water demand at each node. In table 3, earthquake represents the index of failed roads caused by the earthquake, and landslide represents the index of failed roads caused by the landslide. The number of victims at each node is distributed according to the total number of disaster victims [6]. The distribution is based on the weighted average intensity of each node under the 7.0 magnitude earthquake in Jiuzhaigou. The demand for relief supplies at each node is calculated based on the number of people affected by the disaster at each node and the demand for relief materials by a single person.

Results Analysis
The solution of the two-stage model proposed in this paper is programmed using the Python programming language, and NS-based GA. The results of solving the decision variables in the first stage are shown in table 4. A total of 6 warehouses have been established, 2 warehouses have a capacity of 110 cubic meters, 3 warehouses have a capacity of 70 cubic meters, and 1 warehouse has a capacity of 90 cubic meters. This indicates a higher use of small warehouses and a lower use of medium and large warehouses in the results. In addition, the space utilization rate of each warehouse is more than 70 percent, and the highest even reaches 94 percent.
In the second stage of the model for the decision variables, we look primarily at the results for transportation volumes. For each node, it is usually preferred to have its own requirements met by the supplies in the warehouse established by its node. If the node has not established a warehouse or the supplies stored in the node's warehouse cannot meet its demand, the supplies need to be transported from other nodes' warehouses. The results of solving the second stage decision variables for the first five scenarios are shown in table 5. Where arc is used to indicate the route along which the transportation of supplies between the nodes occurs in this scenario. The total objective function value for the model solution results is 89%, meaning that a total of 89% of the supplies requirements are met in the specified time.
Combining the results of the two phases shows that the demand nodes use more relief supplies directly from the warehouse where their node is located, and use less relief supplies transported by other nodes' warehouses. This shows that the model not only improves the efficiency of relief supplies distribution, but also reduces the transportation cost. And by observing the transportation routes in table 5, it can be observed that there is a clear tendency to use directly connected routes when transporting relief supplies. In addition, by looking at the transport volumes in table 5 it can be seen that most vehicles have low-capacity utilization.

Conclusions and Future Research
To improve the resilience of the relief supplies transportation network under sequential disaster scenarios, a two-stage stochastic programming model is proposed in this paper. This model takes into account the uncertainty of sequential hazards using a discrete scenario approach. This paper takes into account road failures due to earthquakes and road failures due to earthquake-induced landslides. During the transportation of relief supplies, cascading failures to avoid failed roads due to sequential disasters are taken into account. Moreover, the effectiveness of the model proposed in this paper is verified through the application of a sequential disaster scenario generated using the Jiuzhaigou 7.0 magnitude earthquake as a prototype. For the solution of this model, this paper uses NS-based GA. However, as this study only considers the maximization of resilience and not vehicle utilization, some vehicles are under-utilized. It is therefore necessary to consider the combination of network resilience and vehicle utilization in subsequent studies. Meanwhile, the priority of each node's demand needs to be considered in the actual transportation. Priority needs to be given to the part of the node with the most pressing needs during the transportation of relief supplies. In addition, this paper has found in case applications that some roads will suffer cascading failures due to overloading when transporting relief supplies. Therefore, when transporting relief supplies, it is important to make reasonable use of road capacity and plan distribution routes rationally to reduce the possibility of cascading failures. For the types of sequential hazards, only earthquakes and earthquake-induced landslides are considered as sequential hazard scenarios in this paper. However, the situation is often more complex in actual disasters, so three or more hazards can be considered as sequential hazard scenarios in subsequent studies.