# Multiobjective Optimization on Hierarchical Refugee Evacuation and Resource Allocation for Disaster Management

Authors: Jian Wang; Danqing Shen; Mingzhu Yu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

This paper studies a location-allocation problem to determine the selection of emergency shelters, medical centers, and distribution centers after the disaster. The evacuation of refugees and allocation of relief resources are also considered. A mixed-integer nonlinear multiobjective programming model is proposed to characterize the problem. The hierarchical demand of different refugees and the limitations of relief resources are considered in the model. We employ a combination of the simulated annealing (SA) algorithm and the particle swarm optimization (PSO) algorithm method to solve the complex model. To optimize the result of our proposed algorithm, we absorb the group search, crossover, and mutation operator of GA into SA. We conduct a case study in a district of Beijing in China to validate the proposed methodology. Some computational experiments are conducted to analyze the impact of different factors, such as the target weight setting, selection of candidate shelters, and quantity of relief resources.


---

## Body

## 1. Introduction

A natural disaster is defined as a natural occurrence or activity that may cause injuries and deaths, loss of property and normal life, public disorder, economic recession, or environmental damage. The growing of population has escalated both the severity and frequency of natural disasters [1]. Research shows that natural disasters impact people with increasing frequency [2]. Site selection of emergency shelters (ESs) and emergency resource allocation have been studied in many disaster risk reduction projects, as they can provide safe havens, medical care, and life-saving resources for refugees [3–5].In recent years, three different categories of evacuation/refuge ESs are widely used after a natural disaster, including immediate shelters, short-term shelters, and long-term shelters. Immediate shelter is used for less than 1 day temporary evacuations, short-term shelter for less than 1 week duration of stay, and long-term shelter for more than 1 week and less than 1 month of stay. This hierarchical design of ESs has been widely used by many research studies [5, 6]. However, the urgency for evacuation and resource demand of refugees in the same period are often different, and refugees should be divided into different priorities [7]. Thus, ESs should be divided into different categories in correspondence with different priorities of refugees in the same period. Furthermore, relief resources are typically not sufficient to meet all demands in the early days of disaster [8]. Hence, the goals of the emergency relief resource allocation planning policy should not only focus on efficiency (e.g., minimizing the evacuation distance and total cost), but also consider the fairness with the limitations of relief resources.In this paper, we consider the facility site selection and relief resource allocation in the response phase of a disaster, especially during the critical 72-h time window after the disaster happens. We propose a hierarchical multiobjective model to determine site selection and the corresponding relief resource allocation. We consider not only traditional objectives, e.g., minimizing the hierarchical evacuation distance, minimizing the total costs (including construction costs of ESs and emergency distribution centers (EDCs), and fixed and transportation costs of emergency resources), and minimizing the total hierarchical unmet resources, but also maximizing the fairness of refugees. A combination of the SA algorithm and the PSO algorithm (SA-PSO) is applied to solve the complex model. The group search, crossover, and mutation operators of GA are combined into SA to optimize and expand the searching space of our proposed algorithm. Through the real-case computational studies, we find that (1) compared to SA-PSO, the improved SA algorithm and the PSO algorithm (ISA-PSO) can optimize the results by expanding its searching space, but population search also increases time consumption; (2) when the objective weight increases, the objective value of the total evacuation distance improves the most; (3) increasing the resources in EDCs can greatly optimize the target of total hierarchical unmet resources with a small increase in the cost; (4) when the evacuation demand of refugees is increased, it is necessary to select more candidate shelters; (5) compared to the results of time-varying demand estimation, emergency evacuation planning under stable demand would cause more evacuation distance and more unmet resources.The main contributions of our research are as follows. First, we formulate a hierarchical multiobjective location-allocation model to determine the selection of hierarchical ESs, emergency medical centers (EMCs), and EDCs, along with the corresponding hierarchical evacuation and relief resource allocation. Second, we proposed a heuristic algorithm to solve the complex nonlinear model. Furthermore, we conduct a case study to validate our model and analyze the impacts of different algorithms and factors. These analyses provide useful suggestions for the government on shelters and resource planning management.This paper is organized as the following. Section2 provides review of the related literature, while Section 3 describes the problem and presents a multiobjective nonlinear programming model. Then, Section 4 proposes two heuristic algorithms to solve this complex problem. We present computational results of the Chaoyang district, Beijing, China, as our case study in Section 5 and conclude in Section 6.


## 2. Literature Review

Emergency disaster management including site selection for emergency facilities and the corresponding resource allocation has been extensively studied in the literature.The site selection of ESs and EMCs has long been a critical and difficult issue. Nowadays, many models have been developed, including the minimal distance model like the P-center problem and the P-median problem and the covering model like the location set covering problem (LSCP) and the maximal coverage location problem (MCLP). Solution of the P-center problem and the P-median problem seeks to set the total number of selected ESs in advance. The objective of the P-center problem is to minimize the longest distance, while the objective of the P-median problem is to minimize the total travel distance from the AAs to the ESs [9, 10]. The LSCP and the MCLP searches for the optimal location of ESs to cover all demands within a certain service distance. The LSCP is to minimize the total number of selected facilities as well as to minimize the construction cost, and the MCLP is to minimize the number of refugees that will not be served within the service distance [11]. However, in an actual response period, the project manager must deal with multiple conflicting objectives in the varying environment [12]. These single-objective models are too simplistic to solve the ES and EMC site selection problems because of the ignorance of important objectives [13]. To solve the complex problem, multiobjective models based on the minimal distance model and the covering model have been developed by many studies. Tsai and Yeh [14] integrated the MCLP with the P-center problem and developed a double-target (minimizing the farthest distance and maximizing the total number of refugees in these candidate ESs) model to decide the site selection problem. Hu et al. [15] considered multiple criteria like the safety of victims, minimization of total evacuation distance, cost reduction, and capacity constraint satisfaction to solve the ES selection problem. Herein, site selection of ESs should consider multiple criteria such as minimizing the total distance/time, minimizing the total cost, and maximizing the coverage.In addition, varying environment corresponds to different levels of service need [5]. Chen et al. [6] developed a three-level hierarchical ES location model based on the varying need of different periods after a disaster, with the objective of minimizing the total hierarchical distance. Previous classification of ESs is based only on the time-varying demand of different periods; however, no study has researched on the varying need of different levels of refugees to different levels of ESs in the same period. And, among these classified earthquake evacuation ESs, according to Xu et al. [16], the immediate shelter is not necessary to determine a specific model for its site selection because it is only used for an immediate evacuation. And, because the long-term shelter is used for a longer period, there are limited candidate locations of these shelters and are normally away from the city center. Thus, the location of a long-term shelter is not a complex problem. But, the location of a short-term shelter must meet multiple criteria and constraints as previously introduced. Since the service demand of refugees varies in the same period, the most critical site selection of the short-term shelter should consider hierarchical division based on different demand of refugees.In the early stage of emergency disaster response, site selection of ESs and allocation of refugees are of great importance to emergency management. However, effective distribution of relief resources is also an indispensable part of this life-saving project [7]. In recent years, some studies focused on simultaneous site selection as well as relief resource allocation during the response phase. A three-stage programming model is designed by Sheu and Pan [17] for disaster management, attempting to identify the integration of the ES location network and the EMC location network. The objective functions of each stage involve not only the most fundamental criteria such as minimizing the total distance and cost, but also consider the psychological cost perceived by refugees in AAs, ESs, and EMCs. Rodríguez-Espíndola and Gaytán [18] developed a multicriteria approach to solve the simultaneous site selection of ESs and EDCs, along with the allocation of human resources and goods by minimizing the total costs and total distance. Dalal and Üster [19] proposed an overall disaster response model including supply (relief) and demand (evacuation) with the trade-offs between the biobjective of evacuation time and system cost, considering uncertainty in refugee demand. They presented a set of disaster scenarios to optimize the average and worst case costs. Furthermore, the demand for emergency resources is uncertain in emergency resource allocation problems [4]. Bozorgi-Amiri and Khorsi [8] pointed out that the relief resources are often not enough to meet the demand and presented a relief logistic planning model making the trade-off between three objectives. Hence, considering the insufficient feature of relief resources and the hierarchical demand of different refugees, the objective of this project should contain not only multiple criteria like cost, distance, and resource satisfaction, but also maximize the fairness of refugees. Fairness can be measured by the hierarchical demand for evacuation and relief resources.Our model differs from most existing studies [4, 7, 16, 17] in two aspects. The first is the formulation of a decision-making project that contemplates a novel model on the assignment of different priority refugees to different levels of ESs, using hierarchy to evaluate the urgent demand of refugees and improve efficiency of management. The other aspect is the allocating of resources considering insufficient relief resources and hierarchical requirement of each demand point.In this paper, we consider the 72-hours gold rescue time after the disaster. Based on these previous research studies, we propose a hierarchical multiobjective nonlinear model, which considers different demands of different levels of refugees, simultaneous site selection of two-level ESs, EMCs, and EDCs, and the allocation of different level refugees and insufficient relief resources. Our study is distinguished from others by bringing a hierarchical design of refugees and ESs to optimize evacuation efficiency, using different hierarchy of refugees and ESs to measure their urgency for evacuation and resource demand. Furthermore, we use hierarchical management to allocate insufficient resources to different priority refugees in AAs, two-level ESs, and EMCs to achieve fairness. In this research, we use GIS to generate distance data to test the effectiveness and efficiency of this model in a case study. Finally, we use ArcGIS and visualize the comparison of target values to examine the trade-offs among three targets and parametric settings.


## 3. Problem Definition and Model Formulation

In this section, we provide the problem definition and model formulation on the emergency evacuation and the resource allocation network. The model considers the specific evacuation and the distribution network problem with hierarchical ESs, EMCs, and AAs (Figure1) in an advanced or a developed city because such a city has basic buildings like public parks that can provide short-term emergency evacuation and adequate hospitals for injuries.Figure 1
Illustration of an evacuation and distribution network.In the evacuation and distribution network, we define the AAs as the source nodes whose populations dictate the required outflow. Another set of nodes of EMCs, A-ESs (A-priority emergency shelters) and B-ESs (B-priority emergency shelters), are the destination nodes. EMCs provide services for injuries, and A-ESs provide higher priority services than B-ESs for refugees. A final set of nodes of EDCs represent the source nodes of resources which provide resources to demand points.


### 3.1. Assumptions

The proposed model is based on the following assumptions:Assumption 1. The time-varying demand of refugees and resources is given. Based on the basic population data of AAs, the number of different priority refugees and injuries is roughly estimated. Besides, the quantity of required resources can also be estimated from the number of people in each demand point.Assumption 2. Refugees are classified as nonevacuated, A-priority evacuated people, B-priority evacuated people, and injured people, and their corresponding places are AAs, A-ESs, B-ESs, and EMCs. We set the order of priority from high to low as injured, A-priority evacuated, B-priority evacuated, and nonevacuated.Assumption 3. The resources are classified into two categories: commodities and medical supplies. Injuries and different hierarchy of refugees require different kinds of resources. The refugees with higher priority can receive more and closer resources.Assumption 4. The EMCs are selected from the existing hospitals. Therefore, regardless of the construction cost of the EMCs, all EMCs can be the selected points. The cost involves the construction cost of A-ESs, B-ESs, and EDCs and the fixed and transportation cost of resources.Assumption 5. The number of evacuation and transportation vehicles is unlimited. Therefore, all refugees, injured, and resources can be transported to the corresponding nodes.


### 3.2. Model Formulation

#### 3.2.1. Sets and Parameters

Sets and parameters of this model, along with their brief description, are given as follows:Sets and indices:AA,a: set and index of the affected areas, a∈AAA-ES,j: set and index of A-priority emergency shelter locations, j∈A−ESB-ES,k: set and index of B-priority emergency shelter locations, k∈B−ESEMC,h: set and index of candidate existing emergency medical centers, h∈EMCEDC,c: set and index of emergency distribution centers, c∈EDCT, t: set and index of the time period, t∈TR, r: set and index of the type of relief resources, as commonly indexed by r; r=1 refers to commodities; r=2 refers to medical supplies, r∈RParameters:Ea: the total number of people in the affected area a, a∈AAEat: the number of refugees in the affected area a at time t, a∈AA, t∈TEAat: the number of A-priority refugees in the affected area a at time t, a∈AA, t∈TEBat: the number of B-priority refugees in the affected area a at time t, a∈AA, t∈TEIat: the number of injured people in the affected area a at time t, a∈AA, t∈TRrc: the quantity of relief resources r in the emergency distribution center c, r∈RMAψ: the maximum capacity (unit: person) of node ψ, ψ∈A−ES ∪ B-ES ∪ EMCCφ: the construction cost (unit: RMB) of node φ, φ∈A−ES ∪ B-ES ∪ EDCCRr: the fixed cost (unit: RMB) of per unit resources r, r∈RCTr: the transportation cost (unit: RMB/km) of per unit resources r, r∈RDaψ: the distance (unit: km) from the affected area a to node ψ, ψ∈A−ES ∪ B-ES ∪ EMCDcΦ: the distance (unit: km) from the emergency distribution center c to node Φ, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMCωΦ: the priority weight of node Φ, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMCδΦr: the number of resources r consumed by per person per unit time in node Φ, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMC, r∈Rα: the ratio of total resources in the emergency distribution center to the total required resources (considering the insufficient feature)


