Multi-objective optimization approach of shelter location with maximum equity: an empirical study in Xin Jiekou district of Nanjing, China

Abstract

Based on the inadequacies and neglect of the equity of refuge resources, refuge demands and the evacuation allocation of traditional methodologies, this study put forwards the multi-objectives layout optimization model of shelters which firstly realizes the maximum equity of shelter location. Our approach has the objectives of maximizing equity, minimizing overall egress time and minimizing the quantity of new shelters. The high-precision population is established through mobile signaling data, while the optimization model adopts a circular circulatory allocation rule derived from a gravity model. The shorter the evacuation time, the larger the shelter capacity and thus more refugees are allocated to the shelter. The evacuation time is determined by the application programming interface (API) of the Baidu Map open platform with Python, which exhibits the authentic evacuation paths and real-time traffic conditions. This study designs a three-stage algorithm ‘genetic algorithm-exhaustive method-evaluation’. The first process of algorithm calculates the minimum quantity of new shelters; the second process selects the feasible layout schemes and determines Pareto optimum solutions; and the third stage evaluates the Pareto optimum solution based on the shelter construction cost and the accessibility from shelters to emergency supply storage points to determine the best location scheme. This study regards Xin Jiekou district in Nanjing as a case area to demonstrate reliability and availability of the proposed methodology.

Keywords:

• Shelter
• multi-objective optimization
• maximum equity
• three-stage optimization algorithm
• evaluation of Pareto optimum solution

Introduction

A shelter is a refuge building/area equipped with emergency infrastructure, support equipment, auxiliary facilities and supplies for post-disaster life support (China MOHURD 2015). It is an important functional element when facing various natural and man-induced calamities. The construction of shelters is a key component of disaster management, and includes determining the location and assigning refugees to shelters. A suitable shelter layout is key to effectively improve the resilience of cities.

Previous optimization models are on the basis of P-center, P-median, and covering models (Gama et al. 2013). P-median and its derivative models are dedicated to determine shelter location with the quantity of ‘P’ to evacuate all refugees as quickly as possible (Tao et al. 2018). The objective of the maximum covering model (Church and Davis 2005, Liu et al. 2023) is to provide refuge service for as many refugees as possible with a limited amount of shelters. The set covering model (Church and Davis 2005) aims to house all refugees with a minimum number of shelters. Furthermore, the objective of the P-center model (Kilci et al. 2015) is to minimize the maximum evacuation cost, resulting in a local (not global) optimum. These inchoate models have just a single objective and ignore constraint conditions.

