Bi-objective robust optimisation on relief collaborative distribution considering secondary disasters

ABSTRACT Developing an effective emergency collaborative distribution system is critical to improve on-time response performance, especially considering secondary disasters. To address the challenge, this paper investigates a bi-objective robust optimisation model on relief collaborative distribution among three echelons of authorities, the province, the municipality, and the county, which aims to minimise the total travel time and the total humanitarian logistics cost simultaneously. The optimal location of relief supply facilities and the relief distribution schemes will be determined by the optimisation model, which considers uncertain demand and travel time. Moreover, two robust optimisation methods are utilised to deduce the robust counterparts of the proposed model. An epsilon-constraint-based approach is used to solve the bi-objective optimisation model. A real-world case study based on an earthquake and aftershocks with different magnitudes in Yunnan Province is provided. The results show that incorporating secondary disaster scenarios contributes to reducing the total travel time and cost. For making full use of emergency resources and preventing situations from worsening, the centralised decision scheme is more effective than the decentralised one. The uncertainty of demand of primary disaster relief has a higher impact on the optimal solution than that of travel time.


Introduction
Disasters incur catastrophes for many areas and cause serious damage to property and even human life.In the first half of 2022, 187 disasters affected over 50 million people and caused approximately $41 million in economic damage (Cred Crunch EM-DAT 2022).Some disasters usually induce relevant secondary disasters (SDs).For example, earthquakes may cause mountain collapses and landslides.Swiss Re Institute announced that in 2021, $111 billion in insured losses was caused by natural disasters worldwide, and 73% of these losses was incurred by SDs.Therefore, it is very important to design an effective humanitarian logistics (HL) system to reduce disaster losses as much as possible.
To design and operate an effective HL system, decision-makers can consider the corresponding factors from a disaster cycle perspective, which generally includes four phases, i.e. mitigation, preparedness, response, and recovery (Gaisie, Han, and Kim 2021).Moreover, it is also vital to transport relief supplies (RSs) to specific facilities throughout the above four phases, CONTACT Shuanglin Li lishuanglin@csu.edu.cnSchool of Traffic & Transportation Engineering, Central South University, Changsha, People's Republic of China Supplemental data for this article can be accessed online at https://doi.org/10.1080/00207543.2023.2217306.which can be regarded as relief distribution (RD) (Ye, Jiao, and Yan 2020).Especially in the response phase, effective and efficient RD will significantly improve emergency management (Sabbaghtorkan, Batta, and He 2020).
However, RD will face challenges due to the disruption of communication and transportation facilities after disasters (Kovacs and Moshtari 2019).The hindered information transfer and the uncertain demand will cause the ripple effect in terms of RD (Ivanov 2020).To address the uncertainties, many optimisation methods have been exploited, such as fuzzy programming, stochastic programming (SP), and robust optimisation (RO).Compared with fuzzy programming and SP, RO uses uncertainty sets to handle uncertain parameters and does not require much historical data.RO can also take decision-makers' risk preferences into consideration (Najafi, Eshghi, and Dullaert 2013).
Furthermore, SDs occurring in a short time will generate extra demand and cause confusion in response decisions related to the primary disaster (PD), which is in line with the phenomenon of the ripple effect in the humanitarian supply chain (Sawik 2022).Dolgui and Ivanov (2021) determined that decision-makers can adopt actions with the interplay of flexibility and redundancy to ensure resilience and ripple effect mitigation.Thus, decisions for PD should be adjusted to prevent a lack of supplies.
Moreover, a single local organisation cannot respond to disasters effectively (Baharmand, Comes, and Lauras 2020).The integration of proactive and reactive strategies from different organisations at different levels can balance supply chain resilience and efficiency (Aldrighetti et al. 2021).Emergency rescue collaboration generally covers two different types: vertical collaboration and horizontal collaboration.The former takes place among governmental departments at different levels, e.g.country-region-province and city-town, while the latter exists among different entities or nonprofit organisations at the same level (Rodríguez-Espíndola, Albores, and Brewster 2018).Kondo and Lizarralde (2021) addressed an optimisation model for the coordination of relief supply facilities (RSFs) at different levels.
Motivated by the above analysis, this paper investigates a robust optimisation problem on relief collaborative distribution considering SDs with uncertain conditions.New settings are adopted to improve the adaptation and viability of the supply chain, especially the humanitarian supply chain in our area of study (Ivanov 2021).A bi-objective optimisation model is presented, which determines the optimal RD scheme from multiple different level facilities to multiple affected areas to achieve the following two conflicting objective functions, i.e. the total HL cost and the total travel time.Robust optimisation methods are adopted to address uncertain demand and travel time.We refer to the above bi-objective robust optimisation model as the BOROM.Furthermore, an epsilon-constraint method is used to determine the location of the RSFs and the RD scheme.Some numerical experiments derived from earthquakes in Yunnan Province, China, illustrate the robustness and advantages of our proposed novel model.A sensitivity analysis of key parameters in the proposed model is carried out, and relevant management insights about relief collaborative distribution are revealed.
The remaining sections of this paper are organised as follows.In Section 2, a review of the related literature is provided.In Section 3, the relief collaborative distribution problem is formulated as a BOROM.Section 4 elaborates on RO and the epsilon-constraint method.A case study, along with several sensitivity analyses of the key parameters, is provided in Section 5. Section 6 draws some managerial insights.Finally, we outline conclusions and directions for future research in Section 7.