#### 3.2.2. Variables

Dependent variables and decision variables in this model are given as follows:Dependent variables:EAjt: the number of A-priority refugees in A-priority emergency shelter j at time t, j∈A−ES, t∈TEBkt: the number of B-priority refugees in B-priority emergency shelter k at time t, k∈B−ES, t∈TEHht: the number of injured people in emergency medical center h at time t, h∈EMC, t∈TRrΦt: the quantity of resources r required by node Φ at time t, r∈R, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMC, t∈TRTr: the total quantity of resources r required by all resource demand points throughout the rescue period of 72h, r∈RRurΦt: the quantity of unmet resources r in node Φ at time t, r∈R, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMC, t∈TFRurΦt: the loss for the unmet resources of node Φ at time t, μ<0,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,t∈TDecision variables:xψ: 1 if node ψ is selected, 0 otherwise, ψ∈A−ES ∪ B-ES ∪ EDC.EAajt: the number of A-priority refugees transported from affected area a to A-priority emergency shelter j at time t, a∈AA, j∈A−ES, t∈TEBakt: the number of B-priority refugees transported from affected area a to A-priority emergency shelter k at time t, a∈AA, k∈B−ES, t∈TEHaht: the number of injured people transported from affected area a to emergency medical center h at time t, a∈AA, h∈EMC, t∈TRrcΦt: the quantity of resources r transported from emergency distribution center c to node Φ at time t, r∈R, c∈EDC, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMC, t∈T


#### 3.2.3. Formulation

The proposed hierarchical multiobjective nonlinear mathematical formulation is presented as follows:Objective function:(1)minff1,f2,f3,(2)f1=min∑t∑a∑jDaj∗EAajt+∑kDak∗EBakt+∑hDah∗EHaht,(3)f2=min∑jxj∗Cj+∑kxk∗Ck+∑cxc∗Cc+∑t∑r∑c∑hCRr∗ωh∗Dch∗Rrcht+∑jCRr∗ωj∗Dcj∗Rrcjt+∑kCRr∗ωk∗Dck∗Rrckt+∑aCRr∗ωa∗Dca∗Rrcat,(4)f3=min∑t∑r∑hωh∗FRurht+∑jωj∗FRurjt+∑kωk∗FRurkt+∑aωa∗FRurat,is subject to(5)EAat=∑jEAajt,a∈AA,j∈A−ES,t∈T,(6)EBat=∑kEBakt,a∈AA,k∈B−ES,t∈T,(7)EIat=∑hEHaht,a∈AA,h∈EMC,t∈T,(8)∑tEAajt≤MAj∗xj,a∈AA,j∈A−ES,t∈T,(9)∑tEBakt≤MAk∗xk,a∈AA,k∈B−ES,t∈T,(10)∑tEHaht≤MAh,a∈AA,h∈EMC,t∈T,(11)∑t∑ΦRrcΦt≤Rrc∗xC,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,c∈EDC,t∈T,(12)EAjt=∑aEAajt,a∈AA,j∈A−ES,t=1,(13)EAjt=EAjt−1+∑aEAajt,a∈AA,j∈A−ES,t≥2,(14)EBkt=∑aEBakt,a∈AA,k∈B−ES,t=1,(15)EBkt=EBkt−1+∑aEBakt,a∈AA,k∈B−ES,t≥2,(16)EHht=∑aEHaht,a∈AA,h∈EMC,t=1,(17)EHht=EHht−1+∑aEHaht,a∈AA,h∈EMC,t≥2,(18)Eat=Ea−EAat−EBat−EIat,a∈AA,t=1,(19)Eat=Eat−1−EAat−EBat−EIat,a∈AA,t≥2,(20)Rrjt=δjr∗EAjt,r∈R,j∈A−ES,t∈T,(21)Rrkt=δkr∗EBkt,r∈R,k∈B−ES,t∈T,(22)Rrht=δhr∗EHht,r∈R,h∈EMC,t∈T,(23)Rrat=δar∗Eat,r∈R,a∈AA,t∈T,(24)RTr=∑t∑ΦRrΦt,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,t∈T,(25)RrΦt≥∑cRrcΦt,RrcΦt≥0,Φ∈AA∪A−ES∪B−ES∪EMC,c∈EDC,r∈R,t∈T,(26)∑cRrc∗xC=α∗RTr,c∈EDC,r∈R,(27)RurΦt=RrΦt−∑cRrcΦt,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,t∈T,c∈EDC,(28)FRurΦt=1−e−μ∗RurΦtμ,μ<0,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,t∈T,c∈EDC.The programming model involves three objective functions (minimization of the total evacuation distance, total cost, and total unmet resources), specified by equations (2)–(4). Constraints (5)–(7) make sure that all refugees requiring hierarchical evacuation in AAs can be evacuated into the corresponding nodes. The maximum capacities of A-ESs, B-ESs, EMCs, and EDCs are limited in constraints (8)–(11). Constraints (12)–(19) reflect the number of people in each node at time t. Constraints (20)–(23) define the quantity of required resources in each node at time t; thus, constraint (24) defines the total quantity of relief resources r required from all resource demand points throughout the rescue period of 72 h. To avoid the wastage of insufficient resources, we define that the quantity of received resources of each node is no more than the quantity of required resources of each node by equation (25). Considering the hierarchical resource demand of each node and insufficient resources of EDCs, to reduce the secondary damage caused by insufficient resources, constraint (26) sets the total quantity of resources in all EDCs (α≥1 means the total quantity of resources is sufficient). Constraint (27) indicates the quantity of unmet resources in each node. Constraint (28) defines the loss for the unmet resources of each node at time t, where the parameter μ is used to show the importance placed on large losses.


## 3.1. Assumptions

The proposed model is based on the following assumptions:Assumption 1. The time-varying demand of refugees and resources is given. Based on the basic population data of AAs, the number of different priority refugees and injuries is roughly estimated. Besides, the quantity of required resources can also be estimated from the number of people in each demand point.Assumption 2. Refugees are classified as nonevacuated, A-priority evacuated people, B-priority evacuated people, and injured people, and their corresponding places are AAs, A-ESs, B-ESs, and EMCs. We set the order of priority from high to low as injured, A-priority evacuated, B-priority evacuated, and nonevacuated.Assumption 3. The resources are classified into two categories: commodities and medical supplies. Injuries and different hierarchy of refugees require different kinds of resources. The refugees with higher priority can receive more and closer resources.Assumption 4. The EMCs are selected from the existing hospitals. Therefore, regardless of the construction cost of the EMCs, all EMCs can be the selected points. The cost involves the construction cost of A-ESs, B-ESs, and EDCs and the fixed and transportation cost of resources.Assumption 5. The number of evacuation and transportation vehicles is unlimited. Therefore, all refugees, injured, and resources can be transported to the corresponding nodes.


## 3.2. Model Formulation

### 3.2.1. Sets and Parameters

Sets and parameters of this model, along with their brief description, are given as follows:Sets and indices:AA,a: set and index of the affected areas, a∈AAA-ES,j: set and index of A-priority emergency shelter locations, j∈A−ESB-ES,k: set and index of B-priority emergency shelter locations, k∈B−ESEMC,h: set and index of candidate existing emergency medical centers, h∈EMCEDC,c: set and index of emergency distribution centers, c∈EDCT, t: set and index of the time period, t∈TR, r: set and index of the type of relief resources, as commonly indexed by r; r=1 refers to commodities; r=2 refers to medical supplies, r∈RParameters:Ea: the total number of people in the affected area a, a∈AAEat: the number of refugees in the affected area a at time t, a∈AA, t∈TEAat: the number of A-priority refugees in the affected area a at time t, a∈AA, t∈TEBat: the number of B-priority refugees in the affected area a at time t, a∈AA, t∈TEIat: the number of injured people in the affected area a at time t, a∈AA, t∈TRrc: the quantity of relief resources r in the emergency distribution center c, r∈RMAψ: the maximum capacity (unit: person) of node ψ, ψ∈A−ES ∪ B-ES ∪ EMCCφ: the construction cost (unit: RMB) of node φ, φ∈A−ES ∪ B-ES ∪ EDCCRr: the fixed cost (unit: RMB) of per unit resources r, r∈RCTr: the transportation cost (unit: RMB/km) of per unit resources r, r∈RDaψ: the distance (unit: km) from the affected area a to node ψ, ψ∈A−ES ∪ B-ES ∪ EMCDcΦ: the distance (unit: km) from the emergency distribution center c to node Φ, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMCωΦ: the priority weight of node Φ, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMCδΦr: the number of resources r consumed by per person per unit time in node Φ, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMC, r∈Rα: the ratio of total resources in the emergency distribution center to the total required resources (considering the insufficient feature)


### 3.2.2. Variables

Dependent variables and decision variables in this model are given as follows:Dependent variables:EAjt: the number of A-priority refugees in A-priority emergency shelter j at time t, j∈A−ES, t∈TEBkt: the number of B-priority refugees in B-priority emergency shelter k at time t, k∈B−ES, t∈TEHht: the number of injured people in emergency medical center h at time t, h∈EMC, t∈TRrΦt: the quantity of resources r required by node Φ at time t, r∈R, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMC, t∈TRTr: the total quantity of resources r required by all resource demand points throughout the rescue period of 72h, r∈RRurΦt: the quantity of unmet resources r in node Φ at time t, r∈R, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMC, t∈TFRurΦt: the loss for the unmet resources of node Φ at time t, μ<0,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,t∈TDecision variables:xψ: 1 if node ψ is selected, 0 otherwise, ψ∈A−ES ∪ B-ES ∪ EDC.EAajt: the number of A-priority refugees transported from affected area a to A-priority emergency shelter j at time t, a∈AA, j∈A−ES, t∈TEBakt: the number of B-priority refugees transported from affected area a to A-priority emergency shelter k at time t, a∈AA, k∈B−ES, t∈TEHaht: the number of injured people transported from affected area a to emergency medical center h at time t, a∈AA, h∈EMC, t∈TRrcΦt: the quantity of resources r transported from emergency distribution center c to node Φ at time t, r∈R, c∈EDC, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMC, t∈T


### 3.2.3. Formulation

The proposed hierarchical multiobjective nonlinear mathematical formulation is presented as follows:Objective function:(1)minff1,f2,f3,(2)f1=min∑t∑a∑jDaj∗EAajt+∑kDak∗EBakt+∑hDah∗EHaht,(3)f2=min∑jxj∗Cj+∑kxk∗Ck+∑cxc∗Cc+∑t∑r∑c∑hCRr∗ωh∗Dch∗Rrcht+∑jCRr∗ωj∗Dcj∗Rrcjt+∑kCRr∗ωk∗Dck∗Rrckt+∑aCRr∗ωa∗Dca∗Rrcat,(4)f3=min∑t∑r∑hωh∗FRurht+∑jωj∗FRurjt+∑kωk∗FRurkt+∑aωa∗FRurat,is subject to(5)EAat=∑jEAajt,a∈AA,j∈A−ES,t∈T,(6)EBat=∑kEBakt,a∈AA,k∈B−ES,t∈T,(7)EIat=∑hEHaht,a∈AA,h∈EMC,t∈T,(8)∑tEAajt≤MAj∗xj,a∈AA,j∈A−ES,t∈T,(9)∑tEBakt≤MAk∗xk,a∈AA,k∈B−ES,t∈T,(10)∑tEHaht≤MAh,a∈AA,h∈EMC,t∈T,(11)∑t∑ΦRrcΦt≤Rrc∗xC,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,c∈EDC,t∈T,(12)EAjt=∑aEAajt,a∈AA,j∈A−ES,t=1,(13)EAjt=EAjt−1+∑aEAajt,a∈AA,j∈A−ES,t≥2,(14)EBkt=∑aEBakt,a∈AA,k∈B−ES,t=1,(15)EBkt=EBkt−1+∑aEBakt,a∈AA,k∈B−ES,t≥2,(16)EHht=∑aEHaht,a∈AA,h∈EMC,t=1,(17)EHht=EHht−1+∑aEHaht,a∈AA,h∈EMC,t≥2,(18)Eat=Ea−EAat−EBat−EIat,a∈AA,t=1,(19)Eat=Eat−1−EAat−EBat−EIat,a∈AA,t≥2,(20)Rrjt=δjr∗EAjt,r∈R,j∈A−ES,t∈T,(21)Rrkt=δkr∗EBkt,r∈R,k∈B−ES,t∈T,(22)Rrht=δhr∗EHht,r∈R,h∈EMC,t∈T,(23)Rrat=δar∗Eat,r∈R,a∈AA,t∈T,(24)RTr=∑t∑ΦRrΦt,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,t∈T,(25)RrΦt≥∑cRrcΦt,RrcΦt≥0,Φ∈AA∪A−ES∪B−ES∪EMC,c∈EDC,r∈R,t∈T,(26)∑cRrc∗xC=α∗RTr,c∈EDC,r∈R,(27)RurΦt=RrΦt−∑cRrcΦt,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,t∈T,c∈EDC,(28)FRurΦt=1−e−μ∗RurΦtμ,μ<0,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,t∈T,c∈EDC.The programming model involves three objective functions (minimization of the total evacuation distance, total cost, and total unmet resources), specified by equations (2)–(4). Constraints (5)–(7) make sure that all refugees requiring hierarchical evacuation in AAs can be evacuated into the corresponding nodes. The maximum capacities of A-ESs, B-ESs, EMCs, and EDCs are limited in constraints (8)–(11). Constraints (12)–(19) reflect the number of people in each node at time t. Constraints (20)–(23) define the quantity of required resources in each node at time t; thus, constraint (24) defines the total quantity of relief resources r required from all resource demand points throughout the rescue period of 72 h. To avoid the wastage of insufficient resources, we define that the quantity of received resources of each node is no more than the quantity of required resources of each node by equation (25). Considering the hierarchical resource demand of each node and insufficient resources of EDCs, to reduce the secondary damage caused by insufficient resources, constraint (26) sets the total quantity of resources in all EDCs (α≥1 means the total quantity of resources is sufficient). Constraint (27) indicates the quantity of unmet resources in each node. Constraint (28) defines the loss for the unmet resources of each node at time t, where the parameter μ is used to show the importance placed on large losses.


