Enhancing evacuation response to extreme weather disasters using public transportation systems: a novel simheuristic approach


 In recent years, there have been an increasing number of extreme weather events that have had major impacts on the built environment and particularly on people living in urban areas. As the frequency and intensity of such events are predicted to increase in the future, innovative response strategies to cope with potential emergency conditions, particularly evacuation planning and management, are becoming more important. Although mass transit evacuation of populations at risk is recognized to play a potentially important role in reducing injury and mortality rates, there is relatively little research in this area. In answering the need for more research in this increasingly important and relatively new field of research, this study proposes a hybrid simulation–optimization approach to maximize the number of evacuees moved from disaster-affected zones to safe locations. In order to improve the efficiency of the proposed optimization approach, a novel multipopulation differential evolution approach based on an opposition-based learning concept is developed. The results indicate that even for large populations the proposed approach can produce high-quality options for decision makers in reasonable computational times. The proposed approach enables emergency decision makers to apply the procedure in practice to find the best strategies for evacuation, even when the time for decision making is severely limited.

During the past five decades, the world's population has doubled and the proportion of people living in urban areas has increased from 36% to 55%. It is estimated by the United Nations that this rate will reach 66% by 2050 (Mandache, 2013). In many countries, a large number of cities have been built or developed in disaster-prone regions (Handayani, Fisher, Rudiarto, Sih Setyono, & Foley, 2019), and in recent decades natural disasters from EWEs such as Hurricanes Katrina, Rita, and Wilma have turned global attention to effective urban evacuation planning and management to reduce potential injury and mortality rates (Yin, Cordahi, Roden, & Wolshon, 2018). The United Nations International Strategy for Disaster Reduction (UNISDR) has defined "Evacuation" as "Moving people and assets temporarily to safer places before, during or after the occurrence of a hazardous event in order to protect them." (UNISDR, 2009). Evacuation management is one of the most important parts of disaster operations management, and is widely recognized as playing a prominent role in reducing mortality and injury rates during such events (Lindell, Kang, & Prater, 2011;Na & Banerjee, 2019;Shahabi & Wilson, 2018;Wei & Lindell, 2017), one of the key targets of the Sendai Framework for Disaster Risk Reduction (Meyer, Mitchell, Purdum, Breen, & Iles, 2018).
Public transportation is one critical element of evacuation management. However, despite its criticality as a resource to evacuation management, it has received relatively little attention from researchers. This is somewhat surprising given rapid changes in mass transportation systems and technologies in many urban areas. For example, during Hurricane Katrina in 2005, the evacuation was only based on private vehicles, and did not adequately consider those without access to personal transportation (Bayram, 2016;Litman, 2006). It is known that this was one of the main reasons for the failure of the evacuation process during Hurricane Katrina (Bayram, 2016;Litman, 2006) where more than 1800 people died. A closer look at the mortality data reveals that a large proportion of the fatalities, near to 70%, were older people who did not have access to private vehicles (Bayram, 2016;Litman, 2006), highlighting the critical need to incorporate public transportation systems in evacuation planning to provide necessary assistance to the aging, disabled, immobilized population or those who do not have access to a personal vehicle. This issue has been largely unaddressed in official evacuation plans (Hess & Gotham, 2007;Urbina & Wolshon, 2003) and there is critical need for both researchers and decision makers to pay more attention to public transport (such as bus transit systems) as an evacuation resource (Hu, Gao, Yu, Liu, & Li, 2016).
Evacuation planning problems are known for their high complexity and the need to produce optimum solutions in short times. In solving many complex optimization problems such as evacuation strategies, metaheuristics have proven successful (Azadeh, Seif, Sheikhalishahi, & Yazdani, 2016;Mostafa Bozorgi & Yazdani, 2019;Yazdani, Aleti, Khalili, & Jolai, 2017;Yazdani, Babagolzadeh, Kazemitash, & Saberi, 2019;Yazdani, Jolai, Taleghani, & Yazdani, 2018;Yazdani, Khalili, Babagolzadeh, & Jolai, 2017;Yazdani, Khalili, & Jolai, 2016). However, one of the main limitations of metaheuristic methods is that parameters are deterministic and conflict with real-world conditions such as those presented by EWE evacuations (Juan, Faulin, Grasman, Rabe, & Figueira, 2015). To address this problem, simheuristics is an innovative and potentially efficient method that integrates simulation into metaheuristicdriven frameworks to take account of uncertainties presented by real-world evacuation problems. However, this method has never been used in an evacuation context. To address this gap in the EWE disaster management and evacuation research literature and find a fast method in order to reach the most appropriate decision in reasonable time, this paper proposes a new approach for evacuation planning by using a novel framework that was proposed by Juan et al. (2015) to take advantages of both simulation and metaheuristic algorithm.
In addition, although there are many studies that utilize a simheuristic approach, the majority of prior research has applied single-solution-based metaheuristics (S-metaheuristics) in their simheuristic framework (e.g. Gruler, Panadero, de Armas, Moreno Pérez, & Juan, 2020;Panadero, Doering, Kizys, Juan, & Fito, 2018). S-metaheuristics apply generation and replacement procedure to a single solution (Gogna & Tayal, 2013), but most of them suffer from fast convergence to local optima, and also they are not able to efficiently explore the search space (Gogna & Tayal, 2013). This gap can be filled by applying population-based metaheuristic algorithms like differential evolution (DE) into simheuristic methods. DE first proposed by Storn and Price (1997) is a potentially powerful population-based metaheuristic algorithm to solve a wide range of complex problems.
Although DE has been known as a powerful populationbased metaheuristic algorithm to solve a wide range of complex problems, it suffers from slow convergence. Therefore, in this study, the basic version of DE has been combined with the opposition-based learning (OBL) concept to increase the exploration of the search space. Using OBL in the proposed method can significantly increase the chance of finding better solutions because current candidate solutions and their corresponding opposite are evaluated simultaneously, and the algorithm continues with more desirable ones. In addition, it has been proven that appropriate mutation strategies and parameters are not constant and can change from one optimization problem to another one (Qin, Huang, & Suganthan, 2009). Furthermore, numerous studies reveal that using different mutation strategies may improve exploitation and exploration capabilities of a DE algorithm (Das & Suganthan, 2011). Therefore, different mutation strategies are used in the proposed DE. Thus, the proposed algorithm uses a set of mutation strategies randomly to maintain its exploitation and exploration capabilities. In addition, previous studies on multipopulation metaheuristic algorithms show that using multipopulation can maintain population diversity (Ali, Awad, & Suganthan, 2015). Therefore, the proposed DE is designed based on a multipopulation strategy with the objective that the outcomes of this research can help decision makers and disaster planners to make better decisions for evacuation in EWE disaster scenarios.
The rest of this paper is structured as follows: Section 2 investigates related works. The problem statement is discussed in Section 3. Section 4 introduces the proposed approach to solving the considered problem. The results are presented in Section 5. Concluding remarks and some future research directions are given and discussed in the last section.