Literature review
A considerable number of papers have been published on the optimisation problem of disaster relief logistics (Cheng, Feng, et al. 2021;Tang et al. 2022a).This paper focuses on the relief collaborative distribution problem considering uncertainties generated by PD and SDs and multiple objectives.Therefore, we divide the related papers into three aspects: (i) works related to RD, (ii) RO methods to address uncertainties, and (iii) HL optimisation with multiple objectives.

Relief distribution
To design a resilient humanitarian logistic supply chain, some researchers have proposed simulation models based on different data inputs, algorithms, and scenarios (Singh et al. 2021;Stewart and Ivanov 2022), while mathematical modelling is usually adopted to solve specific optimisation problems.In the disaster cycle response phase, activities mainly include RD and road network repair (Li and Teo 2019;Li, Ma, and Teo 2020;Ransikarbum and Mason 2016).A timely and effective RD operator can not only reduce property loss but also protect people's lives (Kundu, Sheu, and Kuo 2022).Seraji et al. (2021) located emergency shelters and optimised the vehicle routing of RD.Tang et al. (2022b) took replenishment into consideration and proposed a multi-period model to improve COVID-19 vaccination delivery efficiency.
To provide more efficient humanitarian aid, a collaborative response is valuable (Paciarotti, Piotrowicz, and Fenton 2021).Coordination means interactions among many decentralised agencies, such as authorities, nongovernmental organisations and individual volunteers (Maon, Lindgreen, and Vanhamme 2009).There are many different kinds of decision patterns among several participants, such as cooperation of governments and nongovernment organisations (Sheu and Pan 2014), multiple organisations at different governmental levels (Rodríguez-Espíndola, Albores, and Brewster 2018), vertical and horizontal collaboration of hospitals, emergency shelters and region units (Kamyabniya et al. 2021) or even stations that provide different services (Jiang and Ouyang 2021).
In the real world, SDs can influence the response operations of the first shock.Ahmadi et al. (2020) indicated that considering SDs in relief delivery contributed to higher demand satisfaction than considering PD only.However, there are few papers considering the effect of SDs.Zhang, Li, and Liu (2012) used probability to describe the occurrence of different types of SDs.Conditional probability and probability trees were adopted to address the correlation between PD and SDs (Zhang et al. 2019;Hu et al. 2019).

RO for HL uncertainties
Unpredictable disasters make it difficult for decisionmakers to obtain exact information in the disaster response phase (Yánez-Sandivari, Cortés, and Rey 2021).It is critical to propose a RD scheme considering social factors, timeliness, and cost (Zhang et al. 2020).Scholars have tried to address these uncertain factors in different ways, such as updating information in a timely manner (Ismail 2021;Fang et al. 2021) or using fuzzy methods (Shaw, Das, and Roy 2022), SP (Li, Yu, and Zhang 2021) or even social media technology (Sarma, Das, and Bera 2020).
Because the probability of SP will affect the optimal decision, a RO model was used to describe uncertain parameters in many studies (Sanci and Daskin 2021).Robust optimisation can be classified as robust and distributionally robust optimisation methods by the construction of uncertainty sets (Lu and Shen 2021).Box uncertainty sets are the general methods to address uncertainty parameters (Sun, Wang, and Xue 2021).Dalal and Üster (2021) captured the evacuee, stay-back population, relief demand and disaster duration uncertainty by box uncertainty sets to obtain the optimal time of relief delivery.For distributionally robust optimisation, Wang and Chen (2020) adopted a moment-based ambiguous set to address uncertain distributions of blood demand and optimised the blood supply network design.For easier computation, Wang and Chen (2020) reformulated the distributionally robust model with box and polyhedral ambiguity sets.Furthermore, the combination of machine learning and robust optimisation can also be used to obtain better solutions (Nabavi et al. 2022).
For multi-stage problems, a multi-stage robust framework and dynamic robust programming can also be used to address uncertainties (Yang, Liu, and Yang 2021;Cheng, Adulyasak, and Rousseau 2021;Wang, Yang, and Yang 2022).