## 3.2.1. Sets and Parameters

Sets and parameters of this model, along with their brief description, are given as follows:Sets and indices:AA,a: set and index of the affected areas, a∈AAA-ES,j: set and index of A-priority emergency shelter locations, j∈A−ESB-ES,k: set and index of B-priority emergency shelter locations, k∈B−ESEMC,h: set and index of candidate existing emergency medical centers, h∈EMCEDC,c: set and index of emergency distribution centers, c∈EDCT, t: set and index of the time period, t∈TR, r: set and index of the type of relief resources, as commonly indexed by r; r=1 refers to commodities; r=2 refers to medical supplies, r∈RParameters:Ea: the total number of people in the affected area a, a∈AAEat: the number of refugees in the affected area a at time t, a∈AA, t∈TEAat: the number of A-priority refugees in the affected area a at time t, a∈AA, t∈TEBat: the number of B-priority refugees in the affected area a at time t, a∈AA, t∈TEIat: the number of injured people in the affected area a at time t, a∈AA, t∈TRrc: the quantity of relief resources r in the emergency distribution center c, r∈RMAψ: the maximum capacity (unit: person) of node ψ, ψ∈A−ES ∪ B-ES ∪ EMCCφ: the construction cost (unit: RMB) of node φ, φ∈A−ES ∪ B-ES ∪ EDCCRr: the fixed cost (unit: RMB) of per unit resources r, r∈RCTr: the transportation cost (unit: RMB/km) of per unit resources r, r∈RDaψ: the distance (unit: km) from the affected area a to node ψ, ψ∈A−ES ∪ B-ES ∪ EMCDcΦ: the distance (unit: km) from the emergency distribution center c to node Φ, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMCωΦ: the priority weight of node Φ, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMCδΦr: the number of resources r consumed by per person per unit time in node Φ, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMC, r∈Rα: the ratio of total resources in the emergency distribution center to the total required resources (considering the insufficient feature)


## 3.2.2. Variables

Dependent variables and decision variables in this model are given as follows:Dependent variables:EAjt: the number of A-priority refugees in A-priority emergency shelter j at time t, j∈A−ES, t∈TEBkt: the number of B-priority refugees in B-priority emergency shelter k at time t, k∈B−ES, t∈TEHht: the number of injured people in emergency medical center h at time t, h∈EMC, t∈TRrΦt: the quantity of resources r required by node Φ at time t, r∈R, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMC, t∈TRTr: the total quantity of resources r required by all resource demand points throughout the rescue period of 72h, r∈RRurΦt: the quantity of unmet resources r in node Φ at time t, r∈R, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMC, t∈TFRurΦt: the loss for the unmet resources of node Φ at time t, μ<0,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,t∈TDecision variables:xψ: 1 if node ψ is selected, 0 otherwise, ψ∈A−ES ∪ B-ES ∪ EDC.EAajt: the number of A-priority refugees transported from affected area a to A-priority emergency shelter j at time t, a∈AA, j∈A−ES, t∈TEBakt: the number of B-priority refugees transported from affected area a to A-priority emergency shelter k at time t, a∈AA, k∈B−ES, t∈TEHaht: the number of injured people transported from affected area a to emergency medical center h at time t, a∈AA, h∈EMC, t∈TRrcΦt: the quantity of resources r transported from emergency distribution center c to node Φ at time t, r∈R, c∈EDC, Φ∈AA ∪ A-ES ∪ B-ES ∪ EMC, t∈T


## 3.2.3. Formulation

The proposed hierarchical multiobjective nonlinear mathematical formulation is presented as follows:Objective function:(1)minff1,f2,f3,(2)f1=min∑t∑a∑jDaj∗EAajt+∑kDak∗EBakt+∑hDah∗EHaht,(3)f2=min∑jxj∗Cj+∑kxk∗Ck+∑cxc∗Cc+∑t∑r∑c∑hCRr∗ωh∗Dch∗Rrcht+∑jCRr∗ωj∗Dcj∗Rrcjt+∑kCRr∗ωk∗Dck∗Rrckt+∑aCRr∗ωa∗Dca∗Rrcat,(4)f3=min∑t∑r∑hωh∗FRurht+∑jωj∗FRurjt+∑kωk∗FRurkt+∑aωa∗FRurat,is subject to(5)EAat=∑jEAajt,a∈AA,j∈A−ES,t∈T,(6)EBat=∑kEBakt,a∈AA,k∈B−ES,t∈T,(7)EIat=∑hEHaht,a∈AA,h∈EMC,t∈T,(8)∑tEAajt≤MAj∗xj,a∈AA,j∈A−ES,t∈T,(9)∑tEBakt≤MAk∗xk,a∈AA,k∈B−ES,t∈T,(10)∑tEHaht≤MAh,a∈AA,h∈EMC,t∈T,(11)∑t∑ΦRrcΦt≤Rrc∗xC,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,c∈EDC,t∈T,(12)EAjt=∑aEAajt,a∈AA,j∈A−ES,t=1,(13)EAjt=EAjt−1+∑aEAajt,a∈AA,j∈A−ES,t≥2,(14)EBkt=∑aEBakt,a∈AA,k∈B−ES,t=1,(15)EBkt=EBkt−1+∑aEBakt,a∈AA,k∈B−ES,t≥2,(16)EHht=∑aEHaht,a∈AA,h∈EMC,t=1,(17)EHht=EHht−1+∑aEHaht,a∈AA,h∈EMC,t≥2,(18)Eat=Ea−EAat−EBat−EIat,a∈AA,t=1,(19)Eat=Eat−1−EAat−EBat−EIat,a∈AA,t≥2,(20)Rrjt=δjr∗EAjt,r∈R,j∈A−ES,t∈T,(21)Rrkt=δkr∗EBkt,r∈R,k∈B−ES,t∈T,(22)Rrht=δhr∗EHht,r∈R,h∈EMC,t∈T,(23)Rrat=δar∗Eat,r∈R,a∈AA,t∈T,(24)RTr=∑t∑ΦRrΦt,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,t∈T,(25)RrΦt≥∑cRrcΦt,RrcΦt≥0,Φ∈AA∪A−ES∪B−ES∪EMC,c∈EDC,r∈R,t∈T,(26)∑cRrc∗xC=α∗RTr,c∈EDC,r∈R,(27)RurΦt=RrΦt−∑cRrcΦt,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,t∈T,c∈EDC,(28)FRurΦt=1−e−μ∗RurΦtμ,μ<0,r∈R,Φ∈AA∪A−ES∪B−ES∪EMC,t∈T,c∈EDC.The programming model involves three objective functions (minimization of the total evacuation distance, total cost, and total unmet resources), specified by equations (2)–(4). Constraints (5)–(7) make sure that all refugees requiring hierarchical evacuation in AAs can be evacuated into the corresponding nodes. The maximum capacities of A-ESs, B-ESs, EMCs, and EDCs are limited in constraints (8)–(11). Constraints (12)–(19) reflect the number of people in each node at time t. Constraints (20)–(23) define the quantity of required resources in each node at time t; thus, constraint (24) defines the total quantity of relief resources r required from all resource demand points throughout the rescue period of 72 h. To avoid the wastage of insufficient resources, we define that the quantity of received resources of each node is no more than the quantity of required resources of each node by equation (25). Considering the hierarchical resource demand of each node and insufficient resources of EDCs, to reduce the secondary damage caused by insufficient resources, constraint (26) sets the total quantity of resources in all EDCs (α≥1 means the total quantity of resources is sufficient). Constraint (27) indicates the quantity of unmet resources in each node. Constraint (28) defines the loss for the unmet resources of each node at time t, where the parameter μ is used to show the importance placed on large losses.


## 4. Solution Methodology

### 4.1. Algorithm Frameworks

In the previous section, we proposed the hierarchical multiobjective nonlinear mathematical model. Considering our research on site selection and resource allocation during the critical 72 h after disaster, the general commercial solver like CPLEX and LINGO cannot solve this complex nonlinear problem [5].The SA algorithm is the expansion of the local search algorithm. Different from the local search algorithm, SA accepts a worse solution with a certain probability in iterations to prevent local optimum [20]. Based on this optimization feature, SA has been used in the site selection research [21], manufacturing problem [22], and multiobjective problem [23]. SA can optimize to the global optimum with very slow annealing [24]. Therefore, the SA algorithm can be a good solution to solve the site selection of A-ESs, B-ESs, EMCs, and EDCs.The PSO algorithm is a stochastic evolutionary based on birds’ group behavior. Due to its fast convergence, high efficiency, and robustness, PSO has become a popular tool to solve complex problems in the social network [25] and the computer science [26]. The algorithm has also been used in solving large-scale emergency resource allocation problems and has proved to be time efficient [4]. In our research, PSO is an effective algorithm to solve the resource allocation problem.In this section, we use the combination of two heuristic algorithms to solve our programming model, and outer SA and inner PSO methods are applied to site selection and resource allocation. The flow of our combination algorithms is shown in Figure2. Next, we would explain how the SA and PSO methods are applied to our multiobjective problem with constraints.Figure 2
Flow diagram of combination algorithms.


### 4.2. Fundamental SA-PSO Algorithm

#### 4.2.1. SA Algorithm

The main steps of the SA algorithm are as follows:Step 1. Initialize the parameters: start temperature ts, end temperature te, decay rate of temperature ρ, receive tolerance for inferior solutions ε, initial solution s0, optimal solution unchanged times key, and maximum optimal solution unchanged times keymax.Step 2. Iterate L times in each temperature, generate a new neighbor s′ from the previous solution s based on certain rules in the current temperature tk, according to the Metropolis guidelines [27], and decide if s′ is accepted as the current solution with the following probability:(29)Paccepts′as current solution=exp−fs′−fsε∗tk,fs′−fs>0,1,fs′−fs≤0,wheretk is the current temperature of the kth iteration, tk>0, ∀k, and fs is the target value of solution s.Step 3. Outer loop: if tk is lower than te or key≥keymax, end the algorithm; else, if tk=ρ∗tk, return to step 2.In our combination algorithms of solving site selection and resource allocation problems, we get the optimal objective value through the outer SA algorithm. Thus, the objective functionfs of the SA algorithm is based on all objective functions f1, f2, and f3. By combining the multiobjective into a single-objective function, we use a weighted sum approach that sets the weights of different targets. Therefore, fs can be represented as the following equation, and ∑∇=13ω∇=1:(30)fs=minω1∗f1−f1minf1max−f1min+ω2∗f2−f2minf2max−f2min+ω3∗f3−f3minf3max−f3min.Minimum and maximum objective function values are shown in Table1.Table 1
Minimum/maximum value of each objective function.
Objective functionObjective valueMaximum valueMinimum valuef1(evacuation distance)8.92E + 064.37E + 06f2(cost)1.03E + 099.01E + 08f3(unmet resources)1.47E + 091.79E + 06The random initialization site selection of the A-ES, B-ES, and EMC solution is based on constraints (5)–(10), which ensure that all refugees are evacuated into corresponding nodes with limited capacity. And, the random initialization site selection of the EDC solution is based on constraint (26). The initial site selection solution s0 of the SA algorithm is shown in Figure 3, where xj=1, xk=1, and xc=1 mean the candidates A-ES, B-ES, and EMC are selected and xj=0, xk=0, and xc=0 mean the candidates A-ES, B-ES, and EMC are not selected.Figure 3
Initial site selection solution of the SA algorithm.The SA algorithm randomly selects a new solution from the neighborhood feasible solution spaceNs of the current solution s. In our site selection problem, the method of generating a new solution s′, s′∈Ns, can be divided into three types.Randomly adding an unselected point; randomly removing a selected point on the premise of meeting the above constraints; randomly changing the selected point (randomly adding an unselected point while randomly removing a selected point) on the premise of meeting the above constraints.Taking the location selection of the four candidate points as an example, assume that the initial solutions is as shown in Figure 4, where 1 means the candidate point is selected, 0 otherwise.Figure 4
Initial solution of the SA algorithm.Then, based on the previous method of generating a new solution, the neighborhood feasible solution spaceNs is shown in Figure 5.Figure 5
Neighborhood feasible solution space.


#### 4.2.2. PSO Algorithm

