1 Introduction

Current statistics and information show that unexpected terrestrial, oceanic and atmospheric phenomena are increasing and cause massive human and financial losses worldwide [32]. The great amounts of human and financial losses caused by natural disasters indicate the need to prepare for the disaster [30]. Disasters led to many economic and social challenges and issues. The Haiti earthquake had affected at least three million people and had killed between 217,000 and 230,000 [28]. Disasters like earthquakes impose human, financial and natural losses on society, economy, and nature [9, 48]. Therefore, efficient and comprehensive planning is needed in order to tackle and overcome the damage and losses caused by the earthquake [52]. In order to minimize the human and financial costs of the earthquake, the comprehensive recognition of the relief process and disaster relief supply chain is essential [38].

The amount of storage of relief commodities and the location of facilities and are among the plans that are considered in the pre-disaster preparation phase [35]. Besides, the routing, allocation, and distribution problems of relief commodities are essential and significant in the response phase [26, 51]. Lack of attention to these problems can be the main obstacle to an effective response to the disaster [12, 25]. In this paper, to estimate the demand, first, the urban infrastructure that is affected by the earthquake is identified. Thus, the infrastructures interact with each other using the simulation approach is calculated and finally, the demand is considered as an uncertain parameter.

A new multi-objective mathematical model is developed for pre-disaster and post-disaster phases. In the first phase, the Distribution Centers (DCs), the location, and the capacity of them are determined. In this research, three sorts of DCs are considered in terms of capacity, which include small DC (capacity from 2000 to 30,000), medium DC (capacity from 30,000 to 40,000), and large DC (capacity from 40,000 and 50,000). In the second phase, the allocation of the affected areas to cemeteries for transferring corpses and people to hospitals, routing the items from DCs to shelters and the victims’ evacuation from affected areas to shelters, and the allocation of the hospitals to affected areas for dispatching relief staff are provided. The proposed model integrates these two phases and made simultaneously. A customized Non-Dominated Sorting Genetic Algorithm (NSGA-II), Strength Pareto Evolutionary Algorithm (SPEA-II), and Pareto Envelope based on Selection Algorithm (PESA-II) as well as an efficient Epsilon-constraint approach are utilized for solving the presented mathematical programming model and to find Pareto solutions. A real case study is considered to validate the presented model.

The rest of the paper is addressed as follows: In Section 2, the literature review in the humanitarian relief logistic network field and solution methods are explained. The problem definition and mathematical modeling are discussed in Sections 3 and 4, respectively. In Sections 5 and 6, solution methods and real application are represented. Eventually, the conclusion and future works are reported in Section 7.

2 Literature review

The process of scheduling, organizing, and planning the resources to provide relief to affected humans by disasters is called emergency logistics [41]. Since the 1960s, the number of disasters has grown increasingly, affecting an increasing number of humans [49]. Due to uncertainties, emergency logistics is a challenging problem. Location-allocation-routing problems are considered an important and significant problem in the field of humanitarian logistics. Balcik and Beamon [4] used a covering location model considering capacity and budgetary constraints to provide relief commodities. They presented a multi-commodity scenario-based model in the response phase to maximize the level of demand coverage. Hence, their proposed model for a numerical example is solved by the exact method. The results indicate the efficiency of the presented model. Hong et al. [27] focused on the location of the DCs and the emergency supplies in disaster planning to increase disaster supply chain reliability. The proposed model for the case study in the US Southeast Region is solved using a heuristic approach. Jia et al. [29] formulated models to specify the optimal locations of emergency medical service (EMS) facilities such as air attacks. Authors provided a schedule in the preparation phase based on Petri net for a possible earthquake in Lushan/China. The presented deterministic, single-commodity, and single-period model was solved by the exact method and the results indicated its proper performance. Nedjati et al. [45] presented a deterministic, bi-objective routing and location problem of DCs in disaster preparedness and response phases. Replenishment at Intermediate Depots (CLRPR) was considered. The time window for the transportation of relief commodities, as well as the load and unload time for commodities, were considered in their study. To solve the NP-hard complexity, the NSGA-II algorithm was used.

Balcik [5] designed a deterministic mathematical model for site selection and routing decisions for assessment teams after the occurrence of a disaster (recovery phase). Their proposed model was solved by the Tabu search algorithm. The case study of this research was Van earthquake in Turkey in 2011. Goli et al. [22] presented a framework for supply chain managers in disaster, who face similar problems in other environments. The results indicated that trust and good communications were positive factors in performance of supply chain in the disaster scenario.

Tavana et al. [55] provided a multi-echelon, multi-commodities deterministic model to design the logistics network and location of warehouses in pre-disaster situations. Commodity inventory and routing of perishable commodities were of the problems that have been considered in the post-disaster phase. Their proposed mixed-integer linear programming (MILP) model was solved by utilizing the NSGA-II algorithm and Epsilon-constraint approach for a real case study. Nikoo et al. [46] presented a multi-objective deterministic model to design the transportation network for the response phase under possible earthquake in Tehran Province (Iran). The considered goals included identifying the optimal routes and the length of each route, as well as minimizing the response time. To solve the proposed model, a combination approach of lexicography and weighting was used. The results indicated the importance of decreasing the relief time in the early hours of the disaster due to the coming of voluntary aid to the affected area. Vahdani et al. [59] extended a stochastic mathematical model for the location of the DCs and warehouses of relief commodities. In their research, in the first phase, strategic decisions such as location and the establishment of DCs, and in the second phase, operational decisions such as routing, determining the reliability of the routes, observing the time window of each route were considered. Tlili et al. [58] investigated the routing of relief ambulances during the earthquake. The traffic constraint is considered as one of the contributions of their research. The injured people were divided into two categories of outpatient and severely injured people. The objective of this research was to minimize the relief time for both injured categories in each scenario. The case study consisted of 30 injured people, 7 ambulances, and 4 care centers in Jendouba city of Tunisia. Noham and Tzur [47] offered a scenario-based mathematical model for supply chain network management before and after the disaster. A new Mixed-Integer Linear Programming (MIP) model was formulated along with maximizing the number of allocated commodities per time unit. Also, the p-median approach was used to minimize the coverage radius. The case study considered was 49 potential scenarios for an earthquake in Israel. Maharjan and Hanaoka [43] suggested a model for the location of the temporary hubs under earthquake conditions in order to minimize costs and unsatisfied demands for injured people. To solve this two-objective model, the fuzzy weighting strategy was utilized. The case study was Nepal’s earthquake in 2015.

Liu et al. [42] formulated a mathematical model to minimize costs and unmet demand in post-disaster conditions (prevention phase). Uncertainties in demand and relief time were among the contributions of their study. Authors designed a two-echelon model including temporary facilities and affected areas serviced by helicopters. Accordingly, the proposed model for a case study in Sichuan province in China is provided and is solved by using a precise algorithm. Habibi-Kouchaksaraei et al. [24] presented a bi-objective, multi-echelon, and scenario-based model for supply chain management under disaster conditions (response phase). The proposed supply chain consisted of three levels: supply centers, processing, and DCs. The real case study considered in Mazandaran province/Iran and a robust optimization approach was utilized to solve it. Ghasemi et al. [19, 20] presented a multi-period and multi-objective mathematical model for location and allocation affected areas to hospitals. They considered the types of injured along with the destruction of established centers during an earthquake that was one of their contributions. Minimizing supply chain costs while minimizing the shortage of relief commodities has been their research aim. Their proposed model was solved using NSGA-II and modified multiple-objective particle swarm optimization approaches. The case study was in the city of Tehran and the results indicated the satisfactory performance of their proposed model.

Du et al. [13] provided a reliable p-center facility location to minimize costs in natural disasters (preparedness phase). Two approaches of constraint generation and Bender’s dual cutting plane were used for solving the proposed model. To assess the performance of the stochastic model, the presented approach was compared with a reliable p-median problem and a stochastic p-center problem. Tirkolaee et al. [56] presented a bi-objective mathematical model for disaster rescue units’ allocation and scheduling. Considering the effect of learning in the disaster management problem was one of the contributions of their paper. They applied robust optimization technique and multi-choice goal programming to cope with uncertainty of model. The case study was considered to be Mazandaran/Iran. Khalilpourazari et al. [33] presented an efficient blood supply chain network in disaster. They considered multi-echelon supply chain consist of donors, blood collection centers (permanent and temporary), regional blood centers, local blood centers, regional hospitals, and local hospitals. The main goals of the model were minimizing total transportation time and cost while minimizing unfulfilled demand. The model was solved using neural-learning method. The case study was considered to be earthquake in the Iran–Iraq border in 2017.

