A multi-objective optimization based method for evaluating earthquake shelter location–allocation

ABSTRACT Constructing shelters is an important part of earthquake disaster management which involves selecting the location of them and assigning the evacuees to them. For this work, selecting suitable objectives and solution is important. Thus, in this paper, a multi-objective mathematical model with four groups of the objectives, allied with a modified particle swarm optimization algorithm, has been developed to solve the location–allocation problem for earthquake shelter. The four objective groups are: total shelter number (TSN) and total evacuation distance (TED), TSN and total weighted evaluation time (TWET), total shelter area (TSA) and TED, and TSA and TWET. The solutions of the model include the determination of the shelters from the candidates and how to allocate population to them. Then the solutions of the model with four objective groups are given and compared using safety, capacity and investment evaluation index with the case of Chaoyang district of Beijing, China. Related to government's preferences and future city planning, the most suitable model solutions can be chosen to help decide where it is suitable to construct shelters and how to allocate the population.


Introduction
Earthquake, which is one of the major natural disasters, caused 714,770 deaths worldwide from 2000-2016(EM-DAT 2017. To help reduce the losses of people's lives, a critical and useful part of earthquake disaster management is to instruct effective emergency shelters. One vital factor in the performance of emergency shelters is the location of them at which evacuees can respond to and arrive as soon as possible after an earthquake. It is important that at any time when an earthquake disaster occurred, the emergency shelters have been equipped well and are able to house all of the evacuees and provide relives. The optimum locations of emergency shelters depend on the demand distributions and the budgets of the government. To research the earthquake shelters' location-allocation problems, works about different disasters' shelter can be referenced. A significant bulk of literature has been put out that is on the effective locating of disaster emergency shelters and a wide variety of models have been proposed to solve the shelter selection problem (Saadatseresht et al. 2009;Widener and Horner 2011;Li et al. 2012; CONTACT Wei Xu xuwei@bnu.edu.cn; Juan Du juan.du@bnu.edu.cn Bayram et al. 2015). Even though there are so many diverse models being used, the primary goal of any such model is to determine where emergency shelters can be constructed and how the evacuees can be allocated that best assists all given evacuees. Among the literature on shelter location, much of it applies models derived from the principle of P-median, P-centre and covering models (Hakimi 1964(Hakimi , 1965Church and ReVelle 1974;Berman and Krass 2002;Gama et al. 2013;Bayram et al. 2015;Kilci et al. 2015). Models derived from P-median aim to locate shelters with the number of 'P' so as to allocate all evacuees at the least time. This type is used mostly as the vital objective of evacuees is to arrive at shelters as soon as quickly. Also, models with the principle of covering that aims to cover a set of evacuees with limited shelters are used mostly. The other models that are derivatives of P-centre are used less as their objective is minimizing the maximum evacuation time, which is not the global optimum. Early models with these three types only consider one objective that ignores other limitations. Hierarchical and multi-objective models have been developed which take account of different objectives such as government investment and damage risk (Alçada-Almeida et al. 2009;Doerner et al. 2009;Xu et al. 2018).
With computing's growth and metaheuristics' advent, more complex models have been facilitated to develop to solve location and allocation problems of the disaster shelter. For instance, Kongsomsaksakul et al. (2005) solved a hierarchical model to help determine the location and allocation of flood shelters by applying a genetic algorithm (GA). Also, Doerner et al. (2009) andHu et al. (2014) solved the multi-objective model related to hurricane and earthquake shelter, respectively, by using Non-dominated Sorting GA II. Similarly, a number of other metaheuristics, such as particle swarm optimization algorithms (PSO), have been applied in order to solve the location-allocation problem (Marinakis and Marinaki 2008;Yeh 2009;Ghaderi et al. 2012). Also, some metaheuristics have been combined such as a modified PSO that combine PSO and simulated annealing to optimize the location of earthquake shelters (Hu et al. 2012;Zhao et al. 2015).
Recent developments related to shelter location and allocation research include various aspects such as the use and modification of facility selection model (Sherali et al. 1991;Widener and Horner 2011;Chen et al. 2013) to solve the shelter location problems, especially the development of multiple objective models and the sensitivity analysis of the model (Zhao et al. 2017). Another important development is adding the influence of the disasters in the model. For example, Li et al. (2012) develop a hierarchical model to help determine the location of hurricane disaster shelters and which one should be open within the consideration of the hurricane scenarios and the occurrence time of them. Zhao et al. (2015) take into account the effect of different earthquake scenarios to the number of evacuees when solving the problems of shelter location and evacuee allocation for earthquake shelter.
Notwithstanding a number of advancement has been accomplished in aforementioned works, one limitation of them is that for multi-objective models, only Pareto solutions are given without detailed analysis about how to select the most suitable among them (Hu et al. 2014; Rodr ıguez-Esp ındola and Gayt an 2015). Another limitation is related to the application of different objectives within various shelter location-allocation problems. Taking earthquake shelter as an example, evacuees should be allocated to the shelters as soon as possible. Simultaneously, investment should be taken into account. To express these important objectives, different simplifications and assumptions have been proposed (Ye et al. 2012;Chen et al. 2013;Bayram et al. 2015;Zhao et al. 2015).
The aim of research stated in this paper is to develop different multi-objective models with different objectives to help determine where earthquake emergency shelters are constructed and how the evacuees in communities are allocated to them. Also, a modified particle swarm optimization (MPSO) algorithm (Zhao et al. 2015) is applied to solve these models. Furthermore, in relation to the case study area, Chaoyang district in Beijing, China, a qualitative method is used to evaluate the Pareto results generated by different models.

Methodology
The evaluation of earthquake shelter location-allocation in this paper comprises of two stages: (1) Develop multi-objective models and obtained the Pareto-solutions for these models, and (2) evaluate the Pareto-solutions by considering safety, capacity and investment of shelters.

Multi-objective model
To address the location-allocation problems of earthquake emergency shelters, two objectives should be considered, i.e. evacuation efficiency of evacuees and investment for equipping the shelters. To express these two objectives in mathematical models, different simplified methods have been used. For evacuation efficiency, i.e. the feeling of safety after an earthquake, evacuation distance or evacuation time is used while the number of shelters or the area of shelters is used to present the investment. In order to analyse the influence of different objectives of the model, two different representations of these two objectives respectively are shown in Table 1. All of the models with different objectives are subjected to both capacity and distance constraints as shown in Equations (1) and (2), respectively (Zhao et al. 2015). Also, as staying with acquaintances is another vital thing for evacuees and social structure, all affected people belonging to the same community who are familiar with each other would be allocated to the same shelter. Equation (3) ensures that a community can only be allocated to one shelter.
where M is the community's amount in the study area, N is the candidate shelter's amount selected. d jk denotes the total length of all evacuation paths that compose the shortest route for community j to go to candidate shelter k. Y k indicates if candidate shelter k is selected to be as a shelter (1 if selected, 0 if not selected). The variable S k and C k present the area and capacity of candidate shelter k, respectively. The capacity can be calculated by dividing the shelter area by the minimum required area per capita. The minimum required area per capita is set to 1m2 in this study according to the Beijing Municipal Institute of City Planning & Design (2007). The variable v j is the people's evacuation speed in community j. It can be calculated according to the following equation: where v c , v a and v o are the walking speed of children, adults and elderly people, respectively. These values can be given as shown in the work of Gates (2006). Variables p jc , p ja and p jo are the proportions of them in community j respectively. Then, r is used to adjust the evacuation speed which is set to 1.3 in this paper. Further, B jk indicates if community j is allocated to candidate shelter k. B jk can be 1 only when shelter k is selected as an earthquake shelter and community j is allocated to it, i.e. Y k equals to 1 and community j is allocated to it. B jk will be 0 if community j is not allocated to shelter k. The variable D j is the maximal evacuation distance of community j which can be accomplished in the maximum evacuation time of Tmax j with evacuation speed v j . The variable P j is the number of people in community j that will be evacuated to a shelter. In the formulation of total weighted evaluation time (TWET), W jk is the average width of the evacuation road for evacuees in community j to go to candidate shelter k which include one or more paths.

Results evaluation
There are mainly two approaches related to solving multi-objective problems. Transferring multiobjective to a single objective is a simple one that is accomplished by calculating the weighted sum of all objectives. However, for this approach, because of the difficulty to obtain the prior information of the importance for each objective, it is difficult to assign the weight to them. Besides, this approach can only generate one solution that avoids the potential solutions (Todd 1997). Given this limitation, a new approach that is based on Pareto was proposed (Pareto 1896). A set of optimal and non-dominated solutions can be obtained by using this approach. However, too many solutions are provided to decision-makers, which makes it hard to select the best solution. Therefore, evaluating the Pareto solutions combining different considerable factors for problems becomes an essential work. In this paper, safety, capacity and investment of shelters within each solution are considered and evaluated.

Safety evaluation
For earthquake shelters, their safety is one of the key criteria that will be affected by the landslides, debris flows, etc. (Liu et al. 2011;IFRC 2013;Trivedi and Singh 2017). Besides the distance from earthquake faults, the slope of terrain is also an important factor for the safety of shelter. Thus, the safety of each solution is evaluated by calculating the criteria named SI as shown in Equations (5) and (6). For each solution i in the Pareto set, the parameter si i of a solution is the total slope that all communities should go through to their allocated shelters. The higher the value of si, the higher the risk. Then, according to Equation (5), si i is transferred to the slope indicator -SI i . The value of SI i is from 0-1. In Equation (6), slope ji is the slope of community j of solution i, slope jki is the slope of shelter that community j is allocated in solution i. It should be noted that when slope jki is smaller than slope ji , the result of slope jki minus slope ji is 0. SI should be as low as possible for a shelter location-allocation solution to ensure that the evacuees will evacuate through a flat area.
Simultaneously, evacuation time is important for evacuees to keep safe. As the locations and evacuation roads are different for evacuees in different communities, evacuation time differs significantly for them. If there is no balance between the evacuation time of evacuees in different communities, evacuees in one community may be evacuated after a long time while others go to the shelters quickly. To evaluate the balance of evacuation time between different solutions, the criteria named TI are proposed as shown in Equations (7) and (8). In the equation, TIM ji represents the evacuation time of community j to its selected shelter for solution i, TIM i,max and TIM i,min are maximal and minimal evacuation time among all communities for solution i, respectively. The variable ti i shows the average gap of evacuation times among all communities for the solution i. Then, according to Equation (7), ti i is transferred to theTime indicator, TI i . The value of TI i is from 0 to 1. It means that the value of TI of a shelter location-allocation solution should be with a low value to guarantee the equality.
With the consideration of SI and TI, the Safety index can be obtained according to Equation (9), in which a 1 and b 1 are weights whose sum is equal to 1. Their values are set as 0.5 in this study.

Capacity evaluation
The capacity of the shelters provided to evacuees is also important as there should be enough space to store water, food, personal belongings and other stuff (Nappi and Souza 2014). Thus, the utilization of a shelter is vital as it reflects the robustness of the shelter to meet other needs. Because of the differences of community population and shelter capacities, utilization of the shelters differs significantly. First, the balance of the shelter areas should be ensured as the situation is undesirable, in which some shelter areas are occupied completely while others are utilized less. Second, besides enough capacities for the evacuees, there still needs capacities for managers, volunteers, and other necessary living equipment. Thus, the CI criterion is defined as shown in Equations (10) and (11). In Equation (11), the variable UR ki is the utilized ratio of selected shelter k while UR i,max and UR i,min are maximal and minimal utilized ratio of selected shelters, respectively, for a solution i. The variable n is the number of selected shelters for the solution i, a 2 and b 2 are weights of the balance and robustness of a plan whose sum is equal to 1. Their values can be set according to the government performance, and set as 0.5 in this study. Y ki indicates whether the candidate shelter k is selected or not for solution i (1 if selected; 0 if not). Then, according to Equation (10), ci i is transferred to Capacity indicator (CI i ). The value ofCI i is from 0 to 1. The lower value of CI indicates the better utilized balance and capacity robustness of the plan.

Investment evaluation
For government decision-makers, not only should the safety and shelter service capacity be considered, but also, more importantly the investment efficiency of the relief resources and construction should be taken into account. Thus, for each solution i, the criterion COI i can be calculated by using Equation (12) which shows the investment efficiency of a solution. Its value belongs to 0 and 1, which is a transformation of ce i by considering the maximal value of ce i among all solutionsmax (ce). The value of ce i can be obtained using Equation (13). The variable INV i is the value of government investment for solution i that can be obtained by using Equation (14) while EF i is the value of evacuation efficiency of solution i that can be obtained by using Equation (15). In these two equations, a 3 , b 3 , a 4 and b 4 are weighted values between 0 and 1, which are set as 0.5 in this study.
According to Equation (16), a comprehensive index can be obtained with the weighted sum of Safety index, CI and COI. The values of the three weights, a 5 , b 5 and ' are between 0 and 1, and their sum is 1. In this study, all of them are set to 1/3.

Case study
For the purpose of the research work reported in this paper, Chaoyang district, a central district of Beijing, China is selected as a study area. According to the national standard (AQSIQ and SAC 2008), there are three types of earthquake emergency shelters in China. That is, shelter type I, which can serve evacuees for more than 30 days, shelter type II, which can serve evacuees from 10 to 30 days with better services after an earthquake disaster and shelter type III, which can only provide help to evacuees at most 10 days with basic facilities. Shelters of types I and II also have the function of providing safe places for evacuees soon after an earthquake happens, and can be used as a type III shelter in an emergency situation. In Chaoyang district, there are eight earthquake shelters belonging to types of I and II. Assuming that all evacuees in each community will evacuate to the nearest shelters, the number of evacuees that each shelter should accommodate is presented in Table 2. Also, the ratio of the number of evacuees that each shelter should accommodate to its capacity is also shown in Table 2. It can be seen that except Chaoyang Park, other shelters cannot accommodate all evacuees allocated to them. To accommodate all community residents, parks with area more than 3000 m 2 are also selected as temporary shelters (type III shelter) in this study. Considering the earthquake effect, the locations of candidate earthquake shelters should be at least 500 m away from the earthquake faults (Hu et al. 2014). Thus, in total 72 potential shelters are selected in Chaoyang. Figure 1 shows the residential communities, candidate shelters and evacuation path network. Then, the shortest distance between communities and candidate shelters is calculated by using Dijkstra algorithm (Dijkstra 1959). Also, the population data of 463 communities are provided by the Beijing Bureau of Civil Affairs. In various earthquake scenarios with different earthquake magnitudes, the population that needs to stay in the shelter should be different according to the damage caused by an earthquake (Zhao et al. 2014(Zhao et al. , 2015. However, to meet with the worst situation after an earthquake, we assumed that all of the population would need to be evacuated to the shelter.

Results and discussion
As To solve the complexity shelter location-allocation problem that involves different objectives and constraints, and a large number of communities and candidate shelters, a heuristics named MPSO proposed by Zhao et al. (2015) was used in this study. This algorithm combines particle swarm algorithm and simulated annealing algorithm. Also, it uses both von Neumann and global structures to generate better solutions. By using this algorithm, the Pareto solution sets of the four objective groups (OGs; shown in Table 3) are obtained as shown in Section 4.1. Then, the evaluations of the results are shown in Section 4.2.

Location and allocation results
Four OGs of the models are shown in Table 3. To solve the models in our study, the MPSO was applied. Also, the values of parameters are determined by the premier testing works. Then, by solving the models with these four OGs, 17, 19, 21 and 34 Pareto results are obtained, respectively, for them. The Pareto solutions obtained are shown in Figure 2, which are sorted in ascending order of total shelter number (TSN) and total shelter area (TSA), respectively. As displayed in Figure 2, it can be found that the solutions with lower values of TED or TWET (which mean higher evacuation efficiency) have higher values of TSN or TSA (which mean higher government investment).. The marginal evacuation efficiency of the results has a decreasing trend while the government investment increases. Thus, the decision-makers can determine a suitable solution from all of the Pareto solutions generated according to their necessity and the limitations of investments. Furthermore, for TSN, the ranges of total evacuation distance (TED) and TWET of the results are wider than those of TSA. Besides, there is no solution for TSN from 35 to 40 and from 48 to 58 as shown in Figure 2(a) and Figure 2(b). Similarly, there is no solution for TSA from 20,000,000 to 30,000,000 as shown in Figure 2(c) and Figure 2 In Figure 3, the results with the lowest value of government investment for these four models are shown. Similarly, the results with the lowest value of evacuation efficiency for these four models are shown in Figure 4. In both the figures, the allocations of communities to their shelters are shown using connected lines. It can be observed that all of the lines do not intersect for OG (I) and OG (III) with the objective to minimize TSN. However, for OG (II) and OG (IV), some lines intersect  with each other such as in the bottom left corner of the figure. The reason is that route width is considered into evacuation time. Some routes with longer linear distance but wider width, which makes shorter evacuation time, will be a good choice compared to the routes with shorter distance but narrower width. Furthermore, another interesting finding is that the shelter location allocation results are significantly different in Figure 3 for situations with the same expression of evacuation efficiency. However, the results are same for Figure 4(a,b) and 4(c,d). For OG (I) and OG (II) in Figure 3, shelters with large area are selected to minimize the total shelter number. However, for OG (III) and OG (IV), the shelters with smaller area are selected ignoring the number of them. But in Figure 4, as the investment is enough, people in each community can be allocated to the nearest shelter; thus, results with the same expression of evacuation efficiency are same.

Results evaluation
About 100 Pareto solution results in total are obtained for four OGs. By getting rid of the repeated solutions obtained by all of the four models, there are 85 results totally. In this section, the 85 results are compared by doing safety analysis, capacity analysis and investment analysis. Finally, the best results are given to help the decision-maker determine the final scheme.

Safety analysis
There are two indicators that are considered in our study to do the safety analysis of the solutions. One is named SI which is calculated as in Equations (5) and (6) and the other one is TI which can be obtained by using Equations (7) and (8). In order to obtain the values of SI, the slopes of each community and shelter can be calculated using ArcGIS tools with the Digital Elevation Model (DEM) data. With the slope value of each community and shelter, and the allocation solutions, The SI values of the 85 solutions are shown in Figure 5. In Figure 5, it can be seen that the value of SI is fluctuating between 0.3 and 1. Solutions with lower values of SI are centred around solutions 55-60 and 75-80. As the lower value of SI represents better solutions, it can be found that plan 75 is the best solution within the consideration of slope.
Also, Figure 5 presents the values of TI for all of the 85 solutions. It can be found that the values of all solutions are larger than 0.5 and the differences are insignificant for some solutions such as the solutions 6-20. Then, from the figure, plan 13 can be selected as the best choice as it has the lowest value of TI.
As both of them are important, in our study, their weights are set as 0.5 and the Safety index is obtained and shown in Figure 5. Thus, with the consideration of population safety, the plan 61 with the lowest value of 0.4983 should be the best one.

Capacity analysis
Considering the comfort of the evacuees and the capacity of the shelters to meet with the special situation such as sharp increased number of the evacuees from other districts, the capacity analysis should be taken into account. First, the utilized ratio of each shelters within each solutions is obtained, and then, the index of the capacity, named CI as present in Equations (10) and (11) is calculated and shown in Figure 6. It can be found that its range is around 0.4-1. As described in Section 2, the index CI expresses the capacity difference of selected shelters and the ability of shelters to meet with the requirement of storing relief assets and the increasing population. A high value of it represents that the spare capacity of selected shelters is different seriously, which means that some shelters have limited spaces while other shelters have larger spaces. Thus, lower value of CI is better than higher one. So, the solution 21, with the lowest value, CI = 0.46, is the best choice.

Investment analysis
Besides safety and capacity of the shelter plans, the investment of the government should be the other important factor that affects the choice of the earthquake shelter location and allocation plan. The index of the capacity, named COI, is calculated by using Equation (12) with the values of TWET, TED, TSN and TSA. Then, the value of COI of each solution is shown in Figure 7. It can be found that the solution 73, which has the lowest value of COI = 0.27, is the best choice.
After aforementioned analysis, by considering these three indices with the same weight, the value of a comprehensive index of these three indices is obtained and shown in Figure 8. Figure 8 shows that the solution 57, which has the lowest comprehensive index value of 0.51, is the best solution among the 85 location-allocation solutions. Its values of TWET, TED, TSA and TSN are 159,192,075,1,296,071,37,099,529 and 23,respectively.. Also, the space distribution and residential allocation are shown in Figure 9(a).  In Figure 9(a), it can be found that only five current shelters (Chaoyang Park, Taiyanggong Park, Olympic Park, North River Park and Jingcheng Liyuan) are used for the evacuees to be allocated. The population they should cover are 173,381, 481,113, 62,515, 105,594 and 83,397 respectively, which exceed the current capacities of these shelters. Thus, these five candidate shelters' areas still need to be extended more to construct the earthquake shelter.
Then, within all of these 85 Pareto solutions, no one selected all of these current eight shelters. The largest number of the current shelters selected is 7, for Pareto solutions 20, 21, 54, 55, 71 and 72, of which, the best one is Pareto solution 71 in terms of the three indices. The location and allocation are presented in Figure 9(b), in which it can be found that only one current shelter, Jingcheng  Liyuan, is not selected. The number of evacuees these seven shelters, Yuandadu, Chaoyang Park, Taiyanggong Park, Olympics Park, Xinglong Park, Red Towel Park and Northern River Park, should service is 204802,16351,157848,13193,113688,240131 and 14479, respectively. Among these seven shelters, Chaoyang Park, Olympic Park and Northern River Park have had enough capacity to house the evacuees allocated to them while others did not.

Conclusion
In this paper, a comparison between four models with different objectives for earthquake shelter location-allocation is completed. Importantly, a new evaluation method has been applied to select the 'optimum' solution among all solutions, in which the safety, capacity and investment of selected shelters are taken into consideration.
Models and the evaluation method have been applied to a case study area of Chaoyang district in Beijing, China. In relation to this case study, a total of 85 solutions are obtained with the four models. Following this, solutions obtained are evaluated and compared with different indices. The result shows that solution 57 is the best one although only five current shelters are selected. Among all of these 85 solutions, only seven current shelters can be designated at most and the best one among the solutions selecting seven current shelters is solution 71. However, for the solutions 57 and 71, the current shelters still need to extend their capacities to meet with the necessity of all evacuees.
It has been noted that the six parameters are set equally to 0.5 for simplicity in this study. In real case, local government or planner could adjust them according to regional condition, shelter planning and residents' preferences. Future work could consider the models with more affected factors such as the risk of secondary disasters, the willingness of the evacuees and the collapse of the buildings. Also, the future work could focus on the more detailed evaluation of the final solutions.