The main steps of the PSO algorithm are as follows:Step 1. Initialize the parameters: set N as the particle number, set psi0 as the initial position of the ith particle (i∈1,N), set vi0 as the initial velocity of the ith particle (i∈1,N), and calculate the fitness function fpsi0 of the initial position of the ith particle (i∈1,N).Step 2. Updating the velocity and position of the ith particle in the kth iteration with the following equations:(31)vik+1=w×vik+c1×r1×pbestik−psik+c2×r2×gbestk−psik,(32)psik+1=psik+vik+1.Step 3. Calculating the fitness function fpsik of the ith particle in the kth iteration (k∈1,K), personal best pbestik and global best gbestk are obtained, and if k>K, return gbestk and end the algorithm, else return to step 2.The inner PSO algorithm is to find the best solution to allocate insufficient resources, and the fitness functionfpsik only optimizes the resource allocation part. Thus, based on the objective function of f2 and f3, fpsik can be represented as follows:(33)fpsik=minω2′∗f2′−f′2minf′2max−f′2min+ω3′∗f3−f3minf3max−f3min,f2′=min∑t∑r∑c∑hCRr∗ωh∗Dch∗Rrcht+∑jCRr∗ωj∗Dcj∗Rrcjt+∑kCRr∗ωk∗Dck∗Rrckt+∑aCRr∗ωa∗Dca∗Rrcat.The random initialization resource allocation solution from EDCs to each resource demand points at timet is based on constraints (11) and (25). The initial PSO resource allocation solution of the ith particle (i∈1,N) psi0 can be represented by a four-dimensional array of all Rrcht, Rrcjt, Rrckt, and Rrcat in the model in Section 3.2. The updating method of PSO is based on equations (31) and (32). If psik does not satisfy constraints (11) and (25), Rrc∗t (∗∈EMCs, A-ESs, B-ESs, and AAs) will be set to 0.


### 4.3. Improved SA-PSO Algorithm

The SA algorithm has the advantages of jumping out of local optimal and fast convergence speed. But, its single individual evolution method cannot grasp the entire search space, which makes its search efficiency not high [20, 28]. GA is a global search algorithm based on the theory of evolution and the natural rule of “survival of the fittest” [29]. GA has been widely applied in real-world problems as it is very robust and has strong global search capabilities [30]. Based on the characteristics of SA and GA, we absorb the group search, crossover, and mutation operator of GA into SA to optimize and expand its population diversity. Thus, the outer site selection problem is solved by the ISA algorithm.


#### 4.3.1. GA Algorithm

The main steps of the GA algorithm are as follows:Step 1. Initialize the parameters: population size G, maximal number of generations M, crossover rate pc, and mutation rate pm, set sg0 as initial solution of the gth individual (g∈1,G), calculate the fitness function fsg0 of the gth individual, and find FitnessBest.Step 2. Use roulette wheel selection to select the new population, and use the randomly crossover and mutation to optimize and widen the solution space.Step 3. Calculate the fitness function fsgj of the gth individual in the jth iteration and find FitnessBest (g∈1,G, j∈1,M), and if j>M, return FitnessBest and end the algorithm, else return to step 2.


#### 4.3.2. ISA Algorithm

The main steps of the GA algorithm are as follows:Step 1. Initialize the parameters of SA: start temperature ts, end temperature te, decay rate of temperature ρ, receive tolerance for inferior solutions ε, optimal solution unchanged times key, and maximum optimal solution unchanged times keymax; initialize the parameters of GA: population size G, crossover rate pc, and mutation rate pm. Set sg0 as the initial solution of the gth individual (g∈1,G), calculate the fitness function fsg0 of the gth individual, and find FitnessBest.Step 2. Iterate L times in each temperature, use roulette wheel selection to select the new population Pop′ from prior population Pop, update the new population Pop' by randomly roulette wheel selection, crossover, and mutation, calculate the fitness function fsgl of the gth individual in the lth iteration g∈1,G,l∈1,L, and find FitnessBest′, FitnessBest′=minfsgl. According to the Metropolis guidelines, decide if FitnessBest′ is accepted as the best solution and if Pop′ is accepted as the new population with the following probability:(34)PFitnessBest=FitnessBest′,Pop=Pop′=exp−FitnessBest′−FitnessBestε∗tk,FitnessBest′>FitnessBest,1,FitnessBest′≤FitnessBest.Step 3. Outer loop: if tk is lower than te or key≥keymax, end the algorithm; else if tk=ρ∗tk, return to step 2.Based on the initial site selection and the update solution of the SA algorithm in Section4.2.1, the chromosome, which represents the initial site selection solution of the ISA algorithm, is the same as s0. Combining the characteristics of GA evolution, the update solution of ISA includes selection, crossover, and mutation operation.


## 4.1. Algorithm Frameworks

In the previous section, we proposed the hierarchical multiobjective nonlinear mathematical model. Considering our research on site selection and resource allocation during the critical 72 h after disaster, the general commercial solver like CPLEX and LINGO cannot solve this complex nonlinear problem [5].The SA algorithm is the expansion of the local search algorithm. Different from the local search algorithm, SA accepts a worse solution with a certain probability in iterations to prevent local optimum [20]. Based on this optimization feature, SA has been used in the site selection research [21], manufacturing problem [22], and multiobjective problem [23]. SA can optimize to the global optimum with very slow annealing [24]. Therefore, the SA algorithm can be a good solution to solve the site selection of A-ESs, B-ESs, EMCs, and EDCs.The PSO algorithm is a stochastic evolutionary based on birds’ group behavior. Due to its fast convergence, high efficiency, and robustness, PSO has become a popular tool to solve complex problems in the social network [25] and the computer science [26]. The algorithm has also been used in solving large-scale emergency resource allocation problems and has proved to be time efficient [4]. In our research, PSO is an effective algorithm to solve the resource allocation problem.In this section, we use the combination of two heuristic algorithms to solve our programming model, and outer SA and inner PSO methods are applied to site selection and resource allocation. The flow of our combination algorithms is shown in Figure2. Next, we would explain how the SA and PSO methods are applied to our multiobjective problem with constraints.Figure 2
Flow diagram of combination algorithms.


## 4.2. Fundamental SA-PSO Algorithm

### 4.2.1. SA Algorithm

The main steps of the SA algorithm are as follows:Step 1. Initialize the parameters: start temperature ts, end temperature te, decay rate of temperature ρ, receive tolerance for inferior solutions ε, initial solution s0, optimal solution unchanged times key, and maximum optimal solution unchanged times keymax.Step 2. Iterate L times in each temperature, generate a new neighbor s′ from the previous solution s based on certain rules in the current temperature tk, according to the Metropolis guidelines [27], and decide if s′ is accepted as the current solution with the following probability:(29)Paccepts′as current solution=exp−fs′−fsε∗tk,fs′−fs>0,1,fs′−fs≤0,wheretk is the current temperature of the kth iteration, tk>0, ∀k, and fs is the target value of solution s.Step 3. Outer loop: if tk is lower than te or key≥keymax, end the algorithm; else, if tk=ρ∗tk, return to step 2.In our combination algorithms of solving site selection and resource allocation problems, we get the optimal objective value through the outer SA algorithm. Thus, the objective functionfs of the SA algorithm is based on all objective functions f1, f2, and f3. By combining the multiobjective into a single-objective function, we use a weighted sum approach that sets the weights of different targets. Therefore, fs can be represented as the following equation, and ∑∇=13ω∇=1:(30)fs=minω1∗f1−f1minf1max−f1min+ω2∗f2−f2minf2max−f2min+ω3∗f3−f3minf3max−f3min.Minimum and maximum objective function values are shown in Table1.Table 1
Minimum/maximum value of each objective function.
Objective functionObjective valueMaximum valueMinimum valuef1(evacuation distance)8.92E + 064.37E + 06f2(cost)1.03E + 099.01E + 08f3(unmet resources)1.47E + 091.79E + 06The random initialization site selection of the A-ES, B-ES, and EMC solution is based on constraints (5)–(10), which ensure that all refugees are evacuated into corresponding nodes with limited capacity. And, the random initialization site selection of the EDC solution is based on constraint (26). The initial site selection solution s0 of the SA algorithm is shown in Figure 3, where xj=1, xk=1, and xc=1 mean the candidates A-ES, B-ES, and EMC are selected and xj=0, xk=0, and xc=0 mean the candidates A-ES, B-ES, and EMC are not selected.Figure 3
Initial site selection solution of the SA algorithm.The SA algorithm randomly selects a new solution from the neighborhood feasible solution spaceNs of the current solution s. In our site selection problem, the method of generating a new solution s′, s′∈Ns, can be divided into three types.Randomly adding an unselected point; randomly removing a selected point on the premise of meeting the above constraints; randomly changing the selected point (randomly adding an unselected point while randomly removing a selected point) on the premise of meeting the above constraints.Taking the location selection of the four candidate points as an example, assume that the initial solutions is as shown in Figure 4, where 1 means the candidate point is selected, 0 otherwise.Figure 4
Initial solution of the SA algorithm.Then, based on the previous method of generating a new solution, the neighborhood feasible solution spaceNs is shown in Figure 5.Figure 5
Neighborhood feasible solution space.


### 4.2.2. PSO Algorithm

The main steps of the PSO algorithm are as follows:Step 1. Initialize the parameters: set N as the particle number, set psi0 as the initial position of the ith particle (i∈1,N), set vi0 as the initial velocity of the ith particle (i∈1,N), and calculate the fitness function fpsi0 of the initial position of the ith particle (i∈1,N).Step 2. Updating the velocity and position of the ith particle in the kth iteration with the following equations:(31)vik+1=w×vik+c1×r1×pbestik−psik+c2×r2×gbestk−psik,(32)psik+1=psik+vik+1.Step 3. Calculating the fitness function fpsik of the ith particle in the kth iteration (k∈1,K), personal best pbestik and global best gbestk are obtained, and if k>K, return gbestk and end the algorithm, else return to step 2.The inner PSO algorithm is to find the best solution to allocate insufficient resources, and the fitness functionfpsik only optimizes the resource allocation part. Thus, based on the objective function of f2 and f3, fpsik can be represented as follows:(33)fpsik=minω2′∗f2′−f′2minf′2max−f′2min+ω3′∗f3−f3minf3max−f3min,f2′=min∑t∑r∑c∑hCRr∗ωh∗Dch∗Rrcht+∑jCRr∗ωj∗Dcj∗Rrcjt+∑kCRr∗ωk∗Dck∗Rrckt+∑aCRr∗ωa∗Dca∗Rrcat.The random initialization resource allocation solution from EDCs to each resource demand points at timet is based on constraints (11) and (25). The initial PSO resource allocation solution of the ith particle (i∈1,N) psi0 can be represented by a four-dimensional array of all Rrcht, Rrcjt, Rrckt, and Rrcat in the model in Section 3.2. The updating method of PSO is based on equations (31) and (32). If psik does not satisfy constraints (11) and (25), Rrc∗t (∗∈EMCs, A-ESs, B-ESs, and AAs) will be set to 0.


## 4.2.1. SA Algorithm

The main steps of the SA algorithm are as follows:Step 1. Initialize the parameters: start temperature ts, end temperature te, decay rate of temperature ρ, receive tolerance for inferior solutions ε, initial solution s0, optimal solution unchanged times key, and maximum optimal solution unchanged times keymax.Step 2. Iterate L times in each temperature, generate a new neighbor s′ from the previous solution s based on certain rules in the current temperature tk, according to the Metropolis guidelines [27], and decide if s′ is accepted as the current solution with the following probability:(29)Paccepts′as current solution=exp−fs′−fsε∗tk,fs′−fs>0,1,fs′−fs≤0,wheretk is the current temperature of the kth iteration, tk>0, ∀k, and fs is the target value of solution s.Step 3. Outer loop: if tk is lower than te or key≥keymax, end the algorithm; else, if tk=ρ∗tk, return to step 2.In our combination algorithms of solving site selection and resource allocation problems, we get the optimal objective value through the outer SA algorithm. Thus, the objective functionfs of the SA algorithm is based on all objective functions f1, f2, and f3. By combining the multiobjective into a single-objective function, we use a weighted sum approach that sets the weights of different targets. Therefore, fs can be represented as the following equation, and ∑∇=13ω∇=1:(30)fs=minω1∗f1−f1minf1max−f1min+ω2∗f2−f2minf2max−f2min+ω3∗f3−f3minf3max−f3min.Minimum and maximum objective function values are shown in Table1.Table 1
Minimum/maximum value of each objective function.
Objective functionObjective valueMaximum valueMinimum valuef1(evacuation distance)8.92E + 064.37E + 06f2(cost)1.03E + 099.01E + 08f3(unmet resources)1.47E + 091.79E + 06The random initialization site selection of the A-ES, B-ES, and EMC solution is based on constraints (5)–(10), which ensure that all refugees are evacuated into corresponding nodes with limited capacity. And, the random initialization site selection of the EDC solution is based on constraint (26). The initial site selection solution s0 of the SA algorithm is shown in Figure 3, where xj=1, xk=1, and xc=1 mean the candidates A-ES, B-ES, and EMC are selected and xj=0, xk=0, and xc=0 mean the candidates A-ES, B-ES, and EMC are not selected.Figure 3
Initial site selection solution of the SA algorithm.The SA algorithm randomly selects a new solution from the neighborhood feasible solution spaceNs of the current solution s. In our site selection problem, the method of generating a new solution s′, s′∈Ns, can be divided into three types.Randomly adding an unselected point; randomly removing a selected point on the premise of meeting the above constraints; randomly changing the selected point (randomly adding an unselected point while randomly removing a selected point) on the premise of meeting the above constraints.Taking the location selection of the four candidate points as an example, assume that the initial solutions is as shown in Figure 4, where 1 means the candidate point is selected, 0 otherwise.Figure 4
Initial solution of the SA algorithm.Then, based on the previous method of generating a new solution, the neighborhood feasible solution spaceNs is shown in Figure 5.Figure 5
Neighborhood feasible solution space.