Literature Review
A review of the related literature is presented next in two streams: bus evacuation problems and simheuristic applications. Bish (2011) was among the first researchers who considered bus routing in an evacuation problem context, presenting two mathematical programming formulations for the problem and developing a heuristic algorithm to minimize entire network evacuation time. In the model, buses and shelters have limited capacity, and each pickup location can be visited by several buses to satisfy the evacuation demand. For a variant of this problem, Goerigk, Grün, and Heßler (2013) suggested several methods to find lower and upper bounds and then proposed a branch-and-bound framework based on them for the bus evacuation problem. Sayyady and Eksioglu (2010) also studied the use  of public transport in an emergency evacuation. They proposed a mixed-integer linear programming model to find optimal evacuation routes in order to minimize the total evacuation time and the number of casualties, simultaneously. Due to the high complexity of the problem, a Tabu search algorithm was developed to find evacuation routes for transit vehicles. Dikas and Minis (2016) considered a variant of bus evacuation problem where the buses must return to the single depot at the end of the route. They presented a heuristic and a mixed-integer linear formulation for the problem.  proposed a mixed-integer programming formulation for a bus evacuation problem under uncertainty, where the objective was minimizing the network clearance time. In their model, uncertainty could occur in some parameters such as travel time or capacity in addition to demand. Due to the high complexity of the problem, they presented a Tabu search algorithm in order to find an acceptable quality solution within short computation time. Goerigk, Grün, and Heßler (2014) formulated the integrated bus evacuation problem to determine both the shelter and collection points for evacuating a region using buses. To solve the problem, they proposed a branch-and-price method, where the pricing problem was the shortest path in a round-expanded path, and compared its efficiency using a commercial IP solver. Goerigk, Deghdak, and T'Kindt (2015) considered a different robust version of the problem, which allows buses to change their route as the uncertainty on the number of evacuees diminishes. For this, they proposed a scenario-generation algorithm since the model was too large to be solved. Kulshrestha, Lou, and Yin (2014) presented a mixed-integer linear program to determine the optimal locations of transit pickup points and trip allocations for buses in transit-based evacuation planning under demand uncertainty. The proposed model was solved by a cutting plane scheme. Zheng (2014) studied an evacuation problem based on public transit in which buses run continuously, based upon the spatial and temporal information of evacuee needs. Zheng developed a Lagrangian relaxation-based solution algorithm to minimize the exposed casualty time. Qazi, Nara, Okubo, and Kubota (2017) investigated the effects of introducing demand variations and evacuation route flexibility in a bus-based evacuation problem. They presented a case study for the evacuation of elderly people living in a small town, Kawajima, surrounded by two rivers on its two sides and exposed to flood hazards. Pereira and Bish (2015) studied a bus evacuation problem where evacuees arrive at pickup locations at a constant rate. Minimizing the total waiting time at pickup locations was considered the objective function. Hu et al. (2016) introduced a nonlinear integer programming model and a genetic algorithm to plan a bus bridging evacuation to best plan the bus bridging evacuation in both upstream and downstream directions along a rail transit line during a disruption. Shahabi and Wilson (2018) presented an iterative algorithm for a large-scale evacuation routing, considering unexpected roads' capacity change during evacuation. They used heuristics and dynamic data structures to improve evacuation routing running time. Swamy, Kang, Batta, and Chung (2017) provided a mass evacuation strategy using public transportation before the strike of a hurricane. They assumed that shelter locations, evacuation zones, and the time of strike of the hurricane are predetermined. They developed a simulation tool to model the dispatching of the given number of buses, stochastic arrival of evacuees, queuing effects at the pickup locations, and the transportation of evacuees to the safe regions. Vitali, Riff, and Montero (2017) presented an algorithm based on a greedy randomized adaptive search procedure to solve the bus evacuation problem. Most recently, Downloaded from https://academic.oup.com/jcde/advance-article-abstract/doi/10.1093/jcde/qwaa017/5813879 by Joongbu University user on 03 April 2020
Lakshay and Bolia (2019) developed mathematical models for mass evacuation problems. The model computed the minimum number of buses that are required to carry out the evacuation for any given scenario and determine the optimal trip sequence of each available bus. Despite all of the above-mentioned studies, uncertainty in parameters has only been considered in a limited way and too little attention has been paid to disruption in road networks. A major problem with bus evacuation is associated with computational complexity and to find the best optimum solution in a reasonable time is challenging. Finally, we can find no research that uses an efficient optimization approach or that proposes a hybrid optimization-simulation approach for bus evacuation.

