A Comparative Study of Four Metaheuristic Algorithms, AMOSA, MOABC, MSPSO, and NSGA-II for Evacuation Planning

: Evacuation planning is an important activity in disaster management to reduce the e ﬀ ects of disasters on urban communities. It is regarded as a multi-objective optimization problem that involves conﬂicting spatial objectives and constraints in a decision-making process. Such problems are di ﬃ cult to solve by traditional methods. However, metaheuristics methods have been shown to be proper solutions. Well-known classical metaheuristic algorithms—such as simulated annealing (SA), artiﬁcial bee colony (ABC), standard particle swarm optimization (SPSO), genetic algorithm (GA), and multi-objective versions of them—have been used in the spatial optimization domain. However, few types of research have applied these classical methods, and their performance has not always been well evaluated, speciﬁcally not on evacuation planning problems. This research applies the multi-objective versions of four classical metaheuristic algorithms (AMOSA, MOABC, NSGA-II, and MSPSO) on an urban evacuation problem in Rwanda in order to compare the performances of the four algorithms. The performances of the algorithms have been evaluated based on the e ﬀ ectiveness, e ﬃ ciency, repeatability, and computational time of each algorithm. The results showed that in terms of e ﬀ ectiveness, AMOSA and MOABC achieve good quality solutions that satisfy the objective functions. NSGA-II and MSPSO showed third and fourth-best e ﬀ ectiveness. For e ﬃ ciency, NSGA-II is the fastest algorithm in terms of execution time and convergence speed followed by AMOSA, MOABC, and MSPSO. AMOSA, MOABC, and MSPSO showed a high level of repeatability compared to NSGA-II. It seems that by modifying MOABC and increasing its e ﬀ ectiveness, it could be a proper algorithm for evacuation planning.


Introduction
Natural disasters are threats to human life and the ecosystem in general. Climate changes as well as environmental changes-e.g., deforestation-increase the frequency and intensity of natural disasters such as hurricanes, floods, and landslides [1]. Such extreme catastrophes cause many losses in lives, affect the economy, and leave many damages to the affected area. However, disaster effects can be reduced if the society is prepared and plans-e.g., for evacuation-are in place.

An Overview of Metaheuristic Algorithms
Multi-objective optimization problems (MOOP) involve more than one objective function that is to be minimized or maximized. An answer to these types of problems is to find a set of solutions that define the best tradeoff between conflicting objectives. In recent decades, there has been a trend in the scientific community to solve MOOPs by using metaheuristic methods over exact methods. A metaheuristic is defined as a procedure or technique designed for finding the approximate solution in a short time (low computation time) [23]. Metaheuristic approaches categorized as population-based metaheuristics are emerged to find optimal solutions through the iterative process of generating a new population through natural selection. According to Fister Jr. et al. [24], evolutionary algorithms or bio-inspired-based and swarm-intelligence-based algorithms are the most interesting and widely used approaches in population-based metaheuristics. GA and its variants represent a group of evolutionary algorithms, while ABC, ACO, and PSO are three approaches grouped in swarm-intelligence-based algorithms. Those four algorithms are commonly used to solve real-world problems [15]. Another category of metaheuristics is physics/chemistry-based algorithms, which mimic certain physical and or chemical phenomena, including for instance electrical charges, temperature changes, and gravity or river systems. Such algorithms solve a problem based on the process of improving a single solution. SA is the commonly used algorithm in this category [25,26]. These five metaheuristic algorithms are all global optimization methods and can solve higher-dimensional problems; they are robust with respect to the complexity of the evaluation of functions. They can easily be adjusted to the problem at hand. On the other hand, although a lot of research has used these algorithms, the question of finding which one is the best suited for a specific problem has not been answered satisfactorily. Furthermore, maintaining the diversity of optimal solutions and premature convergence of solutions to local optima are still crucial to population-based algorithms.
In order to evaluate all categories, this study used the multi-objective version of four approaches, that is NSGA-II to represent evolutionary algorithms, MOABC and MSPSO to represent swarm-intelligence-based, and AMOSA to represent physics/chemistry-based algorithms. A brief review of each approach is discussed in the following.

Archive Multi-Objective Simulated Annealing Algorithm
Archive multi-objective simulated annealing (AMOSA) is a global optimization algorithm adapted from the process of annealing in metallurgy. Bandyopadhyay et al. [27] proposed the AMOSA algorithm based on the principle of the original Simulated Annealing (SA) algorithm [28]. In AMOSA, the Pareto dominance approach is adopted and uses the concept of an archive to store all non-dominated solutions. The archive size is limited with two parameters known as hard limit (HL) and soft limit (SL). The HL is the maximum size of the archive on termination, and it is equal to the number of non-dominated solutions required by the user; while SL is the maximum size to which the archive may be filled before clustering is used. The algorithm starts with the set of solutions randomly initialized and refined in Algorithms 2020, 13, 16 4 of 21 the archive by using a hill-climbing technique. A solution is added in the archive if it dominates the previous one and exceeds the HL. If the archive reaches the SL size, then the well-known single-linkage clustering is used to reduce the size of the archive to HL in order to keep a diversity of non-dominated solutions [29]. In the main loop of AMOSA, three cases can occur in dominance:

1.
The current solution dominates the new solution and k points from the archive dominate the new solution. In this situation, a new solution can be accepted as the current solution with a given probability.

2.
The current solution and the new solution are non-dominating with respect to each other. Here, the domination status of a new solution and members of the archive are checked through three situations: when a new solution is dominated by k points in the archive, the new solution is non-dominating with respect to the points in the archive, and when new solution dominates k points of the archive. 3.
The new solution dominates k points of the archive. Here the new solution is selected as the current solution and also added to the archive, while all the k dominated points in the archive are removed. The process in the main loop is repeated through the number of iterations for each temperature, which is reduced to at each iteration using the cooling rate alpha until the minimum temperature is reached. Thereafter, the process stops and the resulting archive contains the final non-dominated solutions.
AMOSA algorithm is capable of solving problems with many objective functions. It has been used to solve medical and engineering-related problems [30,31], but so far there is no literature on AMOSA applied to solve evacuation problems.

Multi-Objective Artificial Bee Colony Algorithm
Akbari et al. [32] proposed a multi-objective artificial bee colony algorithm (MOABC) based on the standard ABC algorithm developed by Karaboga [33]. Recently, a variant version of MOABC developed based on ABC has been used to solve evacuation problems [34,35]. In this study, the MOABC colony consists of three groups of artificial bees: employed, onlookers, and scout bees. This algorithm generates a number of solutions and works through optimizing them. First, a number of scout bees explore the search space of the problem randomly and generate solutions as the initial population. The quality of the solutions is evaluated (fitness value) and the best solutions are stored in the external memory (archive). The scout bees that have high fitness are selected to act as employed bees. Each employed bee explores the neighborhood to update its position. Onlooker bees select a solution with a high amount of fitness from the neighborhood of employed bees. A new scout bee makes a new generation of the solution if the onlooker failed to update the quality of the solution. Then, the fitness values of all bees are compared to select the best solution and store it to the archive. The Pareto-based approach proposed by Deb et al. [36] has been used to rank the non-dominated solutions into Pareto fronts. The archive is updated by non-dominated solutions, at each iteration. The MOABC algorithm terminates when the termination conditions are met, and the archive returns the final best solutions as output.

Multi-Objective Standard Particle Swarm Optimization Algorithm
Particle swarm optimization (PSO) is a population-based metaheuristic algorithm introduced by Kennedy and Eberhart [37]. PSO is a swarm intelligence algorithm inspired by the social behavior of bird flocks, fish school. The algorithm has many variants due to its flexibility and robustness in terms of updating the way the velocity of the particle is updated [38,39]. This velocity is the speed of a particle which is used to find its next position in search space. A particle updates its position through topological relationships in the neighborhood. The links between particles facilitate to share information about the previous best position of particles from one to another. PSO has been adapted in many studies related to evacuation planning [40][41][42][43].
In this study, we used the recent standard PSO (SPSO) that was proposed to provide common procedures and guidance to improve the original PSO [44]. However, the proposed SPSO is not for solving complex problems with many objectives. Therefore, we applied a Pareto-based method to evaluate the two objectives simultaneously, and the algorithm is named multi-objective SPSO (MSPSO). MSPSO starts by initializing a random swarm of particles. Each particle is stored in memory with its position, its fitness, and its initial velocity. Then, at each iteration, the velocity of each particle is re-calculated using an equation that contains: (i) the current position of the particle (pbest); (ii) the current velocity; and (iii) the previous best position in the neighborhood (gbest). The fitness is calculated based on new positions of particles found at each iteration. The algorithm can be stopped if a given maximum number of iterations is met.

