Two-Stage Multi-Objective Stochastic Model on Patient Transfer and Relief Distribution in Lockdown Area of COVID-19

Abstract

The outbreak of an epidemic disease may cause a large number of infections and a slightly higher death rate. In response to epidemic disease, both patient transfer and relief distribution are significant to reduce corresponding damage. This study proposes a two-stage multi-objective stochastic model (TMS-PTRD) considering pre-pandemic preparedness measures and post-pandemic relief operations. The proposed model considers the following four objectives: the total number of untreated infected patients, the total transfer time, the overall cost, and the equity distribution of relief supplies. Before an outbreak, the locations of temporary relief distribution centers (TRDCs) and the inventory levels of established TRDCs should be determined. After an outbreak, the locations of temporary hospitals (THs), the locations of designated hospitals (DHs), the transfer plans for patients, and the relief distribution should be determined. To solve the TMS-PTRD model, we address an improved preference-inspired co-evolutionary algorithm named the PICEA-g-AKNN algorithm, which is embedded with a novel similarity distance and three different tailored evolutionary strategies. A real-world case study of Hunan of China and 18 test instances are randomly generated to evaluate the TMS-PTRD model. The finding shows that the PICEA-g-AKNN algorithm is better than some most widely used multi-objective algorithms.

1. Introduction

From ancient times to the present, epidemics have been threatening the lives of people. Typical outbreaks in recent years include coronavirus disease 2019 (COVID-19) [1], H1N1 influenza [2], and Ebola [3,4]. In 2020, 80,000 people were infected with COVID-19 in Wuhan, the most severely affected city in China. To control the spread of the COVID-19 disease, the urban area of Wuhan was effectively blocked; that is, people and vehicles could not contact the outside world without permission.

In a lockdown area, emergency management must transfer infected patients to hospitals and provide relief supplies to ensure infected patients receive timely medical treatment. Because of the sudden outbreak of an epidemic, lockdown areas are often unable to access relief supplies in time. Given the time required to replenish relief supplies and set up new emergency facilities in affected areas, it is necessary to prepare relief supplies in advance to deal with the early stages of an outbreak. In order to predict the number of infected patients in the future, we can build a time-varying diffusion model. Although diffusion dynamics models have been applied to many studies on emergency management, it is unrealistic to predict with perfect accuracy an epidemic. Because the outbreak of an epidemic disease is highly unpredictable, one of the most important goals is to fully consider the uncertainty of the disaster.

For most disasters, providing relief supplies to injured people or infected patients is not enough. Transferring affected people from disaster areas to hospitals ensures they receive better treatment. During an epidemic outbreak, transferring patients is of great significance for maintaining the stability of affected areas and avoiding further outbreaks of infectious diseases. Although there is much research on relief supplies allocation, the interaction between patient transfer and relief distribution has not been well addressed in epidemic outbreaks. In this study, we not only transfer patients to hospitals but also provide these patients with relief supplies.

Because various stakeholders are concerned about emergency management, one of the major challenges of emergency management is to ensure fair emergency relief for all affected people. On the other hand, the unfair distribution of relief supplies may cause panic and psychological stress among the affected people. Furthermore, if the emergency operations only transfer patients without providing relief supplies, it may lead to the deterioration of infected people after transfer.

This paper focuses on emergency management during an epidemic outbreak to ensure that epidemic patients are treated and relief supplies are provided, considering the uncertainties and the control measures. In this study, coordination between patient transfer and relief distribution considers pre-pandemic preparedness measures and post-pandemic relief operations. At the same time, emergency relief plans provide reasonable decision support to decision makers in terms of equity, efficiency, and economics. The contributions of our study are as follows: first, the TMS-PTRD model is proposed to address the locations of TRDCs and the inventory levels of relief supplies in the preparedness state, and to address the locations of THs, patient transfer, and relief distribution for epidemic disasters; second, a multi-objective programming model for patient transfer and relief distribution problems is proposed to balance the efficiency, the fairness, and the economy of the emergency management system; finally, an improved multi-objective algorithm is proposed for large-scale TMS-PTRD problems, and the insights of emergency management are obtained.

The rest of this article is arranged as follows. The Section 2 reviews the literature and analyzes the gaps in the existing research. The Section 3 introduces the problem and constructs the TMS-PTRD model. The Section 4 designs an improved multi-objective algorithm. The epidemic simulations and numerical experiments are carried out in the Section 5. The last section expresses the conclusions and future research.

2. Literature Review

In order to reduce the damage caused by disasters, theoretical research and practical applications of emergency management have attracted more attention in recent years. Existing research can generally be divided into two categories: relief distribution and patient transfer.

An important measure of emergency management is relief distribution to the affected areas. Effective distribution of relief supplies (such as medicines, blood supplies, and emergency supplies) is a huge challenge for emergency logistics [5]. Considering population distribution and resource constraints, Duhamel et al. presented a location–allocation model to seek optimal solutions for post-disaster operations [6]. Moreno et al. [7] developed a multi-period, multi-commodity relief location-routing model and proposed a customized heuristics model based on time, stage, and scenario. In this research, vehicle reuse was considered to achieve coordination between different participants in emergency logistics. Because of the transportation system caused by an earthquake, a multi-objective, multi-period relief distribution model considering road repair was developed by Vahdani et al. [8]. Ghasemi et al. [9] considered several types of casualties and failure of emergency centers in an earthquake. The total cost and satisfaction of relief materials were considered in the disaster relief network. Elci et al. [10] considered a relief network problem to determine the capacities and locations of emergency facilities in an uncertain post-disaster environment. In general, emergency management consists of four main phases: mitigation, preparedness, response, and recovery [11]. Grass and Fischer [12] considered that two-stage programming is suitable for some emergency management because some decisions must be made before a disaster, and some decisions can only be made after the disaster. Ahmadi et al. [13] proposed a two-stage humanitarian location-routing model that takes into account random travel time to determine the locations of distribution centers and vehicle path planning. In view of its outstanding performance of two-stage planning in emergency preparedness and response, it has received more and more attention in recent years [14]. Vahdani et al. [15] studied a humanitarian logistics model to provide effective relief distribution supplies after an earthquake. In addition, in order to deal with the uncertain factors caused by the disaster, a robust optimization method was used. Considering primary and secondary disasters, Zhang et al. [16] proposed a multi-stage emergency logistics model to address the relief allocation problem. Mansour et al. [17] considered humanitarian logistics regarding pre-position and post-disaster. They presented a two-stage stochastic programming model for the disaster relief problem under such uncertain disaster scenarios.

Another important operation of disaster relief is to transfer the affected patients, which means transporting infected people from affected areas to shelters or hospitals [18,19]. Considering there is a large number of injured people in public emergencies, emergency managers need to formulate strategies for transferring injured people [20,21]. Flores et al. [22] presented a dynamic transfer model, in which affected people were classified. Goerigket [23] proposed a bus transfer model to move affected people to shelters during a disaster. Lim et al. [24] represented a capacity-constrained network model for the short-notice transfer problem. A greedy heuristic model was improved to optimize maximizing total weighted patient outflow. Shahparvari and Abbasi [25] developed an emergency management model to optimize transfer vehicle scheduling in unstable road scenarios after a disaster. Shahparvari et al. [26] developed a short-notice bushfire transfer model to ensure that the maximum number of critical patients can be safely transferred to shelters. As difficult as the relief distribution, patient transfer needs to deal with multiple uncertainties about the affected areas, time, and the scale of a disaster. Bayram and Yaman [27] studied the patient transfer problem considering the uncertainty caused by the destruction of road networks and shelters. Liang et al. [28] considered that transferring victims to shelters during disasters can keep them safe. Considering the uncertainty in the demands of disasters, they proposed a location and evacuation model based on minimizing total transfer time.

Because patient transfer and relief distribution are considered separately, a single emergency operation has a poor result for emergency management [29]. At present, few articles have considered the importance of combining patient transfer with relief distribution for emergency management. Sabah et al. [30] proposed an integrated affected people transfer and relief distribution model to optimize vehicle scheduling. However, the number of the affected people and the relief demand were considered certain in this research. Generally, accurate information is almost impossible to obtain in a disaster. Setiawan et al. [31] proposed a series of humanitarian logistics models to deal with relief distribution and victim transfer problems after a disaster. This research considered coordination among rescue strategies but did not take into account uncertain parameters.

Emergency management should consider not only costs and benefits but also humanitarian issues [32,33]. In addition, non-profit indicators reflect the humanistic care of emergency management [34]. Fereiduni and Shahanaghi [35] presented an emergency network to make optimal choices about location, allocation, and transfer simultaneously. However, the study did not consider humanitarian issues because only the cost objective was considered. Vahdani et al. [36] proposed an emergency location-routing problem to optimize relief supplies delivery and patient transfer. In this study, a fuzzy credibility theory was used to formulate the fuzzy demand of the patients. However, the study only focuses on the setup cost for ELCs, transportation costs, and the fixed cost for vehicles.

Some researchers have taken into account the uncertainty parameters and humanitarian issues in emergency management. Furthermore, a few studies have not obtained an excellent Pareto solution or designed a suitable algorithm for large-scale cases. Considering the demand of different refugees, Wang et al. developed a location–allocation model to optimize total transfer distance, total cost, and total unmet resources. However, this research converted multiple objectives into one objective using a weighting method and did not obtain the actual Pareto solution [37]. Ghasemi et al. [38] developed a victim transfer and relief distribution model to optimize three objectives: the untreated victims, the satisfaction with emergency supplies, and the total cost. However, the research did not consider transfer time, which is very important in emergency management. Mohammadi et al. [39] proposed a multi-objective location–allocation model to optimize victim transfer and relief distribution after a disaster. This study attempts to decide on setting up relief suppliers, selecting locations of relief distribution centers and emergency centers, and transferring injured people and distributing relief.

Based on this review, the conclusions of relevant studies can be drawn as follows:

• (1). Many studies have considered victim transfer and relief distribution separately, but few studies have considered the coordination between the two operations.

• (2). Some studies of emergency management have considered humanitarian issues, but few studies have been able to obtain a real Pareto optimal solution set.

• (3). No study has considered the characteristics of epidemic diseases and the control measures in a multi-objective victim transfer and relief distribution problem.