Simheuristic methods for complex problems under uncertainty
The combination of metaheuristic algorithms with other methodologies is gaining more attention as a method to tackle a wide variety of disciplines, including complex combinatorial optimization problems (Ferreira, 2013). A simulation-optimization algorithm is one such method that has widely been used to solve complex large-scale real-world problems under uncertainty conditions (Rabbani, Heidari, & Yazdanparast, 2019). The general structure of an integrated simulation-optimization approach and their interaction to find near-optimal solutions in search space is presented in Fig. 1. Simheuristics are approaches that benefit both simulation and (meta-)heuristic approaches. In fact, simheuristics incorporate simulation into a (meta-)heuristic-driven framework (Juan et al., 2015). Simheuristics can be applied in solving many combinatorial optimization problems with stochastic components (Rabbani et al., 2019). Recently, Juan et al. (2015) have reviewed and classified the studies on the simheuristic methods proposing a simple but efficient framework where simulation is applied to a subset of the obtained solutions. Over the past few years, Gonzalez-Neira, Ferone, Hatami, and Juan (2017) proposed an approach based on the integration of biased randomization and simulation techniques inside a metaheuristic framework for a distributed assembly permutation flow shop problem with stochastic processing times. Hatami, Calvet, Fernández-Viagas, Framiñán, and Juan (2018) investigated setup starting times in the stochastic parallel flow shop problem. They proposed a simheuristic algorithm based on iter-ated local search metaheuristic with Monte Carlo simulation to minimize makespan-related criteria in a stochastic parallel flow shop problem. Panadero et al. (2018) proposed a simheuristic approach to maximize the net present value in a project portfolio selection problem under uncertainty and rich conditions. Their algorithm integrated Monte Carlo simulation inside a variable neighborhood search framework. Furthermore, their algorithm used other components inspired in simulated annealing and biased randomization techniques. Guimarans, Dominguez, Panadero, and Juan (2018) investigated a 2D vehicle routing problem with stochastic travel times. A hybrid simheuristic algorithm was proposed based on Monte Carlo simulation, an iterated local search framework, and biased randomized routing and packing heuristics. Gruler, Panadero, de Armas, Moreno Pérez, and Juan (2018) presented a simheuristic approach that integrates Monte Carlo simulation within a variable neighborhood search framework to solve the multiperiod inventory rout-Downloaded from https://academic.oup.com/jcde/advance-article-abstract/doi/10.1093/jcde/qwaa017/5813879 by Joongbu University user on 03 April 2020   ing problem with stochastic customer demands. The proposed approach allowed to consider the inventory changes between periods generated by the realization of the random demands in each period, which have an impact on the quantities to be delivered in the next period and therefore on the associated routing plans. In another study, Gruler et al. (2018) presented a variable neighborhood search metaheuristic hybridized with simulation to solve the inventory routing problem with stochastic demands. They claimed that the proposed approach could solve largesized problems for the single-period inventory routing problem with stochastic demands and stock-outs in very short computing times.  (2019) proposed a multi-objective greedy randomized adaptive search procedure metaheuristic procedure coupled with a Monte Carlo simulation and a Pareto archived evolution strategy algorithm to obtain a complete set of Paretooptimal solutions for a bi-objective stochastic permutation flow shop scheduling problem.
Most studies in the field of hybrid optimization-simulation methods have considered that each time the simulation component is performed to evaluate the objective function; however, this procedure is very time consuming and imposes additional efforts to the method to find the suitable solution. In addition, most of the previous studies used single-solution-based metaheuristics such as variable neighborhood search and Tabu search in their framework, which are very simple and fast, but most of them cannot guarantee escaping from local optima and also fast convergence to global optima. Furthermore, almost all the population-based metaheuristic algorithms used in the previous studies are basic methods and are not efficient enough to find good solutions in reasonable time.

