Three-Stage Mixed Integer Robust Optimization Model Applied to Humanitarian Emergency Logistics by Considering Secondary Disasters

In recent years, natural disasters occur frequently, and secondary disasters induced by major disasters will also cause huge losses. The diversity of secondary disasters makes humanitarian emergency logistics (HEL) more challenging but often overlooked by researchers. In order to solve the comprehensive HEL problem of major and secondary disasters, a three-stage mixed integer linear optimization (TS-MILO) model is proposed. Among them, the uncertainty of the demand for relief supplies is also extremely difficult to deal with. In order to resist the interference of uncertainty, based on robust optimization, the TS-MILO model is further transformed into a three-stage mixed integer robust optimization (TS-MIRO) model, which are respectively BTS-MIRO (Box set), PTS-MIRO (Polyhedral set), and ETS-MIRO (Ellipsoid set). The experimental results show that the TS-MILO model can provide the lowest cost but cannot solve the uncertainty problem. The improved TS-MIRO model will pay a robust price (increase by at least 10.05%), but maintains supply stability even in the worst-case scenario. Specifically, ETS-MIRO model has strong robustness, and its cost increase only accounts for 44.66% of BTS-MIRO model in the worst case. The service level of the three TS-MIRO models increases with the safety parameters, among which the service level in the ETS-MIRO model increases significantly from 88.53% to 96.44%. The research results can provide a strong support for the decision making of disaster relief management department.