Due to the complex layout optimization of shelters, scholars are no longer satisfied with establishing a single-objective optimization model, and numerous have proposed multi-objective optimization models through integrating the construction cost of shelters (Panchalee et al. 2021), overall evacuation time (Yu et al. 2018), evacuation distance (Hu et al. 2012), risk on evacuation routes (Alçada-Almeida et al. 2009; Doerner et al. 2009; Zhang et al. 2015) and other optimization objectives (Barzinpour and Esmaeili 2014). For example, Wu and Weng (2011) established an optimization model which minimizes refuge evacuation distance and quantity of shelters by limiting the utilization efficiency of shelters, and compared the differences between multiple schemes via empirical analysis. Chen (2019) and Xiang et al. (2020) considered the influence of uncertainty of road conditions and traffic capacity on shelter layout, which minimizes the refugee evacuation distance and optimizes the conditions of candidate evacuation routes. In general, the methods for solving multi-objective optimization models can be divided into three approaches: (1) transforming multifarious objective function to various single-objective functions through giving different weights to different objectives (Zhang et al. 2020; (2) transforming multi-objective problems to several single-objective problems at different stages (Kongsomsaksakul et al. 2005; Li et al. 2012); and (3) determining Pareto optimum solution via a heuristic algorithm (Doerner et al. 2009; Li et al. 2017).

The development from a single-objective to multi-objective shelter location model has brought about numerous restrictions on the evacuation distance and time (Ding et al. 2013; Ma et al. 2018). The models derived in previous literature consider various shelter hierarchies (Chen et al. 2013), the maximum shelter capacity (Saadatseresht et al. 2009), the behavioral rationality of refugees (Ng et al. 2010; Chen et al. 2019), disaster risk (Fan et al. 2014) and other complicated factors. For instance, Li et al. (2012) developed a hierarchical model to allocate shelters under hurricane hazard based on hurricane scenario and occurrence time. Zhao et al. (2015) considered the influence of diverse seismic scenarios on the quantity of refugees.

The increment on the quantity of constraints and objectives coverts shelter location-allocation problems as non-deterministic polynomial (NP) hard problem (AlShahrani 2019). Numerous heuristic algorithms (Betül et al. 2021; Yildiz and Erdas 2021) have been conducted to solve NP hard problems. Natee et al. (2021) compared the performance of fourteen new and proposed multi-objective metaheuristics. Doerner et al. (2009) and Hu et al. (2014) determined shelter layout scheme with multi-objectives under the events of a hurricane and earthquake using NSGA-II. Kongsomsaksakul et al. (2005) determined optimal hierarchical shelter layout scheme under flood hazard via a genetic algorithm. As well as genetic algorithms, the PSO (Marinakis and Marinaki 2008), simulated annealing algorithm and ant colony algorithm (Arnaout 2013) have also been developed to solve location optimization of shelters. Wu et al. (2018) and Ng et al. (2010) applied the simulated annealing algorithm to the shelter location model to simulate the refugee process from unordered searching to orderly resettlement. Arnaout (2013) minimized the construction costs and adopted the ant colony algorithm to determine layout of shelters and allocation schemes. Zhao et al. (2015) and Hu et al. (2012) combined the simulated annealing and PSO to optimize the layout of seismic shelters. Shankar and Baviskar (2018) proposed evolutionary algorithms with innovative strategies to solve the multi-objective optimization problems which afford faster convergence and better diversity of the Pareto front.

Despite the progress made by the current multi-objective layout optimization model of shelters, they are associated with several limitations. For example, the difference in the population distribution during the day and night is not considered. Current research calculates the refugee demands using permanent population values, which do not simulate the dynamic variations of the population under time-space dimension. In addition, when converting multi-objectives to various single objective, the weight of each objective is difficult to determine. For the conversion to several single-objective problems at different stages, no solution exists when the model has more than two objective functions. If Pareto optimum solutions are calculated via NSGA-II algorithm (Tang et al. 2020), the best location scheme is not provided for decision-makers (Hu et al. 2014). Moreover, the layout of shelters only depends on the quantity of refugees and the spatial relationship between demand points and shelters. Disaster rescue is a complete system, and includes medical rescue, emergency evacuation and emergency supplies. Following a disaster, it is of necessity to distribute emergency supplies to shelters in a timely manner. Present models do not consider the spatial relationship between shelters and emergency supply storage points, and neither to the fully account for different disaster prevention systems. Furthermore, the shelter location model typically assumes that citizens from one demand plot are allocated to a single shelter, which leads to low resource utilization and evacuation efficiencies. Lastly, traditional location models focus on the efficiency of shelter locations, and the optimization objectives include maximizing the number of refugees (maximum coverage model), minimizing the quantity of shelters (set covering model), and minimizing the evacuation cost from the demand points to shelters (P-median model). However, due to the difficulties in modelling and measuring the equity of shelter locations, previous studies generally ignore this variable. Tao et al. (2014) and Li et al. (2017) proposed maximum equity model on residential care facility locations, but the equity on shelter location is not mentioned in previous literature.

Based on the shortcomings of traditional multi-objective shelter location and evacuation allocation models, this study takes equity, efficiency, and economy as the layout optimization principles, and proposes a multi-objective shelter layout optimization method via circular evacuation allocation. The objective functions are minimum total evacuation time, minimum quantity of shelters and maximum equity. In addition, the model is constrained by service range and capacity of shelters. Equity is measured through the difference in accessibility from demand plots to shelters. The difference in accessibility is described through the square deviation of the accessibility from each demand point to shelters.

The proposed approach is summarized as follows. First, the dynamically changing refuge demands during the day and night are determined through mobile signaling data. Second, the optimization model assigns refugees from demand plots to shelters through the gravity model. The shorter the evacuation time, the larger the capacity of shelter, and the greater the amount of refugees allocated. If the number of people receiving in shelter reaches its maximum capacity, the evacuation will stop. Refugees who have not entered to shelters will experience next cycle until all refugees are evacuated. The evacuation time is determined by the application programming interface (API) of the Baidu Map open platform using Python, which can reveal authentic evacuation routes and real-time traffic conditions. Considering the uncertainty in the occurrence of hazards and the evacuation time, the request time for the API is set for both the night and day, which is consistent with the population distribution. Third, a three-stage algorithm is designed to solve the multi-objective optimization layout of shelter. In the first phase, the minimum quantity of new shelters are obtained through genetic algorithm. Feasible schemes are selected from the first stage that satisfy refuge demands under the capacity and service distance constraints in the second phase. The quantity, total evacuation time and Z value (variance between accessibility evaluations from demand points to shelters) of all feasible location schemes are counted to derive Pareto optimum solutions. Based on accessibility evaluation outcomes from shelters to the emergency supply storage points and the number of new shelters, the optimal layout scheme is determined from the Pareto optimum solutions. Furthermore, Xinjiekou District in Nanjing, China, is taken as a case area to evaluate, and verify the proposed methodology.

Compared with traditional shelter location model which concentrates on efficiency and construction cost, this study also focuses on equity of shelter location. The authentic evacuation paths and real-time traffic conditions are considered in OD matrix calculation. An innovative three-stage algorithm ‘genetic algorithm-exhaustive method-evaluation’ and evaluation model for the Pareto optimum solution are designed.

The next section demonstrates the proposed research methodology, including multi-objective layout optimization model of shelters, cyclic assignment model, evaluation model for the Pareto optimum solution and three-stage algorithm. Section 3 describes the empirical area, basic datasets and results of shelter location. Section 4 summarizes the merits and demerits of the proposed methodology. Section 5 illustrates the conclusion with suggestions for future research.

Research methodology

The optimization method includes four key steps: (1) Based on remote sensing images, OpenStreetMap (Haklay and Weber 2008) and field research, the spatial distribution of demand points, shelters, emergency supply storage points, and road grades are determined. The refuge demands of each demand point in the study area are calculated through mobile phone signaling data during the day and night. (2) A multiple objective optimization model for shelter layout is proposed, with the objectives on maximization of equity, minimization of the quantity of shelters and total evacuation time. (3) A three-stage algorithm is constructed to solve the multi-objective optimization model for shelter layout. First, the minimum quantity N of new shelters satisfying refuge demands are determined via genetic algorithm. Second, feasible schemes are selected from ${\mathit{C}}_M^N,$ ${\mathit{C}}_M^{N+1},$ ${\mathit{C}}_M^{N+2}$…, ${\mathit{C}}_M^{M-1},$ ${\mathit{C}}_M^M$ schemes (M represents the quantity of candidate shelters). ${\mathit{C}}_M^N$ is the symbol for permutation and combination, which means selecting N shelters from M candidate shelters. These feasible schemes should meet the refuge demands, maximum capacity and refuge distance limitation. The quantity, total evacuation time and Z values of all feasible location schemes are counted to determine Pareto optimum solutions. Based on accessibility evaluation outcomes from shelters to emergency supply storage points and the number of new shelters, the exclusive optimal layout scheme is chose from Pareto optimum solutions. Figure 1 presents the optimization process.

Figure 1. Multi-objective optimization of shelter location based on circular evacuation allocation.

Basic hypothesis

The evacuation roads are still safe and reliable after the disaster, and blocked roads are not considered.

Multi-objective layout optimization model of shelters

The classical shelter layout and evacuation allocation model mainly focuses on the model efficiency, including maximizing the quantity of the covered population, and minimizing the construction and evacuation costs from the demand points to shelters. Due to difficulties in modelling the equity of shelter layouts, it is not considered in previous research on the optimization of shelters. This study takes maximum equity, high efficiency, and low construction costs as the layout optimization principles, and evacuation time of all refugees, the quantity of newly built shelters and variance of accessibility from demand plots to shelters as the objective functions to propose a three-objective shelter layout optimization model. In order to measure the shelter location equity, the spatial accessibility evaluation model is first established. Accessibility represents the difficulty degree of travelling from the demand point to the shelter, or overcoming spatial obstacles to reach the destination (Panchalee et al. 2021). The shelter location equity is measured by minimizing the differentiation between the accessibility from each demand point to the shelters. The difference in accessibility is reflected in the square deviation of the accessibility of each demand point (Doerner et al. 2009; Wang and Tang 2013). According to the improved two-step floating catchment area method, the objective function of equity maximization is described in Equations (1–4). (1) $${\mathit{A}}_i=0.5\times \sum _{\mathit{j}=1}^{\mathit{n}}\frac{{\mathit{C}}_j\cdot \mathit{f}({\mathit{t}}_{ij,\mathit{daytime}})\cdot {\mathit{y}}_j}{\sum _{\mathit{k}=1}^{\mathit{m}}{\mathit{D}}_{k,\mathit{daytime}}\cdot \mathit{f}({\mathit{t}}_{kj,\mathit{daytime}})}+0.5\times \sum _{\mathit{j}=1}^{\mathit{n}}\frac{{\mathit{C}}_j\cdot \mathit{f}({\mathit{t}}_{ij,\mathit{nighttime}})\cdot {\mathit{y}}_j}{\sum _{\mathit{k}=1}^{\mathit{m}}{\mathit{D}}_{k,\mathit{nighttime}}\cdot \mathit{f}({\mathit{t}}_{kj,\mathit{nightime}})},$$(1) (2) $$\mathit{f}({\mathit{t}}_{ij})=\left\{\begin{array}{c}\frac{{\mathit{e}}^{-1/2{({\mathit{t}}_{ij}/{\mathit{t}}_0)}^2}-{\mathit{e}}^{-1/2}}{1-{\mathit{e}}^{-1/2}}{\mathit{t}}_{ij}\le {\mathit{t}}_0\\ \\ \\ 0{\mathit{t}}_{ij}>{\mathit{t}}_0\mathit{ }\end{array}\right\}\mathit{f}({\mathit{t}}_{kj})=\left\{\begin{array}{c}\frac{{\mathit{e}}^{-1/2{({\mathit{t}}_{kj}/{\mathit{t}}_0)}^2}-{\mathit{e}}^{-1/2}}{1-{\mathit{e}}^{-1/2}}{\mathit{t}}_{kj}\le {\mathit{t}}_0\\ \\ \\ 0{\mathit{t}}_{kj}>{\mathit{t}}_0\mathit{ }\end{array}\right\},$$(2) (3) $$\alpha =\frac{\mathit{C}}{\mathit{D}},$$(3) (4) $$Z=\sum _{i=1}^m{\left(A_{\mathit{i}}-\alpha \right)}^2,$$(4) where i and k are serial numbers of all demand plots; m represents the quantity of demand plots; n is the quantity of all shelters, including current and candidate shelters; t₀ is the time limitation. If we take the target object as the resident emergency aggregate shelter, the maximum service distance limitation as 1500 m (China MOHURD 2015), and the average evacuation speed as 1.27 m/s (Wang 2012), then the time limitation is 1181s. t_ij,daytime and t_ij,nighttime are the evacuation time from plot i to shelter j during day and night; t_kj,daytime and t_kj,nighttime are the evacuation time from plot k to shelter j during day and night; A_i is the measurement of refuge resource accessibility at demand plot i; j represents sequence number of candidate shelters, if shelter j is selected, y_j = 1 and if shelter j is not selected, y_j = 0; C_j denotes the maximal capacity of shelter j; C denotes the total supply scale of refuge resources, which is the sum of the capacities of all selected shelters; D_{k, daytime} and D_{k, nighttime} are the number of refugees during the day and night at demand point k; D is the total demand scales, which is half of the total quantity of refugees during the night and day; α is the proportion of the supply scale of refuge resources to demand scale; f(t_ij) and f(t_kj) represent the attenuation function, where the Gaussian attenuation function is used to evaluate the accessibility of shelters (Dan et al. 2018); and Z is the equity of refugee access to the refuge resources under a group of shelter location schemes, that is, the differentiation in accessibility from each demand plot to shelter.

Equation (1)(1) $${\mathit{A}}_i=0.5\times \sum _{\mathit{j}=1}^{\mathit{n}}\frac{{\mathit{C}}_j\cdot \mathit{f}({\mathit{t}}_{ij,\mathit{daytime}})\cdot {\mathit{y}}_j}{\sum _{\mathit{k}=1}^{\mathit{m}}{\mathit{D}}_{k,\mathit{daytime}}\cdot \mathit{f}({\mathit{t}}_{kj,\mathit{daytime}})}+0.5\times \sum _{\mathit{j}=1}^{\mathit{n}}\frac{{\mathit{C}}_j\cdot \mathit{f}({\mathit{t}}_{ij,\mathit{nighttime}})\cdot {\mathit{y}}_j}{\sum _{\mathit{k}=1}^{\mathit{m}}{\mathit{D}}_{k,\mathit{nighttime}}\cdot \mathit{f}({\mathit{t}}_{kj,\mathit{nightime}})},$$(1) is the accessibility of each demand point to the shelter through improved two-step floating catchment area method. Equation (2)(2) $$\mathit{f}({\mathit{t}}_{ij})=\left\{\begin{array}{c}\frac{{\mathit{e}}^{-1/2{({\mathit{t}}_{ij}/{\mathit{t}}_0)}^2}-{\mathit{e}}^{-1/2}}{1-{\mathit{e}}^{-1/2}}{\mathit{t}}_{ij}\le {\mathit{t}}_0\\ \\ \\ 0{\mathit{t}}_{ij}>{\mathit{t}}_0\mathit{ }\end{array}\right\}\mathit{f}({\mathit{t}}_{kj})=\left\{\begin{array}{c}\frac{{\mathit{e}}^{-1/2{({\mathit{t}}_{kj}/{\mathit{t}}_0)}^2}-{\mathit{e}}^{-1/2}}{1-{\mathit{e}}^{-1/2}}{\mathit{t}}_{kj}\le {\mathit{t}}_0\\ \\ \\ 0{\mathit{t}}_{kj}>{\mathit{t}}_0\mathit{ }\end{array}\right\},$$(2) represents time attenuation function. Equations (3)(3) $$\alpha =\frac{\mathit{C}}{\mathit{D}},$$(3) and (4)(4) $$Z=\sum _{i=1}^m{\left(A_{\mathit{i}}-\alpha \right)}^2,$$(4) represent the equity measurement of refuge resources available for the demand point under one shelter layout scheme.

The spatial accessibility evaluation is the basis for the equity maximization model. Evaluation methods of facility accessibility include the nearest distance method (Guo and Xu 2011), the proportion method, the gravity model (Wu et al. 2017) and the two-step floating catchment area method (Ni et al. 2015). Among them, two-step floating catchment area method is a accessibility measurement approach based opportunity accumulation. It takes the supply and demand points as the center to complete floating search twice (Wang and Luo 2005). In the current study, we improve the traditional two-step floating catchment area approach via three steps to complete the accessibility assessment from each demand point to the shelters under each layout scheme. First, the time attenuation function is introduced on traditional two-step floating catchment area approach, which assumes that the interaction between the demand point and shelter descends with the increasing evacuation time. Second, variations in population during the day and night are considered. Equation (1)(1) $${\mathit{A}}_i=0.5\times \sum _{\mathit{j}=1}^{\mathit{n}}\frac{{\mathit{C}}_j\cdot \mathit{f}({\mathit{t}}_{ij,\mathit{daytime}})\cdot {\mathit{y}}_j}{\sum _{\mathit{k}=1}^{\mathit{m}}{\mathit{D}}_{k,\mathit{daytime}}\cdot \mathit{f}({\mathit{t}}_{kj,\mathit{daytime}})}+0.5\times \sum _{\mathit{j}=1}^{\mathit{n}}\frac{{\mathit{C}}_j\cdot \mathit{f}({\mathit{t}}_{ij,\mathit{nighttime}})\cdot {\mathit{y}}_j}{\sum _{\mathit{k}=1}^{\mathit{m}}{\mathit{D}}_{k,\mathit{nighttime}}\cdot \mathit{f}({\mathit{t}}_{kj,\mathit{nightime}})},$$(1) considers the dynamic variations of population to improve the accessibility assessment accuracy. Equation (1)(1) $${\mathit{A}}_i=0.5\times \sum _{\mathit{j}=1}^{\mathit{n}}\frac{{\mathit{C}}_j\cdot \mathit{f}({\mathit{t}}_{ij,\mathit{daytime}})\cdot {\mathit{y}}_j}{\sum _{\mathit{k}=1}^{\mathit{m}}{\mathit{D}}_{k,\mathit{daytime}}\cdot \mathit{f}({\mathit{t}}_{kj,\mathit{daytime}})}+0.5\times \sum _{\mathit{j}=1}^{\mathit{n}}\frac{{\mathit{C}}_j\cdot \mathit{f}({\mathit{t}}_{ij,\mathit{nighttime}})\cdot {\mathit{y}}_j}{\sum _{\mathit{k}=1}^{\mathit{m}}{\mathit{D}}_{k,\mathit{nighttime}}\cdot \mathit{f}({\mathit{t}}_{kj,\mathit{nightime}})},$$(1) is the accessibility of each demand point to the shelter through improved two-step floating catchment area method, which considers the variations in population during the day and night are considered. The ‘demand’ in the traditional two-step floating catchment area method is calculated by resident population value. Nevertheless, because of the large number of commercial land in the central regions, the population varies greatly between the day and night, while the occurrence time of hazards is uncertain. The accessibility assessment of shelters should thus simultaneously consider both night and daytime conditions. Equation (1)(1) $${\mathit{A}}_i=0.5\times \sum _{\mathit{j}=1}^{\mathit{n}}\frac{{\mathit{C}}_j\cdot \mathit{f}({\mathit{t}}_{ij,\mathit{daytime}})\cdot {\mathit{y}}_j}{\sum _{\mathit{k}=1}^{\mathit{m}}{\mathit{D}}_{k,\mathit{daytime}}\cdot \mathit{f}({\mathit{t}}_{kj,\mathit{daytime}})}+0.5\times \sum _{\mathit{j}=1}^{\mathit{n}}\frac{{\mathit{C}}_j\cdot \mathit{f}({\mathit{t}}_{ij,\mathit{nighttime}})\cdot {\mathit{y}}_j}{\sum _{\mathit{k}=1}^{\mathit{m}}{\mathit{D}}_{k,\mathit{nighttime}}\cdot \mathit{f}({\mathit{t}}_{kj,\mathit{nightime}})},$$(1) considers the dynamic variations of population to improve the accessibility assessment accuracy. 0.5 represents the average value of population during daytime and night.

Third, we improve the two-step floating catchment area method via authentic evacuation time cost. The traditional traffic cost is calculated by the road network distance, which leads to overestimations in the shelter accessibility (Geng et al. 2020; Jing et al. 2020). This study adopts the API (Cheng et al. 2016) of the Baidu Map open platform based on Python to obtain the travel time from every demand plot to every shelter. The transportation mode is set as ‘walking’. Compared with the European and road network distances used in traditional accessibility assessments, the walking time determined here can reflect authentic evacuation routes, as well as dynamic traffic conditions. Considering the uncertain occurrence time of hazards and the evacuation time, the request time in the API is set to both day and night. This is coincident with the difference in the population distribution during the day and night. The time cost input into the two-step floating catchment area method is determined as time cost of citizens from demand area i to shelter j at 4 am and 4 pm.

The equity maximization model is crucial for the layout optimization of shelters and overcomes the absence of equity in traditional models. In particular, it realizes the two policy goals of refuge facilities: efficiency and equity.

Here, we take the total evacuation time, quantity of new shelters and maximal equity of refuge resources as target functions to build a three-objective optimization model:

Objective functions: (5) $$\mathit{min\hspace{0.17em} }{z}_1=\sum _{j\in J}{y}_j,$$(5) (6) $$\mathit{min}{z}_2=0.5\sum _{i\in I}\sum _{j\in J}(x_{ij}·D_{i,\mathit{daytime}}·t_{ij,\mathit{daytime}})+0.5\sum _{i\in I}\sum _{j\in J}(x_{ij}·D_{i,\mathit{nighttime}}·t_{ij,\mathit{nighttime}}),$$(6) (7) $$\mathit{min\hspace{0.17em} }{z}_3=\sum _{i=1}^m{(A_i-\alpha )}^2,$$(7) and constraint formulas: (8) $$\sum _{j\in J}{x}_{ij}\ge 1\forall i\in I,$$(8) (9) $$x_{ij}\le n_{ij}\forall i\in I\forall j\in J,$$(9) (10) $$\sum _{i\in I}(x_{ij}\cdot D_{i,\mathit{daytime}})\le y_j\cdot C_j\forall j\in J,$$(10) (11) $$\sum _{i\in I}(x_{ij}\cdot D_{i,\mathit{nighttime}})\le y_j\cdot C_j\forall j\in J,$$(11) (12) $$\sum _{j=1}^n{p}_{i,j}^{\mathit{daytime}}=D_{i,\mathit{daytime}}\forall i\in I,$$(12) (13) $$\sum _{j=1}^n{p}_{i,j}^{\mathit{night}}=D_{i,\mathit{nighttime}}\forall i\in I,$$(13) (14) $$y_{1,2,..A}=1,$$(14) (15) $$x_{ij}=0or1y_j=0or1,$$(15) where D_i,daytime and D_i,nighttime are the quantity of refugees during the day and night; $p_{i,j}^{\mathit{daytime}}$ and $p_{i,j}^{\mathit{nighttime}}$ are the quantity of refugees transformed from demand plot i to shelter j during day and night; C_j is the maximum capacity limitation of shelter j; J denotes the congregation of shelters; I denotes congregation of the demand plots; n_ij denotes the distance limitation matrix. The research object is the resident emergency aggregate shelter, and the maximal service range limitation is 1500 m (China MOHURD 2015). The average evacuation speed is set as 1.27 m/s (Wang 2012), and thus the time limitation is 1181s. If the evacuation time from demand plot i to shelter j exceed time limitation, the numerical value of n_ij is 1; otherwise, the value is 0. x_ij and y_j are the 0–1 variables. If shelter j is determined as selected shelters, y_j = 1, otherwise y_j = 0; if residents from demand plot i takes refuge at shelter j, then x_ij = 1, otherwise x_ij = 0.

Equations (5–7) are the three target functions, and represent objectives to minimize the square deviation of accessibility for each shelter, the number of new shelters, and total (day and night) evacuation cost.

Equation (8)(8) $$\sum _{j\in J}{x}_{ij}\ge 1\forall i\in I,$$(8) guarantees that all demand plots should be serviced by at least one shelter, while Equation (9)(9) $$x_{ij}\le n_{ij}\forall i\in I\forall j\in J,$$(9) guarantees that the crowd enter shelters within evacuation time limitation. Equations (10(10) $$\sum _{i\in I}(x_{ij}\cdot D_{i,\mathit{daytime}})\le y_j\cdot C_j\forall j\in J,$$(10) and 11)(11) $$\sum _{i\in I}(x_{ij}\cdot D_{i,\mathit{nighttime}})\le y_j\cdot C_j\forall j\in J,$$(11) indicates that the quantity of refugees received by shelter cannot transcend its capacity limitation. Equations (12(12) $$\sum _{j=1}^n{p}_{i,j}^{\mathit{daytime}}=D_{i,\mathit{daytime}}\forall i\in I,$$(12) and 13)(13) $$\sum _{j=1}^n{p}_{i,j}^{\mathit{night}}=D_{i,\mathit{nighttime}}\forall i\in I,$$(13) ensure that all refugees at all demand points are evacuated completely during both the day and night. Equation (14)(14) $$y_{1,2,..A}=1,$$(14) denotes that only current shelters should be selected, where A represents shelters already built and M denotes the candidate shelters. Formula (15)(15) $$x_{ij}=0or1y_j=0or1,$$(15) is a constraint for 0–1 variables.

Cyclic assignment model

In traditional location models, one demand block only corresponds to one shelter, and refugees from one demand block arrive at the same shelter. This results in low evacuation efficiency or a waste of refugee resources due to the mismatch between the population distribution and maximum capacity of shelters. Therefore, this study make the hypothesizes that the multiple demand blocks correspond to multiple shelters. More specifically, refugees from one demand block are assumed to have the authority to choose different shelters rather than strictly following the assignment. The number of refugees from demand block i are assigned to shelters through gravity model within the service radius (Figure 2). We take the assignment process during the day as a case to demonstrate the process.

Figure 2. Schematic diagram of circular distribution model.

According to the distribution principle of the gravity model, the refugee quantity is allocated based on the evacuation time and capacity of the shelter. The shorter the evacuation time, the larger capacity, and thus more refugees will be received by the shelter. Assuming that there are k shelters within the evacuation time limitation of demand plot i, the weighted factor from demand plot i to shelter j is calculated as follows: (16) $${\mathit{w}}_{\text{ij},\text{daytime}}=\frac{{\mathit{C}}_{\mathit{j}}}{{\mathit{t}}_{\mathit{ij},\mathit{daytime}}}/\sum _{\mathit{j}=1}^{\mathit{k}}\left(\frac{{\mathit{C}}_{\mathit{j}}}{{\mathit{t}}_{\mathit{ij},\mathit{daytime}}}\right),$$(16) where C_j is maximum capacity of shelter j; and t_ij,daytime denotes the evacuation time from demand plot i to shelter j during day, derived from Baidu Map API service. The transportation mode is set to ‘walking’. (17) $$p_{ij,\mathit{daytime}}=\left[w_{ij,\mathit{daytime}}{D}_{i,\mathit{daytime}}\right],$$(17) where w_ij,daytime denotes the weight factor for allocation; D_i,daytime represents the quantity of refugees from demand plot i during the day. When the quantity of residents entering shelter j arrives at maximum capacity at time t1 j, shelter j is closed down. At present, refugees who have not completed the evacuation enter the next cycle. This cycle continues until all refugees enter a shelter.

Evaluation model for the Pareto optimum solution

The difficulties in determining the multi-objective shelter layout scheme is described in the introduction. To overcome these bottlenecks, it is necessity of evaluating the Pareto optimum solutions by combining several important factors.

Moreover, the allocation and transportation of emergency supplies is a critical issue in emergency rescue and disaster management. The inefficient distribution of relief supplies will prevent victims from receiving effective assistance. This study introduces the storage points of emergency supplies in the spatial configuration of shelters to establish a post-disaster linked interactive system (Figure 3). We take the accessibility from shelters to the emergency supply storage points and the construction costs in each layout scheme as the evaluation indexes to select the optimal layout scheme from a group of Pareto optimum solutions.

Figure 3. Post-disaster interactive system of shelters and emergency supply centers.

Accessibility assessment from shelters to emergency supply storage points

We improve the traditional two-step floating catchment area method, which has function of measuring spatial accessibility from shelters (based on the Pareto optimum solution in Section multi-objective layout optimization model of shelters) to the emergency supply storage points. In particular, the accessibility evaluation method is improved based on two steps. First, we consider the differences in population between the day and night. The quantity of refugees received by each shelter during the day and night is determined following the method described by mobile signaling data. The equation used to calculate the accessibility assessment based on the day and nighttime populations is then established. Second, evacuation time is a preferable measurement index than the Euclidean distance for accessibility assessment. Distance is a static and unidirectional index. The majority of research obtains travel time through driving speed of different hierarchy on roads and subsequently employing the OD matrix tool through ArcGIS to derive the minimum evacuation time between the starting points and destinations. The evacuation speed is determined through limitation on hierarchy of roads (Jia et al. 2014) or actual travelling speed (Tao et al. 2014). However, the method cannot consider the road congestion varied with time, and only depends on distance and driving speed, which is also a static method. We adopt the Baidu Map APIs with the mode of transportation set as ‘vehicle’ to accurately determine the transportation time from shelters to emergency supply storage points during the day and night. This method considers the real time traffic conditions and adopts a dynamic route selection process.

According to the improved two-step floating catchment area method, the average value of accessibility coefficient from each selected shelter to the emergency supply storage point in a group of shelter layout schemes is taken as the index to evaluate a group of shelter layout schemes. The evaluation index for the accessibility assessment is calculated using Equations (18–20). (18) $${\mathit{A}}_{\mathit{j}\prime }=0.5\times \sum _{\mathit{e}=1}^{\mathit{s}}\frac{{\mathit{C}}_e\cdot \mathit{f}({\mathit{t}}_{\mathit{j}\prime \mathit{e},\mathit{daytime}})}{\sum _{\mathit{k}=1}^{\mathit{n}\prime }{\mathit{DS}}_{\mathit{k}\prime ,\mathit{daytime}}\cdot \mathit{f}({\mathit{t}}_{\mathit{k}\prime \mathit{e},\mathit{daytime}})}+0.5\times \sum _{\mathit{e}=1}^{\mathit{s}}\frac{{\mathit{C}}_e\cdot \mathit{f}({\mathit{t}}_{\mathit{j}\prime \mathit{e},\mathit{nighttime}})}{\sum _{\mathit{k}=1}^{\mathit{n}\prime }{\mathit{DS}}_{\mathit{k}\prime ,\mathit{nighttime}}\cdot \mathit{f}({\mathit{t}}_{\mathit{k}\prime \mathit{e},\mathit{nighttime}})},$$(18) (19) $$\mathit{f}({\mathit{t}}_{\mathit{j}\prime \mathit{e}})=\left\{\begin{array}{c}\frac{e^{-0.5\times {({\mathit{t}}_{\mathit{j}\prime \mathit{e}}/{\mathit{t}}_0^{\prime })}^2}-e^{-0.5}}{1-e^{-0.5}}{\mathit{t}}_{\mathit{j}\prime \mathit{e}}\le {\mathit{t}}_0^{\prime }\\ \\ \\ 0{\mathit{t}}_{\mathit{j}\prime \mathit{e}}>{\mathit{t}}_0^{\prime }\mathit{ }\end{array}\right\}\mathit{f}({\mathit{t}}_{\mathit{k}\prime \mathit{e}})=\left\{\begin{array}{c}\frac{e^{-0.5\times {({\mathit{t}}_{\mathit{k}\prime \mathit{e}}/{\mathit{t}}_0^{\prime })}^2}-e^{-0.5}}{1-e^{-0.5}}{\mathit{t}}_{\mathit{k}\prime \mathit{e}}\le {\mathit{t}}_0^{\prime }\\ \\ \\ 0{\mathit{t}}_{\mathit{k}\prime \mathit{e}}>{\mathit{t}}_0^{\prime }\mathit{ }\end{array}\right\},$$(19) (20) $$\mathit{ZE}=\frac{{\sum }_1^{{\mathit{N}}_{\mathit{j}\prime }}{\mathit{A}}_{\mathit{j}\prime }}{{\mathit{N}}_{\mathit{j}\prime }},$$(20) where j’ and k’ are the serial numbers of the selected shelters from a shelter location scheme; n’ is the quantity of selected shelters; s is the quantity of emergency supply storage points; and ${\mathit{t}}_0^{\prime }$ is the transportation time limitation. Due to the limitations in data collection across the study area, we regard large supermarkets as emergency supply storage points, and directly distribute emergency supplies to shelters after an earthquake. Based on the ‘commercial service facilities in residential areas of the 15 min life circle’ from the Standard for Urban Residential Area Planning and Design (China MOHURD 2018), the threshold of transportation time of large supermarkets is determined, which also denotes the maximum time cost for emergency supplies transported from large supermarkets to selected shelters. We set the ratio of the pedestrian evaluation speed to vehicle speed as 0.2 (Chen et al. 2013). Therefore, the time threshold under the vehicle transportation is determined as 3 min. t_{j’e.daytime} and t_{j’e,nighttime} are the transportation time from selected shelter j’ to supply storage point e; t_{k’e.daytime} and t_{k’e,nighttime} are the transportation time from selected shelter j’ to supply storage point e; A_j’ is the accessibility from the emergency supply resources to selected shelter j’; N_j’ is the quantity of selected shelters; and C_e is the reserve quantity of emergency supply storage point e. It is assumed that emergency supplies include drinking water, cooked food, tents and quilts. We monetize the supply reserve quantity, and take the sum of the monetary values of four types of emergency supplies as the total C_e of each emergency supply storage point e. DS_k’,daytime and DS_{k’,nighttime} are the quantity of refugees received by shelter k’ during the day and night; and f(t_j’e) and f(t_k’e) are the time attenuation function (Wang 2012). According to relevant research on the accessibility evaluation of commercial facilities (Wang et al. 2016; Gan et al. 2022), we employ the Gaussian attenuation function to measure accessibility level from large supermarkets (i.e. emergency material storage sites) to selected shelters. ZE is the average accessibility measurement value from the selected shelters to the emergency supply storage points under a group of shelter layout schemes.

Selection of Pareto optimum solution

Equation (21)(21) $${\mathit{E}}_{\mathit{p}}=\frac{{\mathit{\text{ZE}}}_{\mathit{P}}-\min ({\mathit{\text{ZE}}}_1,{\mathit{\text{ZE}}}_2,\dots ,{\mathit{\text{ZE}}}_{\mathit{Q}})}{\max ({\mathit{\text{ZE}}}_1,{\mathit{\text{ZE}}}_2,\dots ,{\mathit{\text{ZE}}}_{\mathit{Q}})-\min ({\mathit{\text{ZE}}}_1,{\mathit{\text{ZE}}}_2,\dots ,{\mathit{\text{ZE}}}_{\mathit{Q}})}+\frac{{\mathit{J}}_{\mathit{P}}-\min ({\mathit{J}}_1,{\mathit{J}}_2,\dots ,{\mathit{J}}_{\mathit{Q}})}{\max ({\mathit{J}}_1,{\mathit{J}}_2,\dots ,{\mathit{J}}_{\mathit{Q}})-\mathit{min}({\mathit{J}}_1,{\mathit{J}}_2,\dots ,{\mathit{J}}_{\mathit{Q}})},$$(21) is established to determine the optimal layout scheme of shelters: (21) $${\mathit{E}}_{\mathit{p}}=\frac{{\mathit{\text{ZE}}}_{\mathit{P}}-\min ({\mathit{\text{ZE}}}_1,{\mathit{\text{ZE}}}_2,\dots ,{\mathit{\text{ZE}}}_{\mathit{Q}})}{\max ({\mathit{\text{ZE}}}_1,{\mathit{\text{ZE}}}_2,\dots ,{\mathit{\text{ZE}}}_{\mathit{Q}})-\min ({\mathit{\text{ZE}}}_1,{\mathit{\text{ZE}}}_2,\dots ,{\mathit{\text{ZE}}}_{\mathit{Q}})}+\frac{{\mathit{J}}_{\mathit{P}}-\min ({\mathit{J}}_1,{\mathit{J}}_2,\dots ,{\mathit{J}}_{\mathit{Q}})}{\max ({\mathit{J}}_1,{\mathit{J}}_2,\dots ,{\mathit{J}}_{\mathit{Q}})-\mathit{min}({\mathit{J}}_1,{\mathit{J}}_2,\dots ,{\mathit{J}}_{\mathit{Q}})},$$(21) where E_p is the evaluation index of shelter layout scheme group p, and is jointly determined by average accessibility of each shelter to the emergency supply storage points and construction cost of shelters. E_p also represents the evaluation index of the multi-objective shelter layout scheme. Q denotes the total quantity of Pareto optimum solutions; ZE_p denotes the average accessibility of shelters to the emergency supply storage points in layout scheme group p; and J_p represents the number of new shelters in layout scheme group p.

Design of a three-stage solution algorithm

This study proposes a three-stage algorithm to solve the three-objective shelter layout optimization model (Figure 4), and includes four steps: (1) minimizing the quantity of shelters; (2) determining feasible layout schemes; (3) obtaining Pareto optimum solutions; and (4) optimizing the shelter location.

Figure 4. Three-stage solution algorithm of the three-objective optimization model.

First, the minimum quantity of new shelters N is determined via designed genetic algorithm, where M denotes the number of candidate shelters.

In the second phase, feasible layout schemes are selected from ${\mathit{C}}_M^N,$ ${\mathit{C}}_M^{N+1},$ ${\mathit{C}}_M^{N+2}$…, ${\mathit{C}}_M^{M-1},$ ${\mathit{C}}_M^M$ schemes through exhaustive method. The feasible schemes should satisfy all refuge demands during the day and night under ultimate capacity and refuge service range limitations. The total evacuation time, the quantity of new shelters and Z value (standard deviation of accessibility from demand points to shelters) of all feasible location schemes are determined.

In the third stage, we compare the shelter quantity, total evacuation time and Z values of all feasible schemes to screen the Pareto optimum solution groups. The optimal shelter layout scheme is then determined from the Pareto optimum solutions, based on the accessibility evaluation results from the shelters to the emergency supply storage points and the number of selected shelters.

The steps adopted by the genetic algorithm to obtain the minimum quantity of shelters is described in the following.

1. Encoding

   We employ binary coding in this study. The relationship between the layout scheme and gene encodes is based on the assumption that each chromosome denotes a shelter layout scheme group. Each chromosome has m gene points, and each gene point denotes a shelter. The chromosome indicating whether the alternative shelter is selected or not is composed of 0 and 1 gene codes: ‘1’ and ‘0’ indicate that alternative shelter j is selected and not selected.

2. Generation of initial population

   This study adopts the random number approach to produce the chromosomes of the initial population via the rand function in MATLAB to obtain a stochastic number with the interval between 0 and 1. If the numerical value of stochastic number is larger than 0.5, the gene point value is determined as 1; if the number is less than 0.5, the gene point value is set as 0. The stochastic method is then duplicated N times to generate an original population with size of N. A very small initial population size will lead to a local optimum, while if size of population is too large, the number of iterations will increase greatly and the calculation efficiency is reduced. Test calculations determined a suitable population size N of 400.

3. Calculate the fitting function value of layout scheme

   The fitting function is the standard for screening of chromosomes. In the first phase, the fitness function is taken as the quantity of the new shelters. The penalty value is introduced to account for the case when layout scheme (chromosome) does not satisfy the constraints. For this case, the fitness value is taken as the quantity of new shelters plus 200.

4. Determination of genetic algorithm operator

   • (4.1) Selection operation

     The roulette rules and the elite retention strategy are employed to conduct selection operation. The fitting function value of chromosome is calculated as the number of selected shelters. Based on the roulette rules, the probability of choosing chromosome is positively proportional to its fitting function value. Moreover, the roulette rule is conducted to directly copy the chromosome which has optimum fitting value of parent generations to the next generation without crossover and mutation operations. This guarantees that the best-fitting value of the new population is better than that of parent population.

   • (4.2) Crossover operation

     We select the uniform multipoint crossover operation for this step, with fixed probability p_c. In particular, stochastic number r_d∈(0,1) is produced; if r_d is less than p_c, the chromosome begins to conduct a crossover operation. For pairs of chromosomes (parent chromosome A and parent chromosome B) in the parent population, random number r_s∈(0,1) is generated. If r_s is less than 0.5, the gene in the progeny chromosome is inherited from parent chromosome A. If r_s is larger than 0.5, the gene is inherited from chromosome B. This is repeated m times for each gene point to get the progeny chromosome.

   • (4.3) Mutation operation

     In addition to inheriting information of the parent chromosomes, the progeny chromosomes generated in the crossover algorithm also mutate with a certain probability, which is be favorable to maintaining genetic diversity of the population. The progeny chromosomes produced in the crossover operation are mutated with fixed probability p_m. Stochastic number r_m ∈(0,1) is produced for a single chromosome. If r_m is less than p_m, the mutation operation will be conducted using the multi-point uniform mutation operator. Stochastic number r_mp∈(0,1) is produced for the gene spot of a chromosome that performs the mutation operation. If r_m is less than 0.5, the gene spot value is changed as the opposite value; otherwise the gene point value remains unchanged. This process is repeated m times until a new chromosome is produced.

5. Suspension of the algorithm

   The iterative termination condition is determined based on the computation time and accuracy requirements. If the genetic evolution iteration times reaches a maximum value, or the fitness function value of the optimal solution in subsequent generations is not improved, the algorithm is terminated. When the terminal condition is satisfied, the chromosome which has minimum fitting function value is calculated as the solution of the first phase. The genetic point position of the chromosome corresponds to the layout optimization scheme, and its fitness function value is the quantity of the newly built shelter. Figure 4 describes the three-stage solution algorithm established in this study.

The meanings of all symbols are summarized in the following:

D_i,daytime, D_i,nighttime: the quantity of refugees during the day and night;

$p_{i,j}^{\mathit{daytime}},$ $p_{i,j}^{\mathit{nighttime}}:$ the quantity of refugees transformed from demand plot i to shelter j during day and night;

J: the congregation of shelters;

I: the congregation of the demand plots;

i, k: serial numbers of all demand plots;

m: the number of demand plots;

n: the number of all shelters;

j: sequence number of candidate shelters;

j’, k’: the serial numbers of the selected shelters from a shelter location scheme;

n’: the quantity of selected shelters;

s: the number of emergency supply storage points;

J_p: the number of new shelters in layout scheme group p;

N_j’: the number of selected shelters;

C_j: maximum capacity limitation of shelter j;

C: the sum of capacities of all selected shelters;

C_e: the reserve quantity of emergency supply storage point e;

n_ij: distance limitation matrix.

y_j: if shelter j is selected, y_j=1, otherwise y_j=0;

x_ij: if residents from demand plot i takes refuge at shelter j, then x_ij=1, otherwise x_ij=0;

w_ij,daytime: the weight factor for allocation;

t₀: time limitation for refugee service.

${\mathit{t}}_0^{\prime }:$ transportation time limitation from emergency supply storage points to shelters.

t_ij,daytime, t_kj,daytime, t_ij,nighttime, t_kj,nighttime: evacuation time from plot i(k) to shelter j during day and night;

t_{j’e.daytime}, t_{k’e.daytime}, t_{j’e,nighttime}, t_{k’e,nighttime}: transportation time from selected shelter j’ (k’) to supply storage point e;

A_i: measurement of refuge resource accessibility at demand plot i;

A_j’: accessibility from the emergency supply resources to selected shelter j’;

D_{k, daytime}, D_{k, nighttime}: the number of refugees during the day and night at demand point k;

D: the total demand scales;

DS_k’,daytime, DS_{k’,nighttime}: the number of refugees received by shelter k’ during day and night;

α: the proportion of the supply scale of refuge resources to demand scale;

Z: the differentiation in accessibility from each demand plot to shelter;

E_p: the evaluation index of shelter layout scheme group p;

Q: the total quantity of Pareto optimum solutions;

ZE_p: average accessibility of shelters to the emergency supply storage points in layout scheme group p.

Empirical study

General situation of study area

Xinjiekou District of Nanjing city is selected as the study area (Figure 5), with a total area of approximately 13.12 km². Xinjiekou District has both a high population density and high construction intensity, with a resident population density of 31400 people/km², and population density during peak hours of 77200 people/km². The high population concentration intensifies the risk of evacuation. Moreover, the industrial activity of the district is high, with industries including trade, finance, commerce and entertainment. This consequently leads to large differences in the population distribution between day and night.

Figure 5. Location and range of study area.

Basic datasets

The basic data required for the layout optimization of shelters include spatial and population data. The spatial data used here compromises evacuation roads, the spatial distribution and effective refuge areas of current and candidate refuge resources, and the spatial distribution of emergency supply storage points (Figure 6). Following the adaptability evaluation of the candidate sites, 69 candidate plots in Xinjiekou District were selected as suitable for the construction of new resident emergency aggregate shelters.

Figure 6. Spatial data of Xinjiekou District, Nanjing city.

The Baidu Map APIs can plan transportation routes and determine the travel time from demand points to shelters considering real-time road conditions. We use Python to implement the Baidu Map APIs and set the transportation mode as ‘walking’ to obtain the travel time from each demand point to each shelter. Compared with distance accessibility through European and road network distances employed in traditional accessibility assessments, the walking time determined via Baidu Map APIs is more accurate. Considering the uncertainty of the occurrence time of hazards and the subsequent evacuation time, the request time in the Baidu Map APIs is during both the day and night. The time cost used in the two-step floating catchment area method is that at 4 am and 4 pm from demand plot i to shelter j. Table 1 reports the time cost from demand areas to all shelters (including completed and alternative shelters). Similarly, the transportation mode is set as ‘vehicle’ to obtain the transportation time from each shelter to the emergency supply point (Table 2).

Table 1. Travel time matrix from demand points to shelters (unit: seconds).

Demand point	Shelter	1	2	3	…	114	115
1		2072/1856	3551/1424	3421/998	…	2597/1299	574/209
2		1851/1772	1908/1299	1064/942	…	2368/1172	739/185
…		…	…	…	…	…	…
404		1089/224	958/450	1273/991	…	2705/656	3276/1695
405		457/160	2074/571	1861/974	…	1321/725	2615/1765

Note: Shelters 1–45 are current shelters and shelters 46–115 are candidate shelters.

Left value is the time cost at 4 am, and right value is the time cost at 4 am.

Table 2. Travel time matrix from shelters to emergency supply storage points (unit: seconds).

Emergency supply storage point	Shelter	1	2	3	…	114	115
1		1684/357	1371/809	1026/267	…	2530/1650	1672/866
2		1395/251	1865/1650	2068/1196	…	1507/488	2681/1635
…		…	…	…	…	…	…
121		1579/584	2164/572	718/258	…	2161/833	1927/1153
122		1263/621	1308/725	2906/1737	…	/921	/941

Note: Shelters 1–45 are current shelters and shelters 46–115 are candidate shelters.

Left value is the time cost at 4 am, and right value is the time cost at 4 am.

Table 3. Quantity of refugees at each demand plot.

Serial number of demand point	Nighttime refuge demands (person)	Daytime refuge demands (person)
1	1897	4390
2	1543	3642
…	789	1964
404	1614	3832
405	1776	4231

The population census is usually taken as the source of population data in previous studies; however, it only considers data that takes streets as statistics units. In order to obtain more accurate population data that takes plots as statistics units, we determine the quantity of refugees in the study area during a 24-hour period using China mobile phone signaling data from April 12th to 25th, 2021. Xinjiekou District is located at the junction of Gulou District, Qinhuai District and Xuanwu District. Through population data provided by Nanjing Statistical Yearbook (2021), the quantity of users for China mobile communication corporation in Gulou District, Qinhuai District and Xuanwu District accounts for 16.67%, 15.46%, and 15.26% of the resident population in these three districts, respectively. Based on these proportions, we correct the population data derived from the mobile signaling data. The population value at 4 am is used as the resident population and refuge demands at night, while the population value at 12 am is used as refuge demands during the day (Figure 7 and Table 3).

Figure 7. Population variations during a 24-hour period in the study area.

Results

According to the genetic algorithm in former section, the first stage of the algorithm randomly established 400 groups of layout schemes as the original population. The probabilities for crossover and mutation manipulation are calculated as 0.85 and 0.5, respectively. Figure 8 demonstrates the numerical relationship between the iterations quantity and value of fitting function. Following trial calculations with 250 iterations, the value of fitting function value converges to 11. This indicates that at least 11 resident emergency aggregate shelters should be newly constructed.

Figure 8. Numerical relationship between iteration number and fitting function value.

In the second phase, feasible layout schemes are selected from ${\mathit{C}}_{69}^{20}$+${\mathit{C}}_{69}^{21}$+${\mathit{C}}_{69}^{22}$…+${\mathit{C}}_{69}^{68}$+${\mathit{C}}_{69}^{69}$ schemes satisfying all refuge demands under the capacity and service distance limitations through exhaustive method (Figure 9).

Figure 9. Feasible scheme of shelter locations.

The exhaustive method then compares the quantity of new shelters, the total evacuation time and Z values of all feasible schemes in Figure 9 to determine a group of Pareto optimum solutions (Figure 10).

Figure 10. Group of Pareto optimum solutions.

We then compare each group of shelter location schemes in each Pareto optimum solution using Equation (21)(21) $${\mathit{E}}_{\mathit{p}}=\frac{{\mathit{\text{ZE}}}_{\mathit{P}}-\min ({\mathit{\text{ZE}}}_1,{\mathit{\text{ZE}}}_2,\dots ,{\mathit{\text{ZE}}}_{\mathit{Q}})}{\max ({\mathit{\text{ZE}}}_1,{\mathit{\text{ZE}}}_2,\dots ,{\mathit{\text{ZE}}}_{\mathit{Q}})-\min ({\mathit{\text{ZE}}}_1,{\mathit{\text{ZE}}}_2,\dots ,{\mathit{\text{ZE}}}_{\mathit{Q}})}+\frac{{\mathit{J}}_{\mathit{P}}-\min ({\mathit{J}}_1,{\mathit{J}}_2,\dots ,{\mathit{J}}_{\mathit{Q}})}{\max ({\mathit{J}}_1,{\mathit{J}}_2,\dots ,{\mathit{J}}_{\mathit{Q}})-\mathit{min}({\mathit{J}}_1,{\mathit{J}}_2,\dots ,{\mathit{J}}_{\mathit{Q}})},$$(21) . This equation takes the average accessibility measurement from the selected shelters to the emergency supply storage points and construction costs as the evaluation indexes of the multi-objective shelter layout scheme. Based on the two indexes, the exclusive optimal shelter layout scheme is selected from the Pareto optimum solution (Figure 11).

Figure 11. Optimal layout scheme of shelters.

Figure 11 presents the spatial distribution and effective refuge area of the newly built shelters. In particular, there are 2, 13, and 4 shelters providing effective refuge areas between 5000 and 10000 m², 10000 and 20000 m², and 20000 and 25000 m², respectively. Moreover, the effective area of one the shelter reaches 26009 m². Figure 12 depicts the evacuation direction of the demand points, representing the connection between the target shelters and demand points. The average evacuation time for all refuges is determined as 12.4 min, with minimum and maximum evacuation times of 7 and 27.6 min, respectively. All refugees complete the evacuation within 30 min. Figure 13 shows the evacuation time frequency distribution for the demand points. A total of 80.3% of the refugees complete the evacuation within 15 min, while 93.6% complete the evacuation within 20 min. The completion time of the crowd evacuation at two demand points exceeds 25 min.

Figure 12. Evacuation direction of each demand point.

Figure 13. Frequency distribution diagram of evacuation time.

Discussion

This study takes resident emergency aggregate shelters (RS) as the research object to establish a high-precision population model that varies across space and time through mobile signaling data. A multi-objective optimization model for shelters layout, with the optimization targets of maximum equity and quantity of new shelters and total evacuation time, is proposed. A three-stage algorithm, namely the ‘genetic algorithm-exhaustive method-screening’, is designed to determine the solutions of the three-objective model. Compared with the traditional multi-objective layout optimization method, the model established in this study has the following advantages:

1. Prediction for refined refuge demands. The evacuation demands based on the resident population exhibit huge variations in high-density urban areas. The population data used in traditional methods cannot reflect such dynamic variations in the population, and thus is unable to satisfy the refuge demands during peak population periods. Through mobile signaling data, this study establishes an accurate population model that varies with time and space to satisfy the refuge demands in the study area.

2. Refugees at demand points are assigned to shelters through the proposed gravity model, unlike the conventional assumption that refugees from one demand area enter just one shelter. This improves the utilization efficiency of refuge resources.

3. A multi-objective layout optimization model of shelter is proposed. The classical shelter layout and evacuation allocation models concentrate on the high efficiency of facility locations, and the optimization objectives include maximizing the population coverage (maximum coverage model), minimizing the quantity of shelters (set covering model), and minimizing the travel costs from the demand points to shelters (P-median model). Traditional methods do not consider the equity in the layout optimization of shelters due to difficulties in quantifying this variable. The shelter optimization approach proposed in this study takes the maximum equity, high efficiency, low construction costs as principles, and the quantity of newly built shelters and the differences in accessibility from demand points to shelters, overall evacuation time as optimization objectives to build a three-objective shelter layout optimization model. Equity is quantified by the standard deviation of the accessibility results from each demand point to each shelter. The proposed equity maximization model represents key progress in the shelter layout optimization and realizes two important policy objectives, namely, efficiency and fairness.

4. This study adopts the standard deviation of accessibility from each demand point to the shelters to measure equity, improving the traditional two-step floating catchment area method through three steps. First, time attenuation function is introduced, which accounts for the case whereby the interaction between the demand point and shelter decreases with the increasing travel time. Second, population variations during the day and night are considered. Because of the large quantity of commercial land in the central urban area, the population varies greatly across time. Thus, the accessibility assessment of shelters should simultaneously consider both night and day. Third, we improve the two-step floating catchment area method through accurate evacuation time costs. Traffic costs are traditionally calculated by the road network distance, which overestimates the shelter accessibility (Geng et al. 2020; Jing et al. 2020). This study employs the Baidu Map APIs to obtain the evacuation time from demand plot to each shelter. Compared with the European and road network distances in the traditional accessibility assessment, the walking time determined in this study is more accurate and is able to reflect authentic evacuation routes, as well as dynamic traffic conditions.

5. Based on the limitations of the solutions determined by the traditional multi-objectives layout optimization model of shelters, this study proposes a methodology to evaluate the Pareto optimum solutions by combining multiple important factors. In particular, we take the average accessibility from shelters to emergency supply storage points and the construction costs of shelters as the evaluation indicators for the Pareto optimum solution. The optimal layout scheme is then selected from a group of Pareto optimum solutions. This methodology can also be used in another region in the world, not only used in the case study region.

Conclusion

This study takes the resident emergency aggregate shelter as the research object to complete the layout optimization of shelters. Based on the mobile phone signaling data, we establish a population metric that varies with space and time. By taking the total evacuation time, the quantity of new shelters and the maximum equity as the optimization objectives, a multi-objective layout optimization model is proposed. The proposed model makes substantial progress on the calculation of refuge demands, evacuation allocation patterns, and the integration of shelter location equity, the selection of Pareto optimum solutions, and the accessibility evaluation of refuge resources.

We take Xinjiekou District of Nanjing city as the research area and design a genetic algorithm to solve the model. A total of 300 location scheme groups are randomly constructed as the original population. The probabilities for crossover and mutation operation are determined as 0.85 and 0.5, respectively. The value of fitting function converges to 11, indicating that at least 11 new shelters should be constructed. In the second phase, given the minimum quantity of new shelters, the feasible schemes which can satisfy the capacity and service distance limitations under all refuge demands during both the day and night, are selected. For the third stage, the quantity of new shelters, total evacuation time and Z value (standard deviations between accessibility from demand points to shelters) of all feasible schemes are compared to determine Pareto optimum solutions. The unique optimum location scheme is determined from Pareto optimum solutions based on accessibility evaluation results from shelters to emergency supply storage points and the quantity of new shelters. This methodology incorporates the humanitarian from emergency supply storage points into the shelter location model to improve model’s scientificity and efficiency of humanitarian relief.

In the proposed layout scheme of the resident emergency aggregate shelter, 2, 13, and 4 shelters provide effective refuge areas between 5000 and 10000 m², 10000 and 20000 m², and 20000 and 25000 m², respectively, with one shelter reaching an effective area of 26009 m². The average evacuation time for all refuges is 12.4 min. The shortest and longest evacuation times are 7 and 27.6 min, respectively. All refugees complete the evacuation within 30 min; 80.3% complete the evacuation within 15 min; and 93.6% complete the evacuation within 20 min. In addition, the completion time of the crowd evacuation at two demand plots is greater than 25 min.

It is worth mentioning that no changes are required to apply this methodology to another region. Future research will extend the time attenuation function used in the accessibility evaluation model beyond the Gaussian function to adopt various function forms (e.g. the kernel density function). Moreover, disaster process simulations can be embodied into the shelter location model to improve the authenticity and rationality of the results.

Acknowledgments

The authors thank the anonymous reviewers for their helpful suggestions.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are available on request from the corresponding author, Yuwen Lu.

Additional information

Funding

This study was financially supported by the China Postdoctoral Science Foundation (2021M701649), Natural Science Foundation of Jiangsu Province (Grant No. BK20220782), National Natural Science Foundation of China (Grant No. 52108053), and the Social Science Foundation of Jiangsu Province (Grant No. 20ZZC001).