Problem Statement
This section presents a detailed description of bus evacuation planning under disruption in a road network during an EWE. In case of an approaching EWE, the population is exposed to significant risks; thus, authorities broadcast compulsory evacuation warning messages. These messages advise assembly at some organized pickup locations. Then, some limited number of buses are used to move evacuees to safe shelters via available routes. In this case, some parts of road network infrastructures are likely to be disrupted.
The first step in bus evacuation planning is to know all the pickup locations in the endangered area and the number of evac- uees in each pickup location. Then, based on a number of demands and situations, safe shelters are identified. In addition, available routes between pickup locations and selected shelters should be calculated. The next step is confirming the number of available buses at different depots. After selecting an appropriate objective, the optimal evacuation plan can be generated using some advanced optimization method. Figure 2 presents steps of evacuation planning based on buses.
The focus of this research is on the fifth step, "allocation and planning of transport resources," where transport resources in this study are buses. In this problem, first, people are guided by authorities to assemble at the pickup locations i = {1, 2, 3, . . . , I } that are safe and accessible, and then are moved to a set of capacitated safe shelters j = {1, 2, 3, . . . , J } using b = {1, 2, 3, . . . , B} buses. In this problem, there are I pickup locations in a disaster area, and J candidate places for shelters. B is the total number of buses used in the evacuation. Buses are assumed to carry a fixed number of evacuees BC b at once. There might be k different routes between pickup location i and the shelter j. Moreover, as shown in Fig. 3, each route may contain several segments and each of them will have a specific risk profile. Furthermore, evacuation is subjected to disruption of road network infrastructures, and they may not remain accessible and available during an evacuation horizon. The objective is to evacuate the disaster zone in such a way as to maximize the number of evacuees that are moved from a disaster zone to safe shelters in a predefined time window (Fig. 3). In addition, due to the high number of evacuees, it is assumed that each bus travels several times between pickup locations and the shelters within evacuation time window constraint T. Figure 4 illustrates a small sample of the problem. There are two buses and the capacity of each of them is 50. There are two pickup locations and two safe shelters. There are 150 and 100 evacuees in pickup locations 1 and 2, respectively. Also, the capacities of shelters 1 and 2 are 150. There are two routes between each pair of pickup location and a safe shelter. The traveling times between pickup locations and safe shelters for routes 1 and 2 are given in the figure.
Route 2 between pickup point 1 and shelter 1 and route 2 between pickup point 1 and shelter 2 are not available after time 640 due to disruption of the route. Table 1 presents a random feasible solution to the presented example. For example, the first bus goes from yard to pickup location 1 via route 2 and then goes to shelter 2 via route 1, and so on. If it is assumed that window time is 650, then only 200 evacuees can be moved within the evacuation window time using this strategy.