I. INTRODUCTION
In recent years, large-scale natural disasters have frequently occurred, which resulted in a large number of casualties and property losses. In addition to the great damage caused by major disasters, the enormous losses brought about by subsequent secondary disasters cannot be ignored. The Hua-xian earthquake (M = 8 1 ; Jan., 1556) is the largest earth-quake recorded in Shanxi, CHN, in which secondary disa-sters such as landslides and floods are induced by major disasters. According to historical records, 830,000 people were killed [1]. The Xingtai earthquake (M = 6.8; Mar., 1966) occurred in Heibei, CHN. Secondary disasters such as landslides, collapses, flood and water spraying occurred, which The associate editor coordinating the review of this manuscript and approving it for publication was Nagarajan Raghavan . 1 M indicates the magnitude of the earthquake, with destructive and seismic intensity increasing progressively from I to XII. submerged large quantities of farmland and water conservancy facilities [2]. The San Francisco earthquake (M = 8.3; Apr., 1906) led to fire and other secondary disasters [3], [4]. The Kobe earthquake happened in JPN (M = 7.2; Jan., 1995), which resulted in the cracking of natural gas pipelines and further triggered secondary disasters such as fires [5], [6]. Besides, an earthquake (M = 6.8; Jul., 2006) occurred in the southwest of Java, IDN, causing a secondary disaster tsuna-mi [7]. The earthquake in Nepal killed 8,790 people, injured 22,300 people and damaged millions of buildings [8]. Secondary disasters refer to a series of disasters that followed strong earthquakes, such as fire and tsunami. Those events indicate that the losses caused by secondary disasters may even exceed those of major disasters. Con-sequently, we must pay high attention to secondary disa-sters.
HEL plays an indispensable role in post-disaster rescue. At present, scholars usually concentrate on the major disasters but neglect the secondary disasters. As an VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ extraordinary specialized emergency logistics, HEL can satisfy the de-mands of emergency materials like food and medicine in the post-disaster rescue, and play a very important part in reducing casualties and property losses [9]. Bussieres et al. (2017) believe that the implementation of the distribution system can achieve the target of large-scale material supply after major disasters [10]. Benini et al. (2003) used geo-graphic information system to analyze rescue data in major disasters and concluded that the benefits are considerable, but the costs are also high, especially the hidden costs [11]. Olness et al. (2005) consider training programs for rescue personnel to be an effective tool to improve rescue capa-bilities [12]. Chen et al. (2008) used artificial immune algori-thm to confirm the effectiveness of the HEL vehicle path planning model in major disasters [13]. According to the actual situation, there are some challenges in HEL management. Potential secondary disasters can easily inter-fere with the rescue plan, which makes the rescue in the disaster areas more difficult. The previous research often focused on the rescue work after a single stage of major disasters, whereas few scholars paid attention to secondary disasters. Therefore, this paper attempts to meet the rescue demands of major disasters, further planning to satisfy the rescue needs of secondary disasters, and to improve comprehensive rescue capacity of HEL. Considering the impact of secondary disa-sters and the limitations of current research, the purpose of this study is to explore the relationship between major and secondary disasters, and to develop models to respond to uncertain disasters.
In the research of HEL, the previous literature often focused on certain scenarios, while few scholars carried out research on its uncertainty. The uncertainty of the HEL process has a great impact on the rescue plan. It is not only reflected in the uncertainty of transportation demands, but also in that of the supply of relief materials. Because of the uncertainty, it is difficult to obtain accurate information. Secondary disasters are mainly caused by major disasters. Besides, the formers are more difficult to predict in advance and will undoubtedly be more laborious. Some scholars have analyzed and modeled the uncertainty in the actual operation management.  proposed an opportunity con-strained programming model to settle the energy storage problem under uncertainty [14]. Balcik and YanıKoğLu (2020) studied the uncertain routing problem by constructing a robust duration constraint model [15]. These studies have proved that compared with the certainty model, the uncertainty model method is more suitable for the actual situation and has more practical application value.
In view of the uncertainty in the HEL, scholars and experts have carried out research by using heuristics, fuzzy program-ming and stochastic programming. As for heuristics, some scholars have established models to figure out the total costs of all possible factors in HEL. However, these studies ignored the influence of possible secondary disasters. Xu et al. (2011), using ant algorithm under the condition of determined demand, verified the effectiveness of the model with the objective of minimizing supply imbalance [16]. Liu et al. (2016) studied the disruption management issue of post-earthquake HEL system, and worked out the problem via employing hybrid heuristic algorithm [17]. Zhang et al. (2013) proposed a new biome-metic method for the path selection of emergency logistics under the condition of known demand [18]. Zeng et al. (2014) proposed a neighborhood descent algorithm to resolve the problem of cumulative multi warehouse vehicle routing in emergency logistics under the condition of inventory determination [19]. Wang et al. (2015) established a mini-mum objective mathematical model in HEL, and used two-stage heuristics to figure out the problem [20]. And with regard to fuzzy programming and stochastic programming, there are still few scholars who have studied the HEL under uncertain conditions. For the sake of solving the issue of uncertainty, Xu et al. (2016) proposed an ambiguous two-level model and verified the practicability of the model [21]. Yang et al. (2008) proposed a new twostage stochastic programming model to solve the problem of relief allocation [22]. Dai and Ma (2008) proposed hybrid clustering optimization model for large-scale disaster relief with the aim of minimizing the total number of deaths in major disasters [23]. Shen et al. (2019) proposed a triangle model to meet the demand of emergency logistics, and made use of a hybrid two-stage algorithm to solve it [24]. On the whole, the above scholars have studied the uncertainty decision-making problem of HEL from various dimensions. Heuristic algorithm, fuzzy programming and stochastic programming have developed more syste-matically, but there are still some limitations.
• The deterministic linear optimization model only focuses on the single stage decision-making problem.
• It is tough to offer a theoretical exact solution through the heuristic algorithm model, which depends on the selection of termination criteria.
• It is difficult to describe the actual situation by means of the linear chance constrained model. Besides, the solution easily leads to NP-hard problem, is complex.
• It is also hard to describe the membership function via stochastic fuzzy model, which results in lower feasibility of the results.
• Due to the lack of sufficient historical data in stochastic programming model, the blind assumption of probability distribution has high risk. Therefore, it is necessary for us to conduct in-depth research on the uncertainty of HEL and seek a more compati-ble and more practical research methodology.
With the rapid development of theory, robust optimization has been proved to be able to overcome the limitations of stochastic optimization in uncertain problems [25]. Although the research on robust optimization theory is relatively ma-ture, its application still needs to be further enriched, espe-cially in the field of HEL. The advantage of robust optimi-zation is that the model performs reasonably well, even if we estimate the cost at the worst case, and that a robust solution can be acquired. Robust optimization theory was first put forward by Bertsimas et al [26]. As a tool for optimization under uncertain conditions, robust optimization does not depend on probability distribution or membership function, but represents uncertainty through a specific set [27]. Some scholars have used its theory to study decision making problem of path planning. In the case of uncertainty, Zhang proposed a two objective robust optimization model to re-solve the emergency medical service problem [28]. Chen et al. constructed the two-stage robust security-constrained unit commitment model and applied it to a number of research areas [29]- [32]. Through actual cases, it has been proved that it can effectively reduce the conservatism of uncertain sets and ensure the economy of optimization results. Robust optimization can take the uncertainty of parameters and the heterogeneity of risk preference of decision makers into consideration, which will undoubtedly provide valuable empirical suggestions for the research of HEL. Some school-ars have applied the robust optimization to the research of supply chain management [33], [34] and decision planning [35]. Jadidbonab et al. (2020) adopted the non-probability information interval method to optimize the self-scheduling problem under the uncertainty with the goal of maximizing the benefits under the circumstance of uncertainty of renew-able resources [36]. Gholinejad et al. (2019) proposed an energy management system for multiple home energy hubs in a community, and the simulation results have displayed the effectiveness [37]. Nazari-Heris et al. (2019) proposed a multi-objective two-stage random unit model, which was applied to a new flexible energy gas electric networking system [38]. In terms of new energy applications, Mirzaei et al. proposed a great many robust optimization models for uncertain energy demand and variable wind direction, and verified the performance of the robust optimization model [39]- [41]. Al-Sumaiti et al. (2019) proposed a more pro-mising approach to meet the demands of power supply by considering the impact of the probability uncertainty nature of weather [42]. As is mentioned above, robust optimization has been applied to the fields of supply chain management, decision planning and new energy by experts and scholars from all walks of life. However, in the research of HEL, few people use the robust optimization. Therefore, it is innovative to apply the robust optimization theory to researching HEL.
In addition to the application of robust optimization, this paper further extends the single-stage HEL path planning to three-stage. Through careful observation of the actual activities, we can find that in the rescue process, the relief materials are not only transported from the origin to the disaster relief point, but there are also logistics operations of multiple transshipment. However, the single stage model is too idealistic to reflect the reality. Thus, it is necessary to explore and study the multi-stage model. Most of the related researches only involve the design of one or two of the three sub-stages. There is no research on integrating these three sub stages to systematically solve the path planning problem of HEL. In addition, the relationship among the three sub-stages is interdependent, since some inputs of the first stage planning model are the outputs of the second stage model and some inputs of the second stage model are obtained from the third stage. Therefore, the sub-stages must be designed simultaneously with the overall planning model. This method can be used to design a seamless emergency rescue supply network integrating three disaster ma-nagement sub-stages. In the field of application, especially in the HEL, it has not been found whether the two-stage or even multi-stage problem model can be solved by robust optimization.
In summary, the contributions and innovations of this research are as follows. In general, this research can play a guiding role in the construction of the model for HEL, and help decision makers to formulate an appropriate emergency relief strategy to respond to natural disasters.
• This paper not only focuses on the impact of major disasters, but also further analyses the impact of secondary disasters and the internal relationship between them.
• Based on the comprehensive consideration of three substages of disaster management, a three-stage mixed integer linear optimization (TS-MILO) model is proposed.
• Robust optimization theory is applied to the research of three-stage HEL problem. TS-MILO model is further transformed into TS-MIRO model.
• According to the real location data of HEL, computer simulation analysis is carried out. This model provides valuable decision support for government relief departments. The rest of this paper is organized as follows. In Section II, the description of the problems and TS-MILO model are given. In Section III, the TS-MILO model is constructed and it is further transformed into TS-MIRO model in Section IV. The contents in Section V and VI are about the numerical analysis. In Section VII, the conclusion and future research are put forward.