Non-Dominated Sorting Genetic Algorithm-II
The NSGA-II algorithm proposed by Deb et al. [36] is the best known multi-objective optimization genetic algorithm and widely used to solve evacuation planning [11,[45][46][47]. This algorithm belongs to the class of evolutionary algorithms (EA), in the subclass of genetic algorithms (GA), solving the optimization problem through an evolutional process of the population of individuals.
Initially, a random population P t of size N is initialized, evaluated, and sorted on the basis of non-domination. The fitness of each solution is set to a level number; where level 1 is the best, level 2 is the second-best, and so on. The binary tournament selection, crossover, and mutation operators are applied over P t to generate an offspring population of P t with size N. A solution x i of P t wins a tournament with another solution xj if solution xi has a better rank or if it has the same rank but solution xi has better crowding distance than the solution x j . After generating offspring P t , the main loop of NSGA-II starts by combining the two populations R t = P t + P t and sort R t with the size of 2 N on the basis of non-domination. Then, the elitist selection is applied to select the new population with size N from the highest fronts of R t . This main loop is repeated as many times as needed until the satisfaction of an end criterion (i.e., the number of iterations) is reached. NSGA-II has advantages including its low overall complexity of O MN 2 .

Study Area
Kigali is the capital and the most populated city of Rwanda; it accommodates more than 1.135 million inhabitants on an area of 730 km 2 [48]. Due to its geographical location and characteristics, many areas of the city are prone to natural disasters such as floods and landslides. In a study by MIDIMAR [49], the hazard-prone areas in Kigali were highlighted based on the frequency of natural hazards, the topography of the area, and the total damages from the experience of disasters. We selected our case study area from one of the hazardous areas ( Figure 1).

Data Description
To start the research process, the first and most important stage is data collection and compilation. At this stage, all required secondary data including spatial and non-spatial were provided by the city of Kigali. Those data are shapefiles of routes, slope, land use, urban villages, and a boundary map along with the needed attribute data such as population. The National Institute Statistics of Rwanda (NISR) provided documents and the population data from the fourth Population and Housing Census of 2012.
Kigali is the capital and the most populated city of Rwanda; it accommodates more than 1.135 million inhabitants on an area of 730 km 2 [48]. Due to its geographical location and characteristics, many areas of the city are prone to natural disasters such as floods and landslides. In a study by MIDIMAR [49], the hazard-prone areas in Kigali were highlighted based on the frequency of natural hazards, the topography of the area, and the total damages from the experience of disasters. We selected our case study area from one of the hazardous areas ( Figure 1).  Currently, the existing evacuation planning in the city of Kigali does not show a specific place of evacuation. The local authorities are in charge of providing the facilities and locations of safe areas for evacuation and sheltering when a disaster occurs. Thus, the safe areas were selected based on the international standards of evacuation planning for flood and landslide hazards [50]. This includes open spaces, schools, and churches that meet the suitability criteria of being located out of the disaster-prone zones, on gentle slopes and having access to resources, including water sanitation, food, electricity, and toilets.
GIS was used for the preparation and analysis of spatial data. A densely populated region covering an area of 6.9 km 2 , with a population of 176,741, was selected as a study area (Figure 1b). The population data from NISR was aggregated to the level of small blocks. So, there were 1525 small blocks considered as residential/commercial communities and only one population value was considered for each block. A table was created to store each block, its coordinates, the number of evacuees (population size), and the distance to each shelter following the shortest path. Ten shelters were selected and their capacities were calculated based on the crowd density standard considering 3.5 m 2 per person [51]. The table was created to store each candidate shelter, its coordinates, and its capacity to host evacuees. In total, the selected ten shelters have the capacity to host 134,462 evacuated persons. Each shelter might be overloaded due to a large number of evacuees.
A matrix of the shortest distance from each block to each shelter was generated using network analysis in ArcGIS. The origin-destination cost matrix tool was used to compute the minimum distance.

Objective Functions for Evacuation Model
As highlighted earlier evacuation planning, in this study, is a location-allocation problem. We adopted two objective functions proposed by Saadatseresht et al. [11]: 1.
Function to minimize accumulated distance: This objective function aims at allocating each building block to the nearest shelter. 2.
Function to minimize capacity overload: This objective function aims at distributing the overload of the evacuee population among all shelters.
where m represents the number of building blocks; n is the number of safe areas, d ij is the distance between the i th building block and the j th safe area; p ij is the population in the i th building block being evacuated to the j th safe area; and c j is the capacity of the j th safe area for receiving people.

Modeling Metaheuristic Algorithms for Evacuation Planning
This section explains the way the allocation of people from the building blocks (residence, commercial, offices) to the safe areas (shelters) is modeled for each of the four algorithms.
The evacuation problem is solved as an unconstrained problem. The discrete method is used to represent the solution for all four algorithms. Figure 2 shows an example of a discrete encoding of the shelter allocation for a study area consisting of 10 building blocks. In Figure 2, a solution is presented as a list; the size of the list corresponds to the number of building blocks. This list contains elements that correspond to the shelters. Since one shelter can accommodate many people from different places, the elements in the list are repeated (many-to-one assignment).

Modeling Metaheuristic Algorithms for Evacuation Planning
This section explains the way the allocation of people from the building blocks (residence, commercial, offices) to the safe areas (shelters) is modeled for each of the four algorithms.
The evacuation problem is solved as an unconstrained problem. The discrete method is used to represent the solution for all four algorithms. Figure 2 shows an example of a discrete encoding of the shelter allocation for a study area consisting of 10 building blocks. In Figure 2, a solution is presented as a list; the size of the list corresponds to the number of building blocks. This list contains elements that correspond to the shelters. Since one shelter can accommodate many people from different places, the elements in the list are repeated (many-to-one assignment).

Modeling AMOSA
As mentioned above, the AMOSA algorithm was extended from the principal of simulated annealing to handle multiple objectives problems. This extension lies in determining how to calculate the probability of accepting an individual ′ where ( ′ ) is dominated with respect to ( ). The acceptance of new solutions is based on the probability determined by computing the amount of dominance between two solutions and as where M = number of objectives and is the range of the i th objective. A new solution is selected based on the probability computed with the following equation where is the current state and ( , ) and ( , ) are the corresponding energy values of and , respectively [27].
The solutions were generated as demonstrated in Figure 2. Equations (3) and (4) were used to select and sort the non-dominated solutions in the archive. The algorithm stops when the cooling process reaches the predefined low temperature and the maximum number of iterations.

Modeling MOABC
In MOABC, the coding of the population of the bees is equivalent to the coding of the population

Modeling AMOSA
As mentioned above, the AMOSA algorithm was extended from the principal of simulated annealing to handle multiple objectives problems. This extension lies in determining how to calculate the probability of accepting an individual x where f (x ) is dominated with respect to f (x). The acceptance of new solutions is based on the probability determined by computing the amount of dominance between two solutions a and b as where M = number of objectives and R i is the range of the ith objective. A new solution is selected based on the probability computed with the following equation Algorithms 2020, 13, 16 where q is the current state and E(s, T) and E(q, T) are the corresponding energy values of s and q, respectively [27].
The solutions were generated as demonstrated in Figure 2. Equations (3) and (4) were used to select and sort the non-dominated solutions in the archive. The algorithm stops when the cooling process reaches the predefined low temperature and the maximum number of iterations.

Modeling MOABC
In MOABC, the coding of the population of the bees is equivalent to the coding of the population in Figure 2. At the starting stage of the algorithm, a population of scout bees is initialized and each bee represents a food source as an array with the size corresponding to the number of building blocks and composed of 10 repeated indices of the 10 shelters. Modeling the fitness function in the MOABC algorithm is similar to those in AMOSA, using Equations (1) and (2). After the initialization and evaluation of fitness, the best solutions are stored in an external archive (new list). Since the archive contains the best solutions found so far, then each employed bee x id would select a solution from the archive randomly to update it and become v id . The solution is updated through the equations where i represents the food source which is going to be updated, k ∈ {1, 2, . . . , bee}, and d ∈ {1, 2, . . . , D} are randomly chosen indexes. The coefficient w is used to control the influence of the food source k in the production of the new food source.
After evaluating the fitness of employee bees and updating the archive with the best solutions, a roulette wheel selection method is performed to select the onlooker bees for the next generation. The roulette wheel method selects an individual based on the probability p i , found by calculating the proportion of individual fitness f (x i ) in relation to the total fitness of the n population, as shown in Equation (6). Both employed bees and onlooker bees perform the neighborhood search using the expression in Equation (5) [32]. However, this neighborhood search approach is suitable for the continuous problem, and not for the discrete problem. Thus, in this study, we applied a swap method that randomly selects two elements of a solution and interchanges their indexes. Furthermore, a greedy selection method was applied to evaluate the solution with the best value, comparing an existing solution and a new one. By applying this for all employed bees and onlooker bees, a new bee with the best fitness is selected for the next generation. The best solution is stored in the archive at every iteration of the algorithm. Further exploration is carried out by one scout bee that generates a new random solution. The algorithm is terminated when the given termination criterion (maximum iterations) is attained.

Modeling MSPSO
In the MSPSO algorithm used in this research, every possible arrangement of all building blocks to any candidate shelter can be considered as a potential particle in the search space. The MSPSO algorithm looks for a particle location that satisfies the two defined objective functions of evacuation planning. A particle is defined as a solution and initialized randomly (see Figure 2). However, the SPSO algorithm was designed for continuous spaces and with real numbers, while in our case the space of the problem is discrete. To solve this, a rounded value method was used for mapping a discrete problem space to the continuous space and vice versa. The 10 shelters are randomly attributed to the integer values of 1 to 10. The real values generated from updating the position of particles (movements of particles) are rounded in order to obtain integer values between 1 and 10. Figure 3 shows an example of an initial particle in continuous space transformed into discrete space after updating of a particle position. The fitness function is calculated using Equations (1) and (2) Algorithms 2020, 13, 16 9 of 21 and assigned to each particle. A neighborhood topology (ring topology) is used to determine the global best (gbest) for each particle among its neighbors. The algorithm is terminated by attaining the maximum number of iterations.
However, the SPSO algorithm was designed for continuous spaces and with real numbers, while in our case the space of the problem is discrete. To solve this, a rounded value method was used for mapping a discrete problem space to the continuous space and vice versa. The 10 shelters are randomly attributed to the integer values of 1 to 10. The real values generated from updating the position of particles (movements of particles) are rounded in order to obtain integer values between 1 and 10. Figure 3 shows an example of an initial particle in continuous space transformed into discrete space after updating of a particle position. The fitness function is calculated using Equations (1) and (2) and assigned to each particle. A neighborhood topology (ring topology) is used to determine the global best (gbest) for each particle among its neighbors. The algorithm is terminated by attaining the maximum number of iterations.

Modeling NSGA-II
Solving evacuation problem using NSGA-II begins with initializing a population 0 of chromosomes and then initiate the solutions randomly. The coding of a solution in the form of a chromosome is similar to a solution presented in Figure 2. In this study, the number of genes in each chromosome corresponds to the number of building blocks and each gene contains the index of one shelter with repetition. After initialization, the fitness function is evaluated using Equations (1) and (2). The selection of parent chromosomes of the next generation is done using a tournament selection method based on dominance between two individuals. If the two individuals do not inter-dominate,

Modeling NSGA-II
Solving evacuation problem using NSGA-II begins with initializing a population P 0 of chromosomes and then initiate the solutions randomly. The coding of a solution in the form of a chromosome is similar to a solution presented in Figure 2. In this study, the number of genes in each chromosome corresponds to the number of building blocks and each gene contains the index of one shelter with repetition. After initialization, the fitness function is evaluated using Equations (1) and (2). The selection of parent chromosomes of the next generation is done using a tournament selection method based on dominance between two individuals. If the two individuals do not inter-dominate, the selection is made based on crowding distance [52]. This selection technique has also been used by Datta et al. [53] in designing optimal census areas. To generate a new population (offspring), crossover, and mutation operators were applied. The aim of a crossover operator is to exploit the existing best solutions. There are a variety of crossover operators applied in GIS-based genetic procedures [54], and the most used methods include one-point, two-points, and uniform crossover random operators.
Here we used the two-point method. This method randomly selects two crossover points and then swaps the vectors of both parents between the two positions as shown in Figure 4a. The two-point crossover has been applied in [55] for optimizing land use planning.
A mutation operator is used to maintain the diversity from one generation to the next and to prevent the issue of local optimum. Two elements of a chromosome are randomly selected and swapped as showed in Figure 4b. After crossover and mutation operations, the elitism strategy is applied to sort the combined population of parents and offspring using the non-dominated sorting method [36]. the selection is made based on crowding distance [52]. This selection technique has also been used by Datta et al. [53] in designing optimal census areas. To generate a new population (offspring), crossover, and mutation operators were applied. The aim of a crossover operator is to exploit the existing best solutions. There are a variety of crossover operators applied in GIS-based genetic procedures [54], and the most used methods include one-point, two-points, and uniform crossover random operators. Here we used the two-point method. This method randomly selects two crossover points and then swaps the vectors of both parents between the two positions as shown in Figure 4a. The two-point crossover has been applied in [55] for optimizing land use planning. A mutation operator is used to maintain the diversity from one generation to the next and to prevent the issue of local optimum. Two elements of a chromosome are randomly selected and swapped as showed in Figure 4b. After crossover and mutation operations, the elitism strategy is applied to sort the combined population of parents and offspring using the non-dominated sorting method [36].

Comparing and Evaluating the Performances of Algorithms
In this research, the goal of optimization is to find the best combination of building blocks assigned to shelters, with minimum accumulated distance from building blocks to shelters and minimum overload capacity of all shelters. It is assumed that all building blocks will be assigned a safe area.
To evaluate and compare the four algorithms for the given evacuation problem, the different

Comparing and Evaluating the Performances of Algorithms
In this research, the goal of optimization is to find the best combination of building blocks assigned to shelters, with minimum accumulated distance from building blocks to shelters and minimum overload capacity of all shelters. It is assumed that all building blocks will be assigned a safe area.
To evaluate and compare the four algorithms for the given evacuation problem, the different criteria effectiveness, efficiency (convergence trend, execution time), and repeatability were used. The statistical analysis of variance method (Kruskal-Wallis test) [56] was used to allow us to test how each algorithm achieves the best results and to evaluate if there are statistically significant differences between the tested algorithms. The results from the Kruskal-Wallis (KW) test return the Chi-square value and p-value. A high Chi-square indicates the statistical significance of differences, while the p-value determines if the tested hypothesis should be retained or rejected. If the corresponding KW null hypothesis is rejected, the pairwise comparison is done using the Conover-Iman test [57,58].
The effectiveness of the optimization consists of how good the results of each algorithm is. This study compares the effectiveness of four algorithms for evacuation planning to see how each algorithm minimized the two objective functions: the smaller fitness function value, the better performance in terms of effectiveness criteria.
The convergence trend criteria allowed us to evaluate the fitness variation of the algorithm and get information about the speed of the algorithm needed to reach the optimum solution. Execution time helps to evaluate the computational complexity of the algorithm. Since metaheuristic algorithms use randomness to generate initial solutions and to explore search space of feasible solutions, their results are always different from multiple runs. Considering this, to test the repeatability, we run each algorithm thirty times with the same parameters in order to assess their repeatability.

Results of Comparing Algorithms
To compare the performances of the algorithms, we measured and compared four criteria: effectiveness (solution quality), efficiency (convergence speed and execution time), and repeatability.

Parameter Configuration
Initially, each algorithm has a set of parameters that defines the way they perform the optimization. However, to test these parameters is out of the scope of this research and therefore their values were based on the literature. Nonetheless, since each algorithm works in a different way, several pre-runs were executed in order to look for comparable conditions for them. From this exercise, it was noticed that the parameters population size and maximum number of iterations have significant impact on the results and computation time. Therefore, in order to compare the criteria of effectiveness and efficiency, all algorithms were run on an equal number of the population size of 100 and iterated 500. Other parameters were selected based on the original literature of the algorithms [27,32,36,44]. The tested parameters and their initial values are shown in Table 1.
As highlighted in the study by [59], we recommend future studies to further investigate the parameter tuning of the tested algorithms in this study. Parameter tuning analysis aims to obtain the best parameter setting for each algorithm. Also, note that trial and errors along with the experience of the researchers in understanding how these parameter values correlate to the real-world problem being solved are crucial to achieving satisfactory results.

Effectiveness Comparison
This study compares the effectiveness of the four algorithms to see how effective each algorithm optimizes the two defined objective functions. Table A1 shows the average and worst fitness values for both capacity and distance functions (fcapacity, fdistance), as well as the execution times obtained for 30 runs of each algorithm. In all 30 cases, AMOSA was the best one optimizing both objective functions, while MSPSO derived extreme values compare to MOABC and NSGA-II. As shown in Table 2, the Kruskal-Wallis test provided very strong evidence of a difference between the mean ranks of the four methods, optimizing the fcapacity and fdistance. The p-values of both functions are smaller than alpha = 0.05. This means that all algorithms perform differently in terms of optimizing the two objective functions. Table 2. p-values of the Kruskal-Wallis test for evaluating the effectiveness and efficiency of four algorithms in both optimizing capacity and distance function. For effectiveness evaluation in Table 2, the results from the Kruskal-Wallis test shows very strong evidence of a difference (p = 0.000 and p = 0.000) between the mean ranks of the four algorithms in both optimizing fcapacity and fdistance. As shown in Table 3, the pairwise comparison using the Conover-Iman test was carried out to compare the four algorithms and we notice strong evidence of the difference between MSPSO and the other three algorithms, regarding the minimum fitness of capacity as well as distance functions. The asterisk symbol in Table 3 shows where the p-value is less than alpha (α = 0.05), indicating a significant difference between a pair of algorithms in terms of the quality of solutions obtained (see Table A1).

Cost of Fdistance
The box plot presents the average cost of two objectives for the four algorithms ( Figure 5). Figure 5a,b show that AMOSA and MOABC are the algorithms with the minimum average cost for both objectives. NSGA-II is the third in optimizing both objective functions, while MSPSO has a significantly higher average cost.   Figure 6 shows that AMOSA returns the high number of solutions while MOABC returns the small number of solutions in the final Pareto front. AMOSA and NSGA-II effectively converge faster to the minimum fitness compare to MOABC and MSPSO. The large size of AMOSA's solutions is due to its capacity for archiving and clustering solutions that control diversity among the non-dominated solutions. As can be seen in Figure 6, the AMOSA, MSPSO, and NSGA-II show more evenly and smoothly distributed solutions along the Pareto front compared to MOABC.  Figure 6 shows that AMOSA returns the high number of solutions while MOABC returns the small number of solutions in the final Pareto front. AMOSA and NSGA-II effectively converge faster to the minimum fitness compare to MOABC and MSPSO. The large size of AMOSA's solutions is due to its capacity for archiving and clustering solutions that control diversity among the non-dominated solutions. As can be seen in Figure 6, the AMOSA, MSPSO, and NSGA-II show more evenly and smoothly distributed solutions along the Pareto front compared to MOABC. Figure 6 shows that AMOSA returns the high number of solutions while MOABC returns the small number of solutions in the final Pareto front. AMOSA and NSGA-II effectively converge faster to the minimum fitness compare to MOABC and MSPSO. The large size of AMOSA's solutions is due to its capacity for archiving and clustering solutions that control diversity among the non-dominated solutions. As can be seen in Figure 6, the AMOSA, MSPSO, and NSGA-II show more evenly and smoothly distributed solutions along the Pareto front compared to MOABC.

Efficiency Comparison
To evaluate the efficiency of the four algorithms in terms of convergence speed and execution time the Kruskal-Wallis test was used. The convergence speed of an algorithm is evaluated based on the fitness variation. This criterion shows how the algorithm converges toward the optimum solution through a number of iterations, while the execution time reveals how fast the algorithm is in terms of running time. The results of the efficiency criteria are presented in Table 2. The p-values in Table 2 show that there is a very significant difference in convergence speed (fitness variation rate) between the algorithms for both objective functions. This is identified by p-value = 0.000, which is less than alpha = 0.05. The post hoc tests using the Conover test were carried out for pairwise convergence comparison, and the results are presented in Table 4. From Table 4, we found that there are statistically significant differences between all algorithms when optimizing the capacity function (p < 0.05). Regarding fitness variation of distance function, only two paired comparisons of MOABC-MSPSO and AMOSA-NSGA-II did not show significant differences (p > 0.05) among six paired comparisons.   Figure 7, we notice that AMOSA outperforms the other three algorithms, with a minimum average cost for both capacity and distance functions. However, Figure 7c demonstrates that NSGA-II is the fastest algorithm compared to the three others.
This shows that the algorithm with high convergence speed is not always the one with the shortest execution time. The execution time is mostly influenced by the size of the population and the number of iterations. For AMOSA, 500 iterations have increased the computation time compared with when the algorithm runs of 100 iterations.  Figure 8 presents the convergence trends of the four algorithms, for both capacity and distance objective functions. The best fitness values of the two functions were normalized in order to facilitate their comparisons. With the progress of the algorithms, the convergence speed is reduced until the optimal solutions are attained. The mean fitness variation of MSPSO is higher compared to that of AMOSA, NSGA-II, and MOABC as shown in Figure 7a,b. Note that for AMOSA, the number of iterations displayed in Figure 8a did not attain 500 as for other algorithms. This is due to its nested loops that also iterate the cooling rate (from high temperature to lower temperature). To avoid the repetitions of solutions, we only retrieved the minimum fitness value obtained after 500 iterations of every degree of the cooling temperature (set to 100 °C).  Figure 8 presents the convergence trends of the four algorithms, for both capacity and distance objective functions. The best fitness values of the two functions were normalized in order to facilitate their comparisons. With the progress of the algorithms, the convergence speed is reduced until the optimal solutions are attained. The mean fitness variation of MSPSO is higher compared to that of AMOSA, NSGA-II, and MOABC as shown in Figure 7a,b. Note that for AMOSA, the number of iterations displayed in Figure 8a did not attain 500 as for other algorithms. This is due to its nested loops that also iterate the cooling rate (from high temperature to lower temperature). To avoid the repetitions of solutions, we only retrieved the minimum fitness value obtained after 500 iterations of every degree of the cooling temperature (set to 100 • C).
In general, the convergence speed of AMOSA and NSGA-II are higher (better) followed by MOABC and MSPSO. The reason for smoother convergence of NSGA-II is the crossover and mutation operators that influence to obtain the best survivors (offspring) for the next generation. In contrast with that, the neighborhood search strategy in MSPSO does not guarantee a better improvement of the solutions through iterations. The common challenge of this strategy is to deal with local optimums problem.
optimal solutions are attained. The mean fitness variation of MSPSO is higher compared to that of AMOSA, NSGA-II, and MOABC as shown in Figure 7a,b. Note that for AMOSA, the number of iterations displayed in Figure 8a did not attain 500 as for other algorithms. This is due to its nested loops that also iterate the cooling rate (from high temperature to lower temperature). To avoid the repetitions of solutions, we only retrieved the minimum fitness value obtained after 500 iterations of every degree of the cooling temperature (set to 100 °C).

Repeatability Test and Evaluation
A good optimization algorithm is supposed to generate similar results for different runs with the same input parameters. In this section, the repeatability and stability of each algorithm are investigated using the variance of average-normalized fitness values and the average execution time. Each algorithm was implemented five times with the same input data and their results are presented in Table 5. The four algorithms are different regarding repeatability. As shown in Table 5, NSGA-II has the lowest average execution time of 30 runs, followed by AMOSA and MOABC, and then MSPSO. AMOSA, MOABC, and MSPSO have the lowest average-normalized fitness values for both capacity and distance functions. This indicates that in terms of quality of solutions and repeatability, MOABC and AMOSA are to prefer solving evacuation problems.
The box plot in Figure 7c shows that the average execution time of NSGA-II is 363.03, which is less than half the value of MOABC, and less than a third of MSPSO. Although MOABC and MSPSO are both swarm intelligence algorithms, MOABC outperforms MSPSO. The reason for this difference can be related to the time-consuming neighborhood search process by the particles. The main part of the computation is spent on the calculation of neighborhood topology and comparison of local best and global best fitness values of particles. MOABC, on the other hand, performs a quick exploration of scouts and the share of information between employee and onlooker bees. Figure 9 presents the maps of the distribution of the population to shelters as outputs of each algorithm after optimizing the two defined functions. Three solutions are selected from the Pareto front of each algorithm (see Figure 7) by giving higher weight to either minimum objective 1 (Equation (1)) or minimum objective 2 (Equation (2)), or considering the same weight for both objectives. All solutions are optimum and there is a trade-off between them. Meanwhile, decision-makers can select an optimum solution, based on his/her preferences (Figure 9a-c). The lines with different colors represent the allocation of the population from each building block to shelters. The illustrated maps cannot be regarded as an optimal solution for evacuation planning in the city of Kigali. However, decision-makers and planners can use them as input to facilitate the procedure of planning a better distribution of population among the shelters/safe areas. This can be observed in Figure 9 on graphs 1 and 2, where the lines connecting shelters and building blocks look less crowded than graphs 3 and 4.

Conclusion
The objective of this study was to compare the performance of four multi-objective optimization algorithms (AMOSA, MOABC, MSPSO, NSGA-II respectively) for a given spatial problem, namely evacuation planning. In our study, the evacuation problem was aiming to minimize the accumulated distance from high-risk zones to shelters and to minimize the total capacity overload cost of shelters. The higher the minimum fitness values of both capacity and distance are, the better are the obtained alternatives for assigning people to appropriate shelters.
In terms of algorithm performance, all algorithms generated the optimization in a consistent way, and no results were obtained that could suggest that some of them were trapped in a local minimum. By evaluating the convergence speed of the fitness variation of the four algorithms (see Figure 8), we found that AMOSA and NSGA-II followed by MOABC converge faster and smoother towards the final optimal solutions. This justifies not only the competence of NSGA-II, which has been used in the literature to a larger extent than the other algorithms [60]. However, the competence of AMOSA and MOABC shows the capacity of solving multi-objective optimization problems including evacuation problems.
The presented metaheuristic methods and others of its type are not meant to find a 'single perfect solution' but a set of 'good enough' solutions in an efficient way, and therefore, it is possible that a more optimal solution can be achieved by using alternative methods. Decision-makers must be aware of this aspect, in order to properly assess the benefits and limitations of these techniques.
A suggestion for future work, as an alternative approach dealing with this type of spatial multiobjective optimization problems, is to modify the classical algorithms to better fit the problem in hand. For example, based on the results obtained by MOABC and the comparison made to other

Conclusions
The objective of this study was to compare the performance of four multi-objective optimization algorithms (AMOSA, MOABC, MSPSO, NSGA-II respectively) for a given spatial problem, namely evacuation planning. In our study, the evacuation problem was aiming to minimize the accumulated distance from high-risk zones to shelters and to minimize the total capacity overload cost of shelters. The higher the minimum fitness values of both capacity and distance are, the better are the obtained alternatives for assigning people to appropriate shelters.
In terms of algorithm performance, all algorithms generated the optimization in a consistent way, and no results were obtained that could suggest that some of them were trapped in a local minimum. By evaluating the convergence speed of the fitness variation of the four algorithms (see Figure 8), we found that AMOSA and NSGA-II followed by MOABC converge faster and smoother towards the final optimal solutions. This justifies not only the competence of NSGA-II, which has been used in the literature to a larger extent than the other algorithms [60]. However, the competence of AMOSA and MOABC shows the capacity of solving multi-objective optimization problems including evacuation problems.
The presented metaheuristic methods and others of its type are not meant to find a 'single perfect solution' but a set of 'good enough' solutions in an efficient way, and therefore, it is possible that a more optimal solution can be achieved by using alternative methods. Decision-makers must be aware of this aspect, in order to properly assess the benefits and limitations of these techniques.
A suggestion for future work, as an alternative approach dealing with this type of spatial multi-objective optimization problems, is to modify the classical algorithms to better fit the problem in hand. For example, based on the results obtained by MOABC and the comparison made to other algorithms, MOABC could be an interesting algorithm to modify in order to solve complex problems such as evacuation planning. It is also important to consider the use of other methods, such as recoverable robustness, to solve evacuation planning. Iris and Lam [61] proposed a recoverable robust optimization approach for the weekly berth and quay crane planning problem. The results proved the strength of the proposed model for solving a spatial problem.