Solution methodology
Simheuristics have shown good capability in solving a wide range of problems with uncertainty in different fields (Gonzalez-Neira et al., 2017;Guimarans et al., 2018;Hatami et al., 2018;Juan et al., 2015;Rabbani et al., 2019). The proposed framework by Juan et al. (2015) has been modified for this research. In this study, the simheuristic approach assumes that a good solution in the deterministic condition is also likely to be an acceptable solution for its corresponding stochastic condition. However, it does not mean the best solutions for the deterministic and stochastic version are the same (Juan et al., 2015). The proposed integrated framework assumes that an efficient optimization method already exists for the problem under the deterministic condition. Therefore, a novel DE is developed for this problem. The general framework of the proposed simheuristic method is presented in Fig. 5. An initial population is generated in search space, where some parameters are stochastic. In the next step, stochastic values are replaced by their expected values, which are deterministic. Then, the optimization method, which in this study is a novel DE algorithm, is run to find a set of good solutions for the deterministic version (details of this algorithm are discussed in the next sections). After passing a predefined number of iterations, the quality of each of solutions that have been obtained by the proposed DE is evaluated using a simulation component. In this simulation phase, deterministic parameters are replaced by their stochastic values. The obtained solutions through simulation are ranked and then sent to metaheuristic. Once the stopping condition is satisfied, more accurate evaluations can be obtained for obtained solutions by employing intensive simulation. Finally, the solutions are ranked, and the best solution is obtained.
As mentioned earlier, the basic version of DE, in this study, is combined with the OBL concept. In addition, to maintain the population diversity, the multipopulation strategy is applied in the proposed DE algorithm. Moreover, different mutation strategies are used to improve the exploitation and exploration capabilities of the proposed method. Therefore, in the next section, the basic DE algorithm and the most frequently used mutation strategies are discussed.