II. PROBLEM DESCRIPTION
Fire, tsunami and other secondary disasters are mainly brought about by major disasters such as earthquakes, and the occurrence of major disasters is the premise of their generation. HEL planning can be divided into three stages. The HEL rescue process of major disasters includes emer-gency preparation stage (the first stage) and rescue stage (the second and third stage). At the same time, on the basis of meeting the rescue requirements, reasonable planning and allocation of resources can reduce the loss caused by unnecessary operation. Decision making in the first stage is influenced by the results of the second stage. Similarly, the decision making in the second stage is affected by the third stage ( Fig. 1).
In HEL emergency rescue network, there are I original warehouses, J major disaster sites and K secondary disaster sites (Fig. 2). First of all, TS-MILO model is established. The purpose of the model is to explore how to reduce the total costs as many as possible on the basis of meeting the demand to the maximum extent. The first stage is the post-disaster preparation stage. In this stage, under the circumstance of only knowing the location of the disaster, it is necessary to  estimate the demand of disaster relief materials, and then select the location of disaster relief sites. The selection of rescue nodes should meet the demand of the maximum service and the feasibility of practical operation. The second stage is the material distribution stage of the main disaster areas. The goal of this stage is to solve the initial path planning problem of emergency supplies. On the basis of comprehensive consideration of the practical feasibility (cost and capacity), the rescue materials are transported from the original warehouses to the main disaster affected areas, in which the demands of emergency materials are not only affected by their own demands, but also by the demands of the third stage of secondary disaster. The third stage is the material distribution stage of the secondary disaster area. The goal of this stage is to solve the problem of vehicle routing in secondary disaster areas. Considering the potential demand of the secondary disaster areas, the surplus materials should be allocated to the secondary disaster areas after meeting the demand of major disaster.
The three-stage HEL rescue network has the following two functions. On the one hand, overall planning needs to ensure that humanitarian relief needs are satisfied to the maximum extent. It includes not only the demands of the main disaster areas, but also the material demands of the secondary disaster sites. On the other hand, it is also necessary to reduce the costs of operation management through reasonably plan-ning path. The costs are divided into three types: fixed operating costs, vehicle transportation costs and time costs [43], [44].

A. ABBREVIATIONS AND ACRONYMS
The specific parameters are shown in Table 1.