## 4.2.2. PSO Algorithm

The main steps of the PSO algorithm are as follows:Step 1. Initialize the parameters: set N as the particle number, set psi0 as the initial position of the ith particle (i∈1,N), set vi0 as the initial velocity of the ith particle (i∈1,N), and calculate the fitness function fpsi0 of the initial position of the ith particle (i∈1,N).Step 2. Updating the velocity and position of the ith particle in the kth iteration with the following equations:(31)vik+1=w×vik+c1×r1×pbestik−psik+c2×r2×gbestk−psik,(32)psik+1=psik+vik+1.Step 3. Calculating the fitness function fpsik of the ith particle in the kth iteration (k∈1,K), personal best pbestik and global best gbestk are obtained, and if k>K, return gbestk and end the algorithm, else return to step 2.The inner PSO algorithm is to find the best solution to allocate insufficient resources, and the fitness functionfpsik only optimizes the resource allocation part. Thus, based on the objective function of f2 and f3, fpsik can be represented as follows:(33)fpsik=minω2′∗f2′−f′2minf′2max−f′2min+ω3′∗f3−f3minf3max−f3min,f2′=min∑t∑r∑c∑hCRr∗ωh∗Dch∗Rrcht+∑jCRr∗ωj∗Dcj∗Rrcjt+∑kCRr∗ωk∗Dck∗Rrckt+∑aCRr∗ωa∗Dca∗Rrcat.The random initialization resource allocation solution from EDCs to each resource demand points at timet is based on constraints (11) and (25). The initial PSO resource allocation solution of the ith particle (i∈1,N) psi0 can be represented by a four-dimensional array of all Rrcht, Rrcjt, Rrckt, and Rrcat in the model in Section 3.2. The updating method of PSO is based on equations (31) and (32). If psik does not satisfy constraints (11) and (25), Rrc∗t (∗∈EMCs, A-ESs, B-ESs, and AAs) will be set to 0.


## 4.3. Improved SA-PSO Algorithm

The SA algorithm has the advantages of jumping out of local optimal and fast convergence speed. But, its single individual evolution method cannot grasp the entire search space, which makes its search efficiency not high [20, 28]. GA is a global search algorithm based on the theory of evolution and the natural rule of “survival of the fittest” [29]. GA has been widely applied in real-world problems as it is very robust and has strong global search capabilities [30]. Based on the characteristics of SA and GA, we absorb the group search, crossover, and mutation operator of GA into SA to optimize and expand its population diversity. Thus, the outer site selection problem is solved by the ISA algorithm.


### 4.3.1. GA Algorithm

The main steps of the GA algorithm are as follows:Step 1. Initialize the parameters: population size G, maximal number of generations M, crossover rate pc, and mutation rate pm, set sg0 as initial solution of the gth individual (g∈1,G), calculate the fitness function fsg0 of the gth individual, and find FitnessBest.Step 2. Use roulette wheel selection to select the new population, and use the randomly crossover and mutation to optimize and widen the solution space.Step 3. Calculate the fitness function fsgj of the gth individual in the jth iteration and find FitnessBest (g∈1,G, j∈1,M), and if j>M, return FitnessBest and end the algorithm, else return to step 2.


### 4.3.2. ISA Algorithm

The main steps of the GA algorithm are as follows:Step 1. Initialize the parameters of SA: start temperature ts, end temperature te, decay rate of temperature ρ, receive tolerance for inferior solutions ε, optimal solution unchanged times key, and maximum optimal solution unchanged times keymax; initialize the parameters of GA: population size G, crossover rate pc, and mutation rate pm. Set sg0 as the initial solution of the gth individual (g∈1,G), calculate the fitness function fsg0 of the gth individual, and find FitnessBest.Step 2. Iterate L times in each temperature, use roulette wheel selection to select the new population Pop′ from prior population Pop, update the new population Pop' by randomly roulette wheel selection, crossover, and mutation, calculate the fitness function fsgl of the gth individual in the lth iteration g∈1,G,l∈1,L, and find FitnessBest′, FitnessBest′=minfsgl. According to the Metropolis guidelines, decide if FitnessBest′ is accepted as the best solution and if Pop′ is accepted as the new population with the following probability:(34)PFitnessBest=FitnessBest′,Pop=Pop′=exp−FitnessBest′−FitnessBestε∗tk,FitnessBest′>FitnessBest,1,FitnessBest′≤FitnessBest.Step 3. Outer loop: if tk is lower than te or key≥keymax, end the algorithm; else if tk=ρ∗tk, return to step 2.Based on the initial site selection and the update solution of the SA algorithm in Section4.2.1, the chromosome, which represents the initial site selection solution of the ISA algorithm, is the same as s0. Combining the characteristics of GA evolution, the update solution of ISA includes selection, crossover, and mutation operation.


## 4.3.1. GA Algorithm

The main steps of the GA algorithm are as follows:Step 1. Initialize the parameters: population size G, maximal number of generations M, crossover rate pc, and mutation rate pm, set sg0 as initial solution of the gth individual (g∈1,G), calculate the fitness function fsg0 of the gth individual, and find FitnessBest.Step 2. Use roulette wheel selection to select the new population, and use the randomly crossover and mutation to optimize and widen the solution space.Step 3. Calculate the fitness function fsgj of the gth individual in the jth iteration and find FitnessBest (g∈1,G, j∈1,M), and if j>M, return FitnessBest and end the algorithm, else return to step 2.


## 4.3.2. ISA Algorithm

The main steps of the GA algorithm are as follows:Step 1. Initialize the parameters of SA: start temperature ts, end temperature te, decay rate of temperature ρ, receive tolerance for inferior solutions ε, optimal solution unchanged times key, and maximum optimal solution unchanged times keymax; initialize the parameters of GA: population size G, crossover rate pc, and mutation rate pm. Set sg0 as the initial solution of the gth individual (g∈1,G), calculate the fitness function fsg0 of the gth individual, and find FitnessBest.Step 2. Iterate L times in each temperature, use roulette wheel selection to select the new population Pop′ from prior population Pop, update the new population Pop' by randomly roulette wheel selection, crossover, and mutation, calculate the fitness function fsgl of the gth individual in the lth iteration g∈1,G,l∈1,L, and find FitnessBest′, FitnessBest′=minfsgl. According to the Metropolis guidelines, decide if FitnessBest′ is accepted as the best solution and if Pop′ is accepted as the new population with the following probability:(34)PFitnessBest=FitnessBest′,Pop=Pop′=exp−FitnessBest′−FitnessBestε∗tk,FitnessBest′>FitnessBest,1,FitnessBest′≤FitnessBest.Step 3. Outer loop: if tk is lower than te or key≥keymax, end the algorithm; else if tk=ρ∗tk, return to step 2.Based on the initial site selection and the update solution of the SA algorithm in Section4.2.1, the chromosome, which represents the initial site selection solution of the ISA algorithm, is the same as s0. Combining the characteristics of GA evolution, the update solution of ISA includes selection, crossover, and mutation operation.


## 5. Case Study

In this section, we intend to (1) present a case study to illustrate the efficiency of our mathematical model to realize quick location selection and resource allocation during disasters; (2) consider some scenario analyses for further investigation on emergency refugee evacuation and resource allocation decisions; (3) analyze the trade-off between different factors through the result of the computational experiments.Beijing is located in the northern part of the North China Plain, with a latitude of39°26′N to 41°03′N and a longitude of 115°25′E to 117°30′E, mainly in the North China earthquake zone Yin-middle Yanshan seismic zone. The earthquakes with a magnitude greater than 4 in the Beijing area have occurred nearly 200 times, and more than 10 earthquakes with magnitude 5 have occurred. There is a possibility of a moderately devastating earthquake in the Beijing area, which makes the situation of earthquake prevention and mitigation faced by Beijing very serious. Chaoyang district, affiliated to Beijing, located in the eastern part of Beijing, covers an area of 470.8 km2 and is the largest district in the Beijing central city. In recent years, some emergency shelters have been built in Beijing, in which the Chaoyang district has the largest emergency shelters. We choose the Chaoyang district (Figure 6) as our study area to demonstrate the effectiveness and efficiency of our proposed model and algorithms.Figure 6
The map of the Chaoyang district.


### 5.1. Data Processing

We divide the first 3-day gold rescue time after a disaster into 6 hours of time intervals. Thus, the emergency evacuation and resource allocation network are designed in 12 (=72/6) intervals.


#### 5.1.1. Facility Data

According to the database of the Beijing Earthquake Agency, the Chaoyang district has built fifteen earthquake disaster shelters. However, these shelters may not accommodate so many refugees or cause more evacuation distances. Thus, based on the Beijing Central City Earthquake and Emergency Refuge Site Planning Outline and the Chaoyang Government Network, another sixteen candidate shelters were selected. The EMCs need professional medical equipment and nursing staffs. Thus, in our research, we select 26 candidate EMCs from the hospitals located in the Chaoyang district, as shown in Figure7. Since there are no EDCs in the Chaoyang District, Beijing, we assume that candidate EDCs are located at street/district offices. So, there are 42 candidate EDCs.Figure 7
Distribution map of the Chaoyang district.


#### 5.1.2. Population

Chaoyang district governs 24 streets and 19 districts, and the population size of each street/districtEa is extracted from the Chaoyang District Beijing Statistical Yearbook 2017. The number of refugees in each AA of each time interval Eat is roughly estimated [17], and Figure 8 presents the ratio of Eat to Ea. We assume that 10% of the population is injured EIat, 20% of the population is assumed to be A-priority refugees EAat, and 30% is assumed to be B-priority refugees EBat.Figure 8
Number of refugees in AAs over time.


#### 5.1.3. Demand

Quantity of required resources of different refugees is estimated by previous research [17]. Each refugee requires two types of relief resources: commodities and medical supplies. Commodities per unit are herein 1 kg, and the medical supplies per unit are 0.5 kg. Different refugees have different quantity of demand for different resources, and the specific resource requirements are presented in Table 2.Table 2
Quantity of required resources of hierarchical refugees.
δjrδkrδhrδhrr=111.50.52r=21.5120.5


#### 5.1.4. Capacity

We assume that each ES is divided into two parts: A-ES and B-ES, where the capacity of A-ESs is smaller than B-ES. The capacities of A-ES, B-ES, and EMCs are presented in Tables3 and 4.Table 3
Capacity of EMCs.
EMCsMAhCapacityBeijing Obstetrics and Gynecology Hospital, Capital Medical Universityh16600Beijing Obstetrics and Gynecology Hospital, Capital Medical Universityh28000Beijing Chaoyang Hospital, Capital Medical Universityh319000Wangjing Hospital, China Academy of Chinese Medical Sciencesh413000China-Japan Friendship Hospital of the Ministry of Healthh515000People’s Liberation Army No. 306 Hospitalh66000Beijing Anzhen Hospital, Capital Medical Universityh715000Beijing First Hospital of Integrated Traditional Chinese and Western Medicineh816000Beijing Ditan Hospital, Capital Medical Universityh99000Chinese Academy of Medical Sciences Cancer Hospitalh1016000Children’s Hospital of the Capital Institute of Pediatricsh116000Beijing Dangyangliu Hospitalh125600Civil Aviation General Hospitalh136500China Medical University Aviation General Hospitalh1412000Chinese People’s Armed Police Force Beijing Corps Hospitalh1515000Beijing Huaxin Hospitalh167600Coal General Hospitalh175000The Third Affiliated Hospital of Beijing University of Chinese Medicineh185200Beijing Chaoyang District Second Hospitalh194500Shuangqiao Hospital, Chaoyang District, Beijingh205500Beijing Coking Chemical Factory Hospitalh213000Chaoyang District Chinese Medicine Hospital, Beijingh226500Beijing Chaoyang District Maternal and Child Health Hospitalh234500Beijing Tibetan Hospitalh241600Beijing Wuzhou Women and Children Hospitalh255000Beijing Dangyangliu Hospital East Campush264500Table 4
Capacity of ESs.
ESsMAjCapacityMAkCapacityYuan Dadu Archaeological Site Parkj190000k1100000Chaoyang Parkj268000k2182000Sun Palace Parkj320000k390000Olympic Forest Park South Parkj47000k413000Anzhen Yongxi Parkj52600k53000Taiba River Green Beltj630000k650000Xinglong Parkj732800k740000Red Scarf parkj816100k818000Beixiaohe Parkj915700k920000Beijing Liyuanj1050000k1055400Wanghe Parkj1139200k1140000Cuicheng Parkj121950k124000Hongbo Country Parkj13100000k13140000West Dawang Road Community Parkj147200k147000Cuicheng Parkj121950k124000Hongbo Country Parkj13100000k13140000West Dawang Road Community Parkj147200k147000Lishuiqiao Parkj155000k157000Dongba Country Parkj1655000k1680000White Deer Parkj1725000k1740000Jintian Parkj1860000k1890000Dongfeng Parkj1910000k1930000Haitang Parkj2024000k2050000Laojuntang Parkj2148000k2160000Duzhong Parkj2255000k2280000Changying Parkj2377000k2390000Chaolai Forest Parkj2444000k2470000Guta Parkj2555000k2570000Sun River Country Parkj2614000k2630000Dawangjingj2735000k2740000Red Army Parkj3182500k3190000Because of the lack of EDC data in Chaoyang District, Beijing, combining the total demand of resources of all refugees in each node, we assume here that 9 EDCs were selected from 42 candidate EDCs based on previous EDC selection research studies [4, 8, 31–33]. The capacity of each EDC can be derived from the total number of EDCs, total quantity of required resources, and the ratio α.