Basic DE algorithm
In general, DE starts with generating a predefined number (N) of solutions within the search space. At generation G, each individual is represented by a D-dimensional real-valued vector  There are N individuals in the population of DE. Generally, DE starts with a random population, where parameters of the i th individual of the initial population are generated as follows: where D denotes the number of decision variables and rand returns a uniformly distributed random number in [0, 1]. In addition, x low j and x up j are the lower and upper bounds of solutions in the j-dimensional search space, respectively.
Following initialization, DE employs the mutation procedure to create mutant vectors (Huang, Zhang, Song, Zhang, & Shi, 2019). The most frequently used mutation operators are "rand/1: "best/1: "best/2:" "current-to-best/1: where r i 1 = r i 2 = r i 3 = r i 4 = r i 5 are random integer numbers selected from (1, 2, . . . , N). The selected numbers are different from the index i. X G best is the best solution at the current generation. F is a control parameter for scaling the difference vector.
After this stage, the crossover operator is used to mix the target vector and its corresponding mutant vector to yield a trial vector (U G i ): where Cr is crossover rate and is between 0 and 1, and rand is a uniformly distributed random integer between 0 and 1. In addition, j rand is a uniformly distributed random integer in [1, D].
In order to decide whether the target vector (X G i ) or its corresponding trial vector (U G i ) is retained for the next generation, a selection operation is employed between these solutions. For problems where the aim is minimizing, the selection operation is performed as follows: where f (X G i ) and f (U G i ) are the objective functions of the target vector and its corresponding trial vector, respectively.
. , x d are real numbers and x i ∈ [a i , b i ], i = 1, 2, . . . , d, its reflected extended opposition point X reo i (x reo 1 ,x reo 2 , . . . ,x reo d ) is defined as follows:

The proposed DE
In this section, a multipopulation opposition-based DE technique is proposed. The initial individuals are randomly generated in the search space. The structure shown in Fig. 6 is applied to represent solutions. The length of each solution is twice the maximum number of trips a bus can travel between pickup locations and shelters within the maximum available time window. Thus, the maximum available time window is divided by the shortest travel time of the available routes between pickup locations and shelters and the obtained number is multiplied by 2.
As an example, there are two pickup locations and two safe shelters, and between each pair of pickup location and shelter, Downloaded from https://academic.oup.com/jcde/advance-article-abstract/doi/10.1093/jcde/qwaa017/5813879 by Joongbu University user on 03 April 2020 there are three routes. The solution S = (0. 21, 0.71, 0.91, 0.82, 0.45, 0.42) indicates that bus 1 goes from yard to the first pickup location, which here is 2 (because we have two pickup locations, 1 is divided by 2; since 0.71 is larger than 0.5, it means that it is pickup location number 2), via route 1 (because we have three routes, 1 is divided by 3; since 0.21 is smaller than 0.33, it means that it is route number 1), then from pickup point 2 to shelter 2 via route 3, and then from shelter 2 to pickup point 1 via route 2.
Afterward, generated individuals are divided into subpopulations. For each individual, one of the mutation strategies in equations (2)-(6) is selected randomly to generate a mutant vector. Each subgroup is evolved separately for a predefined number of iterations. Then, some solutions are randomly selected to migrate between subgroups. Before inserting selected individuals into new groups, one of the opposition-based strategies is randomly selected and applied to the individual. The process of searching ends once the predefined criterion is satisfied. Pseudocode of the proposed algorithm is presented in Algorithm 1.
Algorithm 1. Pseudocode of the proposed DE algorithm.

1:
The proposed algorithm 2: for i ←1 to NG do ᭠ NG is the number subpopulation.

3:
for j ←1 to NP do ᭠ N P is the number of individuals in each subpopulation.

4:
Generate a solution in the search space (P i j ).

6:
Add solution to the subpopulation i. 7: end for 8: end for 9: δ ← 0 10: while !termination criterion do 11: for i ←1 to NG do 12: for j ←1 to NP do 13: Select P Best and four solutions (P i1 , P i2 , P i3 , and P i4 ) randomly from current subpopulation (i 1 = i 2 = i 3 = i 4 ) 14: R ← rand i(i max ) ᭠ Returns an integer between 1 and i max , which here is 5.

44:
for j ←1 to NP do ᭠ N P is the number of individuals in each subpopulation.

45:
if rand(0, 1) < R OBL then ᭠ Returns random number between 0 and 1 and R OBL is probability of doing OBL.

46:
E ← rand i(i max ) ᭠ Returns an integer between 1 and i max , which here is 5.