Ergün et al. [14] presented a game theoretical model for emergency logistic planning. The main goal of their paper was maximizing the transferred commodity. The case study was considered to be earthquake in Istanbul. The results indicated the proper performance of the proposed model. Alizadeh et al. [2] presented a multi-period model for locating relief facilities in natural disaster situations. The main purpose of their research is to maximize the coverage of hospitals and distribution centers. The Lagrangian approach has been used to solve the proposed model. The results of the case study indicate that increasing demand reduces the level of coverage of the regions. Zhan et al. [60] proposed a mathematical model for locating and allocating relief bases in the face of supply-demand uncertainty. One of the main goals of the research was to minimize the shortage and unsatisfied demand. The case study was for Zhejiang Province in southern China. The Particle Swarm Optimization approach was used to solve their model. The results showed that with increasing number of suppliers, the amount of unsatisfied demand decreases. Cavdur et al. [6] proposed new strategies for distributing relief goods and human resources in times of crisis. Therefore, a two-stage model for locating and allocating temporary relief centers was presented. The first stage was the pre-disaster phase and the second stage was the post-disaster phase. Minimizing chain costs was done in the first phase and allocating resources to centers in the second phase. The results of the sensitivity analysis of the case study indicated that with increasing demand, allocation costs increase exponentially. Ghasemi et al. [18] presented a location-routing-inventory decisions in humanitarian relief chain. The proposed network included affected areas, suppliers, distribution centers, and hospitals. It should be noted that in the second stage, a cooperative game theory of coalition type was considered, which resulted in synergies that minimize the relief golden time. Also, to validate the model, a real case study is provided for a possible earthquake in Tehran.

Hence, the details of the recent papers in the field of location-allocation-routing problems in disaster relief are shown in Table 1.

Table 1 The details and summarizes of the examined papers in the literature review
Full size table

The main contributions of this research in comparison with examined papers are summarized in Table 2.

Table 2 Contributions in this paper
Full size table

3 Problem definition

The main goal of the disaster relief system is to organize injured people and equipment in order to send relief efficiently. The candidate locations of DCs are determined. An optimal number of DCs and shelters should be determined and located. The demand points should be served by them. In this research, first, using the system dynamic approach, the interaction structure of the urban infrastructures in the event of an earthquake is drawn. Using this structure, the demand for relief commodities is estimated. The estimated distribution functions enter the mathematical model as the uncertain parameter in the simulation phase.

In this study, a new mathematical programming is developed to model the aforementioned situation in relief logistics network. Then, the suitable solution for the proposed model in presence of uncertainty is to be found using optimization techniques.

The disaster relief logistics network is assumed to consist of five main parts as 1) areas affected by disaster; 2) hospitals; 3) cemeteries; 4) shelters; and 5) DCs in this research (see Fig. 1).

Fig. 1
figure 1

The schematic view of relief logistic supply chain

Full size image

In this paper, the probability of the routes being destructed or blocked after the earthquake is considered. In the event of an earthquake, depending on its severity, there is a risk of the routes being destructed or blocked due to the traffic during rush hours. The values of this parameter will be changed based on the type of scenario (time and severity of the earthquake), the active fault, and the type of route (Asphalt Street, alley, dirt road, highway, freeway, etc.). For example, the probability of successful evacuation for a 7-magnitude earthquake on the highway is 85% and for the dirt road 45%.

All of the aforementioned tasks should be done within a limited time window after the occurrence of the disaster, which is called the golden time. The relief and rescue operations may become useless after the golden time. Therefore, the proposed model of this study determines the optimal locations for establishing shelters and DCs in the pre-disaster phase. Also, the suggested model determines the victims’ evacuation flow (from affected areas to shelters), relief staff flow, injured people flow (from affected areas to hospitals), corpse flow (from affected areas to cemeteries), relief commodity flow (from DCs to shelters), best routes for the evacuation of victims, distribution of relief items, and the allocation of the vehicles to these flows in the post-disaster phase.

3.1 Research framework

Figure 2 shows the research framework. In the first step, the basic infrastructures are identified in the event of an earthquake and their interactions are determined by the system dynamics. Then, the designed structure is simulated using Enterprise Dynamic simulation software. The output of this step is the estimate of distribution functions of relief commodities containing medicine, tents, food, and water. The distribution functions of relief commodities enter the mathematical model as input. In the second step, the mathematical model is presented for pre- and post-disaster phases. The purpose of this step is to plan location, allocation, distribution, and routing in order to minimize costs and unmet demand and maximize the probability of successful passage through routes. In the next step, the stochastic mathematical model presented by SCCPAis converted to a deterministic model. Finally, in the last step, the model is solved and optimized by three heuristic algorithms.

Fig. 2
figure 2

Research Framework in this paper

Full size image

3.2 Simulation structure

Here, the system dynamic structure of an earthquake is shown in Fig. 3. The considered main branches include the destruction of factory buildings, the destruction of houses, the reduction of the financial utility of companies, the outbreak of infectious diseases such as COVID-19, and the difficult conditions of the suppliers of relief commodities. Other considered infrastructures include transportation, water, fuel, financial and banking infrastructures, and so on. As can be seen, the infrastructure of roads and bridges is destroyed by the occurrence of the earthquake. Destruction of these infrastructures leads to the destruction of the transportation infrastructure and the reduction of the transportation capacity. Reduced transportation capacity leads to reduced port operations and also, logistics delays in the supply of raw materials. Eventually, logistics delays, as well as the destruction of factory buildings, cause damage to industries and products. Finally, the output of the simulation model can estimate the demand for the required relief commodities. The demand parameter is entered into the stochastic mathematical model as the uncertainty parameter after being estimated by the simulation model.

Fig. 3
figure 3

The structure of system dynamic

Full size image

3.3 Simulation model

Enterprise Dynamics (ED) software is a powerful software for simulating discrete event processes that were first developed by the Dutch company INCONTROL in 2003. ED is object-oriented software for modeling, simulating, visualizing, and controlling dynamic processes. The atoms present in this package have the ability to cover all the mentioned processes and due to the flexibility in their performance according to the needs of the model maker, this ability has reached its maximum. ED atoms are designed to implement systems and designs formed in the mind as easily as possible and with the necessary details in the software [53]. In this research, simulation input data were obtained from JICA (Japan International Cooperation Agency). As shown in Figs. 4, 31 server atoms, a source, and a sink are used. The demand server is utilized to estimate the distribution functions of relief commodities. Intended entities include earthquake signals in Richter. AvgContent (cs) is used for performance measurement (PFM).

Fig. 4
figure 4

Structure of simulation model

Full size image

In the simulation model, the inputs for various atoms are different. For example, here the inputs of the two atoms “Electric Grid trans-lines” and “factory building damaged” are described. The settings and inputs of the “factory building damaged” atom are based on Coburn and Spence [10]. This atom estimates the damage caused by the destruction of buildings according to Eq. (1):

$$ HC=R{I}_1.R{I}_2.R{I}_3.R{I}_4.\left( PI+\left(1- PI\right). CI\right) $$
(1)

In (1), HC determines the human casualties in buildings after an earthquake. RI1shows the number of trapped people in demolished buildings. RI2 indicates the percent of the building’s residents. RI3 specifies the number of collapsed buildings. RI4indicates the number of occupants of each building based on statistical data. PI represents the loss ratio after building demolition. Finally, the CI illustrates the proportion of injured who die before rescue teams arrive at the scene.

Also, the HAZUS damage function is used to study the “Electric Grid trans lines” atom. Equation (2) shows the amount of earthquake damage to power lines (RR):

The PGV is the peak ground velocity in centimeters per second and is considered as the input parameter of this atom.

$$ RR=0.3\ast 0.0001\ast PG{V}^{2.25} $$
(2)

Finally, the “Demand” atom has been used to estimate the number of needed commodities. In this study, it is assumed that each homeless person needs 5 l of water, 2.5 kg of food, 0.2 tents, and 0.5 kg of medicine. Moreover, the 4DScript codes written in the “Demand: medicine” atom for relief commodities are as follows (Warm up period = 1000 hours, observation period = 10,000 hours). The warm up period is the time that the simulation will run before starting to collect results. This allows the Queues (and other aspects in the simulation) to get into conditions that are typical of normal running conditions in the simulated system [54].