#### 5.1.5. Distance

We obtain the road networks from the map of Beijing published by the China Cartographic Publishing House in 2016. We calculate the actual distances from AAs to A-ESs, B-ESs, and EMCs and from EDCs to AAs, A-ESs, B-ESs, and EMCs using version 10.1 of ArcGIS.


#### 5.1.6. Cost

The total cost involves the construction cost of candidate A-ESs, B-ESs (existing shelters have no construction cost), and EDCs and the fixed and transportation cost of per unit resourcer. The construction cost of each site is determined by its capacity, which is estimated by the government report and previous research studies [5, 18, 19]. The fixed cost of resources is also estimated by previous research studies, and the transportation cost of resources is estimated by the Beijing logistic price. All these parameters are shown in Table 5.Table 5
Cost parameters.
ParameterEstimationConstruction costCj=40; Ck=15; Cc=15Fixed costCRr=10r=1; CRr=20r=2Transportation costCTr=1


### 5.2. Model Testing and Results

In this section, the proposed model and algorithms are tested with the abovementioned parameter settings. As the priority of refugees in different places is EMCs> A-ESs > B-ESs > AAs and ∑∗ω∗=1 (∗∈EMCs, A-ESs, B-ESs, and AAs), we set the priorities of different places ωh, ωj, ωk, and ωa as following. We consider that the quantity of resources in all EDCs is insufficient and higher priority places have a higher minimum resource satisfaction degree. Different satisfaction degree settings are based on the ratio of total resources in DCs to total required resources. α=0.9 means that the total resources of the emergency distribution center account for 90% of the total resources needed for the entire rescue cycle. μ is a measure of the fairness loss function, and μ<0. When the amount of unsatisfied resources increases, the fairness loss function grows like an exponential function as μ decreases and grows like a linear function as μ increases. For parameters of SA, ω1, ω2, and ω3 are the weights of the three objective functions in the model, and also the weights of the three subobjective functions in the fitness function of the outer simulated annealing algorithm. ts is the start temperature, te is the end temperature, ρ is the decay rate of temperature, ε is the received tolerance for inferior solutions, key is the optimal solution unchanged times, keymax is the maximum optimal solution unchanged times, and L is the iterate times in each temperature. For parameters of PSO, K is the maximum iteration times, N is the particle number, c1 and c2 are the learning factors, c1 indicates the possibility of the particle updating according to its local optimal position, c2 indicates the possibility of the particle updating according to the global optimal position, w is the inertia factor, which reflects the particle’s ability to inherit the previous speed, and ω2′ and ω3′ are the weight setting of the objective functions f2′ and f3. For parameters of GA, G represents the number of groups, pc is the crossover rate, and pm is the mutation rate. According to relevant research [20–22, 25, 26, 29, 30], parameters of the model and algorithms are also shown in Table 6.Table 6
Parameter setting.
ParameterEstimationParameters of the modelωh=0.32;ωj=0.28;ωk=0.22;ωa=0.1;α=0.9;μ=−0.0001Parameters of SAts=100;te=1;ρ=0.9;L=10;ε=0.0002;ω1=1/3;ω2=1/3;ω3=1/3;key=0;keymax=50Parameters of PSON=20;K=100;c1=0.9;c2=0.9;w=0.8;ω2′=0.5;ω3′=0.5Parameters of GAG=50;pc=0.8;pm=0.2By setting the weight of each objective function equal to 1, with the objective of minimizing/maximizingfs of SA, we calculate the minimum/maximum value of each objective of our two test cases (as shown in Table 1).We firstly test the case with these existing shelters as case 1, in which the construction cost of A-ESs and B-ESs set to be zero. Then, these candidate shelters are added in the second test case as case 2, using the parameters shown in Table6.Through 50 experiments, the average value of objective functions of the two test cases is shown in Table7.Table 7
Results of the two test cases.
Objective functionCase 1Case 2Dividend (%)f1(evacuation distance)8.92E + 065.27E + 06−40.92f2(cost)9.05E + 089.40E + 083.85f3(unmet resources)2.13E + 061.97E + 06−7.51The dividend is calculated by the equation=fcase1/fcase2−1, as shown in Table 7, and the results reveal that the total evacuation distance in case 2 decreases by 40.92%, while the total cost increases by 3.85% with selected candidate shelters compared to case 1. And, we see that the total unmet resources in case 2 also decreases. Thus, it is obvious that selecting candidate shelters can reduce the total evacuation distance to a large extent and also reduce the total unmet resources while only adding a small amount of cost.Taking the corresponding refugee evacuation in which the objective function value is the closest to the average objective function value as an example, the allocation results of refugees from AAs to A-ESs and B-ESs at timet=1 of the two cases are shown in Figure 9.Figure 9
Comparison of the two cases of refugee evacuation at timet=1. (a) Refugee evacuation with only the existing shelters. (b) Refugee evacuation with candidate shelters.
(a)(b)One can see from Figure9(a) that many AAs have no nearby shelter, leading to the relatively long evacuation distance of refugees in these AAs. However, as is shown in Figure 9(b), six A-ESs and six B-ESs were selected at time t = 1 (one more A-ES and B-ES were selected at time t=2) and optimized the long-distance allocation result. Through the comparison of Figures 9(a) and 9(b), it is obvious that selecting candidate shelters can greatly reduce the evacuation distance of some AAs.


### 5.3. Scenario Analyses

In this section, five scenarios with added candidate shelters are designed to compare the characteristics of our proposed model and change in disaster situation of refugees. As can be seen from the results of Section5.2, selecting candidate shelters is a more efficient planning method. Therefore, all the scenarios in this section are based on case 2.Scenario 1.
The first scenario is to compare the efficiency of SA-PSO and ISA-PSO. Parameters of the GA are shown in Table8. We analyze these two algorithms by comparing the results of case 2 under SA-PSO and ISA-PSO. The operational iteration graph of the two algorithms is shown in Figure 10, and through 50 experiments, the average value of three objective functions is presented in Table 8.
It can be seen from Figure10 that the iterative result of ISA-PSO is always better than that of SA-PSO. Table 9 presents the objective values calculated by ISA-PSO. The dividend is calculated by the equation =fscenario1/fcase2−1, and it can be seen that compared with SA-PSO, the objective function values by ISA-PSO can achieve Pareto optimization. But, because of the characteristics of population search of ISA-PSO, the time consumption of ISA-PSO is correspondingly longer.
The algorithm used in the following calculations will use ISA-PSO. And, the following comparison of objective values will be based on the results of ISA-PSO.Table 8
Results of Scenario 1.
Objective functionScenario 1Dividend (%)f1(evacuation distance)5.03E + 06−4.55f2(cost)9.34E + 08−0.62f3(unmet resources)1.88E + 06−4.57Figure 10
Iteration graph of SA-PSO and ISA-PSO.Table 9
Results of Scenario 2.
Objective functionS2-1S2-2S2-3f1(evacuation distance)4.46E + 066.08E + 064.43E + 06f2(cost)9.70E + 089.28E + 089.79E + 08f3(unmet resources)1.87E + 061.97E + 061.87E + 06Scenario 2.
The second scenario is to analyze the impact of the objective weight. Thus, we set two of the objective weights to 50% of the original and the other target weight to twice the original: S2-1:ω1=1/6,ω2=1/6,andω3=2/3; S2-2: ω1=1/6,ω2=2/3,andω3=1/6; S1-3: ω1=2/3,ω2=1/6,andω3=1/6. Through 50 experiments, the average value of three objective functions is presented in Table 9.
It can be seen from Table9 that f1 is greatly affected by the objective weight. When the weight increases, the objective value is optimized, while the objective value of f2 and f3 is affected by the weight to a small extent. In addition, it can be seen that f1 and f3 are negatively correlated with f2, and f2 increases as f1 and f3 decrease. The corresponding dividend comparison of three objective values of Scenario 2 is shown in Figure 11.
The dividend is calculated by the equation=fS2−∗/fscenario1−1∗∈1,2,3, and it can be inferred that f1 is the most sensitive to the objective weight. Therefore, provided that the budget is sufficient, increasing the objective weight of f1 appropriately can better optimize the result.Figure 11
Dividend objective values of Scenario 1.Scenario 3.
The third scenario is to analyze the effect of parameterα on three objective functions. The purpose of considering this scenario is to enable managers to improve the efficiency of resource allocation under uncertain or insufficient conditions. In our research, we assume the quantity of relief resources is insufficient, and here we change the ratio α as S3-1: α=0.75, S3-2: α=0.85, and S3-3: α=0.95, respectively. Through 50 experiments, the average value of three objective functions of Scenario 3 is presented in Table 10.
Table10 presents that the change in α has little effect on f1 and f2. It implies that resources and their transportation costs account for only a small part of the total cost (f2). And, obviously, as α increases, f3 decreases. Figure 12 illustrates the corresponding dividend comparison of Scenario 2. The dividend is calculated by the equation =fS3−∗/fscenario1−1∗∈1,2,3. It is observed that as α increases, f3 decreases rapidly and there were almost no change in f1 and f2. This scenario illustrates that adjusting the parameter α significantly changes objective values of f3 in a linear relationship, but does not affect the other objects. This implies that increasing the quantity of resources is an effective method that improves the target value of f3. It can be seen that preparing sufficient resources in advance can play a vital role in the disaster relief process at a small cost.Table 10
Results of Scenario 3.
Objective functionS2-1S2-2S2-3f1(evacuation distance)5.13E + 064.98E + 064.74E + 06f2(cost)9.23E + 089.41E + 089.65E + 08f3(unmet resources)7.89E + 063.35E + 066.81E + 05Figure 12
Dividend objective values of Scenario 3.Scenario 4.
This scenario is to analyze the impact of population in each priority. Here, we increase the number of A-priority and B-priority refugees to be 30% and 40% of the population in S4-1 and 40% and 50% of the population in S4-2.
In this scenario, the total capacity of the existing shelters cannot meet the refugee demand. Thus, the selecting of candidate shelters is essential when the refugee demand increases.
Through 50 experiments, the average value of three objective functions of Scenario 4 is presented in Table11. The dividend is calculated by the equation =fs4−∗/fscenario1−1∗∈1,2. It is obvious that when the number of A- and B-priority refugees increases, the total evacuation distance increases evidently.
Taking the corresponding refugee evacuation in which the objective function value is the closest to the average objective function value as an example, the allocation results of refugees from AAs to A-ESs and B-ESs att = 1 of Scenario 4 are shown in Figure 13. Combining the results in Table 9, all these imply that the construction cost has only very small impact on the total cost f2. Besides, when the number of refugees who need to evacuate increases, the total evacuation distance will also increase even if more candidate shelters are selected. Furthermore, when the evacuation demand of A- and B-priority refugees increases, the number of refugees in A-ESs and B-ESs increases as well. This result demonstrates that refugees in A-ESs and B-ESs are more sensitive to unmet resources than those in AAs.Table 11
Results of Scenario 4.
Objective functionS4-1Dividend (%)S4-2Dividend (%)f1(evacuation distance)7.69E + 0655.049.90E + 0699.50f2(cost)9.51E + 082.429.82E + 085.86f3(unmet resources)2.13E + 0611.262.61E + 0636.65Figure 13
Comparison of the refugee evacuation results of Scenario 4 at timet=1. (a) Refugee evacuation plan of S4-1. (b) Refugee evacuation plan of S4-2.
(a)(b)Scenario 5.
In this scenario, we set the population of each period to be stable and compare the results with time-varying demand and stable demand. Figure14 displays the affected population in each time interval with a stable refuge demand. Through 50 calculations, the average value of three objective functions of Scenario 5 is presented in Table 12. The dividend is calculated by the equation =fscenario5/fscenario1−1. The results in Table 12 show that the total evacuation distance increment is 14.8% higher, while the total cost is only 0.54% less than the result under time-varying demand. The reason of the evacuation distance dividend is that the remaining capacity of A-ESs and B-ESs becomes smaller as time goes by, but the number of A- and B-priority refugees has not become smaller. Therefore, refugees in the later period will have to evacuate to farther A-ESs and B-ESs.
Similarly, most of the resources in EDCs were allocated to the AAs in the early period. The remaining resources in EDCs are getting less and less as time goes by. However, the demand for resources of A-ESs and B-ESs is increasing much faster than the time-varying demand model. Therefore, the resource requirements of A-ESs and B-ESs may not be met in the later period. Recall the implication of Scenario 4, and this might be the reason that the objective value of unmet resources is also 41.88% higher than the results of the time-varying demand. Thus, it is obvious that the time-varying demand estimation is a more efficient evacuation and allocation planning model because it reduces evacuation distance and unmet resources. Furthermore, it fits the actual refuge demand better.Figure 14
Unvarying demand of refugees in AAs over time.Table 12
Results of Scenario 5.
Objective functionS5-1Dividend (%)f1(evacuation distance)5.69E + 0614.80f2(cost)9.23E + 08−0.54f3(unmet resources)2.71E + 0641.88