Results and Discussion
In this section, computational tests performed and results obtained are presented and discussed in order to verify the performance of the proposed novel DE algorithm against traditional DE algorithms in deterministic conditions. We also compare the performance of the proposed simheuristic algorithm against a traditional simulation-optimization method when simulation component is run in each iteration. This section therefore presents the random testing instances generated. Then, parameters of algorithms are tuned and described. Finally, the obtained results are discussed. It should be noted that the algorithms described in the previous sections have been implemented in MATLAB R2017b and all experiments were run on a PC with an Intel QuadCore i7-7700 CPU at 3.60 GHz and 16.0 GB RAM.

Data generation
Since the variant of the bus evacuation problem discussed in this study has not been studied in previous research works, in order to compare the proposed approach with others we had to generate some instances. Therefore, 20 instances are generated with different sizes. Table 2 provides values for the number of pickup points, the number of shelters, and the number of buses. There are three routes between each pair of pickup location and shelter. Stochastic travel time between nodes is generated using a truncated normal distribution, where means are randomly generated between 10 and 200 for route 1, between 20 and 400 for route 2, and between 40 and 800 for route 3. The number of evacuees at each pickup point is generated between 250 and 500. The capacity of each bus is 50. Each route can be disrupted with specific probability; thus, two different probabilities have been considered in the generated instances, 0.1 and 0.3.
Each instance is run 10 times (Hatami et al., 2018). The best results are recorded. In the proposed simheuristic, 100 and 20000 runs are employed during the solution approach and at the end, respectively. In addition, stop condition is set at 60 s, because in the evacuation problem a short time is available to make a good decision; thus, to simulate this condition, the stopping condition is fixed for all instances.

Parameter settings
Choosing the appropriate parameters can help the algorithm to improve its convergence speed as well as the quality of the obtained solutions. The proposed algorithm has several parameters, which are population size (NP), crossover factor (Cr), number of groups, and scaling factor (F). To find the values of parameters, several instances were randomly generated. Some parameters from similar algorithms in the literature were found and Downloaded from https://academic.oup.com/jcde/advance-article-abstract/doi/10.1093/jcde/qwaa017/5813879 by Joongbu University user on 03 April 2020 the algorithm run with their different combinations on the testing set. Table 3 presents the most appropriate parameters for the proposed problem. In addition, the parameters of other DE algorithms are set at a population size (NP) of 100, a scaling factor (F) of 0.5, and a crossover factor (Cr) of 0.9, based on Cai et al. (2017).

Experimental results
In this section, a performance measure is required to represent the performance of algorithms. Therefore, the relative deviation index (RDI) (Pan, Ruiz, & Alfaro-Fernández, 2017) is used to normalize the values obtained. Smaller values of RPD are desirable. RDI is calculated as where Alg sol is the objective value obtained from each algorithm, and Best sol and Worst sol are the best and the worst solutions obtained among all methods, respectively. First, to show the performance of the proposed DE against previous traditional DEs, instances 1-10 are selected and algorithms are run in deterministic condition. The experimental results of different sets are shown in Table 4. In this table, the first column shows the instance number, and the other columns indicate obtained RDI for each algorithm. Furthermore, the last row presents the average values in each column. The best RDI in each row is given in bold.
As it is clear in Table 4, the proposed DE with the average RDI of 12.1% is far superior to other DEs. Furthermore, findings verify that using the OBL concept in the proposed algorithm can significantly improve the performance of the proposed algorithm. Figure 7 presents the performance of the proposed algorithm with and without using the OBL concept in different instances when RDI values of different instances are calculated based on these two methods. As shown in Fig. 7, in the presence of OBL, the proposed algorithm has better outputs in all sizes.
The Wilcoxon signed-rank test as a nonparametric statistical analysis method is used to further verify the performance of the proposed algorithm. The Wilcoxon signed-rank test is used when we compare two related samples or matched samples. Therefore, it can be used to test the significance of the algorithms (Biswas & Biswas, 2017). Wilcoxon's test is conducted with a 1% significance level. The obtained result using the Wilcoxon signed-rank test is presented in Table 5. Based on the obtained results, it is clear that the proposed method outperforms other DEs with the level of significance θ = 0.01 since the P -value of each of the pairs is below 0.01. On the other hand, the z-value exceeds the critical values. This indicates that the comparison is statistically significant.
In the next step, the proposed simheuristic is compared against other simheuristics when their search engines are basic DEs. This comparison is conducted to show the role of the proposed DE in increasing the efficiency of the proposed simheuristic. In both groups of instances, the performances of the proposed algorithm are the same. Table 6 summarizes the results of the experiments and presents the RDI for each algorithm and each instance. The best performing algorithm is the proposed simheuristic with average RDI of 40.4% and 44.7% for two different groups of instances.To check for statistical significance of the results and to confirm which is the best algorithm, we performed the Wilcoxon signed-rank test. According to Table 7, the performance of the simheuristic algorithm is significantly better than those for all problem instances.    In the majority of studies, simulation is run in each iteration, which is a very time-consuming procedure. Therefore, in order to investigate the effects of the proposed simheuristic, the results obtained by the proposed simheuristic are compared with the results when simulation is run in each iteration of the metaheuristic. As it is clear in Table 8, using this simheuristic framework can help the algorithm to find better solutions in different instances. Table 9 summarizes the results of the Wilcoxon signed-rank test. The results reveal that at significance level of 0.01 the proposed approach is statistically superior to the simulationoptimization method when the simulation is run in each iteration. Hence, the proposed approach can help to reduce the computational efforts.