$$ {Server}^{``}{Demand}^{"}: Trigger\ on\ entry=0.5\ast \left( AvgContent\ \left( AtomByName\left(\left[ Medicine\right], Model\right)\right)\right). $$

4 Mathematical formulation

In this section, firstly, the notations of mathematical modeling are introduced. Then, the mathematical model is developed and discussed. The main assumptions of the proposed problem are reported in Table 3.

Table 3 The main assumptions of the proposed problem
Full size table

4.1 Indices, parameters, and decision variables

In this sub-section, the sets and indices are provided in Table 4. Table 5 presents the used parameters in the offered model. In addition, the decision variables in the proposed mathematical modeling are stated in Table 6.

Table 4 Sets and indices
Full size table
Table 5 Parameters and their definitions
Full size table
Table 6 Decision variables and their definitions
Full size table

4.2 Mathematical modeling

The formulation of the multi-objective location-allocation-routing problem mathematical modeling in the relief DCs network is as follows:

$$ {\displaystyle \begin{array}{c}\begin{array}{c}\operatorname{Minimize}\ {f}_1=\left\{\overset{K}{\sum \limits_{k=1}}\overset{S}{\sum \limits_{s=1}} fce{s}_k\times locs{h}_{ks}\right\}+\left\{\overset{K}{\sum \limits_{k=1}}\overset{S}{\sum \limits_{s=1}} vc p\times cap{s}_{ks}\right\}+\left\{\overset{L}{\sum \limits_{l=1}}\overset{J}{\sum \limits_{j=1}}\overset{S}{\sum \limits_{s=1}} fce{d}_{l\mathrm{j}}\times locd{c}_{ljs}\right\}\\ {}+\left\{\overset{S}{\sum \limits_{s=1}}\overset{V}{\sum \limits_{v=1}} vc{d}_v\times pr{o}_s\left[\overset{D}{\sum \limits_{d=1}}\overset{Z}{\sum \limits_{z=1}} disd{q}_{dz}\times vlcd{q}_{dz s}^v+\overset{D}{\sum \limits_{d=1}}\overset{K}{\sum \limits_{k=1}}\overset{E}{\sum \limits_{e=1}}\mathit{\ln}{g}_{dk}^e\times vlcs{f}_{dk s}^v+\overset{D}{\sum \limits_{d=1}}\overset{B}{\sum \limits_{b=1}} disd{g}_{dz}\times vlcd{c}_{db s}^v+\overset{Z}{\sum \limits_{z=1}}\overset{D}{\sum \limits_{d=1}} disd{q}_{dz}\times vlcq{d}_{zds}^v+\overset{O}{\sum \limits_{l=1}}\overset{P}{\sum \limits_{k=1}}\overset{T}{\sum \limits_{t=1}} disl{f}_{lk}\times vlcm{f}_{lk s}^{vt}\right]\right\}\end{array}\\ {}\left\{\overset{S}{\sum \limits_{s=1}}\overset{V}{\sum \limits_{v=1}}f{c}_v\times pr{o}_s\left[\overset{D}{\sum \limits_{d=1}}\overset{Z}{\sum \limits_{z=1}} vlcd{q}_{dz s}^v+\overset{D}{\sum \limits_{d=1}}\overset{K}{\sum \limits_{k=1}} vlcs{f}_{dk s}^v+\overset{D}{\sum \limits_{d=1}}\overset{B}{\sum \limits_{b=1}} vlcd{c}_{db s}^v+\overset{Z}{\sum \limits_{z=1}}\overset{D}{\sum \limits_{d=1}} vlcq{d}_{zds}^v+\overset{O}{\sum \limits_{l=1}}\overset{P}{\sum \limits_{k=1}}\overset{T}{\sum \limits_{t=1}} vlcm{f}_{lk s}^{vt}\right]\right\}+\\ {}\begin{array}{l}\left\{\underset{d,s,h}{\operatorname{Maximize}}\left\{0,\left[\left(\frac{nj{p}_{hds}}{\sum_{f=1} nd{c}_{1 ds}\ast ns{d}_{hf}}\right)+{\sum}_{Z=1}{\sum}_{e=1} path{\prime}_{dz s}^e\times {t}_{dz}\right]-72\right\}\times {w}_1\times \mathit{\cosh}\right\}\\ {}+\left\{\underset{d,s,g}{\operatorname{Maximize}}\left\{0,\left[\left(\frac{nh{p}_{ds}}{\sum_{n=1} nd{c}_{2 ds}\ast ns{n}_{gn}}\right)+{\sum}_{k=1}{\sum}_{e=1}t{\prime}_{dk}\times pat{h}_{dk s}^e\right]-100\right\}\times {w}_2\times cosg\right\}\\ {}+\left\{\underset{d,s,a}{\operatorname{Maximize}}\left\{0,\left[\left(\frac{ncrp{s}_{ds}}{\sum_{q=1} nd{c}_{3 ds}\ast ns{\mathrm{c}}_{aq}}\right)+{\sum}_{b=1}{\sum}_{e=1}{t}^{\prime {\prime}_{db}}\times pat{h}^{\prime {\prime}_{db s}^e}\right]-120\right\}\times {w}_3\times cosb\right\}\end{array}\end{array}} $$
(3)
$$ \operatorname{Minimize}\ {f}_2=\overset{S}{\sum \limits_{s=1}}\mathit{\Pr}{o}_s\times \overset{R}{\sum \limits_{r=1}}\underset{d}{\mathit{\operatorname{maximize}}}\ \left\{ sh{s}_{rds}\right\} $$
(4)
$$ \mathit{\operatorname{Minimize}}\ {f}_3=\overset{S}{\sum \limits_{s=1}}\overset{E}{\sum \limits_{e=1}}\overset{D}{\sum \limits_{d=1}}\overset{K}{\sum \limits_{k=1}}\overset{Z}{\sum \limits_{z=1}}\overset{B}{\sum \limits_{b=1}}\ \left\{\mathit{\Pr}{o}_s\times \left[\left(1- pr{s}_{dks}^e\right)\times pat{h}_{dks}^e+\left(1- pr s{\prime}_{dzs}^e\right)\times pat h{\prime}_{dzs}^e+\left(1- pr{s}^{\prime {\prime}_{dbs}^e}\right)\times {path}^{\prime }{\prime}_{dbs}^e\right]\right\} $$
(5)

Constraints:

$$ nh{l}_{dks}\le {M}_{Big}\times \overset{E}{\sum \limits_{e=1}} pat{h}_{dks}^e\times pr{s}_{dks}^e\kern1.00em \forall d,k,s $$
(6)
$$ nw{p}_{hdzs}\le {M}_{Big}\times \sum \limits_{e=1} path{\prime}_{dzs}^e\times prs{\prime}_{dzs}^e\kern1.00em \forall d,h,z,s $$
(7)
$$ ncrp{s}_{ds}\le {M}_{Big}\times \sum \limits_{e=1}{path}^{\prime }{\prime}_{dbs}^e\times {prs}^{\prime }{\prime}_{dbs}^e\kern1.00em \forall d,b,s $$
(8)
$$ {\sum}_{e=1} pat{h}_{dks}^e\le 1\kern0.75em \forall d,k,s $$
(9)
$$ {\sum}_{e=1} path{\prime}_{dzs}^e\le 1\kern0.75em \forall d,z,s $$
(10)
$$ \sum \limits_{e=1}{path}^{\prime }{\prime}_{dbs}^e\le 1\kern0.75em \forall d,b,s $$
(11)
$$ \overset{D}{\sum \limits_{d=1}} nh{l}_{dks}\le cap{s}_{ks}\kern0.5em \forall k,s $$
(12)
$$ \left\{\left(\overset{D}{\sum \limits_{d=1}} nh{l}_{dks}\right)\times {\overset{\sim }{nd}}_{its}\right\}= dm{r}_{ikts}\kern0.5em \forall k,t,i,s $$
(13)
$$ nd{c}_{rds}- nr{s}_{rds}= srp{d}_{rds}- sh{s}_{rds}\kern0.5em \forall d,r,s $$
(14)
$$ \overset{t}{\sum \limits_{t=1}}\overset{L}{\sum \limits_{l=1}} ri{t}_{ilks}^t-\overset{t}{\sum \limits_{t=1}} dm{r}_{ikts}= srp{f}_{ikts}- sh{c}_{ikts}\kern1.50em \forall k,i,t,s $$
(15)
$$ \overset{D}{\sum \limits_{d=1}} nd{c}_{rds}\le st{f}_{rzs}\kern0.5em \forall r,z,s $$
(16)
$$ \overset{I}{\sum \limits_{i=1}} ri{t}_{ilks}^t\le \left(\overset{J}{\sum \limits_{j=1}} capd{c}_j\times locd{c}_{ljs}\times c{o}_i\right)\kern0.5em \forall l,k,t,s\kern0.5em $$
(17)
$$ \left(\overset{Z}{\sum \limits_{z=1}} nw{p}_{hdzs}\right)+ sh{w}_{hds}= nj{p}_{hds}\kern0.5em \forall d,h,s\kern0.75em $$
(18)
$$ \overset{Z}{\sum \limits_{z=1}} nw{p}_{hdzs}\le nj{p}_{hds}\kern0.5em \forall d,h,s\kern0.75em $$
(19)
$$ \overset{K}{\sum \limits_{k=1}} nh{l}_{dks}= nh{p}_{ds}\kern0.5em \forall d,s\kern0.5em $$
(20)
$$ \overset{B}{\sum \limits_{b=1}} cr{p}_{dbs}= ncrp{s}_{ds}\kern0.75em \forall d,s\kern0.5em $$
(21)
$$ \overset{J}{\sum \limits_{j=1}} locd{c}_{ljs}\le 1\kern0.5em \forall l,s $$
(22)
$$ \overset{E}{\sum \limits_{e=1}} pat{h}_{dks}^e\le locs{h}_{ks}\kern0.5em \forall d,k,s\kern0.5em $$
(23)
$$ cap{s}_{ks}\le \left({M}_{Big}\times locs{h}_{ks}\right)\kern1.00em \forall k,s $$
(24)
$$ \left(\overset{I}{\sum \limits_{i=1}} ri{t}_{ilks}^t\times vl{u}_i\right)\le \left(\overset{V}{\sum \limits_{v=1}} vlcm{f}_{lks}^{vt}\times cp{v}^v\right)\kern0.5em \forall l,k,t,s\kern0.5em $$
(25)
$$ \left(\overset{I}{\sum \limits_{i=1}} ri{t}_{ilks}^t\times wu{t}_i\right)\le \left(\overset{V}{\sum \limits_{v=1}} vlcm{f}_{lks}^{vt}\times cp{w}^v\right)\kern0.5em \forall l,k,t,s\kern0.5em $$
(26)
$$ \overset{D}{\sum \limits_{d=1}} nw{p}_{hdzs}\le be{d}_{hzs}\kern0.5em \forall h,z,s\kern0.75em $$
(27)
$$ cr{p}_{db s}\le {M}_{Big}\times al{c}_{db}\kern0.5em \forall d,b,s\kern0.5em $$
(28)
$$ cr{p}_{dbs}\le \left(\overset{V}{\sum \limits_{v=1}} vlcd{c}_{dbs}^v\times cp{c}^v\right)\kern0.5em \forall d,b,s $$
(29)
$$ \overset{H}{\sum \limits_{h=1}} nw{p}_{hdzs}\le \left(\overset{V}{\sum \limits_{v=1}} vlcd{q}_{dzs}^v\times cp{h}^v\right)\kern0.75em \forall d,z,s\kern0.5em $$
(30)
$$ nh{l}_{dks}\le \left(\overset{V}{\sum \limits_{v=1}} vlcs{f}_{dks}^v\times cp{p}^v\right)\kern0.5em \forall d,k,s\kern0.5em $$
(31)
$$ \overset{R}{\sum \limits_{r=1}} nd{c}_{rds}\le \left(\overset{V}{\sum \limits_{v=1}} vlcq{d}_{zds}^v\times cp{m}^v\right)\kern0.5em \forall d,z,s $$
(32)
$$ {\displaystyle \begin{array}{c}\begin{array}{l} dm{r}_{ikts}, sh{c}_{ikts}, sh{s}_{rds}, sh{w}_{hds}, srp{f}_{ikts}, srp{d}_{rds}, vlcs{f}_{dks}^v,\\ {} vlcd{c}_{dbs}^v, vlcd{q}_{dzs}^v, vlcq{d}_{zds}^v, vlcm{f}_{ops}^{vt}, ri{t}_{ilks}^t, nd{c}_{rds}, nw{p}_{hdzs},\kern0.5em \end{array}\\ {} nh{l}_{dks}, cr{p}_{dbs}, cap{s}_{ks}\in \left\{{R}^{+}\right\}\end{array}} $$
(33)
$$ locs{h}_{ks}, locd{c}_{ljs}, pat{h}_{dks}^e, pat h{\prime}_{dzs}^e,{path}^{\prime }{\prime}_{dbs}^e\in \left\{0,1\right\} $$
(34)

The objective function (3) minimizes the expected value of the total costs of the relief supply chain. It includes the total costs of pre-disaster activities (including fixed costs and variable costs of establishing shelters and DCs) and the costs of post-disaster activities (including fixed costs and variable costs of transportation vehicles). The second part of the objective function is to decrease the maximum cost of relief. The expression \( {\mathrm{njp}}_{\mathrm{hds}}/\left(\sum \limits_{\mathrm{f}=1}{\mathrm{ndc}}_{1\mathrm{ds}}.{\mathrm{nsd}}_{\mathrm{hf}}\right) \)specifies the time of serving the injured people. The expression \( {\sum}_{\mathrm{z}=1}{\sum}_{\mathrm{e}=1}\mathrm{path}{\prime}_{\mathrm{dz}\mathrm{s}}^{\mathrm{e}}.{\mathrm{t}}_{\mathrm{dz}} \) calculates the total transportation time of the injured people. The difference of this amount from golden time (72 hours) indicates the amount of time violation from the relief. Multiplying the amount of time violation in the priority factor and the cost of not giving relief w1coshto the injured person determines the cost of serving the injured person. Similarly, subsequent expressions show the cost of not serving the homeless people and the corpses. The amount of golden time for serving the homeless is 100 hours and for the corpses, 120 hours. The priority weighting factor of injured people determines the priority level of people in relief. Accordingly, the priority factor for serious injury, homelessness and death is 0.6, 0.3 and 0.1, respectively. These values are based on Ghasemi et al. [19]. Accordingly, the first objective function also guarantees fairness in evacuation.

The objective function (4) shows to decrease the maximum number of unsatisfied demands for relief staff in each affected area under all scenarios. The objective function (5) minimizes the total probability of unsuccessful evacuation in routes. It should be noted that the probability of successful evacuation in routes has been estimated by road construction experts and earthquake engineers who have sufficient experience in route strength. Therefore, the third objective function seeks to reduce the probability of crossing routes with less reliability. For example, if in scenario S there are two routes among the affected area and the shelter with probability of successful evacuation \( \left({prs}_{dks}^e\right) \) equal to 0.3 and 0.6, then the model selects the route with the highest probability of successful evacuation. In fact, the third objective function minimizes the possibility of an unsuccessful evacuation.

According to the parameterspros,\( {prs}_{dks}^e \),\( {prs}_{dzs}^{\prime e} \) and \( {prs}_{dbs}^{\prime \prime e} \)are of the probability type and the binary variables \( {path}_{dks}^e \), \( {path}_{dzs}^{\prime e} \), and \( {path}_{dbs}^{\prime \prime e} \)are scalable, so the objective function is of the probability type.

Constraints (6–8) prevent the use of paths between two nodes that are not exist. Constraints (9)–(11) guarantee only one evacuation path is chosen between two nodes. Given the constraints (9)–(11) and the structure of the model, transportation in this problem is possible only between two points; for example, between hospital and affected area, shelter and affected area, cemetery, and affected area. So, given that the allocation of these locations to vehicles is one-to-one, and knowing that the creation of a sub-tour requires at least three nodes; there is no need to write the sub-tour constraint separately in the problem. So, given the nature of the problem, and as it is clear in the results, no loop formation is possible in the route of vehicles. Constraint (12) shows the capacity constraint for shelters. Constraint (13) satisfies the demand of each shelter. The unsatisfied demand for relief staff at affected areas is investigated in constraint (14). Constraint (15) guarantees unsatisfied demand for commodities at shelters. Constraint (16) shows the dispatched relief staff does not exceed the available relief staff of the current hospital. Constraint (17) indicates the dispatched commodities do not exceed the available commodities of current DC. Constraint (18) determines unserved injured people in affected areas. Constraint (19) illustrates the amount of injured people dispatched from a specific affected area does not exceed the number of injured people from that affected area. Constraint (20) ensures that all the homeless people are evacuated from affected areas. Constraint (21) shows all corpses are transferred from affected areas. Constraint (22) prevents establishing more than one DC at any node. Constraint (23) indicates the establishment of a new shelter before transferring homeless people. Constraint (24) represents a new shelter should be opened before it is used. Constraint (25) ensures that the volume capacity of vehicles. Constraint (26) shows the weight capacity of vehicles. Constraint (27) displays the maximum allowed capacity of the hospitals. Constraint (28) ensures that the corpses of each affected area are carried to allowed cemeteries. The number of corpses considering the capacity of vehicles, the number of injured people based on the capacity of vehicles, the number of homeless people considering the capacity of vehicles, and the number of relief staff considering the capacity of vehicles, respectively are shown in constraints (29)–(32). Constraints (33)–(34) indicate the type of decision variables.

4.3 Stochastic chance constraint programming

In this study, Stochastic Chance Constraint Programming Approach (SCCPA) has been used to convert stochastic model to deterministic model. This approach was stated by Charnes and Cooper [7]. Numerous and successful applications of this approach have been reported. Many papers have been used the SCCPA to convert stochastic technique to deterministic model [1, 21, 39, 40]. Assume that ~is the symbol of an uncertain parameter and k is the number of objective functions. Also assume at least one of the parameters of aij, hi or cjk defined as stochastic parameter. Therefore, the uncertain model is considered as follows:

$$ \mathit{\min}{f}_k=E\left(\overset{n}{\sum \limits_{j=1}}{c}_{kj}^{\sim }{z}_j\ge {h}_i^{\sim}\right)\kern0.5em k=1,\dots, K;i=1,2,\dots, m $$
(35)
$$ i=1,2,\dots, m\kern0.5em p\left(\overset{n}{\sum \limits_{j=1}}{a}_{ij}^{\sim }{z}_j\ge {h}_i^{\sim}\right)\ge {\alpha}_i $$
(36)
$$ z=\left({z}_1,\dots, {z}_n\right) $$
(37)
$$ z\in S $$
(38)
$$ z\ge 0 $$
(39)

where ckj shows the benefit ratio of jth decision variable in kth objective function. aij, hi, and yj indicate the technology coefficient of jth decision variable in ith constraint, the right hand side of constraint i, and jth decision variable, respectively.

Finally, the deterministic result of the general model for the maximum and minimum states is as Eqs. (40)–(42):

$$ E\left(\overset{n}{\sum \limits_{j=1}}{d}_{kj}^{\ast }{z}_j-{f}_k^{-}\right)-{\varphi}^{-1}\left({\alpha}_k\right)\sqrt{Var\Big(}\overset{n}{\sum \limits_{j=1}}{d}_{kj}^{\ast }{z}_j-{f}_k^{-}\Big)\ge 0\kern0.5em k=1,\dots, K $$
(40)

(Therefore,)

$$ {f}_k^{-}=\min \overset{n}{\sum \limits_{j=1}}{d}_{kj}^{\ast }{z}_j $$
$$ E\left(\overset{n}{\sum \limits_{j=1}}{d}_{kj}^{\ast }{z}_j-{f}_k^{+}\right)+{\varphi}^{-1}\left({\alpha}_k\right)\sqrt{Var\Big(}\overset{n}{\sum \limits_{j=1}}{d}_{kj}^{\ast }{z}_j-{f}_k^{+}\Big)\le 0\kern0.5em k=1,\dots, K $$
(41)

(Moreover,)

$$ {f}_k^{+}=\max {d}_{kj}^{\ast }{z}_j $$
$$ E\left(\overset{n}{\sum \limits_{j=1}}{a}_{ij}^{\sim }{z}_j-{h}_i^{\sim}\right)-{\varphi}^{-1}\left(1-{\alpha}_i\right)\sqrt{Var\Big(}\overset{n}{\sum \limits_{j=1}}{a}_{ij}^{\sim }{z}_j-{h}_i^{\sim}\Big)\ge 0\kern0.5em i=1,2,\dots, m $$
(42)

Based on Eq. 13, the multi-objective chance constraint model can be converted as a deterministic model at α%level as constraint (43):

$$ \left(\overset{\mathrm{D}}{\sum \limits_{\mathrm{d}=1}} nh{l}_{\mathrm{d}\mathrm{ks}}\right).E\left(\overset{\sim }{n{d}_{its}}\right)+{\varphi}^{-1}\left(1-{\alpha}_i\right).\sqrt{\operatorname{var}\overset{\sim }{n{d}_{its}}}= dm{r}_{ikts}\kern0.5em \forall k,t,i,s $$
(43)

where \( E\left({\overset{\sim }{nd}}_{its}\right) \) and \( Var\left({\overset{\sim }{nd}}_{its}\right) \) illustrate mean and variance of Normal distribution function estimated by simulation model (see Table 20 in section 7.1.1), α shows normal distribution at confidence level (1 -α)%, and φindicate standard normal distribution, with zero mean and unit variance.

5 Proposed solution methods

Since the proposed problem is NP-hard, to solve mathematical modeling, the meta-heuristic algorithms are employed in this paper. The NSGA-II, SPEA-II, and PESA-II algorithm and Epsilon-constraint method are utilized to solve the model in different sizes.

5.1 Proposed NSGA-II

The NSGA-II is one of the methods widely utilized by managers and engineers for solving multi-objective problems [36]. This algorithm has shown great accuracy in finding solutions to large-scale problems by examining most of the possible states which by humans is almost impossible.

5.1.1 Chromosome representation

In this research, there are h (h = 1,…,H) injured people, v (v = 1,…,V) vehicles, and finally m (m = 1, …,M) number of predetermined points and a candidate for the establishment of DCs and shelters for distribution of commodities. There is also a chromosome with (2 V + H) gene. This chromosome is shown in Fig. 5.

Fig. 5
figure 5

Chromosome representation

Full size image

Designed chromosome consists of 3 sections. The number of genes in the first section of this chromosome is V, and is filled with random numbers between 1 and M. Repeated numbers are allowed in the first section, which means that several vehicles are allocated to a DC and a shelter (see Fig. 5). The second section of the chromosome is designed to have V genes, and its alleles are integer and non-repetitive random numbers from 1 to M. Non-repeatability of alleles of the second section means that injured people can’t be served more than once by a vehicle. For example, as shown in Fig. 5, the second section of the chromosome has 4 genes and is completed with numbers from 1 to 9. According to the second section of the chromosome, the first vehicle visits the fourth injured person on its first move, the second vehicle visits the fifth injured person on its first move, and eventually, the fourth vehicle visits the second injured person on its first move. As it is obvious, the second section of the chromosome is relevant to the first injured person that each vehicle visits. The third section of the chromosome has H genes, which are shown with non-repetitive numbers like the first section. This section shows the serving sequence to the injured people.

5.1.2 Crossover operator

A double-point crossover operator is designed for generating off-springs using three-part chromosome on the basis of a crossover probability. Figure 6 indicates the schematic view of the proposed crossover operator.

Fig. 6
figure 6

Double-point crossover

Full size image

5.1.3 Mutation operator

As shown in Fig. 7, the mutation operator performed its duty according to mutation rate. The two genes are considered in the third section of the designed chromosome (for example, genes between 2 V + 1 and 2 V + H). Then, as shown in Fig. 7, the swap operator changes the alleles of the desired genes.

Fig. 7
figure 7

Mutation operator

Full size image

5.1.4 Stopping criterion

If any of the following two conditions are met, the algorithm stops:

1-If a certain number of iterations are satisfied.

2-If there is no improvement after a certain number of iterations.

5.2 Proposed SPEA-II

SPEA-II is an efficient algorithm that employs an external archive to store the non-dominated solutions which is found while searching the algorithm. Hence, Zitzler et al. [63] proposed the SPEA-II algorithm for the first time. The framework of this algorithm is described below:

  1. Step 1.

    Generating an initial population of solutions P0and set E0 = ∅

  2. Step 2.

    Computing the fitness of each solution i in set (Pt ∪ Et) according to the Eq. (44).

  • First, calculate an initial fitness of solution i based on Eq. (44):

$$ R(i)=\sum \limits_{j\in {P}_t\cup {E}_t\&j>i}s(j) $$
(44)

where j > i means that the j solution overcomes on the i solution. Also, s(j) shows the strength of the solutions, which is obtained from Eq. (45). Indeed, calculate the number of solutions that are dominated by i.

$$ s(i)=\left|\left\{j|j\in {P}_t\cup {E}_t\&i>j\right\}\right| $$
(45)

  • Calculate the congestion of the i solution according to the Eq. (46):

$$ D(i)=\frac{1}{\sigma_i^k+2} $$
(46)

where \( {\sigma}_i^k \) is the distance among i solution and kth neighborhood is close to it.

  • Finally, the fitness value is obtained from the sum an initial fitness value and congestion of the i solution:

$$ F(i)=R(i)+D(i) $$
(47)
  1. Step 3.

    Copying all non-dominated solutions in the set (Pt ∪ Et) to Et + 1.

  2. Step 4.

    If provided stop conditions, the algorithm stops and returns Et + 1 solutions.

  3. Step 5.

    Using the binary tournament method, the parent is selected of sets Et + 1.

  4. Step 6.

    Employing crossover and mutation operators on parent and generating offspring as many as NP. Offspring copy to the set Pt + 1 and add to the counter value of one unit and go to step 2.

It should be noted that this algorithm utilizes the same approach of crossover and mutation as described in Sections 5.1.3 and 5.1.4.

5.3 Proposed PESA-II

One of the most well-known multi-objective algorithms is the second version of the Pareto envelops based selection algorithm (PESA-II), which utilizes genetic algorithm operators to generate new solutions. The PESA-II was introduced by Corne et al. [11]. In the following, the steps of this algorithm describe:

  1. Step 1.

    Starts with a random initial population (P0)and sets the external archive (E0)equal to empty and the t = 0 counter.

  2. Step 2.

    Divides the solution space to hypercube nk, so n is the number of networks in each axis of the objective functions and k is the number of objectives.

  3. Step 3.

    Combines non-dominated solutions archive Et with new solutions of Pt.

  4. Step 4.

    If the stop condition is met, stop and indicate the final Et.

  5. Step 5.

    By setting Pt = ∅, solutions of Et is select for combining and mutation based on congestion information of hypercube. This selection is made utilizing the roulette wheel so that the probability of choosing hypercubes with a smaller population is more. Use the combination and mutation to generate offspring Np and copy it into Pt + 1.

  6. Step 6.

    Set t to t + 1 and go to the step 3.

For further details of PESA-II and SPEA-II algorithms could refer to references Gadhvi et al. [17]; Anusha and Sathiaseelan [3]; Kumar and Guria [37]; Moayedikia [44]; Chen and Li [8].

5.4 Epsilon-constraint method

The Epsilon-constraint method is one the most successful exact multi-objective optimization techniques. In the Epsilon-constraint method, one of the objective functions will be selected to be optimized and the other objective functions will act as the constraints [15].

In mathematical terms, If fj(x)j ∈ {1, …, k} the jth objective function chosen to be optimized, the multi-objective optimization is replaced with the following single objective optimization problem:

$$ {\displaystyle \begin{array}{l}\operatorname{Minimize}\ {f}_j(x)\ j\in \left\{1,\dots, k\right\}\\ {}s.t.\\ {}\kern1em {f}_i(x)\le {\varepsilon}_i,\kern0.5em \forall i\in \left\{1,\dots, k\right\},i\ne j\\ {}\kern1em x\in S\end{array}} $$
(48)

where, S is the feasible solution space, and εi is the upper bound of ith objective function. One advantage of the Epsilon-constraint approach is that it can be achieved effective points in a non-convex Pareto curve. Therefore, the decision maker can change the upper bounds εi to achieve weak Pareto optima [31]. There are several successful applications of the Epsilon-constraint approach (see [50]).

6 Computational experiments

According to the problem is NP-hard, the Epsilon-constraint approach is not capable for solving the proposed problem on a large-scale. According to Table 7, sample problems are classified into two groups of small and medium. In addition, the value of the proposed parameters is reported in Table 8. Therefore, after proving the efficiency of proposed meta-heuristic algorithms for small- and medium-scale problems, the proposed meta-heuristic algorithms are used to find the Pareto front of large-scale problems.

Table 7 The different size of the problems
Full size table
Table 8 The value of the proposed parameters
Full size table

The Epsilon-constraint method is coded in GAMS 24 and three meta-heuristic algorithms are coded in MATLAB R2020a software on a 2.3 GHz laptop computer with 4 GB of RAM using Win 7, 64bit. Due to the nonlinearity of the proposed model for small- and medium-scale, the model is solved with BARON solver. The proposed meta-heuristic algorithms codes were run using several combinations of parameters to obtain the optimum configuration of parameters. The proper configurations of parameters for NSGA-II, SPEA-II, and PESA-II are listed in Table 9.

Table 9 Fitted parameters of the suggested algorithms
Full size table

Tables 10 and 11 show the efficiency of proposed Epsilon-constraint and meta-heuristic algorithms in small and medium-scale problems. The first column of the table indicates the number of problems. The first 5 problems are small-scale problems and the problems of 5 to 10 are medium-scale problems.

Table 10 The results of Epsilon-constraint for small and medium size test problems
Full size table
Table 11 The results of meta-heuristic algorithms for small and medium-size test problems
Full size table

As shown in Table 11, the error average of the objective functions for the NSGA-II is 0.03, 0.63, and 0.6, for the SPEA-II approach is 0.06, 1.63, and 2.13, and for the PESA-II approach is 0.04, 1.29, and 1.58, respectively. As is known, the error rate of the NSGA-II is less than SPEA-II and PESA-II algorithms.

According to the obtained error values, efficiency and reliability of NSGA-II algorithm for small and medium-scale problems are proven. The maximum mean error in the NSGA-II algorithm is 0.63%. So, with this argument, the proposed NSGA-II method can be trusted to solve large-scale problems.

Figure 8 indicates the outcomes of CPU time for small- and medium-scale problems according to the presented methods. The comparison between the CPU time of NSGA-II, PESA-II, and SPEA –II methods shows that the NSGA-II CPU time is lower. Also, the CPU time of the Epsilon-constraint approach for small- and medium-size problems has increased exponentially, which proves the NP-hardness of the proposed model.

Fig. 8
figure 8

The results of CPU time for small- and medium-scale problems of the proposed methods

Full size image

6.1 Taguchi method

This sub-section sets out an, experimental design for controlling and tuning the algorithm’s parameters. In this paper, the Taguchi approach is used to tune the parameters of the algorithms set out. For more details about this method, interested researchers can examine a few papers related to the Taguchi method, such as Tirkolaee et al. [57]; Goodarzian et al. [23], etc.

The Signal-to-Noise ratio (S/N) is used to analyze the experiment. The S/N value expresses the degree of scatter around a certain value or, in other words, how our solutions have changed across several experiments. To obtain the S/N value, there are three equations, each of which applies to specific conditions.

$$ SB=\frac{1}{n}\sum {\left({y}_i\right)}^2 $$
(49)
$$ LB=\frac{1}{n}\sum {\left(\frac{1}{y_i}\right)}^2 $$
(50)
$$ NB=\frac{1}{n}\sum {\left({y}_i-{y}_0\right)}^2 $$
(51)

Moreover, the Taguchi approach mechanism focuses on the sort of solution. Therefore, the solution obtained that is relevant to three Eqs. (49)–(51) of this method is divided into three groups including the smaller is better (see Eq. (49)), larger is better (see Eq. (50), and the nominal is better (see Eq. (51)). In the proposed model, objective functions are a type of minimization, to control the parameters of the algorithm set out, the smaller better is used. Thus, Eq. (52) shows the S/N ratio value set out in this paper.

$$ \frac{S}{N}=-10\times \mathit{\log}\left(\frac{\sum_{i=1}^n{Y}_i^2}{n}\right) $$
(52)

where n shows the orthogonal array as well as yirepresents the value of the solution for the orthogonal array.

Since the scale of objective functions in each example is various, they could not be used directly. Accordingly, the Relative Percent Deviation (RPD) is used for each example to solve this problem. The RPD value for the data is obtained using Eq. (53).

$$ RPD=\frac{Alg_{sol}-{\mathit{\operatorname{Min}}}_{sol}}{{\mathit{\operatorname{Min}}}_{sol}}\times 100 $$
(53)

where Minsol and Algsol show the achieved best solution and the values of the achieved objective for each iteration of the experiment in a provided example, respectively. Therefore, the mean RPD is computed for each experiment after transforming the values of the objective to RPDs.

The provided levels with the factors (the algorithm’s parameters) are reported in Table 12, which for algorithm’s factors, three levels are provided. It should be noted that the Taguchi approach reduces the total number of tests by presenting a set of orthogonal arrays to control the proposed algorithms in a reasonable time. This approach suggests L9 for three VNS, GWO, and MGWO (see Table 12) obtained by Minitab Software.

Table 12 The orthogonal array L9 for NSGA-II, SPEA-II, and PESA-II
Full size table

Therefore, the output of the S/N ratio should be analyzed by Minitab Software to detect the best levels for each algorithm (see Figs. 9, 10 and 11). Where the S/N index has reached its minimum, the levels can be selected as the optimal levels.

Fig. 9
figure 9

Minitab Software output for the Taguchi method of the NSGA-II algorithm

Full size image
Fig. 10
figure 10

Minitab Software output for the Taguchi method of the SPEA-II algorithm

Full size image
Fig. 11
figure 11

Minitab Software output for the Taguchi method of the PESA-II algorithm

Full size image

The RPD is also used to confirm the selected best factors based on S/N ratios. Figure 9 shows shows the outcomes of RPD for each parameter level. It is clear that the RPD shows the best factors in Fig. 12, which confirm the same outcomes as S/N ratios.

Fig. 12
figure 12

Mean RPD plot for each level of the factors

Full size image

6.2 Metrics to evaluate algorithm efficiency

In order to assess the efficiencies of the multi-objective meta-heuristic algorithms, two assessment metrics are utilized as follows.

  • Spacing: measures the standard deviation of the distances between solutions of the Pareto front [62].

  • Mean ideal distance (MID): measures the convergence rate of Pareto fronts up to a certain point (0, 0) [61].

Results of the comparison of four methods of epsilon-constraint, NSGA-II, SPEA-II, and PESA-II methods using the metrics of MID and SM are shown in Table 13.

Table 13 The results of the obtained assessment metrics of the proposed algorithms
Full size table

As shown in Table 13, samples 1 to 5 are small-scale samples and samples 6 to 10 are in medium-scale. In each group of samples, the problem scale increases gradually with increasing the sample number. Gradually, with the increase in the scale of the problem to sample number 5, the difference remains almost the same. But from the sample 5, due to the increase in the scale of the problem from small to medium, significant changes are made in the values ​​of the metrics. For example, the SM value of the NSGA-II method increases from 0.383 to 0.426, and the MID for NSGA-II increases from 6.44 to 6.63.

The mean value of SM for the NSGA-II method is 0.402, for the SPEA-II algorithm is 0.4944, and for the PESA-II algorithm is 0.5106. This suggests that Spacing Metric for NSGA-II is slightly better than the two presented meta-heuristic algorithms. Moreover, it can be concluded that the outcomes of MID and SM for NSGA-II are slightly better than the other two presented meta-heuristic algorithms for small- and medium-scale problems.

In order to use the solution space in a better way, the number of Pareto solutions is considered to be 50. If the number of calculated Pareto solutions is more than 50, only the best 50 solutions are selected by the crowding distance operator.

Moreover, a set of statistical comparisons between the proposed meta-heuristic algorithms according to the assessment metrics are performed to realize which the presented meta-heuristic is the best. Furthermore, the attained outputs of the suggested problems are transformed to Relative Deviation Index (RDI) that is calculated based on Eq. (54):

$$ RDI=\frac{\left|{Alg}_{sol}-{Best}_{sol}\right|}{{\mathit{\operatorname{Max}}}_{sol}-{\mathit{\operatorname{Min}}}_{sol}}\times 100 $$
(54)

where Bestsol shows the best solution between the algorithms. Algsol indicates the value of the objective function. Minsol and Maxsol illustrate the minimum and maximum values of the proposed assessment metrics. Accordingly, the confidence interval of 95% for the assessment metrics in the algorithms set out is conducted to analyze the efficiency of the presented meta-heuristics statistically. Figure 13 depicts the means plot and LSD intervals for the proposed meta-heuristics.

Fig. 13
figure 13

The means plot and LSD intervals for the suggested algorithms

Full size image

It should be noted that the lower value of RDI is better. In Fig. 13, in terms of the spacing and MID metrics, the value of RDI for NSGA-II is lower than the PESA-II and SPEA-II algorithms and indicates has the best statistical performance. But, PESA-II has the high value of the RDI and shows the worst performance.

7 Real application

In this section, a real case study is introduced and discussed according to the proposed problem. The proposed model and solution procedures are applied. The results and some analysis of the results are further investigated.

7.1 Case study

Iran is one of Asia’s most populous and developing countries. Tehran is the capital of Iran and has about 9 million people. Tehran is among the rapid growing cities in Asia. Tehran is divided into 22 regions, 134 areas and 374 neighborhoods. The 2017 census shows that region 6 has the highest density with 434 people per hectare. Region 7 with 402 people per hectare is next. Regions 1, 4 and 8 with an average density of 250 people per hectare are among the relatively densely populated regions. Figure 14 shows the density ratio of buildings in the 22 regions of Tehran city on its faults. Due to the location of these two regions on the faults of Tehran city, the high population density in these two regions and the location of many organizations and schools in them, they are among the most vulnerable regions of Tehran. Therefore, due to the explanations presented; these two regions are considered as case studies.

Fig. 14
figure 14

Density ratio of buildings in 22 regions of Tehran on faults (Firuzi et al. [16])

Full size image

Regions 6 and 7 are the most populated and most crowded regions of Tehran including 36 hospitals complexes. Map of regions 6 and 7 of Tehran is presented in Fig. 15, which includes 35 candidate locations for establishing shelters, 45 candidate locations for establishing DCs, 36 hospitals, 30 affected areas and 20 cemeteries.

Fig. 15
figure 15

Case Study Map

Full size image

It is clear that in Fig. 15, the H points indicate the location of the hospitals, the D points display candidate locations for setup DCs, C points indicate the cemetery locations, A points show affected areas and S points show the candidate location for setup shelters.

Table 14 shows the scenarios of this research and their probability of occurrence. The probability of occurrence of each scenario is obtained by examining the time and severity of earthquakes in Tehran over the past 200 years [34].

Table 14 The scenarios investigated in the case study
Full size table

This research has considered six scenarios based on the severity of the earthquake and the time of its occurrence. As is clear, disaster happening in the rush hour causes more financial costs and human casualties than happening in a low-traffic and off-peak hour, also the capacity of hospitals and relief centers in high-traffic and rush hours are lesser due to the possibility of the occurrence of various events in the city. Given the large volume of input data of the problem, the values of some parameters, only for the first and second scenarios, are as follows: The numbers of homeless and injured people and corpses of each affected area are reported in Table 15. The numbers of required relief staff for each affected area are shown in Table 16. The capacity of each hospital for dispatching each kind of relief staff is provided in Table 17. The capacity of each hospital for acceptance of each kind of injuries is represented in Table 18.

Table 15 The number of homeless and injured people and corpses (person)
Full size table
Table 16 The number of required relief personnel (person)
Full size table
Table 17 The capacity of hospitals for dispatching relief personnel (person)
Full size table
Table 18 The capacity of hospitals for acceptance of injured people (person)
Full size table

7.1.1 The results of the case study

Table 19 is presented to better understand the distribution function estimation process. Also, the estimated value of the tent after 70 runs for 10,000 hours is shown in Table 19.

Table 19 Estimating the amount of required tent in scenario 1
Full size table

Next, the calculated mean distribution functions by Minitab 19.2020.1 software and the function fitting process are estimated. Figure 16 shows the simulation results. According to this figure, the distribution functions of demand including food, water, medicine and tents are estimated. The estimated distribution functions for food, water, medicine, and tent goods are normal, with correlation coefficients of 0.985, 0.990, 0.950, and 0.975, respectively.

Fig. 16
figure 16

Estimating the distribution functions

Full size image

After estimating the distribution functions of relief commodities, its mean and standard deviation are determined by the function fitting approach in Minitab 19.2020.1 software. Table 20 shows the simulation output. In this table, the average and standard deviation of relief commodities in each scenario are specified. For example, in the second scenario, the amount of required medicine follows a Normal distribution with an average of 156,208 kg and a standard deviation of 562.5.

Table 20 The results of simulation
Full size table

The validation of a simulation model indicates whether the model shows a true reflection of real-world performance. To prevent variability, the simulation model was run 100 times for 10,000 hours and the average performance of the model was reported. Figure 17 shows a comparison between the performance of the simulation model and the real-world. The results of the real-world are taken from the historical data of similar earthquakes in the region.

Fig. 17
figure 17

Comparison between results obtained by simulation and real system

Full size image

Figure 17 shows the validation of the simulation model. The vertical axis shows the estimated value of the relief commodities. It can be said, with a 95% confidence interval, that the simulation model has had a real estimate of the real-system.

In Table 21, there are 10 optimal solutions to illustrate the effective performance of the NSGA-II, SPEA-II, and PESA-II.

Table 21 The results of the suggested case study
Full size table

As it is clear that the quality of the solutions of the NSGA-II and the CPU time of this algorithm are better than the other two meta-heuristic algorithms. Figure 18 shows the results of CPU time for large-scale problems (case study) according to the three meta-heuristic algorithms.

Fig. 18
figure 18

The results of CPU time for large-scale problems of the presented meta-heuristic algorithms

Full size image

According to the stated arguments, the NSGA-II is better and more robust than the other two meta-heuristic algorithms. In the following, the case study will be solved with this algorithm.

Figure 19 shows the number of homeless people allocated to each shelter under various scenarios. The labels of each graph show the number of the affected area. As can be seen, the number of shelters established in the first and second scenarios is 20 and 24, for the third and fourth scenarios, 26 and 31, and for the fifth and sixth scenarios, 34. So, with the increase in severity of the earthquakes, the need for shelters increases as well. It is also clear that the number of people transported to shelters in scenarios 5 and 6 is higher than in other scenarios. Shelters number 16 and 18 are the most active and most important shelters as all of the scenarios transfer 276 and 251 homeless people to these two shelters, respectively.

Fig. 19
figure 19

Homeless people allocated to each shelter under each scenario

Full size image

Figure 20 shows how DCs are established based on their capacity in each scenario. In the first and second scenarios, 31 and 29 DCs are established, respectively. In the first scenario, 3, 17, and 11 DCs are established in small-, medium- and large-scale, respectively, and in the second scenario, 5, 14, and 9 DCs, in small-, medium- and large-scale, respectively. In the fifth scenario, 41 DCs are established, which include 23, 11, and 7 centers, in small-, medium- and large-scale, respectively. Also, in the sixth scenario, 43 ​​DCs are established, which include 24, 15, and 4 centers, in small-, medium- and large-scale, respectively. As you can see, the number of DCs established in scenarios with more severe earthquakes (scenarios 5 and 6) is far greater than their number for scenarios with less severe earthquakes (scenarios 1 and 2). Also, in scenarios with more severe earthquakes, there is a great tendency to establish more DCs with small capacity, and in scenarios with less severe earthquakes, there is a tendency to establish fewer DCs with large capacity.

Fig. 20
figure 20figure 20

The process of establishing DCs based on the capacity of DCs in each scenario

Full size image

7.1.2 Sensitivity analysis

According to Fig. 21, as the capacities of the DCs increase, the cost will increase, too. The increase in the size of the capacity causes not establishing the new DCs. Hence, the fixed cost of the establishment will be decreased. This reduction continues up to a certain point. As the capacities increases, the numbers of established DCs are reduced and transportation costs are increased.

Fig. 21
figure 21

Relationship between capacities of DCs and first objective function

Full size image

According to this figure, scenario 6 (earthquake of 7–8 magnitude during rush hour) will cost more than other scenarios. According to the explanations given in this scenario, there is a tendency to establish more DCs with small and medium capacities. The same behavior is also seen in scenario 5. In the first scenario (earthquake of 5–6 magnitude), there is a tendency to establish fewer DCs with large capacity. In these two scenarios, the cost slope is less than other scenarios.

The probability of existence of the routes and its impact on the third objective function is analyzed and shown in Fig. 22.

Fig. 22
figure 22

Relationship between probability of paths and third objective function

Full size image

According to Fig. 22, as the probability of the successful evacuation (the probability of routes being safe) increases, the value of third objective function (i.e., risk) decreases. Figure 22 shows the behavior of the third objective function for different scenarios when the probability of successful evacuation increases. Here, for all scenarios, three probabilities of successful evacuation are considered between 0 and 0.3, 0.3–0.6 and 0.6–0.9. In general, in all scenarios, with the increase of probability of successful evacuation, the value of the third objective function (risk of the routes) decreases. Also, the effect of increasing the probability of successful evacuation in reducing the value of the objective function in scenarios 5 and 6 is higher than other scenarios. This means that the effect of safe routes in more severe earthquakes is more than their effect in less severe earthquakes (scenarios 1 and 2).

8 Conclusion and future works

A new integrated stochastic multi-objective location-allocation-routing supply chain logistic programming model in disaster situations was proposed which has helped to design both pre-disaster and post-disaster plans. In order to integrate strategic and operational decisions, a new multi-objective mathematical model was presented for disaster relief management. The proposed model was used in a case study in two regions of Tehran Province in Iran. In this study, disaster relief decisions are considered at both strategic and operational levels for managers. Based on the decisions on the location of shelters, the location of DCs, and the allocation of capacity to them are stated as strategic decisions. Decisions on routing, the number of injured people, homeless people, and corpses who were transferred to the hospitals, shelters, and cemeteries, the number of commodities sent from the DCs to the shelters, and so on are selected as operational decisions. The simulation approach has been used to estimate the demand for relief commodities. In this way, the basic urban infrastructure has been identified and their interaction with each other during the earthquake has been calculated. The simulation results show that the demand for food, water, medicine, and tents follows the Normal distribution function with correlation coefficients of 0.985, 0.990, 0.950, and 0.975. Also, the validation results of the simulation model indicate that with a 95% confidence interval, it will have an accurate estimate of the real system. The estimated value is entered into the mathematical model and converted to a deterministic model by using the stochastic chance constraint programming approach.

The results of strategic and operational decisions extracted from the case study are presented in the following. Decisions on the number of established DCs and their capacity determination are of the strategic decisions. Also, decisions on how to transfer the homeless people from the affected areas to the shelters are considered as the operational decisions. The proposed model included three objectives as minimization of the expected value of the total costs of the relief supply chain, minimization of the maximum amount of unsatisfied demand for relief staff, minimization of the total probability of unsuccessful evacuation in routes. The destruction of facilities, infrastructures, and routes, also the amount of road traffic due to the disaster, and the severity of the disaster are considered. Four multi-objective algorithms, namely, NSGA-II, PESA-II, SPEA-II, and Epsilon-constraint methods were utilized to solve the suggested model. The outcomes of the solution in small and medium-sized problems indicate that the quality of the solutions and the average CPU time for the NSGA-II is better and more powerful than the other two meta-heuristic algorithms. Therefore, the case study has been solved using this algorithm.

A sensitivity analysis was performed to investigate the efficiency of the proposed model against the volatility of important parameters. The outcomes of the sensitivity analysis illustrated by increasing the severity of the earthquake (scenarios 5 and 6), more DCs with less capacity will be needed. This causes the significant spacing of DCs in the region. Also, less severe earthquakes (scenarios 1 and 2) required the establishment of fewer DCs with more capacity.

This research will help managers and decision-makers to simultaneously make strategic decisions before a disaster occurs and operational decisions during a disaster. Hence, determining the optimal capacity of established shelters and DCs helps managers to deal with possible shortages during the disaster. In addition, estimating the demand for relief commodities has always been one of the major concerns of disaster relief managers and decision-makers. Then, the proposed simulation approach and determining the demand distribution function of relief commodities such as water, food, tents, and medicine can increase the readiness of managers to cope with the inherent uncertainty of earthquakes.

According to the research results, the higher the magnitude of the earthquake, the higher the costs and loss of life, and the more relief commodities are needed. Moreover, managers are advised to be properly prepared for relief operations in such situations, especially during high traffic, which increases the severity of the damage. It is also essential for managers to know that as the magnitude of an earthquake increases, especially during high traffic, the capacity of established DCs increases dramatically.

There is also the possibility of destruction and blockage of roads during high-magnitude earthquakes. Therefore, managers are suggested to strengthen the main urban routes. The location of potential points for the establishment of shelters and DCs should also be given the strength of the paths leading to it. This can greatly reduce the likelihood of an unsuccessful evacuation.

Finally, in order to observe justice in the evacuation, fast triage operations are recommended to decision-makers. Depending on the different evacuation priorities for the injured, the homeless, and the dead, a quick triage can speed up the evacuation operation and prevent further casualties.

The following points are recommended for future research:

  • considering other phases of disaster relief management: for example, incorporating uncertainty in recovery operations and waste collection with disaster remnants,

  • considering a cooperative game for relief groups so that with their cooperation, the golden time of rescue and its costs are reduced.

  • using the other meta-heuristic algorithms such as multi-objective simulated annealing, multi-objective particle swarm optimization, etc.

The research limitations are as follows:

  • As there was no official database for some parts of cost elements, the expert’s estimations were asked to help. The questions about the cost for establishing distribution centers and shelters have been categorized and the estimated costs have been entered into the mathematical model.

  • Usually, numerous runs are required for each simulation model, and this can lead to high costs for using a computer. Simulation also requires access to a computer system equipped with features such as high RAM and CPU.