## 5.1. Data Processing

We divide the first 3-day gold rescue time after a disaster into 6 hours of time intervals. Thus, the emergency evacuation and resource allocation network are designed in 12 (=72/6) intervals.


### 5.1.1. Facility Data

According to the database of the Beijing Earthquake Agency, the Chaoyang district has built fifteen earthquake disaster shelters. However, these shelters may not accommodate so many refugees or cause more evacuation distances. Thus, based on the Beijing Central City Earthquake and Emergency Refuge Site Planning Outline and the Chaoyang Government Network, another sixteen candidate shelters were selected. The EMCs need professional medical equipment and nursing staffs. Thus, in our research, we select 26 candidate EMCs from the hospitals located in the Chaoyang district, as shown in Figure7. Since there are no EDCs in the Chaoyang District, Beijing, we assume that candidate EDCs are located at street/district offices. So, there are 42 candidate EDCs.Figure 7
Distribution map of the Chaoyang district.


### 5.1.2. Population

Chaoyang district governs 24 streets and 19 districts, and the population size of each street/districtEa is extracted from the Chaoyang District Beijing Statistical Yearbook 2017. The number of refugees in each AA of each time interval Eat is roughly estimated [17], and Figure 8 presents the ratio of Eat to Ea. We assume that 10% of the population is injured EIat, 20% of the population is assumed to be A-priority refugees EAat, and 30% is assumed to be B-priority refugees EBat.Figure 8
Number of refugees in AAs over time.


### 5.1.3. Demand

Quantity of required resources of different refugees is estimated by previous research [17]. Each refugee requires two types of relief resources: commodities and medical supplies. Commodities per unit are herein 1 kg, and the medical supplies per unit are 0.5 kg. Different refugees have different quantity of demand for different resources, and the specific resource requirements are presented in Table 2.Table 2
Quantity of required resources of hierarchical refugees.
δjrδkrδhrδhrr=111.50.52r=21.5120.5


### 5.1.4. Capacity

We assume that each ES is divided into two parts: A-ES and B-ES, where the capacity of A-ESs is smaller than B-ES. The capacities of A-ES, B-ES, and EMCs are presented in Tables3 and 4.Table 3
Capacity of EMCs.
EMCsMAhCapacityBeijing Obstetrics and Gynecology Hospital, Capital Medical Universityh16600Beijing Obstetrics and Gynecology Hospital, Capital Medical Universityh28000Beijing Chaoyang Hospital, Capital Medical Universityh319000Wangjing Hospital, China Academy of Chinese Medical Sciencesh413000China-Japan Friendship Hospital of the Ministry of Healthh515000People’s Liberation Army No. 306 Hospitalh66000Beijing Anzhen Hospital, Capital Medical Universityh715000Beijing First Hospital of Integrated Traditional Chinese and Western Medicineh816000Beijing Ditan Hospital, Capital Medical Universityh99000Chinese Academy of Medical Sciences Cancer Hospitalh1016000Children’s Hospital of the Capital Institute of Pediatricsh116000Beijing Dangyangliu Hospitalh125600Civil Aviation General Hospitalh136500China Medical University Aviation General Hospitalh1412000Chinese People’s Armed Police Force Beijing Corps Hospitalh1515000Beijing Huaxin Hospitalh167600Coal General Hospitalh175000The Third Affiliated Hospital of Beijing University of Chinese Medicineh185200Beijing Chaoyang District Second Hospitalh194500Shuangqiao Hospital, Chaoyang District, Beijingh205500Beijing Coking Chemical Factory Hospitalh213000Chaoyang District Chinese Medicine Hospital, Beijingh226500Beijing Chaoyang District Maternal and Child Health Hospitalh234500Beijing Tibetan Hospitalh241600Beijing Wuzhou Women and Children Hospitalh255000Beijing Dangyangliu Hospital East Campush264500Table 4
Capacity of ESs.
ESsMAjCapacityMAkCapacityYuan Dadu Archaeological Site Parkj190000k1100000Chaoyang Parkj268000k2182000Sun Palace Parkj320000k390000Olympic Forest Park South Parkj47000k413000Anzhen Yongxi Parkj52600k53000Taiba River Green Beltj630000k650000Xinglong Parkj732800k740000Red Scarf parkj816100k818000Beixiaohe Parkj915700k920000Beijing Liyuanj1050000k1055400Wanghe Parkj1139200k1140000Cuicheng Parkj121950k124000Hongbo Country Parkj13100000k13140000West Dawang Road Community Parkj147200k147000Cuicheng Parkj121950k124000Hongbo Country Parkj13100000k13140000West Dawang Road Community Parkj147200k147000Lishuiqiao Parkj155000k157000Dongba Country Parkj1655000k1680000White Deer Parkj1725000k1740000Jintian Parkj1860000k1890000Dongfeng Parkj1910000k1930000Haitang Parkj2024000k2050000Laojuntang Parkj2148000k2160000Duzhong Parkj2255000k2280000Changying Parkj2377000k2390000Chaolai Forest Parkj2444000k2470000Guta Parkj2555000k2570000Sun River Country Parkj2614000k2630000Dawangjingj2735000k2740000Red Army Parkj3182500k3190000Because of the lack of EDC data in Chaoyang District, Beijing, combining the total demand of resources of all refugees in each node, we assume here that 9 EDCs were selected from 42 candidate EDCs based on previous EDC selection research studies [4, 8, 31–33]. The capacity of each EDC can be derived from the total number of EDCs, total quantity of required resources, and the ratio α.


### 5.1.5. Distance

We obtain the road networks from the map of Beijing published by the China Cartographic Publishing House in 2016. We calculate the actual distances from AAs to A-ESs, B-ESs, and EMCs and from EDCs to AAs, A-ESs, B-ESs, and EMCs using version 10.1 of ArcGIS.


### 5.1.6. Cost

The total cost involves the construction cost of candidate A-ESs, B-ESs (existing shelters have no construction cost), and EDCs and the fixed and transportation cost of per unit resourcer. The construction cost of each site is determined by its capacity, which is estimated by the government report and previous research studies [5, 18, 19]. The fixed cost of resources is also estimated by previous research studies, and the transportation cost of resources is estimated by the Beijing logistic price. All these parameters are shown in Table 5.Table 5
Cost parameters.
ParameterEstimationConstruction costCj=40; Ck=15; Cc=15Fixed costCRr=10r=1; CRr=20r=2Transportation costCTr=1


## 5.1.1. Facility Data

According to the database of the Beijing Earthquake Agency, the Chaoyang district has built fifteen earthquake disaster shelters. However, these shelters may not accommodate so many refugees or cause more evacuation distances. Thus, based on the Beijing Central City Earthquake and Emergency Refuge Site Planning Outline and the Chaoyang Government Network, another sixteen candidate shelters were selected. The EMCs need professional medical equipment and nursing staffs. Thus, in our research, we select 26 candidate EMCs from the hospitals located in the Chaoyang district, as shown in Figure7. Since there are no EDCs in the Chaoyang District, Beijing, we assume that candidate EDCs are located at street/district offices. So, there are 42 candidate EDCs.Figure 7
Distribution map of the Chaoyang district.


## 5.1.2. Population

Chaoyang district governs 24 streets and 19 districts, and the population size of each street/districtEa is extracted from the Chaoyang District Beijing Statistical Yearbook 2017. The number of refugees in each AA of each time interval Eat is roughly estimated [17], and Figure 8 presents the ratio of Eat to Ea. We assume that 10% of the population is injured EIat, 20% of the population is assumed to be A-priority refugees EAat, and 30% is assumed to be B-priority refugees EBat.Figure 8
Number of refugees in AAs over time.


## 5.1.3. Demand

Quantity of required resources of different refugees is estimated by previous research [17]. Each refugee requires two types of relief resources: commodities and medical supplies. Commodities per unit are herein 1 kg, and the medical supplies per unit are 0.5 kg. Different refugees have different quantity of demand for different resources, and the specific resource requirements are presented in Table 2.Table 2
Quantity of required resources of hierarchical refugees.
δjrδkrδhrδhrr=111.50.52r=21.5120.5


## 5.1.4. Capacity

We assume that each ES is divided into two parts: A-ES and B-ES, where the capacity of A-ESs is smaller than B-ES. The capacities of A-ES, B-ES, and EMCs are presented in Tables3 and 4.Table 3
Capacity of EMCs.
EMCsMAhCapacityBeijing Obstetrics and Gynecology Hospital, Capital Medical Universityh16600Beijing Obstetrics and Gynecology Hospital, Capital Medical Universityh28000Beijing Chaoyang Hospital, Capital Medical Universityh319000Wangjing Hospital, China Academy of Chinese Medical Sciencesh413000China-Japan Friendship Hospital of the Ministry of Healthh515000People’s Liberation Army No. 306 Hospitalh66000Beijing Anzhen Hospital, Capital Medical Universityh715000Beijing First Hospital of Integrated Traditional Chinese and Western Medicineh816000Beijing Ditan Hospital, Capital Medical Universityh99000Chinese Academy of Medical Sciences Cancer Hospitalh1016000Children’s Hospital of the Capital Institute of Pediatricsh116000Beijing Dangyangliu Hospitalh125600Civil Aviation General Hospitalh136500China Medical University Aviation General Hospitalh1412000Chinese People’s Armed Police Force Beijing Corps Hospitalh1515000Beijing Huaxin Hospitalh167600Coal General Hospitalh175000The Third Affiliated Hospital of Beijing University of Chinese Medicineh185200Beijing Chaoyang District Second Hospitalh194500Shuangqiao Hospital, Chaoyang District, Beijingh205500Beijing Coking Chemical Factory Hospitalh213000Chaoyang District Chinese Medicine Hospital, Beijingh226500Beijing Chaoyang District Maternal and Child Health Hospitalh234500Beijing Tibetan Hospitalh241600Beijing Wuzhou Women and Children Hospitalh255000Beijing Dangyangliu Hospital East Campush264500Table 4
Capacity of ESs.
ESsMAjCapacityMAkCapacityYuan Dadu Archaeological Site Parkj190000k1100000Chaoyang Parkj268000k2182000Sun Palace Parkj320000k390000Olympic Forest Park South Parkj47000k413000Anzhen Yongxi Parkj52600k53000Taiba River Green Beltj630000k650000Xinglong Parkj732800k740000Red Scarf parkj816100k818000Beixiaohe Parkj915700k920000Beijing Liyuanj1050000k1055400Wanghe Parkj1139200k1140000Cuicheng Parkj121950k124000Hongbo Country Parkj13100000k13140000West Dawang Road Community Parkj147200k147000Cuicheng Parkj121950k124000Hongbo Country Parkj13100000k13140000West Dawang Road Community Parkj147200k147000Lishuiqiao Parkj155000k157000Dongba Country Parkj1655000k1680000White Deer Parkj1725000k1740000Jintian Parkj1860000k1890000Dongfeng Parkj1910000k1930000Haitang Parkj2024000k2050000Laojuntang Parkj2148000k2160000Duzhong Parkj2255000k2280000Changying Parkj2377000k2390000Chaolai Forest Parkj2444000k2470000Guta Parkj2555000k2570000Sun River Country Parkj2614000k2630000Dawangjingj2735000k2740000Red Army Parkj3182500k3190000Because of the lack of EDC data in Chaoyang District, Beijing, combining the total demand of resources of all refugees in each node, we assume here that 9 EDCs were selected from 42 candidate EDCs based on previous EDC selection research studies [4, 8, 31–33]. The capacity of each EDC can be derived from the total number of EDCs, total quantity of required resources, and the ratio α.


## 5.1.5. Distance

We obtain the road networks from the map of Beijing published by the China Cartographic Publishing House in 2016. We calculate the actual distances from AAs to A-ESs, B-ESs, and EMCs and from EDCs to AAs, A-ESs, B-ESs, and EMCs using version 10.1 of ArcGIS.


## 5.1.6. Cost

The total cost involves the construction cost of candidate A-ESs, B-ESs (existing shelters have no construction cost), and EDCs and the fixed and transportation cost of per unit resourcer. The construction cost of each site is determined by its capacity, which is estimated by the government report and previous research studies [5, 18, 19]. The fixed cost of resources is also estimated by previous research studies, and the transportation cost of resources is estimated by the Beijing logistic price. All these parameters are shown in Table 5.Table 5
Cost parameters.
ParameterEstimationConstruction costCj=40; Ck=15; Cc=15Fixed costCRr=10r=1; CRr=20r=2Transportation costCTr=1


## 5.2. Model Testing and Results