HL optimisation with multiple objectives
Regarding the optimisation problem of HL, the objective functions can be divided into three categories: efficiency, effectiveness, and equity (Gralla, Goentzel, and Fine 2014).In terms of equity, there are two policies for decisions, egalitarian and utilitarian.The objective of the egalitarian objective function is to optimise the worst-case performance, such as maximising the minimum total cost (Zhang et al. 2023) or minimising the maximum number of locations (Jenkins, Lunday, and Robbins 2020).The utilitarian objective aims to optimise the entire disaster chain, minimising the total cost (Aliakbari et al. 2022) or maximising the profits obtained from the drone nodes and links (Zhang et al. 2021).
In the decision-making process, HL optimisation always has multiple objectives to meet different kinds of requirements at the same time (Jenkins, Lunday, and Robbins 2020).Mosallanezhad et al. (2021) compared different kinds of personal protection equipment and obtained a distribution scheme to minimise the total cost, as well as the unmet demand.Wang et al. (2022) addressed the RD optimisation problem, aiming to simultaneously minimise the expected total cost and maximise the supply shortage rate of a health care coalition.
Some scholars have adopted the fuzzy method and epsilon-constraint method to obtain the equivalent single-objective model (Akbarpour, Torabi, and Ghavamifar 2020;Maharjan and Hanaoka 2020).For large-scale multi-objective optimisation problems, NSGA-II is regarded as an effective method (Ransikarbum and Mason 2022).There are also some other metaheuristic approaches, such as multi-objective vibration damping optimisation (Seraji et al. 2021).

Comparisons with the state of the art
As shown Table 1, the main related studies can be divided into two groups based on the five criteria, disaster response operations and modelling.The literature on disaster response operations is further categorised based on whether SDs and collaboration are considered.For the second group, there are the following three criteria: the number of objective functions, the kinds of uncertainties, and the optimisation methods.
The following conclusions are obtained based on the above review of related papers.
(i) Few studies have considered the impact of SDs or cooperation on decision-making.(ii) Most studies cope with uncertainty with multi-stage SP or RO.(iii) There is no study that considers SDs, collaborative responses, and multiple objectives in RD plans at the same time.
To fill the above research gap, we address a BOROM on the relief collaborative distribution problem considering the PD and SDs with the uncertainty of demand and travel time, which aims to minimise the total travel time and the total HL cost simultaneously.
Our main contributions are fourfold.(i) A BOROM on relief collaborative distribution with SDs is proposed, which minimises the total HL cost and total travel time simultaneously.(ii) The robust counterparts are deduced to address uncertain relief demand and travel time.The epsilon-constraint method is utilised to obtain the Pareto front.(iii) Numerical experiments derived from

Real-life RD system description
This paper was motivated by the current emergency distribution system in China, which is a large country with frequent natural disasters.To establish a robust emergency management system that can respond to disasters rapidly and effectively, the State Council promulgated the 14th five-year plan of the national emergency system.Based on the centralised political system with several echelon government authorities, the plan suggested that Chinese emergency management adopt the form of territorial management and echelon coordination.At present, central-level relief material reserves have been established in 20 cities to meet the demand of major catastrophes.Supply facilities at the province, municipality, and county levels are only able to meet the RS support needs for initiating level II emergency response.
There are some factors disrupting RD plans in HL: (i) Uncertain demand and travel time: The types of disasters, time and location of occurrence make the demand of RSs uncertain.Some catastrophes will damage the transportation network and delay delivery.(ii) Conflict objectives: Due to the limitations of RSs and budgets, decision-makers should consider the cost of the optimal solution.However, timeliness also plays a huge role in RD.
(iii) SDs: Subsequent shocks occur with different probabilities and lead to uncertain disaster response duration.The RD will be disrupted due to uncertain road network conditions and RS demand (Li et al. 2023).(iv) Collaboration modes: In China, once the PD occurs or even when it is approaching, the local government will complete the RD alone first, and then the government at higher levels will deliver the cannibalised RSs to meet all requirements.The process of seeking for help will cost extra time and money.
Therefore, it is a core mission of decision-makers to transport relief from different levels of supply facilities to affected areas while considering SDs.