B. ASSUMPTION
Combined with the above analysis, following assumptions are made for the three-stage HEL problems.
• It is assumed that in addition to meeting their own needs, the major sites can also supply materials to multiple second-ary disaster relief sites.
• The rescue demands of all sites must be met. And any rescue demand site has at least a supply site to supply.
• The vehicles used at the same stage are of the same type with the same fuel consumption and load capacity.
• The geographical location and time window of candidate initial disaster relief sites are known.
• The average speed of vehicles in different regions are different, which depends on road conditions and time period in the area.

III. MODEL CONSTRUCTION
In this section, the basic three-stage mixed integer linear optimization (TS-MILO) model is constructed to analyze the three-stage HEL problem deeply.

A. TS-MILO MODEL
According to the actual situation of HEL, the specific costs involved include: the fixed cost, handling cost, transportation cost and path time cost [45], [46]. TS-MILO model is constructed, which aims to meet the demands of human-itarian relief, achieve the goal of selecting material storage sites and calculating the minimization of overall costs. According to the real situation, the first stage of TS-MILO is shown below. (1) Among them, the first item of (1) is fixed cost, which is infrastructure investment cost, including office equipment loss cost and basic water and electricity cost. Fixed cost has nothing to do with vehicle routing planning. The second is the handling cost, which is affected by the demand. The third cost is affected by the second stage of vehicle routing planning.
Specific constraints: (2) represents that the total handling capacity cannot be higher than the total demand, and there is no other outflow part; (3) shows that the maximum handling capacity is less than its maximum capacity; (4) only the selected initial sites participate in the corresponding logistics operation; constraint (5) is related variable constraint, in which only participates in the corresponding emergency material distribution when x j = 1.
In the second stage of TS-MILO, the transportation cost is minimized on the basis of meeting the demand of initial sites. Among them, the first item of (6) is the vehicle transportation cost of the second stage. The second is the time cost, and the third one is related to the path planning in the third stage. Specific constraints: (7) represents that the total distribution volume cannot be higher than the total demand; (8) represents the maximum storage capacity; (9) represents time constraint; (10) only the selected initial node participates in the corresponding emergency logistics operation; (11) is related variable constraint.
In the third stage of TS-MILO, the purpose is to minimize the total cost on the basis of meeting the needs of secondary disaster relief sites to the maximum extent. The first term of (12) is the time penalty cost of the third stage, and the second term is the vehicle transportation cost. Specific constraints: (13) represents the total distribution volume unable to be higher than the total material supply; (14) represents the maximum storage capacity; (15) represents time constraint; (16) is related variable constraint.

IV. CONSTRUCTING TS-MIRO MODEL
Based on the TS-MILO model, this section further constructs three robust optimization models under uncertain conditions, which are three-stage mixed integer robust optimization model based on box set (BTS-MIRO model), three-stage mixed integer robust optimization model based on polyhedron set (PTS-MIRO model), and three-stage mixed integer robust optimization model based on ellipsoid set (ETS-MIRO model).

A. UNCERTAINTY SCENARIOS
In the process of post disaster rescue, there is no ideal model. On the contrary, the external environment is full of complexity and uncertainty. It is very difficult to obtain the precise value or probability distribution of key parameters, especially the demand parameters [47]. The idealization of theoretical research may lead to low degree of feasibility in reality, in other words, the stability of deterministic model under uncertain conditions is infeasible. Therefore, we intro-duce the concept of robust optimization. In reality, the problem to be solved is uncertain, which is reflected in the uncertainty of parameters. Robust model can provide an effective uncertainty measurement method, so the research of robust optimization has high applicability. Applying the theory of robust optimization [48], the TS-MILO model is transformed into TS-MIRO model. The uncertain parameters change in uncertain sets, so it can be studied without depending on the probability distribution. The greater the volatility of the first type of initial rescue sites demand is, the greater the uncertainty and uncontrollability will be. The uncertain parameter Dj 0 = Dj +Dj, where, Dj 0 is nominal value,Dj = εDj 0 is fluctuation demand and ε is disturbance proportion. The demand fluctuation of secondary disaster stations is high and uncontrollable. It is defined as a random demand parameter k is demand fluctuation and ξ is disturbance proportion. Therefore, BTS-MIRO, PTS-MIRO and ETS-MIRO models are established, respectively [49].
Theorem 1: Under uncertain situation, when the box uncertainty set is empty, the BTS-MIRO model will degrade to TS-MILO model. When the box uncertainty parameter set is non-empty, (1) - (16) robust equivalent model is (17) - (34).
The first stage of BTS-MIRO is (17) - (22), and its goal is to minimize the total cost under uncertain conditions.
x j ∈ {0, 1}, y ij ≥ 0, ∀i ∈ I , ∀j ∈ J The second stage of BTS-MIRO is (23) - (28), and its goal is to minimize the distribution cost of emergency materials on the basis of maximizing the demand of major disaster relief sites.
The third stage of BTS-MIRO is (29) - (34), and its goal is to minimize the cost of secondary disaster rescue on the basis of maximizing the demand of secondary disaster sites.

