A hierarchical mathematical model of the earthquake shelter location-allocation problem solved using an interleaved MPSO–GA

Abstract Earthquake disaster management involves determining locations in which to construct shelters and how to allocate the affected population to them. A multi-objective, hierarchical mathematical model, allied with an interleaved modified particle swarm optimization algorithm and genetic algorithm (MPSO–GA), have been developed to solve the earthquake shelter location-allocation problem. From a set of candidate shelter locations, the model first determines which of these should act as emergency shelters and then which should be used as long-term shelters, while simultaneously optimizing the allocation of a population to them. Damage caused to evacuation routes is considered in addition to the number of evacuees and shelter capacity. In terms of the model’s emergency and long-term shelter stages, the objectives are to minimize (i) total weighted evacuation time, and (ii) total shelter area used. An interleaved MPSO–GA applied to the model yielded better results than achieved using MPSO or GA in isolation. For a case study with an earthquake affecting the area of Jinzhan within Beijing’s Chaoyang district in China, results generated present government with a range of solution options. Thus, based on government preferences, choices can be made regarding the locations in which to construct shelters and how to allocate the population to them.


Introduction
Since 1940, the number of floods and storms has followed an upward trend whereas drought and earthquakes has been relatively stable according to the record of EM-DAT (2017). Furthermore, since 1980, the damage caused by these events has increased. Indeed, natural disasters such as earthquakes, floods, storms and hurricanes can cause significant loss of human life and serious injury to people, along with damage, disruption and economic losses. Further, amongst the major types of natural disasters, earthquakes tend to cause the most damage despite occurring less frequently than other types of disaster.
A number of engineering techniques exist to enhance the resilience of buildings and reduce the damage to them caused by earthquakes such as seismic base isolation (Datta 2010) and seismic shock absorption (Lin et al. 2016). However, in some cases, due to the severity of the earthquake, buildings cannot offer protection to people. Decisions as to where to house displaced people and provide them with sufficient provisions to ensure their safety and survival is an important problem to be solved. In order to provide assistance to government decision makers, much research has been carried out related to determining the optimized position of disaster shelters. Most of the studies conducted have solved this problem by modifying site selection models first proposed between 1960-70, such as the P-median model (Hakimi 1964), P-centre model (Hakimi 1965) and covering model (Toregas et al. 1971). These models have been used widely in disaster shelter location problems (Sherali et al. 1991;Gama et al. 2013;Bayram et al. 2015;Kilci et al. 2015). Based on the three site selection models mentioned, hierarchical models (Chang et al. 2007;Li et al. 2011Li et al. , 2012Widener and Horner 2011;Chen et al. 2013) and multi-objective models (Doerner et al. 2009;Saadatseresht et al. 2009; Barzinpour and Esmaeili 2014, Rodr ıguez-Esp ındola and Gayt an 2015) have been developed to solve the shelter location and population allocation problems with each having different objectives. Hierarchical models have two main types, namely bi-level model and general hierarchical model. Bi-level models have been used widely to determine the shelter location and evacuee allocation before hurricanes and floods occur. For example, to solve a flood shelter selection problem, Kongsomsaksakul et al. (2005) proposed a bi-level model to solve the flood shelter location-allocation problem with an upper level objective to minimize total evacuation time lower level objective to ensure each evacuee travels to a shelter as soon as possible. Ng et al. (2010) proposed a hybrid bi-level model for hurricane shelter determination and evacuees' allocation in which the upper level selects shelter locations and the lower level gives the evacuation paths that are selected freely by evacuees. Li et al. (2011) used a bi-level model to solve the hurricane shelter location-allocation problem with a case study in the Gulf Coast region, USA. Similarly, in other studies, in the upper level (called the preparedness level) the location of shelters is determined. Subsequently, in the lower level (called the response level) both evacuees and resources are distributed to shelters. Li et al. (2012) presented a stochastic bi-level model to solve the shelter location-allocation problem in a hurricane scenario. The upper level is divided into two stages in which the locations of shelters are established followed by the determination of opened shelters which are outside the area affected by the hurricane. Similarly, with the work of Ng et al. (2010) and Li et al. (2011), in the second level, the evacuees select their evacuation paths according to the result of the first level in which locations of shelters are determined. General hierarchical models are used to determine the location of shelters of different types. For example, Widener and Horner (2011) developed a hierarchical model with the objective of minimizing the distance from all demand points to their assigned facilities for the hurricane relief point location-allocation problem. In the model, the lower level provides basic relief goods, whereas the upper level provides special relief goods. For an earthquake disaster, Chen et al. (2013) allocate people to three types of shelter with a hierarchical model according to a single objective, i.e. minimizing total distance travelled by evacuees.
Notwithstanding the importance of the aforementioned research, there remains scope for further advances to be made in relation to the earthquake shelter locationallocation problem. Specifically, opportunities exist to make original contributions in terms of model development. Many mathematical models of the earthquake shelter location-allocation problem do not consider the effect of earthquake intensity on evacuation routes in that those in closer proximity to the epicentre will suffer more damage than those further away. Consequently, the rate at which people can evacuate via damaged routes is lower than via those undamaged by the earthquake. For example, the models proposed by Chen et al. (2013) and Hu et al. (2012Hu et al. ( , 2014 only consider evacuation route length while ignoring the damage caused to them by an earthquake. Salman and Y€ ucel (2015) do consider path damage when solving the shelter location-allocation problem with a path having one of two possible states, i.e. failure or available. However, in real earthquake situations, paths may not be damaged entirely. Some researchers have studied road damage after an earthquake, such as Haghighattalab et al. (2010), but their work focussed on the damage assessment using satellite images. In the context of mathematical modelling, road damage estimation has received limited consideration. Another limitation of current models relates to the types of shelter considered in earthquake evacuation and how people are allocated to them. According to planning of the Beijing Municipal Institute of City Planning & Design (2007), an earthquake shelter should be defined as being either an emergency shelter (EMS) or a long-term shelter (LTS). As mentioned previously, although there are some works regarding different types of hurricane shelters, such as the work of Widener and Horner (2011) that determines the placement of different hurricane relief goods within a hierarchical model, the authors only take into account the objective of minimizing the travel distance. Also, an earthquake disaster is different from a hurricane which means the findings obtained from the research of hurricane shelter location-allocation problem may not be completely applicable in the shelter location-allocation problem of an earthquake disaster. For the earthquake disaster, research considering different types of shelters is limited. All except Chen et al. (2013) do not account for different earthquake shelter types, but Chen et al. (2013) allocate people to different types of shelter according to a single objective, i.e. minimizing total distance travelled by evacuees. However, there is a need to consider other objectives such as minimizing total evacuation time and total shelter area used.
Recognizing the scope for improvement in existing earthquake shelter locationallocation models, the research reported in this article is aimed at developing and using a new multi-objective hierarchical model with two stages named the EMS stage and LTS stage to determine the location of EMSs and LTSs from a set of candidate shelters, along with allocating the affected population to them, taking into account the damage caused by an earthquake to evacuation routes based on proximity to the epicentre. Also, the population of a community is divided into sub-communities and the number of people in an EMS is divided into different groups to be evacuated. This model will be applied with a case study of the Jinzhan area within the Chaoyang district of Beijing in China.
In the research presented in this article, it is important to note that although the area considered is 50 km 2 , which may be viewed as relatively small -scale, the problem includes detailed consideration of path damage and dividing communities into sub-communities. That is, the problem considered is complex in that it involves 64 sub-communities in the EMS stage and 61 evacuee groups (EGs) in the LTS stage. Moreover, the model involves two different objectives, namely to minimize total weighted evacuation time and total shelter area, and includes capacity constraints and distance constraints. Thus, heuristic optimization algorithms may be applicable to solve the problem within a reasonable time such as genetic algorithms (GAs) (Goldberg 1989), particle swarm optimization (PSO) (Kennedy and Eberhart 1995) and simulated annealing (SA) (Kirkpatrick et al. 1983). In terms of their application, GAs have been used to solve a variety of different problems including the shelter location-allocation problem (Kongsomsaksakul et al. 2005;Doener et al. 2009;Hu et al. 2014). PSO has also been used in many fields (Jin et al. 2007;Shen et al. 2007;Ai and Kachitvichyanukul 2009) including to solve the location-allocation problem (Marinakis and Marinaki 2008;Yeh 2009;Ghaderi et al. 2012). SA has also been used in fields such as routing (Yu and Lin 2015) and packing problems (Gao 2015); however, in addition, SA has been used with other algorithms to help avoid premature convergence (Ahonen et al. 2014;Mousavi and Tavakkoli-Moghaddam 2013). In terms of solution algorithm development, scope exists to investigate how heuristic optimization algorithms can be used together in order to establish if better solutions to the earthquake shelter location-allocation problem can be obtained than is possible if used in isolation.
To solve the model proposed in this article, an optimization approach has been developed that interleaves the execution of a GA and a modified PSO (MPSO) algorithm, which incorporates SA, thus allowing better solutions to be obtained for the earthquake shelter location-allocation problem than is possible using these algorithms individually. In relation to a case study of the Jinzhan area within the Chaoyang district of Beijing in China, results are presented which provide government with a range of solutions to the earthquake shelter location-allocation problem.

Case study area
The area of Jinzhan within the Chaoyang district of Beijing in China is considered in this case study. The area of Jinzhan is around 50 km 2 , which contains 58,000 population. The location of Beijing in China and the location of Jinzhan within the Chaoyang District of Beijing are indicated in Figure 1(a,b), respectively.
Maps showing the evacuation path network and the locations of Jinzhan's 15 communities and 10 candidate shelters are presented in Figure 5, which was provided by the Key Laboratory of Environmental Change and Natural Disaster of Ministry of Education, Beijing Normal University. In Figure 5, in consideration of the potential damage that can be caused by an earthquake, the locations of shelters should be at least a distance of 500 m from the earthquake faults (Hu et al. 2014).
The road network is as shown in Figure 2a that has 740 paths and 521 path nodes. Based on this road network, ArcGIS mapping software (ESRI 2011) was used to determine the length of each path and subsequently, via Dijkstra's algorithm (Dijkstra 1959), the length of the shortest evacuation route, d kij , from sub-community i of community j to candidate shelter k within the EMS stage of the model (see Appendix I), and from EG i of EMS j to candidate shelter k within the LTS stage of the model (see Appendix II). ArcGIS mapping software was also used to determine the area of each of the 10 candidate shelters shown in Figure 2 (see Appendix I). The width of each path in the case study area was obtained using Google Earth. Further, the adjustment factor (see Equation (4) in Section 3.1.1) associated with each path was determined in relation to the earthquake scenario considered, which is described in Section 2.2.
In the EMS stage of the model, the number of people within each community was divided into a set of sub-communities consisting of up to 1000 people. For example, community '1' consisting of 3848 people was divided into four sub-communities; three of these with 1000 people and the other with 848. Appendix I presents a table in which the first column indicates the index of each community and the number of sub-communities within them (in brackets); the second column indicates the number of people in each of the 15 communities. The same approach as described for the model's EMS stage was used in the LTS stage in relation to dividing the number of people in each EMS into EGs, which may consist of people from different sub-communities.

Earthquake scenario
In this case study, a 6.5 magnitude earthquake with its epicentre in the Tongzhou district of Beijing is considered. In Figure 1(b), the Tongzhou district, in which an earthquake of this magnitude occurred in 1665, is south-east of the Chaoyang district. Figure 3 indicates the location of the earthquake's epicentre in relation Jinzhan.
In terms of earthquake intensity, Equation (1) states the specific form of Equation (2) in Appendix III for Jinzhan obtained using historical data (Zhang et al. 2009).
With a magnitude of 6.5 at its epicentre, the earthquake's intensity at this location is 6.760. At the farthest point (15.537 km) from the epicentre in the affected case study area, the intensity is 6.595. Correspondingly, via Equation (5), the damage ratio varies from 0.519 at the farthest point to 0.552 at the epicentre as shown in Appendix IV. Thus, the upper index of damage level r u , was determined as 6 using the iterative, three-step method described in Appendix III.

Problem formulation and mathematical model
An earthquake shelter can be defined as being either an EMS or a LTS (Beijing Municipal Institute of City Planning & Design 2007). EMSs are equipped to accommodate people for the first day after an earthquake whereas LTSs can be used to house people for up to one month or longer. To determine the number and location of EMSs and LTSs from a set of candidate shelters, as well as the allocation of evacuees to them in the first day and then beyond, a hierarchical mathematical model has been developed and is presented in this article.
An overview of the earthquake EMS and LTS location-allocation problem considered in this article is illustrated in Figure 4.
The hierarchical model, allied with an optimization algorithm, leads to the selection of n ems EMS locations from a set of N candidate shelter (abbreviated as CS) locations; these EMSs are designated as shelters to which evacuees from all subcommunities (abbreviated as SC) of each of M communities (abbreviated as C) are allocated. Initially, while all communities have different locations, all sub-communities within a community have the same location. Also, people within the same subcommunity will be allocated to the same EMS. Subsequently, the locations of n lts LTSs are determined, some of which may have been initially selected as EMSs. Simultaneously, evacuees initially housed in EMSs are divided into groups of evacuees and allocated to LTSs. The hierarchical model is defined by the equations and constraints presented in Equations (1-12). In relation to these equations and constraints, a nomenclature is given in Appendix V.

EMS stage of the model
The two objectives related to the EMS stage of the model are to minimize: (i) total weighted evacuation time for sub-communities to travel from their respective community's location to EMSs (TWET EMS ), f 1 (see Equation (1)); (ii) total shelter area of EMSs (TSA EMS ) to which sub-communities are allocated, f 2 (see Equation (5)).
3.1.1 Total weighted evacuation time for sub-communities to travel from their respective community's location to EMSs (TWET EMS ) The TWET EMS objective function is defined as where n sc;j is the number of sub-communities in community j. In Equation (2), d kij is the length of the shortest evacuation route, which is made up of one or more paths, from sub-community i of community j to candidate shelter k. In this article, as a simplification, all of the people within the same sub-community travel to the shelters to which they are allocated at an average evacuation speed, v ij , which is calculated via Equation (3), where p jc , p ja and p je are the proportions of community j's children, adults and elderly people, respectively, and v c , v a and v e represent the speed of children, adults and elderly people in all communities. In Equation (3), it is noted that an adult will help a child and thus the speed of this adult is same as that of the child. In the case study presented in this article, the speed of children, v c , adults, v a and elderly people, v e , is 1.05, 1.27 and 1.12 m/s, respectively (Gates et al. 2006). Further, the proportions of children, adults and elderly people in each of Jinzhan's 15 communities was determined using population data provided by the Beijing Bureau of Civil Affairs. Due to the limited data available, these proportions were taken to be constant for all communities, i.e. p jc ¼ 0.025, p ja ¼ 0.928 and p je ¼ 0.047. In Equation (2), the quotient involving P ij and W kij is included to adjust the evacuation time due to the congestion of evacuation paths. The parameter P ij represents the number of people within sub-community i of community j. The parameter W kij is the weighted mean width of the evacuation paths that form the entire route taken by sub-community i of community j to candidate shelter k; the term weighted is used to indicate that the length, as well as width, of each path that forms the evacuation route is considered in determining W kij as shown in Equation (4).
In Equation (4), the parameter Q is the total number of paths in an earthquake affected case study area and the variable H qkij indicates whether or not path q forms part of the evacuation route taken by sub-community i of community j to shelter k, i.e. H qkij ¼ 1 if path q forms part of the evacuation route whereas H qkij ¼ 0 if not. Furthermore, w q represents the width of path q which forms part of the evacuation route. To account for earthquake damage to each path making up an evacuation route, the parameter w q is modified using an adjustment factor, a q , which is calculated using Equation (5), where r is the damage level index associated with an earthquake affected case study area ranging from 1 to an upper value, termed r u ; lower values of r signify areas closer to the earthquake's epicentre as shown in Figure 5. The upper damage level index, r u , is determined via an iterative, three-step method as shown in Appendix III. The index l qr and l q represent the length of the path in the damage level area and total path length. The index DR i r and DR o r are inner boundary and outer boundary of damage level area r.
Returning to Equation (2), the variable B kij indicates whether or not sub-community i of community j is allocated to candidate shelter k (providing it has been selected as a shelter), i.e. B kij ¼ 1 if allocated whereas B kij ¼ 0 if not.
3.1.2 Total shelter area of EMSs (TSA EMS ) to which sub-communities are allocated The TSA EMS objective function is defined as where X k indicates whether or not candidate shelter k is selected as an EMS, i.e. X k ¼ 1 if selected, otherwise X k ¼ 0. Further, the parameter S k indicates the available area of candidate shelter k. In the case study described in Section 3, S k is defined as 60% of a shelter's area due to only this proportion being able to be used to house evacuees whereas the remaining 40% is unsuitable (Beijing Municipal Institute of City Planning & Design 2007).

Constraints of EMSs
The constraints associated with the model's EMS stage are related to evacuation time (see Constraint (7)), shelter capacity (see Constraint (8)), and ensuring that a subcommunity of a community can be allocated to only one shelter (see Constraint (9)).
d kij Â B kij ÀD ij 0; 8i ¼ 1; 2; :::; n sc;j ; 8j ¼ 1; 2; :::; n ems ; 8k ¼ 1; 2; :::; N (7) 8i ¼ 1; 2; :::; n sc;j ; 8j ¼ 1; 2; :::; n ems (9) In Constraint (7), the parameter D ij is the maximum evacuation distance that subcommunity i of community j can travel in t ij seconds if moving at speed v ij as defined in Equation (3). The parameter t ij is set to 4600 s in order to ensure that each subcommunity in the case study area can reach at least two candidate shelters. In Constraint (8), the parameter C k is the accommodation capacity of candidate shelter k as an EMS, i.e. the number of evacuees that can be housed in shelter k. This parameter can be determined by dividing S k (see Equation (6)) by the average area occupied by a single person. For an EMS, given that evacuees will stay for only a short period of time, a small area, 1 m 2 , is needed per person (Beijing Municipal Institute of City Planning & Design 2007).

LTS stage of the model
The results of the EMS stage of the model provide the input for the LTS stage, which is then solved to determine the location of LTSs and how the evacuees are allocated to them from EMSs. The model's LTS stage includes two objectives which are to minimize: (i) total weighted evacuation time for groups of evacuees, potentially consisting of people from different sub-communities, to travel from EMSs to LTSs (TWET LTS ), f 3 (see Equation (10)); (ii) total shelter area of LTSs (TSA LTS ) to which groups of evacuees are allocated, f 4 (see Equation (11)).

Total weighted evacuation time for groups of evacuees to travel from EMSs to LTSs (TWET LTS )
The TWET LTS objective function is defined as Equation (10) is similar to Equation (2); however the summation limits differ in that: (a) i varies from 1 to the number of groups of evacuees situated in EMS j, n eg,j , each of which may consist of people from the same or different sub-communities (as opposed to the number of sub-communities of community j, n sc,j , as seen in Equation (2)); (b) j varies from 1 to the number of EMSs, n ems , (rather than the number of communities, M, as seen in Equation (2)).

Total shelter area of LTSs (TSA LTS ) to which groups of evacuees are allocated
The TSA LTS objective function is defined as Equation (11) is similar to Equation (6); however the variable X k indicates whether or not candidate shelter k is selected as an LTS (rather than selected as an EMS as in Equation (6)).

Constraints of LTSs
The constraints associated with the model's EMS stage are related to the shelter capacity constraint (see Constraint (12)), and ensuring that an EG from an EMS can be allocated to only one LTS (see Constraint (13)).
X n ems j¼1 X n eg;j i¼1 P ij Â B kij À C k Â X k À Á 0; 8k ¼ 1; 2; :::; N (12) X N k¼1 B kij Â X k ð Þ¼ 1; 8i ¼ 1; 2; :::; n eg;j ; 8j ¼ 1; 2; :::; n ems (13) Constraint (12) is similar to Constraint (8); however, the summation limits differ as stated in Section 2.2.1, and C k is the accommodation capacity of candidate shelter k as an LTS (rather than the capacity of the shelter as an EMS as in Constraint (8)). The parameter C k can be determined by dividing S k (see Equation (11)) by the average area occupied by a single person in an LTS, which is defined as 3 m 2 (Beijing Municipal Institute of City Planning & Design 2007).

Interleaved MPSO-GA
In this research, due to the multi-objective nature of each stage of the mathematical model, coupled with these objectives being conflicting, an interleaved MPSO-GA combining a GA and MPSO algorithm has been implemented to solve the model. A nomenclature related to these algorithms is given in Appendix VI.
The interleaved MPSO-GA (see Figure 6(a)) begins with an initial population of size 200 being randomly generated, via the INITIALIZE function, followed by the MPSO algorithm (see Figure 6(b)) being executed to solve the location-allocation problem. In each iteration, the solution is compared with the solution generated in the previous iteration, and the non-dominated solution PS is updated with the function NONDOMINATED. After the first one hundred iterations of the MPSO algorithm, and each subsequent iteration, the current Pareto set is assessed against the previous fifty Pareto sets such that if no difference exists between them, i.e. convergence it taken as having occurred, then the execution passes from the MPSO algorithm to the GA (see Figure 6(c)). The GA continues to be executed until no difference exists between the Pareto sets in the same way as described for the MPSO algorithm. Also, at this point, execution passes from the GA to the MPSO algorithm. This process of interleaving the execution of the MPSO algorithm and GA continues until the convergence of the Pareto sets is met simultaneously by both algorithms. It is noted that when changing from one algorithm to the other, the population generated in the last iteration is taken as the initial population of the algorithm to be executed. The decision to compare the current Pareto set with the previous fifty Pareto sets in order to establish if convergence had occurred was determined via experimentation. Figure 6(b) presents the MPSO algorithm (Zhao et al. 2015). With a population size of 200, each iteration, the movement of any particle, u, through the search space is informed by its knowledge of the best position it has occupied so far in terms of objective value(s), p best,u , and the position of the particle with the best objective value(s) so far within (a) neighbouring particles, n best,u , (von Neumann topology) or (b) all particles, g best , (global topology).
For each iteration of the MPSO algorithm, the TWET and TSA objective values associated with each particle are evaluated using the COMPUTEOV function that uses Equations (1) and (5) for the EMS stage and Equations (9) and (10) for the LTS stage. Next, the velocity and position of each particle are updated using the functions UPDATE_v and UPDATE_p. For velocity updating, in this research, based on experimentation, to achieve a balance between exploration and exploitation, the von Neumann topology is used in the MPSO algorithm's first one hundred iterations and thereafter the global topology is used. Subsequently, for each particle, if its current position, p current , is better than its best position so far, p best , then p current becomes p best . However, if p current is worse than p best , then in contrast to a general PSO algorithm, SA is applied to update p best such that a worse position has the potential of being accepted, albeit with a lower probability than a better position. Following the update of p best , Algorithm 2 determines the position of the particle with the best objective value(s) so far within (a) neighbouring particles, n best , (von Neumann topology) or (b) all particles, g best , (global topology). Here, p nns is updated via the UPDATE_n function, which compares the positions of neighbours of a particle, p n , and the positons of particles in the non-dominated set obtained by neighbours so far. Similarly, p gs is updated with via the UPDATE_g function, which compares the positions of all particles, p g , and the positons of particles in the Pareto set.
In relation to n best and g best , the particle selected, via the SELECT function, is that with the largest rectangular area unoccupied by any other solutions as shown in Figure 7's visual representation of the search space. Selecting n best and g best in this manner facilitates a search within its local proximity that may lead to a better solution being found. This selection approach is adapted from the use of crowding distance by Deb et al. (2002). Figure 6(c) presents the GA implemented in this research, which was developed via experimentation. In the GA, with the initial population of size 200 being taken from the last iteration of the MPSO algorithm, each iteration of the GA uses a COMPUTEFITNESS function to evaluate the fitness of each individual in the current population. The fitness of individual u, f u , is calculated using Equation (20), where n is the number of individuals in the population and R u is the rank of individual u based on dominance in relation to the TWET and TSA objectives (for both the model's EMS and LTS stages). In the GA, once each individual's fitness has been evaluated, the fittest 5% of individuals in the current generation, i.e. iteration of the GA, are preserved as P elite using the ELITE function; later in the GA these replace the worst 5% of individuals in the next generation. The procedure of determining the next generation involves using selection, crossover and mutation operations. The selection of individuals from the population, via the SELECTION_p function, involves the use of a fitness sharing method and a roulette wheel based approach (Goldberg 1989). Next, depending on the crossover probability, c, the selected individual can either (a) mate with another individual, via the CROSSOVER function, to produce offspring to be included in the next generation, or (b) be copied directly to the next generation. In the GA developed in this research, experimentation has established that a good value of c is 0.95. In the SELECTION_p mate function, a strategy is applied such that only a sufficiently different individual within a specified proximity in the search space can be selected as a mate for another individual. Again, via experimentation, it was determined that an individual can only be selected as the mate of another individual providing their respective chromosomes are at least 30% different and they are a Euclidean distance of <5,000,000 apart in the TWET-TSA search space. In relation to the CROSSOVER function, a method proposed by Haupt and Haupt (1998) is used in which uniform point crossover is combined with a blend of the genes of two parents to produce two offspring. Once the next generation is fully populated, each individual can be mutated, via the MUTATE function, according to the mutation probability, m, which in this implementation of a GA is set to 0.04 based on an experimental analysis. Furthermore, the fittest 1% of individuals is immune from mutation. As referred to earlier, within the next generation, P next , the worst 5% of individuals, P worst , selected via the WORST function are replaced by the fittest 5% of individuals preserved from the last generation, P elite .

Results and discussion
Initially, this section presents a comparison of the algorithms mentioned in Section 4 demonstrating that the interleaved MPSO-GA yields better solutions to the earthquake shelter location-allocation problem than if the MPSO algorithm or GA were used in isolation. Following this comparison, in the context of the case study area of Jinzhan within the Chaoyang district of Beijing in China, results are presented from the application of the interleaved MPSO-GA to the model's EMS and LTS stage, respectively.

Comparison of the MPSO algorithm, GA and interleaved MPSO-GA
In order to compare the algorithms discussed in Section 4, the EMS stage of the mathematical model was used in which the location of the EMSs is to be determined along with the allocation of sub-communities to them. Fifteen replicates were performed using each algorithm; however, it was observed that no additional Pareto solutions were found beyond seven executions of the MPSO algorithm and GA, and beyond two executions of the interleaved MPSO-GA. For the MPSO algorithm and GA, 2500 iterations were performed each execution as this was found, via experimentation, to give convergence. However, in using the interleaved MPSO-GA, an average of only 1135 iterations was required before convergence of the Pareto sets of the MPSO algorithm and GA was met simultaneously. Figure 8 shows the Pareto solutions obtained using each algorithm.
In Figure 8, it can be seen that predominantly the interleaved MPSO-GA outperforms both the MPSO algorithm and GA in terms of generating Pareto optimal solutions that minimize the TWET EMS and TSA EMS objectives. It is noted that the units  of the TWET EMS objective is weighted seconds, which not only accounts for the distance travelled along the shortest evacuation route at a particular speed, but also considers the damage to these routes and the number of people moving along them. The Pareto optimal set obtained by the GA is smaller than that obtained by the other two algorithms. In contrast, the Pareto optimal set obtained using the MPSO algorithm spans the same solution space, in terms of TWET EMS and TSA EMS , as the interleaved MPSO-GA; although the number of Pareto solutions is less (12 compared against 16) and these are dominated by those obtained using the interleaved MPSO-GA.
Based on the findings of the comparison presented, the interleaved MPSO-GA was used to solve the hierarchical mathematical model of the earthquake shelter locationallocation problem.

Results of the earthquake shelter location-allocation model
This section presents the results obtained by the interleaved MPSO-GA in solving the mathematical model's EMS stage and then the LTS stage. Figure 9 shows the Pareto optimal set obtained which consists of 16 solutions. Figure 9 indicates that if sub-communities are to be evacuated in less time than more EMSs, or larger EMSs nearer to sub-communities, need to be constructed. Taking solutions on the Pareto front marked 'A', 'B' and 'C' in Figure 9 as examples, the location of selected candidate shelters to be used as EMSs and how the sub-communities are allocated to them are shown in Figure 10(a-c), respectively. In Figure  10, it can be seen that all sub-communities belonging to the same community are allocated to the same EMSs. For example, in Figure 10(a), all sub-communities from 10 communities are allocated to candidate shelter 8 whereas those from five communities are allocated to candidate shelter 9. As indicated in Figure 10(b,c), three (1,8,9) and five (1,2,6,8,9) candidate shelters are used as EMSs.
In Figure 11, it can be seen that for Pareto solutions 'A', 'B' and 'C', shelters 8 and 9 are always selected as EMSs. Further, all sub-communities from communities 1, 5, 6, 7, 9 and 11 are always allocated to shelter 8 whereas all sub-communities from communities 8 and 15 are always allocated to shelter 9. As shown in Figure 9, the TWET EMS objective decreases from Pareto optimal solution 'A' to 'B'. The difference between these two solutions, as illustrated in Figure 10(a,b), is that solution 'B' also includes shelter 1 as an EMS. Further, all sub-communities from communities 2, 3 and 13 are allocated to this shelter, rather than shelter 9 as for solution 'A', since it is nearer to them. Similarly, sub-communities from community 4 are allocated to shelter 1 rather than shelter 8. In consideration of the Pareto optimal solution labelled 'C' in Figure 9, shelters 2 and 6 are added to 1, 8 and 9 as EMSs as shown in Figure  10(c). Also, all sub-communities from communities 10, 12 and 14 are allocated to shelters 2 and 6. The solution to the EMS stage of the model forms the input for the LTS stage. Thus, each of the sixteen Pareto solutions for the model's EMS stage shown in Figure 9 will lead to a different Pareto optimal set of solutions to the LTS stage. In this research, the Pareto solution with the lowest TWET EMS for the EMS stage, labelled 'C' in Figure  9 and illustrated in Figure 10(c), was selected as the input for the LTS stage. In this solution, sub-communities from communities are allocated to EMSs as indicated in Table 1. Within the LTS stage of the model, the sub-communities allocated to the five EMS are divided into 61 EGs in each of which the number of people is up to 1000. The number of EGs in each EMS to be allocated to LTSs is shown in Table 1.
In the LTS stage of the model, all 10 original candidate shelters can potentially become LTSs, five of which were selected as EMSs after solving the model's EMS stage. Thus, on solving the LTS stage of the model, an EMS may become an LTS. On applying the interleaved MPSO-GA to solve the LTS stage of the model, the Pareto optimal set containing 19 solutions was obtained as shown in Figure 10.
In Figure 11, the solutions labelled 'A' and 'C' signify those at either end of the Pareto front, and the solution labelled 'B' is located approximately at the mid-point. In comparing solutions 'A' and 'B', both with a TSA LTS <500,000 m 2 , TWET LTS decreases sharply from 26.4 to 6.9 million weighted seconds. However, in comparing solutions 'B' and 'C', as TSA LTS increase from 500,000 m 2 to approximately 1,200,000 m 2 , TWET LTS decreases gradually from 6.9 million weighted seconds to zero. A TWET LTS of zero corresponds with a solution to the LTS stage of the model that is the same as that obtained at the EMS stage; thus EGs are not reallocated from EMSs to LTSs.
For the solutions on the Pareto front in Figure 11 marked 'A', 'B' and 'C', Figure  12(a-c) shows the candidate shelters to be used as LTSs and includes lines indicating which EGs are allocated to them.
In Figure 12, it can be seen that all EGs located in an EMS can be allocated to different LTSs or the same LTS. For example, in Figure 12(a), EGs from EMS 1 are allocated to LTS 3 and 10 whereas all EGs from EMS 2 are allocated to LTS 10. As indicated in Figures 12(b), shelters 6, 8, 9 initially selected as EMSs go on to be selected as LTSs. Also, as shown in Figure 12(c), all shelters initially selected as EMSs go on to be selected as LTSs. Thus, evacuees allocated to shelters in solving the EMS stage of the model are not reallocated in the LTS stage. In contrast, in Figure 12(b), all EGs allocated to EMS 2 are reallocated to LTS 6. Similarly, some EGs in EMS 1 are reallocated to LTSs 6 and 8. Table 2 presents details of the three Pareto solutions discussed. In addition to indicating the index of shelters selected as EMSs and LTSs, Table 2 specifies how many evacuees are allocated to the LTSs and the time taken for them to evacuate from their designated EMS to LTS. It is noted that the evacuation times stated represent the actual time for EGs to travel from their respective EMS to their designated LTS, taking into account the number of evacuees and the earthquake damage to evacuation routes. In Pareto solution A, only two LTSs are selected, namely 3 and 10, to which all evacuees in EMSs are reallocated. All evacuees in EMS 2 and 6 are reallocated to LTS 10, and most evacuees in other EMSs are allocated to LTS 10 but some to LTS 3. There are two advantages of Pareto solution A compared with solutions B and C: (i) the value of TSA is less than that of Pareto solutions B and C, which represents the total cost for LTS is less assuming the constructing cost of shelters is same per square metre; (ii) the two LTSs, i.e. 3 and 10, are located far from the epicentre of the earthquake, which means LTSs of solution A is safer than LTSs of solutions B and C. However, it has a disadvantage in that all of evacuees in EMSs should be reallocated and thus the total evacuation time is longer than for Pareto solutions B and C. For Pareto solution B, only evacuees in EMS 1 and 2 need to be reallocated to LTS 6 and 8. The advantages of this solution are that it provides a balance between evacuation time and shelter area. Further, the evacuation time is less than that for Pareto solution A; however, it is more than for Pareto solution C. In Pareto solution C, all EMSs are assigned as LTSs meaning that all evacuees can remain in their initial shelters. Thus, the advantage of this solution is that the evacuees do not need to reallocated. However, this solution's disadvantages are that the value of TSA, representing of total cost of LTSs construction, is more than that of Pareto solutions A and B and there are two LTSs near to the earthquake's epicentre, i.e. shelter 1 and 9, which will accommodate 26,411 people.

Conclusions
This article presents a new multi-objective, hierarchical mathematical model of the earthquake shelter location-allocation problem. Importantly, the model accounts for damage caused by an earthquake to evacuation routes and the effect of this on the time taken to evacuate from communities to EMSs and, subsequently, from EMSs to LTSs. Furthermore, an interleaved MPSO-GA has been used to solve the model for a particular case study in order to determine the location of EMSs and LTSs along with how evacuees should be allocated to them in the aftermath of an earthquake. In relation to the interleaved MPSO-GA developed, this has been demonstrated to yield better solutions to the earthquake shelter location-allocation problem than using either the MPSO or GA in isolation. The model and interleaved MPSO-GA have been applied to an earthquake scenario in the case study area of Jinzhan within the Chaoyang district of Beijing in China. For this case study, in solving the model's EMS stage, a set of sixteen Pareto solutions was obtained. Following this, taking the Pareto solution with the lowest TWET EMS as input, the LTS stage of the model was solved yielding a set of nineteen Pareto solutions. These solutions present government with a range of options, each of which offers variations in terms of values for the TWET LTS and TAS LTS objectives. Thus, based on government preferences, choices can be made regarding the locations in which to construct earthquake shelters. In this article, although the model is tested with a relatively small scale problem, namely Jinzhan, in order to demonstrate the model and methods used, the model is also suitable for larger scale problems.
Future work could consider the effect of earthquake damage in more detail such as the damage caused to buildings and the effect of this in terms of obstructing evacuation routes. Another aspect of future work could focus on the time of day at which an earthquake occurs given that the distribution of the population in a given geographical area will vary as will the density of traffic. Furthermore, the adherence of evacuees in using evacuation route determined by the authorities should be taken account. Universities (No. B08008) and the National Natural Science Foundation of China (No. 41201547). Length of path q l qr Length of the path q in the damage level area r M e Earthquake magnitude at the epicentre P ij The number of people within sub-community i of community j in EMS level and within EG i of EMS j in LTS level p ja Proportions of adults in community j for EMS level and in EMS j for LTS level p jc Proportions of children in community j for EMS level and in EMS j for LTS level p je Proportions of elderly people in community j for EMS level and in EMS j for LTS level t ij Travel time of maximum distance v a Speed of adults v c Speed of children v e Speed of elderly people v ij Evacuation speed of the people in sub-community i of community j in EMS level and in EG i of EMS j in LTS level W kij The weighted mean width of the evacuation paths that form the entire route taken by sub-community i of community j to candidate shelter k in EMS level and taken by evacuees' group i of EMS j to candidate shelter k in LTS level w q The width of path q a q Adjustment factor using to modified w q a diff A prescribed value D A set fa ru¼nÀ1 Cognitive acceleration coefficient c 2 Social acceleration coefficient f u Fitness of individual u g best The position of particle with the best objective value(s) so far within all particles i Counter for successive interleaved MPSO-GA convergence n The number of individuals in the population n best,k The position of particle k with the best objective value(s) so far within neighbouring particles P Population P elite Fittest individuals P next Next generation P offspring Offspring P worst Worst of individuals p Position of particle in MPSO and individual in GA p best,k The best position particle k has occupied so far in terms of objective value (