RD network illustration
A relief collaborative distribution network is illustrated in Figure 1.There are two types of facilities in this network: relief supply facilities (RSFs) and demand points (DPs).
(i) RSFs: These facilities can store and supply various RSs.Based on administrative institutions, RSFs belong to three echelons of government: province, municipality, and county.The lower the echelon is, the lower its capacity and fixed cost are.Only the local government is able to manage RSF inventory.
There are administrative boundaries in RSFs.(ii) DP: Once a PD occurred, all of these places were attacked by different types of SDs, causing the relief requirements of the DPs to vary.
The left side of Figure 1 shows the existing scheme in China, which is called a decentralised decision scheme.The local government in the affected area will respond to the disaster independently at first.Only when the local RSs nearly run out will the government search for additional help through successively higher levels of government.Then, the RSFs at higher levels will transport relief to the RSFs at lower levels.All the RSs are distributed from local RSFs.
In this paper, we propose a centralised decision scheme, which can be seen on the right side of Figure 1.The provincial authority will command the distribution.The individual RSFs at different levels can serve DPs at the same time.In the collaborative RD, RSFs and DPs are interconnected.Thus, the network can be called an interwind supply network considering viability (Ivanov and Dolgui 2020).Because the disruption duration is uncertain, relief transportation can be divided into two parts: (i) meeting the post-disaster demand immediately and (ii) keeping SDs, which can also be called safety stock.
Based on the emergency management system, a BOROM is proposed to minimise the total travel time and the total HL cost, which allows efficient and effective response operations in the PD and SD response phases (Aldrighetti, Battini, and Ivanov 2023).Once the PD occurs, the provincial authority decision-makers will obtain the RSF location and RS allocation plan under demand and travel time uncertainties considering a trade-off between efficiency and resilience.

Mathematical formulation
Assumptions, definitions of sets, parameters, and decision variables are given as follows.

S jm
Amount of safety stock of RS m at DP j for SDs, which is held for SPs in advance.
Based on the abovementioned descriptions and definitions, a collaborative RD nominal model is proposed. min Equation ( 1) minimises the maximum total travel time in the PD and SDs, which is calculated by summing the weighted required travel time from the RSFs to the DPs of each trip.Equation ( 2) minimises the total HL cost, including the fixed cost for opening RSFs at different levels, the distribution cost of RSs from RSFs to DPs, and the penalty cost of unsatisfied demand in the PD and SDs.Constraint (3) requires that one candidate RSF can be established at only one level.Constraints (4) and ( 5) indicate that the number of RSs transported from all RSFs and the unmet demand should be greater than or equal to the requirement at each DP.Constraints ( 6) and ( 7) ensure that the number of RSs delivered to each DP should satisfy the lowest demand of the DP.In particular, SDs lead to extra demand in the PD response phase; thus, we introduce S jm to describe the safety stock for SDs.Constraint (8) establishes that the number of RSs delivered from one RSF cannot exceed its capacity.Constraint (9) enforces the binary restriction on the corresponding decision variable.Constraint (10) defines the nonnegative and integer constraints of decision variables.

Solution approach
In this paper, we investigate a BOROM with uncertainties, which includes two difficulties.One is that the relief requirement and required travel time between each pair of nodes are both uncertain.The other is that the optimal solution should be performed on more than one objective function.Therefore, to solve this model, we address the uncertainties by RO and then develop the equivalent single-objective model by the epsilon-constraint method.