C. PTS-MIRO MODEL
In the PTS-MIRO model, the uncertain set of demand sites is polyhedron set, which is defined by l 1 norm: D :  (52).
The first stage of PTS-MIRO is (35) - (40), and its goal is to minimize the overall rescue cost.
The second stage of PTS-MIRO Model is (41) - (46), and its goal is to minimize the initial rescue cost.

D. ETS-MIRO MODEL
In ETS-MIRO model, the uncertain parameters float in the ellipsoid set. Through the l 2 norm, U E1,E2 are defined as: {ς, ζ } : ε 2 , ξ 2 ≤ j,k = ε J ε 2 j , ξ K ξ 2 k ≤ j,k , in which j,k is the adjustable SP and also the spherical diameter of uncertain set [56].
the above formula is converted to C(D 0 j,k ) + j,k P j,k ≤ Z E . Since the goal is to minimize the cost, relaxation constraint

Theorem 3: Under uncertain situation, when the ellipsoid uncertainty is empty, the ETS-MIRO model will degrade to TS-MILO model. When the box uncertainty parameter set is non-empty, (1)-(16) robust equivalent model is (53) -(76).
The first stage of ETS-MIRO is (53) -(60), and its goal is to minimize the total cost under uncertain conditions.
The second stage of ETS-MIRO is (61) -(69), and its goal is to minimize the cost of major disaster relief.

E. ALGORITHM DESIGN
Through MATLAB (R2016a) as programming platform, the algorithm process of model is shown (Table 2). All data analysis is carried out on the same PC, which has 2.44 GHz 4-core CPU (Inter Core), 256GB SSD and 8GB RAM.

V. NUMERICAL EXPERIMENT
In this section, a real case is used to verify the effectiveness of the models in HEL. The model selects Yushu Earthquake (M = 7.1, Apr., 2010) as an example (Fig. 3). In actual rescue process, the Rescue Department is faced with the HEL problem. The first stage is to determine the location of major disaster relief sites, which are determined according to the disaster situation. Taking Jiegu street as the center, the location of the main disaster relief sites is determined first. The goal is to determine the location and minimize the overall total cost. These major sites offer dual functions. One is to provide materials for their own disaster areas, and the other is to provide materials for the secondary disaster areas. The original center of relief materials is determined as O 1 of Batang Airport (Cargo mark), and Jiegu, et al. 2 towns are selected as the major disaster relief sites (Red Cross), which are M 1 , M 2 , · · · , M 5 . The second stage is the path planning from the original center to the major sites. There are five major sites, which are selected among the alternative sites in the first stage. The aim of the second stage is to minimize the initial distribution costs. The third stage is the path planning from the major sites to the secondary sites. Secondary disasters are often caused by the follow-up effects of major disasters, such as small-scale aftershocks, landslides, dammed lakes, floods and so on. There are 14 secondary sites (Black solid), namely as Lixin, et al., 3 which are represented by S 1 , S 2 , · · · , S 14 . The aim of the third stage is to minimize the costs from major sites to the second sites.
In the process of route planning, we remove the interference factors and get the relative position of the sites (Fig. 4). On the basis of maximizing the satisfaction of demand, the minimum costs are pursued. In the HEL system, there is one original center, 5 major sites and 14 secondary sites. Any alternative path corresponds to different costs. The transportation costs are determined by the comprehensive calculation of real-time oil price, actual distance, and traffic congestion and time limit. On the basis of comprehensive consideration of relevant costs, the following case simulation is carried out.

A. RELATED DATA
The data are mainly from China Seismological Bureau and National Bureau of Statistics. In this paper, the input parameters are estimated using historical data of a disaster-prone area (Yushu, Qinghai Province). The basic information includes information of vehicle (Table 3). The fixed operating costs, demand for major sites, demand for secondary sites, and vehicle speed ( Table 4). The distance between sites can be obtained through Google map [57], [58] as is shown in Table 5.    Fig. 5. We can find that the all the sites are selected as storage warehouses for major relief materials. In the second stage, materials will be distributed in the main disaster relief areas, i.e. from the initial center (Airport) to major disaster relief sites planned according to the material transportation. The transportation of relief materials mainly involves major site M 4 and M 5 , which respectively account for 32.91% and 32.52% of the total demands. Both of them are responsible for providing major relief supplies. The third stage is aimed at offering the material distribution services of the secondary disaster sites. On the premise of meeting the demands of each node, the distribution path of TS-MILO model almost traverses all feasible paths. Through careful analysis, we can find that although the distribution path can guarantee the supply of basic materials via the TS-MILO model, there are still some problems in the specific service process. For example, the long-distance transportation in the route planning will increase the fuel cost. The time costs caused by circuitous transportation in route planning will increase the suffering of the rescued. The unreasonable use of major disaster relief sites will lead to the increase of costs of subsequent re transportation. Once there is uncertainty in the actual rescue process, for example, the fluctuation of demands in the target area will increase. These uncertainties will make the stability and sustainability of TS-MILO model unable to be guaranteed, which makes the relief material logistics service face some challenges and difficulties. As a result, in the rescue process, we must plan reasonably, explore more optimized improvement strategies, and seek more robust path planning scheme. According to the above analysis, we must optimize the distribution path.