Therefore, we propose a two-stage multi-objective stochastic model (TMS-PTRD) considering patient transfer and relief distribution to deal with an epidemic disease. In addition, emergency management is a typical NP-hard problem involving the location of emergency facilities, personnel transfer, material distribution, and vehicle routing. So, it poses a great challenge for multi-objective optimization algorithms. Exact algorithms have been shown to solve small-scale emergency management problems. However, evolutionary algorithms such as NSGA-II [40] and MOPSO [41] have been usually used to solve large-scale multi-objective emergency management problems. When solving problems with more than three objectives, some classical multi-objective optimization algorithms are likely to have poor performance [42]. In order to solve the problem with more objectives, researchers have developed some many-objective algorithms such as MOEA/D [43] and HypE [44]. The principle of the above many-objective algorithms is based on the relationship between multiple objectives rather than the dominant relationship between solutions. A novel concept of decision-maker preference co-evolution has been designed for the many-objective problem (excess of three objectives) [45]. Based on the concept, a preference-inspired co-evolutionary algorithm (PICEA-g) was developed by Wang et al. [46,47]. Previous studies have shown that one heuristic algorithm is usually not suitable for different models [7]. Hence, it is important to design a tailored algorithm for the TMS-PTRD model proposed in this paper.

3. Problem Definition and Formulation

The TMS-PTRD model addresses patient transfer and relief distribution during an outbreak. We consider the problem (see Figure 1) as an integrated emergency network, consisting of TRDCs, fever clinics, THs, and DHs. In the preparation phase, the locations of the TRDCs and the inventory levels of relief supplies are determined. After the outbreak, emergency decisions include the locations of the THs, the locations of the DHs, the relief distribution for the fever clinics, and the transfer plan for different types of infected people.

In Figure 1, black arrows represent relief distribution. The blue arrows show the flows of infected people, and the red arrows show the flows of critical patients. The main assumptions are as follows:

• (1). Each community has a fever clinic, in which the general physicians in the community diagnose illnesses and transfer confirmed cases to the appropriate hospitals;

• (2). The confirmed cases consisting of the infected people and critical patients can be estimated by the modified SEIR model;

• (3). We assume that untransferred confirmed cases do not affect the forecast results, because only a small number of confirmed cases are not transferred in time and will be preferentially transferred in the next period;

• (4). The infected people need to be transferred to the THs and the critical patients need to be transferred to the DHs;

• (5). The fever clinics provide a set of medical and ancillary supplies (such as an N95 mask and an additional protective suit) for each infected patient and critical patient;

• (6). Each established TRDC can distribute relief supplies to each fever clinic using a vehicle with a homogeneous capacity.

Note that the planning horizon is a period of time (such as a week) at the beginning of the epidemic, after which the affected areas can get external relief supplies and build more emergency facilities. How measures are implemented goes beyond our research scope.

3.1. Time-Varying Epidemic Prediction Model

The relief demand depends on the number of confirmed cases during an epidemic. When making emergency management plans, it is necessary to assess the development of the epidemic in advance. The SEIR models were applicable to epidemic diffusion research [48,49]. Considering epidemic outbreaks are usually short-lived, we assume that the population in the affected areas will not change significantly. In the SEIR model, we assume that the recovered individuals receive lifelong immunity. Hence, ordinary differential equations can be used to describe the epidemic spread as follows.