Formulation of the robust model
Because the damage caused by natural disasters and the road network conditions in the affected areas are unpredictable, the demand and required travel time for delivering RSs are uncertain parameters, which may cause ripple effects (Sawik 2023).However, SP needs sufficient historical data to estimate the probability distribution (Tirkolaee et al. 2020).Moreover, the more scenarios that are taken, the more precise the results will be, and too many scenarios for incorporating uncertainty will increase the complexity of computation (Gao and Cao 2020).Compared to SP, RO introduces uncertainty sets for modelling uncertainties, which provides a series of solutions depending on the risk preferences of decision-makers.
There is an uncertain parameter, travel time tkij , in Objective Function (1).Based on Bertsimas and Sim (2004), we introduce an auxiliary variable f to convert Function (1) to Equation (11) and Constraint (12) equivalently. min Furthermore, Constraint (12) can be rewritten as: In Constraint (13), tkij (j ∈ J) takes values in [ tkij − tkij , tkij − tkij ], where tkij denotes the nominal required travel time and tkij represents the maximum deviation of tkij .In practice, not all required travel times between RSFs and DPs are uncertain.Parameter ki ∈ [0, |J|] denotes the level of conservatism.If ki is an integer, ki of all required travel times will fluctuate within the value range, and the maximum deviation value will be taken.If it is not an integer, ki of all required travel times will fluctuate to the maximum deviation value, and only one required travel time is ( ki − ki ) tkij , while the solution remains feasible.The larger ki is, the less risk decision-makers will take, and the more conservative the solution is.However, if more uncertain coefficients deviate from their nominal values within the predefined range, the probability that the robust solution remains feasible satisfies the following: where Y α * kijm and Y b * kijmω are the optimal solutions of the model and (θ) is the cumulative distribution function of a standard normal distribution.
We introduce a protection function β(y, ki ) for Constraint (13) as follows: Then, Constraint ( 13) is equivalent to: β(y, ki ) can be transformed into the following linear programming problem: The dual problem of the linear programming problem is: in which p kij and z ki are the dual variables.

Transformation of the bi-objective model
It is impossible to find a solution to optimise two conflicting objectives at the same time.The weighted sum method is usually used to address multi-objective optimisation problems.Although the method is very simple, the weight coefficient has a great influence on the solution quality (Mohamadi, Yaghoubi, and Pishvaee 2019).The epsilon-constraint method controls the effective solutions by relaxing each objective as a constraint with different thresholds (Mavrotas 2009).Therefore, we choose the epsilon-constraint method to solve the problem with two conflicting objective functions.The specific form is as follows.

Case study
To investigate the practicality of BOROM and the algorithm, we conduct numerical tests based on the RD problem in the Yunnan earthquake.We used MAT-LAB software coupled with CPLEX 12.9 as the development environment.The tests ran on a laptop with an Intel(R)Core (TM) i5-1135G7 2.40 GHz CPU and 16 GB RAM.

Case study description
Complex landforms and climates cause various natural disasters, such as earthquakes, droughts, and floods, which frequently occur in Yunnan Province, China.We take earthquakes and aftershocks with different magnitudes in this location as the research background.Figure 2 denotes the map of RSFs and DPs.
There are 16 cities and 112 counties, where RSFs are divided into the province, municipality, and county levels.Suppose aftershocks of different magnitudes represent the SDs ω 1 , ω 2 , ω 3 , ω 4 , ω 5 .The unit penalty cost of food, water, and medicine is 12 times the distribution cost, while the unit penalty cost of bedding is the same as the distribution cost (Wei, Hou, and Guo 2019;Cai, Zhu, and Wang 2017).RSFs at different levels have different fixed costs and different capacities.Some parameters of the case study can be found in the supplementary materials.

Simulation results in the centralised decision scheme
To discuss the influences of uncertain demand (including PD demand and SD demand) and travel time on the main objective function, suppose that parameter ∈ [0, |J|] denotes the conservation level of the uncertain parameters.Thus, α jm = Table 2 gives the optimal RD plan in the centralised decision scheme ( ε 2 = 1E + 10, = 5, γ = 10% ).The last column shows when the RSF will open, including the PD and SD response phases ( = {ω 1 , ω 2 , ω 3 , ω 4 , ω 5 }.)The minimum total travel time is 1.0222E + 06.As shown in Table 2, all RSFs in the DPs are chosen, which is consistent with the real-world RD.For RSs of large size, the demands of DPs must be supplied by other RSFs, such as Yongsheng, Fumin, Jianshui, and Nanjian counties.Some RSFs only open for some SDs to supply the bedclothes, such as Lijiang Municipality, Shiping, Midu and Weishan counties.In addition, regardless of the type of SD, the RSF response strategy will not change.Moreover, some RSFs have strong synergies, such as Nujiang-Fugong-Gongshan and Gengma-Cangyuan-Shuangjiang.In most cases, these RSFs will serve the same DPs.
Table 3 shows the safety stock of each RS for SDs, which are transported to demand points after a PD.For example, J 1 will prepare 1757 food resources, 21,957 water resources, 983 medicine resources and 37,831 bedclothes resources for SDs.It is obvious that all demand points need safety stock but not in all commodities.J 4 and J 8 keep only water, medicine, and bedclothes.