C. BTS-MIRO Model
In the BTS-MIRO model, the effect of j,k on the total costs is constantly changing. The results of the BTS-MIRO model are shown in Table 6 below. The total costs of HEL increase with the increasing of SP. According to the Theorem 1, when j,k = 0, the BTS-MIRO model corresponds to the TS-MILO model and the total costs are 8.5274 E+05 CNY.
In BTS-MIRO model, the rate rises from 8.5274 E+05 to 9.3323 E+05, with an increase of 9.44%.  Fig. 6, the BTS-MIRO model is compared with the TS-MILO model: In the first two stages, three major sites M 3 , M 4 , M 5 bear more than 81.31% of the material supply of secondary disaster relief services. In the third stage of secondary disaster relief logistics planning, the main disaster relief site M 3 undertakes a larger distribution service than before, with the transshipment proportion increasing from 11.8% to 28.57%. The transshipment proportion involving M 4 , M 5 declined, but they still account for 27.51% and 29.87% of the main relief sites respectively. Although the distribution line is centralized to M 3 , and the route planning is relatively reasonable, there are still some circuitous transportation problems in other lines, so there is still room to optimize the BTS-MIRO model.

D. PTS-MIRO MODEL
In the PTS-MIRO model, the total costs vary with SP. The results of PTS-MIRO model are shown in Table 7. The total costs of emergency relief logistics increase with the increase of SP. According to the Theorem 2, when j,k = 0, the PTS-MIRO model is equivalent to the TS-MILO model.  In BTS-MIRO model, the rate rises from 8.5274 E+05 to 1.0205 E+06CNY, with an increase of 19.67%. It can be seen that compared with BTS-MILO model, PTS-MILO model requires higher costs.
As is shown in Fig. 7, the PTS-MIRO model is compared with the BTS-MIRO and TS-MILO models: In the initial route planning of first two stages, M 3 , M 4 , M 5 are still the main material transfer centers. From the comparison of relative location sites, it can be found that the transfer of transport transfer center makes the route planning more reasonable and reduces the circuitous transportation in secondary stage. Compared with BTS-MIRO model and TS-MILO model, PTS-MIRO model reduces the proportion of long-distance line transportation and increases the proportion of shortdistance line transportation in the third stage of secondary disaster relief path planning, which can provide specific reference for quickly reaching destination, which has a positive role in promoting path optimization.

E. ETS-MIRO MODEL
The results of the ETS-MIRO model are shown in Table 8. Similarly, the total costs increase with the increase of SP of  the model. According to Theorem 3, when j,k = 0, ETS-MIRO model is equivalent to TS-MILO model, and the cost is 8.5274 E+05CNY. In ETS-MIRO model, the rate rises from 8.5274 E+05 to 8.8568 E+05CNY, only with an increase of 3.87%. It can be seen that compared with BTS-MILO and PTS-MILO model, ETS-MIRO shows higher robustness.
As is shown in Fig. 8, the ETS-MIRO model is compared with the PTS-MIRO, BTS-MIRO and TS-MILO models: In the initial planning of the first two stages, the proportion of transport transfer is relatively evenly distributed in the major transport centers. The transit capacity and load pressure of each major disaster relief point are relatively balanced. The sites have penetrated into the hinterland of the disaster area and is closer to the affected people. In consequence, increasing the storage ratio of M 2 makes the path planning more reasonable. In the third stage of planning, the ratio of PTS-MIRO, BTS-MIRO and TS-MILO models is relatively higher. The proportion of long-distance transportation is further reduced, while the proportion of short-distance transportation is increased in the ETS-MIRO model.
Especially after making full use of this node which goes deep into the disaster area, the performance is more obvious. In comparison, the proportion of services in each path tends to be short path. The costs of short distance route are fewer, and its proportion in the supply of materials shows an upward trend. Due to the higher efficiency of vehicle mileage and more accurate and quick distribution route, TS-MIRO has shown better optimization performance.

VI. MODEL SENSITIVITY ANALYSIS
This section compares and analyses the performance of each model, including operational efficiency, uncertainty, and demand fluctuation.