(1)$\{\begin{array}{c}\raisebox{1ex}{$dS\left(t\right)$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=-r\beta I\left(t\right)S\left(t\right)\hfill \\ \raisebox{1ex}{$dE\left(t\right)$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=r\beta I\left(t\right)S\left(t\right)-\omega E\left(t\right)\hfill \\ \raisebox{1ex}{$dI\left(t\right)$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=\omega E\left(t\right)-\gamma I\left(t\right)\hfill \\ \raisebox{1ex}{$dR\left(t\right)$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=\gamma I\left(t\right)\hfill \end{array}$

where N is the total population of the infected areas, S(t) is the number of susceptible people in an outbreak, E(t) is the number of exposed people in contact with the infection, I(t) is the number of infected people with all symptoms, R(t) is the number of recovered people, $r$ is the exposure ratio, $\beta$ is the transmission probability, $\omega$ is the morbidity of patients in the latent period, and $\gamma$ is the recovery rate.

For many diseases such as SARS and COVID-19, patients in the latent period can also spread to susceptible people. In practice, transferring confirmed cases to hospitals can control these patients from spreading the disease. Considering the interventions of epidemic diseases, the susceptible people in isolation (SI), the exposed people in isolation (EI), and the critical patients (H) in urgent need of hospitalization are added to the modified model. Based on the above point of view, we use a modified SEIR model to evaluate the epidemic [50].

(2)$\{\begin{array}{c}\raisebox{1ex}{$\mathrm{d}S\left(t\right)$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=-\left[r\beta +rq\left(1-\beta \right)\right]S\left(t\right)\left[I\left(t\right)+\theta E\left(t\right)\right]+\lambda S_q\left(t\right)\hfill \\ \raisebox{1ex}{$\mathrm{d}E\left(t\right)$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=r\beta \left(1-q\right)S\left(t\right)\left[I\left(t\right)+\theta E\left(t\right)\right]-\omega E\left(t\right)\hfill \\ \raisebox{1ex}{$\mathrm{d}I\left(t\right)$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=\omega E\left(t\right)-\left({\delta }_I+\kappa +{\gamma }_I\right)I\left(t\right)\hfill \\ \raisebox{1ex}{$d\mathrm{S}\mathrm{I}\left(t\right)$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=rq\left(1-\beta \right)S\left(t\right)\left[I\left(t\right)+\theta E\left(t\right)\right]-\lambda \mathrm{S}\mathrm{I}\left(t\right)\hfill \\ \raisebox{1ex}{$d\mathrm{E}\mathrm{I}\left(t\right)$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=r\beta qS\left(t\right)\left[I\left(t\right)+\theta E\left(t\right)\right]-{\delta }_q\mathrm{E}\mathrm{I}\left(t\right)\hfill \\ \raisebox{1ex}{$\mathrm{d}H\left(t\right)$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.={\delta }_{I}I\left(t\right)+{\delta }_{q}E\left(t\right)-\left(\kappa +{\gamma }_H\right)H\left(t\right)\hfill \\ \raisebox{1ex}{$dR\left(t\right)$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.={\gamma }_{I}I\left(t\right)+{\gamma }_{H}H\left(t\right)\hfill \end{array}$

where $q$ is the isolation ratio, $\theta$ is the transmission coefficient in the latent period, $\lambda$ is the rate of release from quarantine, $\kappa$ is the death rate, ${\delta }_q$ is the morbidity of quarantined patients in the latent period, ${\delta }_I$ is the probability that infected people will be isolated for treatment, ${\gamma }_H$ is the recovery rate of hospitalized patients, and ${\gamma }_I$ is the recovery rate of infected people.

Given the initial values and related parameters, the modified SEIR model can be used to assess the time series data of patients. We define the relief supplies required by each patient in a day as θ. The infected people need to be transferred to the THs for isolation and treatment, and the critical patients need to be transferred to the DHs for specialty treatment [50]. Before reaching the hospitals, relief supplies must be provided to the infected people and critical patients to ensure that the disease is not spread during transfer. So, the demand for relief supplies is formulated as:

(3)$d=\theta \left[I\left(t\right)+H\left(t\right)\right]$

3.2. Notation

Sets and Index:

$S$:	Set of scenarios, $s\in S$;
$T$:	Set of planning periods, $t\in T$;
$I$:	Set of TRDCs, $i\in I$;
$C$:	Set of fever clinics, $\mathrm{c}\in C$;
$J$:	Set of THs, $j\in J$;
$H$:	Set of DHs, $h\in H$.

Model Parameters:

$P_s$:	Probability of scenario s;
$t_{nm}^p$:	Transportation time from the node n to the node m;
$F{c}_i$:	Fixed cost of the TRDC i;
$F{t}_j$:	Fixed cost of the TH j;
$F{d}_h$:	Fixed cost of the DH h;
$c_1$:	Holding cost of a relief supply;
$c_2$:	Penalty cost of a unit of relief supply;
$C_{cj}^M$:	Transfer cost for the infected people from the fever clinic c to the TH j;
$C_{ch}^S$:	Transfer cost for a critical patient from the fever clinic c to the DT h;
$C{V}_{ic}^R$:	Transport fees for one vehicle from the TRDC i to the fever clinic c;
$C{R}_i$:	Maximum capacity of the TRDC i;
$Cav$:	Capacity of a vehicle to carry the relief supplies;
$C{m}_j$:	Maximum capacity of the TH j;
$C{s}_h$:	Maximum capacity of the DH h;
$\tilde{I}{p}_c^{ts}$:	Number of infected people at the fever clinic c in scenario s;
$\tilde{C}{p}_c^{ts}$:	Number of critical patients at the fever clinic c in scenario s;
${\varphi }^S$:	Priority of satisfying the critical patient;
$\alpha$:	Confidence levels of the chance-constrained model;
$M_{big}$:	A large positive number.

Decision Variables:

$Y_i$:	1, if the TRDC i is established, 0 else;
$Z_j^s$:	1, if the TH j is established under scenario s, 0 else;
$K_h^s$:	1, if the designated hospital h be established under scenario s, 0 else;
$V{N}_i$:	Amount of inventory at the TRDC i;
$X_{ic}^{ts}$:	Number of relief supplies transported from the TRDC i to the fever clinic c in scenario s.

3.3. Model Formulation

Equation (4) minimizes the number of patients not transferred at different periods. Equation (5) minimizes the transfer time. Equation (6) minimizes the fixed cost of the TRDCs and pre-stored relief supplies in the first phase. In the second phase, the objective function minimizes the fixed THs, the fixed DHs, the cost of distributing relief supplies, and the cost of transferring patients. Equation (7) minimizes the dissatisfaction with relief service. Equations (8)–(10) indicate only one emergency facility can be established at each point. Equation (11) describes that only established TRDCs can be applied. Equation (12) describes the amount of unused relief supply. Equation (13) describes the carrying capacity of the transportation vehicles. Equation (14) describes the capacity of established THs. Equation (15) describes the capacity of the designated hospitals. Equation (16) defines the unserved confirmed patients in the fever clinics. Equation (17) defines a satisfaction function for the fever clinics. Equations (18) and (19) describe the binary variables and the non-negative variables.

(4)$\min f_1={\displaystyle \sum _{s\in S}{\displaystyle \sum _{t\in T}{\displaystyle \sum _{c\in C}{P}_s\cdot \left[\left({\displaystyle \sum _{\tau =1}^t\tilde{I}{p}_c^{\tau s}}-{\displaystyle \sum _{\tau =1}^t{\displaystyle \sum _{j\in J}P{M}_{cj}^{\tau s}}}\right)+{\varphi }^S\left({\displaystyle \sum _{\tau =1}^t\tilde{C}{p}_c^{\tau s}}-{\displaystyle \sum _{\tau =1}^t{\displaystyle \sum _{h\in H}P{S}_{ch}^{\tau s}}}\right)\right]}}}$

(5)$\min f_2={\displaystyle \sum _{s\in S}{P}_s}\cdot \left({\displaystyle \sum _{t\in T}{\displaystyle \sum _{c\in C}{\displaystyle \sum _{j\in J}{t}_{cj}P{M}_{cj}^{ts}}+{\varphi }^S{\displaystyle \sum _{t\in T}{\displaystyle \sum _{c\in C}{\displaystyle \sum _{h\in H}{t}_{ch}P{S}_{ch}^{ts}}}}}}\right)$

(6)$\begin{array}{l}\min f_3={\displaystyle \sum _{i\in I}F{c}_i\cdot Y_i+{\displaystyle \sum _{i\in I}{c}_1\cdot V{N}_i}}+{\displaystyle \sum _{s\in S}{P}_s}\cdot \left({\displaystyle \sum _{j\in J}F{t}_j\cdot Z_j^s}+{\displaystyle \sum _{h\in H}F{d}_h\cdot K_h^s}+{\displaystyle \sum _{t\in T}{\displaystyle \sum _{i\in I}{\displaystyle \sum _{c\in C}C{V}_{ic}^{R}V{R}_{ic}^{ts}}}}+\right.{\displaystyle \sum _{t\in T}{\displaystyle \sum _{c\in C}{\displaystyle \sum _{j\in J}{C}_{cj}^{M}P{M}_{cj}^{ts}}+}}\\ \left.{\displaystyle \sum _{t\in T}{\displaystyle \sum _{c\in C}{\displaystyle \sum _{h\in H}{C}_{ch}^{S}P{S}_{ch}^{ts}+{\displaystyle \sum _{i\in I}{c}_2\cdot UV{N}_i^s}}}}\right)\end{array}$

(7)$\min f_4={\displaystyle \sum _{s\in S}{\displaystyle \sum _{t\in T}{\displaystyle \sum _{c\in C}{P}_s\cdot \phi \left(U_c^{ts}\right)}}}$

(8)$Y_i\le 1,\forall i\in I$

(9)$Z_j^s\le 1,\forall j\in J$

(10)$K_h^s\le 1,\forall h\in H$

(11)$V{N}_i\le C{R}_i\cdot Y_i,\forall i\in I$

(12)$V{N}_i-{\displaystyle \sum _{t\in T}{\displaystyle \sum _{c\in C}{X}_{ic}^{ts}}}\le UV{N}_i^s,\forall i\in I,s\in S$

(13)$X_{ic}^{ts}\le cav\cdot V{R}_{ic}^{ts},\forall c\in C,i\in I,t\in T,s\in S$

(14)${\displaystyle \sum _{t\in T}{\displaystyle \sum _{c\in C}P{M}_{cj}^{ts}}}\le C{m}_j{Z}_j^s,\forall j\in J,s\in S$

(15)${\displaystyle \sum _{t\in T}{\displaystyle \sum _{c\in C}P{S}_{ch}^{ts}}}\le C{s}_h{K}_h^s,\forall h\in H,s\in S$

(16)$U_c^{ts}\le 1-\frac{{\displaystyle \sum _{i\in I}{X}_{ic}^{ts}}}{\tilde{I}{p}_c^{ts}+\tilde{C}{p}_c^{ts}},\forall c\in C,t\in T,s\in S$

(17)$\phi \left(U_c^{ts}\right)=e^{\frac{-U_c^{ts}}{1-U_c^{ts}}},\forall c\in C,t\in T,s\in S$

(18)$Y_i,Z_j^s,K_h^s\in \left\{0,1\right\},\forall i\in I,j\in J,h\in H,s\in S$

(19)$V{N}_i,X_{ic}^{ts},V{R}_{ic}^{ts},UV{N}_i^s,P{M}_{cj}^{ts},P{S}_{ch}^{ts}\ge 0\&Integer,\forall c\in C,i\in I,j\in J,h\in H,t\in T,s\in S$

3.4. A Chance-Constrained Model for Patient Transfer and Relief Distribution

Considering the results predicted by the modified SEIR model will partially deviate from the actual data, this may have stochastic characteristics in the two-stage patient transfer and relief distribution problem. The stochastic model is provided with a predetermined reliability level (α) by using chance constraints. So, we used the chance-constraint model to define reliable sets as follows [51].

(20)$\underset{x\in X}{\min }{f}_k\left(x\right)={\displaystyle \sum _{j=1}^n{\tilde{c}}_{kj}}{x}_j,k=1,2,...,K$

(21)$\mathrm{Pr}\left\{{\displaystyle \sum _{j=1}^N{a}_{ij}{x}_j\ge {\tilde{b}}_i}\right\}\ge {\alpha }_i,i=1,2,...,M$

(22)$x\in X$

where $X$ is the deterministic feasible region, f(x) is the value of the minimized function, and α is the confidence level of the chance-constrained optimization. According to the inverse function theorem, the chance constraint is rewritten are as follows:

(23)$\underset{x\in X}{\min }{f}_k\left(x\right)=E\left({\displaystyle \sum _{j=1}^N{\tilde{c}}_{kj}}{x}_j\ge {\tilde{b}}_i\right),k=1,2,...,K;i=1,2,...,M$

(24)$\mathrm{Pr}\left\{{\displaystyle \sum _{i=1}^I{\tilde{a}}_{ij}{x}_i\ge {\tilde{b}}_i}\right\}\ge {\alpha }_i,i=1,2,...,M$

(25)$x=\left(x_1,x_2,...,x_n\right)$

(26)$x\in X$

Based on the above theory, the summary can be rewritten as follows.

(27)$E\left({\displaystyle \sum _{j=1}^n{c}_{kj}^{\ast }{x}_j}-f_k^{-}\right)-{\phi }^{-1}\left({\alpha }_k\right)\sqrt{Var}\left({\displaystyle \sum _{j=1}^n{c}_{kj}^{\ast }{x}_j}-f_k^{-}\right)\ge 0,k=1,2,...,K$

(28)$E\left({\displaystyle \sum _{j=1}^n{c}_{kj}^{\ast }{x}_j}-f_k^{+}\right)+{\phi }^{-1}\left({\alpha }_k\right)\sqrt{Var}\left({\displaystyle \sum _{j=1}^n{c}_{kj}^{\ast }{x}_j}-f_k^{+}\right)\le 0,k=1,2,...,K$

where $f_k^{-}=\min c_{kj}^{\ast }{x}_j$ and $f_k^{+}=\max c_{kj}^{\ast }{x}_j$.

(29)$E\left({\displaystyle \sum _{j=1}^n{\tilde{a}}_{ij}{x}_j}-{\tilde{b}}_i\right)-{\phi }^{-1}\left(1-{\alpha }_i\right)\sqrt{Var}\left({\displaystyle \sum _{j=1}^n{\tilde{a}}_{ij}{x}_j}-{\tilde{b}}_i\right)\ge 0,i=1,2,...,M$

Making use of the auxiliary variables $W_{stc}^1$ and $W_{stc}^2$, we rewrite the first objective function as a deterministic model at the confidence level $\alpha$ as follows:

(30)$\min f_1={\displaystyle \sum _{s\in S}{\displaystyle \sum _{t\in T}{\displaystyle \sum _{c\in C}{P}_s\cdot \left(W_{stc}^1+{\varphi }^S{W}_{stc}^2\right)}}}$

(31)$E\left({\displaystyle \sum _{\tau =1}^t\overline{I}{p}_c^{ts}}\right)+{\phi }^{-1}\left(1-{\alpha }_j\right)\sqrt{Var}\left({\displaystyle \sum _{\tau =1}^t\tilde{I}{p}_c^{ts}}\right)-{\displaystyle \sum _{\tau =1}^t{\displaystyle \sum _{j\in J}P{M}_{cj}^{\tau s}}}\le W_{stc}^1$

(32)$E\left({\displaystyle \sum _{\tau =1}^t\overline{C}{p}_c^{\tau s}}\right)+{\phi }^{-1}\left(1-{\alpha }_j\right)\sqrt{Var}\left({\displaystyle \sum _{\tau =1}^t\tilde{C}{p}_c^{\tau s}}\right)-{\displaystyle \sum _{\tau =1}^t{\displaystyle \sum _{h\in H}P{S}_{ch}^{\tau s}}}\le W_{stc}^2$

(33)$W_{stc}^1,W_{stc}^2\ge 0$

Similarly, Equation (16) is modified as follows:

(34)$U_c^{ts}\le 1-\frac{{\displaystyle \sum _{i\in I}X{S}_{ic}^{ts}}}{\left(\overline{I}{p}_c^{ts}+\overline{C}{p}_c^{ts}\right)+{\phi }^{-1}\left(1-{\alpha }_j\right)\sqrt{Var}\left(\tilde{I}{p}_c^{ts}+\tilde{C}{p}_c^{ts}\right)},\forall c\in C,t\in T,s\in S$

Finally, the chance-constrained model for the TMS-PTRD problem is represented as follows.

Min. (30), (5), (6), and (7)

S.t. (8)~(15), (17)~(19), (31)~(34).

4. Solution Algorithm

4.1. ε-Constraint Algorithm

In the ε-constraint algorithm, the highest priority objective is selected based on the preferences of the decision maker, and other objectives are converted to constraints. The Pareto front is obtained by changing the upper and lower boundaries of these constraints [52]. The ε-constraint algorithm is as follows:

(35)$\min f_1\left(X\right)$

S.t.

(36)$f_2\left(X\right)\le {\epsilon }_2$

(37)$f_m\left(X\right)\le {\epsilon }_m$

where X is a feasible set. We use the epsilon parameters to show the range of each objective function. The upper and lower bounds of the objective function are shown in Equations (38) and (39).

(38)${\overline{\epsilon }}_i=\min \left\{f_i\left(x\right)\left|x\in X\right.\right\}$

(39)${\underset{¯}{\epsilon }}_i=\max \left\{f_i\left(x_1^{\ast }\right),f_i\left(x_2^{\ast }\right),...f_i\left(x_m^{\ast }\right)\right\}$

We use Equation (40) to calculate the ε parameters.

(40)${\epsilon }_{i{j}_i}={\underset{¯}{\epsilon }}_i-\left(\frac{{\underset{¯}{\epsilon }}_i-{\overline{\epsilon }}_i}{q_i}\right)\cdot j_i,j_i=0,1,...,q_i$

where $q_i$ represents the number of intervals between the upper and lower bounds. In this algorithm, a payoff constrain table of the ε-constraint algorithm is as follows:

• Step 1:. Calculate a payoff table

  $\left(\begin{array}{c}\begin{array}{cccc}{f}_2\left(x_1^{\ast }\right)& f_3\left(x_1^{\ast }\right)& ...& f_m\left(x_1^{\ast }\right)\end{array}\\ ...\\ \begin{array}{cccc}{f}_2\left(x_2^{\ast }\right)& f_3\left(x_3^{\ast }\right)& ...& f_m\left(x_m^{\ast }\right)\end{array}\end{array}\right)$

  $x_i^{\ast }={\arg }_{x\in X}\min f_i\left(x\right),i=1,2,...,m$

• Step 2:. Set ${\overline{\epsilon }}_i=f_i\left(x_i^{\ast }\right),{\underset{¯}{\epsilon }}_i=f_i\left(x_1^{\ast }\right),i=1,2,...,m$ a

• Step 3:. $P=\left\{x_1^{\ast },x_2^{\ast },...,x_m^{\ast }\right\},F=\left\{f_1\left(x_1^{\ast }\right),f_2\left(x_2^{\ast }\right),...,f_m\left(x_m^{\ast }\right)\right\}$

• Step 4:. Solve $x^{\ast }=opt\left(f_1,{\epsilon }_{2j_2},{\epsilon }_{3j_3},...,{\epsilon }_{m{j}_m}\right)$

  $$\begin{gathered} f_2\left(X\right)\le {\epsilon }_{2j_2}\mathrm{for}{j}_2=0:1:q_2 \\ {f}_3\left(X\right)\le {\epsilon }_{3j_3}\mathrm{for}{j}_3=0:1:q_3 \\ \dots \\ {f}_m\left(X\right)\le {\epsilon }_{m{j}_m}\mathrm{for}{j}_m=0:1:q_m \end{gathered}$$

• Step 5:. $P=P\cup \left\{x^{\ast }\right\},F=F\cup \left\{f_1\left(x^{\ast }\right),f_2\left(x^{\ast }\right),...,f_m\left(x^{\ast }\right)\right\}$

• Step 6:. Return P and F

Given the priority of the objective functions, the ε-constraint algorithm can obtain the Pareto front of multi-objective problems.

4.2. An Improved PICEA-g Algorithm

Evolutionary algorithms have been applied to solve multi-objective problems. Numerous studies have shown that utilizing neighborhood information can make significant progress in algorithm research [53,54]. Considering it is very difficult to define the neighborhood of the TMS-PTRD model using Euclidean distance, a novel similarity distance is designed to construct the K-nearest neighborhood. The PICEA-g-AKNN algorithm is proposed by integrating the presented K-nearest neighborhood into the framework of the preference-inspired co-evolutionary algorithm (PICEA-g). At the same time, we use an evaluation method according to the fitness of the solutions and assign a tailored evolutionary strategy for each type of chromosome to improve the efficiency of the algorithm.

4.2.1. An Adaptive K-Nearest Neighborhood Method Based on a Novel Similarity Distance

We introduce an intuitive definition for similarity distance to determine the neighborhood for each chromosome. The two types of similarities between chromosomes are defined as follows:

• The similarity of the locations of emergency facilities

This type of similarity between $X$ and $X^{\prime }$ are measured by the location of the emergency facilities. Considering two factors (locations and size), the similarity of the emergency facilities is formulated as:

(41)$S{D}_1(X,X^{\prime })=\frac{1}{\left|I\right|\left|L\right|}{\displaystyle \sum _{i\in I}{\displaystyle \sum _{l\in L}{\phi }_1\left(x_{il},x_{il}^{\prime }\right)}}$

(42)${\phi }_1\left(x_i,x_i^{\prime }\right)=\left\{\begin{array}{c}1,if{x}_{il}=x_{il}^{\prime }\\ 0,otherwise\end{array}\right.$

where $x_{il}$ is one if emergency facility I is established; otherwise, $x_{il}$ is zero. The example in Figure 2 illustrates the encoding process for emergency facility locations.

In Figure 2, Chromosome (a) chooses to establish the second TRDC at Size 1 and the fourth TRDC at Size 3. Chromosome (b) chooses to establish the second TRDC at Size 1, the third TRDC at Size 2, and the fourth TRDC at Size 2.

• 2. The similarity between distribution planning and transfer planning

The similarity of the patient flows between solution $\mathrm{X}$ and solution $X^{\prime }$ can be measured by the number of patients transferred. Similarly, the similarity of the relief flows is measured by the transportation quantity of relief supplies.

(43)$S{D}_2(X,X^{\prime })={\displaystyle \sum _{i\in N}{\displaystyle \sum _{j\in M}\frac{1}{\left|N\right|\left|M\right|}}{\phi }_2\left({\alpha }_{stij},{\alpha }_{stij}^{\prime }\right)}$

At the decision period $t$, ${\alpha }_{stij}$, and ${\alpha }_{stij}^{\prime }$ are the number of relief supplies transported from the TRDC $i$ to the clinic $j$ in scenario $s$.

(44)${\phi }_2\left({\alpha }_{stij},{\alpha }_{stij}^{\prime }\right)=\left\{\begin{array}{c}{\alpha }_{stij}/{\alpha }_{stij}^{\prime },if{\alpha }_{stij}^{\prime }>{\alpha }_{stij}\\ {\alpha }_{stij}^{\prime }/{\alpha }_{stij},if{\alpha }_{stij}>{\alpha }_{stij}^{\prime }\\ 1,if{\alpha }_{stij}={\alpha }_{stij}^{\prime }\end{array}\right.$

The example in Figure 3 illustrates the encoding process of the patient flows and the relief flows.

From Figure 3, 100 relief supplies are transported from the second TRDC to Clinic 2, and 300 relief supplies are transported from the second TRDC to Clinic 4 in Chromosome (a). Similarly, 300 relief supplies are transported to Clinic 2 and Clinic 3.

Considering the two-stage stochastic model, we construct the similarity distance between two chromosomes as follows:

(45)$SD(X,X^{\prime })=S{D}_1(X,X^{\prime })+\frac{1}{\left|S\right|\left|T\right|}{\displaystyle \sum _{s\in S}{\displaystyle \sum _{t\in T}S{D}_2(X,X^{\prime })}}$

The K-nearest neighborhood has been proven useful in the field of computer science and engineering [55,56]. Generally speaking, a smaller neighborhood size can enhance local search ability, while a larger neighborhood size can improve global search ability [57]. So, an adaptive K-nearest neighborhood method is used to dynamically adjust the neighborhood size to balance exploitation ability and exploration ability. By evaluating the evolutionary state of the solutions, the adaptive K-nearest method decreases the neighborhood radius k to enhance the exploitation ability when the solution is in a continuous evolution state and increases the neighborhood radius k to improve the exploration ability when the solution is in a state of stagnant evolution. The adaptive K-nearest construction method is as Algorithm 1.

Algorithm 1: Framework of the K-nearest neighborhood method

Input: solution S, neighborhood size K1 and K2

1. set the current neighborhood size Kc and K1;

2.for$gen=1:genMax$;

3.          Find non-dominated solutions NS in S;

4.            Generate offspring solutions OS based on neighborhood size Kc;

5.            if$\mathrm{mod}(gen,5)==0$

6.                  Find non-dominated solutions NOS in OS;

7.                  if$x^{\prime }\prec x,x\in S,x^{\prime }\in OS$

8.                      $K_c=K_c-1$;

9.                  else

10.                    $K_c=K_c+1$;

11.                  end if

12.              end if

13.              if$K_c>K_1$

14.                $K_c=K_1$;

15.                else

16.                $K_c=K_2$

17.              end if

18. end for

4.2.2. The PICEA-g-AKNN Algorithm

In the presented algorithm, all solutions are divided into three evolutionary states by an assessment method. In addition, each state is assigned a customized evolutionary strategy. The assessment method and crossover operation for individuals are as follows.

• 1. We define a dominant individual if the solution dominates all of the neighborhoods. It is probable that the dominant individual is close to the Pareto front. Therefore, an SBX local search strategy is used to modify this individual. The SBX local search strategy not only enables the dominant individual to move close to the Pareto fronts but also does not cause a large disturbance for the outstanding individual. Let $Pc$ is the probability of executing SBX and the SBX is described as Equation (46).

  (46)$c{x}_{il}=\{\begin{array}{c}0.5\cdot \left[\left(1+\beta \right)x_{il}+\left(1-\beta \right)x_{jl}\right],ifrc\le Pc\\ 0.5\cdot \left[\left(1+\beta \right)x_{il}-\left(1-\beta \right)x_{jl}\right],otherwise\end{array}$

  where $rc$ is a randomly generated number in the range [0, 1]. Note that if $c{x}_{il}\le 0$, we assign that $c{x}_{il}$ is equal to zero. If $x_{il}=x_{jl}$, then assign $c{x}_{il}$ to $x_{il}$, where $c{x}_{il}$ is the offspring $x_i$ in the $l$ dimension, $x_j$ is a parent chromosome selected at random from the neighborhood, and $\beta$ is a random number generated as follows according to [58].

  (47)$\beta =\{\begin{array}{c}{\left(2\cdot rand\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$1+\eta $}\right.},ifrand\le 0.5\\ {\left[1/\left(2-2\cdot rand\right)\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$1+\eta $}\right.},otherwise\end{array}$

  where $\eta$ is a parameter that represents the degree of learning from the parent individual.

• 2. A solution is defined as an exploring individual if it is not dominated by any neighborhood and there are other non-dominant individuals in the neighborhood. The exploring individual is likely to obtain useful information from the non-dominant individual in the neighborhood. We use a classic and robust DE operator named “DE/rand/1” to generate offspring. Because more parents can be referenced in the DE, children individuals can exchange excellent information with individuals in the neighborhood [59]:

  (48)$c{x}_{il}=x_{il}+F\cdot \left(x_{r1,l}-x_{r2,l}\right)$

  where $c{x}_i$ is the children individual of $x_i$, and $x_{r1}$ and $x_{r2}$ are randomly selected from the non-dominant individuals in the neighborhood. When there is only one non-dominant solution of $x_i$ in the neighborhood, the other parent individual is selected randomly from the other individuals in the neighborhood. $F\in \left(0,2\right)$ is the mutation factor. We assign $c{x}_{il}$ is equal to zero if $c{x}_{il}\le 0$, and assign $c{x}_{il}$ to $x_{il}$ if $x_{il}=x_{jl}$.

• 3. A solution is defined as a learning individual if it is dominated by individuals in the neighborhood. For the learning individual, it is necessary to try to learn from the outstanding individuals in the neighborhood. A directional search DE named “DE/current-to-dominance” has good global search ability and is used to generate offspring [60].

  (49)$c{x}_{il}=x_{il}+K\cdot \left(x_{dl}-x_{il}\right)+F\cdot \left(x_{r1,l}-x_{r2,l}\right)$

  where $x_d$ is an individual that dominates the current individual in the neighborhood, $K\in \left(0,1\right)$ is another mutation factor, and $x_{r1}$ and $x_{r2}$ are randomly selected from the individuals that dominate $x_i$. If only one solution dominates $x_i$, the other parent individual can choose from the non-dominated solution of $x_i$. We assign $c{x}_{il}$ when it is equal to zero if $c{x}_{il}\le 0$, then assign $c{x}_{il}$ to $x_{il}$ if $x_{il}=x_{jl}$.

In the PICEA-g-AKNN algorithm, the fitness value of the candidate solution is positively correlated with the number of goal vectors dominated by this candidate solution, but the fitness value of the target vector is inversely proportional to the number of candidate solutions dominating this target vector. The PICEA-g-AKNN algorithm is shown in Figure 4.

5. Computational Examinations

For validating the TMS-PTRD model, a case study is presented based on coronavirus disease 2019 (COVID-19) in Wuhan, Hubei, China. Moreover, we design a series of TMS-PTRD models to verify the PICEA-g-AKNN algorithm.

5.1. Simulation of Forecasting Phase

This part uses a modified SEIR epidemic dynamics model to evaluate COVID-19 in Wuhan. Wuhan, which covers an area of 8483 km² and has more than 10 million permanent residents, was one of the most seriously affected areas during the COVID-19 outbreak. Since 23 January 2020, the urban area of Wuhan has been effectively blocked. According to the Wuhan Novel Pneumonia Prevention and Control Headquarters, some temporary hospitals were established to receive infected people with mild symptoms, and some designated hospitals were selected to treat crucial patients. In Wuhan city, the population size within the lockdown area has not changed significantly. Because of the intervention measures, people were asked to stay at home to minimize contact with pedestrians.

The actual data were selected from the daily data of the COVID-19 epidemic in Wuhan from 21 January to 21 April 2020 (wjw.hubei.gov.cn/, accessed on 25 April 2020). The initial values were chosen as: N = 1,2320,000, E(0) = 149, I(0) = 440, R(0) = 28, H(0) = 102, SI(0) = 2197, EI(0) = 1394, S(0) = N−E(0)−I(0)−R(0)−H(0)−SI(0)−I(0). In line with Chao et al. [50], the parameters are shown in Table 1.

Note that because of changes in statistical methods, official figures of the cases have risen sharply since 12 February, and the model parameters have been modified appropriately. Figure 5 shows that the number of confirmed patients predicted by the SEIR model reaches a maximum at around 19 February. Considering the results show that the correlation coefficient between the predicted number and the actual number is more than 90%, there is no statistically significant difference between the actual number and the predicted number. So, this forecasting accuracy is enough for emergency management decision making.

5.2. Simulation of the Second Stage

Suppose that the lockdown city has 10 affected areas, and each affected area has a fever clinic where the general physicians diagnose illnesses and transfer confirmed cases to the appropriate hospitals. There are three TRDCs, three THs, and three DHs. We set the planning horizon to one week, taking into account the availability of external relief supplies and the additional emergency facilities in the affected areas for some time after the epidemic outbreak. Because it is difficult to obtain the start date of the epidemic, a certain time in the past is selected as the start date of the planning period. In a scenario, the parameters of the affected areas are the same as Table 1, and the initial values are shown in Table 2.

The parameters of the above scenarios are based on the implementation of strict control measures during the COVID-19 outbreak in Wuhan. So, hospital patients do not cause infection, and the exposure ratio of the patient is zero. Because the exposure ratio is likely to increase under conditions of inadequate intervention measures, we attempt to assess scenarios with inadequate measures. In addition, there is a relationship between the value of $q$ and the number of infected cases. To analyze the quarantine during the COVID-19 outbreak, different isolation ratio scenarios are considered in the TMS-PTRD model. We consider the case study consisting of four scenarios: S₁ ($r1=r$, $q1=q$), S₂ ($r1=2r$, $q1=q$), S₃ ($r1=r$, $q1=0.8q$), and S₄ ($r1=2r$, $q1=0.8q$). Each scenario is assumed to occur with equal probability ($p_s=\frac{1}{\left|S\right|}$). Within the next week, the number of infected people and critical patients in Disaster Area 1 is calculated by the modified SEIR model, as shown in Table 3.

Assume that the transportation cost ($c_{nm}$) and the delivery time ($t_{nm}$) from point n to point m are linear. The information is given in Table 4.

A standard carton is used as a transportation unit, including a set of medical and ancillary supplies. The loading capacity of the transportation vehicle is 50. Table 5 and Table 6 describe the parameters of the hospitals and the TRDCs, respectively.

In this section, experimental results are provided. The initial values of the affected areas were set as in Table 2, and the predicted values are shown in Table 3. In this case, the key parameters are shown in Table 4, Table 5 and Table 6. The model was coded using CPLEX and MATLAB software version R2019b, and the numerical experiments were run on a PC with 8 GB of RAM. In this study, we set the priority of critical patients to five. Considering the correlation coefficient between the predicted number and the actual number is more than 90%, the variation domain/perturbation level is set to 10% and the budget parameter of the confidence level is set to 95%.

Figure 6 and Figure 7 show the results of solving a small-scale multi-objective problem by the two methods of the ε-constraint algorithm and the PICEA-g-AKNN algorithm, respectively. In the test problem, the first objective function is on a scale of 10², and the second function is on a scale of 10³. Similarly, the third objective function is on a scale of 10⁵, and the fourth function is on a scale of 10⁻¹.

From Figure 6 and Figure 7, we calculated that the mean error value of the two algorithms in the first objective was 0.67%. Similarly, the mean error values of the second, third, and fourth objectives were 0.55%, 1.02%, and 0.78%, respectively. So, the PICEA-g-AKNN algorithm is credible in solving TMS-PTRD problems.

To analyze the uncertain parameters, the TMS-PTRD problems were tested under different uncertain conditions, and the output of the TMS-PTRD implementation is shown. In Figure 8, the main objective function increases when the uncertain environment gets worse. The results also show that an accurate assessment of the uncertainty can positively impact the results, where $a$ is the variation domain/perturbation level, and $\alpha$ is the budget parameter of the confidence level.

Table 7 shows the relief flows sent from the established TRDCs to the fever clinics for all periods ($a=0.1,\alpha =0.95$). For example, TRDCs 1 and 2 are established, and relief supplies are distributed to 10 fever clinics through these two established TRDCs.

Table 8 and Table 9 show the suitable locations for the established THs and DHs in each scenario. For example, all of the THs (three temporary hospitals) are established under Scenario 4, and we just need to establish two THs under Scenarios 1, 2, and 3. In addition, Table 8 and Table 9 determine the number of infected people and critical patients transferred from each fever clinic to hospitals for all periods ($a=0.1,\alpha =0.95$).

5.3. Computational Performance Analysis

We consider the 20, 30, and 50 affected areas, that is, the number of fever clinics that should be serviced. For each affected area, the parameters of the modified SEIR were randomly selected from an interval. Eighteen instances were randomly generated to test the performance of the PICEA-g-AKNN algorithm. In Table 10, for example, T-4-20-4-4-7-4 means that the case with three periods of two uncertain scenarios has 4 TRDCs, 20 fever clinics, 4 THs, and 4 DHs. Initial values of the disaster areas were randomly generated in the intervals shown in Table 11.

Without loss of generality, we set the variation perturbation level at 10% and the budget parameter of the confidence levels at 95%. Experimental results were obtained by a personal computer with 2.50 GHz and 8 GB of RAM. Each test instance was run 10 times to obtain the approximate solutions. The experimental results of the PICEA-G-AKNN algorithm and PICEA-g, MOEA/D, and NSGA-II were compared. The parameters of the algorithms are as Table 12.

Set coverage (C-metric) is used to compare the efficiency of the multi-objective algorithms [61]. In a minimization problem, the C-metric can be represented as follows.

(50)$C\left(A,B\right)=\frac{\left\{y\in B\left|\exists x\in A:x\prec y\right.\right\}}{\left|B\right|}$

where A is Pareto solutions obtained by one multi-objective algorithm, and B is Pareto solutions obtained by another multi-objective algorithm. The C-metric mean values of the numerical experiments are shown in Table 13.

From the results, the PICEA-g-AKNN algorithm outperforms the other three multi-objective algorithms in solving the TMS-PTRD problems. In particular, the PICEA-g-AKNN algorithm performs well in solving the TMS-PTRD problems with a large number of scenarios and periods. Considering the characteristics of the multi-scenario and multi-period problems, the novel neighborhood method proposed in Section 4.2.1 can reasonably construct a neighborhood of the TMS-PTRD problems. In addition, in the PICEA-g-AKNN algorithm, an appropriate method is used to evaluate the multi-objective solutions, and different evolutionary strategies are assigned to different types of chromosomes, which can effectively improve the efficiency of the multi-objective algorithm. We believe the proposed algorithm has advantages in solving multi-objective problems with multi-scenario and multi-period problems.

To analyze the epidemic interventions in the lockdown area of COVID-19, we considered different levels of epidemic interventions in the lockdown area. When strict intervention is used, the isolation ratio $q$ increases but the exposure rate $r$ decreases. Conversely, lax intervention leads to a decrease in the isolation ratio and an increase in the exposure rate. We tested T-8-50-8-8-14-12 by the PICEA-g-AKNN algorithm and showed the results as follows.

From Figure 9, Figure 10, Figure 11 and Figure 12, we can see that with the increase in the exposure rate $r$, the first three objectives increase significantly. However, the increase in the isolation ratio $q$ has no obvious effect on the results. We can conclude that reducing the average number of daily contacts and increasing the proportion of quarantined can reduce the expectation of the objectives. The average number of contacts has a more important impact on the development of the epidemic than the proportion quarantined. So, the first three objectives are better when strict intervention has been used. As can be seen from the fourth objective, the presented research can provide reliable equitable relief distribution under different epidemic interventions.

In this study, the TMS-PTRD model and the PICEA-g-AKNN algorithm have the advantages in tackling public health emergencies: (1) A reservation mechanism of relief supplies in densely populated urban areas is very important for epidemic intervention; (2) A relatively accurate method of epidemic prediction can support the decision of emergency management, for instance, setting up a reasonable number of temporary hospitals in different scenarios; (3) The TMS-PTRD model can be used to coordinate patient transfer and relief distribution problems, and multi-objective results can support decision makers to make decisions in their preferences and actual situations; (4) A multi-objective algorithm developed according to the TMS-PTRD model structure is helpful to obtain better Pareto solutions; (5) Reducing person-to-person contact has a positive significance for emergency management in the epidemic.

6. Conclusions

In response to public health emergencies, scientific and reasonable emergency relief can reduce epidemic disaster losses. Because of the suddenness and unpredictability of an epidemic outbreak, the TMS-PTRD model has advantages in pre-disaster preparation and post-disaster rescue. Considering that humanitarian factors cannot be ignored in emergency relief, timeliness and fairness are even more important than economy in emergency rescue. In order to ensure that the study can be more realistic, uncertainty factors and epidemic control measures were considered. This paper investigated the patient transfer and relief distribution problem in lockdown areas of COVID-19, which was depicted by a two-stage multi-objective stochastic model. The proposed model considers the following objectives: the total number of untreated infected patients at different periods, the total transfer time, the overall cost, and the equity distribution of relief supplies. To better cope with epidemic disasters, this paper considered setting up emergency facilities and storing a certain number of emergency supplies before an epidemic outbreak occurs. After an outbreak, the transfer plans of patients and the relief distribution are determined. Considering stochastic parameters, this study used chance-constraint programming to define the condition of a reliable set. To efficiently solve the large-scale TMS-PTRD problem, we proposed the PICEA-g-AKNN algorithm based on a novel similarity distance and tailored evolutionary strategies.

A real-world case study of Hunan of China and 18 test instances were generated to evaluate the TMS-PTRD model and the proposed algorithm. Numerical experiments show that the proposed method can well meet the needs of infected patients in an epidemic outbreak. In other words, hospital beds and relief supplies can be provided to infected patients in a timely manner. At the same time, numerical experiments showed that exposure reduction is more effective than other control measures. In this paper, the improved multi-objective algorithm can effectively solve large-scale TMS-PTRD problems. Experimental results show that the proposed PICEA-g-AKNN algorithm outperforms the PICEA-g, MOEA/D, and NSGA-II algorithms in solving the TMS-PTRD model. Moreover, the proposed TMS-PTRD model and PICEA-g-AKNN algorithm can effectively deal with the transfer of epidemic patients and the relief distribution.

In the future, the following extension issues can be considered. Patient migration and psychological panic should be considered in emergency management. Second, it is necessary to consider the routing and scheduling of delivery vehicles for the patients and the relief supplies. Third, more epidemiological interventions in lockdown areas need to be considered in emergency management.

Author Contributions

Conceptualization, S.L. (Shengjie Long) and D.Z.; Methodology, S.L. (Shengjie Long); Software, S.L. (Shengjie Long) and S.L. (Shuanglin Li); Validation, D.Z. and S.L. (Shengjie Long); Formal analysis, S.L. (Shuanglin Li); Investigation, S.L. (Shuanglin Li); Resources, D.Z.; Data curation, D.Z.; Writing—original draft preparation, S.L. (Shengjie Long); Writing—review and editing, S.L. (Shengjie Long) and S.L. (Shuanglin Li); Visualization, S.L. (Shuangyan Li); Supervision, D.Z. and S.L. (Shuangyan Li); Project administration, S.L. (Shuangyan Li); Funding acquisition, D.Z. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.

Acknowledgments

The authors would like to sincerely thank the editors and the anonymous reviewers for their constructive comments that greatly contributed to improving the manuscript.

Conflicts of Interest

The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Patient transfer and relief distribution network in lockdown area of epidemic.

Figure 2. An encoded chromosome for emergency facility location.

Figure 3. An encoded chromosome for relief planning.

Figure 4. PICEA-g-AKNN framework.

Figure 5. Theoretical estimates of infected patients and the critical patients by the modified SEIR model.

Figure 6. Results for a small test problem by an ε-constraint algorithm.

Figure 7. Results for a small test problem by the PICEA-g-AKNN algorithm.

Figure 8. Results in different uncertainty parameters.

Figure 9. Boxplot chart of the total number of untreated infected patients.

Figure 10. Boxplot chart of the total transfer time.

Figure 11. Boxplot chart of the overall cost.

Figure 12. Boxplot chart of criteria of equitable distribution.

Table 1

Parameters in modified SEIR model.

Parameters	r	β	q	α	δI	δq	γI	γH
	2	0.045	1.0 × 10−6	2.7 × 10−4	0.13	0.13	0.007	0.014

Table 2

Initial values of affected areas.

Affected Areas	1	2	3	4	5	6	7	8	9	10
N	150,005	113,553	175,424	226,502	200,740	164,807	187,551	159,922	210,337	209,542
E(0)	125	131	320	225	477	205	320	103	98	155
I(0)	5	3	9	8	5	7	9	4	7	3
R(0)	37	41	33	36	40	23	35	20	15	32
H(0)	1	3	5	4	2	5	2	3	0	1
SI(0)	219	178	409	291	611	321	400	184	129	233
EI(0)	139	171	230	180	320	189	306	99	102	174

Table 3

Predicted information of Disaster Area 1.

dI(t)/dH(t)	Period
Scenario	1	2	3	4	5	6	7
1	15/4	28/8	38/11	46/13	51/15	55/16	57/17
2	15/5	31/10	45/14	56/18	64/21	70/23	77/25
3	15/4	30/9	43/13	54/16	61/18	67/19	70/20
4	25/7	47/14	65/18	78/22	88/25	94/27	98/29

Table 4

Transportation cost and time distance.

${\mathit{c}}_{\mathit{n}\mathit{m}}/{\mathit{t}}_{\mathit{n}\mathit{m}}$	Fever Clinic
		1	2	3	4	5	6	7	8	9	10
TRDC	1	5/1.2	5/1.1	6/2.0	7/1.6	3/1.1	6/1.5	4/2.2	7/3.0	9/2.5	4/1.3
	2	7/2.0	4/2.2	5/2.3	3/1.5	2/1.0	5/1.3	5/1.3	6/1.5	4/1.8	6/2.4
	3	4/1.5	5/1.7	3/1.0	5/1.5	6/2.0	3/1.1	7/2.6	4/1.6	6/1.7	7/1.9
TH	1	6/1.6	7/1.9	5/1.5	4/1.7	6/2.4	5/2.2	8/2.5	3/1.0	5/1.1	7/1.4
	2	7/2.3	6/2.2	6/3.0	3/1.2	4/1.2	6/1.8	5/2.0	7/2.0	5/1.5	6/1.4
	3	5/1.5	8/2.0	7/2.2	8/2.7	4/2.0	6/1.8	7/2.1	5/1.8	9/2.9	7/3.0
DH	1	5/1.2	6/2.0	3/1.4	6/1.8	4/1.1	7/1.5	5/1.1	8/2.4	3/1.1	8/2.0
	2	7/1.9	5/1.5	7/2.6	6/2.5	5/1.4	8/2.0	8/3.2	6/1.5	3/2.0	7/2.3
	3	5/2.4	6/2.1	4/1.5	6/1.9	7/2.0	3/1.1	6/1.6	9/2.0	5/1.8	4/1.1

Table 5

Hospital parameters.

		Capacity	Fixed Cost
TH	1	2000	30,000
	2	3000	40,000
	3	3000	44,000
DH	1	900	30,000
	2	1000	50,000
	3	800	40,000

Table 6

TRDC parameters.

		Fixed Cost	Variable Cost
TRDC	1	30,000	1.2
	2	25,000	1.5
	3	20,000	1.7

Table 7

Relief commodities sent from established TRDCs to fever clinics.

Scenario	TRDC	Fever Clinic
		1	2	3	4	5	6	7	8	9	10
1	1	336		29			15	837			460
	2	35	456	652	521	931	513		279	274
	3
2	1	263	51	46	18	69	41	713	57	26	393
	2	160	380	630	527	820	495	67	252	259	71
	3
3	1	347	36	77		32	72	620	38		409
	2	56	390	710	499	861	441	181	237	260	24
	3
4	1	403	17		21		114	666	94	19	393
	2	102	389	766	559	828	417	58	192	258	83
	3

Table 8

Infected people transferred from each fever clinic to established THs.

Scenario	TH	Fever Clinic
		1	2	3	4	5	6	7	8	9	10
1	1	311	377	559					240	234	382
	2				433	760	441	469
	3
2	1	382	398	594					288	263
	2				507	829	496	719			429
	3
3	1	362	390	756					248	244
	2				470	791	478	714			394
	3
4	1		418	781					295	274
	2				541	866	266	752			435
	3	523					283

Table 9

Critical patients transferred from each fever clinic to established DHs.

Scenario	DH	Fever Clinic
		1	2	3	4	5	6	7	8	9	10
1	1	98	80	167	44	229		207		75
	2
	3		38	4	91		138		76		119
2	1	134		212		251		221		82
	2
	3		125	19	155		155		88		134
3	1	115	12	182	54	240		218		79
	2
	3		112		90		147		79		124
4	1	161		165		262		227		85
	2		131						90
	3			73	195		168				135

Table 10

Parameters of the test instances.

Instance	|I|	|C|	|J|	|H|	|T|	|S|
T-4-20-4-4-7-4	4	20	4	4	7	4
T-4-20-4-4-7-8	4	20	4	4	7	8
T-4-20-4-4-7-12	4	20	4	4	7	12
T-4-20-4-4-14-4	4	20	4	4	14	4
T-5-20-4-4-14-8	4	20	4	4	14	8
T-5-20-4-4-14-12	4	20	4	4	14	12
T-5-30-5-5-7-4	5	30	5	5	7	4
T-5-30-5-5-7-8	5	30	5	5	7	8
T-5-30-5-5-7-12	5	30	5	5	7	12
T-5-30-5-5-14-4	5	30	5	5	14	4
T-5-30-5-5-14-8	5	30	5	5	14	8
T-5-30-5-5-14-12	5	30	5	5	14	12
T-8-50-8-8-7-4	8	50	8	8	7	4
T-8-50-8-8-7-8	8	50	8	8	7	8
T-8-50-8-8-7-12	8	50	8	8	7	12
T-8-50-8-8-14-4	8	50	8	8	14	4
T-8-50-8-8-14-8	8	50	8	8	14	8
T-8-50-8-8-14-12	8	50	8	8	14	12

Table 11

Initial values of the starting time of the disaster areas.

Initial Values	N	E(0)	I(0)	R(0)	H(0)	Sq(0)	Eq(0)
	[15,000~30,000]	[100~200]	[0~10]	[0~10]	[0~10]	[100~200]	[100~200]

Table 12

Algorithm parameters.

Parameters	PICEA-g-ANK	PICEA-g	MOEA/D	NSGA-II
Maximum generations maxGenPopulation size NNumber of goal vectors NgProbability pcProbability pm	50001001000.80.2	50001001000.80.2	5000100-0.80.2	5000100-0.80.2
Neighborhood size T	N/20-N/5	-	20	-

Table 13

Results of the numerical experiments (A: PICEA-g-AKNN, B: PICEA-g, C: MOEA/D, D: NSGA-II).

Instance	$\mathit{C}\left(\mathit{A},\mathit{B}\right)$	$\mathit{C}\left(\mathit{B},\mathit{A}\right)$	$\mathit{C}\left(\mathit{A},\mathit{C}\right)$	$\mathit{C}\left(\mathit{C},\mathit{A}\right)$	$\mathit{C}\left(\mathit{A},\mathit{D}\right)$	$\mathit{C}\left(\mathit{D},\mathit{A}\right)$
T-4-20-4-4-7-4	0.213	0.253	0.202	0.183	0.212	0.013
T-4-20-4-4-7-8	0.266	0.135	0.214	0.288	0.302	0.044
T-4-20-4-4-7-12	0.184	0.117	0.25	0	0.360	0.161
T-4-20-4-4-14-4	0.20	0.304	0.202	0.104	0.105	0.255
T-5-20-4-4-14-8	0.496	0.121	0.24	0.166	0.406	0
T-5-20-4-4-14-12	0.211	0.058	0.057	0.204	0.237	0.187
T-5-30-5-5-7-4	0.419	0.013	0.358	0.072	0.023	0.303
T-5-30-5-5-7-8	0.168	0.263	0.072	0.227	0.042	0.164
T-5-30-5-5-7-12	0.283	0.105	0.299	0.056	0.183	0.133
T-5-30-5-5-14-4	0.218	0.152	0.195	0.012	0.153	0.389
T-5-30-5-5-14-8	0.285	0.151	0.158	0.221	0.248	0.242
T-5-30-5-5-14-12	0.480	0	0.184	0.118	0.223	0.294
T-8-50-8-8-7-4	0.515	0.124	0.231	0.132	0.171	0.40
T-8-50-8-8-7-8	0.317	0.15	0.147	0.294	0.297	0.145
T-8-50-8-8-7-12	0.275	0.164	0	0.412	0.194	0.290
T-8-50-8-8-14-4	0.360	0.218	0.159	0.113	0.315	0.047
T-8-50-8-8-14-8	0.581	0	0.125	0.177	0.307	0.031
T-8-50-8-8-14-12	0.338	0.081	0.149	0.033	0.238	0.155
Average	0.3629	0.1506	0.2026	0.1758	0.2510	0.2031

References

1. Verity, R.; Okell, L.C.; Dorigatti, I.; Winskill, P.; Whittaker, C.; Imai, N.; Cuomo-Dannenburg, G.; Thompson, H.; Walker, P.G.; Fu, H. Estimates of the severity of coronavirus disease 2019: A model-based analysis. Lancet Infect. Dis.; 2020; 20, pp. 669-677. [DOI: https://dx.doi.org/10.1016/S1473-3099(20)30243-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32240634]

2. Scalera, N.M.; Mossad, S.B.J.P.m. The first pandemic of the 21st century: Review of the 2009 pandemic variant influenza A (H1N1) virus. Postgrad. Med.; 2009; 121, pp. 43-47. [DOI: https://dx.doi.org/10.3810/pgm.2009.09.2051] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19820273]

3. Chowell, G.; Nishiura, H. Transmission dynamics and control of Ebola virus disease (EVD): A review. BMC Med.; 2014; 12, 196. [DOI: https://dx.doi.org/10.1186/s12916-014-0196-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25300956]

4. Büyüktahtakın, İ.E.; des-Bordes, E.; Kıbış, E.Y. A new epidemics–logistics model: Insights into controlling the Ebola virus disease in West Africa. Eur. J. Oper. Res.; 2018; 265, pp. 1046-1063. [DOI: https://dx.doi.org/10.1016/j.ejor.2017.08.037]

5. Cozzolino, A. Humanitarian logistics and supply chain management. Humanitarian Logistics; Springer: Berlin/Heidelberg, Germany, 2012; pp. 5-16.

6. Duhamel, C.; Santos, A.C.; Brasil, D.; Châtelet, E.; Birregah, B. Connecting a population dynamic model with a multi-period location-allocation problem for post-disaster relief operations. Ann. Oper. Res.; 2016; 247, pp. 693-713. [DOI: https://dx.doi.org/10.1007/s10479-015-2104-1]

7. Moreno, A.; Alem, D.; Ferreira, D. Heuristic approaches for the multiperiod location-transportation problem with reuse of vehicles in emergency logistics. Comput. Oper. Res.; 2016; 69, pp. 79-96. [DOI: https://dx.doi.org/10.1016/j.cor.2015.12.002]

8. Vahdani, B.; Veysmoradi, D.; Shekari, N.; Mousavi, S.M. Multi-objective, multi-period location-routing model to distribute relief after earthquake by considering emergency roadway repair. Neural Comput. Appl.; 2018; 30, pp. 835-854. [DOI: https://dx.doi.org/10.1007/s00521-016-2696-7]

9. Ghasemi, P.; Khalili-Damghani, K.; Hafezalkotob, A.; Raissi, S. Uncertain multi-objective multi-commodity multi-period multi-vehicle location-allocation model for earthquake evacuation planning. Appl. Math. Comput.; 2019; 350, pp. 105-132. [DOI: https://dx.doi.org/10.1016/j.amc.2018.12.061]

10. Elçi, Ö.; Noyan, N. A chance-constrained two-stage stochastic programming model for humanitarian relief network design. Transp. Res. Part B Methodol.; 2018; 108, pp. 55-83. [DOI: https://dx.doi.org/10.1016/j.trb.2017.12.002]

11. Banomyong, R.; Varadejsatitwong, P.; Oloruntoba, R. A systematic review of humanitarian operations, humanitarian logistics and humanitarian supply chain performance literature 2005 to 2016. Ann. Oper. Res.; 2019; 283, pp. 71-86. [DOI: https://dx.doi.org/10.1007/s10479-017-2549-5]

12. Grass, E.; Fischer, K.; Science, M. Two-stage stochastic programming in disaster management: A literature survey. Surv. Oper. Res.; 2016; 21, pp. 85-100. [DOI: https://dx.doi.org/10.1016/j.sorms.2016.11.002]

13. Ahmadi, M.; Seifi, A.; Tootooni, B. A humanitarian logistics model for disaster relief operation considering network failure and standard relief time: A case study on San Francisco district. Transp. Res. Part E Logist. Transp. Rev.; 2015; 75, pp. 145-163. [DOI: https://dx.doi.org/10.1016/j.tre.2015.01.008]

14. Oksuz, M.K.; Satoglu, S.I. A two-stage stochastic model for location planning of temporary medical centers for disaster response. Int. J. Disaster Risk Reduct.; 2020; 44, 101426. [DOI: https://dx.doi.org/10.1016/j.ijdrr.2019.101426]

15. Vahdani, B.; Veysmoradi, D.; Noori, F.; Mansour, F. Two-stage multi-objective location-routing-inventory model for humanitarian logistics network design under uncertainty. Int. J. Disaster Risk Reduct.; 2018; 27, pp. 290-306. [DOI: https://dx.doi.org/10.1016/j.ijdrr.2017.10.015]

16. Zhang, J.; Liu, H.; Yu, G.; Ruan, J.; Chan, F.T. A three-stage and multi-objective stochastic programming model to improve the sustainable rescue ability by considering secondary disasters in emergency logistics. Comput. Ind. Eng.; 2019; 135, pp. 1145-1154. [DOI: https://dx.doi.org/10.1016/j.cie.2019.02.003]

17. Doodman, M.; Shokr, I.; Bozorgi-Amiri, A.; Jolai, F. Pre-positioning and dynamic operations planning in pre-and post-disaster phases with lateral transhipment under uncertainty and disruption. J. Ind. Eng. Int.; 2019; 15, pp. 53-68. [DOI: https://dx.doi.org/10.1007/s40092-019-0317-7]

18. Saadatseresht, M.; Mansourian, A.; Taleai, M. Evacuation planning using multiobjective evolutionary optimization approach. Eur. J. Oper. Res.; 2009; 198, pp. 305-314. [DOI: https://dx.doi.org/10.1016/j.ejor.2008.07.032]

19. Coutinho-Rodrigues, J.; Tralhão, L.; Alçada-Almeida, L. Solving a location-routing problem with a multiobjective approach: The design of urban evacuation plans. J. Transp. Geogr.; 2012; 22, pp. 206-218. [DOI: https://dx.doi.org/10.1016/j.jtrangeo.2012.01.006]

20. Goodarzian, F.; Taleizadeh, A.A.; Ghasemi, P.; Abraham, A. An integrated sustainable medical supply chain network during COVID-19. Eng. Appl. Artif. Intell.; 2021; 100, 104188. [DOI: https://dx.doi.org/10.1016/j.engappai.2021.104188] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33619424]

21. Mansoori, S.; Bozorgi-Amiri, A.; Pishvaee, M.S. A robust multi-objective humanitarian relief chain network design for earthquake response, with evacuation assumption under uncertainties. Neural Comput. Appl.; 2020; 32, pp. 2183-2203. [DOI: https://dx.doi.org/10.1007/s00521-019-04193-x]

22. Flores, I.; Ortuño, M.T.; Tirado, G.; Vitoriano, B. Supported evacuation for disaster relief through lexicographic goal programming. Mathematics; 2020; 8, 648. [DOI: https://dx.doi.org/10.3390/math8040648]

23. Goerigk, M.; Grün, B.J.O.S. A robust bus evacuation model with delayed scenario information. Or Spectr.; 2014; 36, pp. 923-948. [DOI: https://dx.doi.org/10.1007/s00291-014-0365-8]

24. Lim, G.J.; Zangeneh, S.; Baharnemati, M.R.; Assavapokee, T. A capacitated network flow optimization approach for short notice evacuation planning. Eur. J. Oper. Res.; 2012; 223, pp. 234-245. [DOI: https://dx.doi.org/10.1016/j.ejor.2012.06.004]

25. Shahparvari, S.; Abbasi, B. Robust stochastic vehicle routing and scheduling for bushfire emergency evacuation: An Australian case study. Transp. Res. Part A Policy Pract.; 2017; 104, pp. 32-49. [DOI: https://dx.doi.org/10.1016/j.tra.2017.04.036]

26. Shahparvari, S.; Abbasi, B.; Chhetri, P.; Abareshi, A. Fleet routing and scheduling in bushfire emergency evacuation: A regional case study of the Black Saturday bushfires in Australia. Transp. Res. Part D Transp. Environ.; 2019; 67, pp. 703-722. [DOI: https://dx.doi.org/10.1016/j.trd.2016.11.015]

27. Bayram, V.; Yaman, H. Shelter location and evacuation route assignment under uncertainty: A benders decomposition approach. Transp. Sci.; 2018; 52, pp. 416-436. [DOI: https://dx.doi.org/10.1287/trsc.2017.0762]

28. Liang, B.; Yang, D.; Qin, X.; Tinta, T. A risk-averse shelter location and evacuation routing assignment problem in an uncertain environment. Int. J. Environ. Res. Public Health; 2019; 16, 4007. [DOI: https://dx.doi.org/10.3390/ijerph16204007]

29. Sheu, J.-B.; Pan, C. A method for designing centralized emergency supply network to respond to large-scale natural disasters. Transp. Res. Part B Methodol.; 2014; 67, pp. 284-305. [DOI: https://dx.doi.org/10.1016/j.trb.2014.05.011]

30. Sabouhi, F.; Bozorgi-Amiri, A.; Moshref-Javadi, M.; Heydari, M. An integrated routing and scheduling model for evacuation and commodity distribution in large-scale disaster relief operations: A case study. Ann. Oper. Res.; 2019; 283, pp. 643-677. [DOI: https://dx.doi.org/10.1007/s10479-018-2807-1]

31. Setiawan, E.; Liu, J.; French, A. Resource location for relief distribution and victim evacuation after a sudden-onset disaster. IISE Trans.; 2019; 51, pp. 830-846. [DOI: https://dx.doi.org/10.1080/24725854.2018.1517284]

32. Dasaklis, T.K.; Pappis, C.P.; Rachaniotis, N.P. Epidemics control and logistics operations: A review. Int. J. Prod. Econ.; 2012; 139, pp. 393-410. [DOI: https://dx.doi.org/10.1016/j.ijpe.2012.05.023]

33. Bozorgi-Amiri, A.; Khorsi, M. A dynamic multi-objective location–routing model for relief logistic planning under uncertainty on demand, travel time, and cost parameters. Int. J. Adv. Manuf. Technol.; 2016; 85, pp. 1633-1648. [DOI: https://dx.doi.org/10.1007/s00170-015-7923-3]

34. Caunhye, A.M.; Nie, X.; Pokharel, S. Optimization models in emergency logistics: A literature review. Socio-Econ. Plan. Sci.; 2012; 46, pp. 4-13. [DOI: https://dx.doi.org/10.1016/j.seps.2011.04.004]

35. Fereiduni, M.; Shahanaghi, K. A robust optimization model for distribution and evacuation in the disaster response phase. J. Ind. Eng. Int.; 2017; 13, pp. 117-141. [DOI: https://dx.doi.org/10.1007/s40092-016-0173-7]

36. Vahdani, B.; Veysmoradi, D.; Mousavi, S.; Amiri, M. Planning for relief distribution, victim evacuation, redistricting and service sharing under uncertainty. Socio-Econ. Plan. Sci.; 2022; 80, 101158. [DOI: https://dx.doi.org/10.1016/j.seps.2021.101158]

37. Wang, J.; Shen, D.; Yu, M. Multiobjective optimization on hierarchical refugee evacuation and resource allocation for disaster management. Math. Probl. Eng.; 2020; 2020, 8395714. [DOI: https://dx.doi.org/10.1155/2020/8395714]

38. Ghasemi, P.; Khalili-Damghani, K.; Hafezalkotob, A.; Raissi, S. Stochastic optimization model for distribution and evacuation planning (A case study of Tehran earthquake). Socio-Econ. Plan. Sci.; 2020; 71, 100745. [DOI: https://dx.doi.org/10.1016/j.seps.2019.100745]

39. Mohammadi, S.; Darestani, S.A.; Vahdani, B.; Alinezhad, A. A robust neutrosophic fuzzy-based approach to integrate reliable facility location and routing decisions for disaster relief under fairness and aftershocks concerns. Comput. Ind. Eng.; 2020; 148, 106734. [DOI: https://dx.doi.org/10.1016/j.cie.2020.106734]

40. Deb, K.; Pratap, A.; Agarwal, S.; Meyarivan, T. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans. Evol. Comput.; 2002; 6, pp. 182-197. [DOI: https://dx.doi.org/10.1109/4235.996017]

41. Mostaghim, S.; Teich, J. Strategies for finding good local guides in multi-objective particle swarm optimization (MOPSO). Proceedings of the 2003 IEEE Swarm Intelligence Symposium. SIS’03 (Cat. No. 03EX706); Indianapolis, IN, USA, 26–26 April 2003; pp. 26-33.

42. Ishibuchi, H.; Tsukamoto, N.; Nojima, Y. Behavior of evolutionary many-objective optimization. Proceedings of the Tenth International Conference on Computer Modeling and Simulation (uksim 2008); Cambridge, UK, 1–3 April 2008; pp. 266-271.

43. Liu, H.-L.; Chen, L.; Zhang, Q.; Deb, K. An evolutionary many-objective optimisation algorithm with adaptive region decomposition. Proceedings of the 2016 IEEE Congress on Evolutionary Computation (CEC); Vancouver, BC, Canada, 24–29 July 2016; pp. 4763-4769.

44. Bader, J.; Zitzler, E. HypE: An algorithm for fast hypervolume-based many-objective optimization. Evol. Comput.; 2011; 19, pp. 45-76. [DOI: https://dx.doi.org/10.1162/EVCO_a_00009]

45. Purshouse, R.C.; Jalbă, C.; Fleming, P.J. Preference-driven co-evolutionary algorithms show promise for many-objective optimisation. Proceedings of the International Conference on Evolutionary Multi-Criterion Optimization; Ouro Preto, Brazil, 5–8 April 2011; pp. 136-150.

46. Wang, R.; Fleming, P.J.; Purshouse, R.C. General framework for localised multi-objective evolutionary algorithms. Inf. Sci.; 2014; 258, pp. 29-53. [DOI: https://dx.doi.org/10.1016/j.ins.2013.08.049]

47. Wang, R.; Purshouse, R.C.; Fleming, P.J. Preference-inspired coevolutionary algorithms for many-objective optimization. IEEE Trans. Evol. Comput.; 2012; 17, pp. 474-494. [DOI: https://dx.doi.org/10.1109/TEVC.2012.2204264]

48. Liu, M.; Zhang, D. A dynamic logistics model for medical resources allocation in an epidemic control with demand forecast updating. J. Oper. Res. Soc.; 2016; 67, pp. 841-852. [DOI: https://dx.doi.org/10.1057/jors.2015.105]

49. Mishra, B.K.; Saini, D.K. SEIRS epidemic model with delay for transmission of malicious objects in computer network. Appl. Math. Comput.; 2007; 188, pp. 1476-1482. [DOI: https://dx.doi.org/10.1016/j.amc.2006.11.012]

50. Cao, S.; Feng, P.; Shi, P. Study on the epidemic development of COVID-19 in Hubei province by a modified SEIR model. J. Zhejiang Univ.; 2020; 49, pp. 178-184.

51. Charnes, A.; Cooper, W.W. Chance-constrained programming. Manag. Sci.; 1959; 6, pp. 73-79. [DOI: https://dx.doi.org/10.1287/mnsc.6.1.73]

52. Ghomi-Avili, M.; Naeini, S.G.J.; Tavakkoli-Moghaddam, R.; Jabbarzadeh, A. A fuzzy pricing model for a green competitive closed-loop supply chain network design in the presence of disruptions. J. Clean. Prod.; 2018; 188, pp. 425-442. [DOI: https://dx.doi.org/10.1016/j.jclepro.2018.03.273]

53. Tang, M.; Ji, B.; Fang, X.; Yu, S.S. Discretization-Strategy-Based Solution for Berth Allocation and Quay Crane Assignment Problem. J. Mar. Sci. Eng.; 2022; 10, 495. [DOI: https://dx.doi.org/10.3390/jmse10040495]

54. Zhang, S.X.; Zheng, L.M.; Liu, L.; Zheng, S.Y.; Pan, Y.M. Decomposition-based multi-objective evolutionary algorithm with mating neighborhood sizes and reproduction operators adaptation. Soft Comput.; 2017; 21, pp. 6381-6392. [DOI: https://dx.doi.org/10.1007/s00500-016-2196-9]

55. Guo, W.; Liu, T.; Dai, F.; Xu, P. An improved whale optimization algorithm for feature selection. Comput. Mater. Contin.; 2020; 62, pp. 337-354. [DOI: https://dx.doi.org/10.32604/cmc.2020.06411]

56. Epitropakis, M.G.; Plagianakos, V.P.; Vrahatis, M.N. Finding multiple global optima exploiting differential evolution’s niching capability. Proceedings of the 2011 IEEE Symposium on Differential Evolution (SDE); Paris, France, 11–15 April 2011; pp. 1-8.

57. Zhang, Y.-H.; Gong, Y.-J.; Gao, Y.; Wang, H.; Zhang, J. Parameter-free voronoi neighborhood for evolutionary multimodal optimization. IEEE Trans. Evol. Comput.; 2019; 24, pp. 335-349. [DOI: https://dx.doi.org/10.1109/TEVC.2019.2921830]

58. Deb, K.; Beyer, H.-G. Self-adaptive genetic algorithms with simulated binary crossover. Evol. Comput.; 2001; 9, pp. 197-221. [DOI: https://dx.doi.org/10.1162/106365601750190406] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11382356]

59. Storn, R.; Price, K. Differential evolution—A simple and efficient heuristic for global optimization over continuous spaces. J. Glob. Optim.; 1997; 11, pp. 341-359. [DOI: https://dx.doi.org/10.1023/A:1008202821328]

60. Qin, A.K.; Huang, V.L.; Suganthan, P.N. Differential evolution algorithm with strategy adaptation for global numerical optimization. IEEE Trans. Evol. Comput.; 2008; 13, pp. 398-417. [DOI: https://dx.doi.org/10.1109/TEVC.2008.927706]

61. Lei, H.; Wang, R.; Laporte, G.J.C.; Research, O. Solving a multi-objective dynamic stochastic districting and routing problem with a co-evolutionary algorithm. Comput. Oper. Res.; 2016; 67, pp. 12-24. [DOI: https://dx.doi.org/10.1016/j.cor.2015.09.002]