Robustness analysis
The optimal solution is investigated for the BOROM with different levels of uncertainty.The required travel time ki is generated from the interval [0, |J|].PD demands and SD demands are uncertain, and their budget of uncertainty takes values within the interval [0.1].To analyse the impact of uncertainty, suppose that uncertain parameters deviate 5%, 10%, and 15% from their nominal values, which means data variability γ = 5%, 10%, 15%.Only one uncertainty is analysed at a time, and the other parameters take the nominal value.The incremental variation of the main objective function is calculated by subtracting the objective function value of the robust model from that of the nominal model.( ω 1 , ω 2 , ω 3 , ω 4 , ω 5 represent different subsequent shocks).
Since the total humanitarian cost is converted into a budget constraint in Section 4.2, we propose a novel criterion, rescue quality (RQ), to estimate the emergency management quality of the cost.The total humanitarian cost includes two parts, i.e. the rescue cost and the penalty cost.The former is actually used in the disaster relief process, while the latter occurs due to unmet needs.RQ is the proportion of the total cost that is associated with RSF location and RD.The higher the proportion is, the more demand will be satisfied.RQ = (fixed cost for opening RSFs + distribution cost in a PD and SDs)/total humanitarian cost.
Figures 3(a) and 4(a) shows that the total travel time and RQ increase with the increase in the budget of PD demand uncertainty ( ε 2 = 1E + 10 ).The higher the data variability is, the faster the rate of growth for the total travel time is, as well as RQ.In Figure 3(a), the increase in the total travel time exhibits linear growth in all data variability.However, Figure 4(a) shows that the increase in RQ increases all data variability.Now, we address the effects of SD demand uncertainty on total travel time and RQ. Figure 3(b) shows the same trend ( ε 2 = 1E + 10 ) as Figure 3(a).The greater the demand uncertainty is, the greater the impact is on the value of the total travel time.Figure 3(b) shows that the impact generated by the uncertain demand changes in the PD for the main objective function is much greater than that of the SDs under the same uncertainty level.The budget of SD demand uncertainty can lead to an increase in RQ but shows no relationship with the growth of increases (see Figure 4(b)).
Figures 3(c) and 4(c) indicate the effects of PD and SD demand uncertainty on the total travel time and RQ ( ε 2 = 1E + 10 ).The growth of total travel time increases monotonically with the budget of uncertainty and the data variability at the same time (see Figure 3(c)).Figure 4(c) shows that the more uncertain  the parameters are, the higher the RQ of the distribution scheme is, which indicates the robustness of the proposed model.
Figure 3(d) depicts the change in the total travel time under different budgets of required travel time ( ε 2 = 1E + 10 ).The uncertainty of the travel time affects the growth of the total travel time only within a specific limit.If the budget of required travel time is equal to a certain value (e.g. 2 in this study), the increase in the value of the main objective will remain constant.In addition, the maximum budget of required travel time is much higher than the specific value.The uncertainty of travel time has no influence on the increase in RQ, as illustrated in Figure 4(d).The sensitivity of the robust model to variation in travel time can be found in the supplementary materials.
Comparison between robust model and other different models can be found in supplementary materials.

Safety stock analysis
Because the safety stock is prepared for SDs in the PD response phase, we utilise the BOROM with different uncertain budgets and data variabilities to find the rules of safety stock changes ( ε 2 = 2E + 10 ).The results are shown as follows: As shown in Figure 5(a-1) and (a-2), the ratio of water and bedclothes safety stock to both PD and SD nominal demand maintains the same trend.With the increase in data variability, the ratio of food and medicine safety stock to PD nominal demand increases slightly.Regarding the ratio of safety stock to SD nominal demand, that of medicine, water, and food increased linearly.Obviously, data variability has a greater impact on the ratio of safety stock to SD nominal demand.It can be seen from Figure 5(b-1) and (b-2) that with the increase in budget, the safety stock of water, food, and medicine increases.However, the ratio of bedclothes safety stock to PD and SD nominal demand decreases linearly.

Comparison with and without secondary disaster effect
If SDs are not considered in advance, the RS dispatching problem will be divided into the PD and SD response phases.RD after a PD will only meet the demand in time, and nothing else will be prepared for SDs.The corresponding mathematical model can be found in the supplementary materials.
Figure 6 reports the component of the total humanitarian cost of the two cases.Obviously, if decision-makers consider the SDs, the total fixed cost will decrease by 1.12E + 07 CNY for fewer SPs.The RSs are transported to affected areas to meet the PD demand and prepare for SDs.Therefore, the distribution cost in the PD response phase is 3.27% higher than that in the SD response phase, while the distribution cost in the SD response phase is 2.37E + 05 CNY lower.Therefore, the total penalty cost is 1.36% lower than that in the system without the SD effect.Considering SDs in HL will lead to shorter travel times and lower costs.