A. OPERATING PERFORMANCE COMPARION
This section compares the operation efficiency of each model by observing the operation time of the models. In order to ensure the effectiveness of the results, the following settings are made: running in the same computer environment, setting the SP as the only variable. We use MATLAB as the platform. The solvers GUROBI and CPLEX are used to solve the above models respectively, and the results are shown in Table 9. In the total costs comparison of the model, the costs in the calculation results of the model is not affected by the algorithm. Gap%<0.1 is found in the error comparison of the calculation, which indicates that the model is effective without significant difference in different solvers. From com-parison of the calculation performance of the model, it can be seen that those two solvers are used to solve the models. GUROBI can save at least 38% of the time consumption compared with VOLUME 8, 2020 the solvers CPLEX. Obviously, the algorithm of GUROBI is better than that based on CPLEX. Fig. 9 show that, the TS-MILO model has the highest running efficiency and the overall operation time is much shorter than that of MIRO model (≤2600 ms). Detailed comparison of the four models shows that, when SP ≤ 3, BTS-MIRO>ETS-MIRO>PTS-MIRO, which means that the BTS-MIRO model is the most efficient. When 3 ≤ SP ≤ 8, PTS-MIRO>ETS-MIRO>BTS-MIRO, which means that the PTS-MIRO model is the most efficient. When 9 ≤ SP, ETS-MIRO>PTS-MIRO>BTS-MIRO, which means that the ETS-MIRO model is the most efficient with a large number of uncertain parameters. Due to small scale of this HEL problem, there is little difference in time. However, when the constraints and variables in the model increase to tens of thousands or even tens of millions, the operational efficiency of the model solution will be significantly different. There-fore, when solving practical problems, we can build an ap-propriate model according to real data.

B. IMPACT OF DEMAND VOLATILITY AND SAFETY PARAMETERS
In this section, we analyze the impact of demand volatility on the total cost. The control variable method is used to analyze the effect of the fluctuation of random parameters on the total cost, where ( j , j , j = 3; k , k , k = 7). The calculation results are shown in Figure 10. Although TS-MIRO model will pay a certain cost, that is, the total cost is higher than TS-MILO model, it can still give path planning scheme even if in uncertain situation. The growth rate of costs is slightly different. The costs of PTS-MIRO model increase sharply, while that of ETS-MIRO model increases slowly. In the PTS-MIRO model, when the demand volatility increases from 0 to 30%, the total costs increase from 8.5512 E+05CNY to 9.6885 E+05CNY, with an increase of 13.6%. In the ETS-MIRO model, when the demand volatility increases from 0.00% to 30%, the total costs increase from 8.5512 E+05CNY to 8.8907 E+05CNY, with an increase

C. SERVICE LEVEL AND ITS RESPONSIVENESS
The performance is analyzed through Level of Service (SL). Due to high requirement for timeliness in HEL, this section  compares the SL of models through time differences, and analyses the advantages and disadvantages of the different models. The SL is calculated by following formula. Com-puter simulation results under different parameters are shown in Fig. 12 and 13.
where, I , J represents the number of arcs in the model. Fig. 12 shows the influence of volatility on SL, under fixed SP ( j , j , j = 3; k , k , k = 7). Through analysis of the influence of the volatility on the model, the following three conclusions can be drawn. The SL of TS-MILO model is not affected by fluctuation, and of course, the problem under uncertain conditions is unable to solve. Due to the determined data, the SL is also the highest (SL=92.87% of the second stage and 86.97% of the third stage). The SL of three TS-MIRO models shows downward trend with the increase of volatility. The greater the amplitude of volatility is, the lower SL is. The impact of demand volatility on SL in the second stage is greater than that in the third stage. The reason is that the uncontrollability of rescue work aggravates the uncertainty in the process of HEL. For the major disaster areas, the uncertainty of secondary disasters should be considered in addition to the impact of uncertainty fluctu-ation in their own regions. Due to the lack of information and low efficiency of resource turnover for secondary disaster areas, the cumulative uncertainty has a more significant influence on the second stage (SL downward). Volatility has a higher influence on SL in the second stage than that on the third stage. The reason is that in the planning of major disaster relief material dispatch, besides the impact of uncer-tainty itself, the influence of uncertainty on secondary disa-ster needs considering. It can be clearly seen from the com-parison that the ETS-MIRO model has strong robustness. The BTS-MIRO model is centered and the PTS-MIRO model is most affected by volatility, which means that PTS-MIRO model has the weakest ability to resist uncertainty, while ETS-MIRO model performs better.
Fortunately, Fig. 13 illustrates the impact of SP on the SL under the condition of fixed volatility (ε = 0.15). We can see that the SL tends to increase in SP. SP has a great effect on the improvement of service quality, and it performs well both in the second and the third stage. To a certain extent, this makes up for the robustness cost caused by uncertainty and also reduces the loss of SL owing to volatility. Careful comparison shows that ETS-MIRO model has strong robustness. When SP increases from 1 to 14, SL will increase from 88.53% to 96.44% in the second stage and from 78.74% to 89.61% in the third stage. In the process of disaster relief, managers must pay attention to the rapid responsiveness of rescue. Considering uncertainties, although each TS-MIRO model can give path planning plans, the performance and application scope of each plan are also different. HEL decision makers must review the situation and make the most reasonable plan according to local conditions.