In this section, the proposed model and algorithms are tested with the abovementioned parameter settings. As the priority of refugees in different places is EMCs> A-ESs > B-ESs > AAs and ∑∗ω∗=1 (∗∈EMCs, A-ESs, B-ESs, and AAs), we set the priorities of different places ωh, ωj, ωk, and ωa as following. We consider that the quantity of resources in all EDCs is insufficient and higher priority places have a higher minimum resource satisfaction degree. Different satisfaction degree settings are based on the ratio of total resources in DCs to total required resources. α=0.9 means that the total resources of the emergency distribution center account for 90% of the total resources needed for the entire rescue cycle. μ is a measure of the fairness loss function, and μ<0. When the amount of unsatisfied resources increases, the fairness loss function grows like an exponential function as μ decreases and grows like a linear function as μ increases. For parameters of SA, ω1, ω2, and ω3 are the weights of the three objective functions in the model, and also the weights of the three subobjective functions in the fitness function of the outer simulated annealing algorithm. ts is the start temperature, te is the end temperature, ρ is the decay rate of temperature, ε is the received tolerance for inferior solutions, key is the optimal solution unchanged times, keymax is the maximum optimal solution unchanged times, and L is the iterate times in each temperature. For parameters of PSO, K is the maximum iteration times, N is the particle number, c1 and c2 are the learning factors, c1 indicates the possibility of the particle updating according to its local optimal position, c2 indicates the possibility of the particle updating according to the global optimal position, w is the inertia factor, which reflects the particle’s ability to inherit the previous speed, and ω2′ and ω3′ are the weight setting of the objective functions f2′ and f3. For parameters of GA, G represents the number of groups, pc is the crossover rate, and pm is the mutation rate. According to relevant research [20–22, 25, 26, 29, 30], parameters of the model and algorithms are also shown in Table 6.Table 6
Parameter setting.
ParameterEstimationParameters of the modelωh=0.32;ωj=0.28;ωk=0.22;ωa=0.1;α=0.9;μ=−0.0001Parameters of SAts=100;te=1;ρ=0.9;L=10;ε=0.0002;ω1=1/3;ω2=1/3;ω3=1/3;key=0;keymax=50Parameters of PSON=20;K=100;c1=0.9;c2=0.9;w=0.8;ω2′=0.5;ω3′=0.5Parameters of GAG=50;pc=0.8;pm=0.2By setting the weight of each objective function equal to 1, with the objective of minimizing/maximizingfs of SA, we calculate the minimum/maximum value of each objective of our two test cases (as shown in Table 1).We firstly test the case with these existing shelters as case 1, in which the construction cost of A-ESs and B-ESs set to be zero. Then, these candidate shelters are added in the second test case as case 2, using the parameters shown in Table6.Through 50 experiments, the average value of objective functions of the two test cases is shown in Table7.Table 7
Results of the two test cases.
Objective functionCase 1Case 2Dividend (%)f1(evacuation distance)8.92E + 065.27E + 06−40.92f2(cost)9.05E + 089.40E + 083.85f3(unmet resources)2.13E + 061.97E + 06−7.51The dividend is calculated by the equation=fcase1/fcase2−1, as shown in Table 7, and the results reveal that the total evacuation distance in case 2 decreases by 40.92%, while the total cost increases by 3.85% with selected candidate shelters compared to case 1. And, we see that the total unmet resources in case 2 also decreases. Thus, it is obvious that selecting candidate shelters can reduce the total evacuation distance to a large extent and also reduce the total unmet resources while only adding a small amount of cost.Taking the corresponding refugee evacuation in which the objective function value is the closest to the average objective function value as an example, the allocation results of refugees from AAs to A-ESs and B-ESs at timet=1 of the two cases are shown in Figure 9.Figure 9
Comparison of the two cases of refugee evacuation at timet=1. (a) Refugee evacuation with only the existing shelters. (b) Refugee evacuation with candidate shelters.
(a)(b)One can see from Figure9(a) that many AAs have no nearby shelter, leading to the relatively long evacuation distance of refugees in these AAs. However, as is shown in Figure 9(b), six A-ESs and six B-ESs were selected at time t = 1 (one more A-ES and B-ES were selected at time t=2) and optimized the long-distance allocation result. Through the comparison of Figures 9(a) and 9(b), it is obvious that selecting candidate shelters can greatly reduce the evacuation distance of some AAs.


## 5.3. Scenario Analyses

In this section, five scenarios with added candidate shelters are designed to compare the characteristics of our proposed model and change in disaster situation of refugees. As can be seen from the results of Section5.2, selecting candidate shelters is a more efficient planning method. Therefore, all the scenarios in this section are based on case 2.Scenario 1.
The first scenario is to compare the efficiency of SA-PSO and ISA-PSO. Parameters of the GA are shown in Table8. We analyze these two algorithms by comparing the results of case 2 under SA-PSO and ISA-PSO. The operational iteration graph of the two algorithms is shown in Figure 10, and through 50 experiments, the average value of three objective functions is presented in Table 8.
It can be seen from Figure10 that the iterative result of ISA-PSO is always better than that of SA-PSO. Table 9 presents the objective values calculated by ISA-PSO. The dividend is calculated by the equation =fscenario1/fcase2−1, and it can be seen that compared with SA-PSO, the objective function values by ISA-PSO can achieve Pareto optimization. But, because of the characteristics of population search of ISA-PSO, the time consumption of ISA-PSO is correspondingly longer.
The algorithm used in the following calculations will use ISA-PSO. And, the following comparison of objective values will be based on the results of ISA-PSO.Table 8
Results of Scenario 1.
Objective functionScenario 1Dividend (%)f1(evacuation distance)5.03E + 06−4.55f2(cost)9.34E + 08−0.62f3(unmet resources)1.88E + 06−4.57Figure 10
Iteration graph of SA-PSO and ISA-PSO.Table 9
Results of Scenario 2.
Objective functionS2-1S2-2S2-3f1(evacuation distance)4.46E + 066.08E + 064.43E + 06f2(cost)9.70E + 089.28E + 089.79E + 08f3(unmet resources)1.87E + 061.97E + 061.87E + 06Scenario 2.
The second scenario is to analyze the impact of the objective weight. Thus, we set two of the objective weights to 50% of the original and the other target weight to twice the original: S2-1:ω1=1/6,ω2=1/6,andω3=2/3; S2-2: ω1=1/6,ω2=2/3,andω3=1/6; S1-3: ω1=2/3,ω2=1/6,andω3=1/6. Through 50 experiments, the average value of three objective functions is presented in Table 9.
It can be seen from Table9 that f1 is greatly affected by the objective weight. When the weight increases, the objective value is optimized, while the objective value of f2 and f3 is affected by the weight to a small extent. In addition, it can be seen that f1 and f3 are negatively correlated with f2, and f2 increases as f1 and f3 decrease. The corresponding dividend comparison of three objective values of Scenario 2 is shown in Figure 11.
The dividend is calculated by the equation=fS2−∗/fscenario1−1∗∈1,2,3, and it can be inferred that f1 is the most sensitive to the objective weight. Therefore, provided that the budget is sufficient, increasing the objective weight of f1 appropriately can better optimize the result.Figure 11
Dividend objective values of Scenario 1.Scenario 3.
The third scenario is to analyze the effect of parameterα on three objective functions. The purpose of considering this scenario is to enable managers to improve the efficiency of resource allocation under uncertain or insufficient conditions. In our research, we assume the quantity of relief resources is insufficient, and here we change the ratio α as S3-1: α=0.75, S3-2: α=0.85, and S3-3: α=0.95, respectively. Through 50 experiments, the average value of three objective functions of Scenario 3 is presented in Table 10.
Table10 presents that the change in α has little effect on f1 and f2. It implies that resources and their transportation costs account for only a small part of the total cost (f2). And, obviously, as α increases, f3 decreases. Figure 12 illustrates the corresponding dividend comparison of Scenario 2. The dividend is calculated by the equation =fS3−∗/fscenario1−1∗∈1,2,3. It is observed that as α increases, f3 decreases rapidly and there were almost no change in f1 and f2. This scenario illustrates that adjusting the parameter α significantly changes objective values of f3 in a linear relationship, but does not affect the other objects. This implies that increasing the quantity of resources is an effective method that improves the target value of f3. It can be seen that preparing sufficient resources in advance can play a vital role in the disaster relief process at a small cost.Table 10
Results of Scenario 3.
Objective functionS2-1S2-2S2-3f1(evacuation distance)5.13E + 064.98E + 064.74E + 06f2(cost)9.23E + 089.41E + 089.65E + 08f3(unmet resources)7.89E + 063.35E + 066.81E + 05Figure 12
Dividend objective values of Scenario 3.Scenario 4.
This scenario is to analyze the impact of population in each priority. Here, we increase the number of A-priority and B-priority refugees to be 30% and 40% of the population in S4-1 and 40% and 50% of the population in S4-2.
In this scenario, the total capacity of the existing shelters cannot meet the refugee demand. Thus, the selecting of candidate shelters is essential when the refugee demand increases.
Through 50 experiments, the average value of three objective functions of Scenario 4 is presented in Table11. The dividend is calculated by the equation =fs4−∗/fscenario1−1∗∈1,2. It is obvious that when the number of A- and B-priority refugees increases, the total evacuation distance increases evidently.
Taking the corresponding refugee evacuation in which the objective function value is the closest to the average objective function value as an example, the allocation results of refugees from AAs to A-ESs and B-ESs att = 1 of Scenario 4 are shown in Figure 13. Combining the results in Table 9, all these imply that the construction cost has only very small impact on the total cost f2. Besides, when the number of refugees who need to evacuate increases, the total evacuation distance will also increase even if more candidate shelters are selected. Furthermore, when the evacuation demand of A- and B-priority refugees increases, the number of refugees in A-ESs and B-ESs increases as well. This result demonstrates that refugees in A-ESs and B-ESs are more sensitive to unmet resources than those in AAs.Table 11
Results of Scenario 4.
Objective functionS4-1Dividend (%)S4-2Dividend (%)f1(evacuation distance)7.69E + 0655.049.90E + 0699.50f2(cost)9.51E + 082.429.82E + 085.86f3(unmet resources)2.13E + 0611.262.61E + 0636.65Figure 13
Comparison of the refugee evacuation results of Scenario 4 at timet=1. (a) Refugee evacuation plan of S4-1. (b) Refugee evacuation plan of S4-2.
(a)(b)Scenario 5.
In this scenario, we set the population of each period to be stable and compare the results with time-varying demand and stable demand. Figure14 displays the affected population in each time interval with a stable refuge demand. Through 50 calculations, the average value of three objective functions of Scenario 5 is presented in Table 12. The dividend is calculated by the equation =fscenario5/fscenario1−1. The results in Table 12 show that the total evacuation distance increment is 14.8% higher, while the total cost is only 0.54% less than the result under time-varying demand. The reason of the evacuation distance dividend is that the remaining capacity of A-ESs and B-ESs becomes smaller as time goes by, but the number of A- and B-priority refugees has not become smaller. Therefore, refugees in the later period will have to evacuate to farther A-ESs and B-ESs.
Similarly, most of the resources in EDCs were allocated to the AAs in the early period. The remaining resources in EDCs are getting less and less as time goes by. However, the demand for resources of A-ESs and B-ESs is increasing much faster than the time-varying demand model. Therefore, the resource requirements of A-ESs and B-ESs may not be met in the later period. Recall the implication of Scenario 4, and this might be the reason that the objective value of unmet resources is also 41.88% higher than the results of the time-varying demand. Thus, it is obvious that the time-varying demand estimation is a more efficient evacuation and allocation planning model because it reduces evacuation distance and unmet resources. Furthermore, it fits the actual refuge demand better.Figure 14
Unvarying demand of refugees in AAs over time.Table 12
Results of Scenario 5.
Objective functionS5-1Dividend (%)f1(evacuation distance)5.69E + 0614.80f2(cost)9.23E + 08−0.54f3(unmet resources)2.71E + 0641.88


## 6. Conclusions

This paper studies decision-making on the dynamic emergency shelter and resource allocation for disaster management, considering the hierarchical requirement of refugees and allocation of insufficient relief resources. The planning of evacuating hierarchical refugees to hierarchical shelters and EMCs along with the corresponding resource allocation to these refugees is a difficult problem with multicriteria restrictions, especially when the quantity of relief resources is insufficient. Thus, a multiobjective location-allocation model of minimizing the total hierarchical evacuation distance, total cost, and total hierarchical unmet resources is proposed to characterize the problem. It considers both traditional objectives and fairness perceived by refugees. The SA algorithm was applied to site selection in the outer layer, and the PSO algorithm in the inner layer is to solve the corresponding resource allocation. And, an improved SA-PSO is proposed in order to expand the local exploration scope of SA and enhance its population diversity.A case study of Chaoyang District, Beijing, is conducted to verify the capability of the proposed model and algorithms. The computational experiment results indicate that selecting candidate shelters for refugees is obviously more worthwhile and more effective than evacuating refugees to existing shelters. The results of this study can be considered as a shelter selection proposal for Chaoyang District, Beijing, emergency facility planning. Further, some scenario analyses are considered for further investigation on emergency refugee evacuation and resource allocation decisions. The scenario analyses not only provide useful suggestions on site selection of candidate emergency facilities and corresponding resource allocation, but also perceive a more cost-effective and realistic estimation on the refugee demand.For future research, the real-time road network like congestion and destruction should be considered in the refugee evacuation planning. Thus, the evacuation planning and resource allocation planning would be different over time, with the change in the road distance. In addition, integration with physicians and nurse allocation and multitype vehicle allocation could provide a more comprehensive disaster reduction system.



---

*Source: 8395714-2020-08-19.xml*