Comparison with and without collaborative distribution
According to the current administrative system of China, when a disaster occurs in an area, the local government in the affected area will orchestrate rescue operations alone.If disaster damage is so serious that local RSs nearly run out, the government will ask for additional help from successively higher levels of government.Relief demand will be satisfied until the RSs are out of stock all over the country.Suppose ε 2 = 1E + 10, the budget of uncertainty = 5, and data variability γ = 10%.If RSs are delivered independently, the HL cost is 9.30E + 09 CNY, and the travel time is at least 1.04E + 06 h.For the centralised decision-making scheme, the total cost of distributing relief collaboratively is 1.00E + 10 CNY, and it requires only 1.02E + 06 h.
Although centralised decision-making improves humanitarian logistics costs, it can also effectively improve the allocation efficiency of RSs and has practical application value.As shown in Figure 7, DPs can be served by RSFs at all levels at the same time, rather than facing RS shortages and having to wait for more material stored in higher-level RSFs.Obviously, if distributing relief collaboratively, DPs will be served by many more RSFs that are located close by.For the decentralised decision scheme, only a few municipalities and provincelevel RSFs are chosen, which may lead to a lower cost but take a long time to complete transportation.Then, some RSFs can form coalitions to manage relief inventory together, regardless of whether disasters occurred in the location of the RSFs.On the one hand, cooperation can achieve economies of scale to decrease purchase cost, delivery cost, holding cost, and so on.On the other hand, the involvement of more than one RSF will provide a robust and adjustable RD plan.In summary, decisionmakers should apply the collaborative mechanism in HL management.

Managerial insights
Based on the numerical results in Section 5, we obtained the following generalisable findings, which benefit decision-makers in improving their post-disaster RD considering uncertainties and SDs.
(i) According to the sensitivity analysis in Section 5.3, demand and travel time uncertainties (Figures 3 and 4) will increase the total travel time and RQ.Thus, novel approaches and technology should be adopted to predict reliable emergency demand and travel time.It makes sense to use GIS, GPS, and satellite big data analytics to obtain the state of transport infrastructure and repair damaged infrastructure for more real-time information (Nagendra, Narayanamurthy, and Moser 2022).Data-driven approaches can also improve management with tractable decisions (Chen et al. 2022), while reinforcement learning can learn from historical data to forecast and propose optimisation solutions under uncertainty (Yan et al. 2022).(ii) As shown in Sections 5.4 and 5.5, considering SDs reduces the total travel time and cost in the HL system (Figure 6).Therefore, it is significant for decision-makers to consider the SD effect in PD response operations.RSFs can preposition more RSs to reduce the penalty cost and improve demand satisfaction (Li, Zhang, and Yu 2020).RSs and safety stock can be transported together to achieve economies of scale.(iii) Regarding the findings in Section 5.6, Figure 7 illustrates that centralised decision-making is a more effective design for RD systems to make full use of resources than decentralised decision-making.HL systems can adopt centralised decision-making schemes.It is necessary to set up disaster response organisations that consist of various government departments.The leader groups contain decision-makers from different DPs and command emergency management operations both pre-and post-disaster, e.g.facility location, RD, and asset prepositioning.

Conclusion and future research
This study investigated a BOROM to locate RSFs at different levels and allocate RSs to DPs considering PD and SDs at the same time to minimise the total travel time and the total HL cost.The uncertainties of demand during and after DPs and required travel time between origindestination pairs are taken into consideration.To solve this model, first, we develop a nominal model for designing the relief collaborative distribution scheme.Second, we utilise RO methods to address left-hand side and right-hand side uncertainties and derive robust counterparts.Third, the epsilon-constraint method is used to solve the bi-objective problem and provides the Pareto front.A case study derived from the relief collaborative distribution in a Yunnan Province earthquake reveals some important practical insights and managerial implications.The proposed model adopted a collaborative mechanism considering subsequent shocks to make the network more robust to external uncertainties, which can be used in other cases to improve supply chain viability (Ivanov 2022), such as the ripple effect analysis in epidemic disruptions (Ivanov and Keskin 2023;Ivanov et al. 2023).
Based on the limitations of this paper, we suggest further research directions as follows.(i) The characteristics of items should be taken into consideration, such as replenishment or relief kits (Tang et al. 2022b;Zhang et al. 2023).(ii) To improve the model performance, novel approaches should be considered, such as distributionally robust optimisation, viable supply chain modelling with ripple effect, and reinforcement learning algorithms (Zekhnini et al. 2022;Brusset et al. 2022;Rolf et al. 2022).(iii) The proposed model provides only single-period decisions.Uncertain disruption duration lengths can be divided into several periods to capture dynamic disaster response operations (Li et al. 2023).