VII. CONCLUDING REMARKS
Humanitarian emergency logistics (HEL) management is a major decision-making related to social economy and people's life safety, which is highly valued by the national government. Most of the previous studies focused on major disasters, but the impact of secondary disasters was easily ignored. On the basis of fully analyzing the internal relationship between primary and secondary disasters, this paper establishes a three-stage mixed integer linear optimization (TS-MILP) model based on inventory, vehicle load and time constraints. Under the condition that all parameters are determined, this paper compiles the algorithm based on GUROBI, and gives the basic feasible path planning scheme. In the previous research, it was often carried out under certain circumstances, while the actual rescue process is uncertain. Therefore, the management of HEL under uncertain conditions has higher research value. In order to VOLUME 8, 2020 prevent the dissimilation effect of uncertainty, robust optimization is introduced. By constructing the corresponding uncertain set, the TS-MILO model is transformed into a three-stage mixed integer robust optimization (TS-MIRO) model. Finally, the practicability is verified by experiments.
The main contributions of this article are: First, the traditional HEL problem in major and secondary disasters is subdivided, and the single-stage problem is subdivided into three-stage problem. Secondly, the robust optimization theory is introduced, and the deterministic optimization model is further extended to the robust optimization model with uncertainty in mind. Finally, experiments are carried out with real data, and a feasible route planning scheme is given based on the model calculation results.
Through the study of this paper, we can draw the following conclusions. On account of the idealization of parameter data, the three-stage mixed integer linear optimization model can give the lowest total costs and the highest service level. However, in reality, such an ideal data set is extremely difficult to obtain, even unavailable, so its practical feasibility is not high. On the whole, the three-stage mixed integer robust optimization model will pay a certain robust cost, such as the increase of costs and the increase of complexity of calculation, due to the influence of stochastic demand fluctuation. But they can solve uncertain humanitarian logistics problems and maintain robustness. From the details, the service level of the improved robust model has been significantly improved. The ellipsoid set robust optimization model performs best, while the polyhedron set robust optimization model has the weakest ability to resist environmental changes. The robust optimization model proposed in this paper can better solve the uncertainty of relief work.
The proposed model can not only be used to solve the problem of HEL in the earthquake areas, but also has great expansibility. It can also be applied to other problems, such as supply chain management with secondary replenishment, vehicle routing planning with secondary transshipment, home medical service with personnel transfer, etc. These problems have similar characteristics, that is, in addition to the main demand or supply, there is also a demand for secondary. On the technical level, the programming algorithm used in this paper still has plenty of room for improvement. Through the application of new technology, the proposed model will be more accurate and more efficient. In the future, our research will explore the application mode of science and technology. It can be predicted that the informatization and intelligence of emergency rescue management is the future direction.
JIANHUI DU is currently the Director of Guangdong Maritime Chamber of Commerce, Guangdong Province. He has unique views on supply chain management, coordinated sustainable development of economic and social resources, and mathematical modeling and optimization. He has strong scientific research and innovation ability. His main research interests include supply chain management, game theory and optimization, and robust optimization. VOLUME 8, 2020 YING JI was a Postdoctoral Candidate with the Asia-Pacific Logistics Research Institute, National University of Singapore, a Visiting Scholar with the Business School of National University of Singapore, and a Postdoctoral Candidate with the Management Science and Engineering Mobile Station, Harbin University of Technology. She was also an Associate Professor with the Institute of Basic and Interdisciplinary Sciences, Harbin University of Technology. She is currently a Professor and a Doctoral Supervisor. Her main research interests include supply chain finance, interest rate optimization, and optimization theory and application. Her current research interests include data analysis, portfolio optimization, and decision science.
DEQIANG QU has been engaged in the research of physical geography and computer software for a long time, and has made outstanding contributions in the field of GIS management. He also has in-depth research on the application of the Internet of Things technology and an understanding of frontier knowledge. In addition, he has extensive research on industrial manufacturing and deep insight in global urban manufacturing.
XIAOQING WU received the bachelor's degree in management from Guangzhou University, Guangzhou, China, in 2020. She is currently pursuing the master's degree in social work with Hohai University, Nanjing, China. Her current courses include principles of management, management psychology, urban politics, social work theory, and advanced social work practice. Her main research interests include social management, social work, humanitarian relief, and migration.
DAN YANG received the bachelor's degree (Hons.) from the Guangzhou Maritime Institute and the M.Sc. degree (Hons.) from Cranfield University, U.K. He has a deep understanding of mathematical modeling and is good at Machine Learning and Big Data analysis. His main research interest inlcudes supply chain network design and optimization.