Conclusions
The aim of this paper was to propose a hybrid simulationoptimization approach to maximize the number of evacuees moved from disaster-affected zones to safe locations. Addressing a gap in both EWE disaster management research and operation research to solve complex evacuation planning using public transport systems, this study developed an innovative solution to optimize the available buses to maximize the number of evacuees in the affected area moved via the available emergency evacuation routes to the available shelters during a specific period and considering the evacuation route disruption risk. Many cities across the world are prone to EWEs such as major floods. In case of approaching some extreme natural hazards, the population at risk is advised to evacuate, but some evacuees do not have access to personal transportation means. Therefore, buses should be provided for these evacuees in order to move them to safe shelters. In many situations, there are a limited number of buses; therefore, they need to be run continuously during the planning horizon. In addition, evacuation planning is associated with considerable uncertainties such as travel time of buses through the road network. Furthermore, the road network can be affected by disasters, and it may affect travel time. Also, evacuees with personal vehicle tend to use some short routes, which may further cause route congestion. Considering these uncertainties in evacuation planning, finding an appropriate solution is challenging and complex.
To address the challenges associated with the emergency bus evacuation, this study focused on route planning of the available buses to maximize the number of evacuees in the affected area moved via the available emergency evacuation routes to the available shelters during a specific period and considering the evacuation route disruption risk. This study presented an efficient method based on combining a metaheuristic with a simulation method. To increase the chance of finding better solu-tions using the proposed method, a novel DE was developed for the search engine of the approach. In the proposed DE, some OBL strategies were used to evaluate current candidate solutions and their corresponding opposite at the same time. In addition, a random strategy was applied to use different mutation strategies in order to improve the exploitation and exploration capabilities of the proposed DE algorithm. Furthermore, to maintain population diversity, the proposed DE was designed based on the multipopulation strategy. In order to assess the performance of the proposed solution approaches, a set of numerical experiments was conducted to evaluate the proposed approach in terms of solution quality. First, the proposed novel DE algorithm was found to be more promising than traditional DE algorithms in deterministic condition. Furthermore, the proposed hybrid approach demonstrated a good performance when the proposed DE algorithm is used as the search engine. Finally, the performance of the proposed simheuristic algorithm was compared against a traditional simulation-optimization method when the simulation component is run in each iteration. The results indicate that the proposed method can achieve more desirable solutions. The proposed hybrid simulation-optimization approach can be used as efficient, practical tools by emergency management authorities in improving the utilization of emergency evacuation routes and emergency shelters and reducing the travel time of evacuees during an emergency evacuation. Moreover, the proposed approach is expected to improve the overall effectiveness of the bus emergency evacuation process and ensure the safety of evacuees, including aging and disabled people.
One of the future research directions of this work is considering demand uncertainty in the problem. In addition, the exact method could be used to obtain optimal solutions or lower bounds in small-scale scenarios. Considering other objectives for the problem can be an interesting future research direction.