Disclosure statement
No potential conflict of interest was reported by the author(s).

Figure 1 .
Figure 1.Relief distribution network under a decentralised decision scheme (not collaborative) and a centralised decision scheme (collaborative).
Assumptions(i) The RS requirements of DPs and the required travel time between RSFs and DPs are both influenced by the scale of disasters.A higher scale of disaster incurs higher requirements.Therefore, to keep the emergency managerial solution more robust, we can assume that these parameters can be presented by uncertain intervals.Moreover, based on the information, the intervals can be determined and set as uncertain parameters (Sun, Wang, and Xue 2021).(ii) Because the environment of DPs is different, the type of damage in DPs caused by the PD and the corresponding SDs are usually also different.Therefore, DPs may require different kinds of RSs, such as food, water, and medical materials (Demirbas and Ertem 2021).(iii) Different kinds of RSs have different distribution costs and minimum satisfaction, which has an impact on the RD solution, creating different objective function values.

Figure 3 .
Figure 3. Effects of uncertainty on total travel time.

Figure 4 .
Figure 4. Effects of uncertainty on RQ.

Figure 5 .
Figure 5. Effects of data variability and budget of uncertainty on safety stock.

Figure 6 .
Figure 6.Comparison of the cost component of optimal solutions with and without the SD effect.

Figure 7 .
Figure 7. Maps of RSFs under a decentralised decision scheme (not collaborative) and a centralised decision scheme (collaborative).
The work was supported by the National Natural Science Foundation of China:[Grant Number 72074073]; the High-end Think Tank Project of Central South University:[Grant Number 2021znzk08]; Natural Science Foundation of Hunan Province, China: [Grant Number 2021JJ30857, 2021JJ31167]; Hunan Social Science Foundation: [Grant Number 19YBA378]; and Science and Technology Research and Development Plan Project of China National Railway Group Co.: [Grant Number 2023X022].Notes on contributorsDezhi Zhang is a professor at the School of Traffic & Transportation Engineering, Central South University.His research interests are green logistic network design and optimisation, green vehicle routing problem and schedule, optimisation of green supply chain network design, and optimisation of multimodal freight transport based on low-carbon.Yarui Zhang received the B.E. degrees from the School of Transportation and Logistics, Southwest Jiaotong University, China, in 2020, where she is currently pursuing the master's degree with the School of Traffic & Transportation Engineering, Central South University, China.Her research interest is emergency logistics system optimisation.Shuanglin Li received his Ph.D. in Logistics Engineering from Southwest Jiaotong University, Chengdu, China.He is currently an associate professor with the Central South University, Changsha, China.His research interests include emergency management, emergency logistics, logistics system optimisation, and operations management.Shuangyan Li is currently an associative professor with the College of Logistics and Transportation, Central South University of Forestry and Technology, China.Her research interests include green logistic network design and optimisation, green vehicle routing problem biomass supply chain network design and logistics system optimisation.Wanru Chen received the B.E. degrees from the School of Traffic & Transportation Engineering, Central South University, China, in 2020, where she is currently pursuing the doctorate's degree with the School of Traffic & Transportation Engineering.Her research interests are logistics system optimisation and intelligent optimisation algorithm.

Table 1 .
A brief review of humanitarian logistics system optimisation.
earthquakes in Yunnan Province, China, illustrate the robustness and advantages of the novel model.(iv) A sensitivity analysis of the budget of uncertainty, safety stock, is carried out to obtain valuable managerial insights.
The appropriate parameters ε 2 , ε 3 , . .., ε n should be determined in advance to obtain optimal solutions of the model.Constraints remain the same, and then a single objective optimisation model Q 1 with the optimal solution Z *

Table 2 .
Part of scheme of bi-objective robust model.

Table 3 .
Relief supply safety stock for secondary disasters.