{'fraction_non_alphanumeric': 0.03685814957045014, 'fraction_numerical': 0.00991725721525001, 'mean_word_length': 4.341404102347219, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Scalable Data-Driven Technique for Joint Evacuation Routing and Scheduling Problems

Evacuation planning is a crucial part of disaster management where the goal is to relocate people to safety and minimize casualties. Every evacuation plan has two essential components: routing and scheduling. However, joint optimization of these two components with objectives such as minimizing average evacuation time or evacuation completion time, is a computationally hard problem. To approach it, we present MIP-LNS, a scalable optimization method that combines heuristic search with mathematical optimization and can optimize a variety of objective functions. We use real-world road network and population data from Harris County in Houston, Texas, and apply MIP-LNS to find evacuation routes and schedule for the area. We show that, within a given time limit, our proposed method finds better solutions than existing methods in terms of average evacuation time, evacuation completion time and optimality guarantee of the solutions. We perform agent-based simulations of evacuation in our study area to demonstrate the efficacy and robustness of our solution. We show that our prescribed evacuation plan remains effective even if the evacuees deviate from the suggested schedule upto a certain extent. We also examine how evacuation plans are affected by road failures. Our results show that MIP-LNS can use information regarding estimated deadline of roads to come up with better evacuation plans in terms evacuating more people successfully and conveniently.


I. INTRODUCTION
Evacuation plans are essential to ensure the safety of people living in areas that are prone to potential disasters such as floods, hurricanes, tsunamis and wildfires. Large-scale evacuations have been carried out during the past hurricane seasons in Florida, Texas, Louisiana, and Mississippi. Examples of hurricanes when such evacuations were carried out include, Katrina & Rita (2005), Ike & Gustav (2008), Irma & Harvey (2017), Laura (2020), and Ida (2021). The recent category four hurricane, Ida, caused a total of $75 billion in damages and 55 deaths in the United States alone [1]. We are also anticipating that the 2022 hurricane season will have abovenormal activity [2]. To give a sense of the scale of evacuations due to such hurricanes, about 2.5 million individuals were evacuated from the coastal areas of Texas [3] before the landfall of Hurricane Rita. At such scale, it is essential to have an evacuation plan to ensure that people can evacuate in a safe and orderly manner. Any such plan needs to have two essential components: (i) Evacuation Routes, which are paths that the evacuees will take to egress out of the area under danger, and (ii) Evacuation Schedule which dictates when people should leave from different regions. The goal is to find routes and schedule that optimizes a desired objective such as average evacuation time, evacuation completion time.
The focus of our paper is on jointly optimizing the routes and schedules. Informally, the idea is to find a time to schedule when individuals can begin evacuation (within a given time window) and a route that would be used to evacuate, so as to minimize the objective functions capturing the system level evacuation time (see Section III for formal definition of the problems). Jointly optimizing over the routes and schedule is significantly harder from a computational standpoint (See Section IV for hardness results). Existing methods, even those designed to find bounded sub-optimal solutions, do not scale to city or county level planning problems. Thus, finding good evacuation routes and schedule within a reasonable amount of time, for a city or county with a large population, remains an open problem. Moreover, in real-life scenarios, evacuees do not follow a prescribed evacuation schedule exactly, which may negatively impact the performance of an evacuation plan. Therefore, sensitivity of an evacuation plan to such non-compliance should be examined carefully. Furthermore, the availability of all road segments throughout the entire evacuation time period is not guaranteed. For instance, before a hurricane makes a landfall, roads in low lying areas can become flooded due to heavy rainfall. Thus, a method for designing evacuation plans should be able to use information regarding the availability of roads at different times.
Our Contributions As our first contribution, we present MIP-LNS, a scalable optimization method that can find solutions to a class of evacuation planning problems, while optimizing for a variety of objectives. MIP-LNS is designed based on the well known Large Neighborhood Search (LNS) framework. It combines the idea of heuristic search with mathematical optimization (Section III). In this paper, we focus on three objectives: minimizing average evacuation time, minimizing evacuation completion time, and minimizing the average evacuation time of 'non-outlier' evacuees (Section III). We show that all of these three optimization problems are hard to approximate within a logarithmic factor (Section IV). In the same section, we also show that even if the underlying graph is a subgraph of a grid, these problems remain NP-hard.
Second, we choose Harris county in Houston, Texas as our study area and apply MIP-LNS. The county has about 1.5 million households and spans an area of 1, 778 square miles. We have used real-world road network data from HERE maps [4] and population data generated by Adiga et al. [5] to construct a realistic problem instance. In terms of size, our road network is at least five times larger than the networks considered in earlier papers. The study area has been affected by major hurricanes in the past (e.g. Rita, Ike, Harvey, Laura). Using MIP-LNS, we calculate evacuation routes and schedule. We show that, within a given time limit, MIP-LNS finds solutions that are on average 13%, 20.7% and 58.43% better than the existing method by Hasan and Van Hentenryck et al. [6] in terms of average evacuation time, evacuation completion time and optimality guarantee of the solution, respectively (Section VI-B). Through a scalability study, we also show that for smaller problem instances, MIP-LNS can find near optimal solutions very quickly.
Third, we demonstrate the robustness of our solution by considering scenarios where evacuees do not follow their prescribed schedule exactly. Instead, they either deviate slightly from their prescribed departure time or they choose their departure time completely at random in the allotted time. We simulate these scenarios using the agent-based simulator QueST [7] and show that in the former case, the original evacuation plan still remains almost as effective as before, showing our solution is robust to perturbations (Section VI).
Finally, we examine scenarios where certain roads fail and observe how that affects the efficacy of an evacuation plan. Our experimental results indicate that a risk model, which can predict when roads are likely to become unavailable, can help to design better evacuation plans in terms of convenience and successfully evacuating everyone. We show that our algorithm can easily use information from such models, making it a versatile, efficient and effective tool for evacuation planning.

II. RELATED WORK
Researchers have approached the evacuation planning problem in different ways in the past. Hamacher and Tjandra [8] formulated it as a dynamic network flow optimization problem and introduced the idea of time expanded graphs to solve it using mathematical optimization methods. However, the computational cost of their proposed methodology was prohibitively expensive. This led to several heuristic methods [9], [10] that are capable of working with larger networks. However, these methods are designed to solve the routing problem only and they either do not consider the scheduling problem at all or propose simple schemes such as letting evacuees leave at a constant rate. On the other hand, in a series of research works, Even and Pillac et al. [11], Romanski and Van Hentenryck et al. [12], and Hasan and Van Hentenryck [6], [13] considered the joint optimization problem of routing and scheduling. They formulated the problem as Mixed Integer Programs and used decomposition techniques [14], [15] to separate the route selection and scheduling process. A review of existing works on evacuation planning and management can be found in the survey paper by Bayram [16].
We use the most recent method by Hasan and Van Hentenryck [6] as our baseline and show that MIP-LNS finds better solutions in terms of average evacuation time, evacuation completion time and optimality guarantee of the solution. In addition to the minimizing average evacuation time objective, we provide direct MIP formulations for two other objectives: minimizing the average evacuation time for 'nonoutlier' evacuees and minimizing the evacuation completion time. MIP-LNS is able to optimize all three objectives without needing any modifications. Finally, in the evaluation with agent-based simulation phase, we consider non-compliance of evacuees to the prescribed schedule. We show that the original plan remains effective even if there is a certain amount of randomness in the departure time of the evacuees.
The use of convergent evacuation routes has been explored in the literature [6], [11]- [13], where all evacuees coming to an intersection follow the same path afterwards. This is also known as confluent flow [17]. Golin et al. [18] investigated the single-sink confluent quickest flow problem where the goal is minimizing the time required to send supplies from sources to a single sink. They showed that the problem cannot be approximated in polynomial time within a logarithm approximation factor. In this paper, we prove that finding evacuation routes and schedule that minimizes average evacuation time is also hard to approximate.
Heuristic search methods and meta-heuristics are generally applied to problems that are computationally intractable. The goal is to find good solutions in a reasonable amount of time. The Large Neighborhood Search (LNS) framework [19] has been successfully applied to various hard combinatorial optimization problems in the literature [20]. Very recently, Li et al. [21] proposed MAPF-LNS, where the authors applied the LNS framework to find solutions for the Multi-Agent Path Finding Problem. The joint routing and scheduling problem considered here is significantly different than their problem.

III. PROBLEM FORMULATION
In this section, we introduce some preliminary terms that we use in our problem formulation. Then we define three different objective functions for the evacuation planning problem. This allows us to formally define the optimization problems A-DCFP, CT-DCFP, and O-DCFP. Next, we describe how we construct time expanded graphs to model the flow of evacuees over time using a sample problem instance. Finally, we present a Mixed Integer Program (MIP) that represents a class of evacuation planning problems and show how we use it to formulate the above mentioned optimization problems.
Definition III.1. A road network is a directed graph G = (N , A) where every edge e ∈ A has (i) a capacity c e , representing the number of vehicles that can enter the edge at a given time and (ii) a travel time T e representing the time it takes to traverse the edge.
Definition III.2. Given a road network, a single dynamic flow is a flow f along a single path with timestamps a v , representing the arrival time of the flow at vertex v that obeys the travel times. In other words, a v − a u ≥ T uv . A valid dynamic flow is a collection of single dynamic flows where no edge at any point in time exceeds its edge capacity.
Definition III.3. An evacuation network is a road network that specifies E, S, T ⊂ N , representing a set of source, safe and transit nodes respectively. Furthermore, for each source node k ∈ E, let W (k) and d k represent the set of evacuees and the number of evacuees at source k respectively. Let W denote the set of all evacuees.
For the purpose of scheduling an evacuation, we observe that once an evacuee has left their home, it is difficult for them to pause until they reach their desired destination. We also assume that people from the same location evacuate to the same destination. Similarly, we assume that if two evacuation routes meet, they should both be directed to continue to the same location.
Definition III.4. Given an evacuation network, we say a valid dynamic flow is an evacuation schedule if the following are satisfied: • all evacuees end up at some safe node, • no single dynamic flow has any intermediary wait-time (i.e. a v − a u = T uv and, • the underlying flow (without considering time) is confluent, where if two single dynamic flows use the same vertex (possibly at different times), their underlying path afterwards is identical.
A natural objective to optimize for during evacuation planning is the metric average evacuation time of the evacuees. Optimizing for this objective ensures that the evacuation time of the evacuees, on average, is as small as possible. To define it formally, let t i denote the evacuation time of evacuee i. We can then formally define the following problem: Problem 1. Average Dynamic Confluent Flow Problem (A-DCFP). Given an evacuation network, let T max represent an upper bound on evacuation time. Find an evacuation schedule such that all evacuees arrive at some safe node before time T max while minimizing 1 |W | i∈W t i . We formally define two other planning problems. Minimizing Evacuation Completion Time: i.e. time when the last evacuee reaches safety, is another natural objective to optimize during evacuation planning. Formally:

Problem 2. Completion Time Dynamic Confluent FLow
Problem (CT-DCFP). Given an evacuation network, find an evacuation schedule such that all evacuees arrive at some safe node while minimizing max i∈W t i .
Minimizing the Average/Total Evacuation Time of pfraction of the evacuees: In evacuation scenarios, some evacuees may be in such a situation that the cost of evacuating them may dramatically increase the overall evacuation objective, e.g. the Average Evacuation Time. We consider such evacuees as 'outliers'. A common way to handle outliers is to optimize the desired objective for the 'non-outlier' evacuees, while  taking into consideration that the 'outlier' evacuees will also evacuate and use the same road network. We formally define the problem as follows:

A. Time Expanded Graph for Capturing Flow Over Time
Joint routing and scheduling over networks requires one to study flows over time. However, capturing the flow of evacuees over time using only static flows on the road network is challenging. For instance, let us consider the sample evacuation network shown in Figure 1a. All three edges in this network have a capacity of 1, which means 1 car can enter the link in a single timestep. However, sending flow from both sources (0 and 1) at a rate of 1 evacuee per timestep will not work because then two evacuees will reach node 2 at the same timestep but only one evacuee will be able to enter edge (2, A). The main issue here is that we need to keep track of available capacity on the edges at different points in time. We cannot do it using static flows because static flows have the underlying assumption that flows travel instantaneously. To address this issue, researchers have defined dynamic flows ( [22]- [24]) and used time expanded graphs to solve dynamic flow problems ( [6], [8], [11], [12]). In this paper, we also use a time expanded graph (TEG) to capture the flow of evacuees over time.
Time expanded graph is a directed graph denoted by G To construct it, we first fix a time horizon H and discretize the temporal domain into discrete timesteps of equal length. Then we create copies of each node at each timestep within H. After that, for each edge e(u, v) in the road network, we create edges in the TEG as e t (u t , v t+Te ) for each t ≤ H − T e where the edges e t have the same flow capacity as e. Finally, we add a super sink node v t that connects to the nodes u t for each u ∈ S and each t ≤ H. Edges to the super sink node are assigned an infinite amount of capacity. Note that, when creating the time expanded graph, we are adding an additional dimension (i.e. time) to the road network. The size of the TEG is about H times as large as the road network in terms of number of nodes and edges. -yielding a substantially larger problem representation. For instance, the original Houston network with about 1330 nodes, a discretization step of 2 minutes and an evacuation horizon of 15 hours yields a time expanded graph with about 600,000 nodes -yielding a 400 times larger TEG. This time expanded representation is one important reason for the underlying computational space and time complexity of the problem and motivates the need for efficient heuristics.
A sample evacuation network and its corresponding TEG with time horizon H = 3 are shown in Figure (1a-1b). The source, safe and transit nodes are denoted by squares, triangles, and circles respectively. In the TEG, there may be some nodes that are (i) not reachable from the source nodes, or (ii) no safe node can be reached from these nodes within the time horizon. These nodes are greyed out in Figure 1b.
An optimal solution of A-DCFP (and CT-DCFP) for this sample problem instance is to use the routes 0 → 2 → A from source node 0 and 1 → 2 → A from source node 1, where the evacuee at source node 0 and 1 leave at timestep 1 and 0 respectively.

B. Mixed Integer Program (MIP) Model
In this section, we present the Mixed Integer Program (1-8) that we use to represent a class of evacuation planning problems. As shown in objective (1), we can have different objectives in the program, each representing a certain planning problem. Here, we use two types of variables: (i) Binary variable x e , ∀e ∈ A, which will be equal to one if and only if the edge e is used for evacuation. Otherwise, it will be zero. (ii) Continuous variable φ e , ∀e ∈ A x , which denotes the flow of evacuees on edge e.

Objective to Optimize
(1) The constraints of the model are explained in Table I. The constraint that evacuation completion time needs to be less Explanation Constraint (2) Ensures that there is only one outgoing edge from each evacuation node.

Constraint (3)
Ensures that at each transit node, there is at most one outgoing edge. This is necessary for convergent routes.
Constraint (4) Ensures that the total flow coming out of every evacuation node is equal to the number of evacuees at the corresponding node.
Constraint (5) Flow conservation constraint, incoming flow equals outgoing flow. δ − (i) and δ + (i) denote the set of incoming and outgoing edges to/from node i, respectively.
Constraint (6) Flow capacity constraint, flow on an edge will not exceed its capacity. Also, flows are only allowed on assigned edges.
Constraint (8) Binary edge assignment variables. To solve A-DCFP using model (1)(2)(3)(4)(5)(6)(7)(8), we need to minimize the average evacuation time over all evacuees. To do that, we first represent the total evacuation time over all evacuees using the variables x e and φ e as follows: Here, t s (e) denotes the timestep of the starting node of edge e. Dividing (9) by the total number of evacuees will give us the average evacuation time. Note that, minimizing the average evacuation time and the total evacuation time are equivalent as the total number of evacuees is a constant. So, the A-DCFP objective in our MIP model would be: Here, we have just provided details on how to formulate A-DCFP. Details on how CT-DCFP and O-DCFP are formulated as MIPs are provided in the Appendix of the full version [25]. Each of the three problems can be solved using our proposed algorithm MIP-LNS.

IV. INAPPROXIMABILITY RESULTS
In this section, we show that the problems we consider are not only NP-hard but also hard to approximate. Even when we consider special planar graphs that perhaps is closer to a city's road network where G is a subgraph of a grid and all destinations are along the border, these problems remain NPhard. A summary of the hardness results is found in Table II.   Theorem 3. For A-DCFP and O-DCFP with many sources and one safe node, it is NP-hard to approximate within a factor of O(log n).
To prove Theorem 1, we rely on the general k-Node-Disjoint Path Problem: Problem 4 (k-Node-Disjoint Path Problem (kNDP)). Given a graph G, a set of k source-sink pairs s i , t i , find node-disjoint paths from each source s i to sink t i .
The above problem was proven to be hard even when G is a subgraph of a grid [26].
Theorem 4. The k-Node-Disjoint Path Problem is hard to approximate to a factor of 2 Ω( √ log n) even when the graph is a subgraph of a grid and all sources lie on the outer boundary.
The proofs for Theorem 1, 2, and 3 are provided in the Appendix of the full version [25].

V. HEURISTIC OPTIMIZATION
As shown in Section IV, solving A-DCFP, CT-DCFP, and O-DCFP is computationally hard. For this reason, we present the scalable method MIP-LNS where we use MIP solvers in conjunction with combinatorial methods.
In MIP-LNS, we first calculate an initial feasible solution in two steps: (i) calculating an initial convergent route set, and (ii) calculating the schedule that minimizes the target objective using the initial routeset. To calculate the initial route set we take the shortest path from each source to its nearest safe node by road. To calculate the schedule, we use the just calculated route set to fix the binary variables x e in model (1)(2)(3)(4)(5)(6)(7)(8). This gives us a linear program that can be solved optimally to get the schedule.
Next, we start searching for better solutions in the neighborhood of the solution at hand (Algorithm 1). Here, we run n iterations. In each iteration, we select q = (100 − p)% of source locations uniformly at random and keep their routes fixed. This reduces the size of the MIP as we have fixed values for a subset of the variables. We then optimize this 'reduced' MIP model using the Gurobi [27] MIP solver. Essentially, we are searching for a better solution in the neighborhood where the selected q% routes are already decided. Any solution found in the process will also be a feasible solution for the original problem. If we find a better solution with an evacuation completion time T that is less than the current time horizon (T ), then we also update the model by setting the time horizon to T . When resetting the time horizon, we (i) remove edges in the time expanded graph whose start or end node have a time (ii) nodes from which none of the safe nodes can be reached within time horizon T 10 p ← p + p inc 11 return sol stamp greater than T , and (ii) we prune the TEG by removing nodes that are unreachable from the evacuation nodes, and nodes from which none of the safe nodes can be reached within time T . This pruning process reduces the number of variables in the MIP model and simplifies the constraints. At the end of each iteration, we increase the value of p by p inc amount. Note that, when p = 100, we will be solving the original optimization problem. In our experiments, we set the initial value of p to 75 and used p inc = 0.5.
When solving the reduced problem in each iteration (line 4), we use (i) a time limit, and (ii) a parameter threshold gap to decide when to stop. MIP solvers keep track of an upper bound (Z U ) (provided by the current best solution) and a lower bound (Z L ) (obtained by solving relaxed LP problems) of the objective value. We stop the optimization when the relative gap (Z U − Z L )/Z U becomes smaller than the threshold gap. In our experiments, we set this threshold to 5%. In total, MIP-LNS has four parameters: n, p, p inc , and threshold gap. In some iterations, it may happen that the current solution is already within the threshold gap. In that case, the algorithm will simply continue to the next iteration.

VI. EXPERIMENTS
In this section, we first present details of our problem instance (Section VI-A). Then, we present our experiment results with A-DCFP and compare the solutions of MIP-LNS with the baseline method (Section VI-B). Section (VI-C) contains simulation results showing the efficacy and robustness of our solution. We also present the characteristics of the solutions found by solving CT-DCFP and O-DCFP in Section (VI-D). Experiment results showing the scalability of MIP-LNS are provided in Section (VI-E). Finally, Section (VI-F) examines the effect of road failure on the performance of evacuation plans.

A. Problem Instance
In our experiments, we used Harris County in Houston, Texas as our study area. It is situated on the Gulf of Mexico and is prone to hurricanes every year during hurricane season. We have used data from HERE maps [4] to construct the road network. The network contains roads of five (1 to 5) different function classes, which correspond to different types of roads. For instance, function class 1 roads are the major highways or freeways, and function class 5 roads are roads in residential areas. The road network has a hierarchical structure where lower-level roads (e.g., function class 3/4) and the higher-level roads (e.g., function class 1/2) are connected through entrance and exit ramps.
For the purposes of our experiments, we consider the nodes (of the road network) which connect and lead from function class 3/4 roads to function class 1/2 roads as the start/source locations of the evacuees (i.e., the evacuation nodes). We then consider the problem of (i) when should evacuees target to enter the function class 1/2 roads and (ii) how to route them through the function class 1/2 roads to safety. As safe locations, we selected eight locations at the periphery of Harris County which are on major roads. A visualization of the dataset is presented in Figure 2. Additional details regarding the road network, the study area, the time expanded graph and the corresponding MIP model are provided in Table III  We use a synthetic population (as described by Adiga et al. [5]) to determine the location of the households. We consider that one vehicle is used per household for evacuation. From this data, we first extract the location of each household. Then, we assign the nearest exit ramp to each household as their source location.

B. Algorithm Execution and Comparison with Baseline
We performed all our experiments and subsequent analyses on a high-performance computing cluster, with 128GB RAM Fig. 2: Harris County Problem Instance and 4 CPU cores allocated to our tasks. In addition to MIP-LNS, we used two more methods to solve A-DCFP. We used a time limit of one hour for each method and compared the best solutions found within this time. The three methods we experimented with are: 1) Gurobi MIP solver to solve A-DCFP directly.

2) Benders decomposition method, designed by Hasan and
Van Hentenryck [6]. We repurposed this method to solve A-DCFP. To the best of our knowledge, there is no publicly available implementation of the method (or the dataset used by the authors). We, therefore, implemented it and plan to make it public.
2) Benders Decomposition: Hasan and Hentenryck [6] presented Benders decomposition method to solve the 'Convergent Evacuation Planning' problem. Their problem is similar to A-DCFP, differing in the objective function, which is maximizing flow of evacuees instead of minimizing average evacuation time. We repurposed their method to solve A-DCFP and used it as our baseline. At the end of the one hour time limit, we got a solution with average evacuation time of 2.54 hours and evacuation completion time of 7.83 hours.
An advantage of using Benders method is that, in addition to finding solutions to a problem, it can provide lower bounds for the objective value. However, when we used it for our problem instance of A-DCFP, the method was not able to come up with a lower bound better than the trivial value of zero, due to the size of the problem instance. We, therefore, used the lower bound found by Gurobi to calculate the optimality guarantee of the solution. The optimality guarantee was ∼ 20.47%.
3) MIP-LNS: In our experiments with MIP-LNS, for A-DCFP, we used thirty iterations with a total time limit of one hour (same as the baseline). Also, since we have a random selection process within MIP-LNS, we ran ten experiment runs with different seeds. To compare the quality of our solutions with the baseline, we use three different metrics: average evacuation time, evacuation completion time, and optimality Metric Baseline [6] MIP   [6] in terms of three metrics: average evacuation time, evacuation completion time and optimality guarantee. Even the worst solution from MIP-LNS outperforms the baseline in terms of all three metrics. On average, we see a ∼ 13%, ∼ 21%, and ∼ 58% improvement in the three metrics respectively.  guarantee. Optimality guarantee is defined as follows: let the objective value of the solution sol be z sol and the optimal objective be z opt . Then, the optimality guarantee of sol is (z sol − z opt )/z sol , i.e. the smaller the value of optimality guarantee, the better. Table IV shows a comparison of our solutions with the baseline ones in terms of the three metrics. We observe that even the worst solution from MIP-LNS, over the ten experiment runs, is better than the baseline in terms of all three metrics. The best and the average result from MIP-LNS, therefore, also outperform the baseline. Let the value of a metric m for the baseline solution and the MIP-LNS solution be m base and m lns respectively. Then, we quantify the improvement over the baseline as (m base − m lns )/m base . On average, we see an improvement of 13%, 21%, and 58% over the baseline in the three above-mentioned metrics respectively. This indicates that MIP-LNS finds better solutions than the baseline within the given time limit.
To visualize the progress of MIP-LNS over the iterations, we look at the experiment run that returned the best solution.

C. Evaluating the solutions using an agent-based simulation
To evaluate the results returned by MIP-LNS, we run an agent-based simulation of the evacuation. We use the QueST simulator presented by Islam et al. [7] for this purpose. It is an agent-based queuing network simulation system where the roads in the road network are represented by queues and evacuees are agents that traverse these queues. Evacuation routes and schedule (i.e. departure time of evacuees from their initial location) are provided as input to the simulator. We use QueST for the following two purposes: 1) Compare the MIP-LNS and the baseline solution 2) Examine the the effect of non-compliance of evacuees to their designated schedule

1) Comparison of the MIP-LNS and the Baseline Solution:
In this experiment, we performed two simulation runs: one with the routes and schedule from the best MIP-LNS solution, and the other with the routes and schedule from the baseline solution. Figure 3b shows a comparison of the results from the two simulation runs in terms of percentage of people evacuated with time. We see that the blue curve (corresponding to MIP-LNS) is steeper than the orange (corresponding to the baseline). This indicates that the MIP-LNS solution induces a higher rate of evacuation than the baseline.
We also see the difference in evacuation completion time, i.e. ∼ 5.91 hours for the MIP-LNS solution and ∼ 7.88 hours for the baseline solution. Note that, the completion times reported by the simulator are slightly different from the values reported by the solution methods. This is because, in our problem formulation, we discretize the temporal domain using a time unit of two minutes and therefore we lose some precision. The QueST simulator, on the other hand has arbitrary precision which makes it more accurate.
2) Robustness of the solution wrt. Prescribed Departure Time: In real-life evacuation scenarios, even if a schedule is prescribed, it is expected that not everyone will follow it exactly. This may adversely affect the performance of an evacuation plan. To test this, we ran simulations of evacuation in our study area, where people may decide not to start evacuation at their prescribed time.
We use two types of random sampling to decide the departure time evacuees: (i) Normal Distribution, and (ii) Uniform Random Distribution. Let µ i denote the prescribed departure time for evacuee i. Then, in the first sampling method, evacuee i can decide to depart at any random time, sampled from the normal distribution N (µ i , σ 2 ). Here, µ i and σ are the mean and standard deviation respectively. We have experimented with three different values for σ: 0.1, 0.5, and 1 hour. We performed one experiment run for each choice of σ. Fig. 4: Box-plots showing the departure time of evacuees in different schedules. We see that in the random schedules (four boxplots on the right) some evacuees start their evacuation later than their prescribed departure time. We have not allowed evacuees to leave before the declared evacuation start time (i.e. hour zero). In the uniform random sampling method, we sample the departure time of each evacuee from the uniform distribution over the range of 0 to 6 hours. We chose this range because in the original plan it takes about 6 hours for everyone to reach safety. Figure 4 shows the departure time of evacuees in the different schedules using boxplots. Figure 5 shows the arrival time of evacuees at safe locations when they follow different schedules. We observe that: 1) The estimation of our MIP formulation for arrival time of evacuees at safe locations very closely matches the result of our simulation (first two boxplots from the left). 2) We see that when evacuees pick their departure time by random sampling from normal distributions, the evacuation rate remains almost same; i.e. the four boxplots in the middle match very closely. However, we do see some outlier evacuees who reach safety quite late.
3) The uniform random schedule induces a slower rate of evacuation with a median value of ∼ 3.36 hours (compared to the median value of ∼ 2 hours for the other schedules).
The above observations suggest that deviation (within lim-  Figures 6a, 6b, 6c show the comparison in terms of the three metrics average evacuation time, evacuation completion time, and non-outlier average evacuation time, respectively. We observe that the A-DCFP, CT-DCFP, and O-DCFP solutions are superior compared to the other solutions in terms of the objective they are designed to optimize, respectively. its) from the prescribed schedule degrades the performance of the prescribed plan minimally, in terms of evacuation rate. In other words, the resulting schedule and the routes are robust to reasonble perturbations. On the other hand, when individuals choose their departure time uniformly at random, we do see a degradation in overall performance. This is somewhat expected as their departure time does not take the collective behavior of the individuals into account.

D. CT-DCFP and O-DCFP Solutions
In this section, we provide experiment results on the CT-DCFP and O-DCFP problems.
We ran ten experiment runs (with different random seeds) of MIP-LNS for both CT-DCFP and O-DCFP. We again used the one hour time limit for these experiment runs. Figures  6a, 6b, 6c shows a comparison of the A-DCFP, CT-DCFP, and O-DCFP solutions in terms of the three metrics: average evacuation time, evacuation completion time, and non-outlier average evacuation time. Note that, these metrics are also the objectives of A-DCFP, CT-DCFP, and O-DCFP respectively.
From the box-plots we observe that, in general, A-DCFP, CT-DCFP, and O-DCFP solutions are superior compared to the other solutions in terms of the objective they are designed to optimize, respectively. For instance, A-DCFP solutions are, in general, superior than CT-DCFP and O-DCFP solutions in terms of average evacuation time, but worse in terms of the other two objectives. This shows the effectiveness of MIP-LNS and of our formulation in solving evacuation planning problems, while optimizing for different objectives.

E. MIP-LNS Scalability Study
To test the scalability of MIP-LNS we constructed two additional problem instances of varying sizes as shown in Figure 7. The size of these problem instances are given in Table V. Details of the original problem instance are also provided for reference. Our goal is to run MIP-LNS on these three instances and observe how the run time as well as the quality of the solutions found change.
We ran our experiments on A-DCFP and performed ten experiment runs for each of the three problem instances. Fig. 7: Problem instances for our scalability study. Instance 1 consists of the area surrounded by the smallest polygon at the center, the source nodes within it and the blue safe nodes. Instance 2 consists of the orange polygon, the source nodes within and the orange safe nodes. Instance 3 is the same as our original problem instance.    Figure (8a) shows the average run time of MIP-LNS for each problem instance vs the number of binary variables in the those instances. We observe that, starting from problem instance 1, average run time increases almost linearly (with a slope of ∼ 2) with the number of binary variables.
Figure (8a) also shows the run time of MIP-LNS, and Figure (8b), the optimality guarantee of the solutions in each experiment run. We observe that the run time of MIP-LNS and the optimality guarantee of the solution found for problem instance 1 are both very small over all experiment runs. In fact, the optimality guarantee is 0.002% in all runs which means the solutions found are near optimal. This suggests that for smaller instances, MIP-LNS can find almost optimal solutions and it can do so very quickly. We do see more variance in run time and in optimality guarantee for problem instance 2 and 3. However, optimality guarantee is less than 6% over all runs and all problem instances.
In summary, the results in this study suggest that MIP-LNS can find near optimal solutions for small problem instances very quickly. The run time of MIP-LNS increases almost linearly (with a slope of ∼ 2) with the number of binary variables in the MIP formulation, when we start from our (a) The blue plot shows average run time of MIP-LNS vs the NO. of binary variables in the three problem instances. After instance 1, the run time increases linearly (with slope ∼ 2) with the NO. of binary variables. The box-plots show MIP-LNS run time over ten experiment runs for each problem instance.
(b) Boxplot showing optimality guarantee of the final solution of MIP-LNS over ten experiment runs for each problem instance. We observe that an exact optimal solution for instance 1 is found in all ten experiment runs. Optimality guarantee is less than 6% over all runs and all problem instances. smallest problem instance. However, the optimality guarantee does remain reasonably good. We would like to mention that we can always tune the four parameters of MIP-LNS and control the trade-off between run time and solution quality.

F. Effect of Road Failures
So far in our methodology, we have not considered failure of edges. Here, we experiment with the failure of a certain edge and observe its impact on evacuation plans. Specifically, we consider the failure of the Fred Hartman bridge (that spans the Houston ship channel), two hours from the beginning of the evacuation. We selected this bridge because it is used by about 25,000 people in our prescribed evacuation plan.
We consider three different evacuation plans; (i) The original plan we get from MIP-LNS, (ii) Calculate new plan where we avoid the bridge entirely, and (iii) Calculate new plan assuming that the bridge will be unavailable after two hours, this can be done by deleting copies of the bridge from the time expanded graph after the two hour mark. Results of following these plans are shown in Table VI.
The results show that if we have information on the risk of the roads at different times, then MIP-LNS can use it to design better evacuation plans than avoiding roads completely in anticipation of failure. For instance, Agarwal et al. [28] presented a model where the probability of a resource being damaged depends on the spatial distance between the location of disaster and the location of the resource. However, they did not consider the temporal domain, i.e., when the resource is likely to become damaged. A natural extension of their model was proposed by Islam et al. [29] where the authors proposed a time-varying risk model.

VII. CONCLUSION
In this paper, we have proposed a scalable general-purpose optimization method MIP-LNS to solve a class of evacuation planning problems. We showed that it finds near optimal solutions very quickly for small problem instances. We then  demonstrated the scalability of our method by applying it on a problem instance that is at least five times larger than instances previously considered. We also showed that, for our problem instance, MIP-LNS finds solutions that are on average 13%, 21% and 58% better than the baseline method in terms of average evacuation time, evacuation completion time and optimality guarantee of the solution respectively. Using agent-based simulations, we demonstrated that, to a certain extent, our solution is robust to non-compliance of evacuees in following their prescribed schedule. We have also showed how our method can use information regarding failure of roads to come up with better evacuation plans.

A. Formulating CT-DCFP and O-DCFP as MIPs
In this section, we present the details of how we formulate CT-DCFP and O-DCFP as MIPs. Note that, both of these problems can be solved using MIP-LNS without making any changes to the algorithm.
1) CT-DCFP: Evacuation completion time cannot be represented as a linear function of the variables x, φ. To optimize for this objective we modify the time expanded graph as follows: we add a new node z t for each timestep t ≤ H. Then, we add an edge from each safe node (node, t) ∈ S x to the node z t . Finally, we add add an edge from node z t , ∀t ≤ H to the super-sink (v t ). Figure 9 shows an example with two safe nodes on how to do this construction.   φ e > 0 ⇐⇒ y e = 1. We can do it by adding the following constraints to model (1)(2)(3)(4)(5)(6)(7)(8): Here, is a small positive constant (e.g. 0.001) and M is a positive constant which is equal to the total number of evacuees. With these new variables and constraints, we can now represent ECT as: By minimizing (13) with model (2)(3)(4)(5)(6)(7)(8)(11)(12) we can achieve the goal of minimizing evacuation completion time.
2) O-DCFP: To optimize the objective of O-DCFP, we extend the time expanded graph (TEG) as follows: First, we add a new 'bypass' node (A , 0) to the TEG. Then, from each safe node s ∈ S x , we add an edge to (A , 0) with infinite capacity. Finally, we add an edge from (A , 0) to the supersink (v t ) with capacity ( k∈E d k ) * (100 − p)%. Figure 10 shows the newly constructed time expanded graph for our sample problem instance. With this new TEG, minimizing objective (9) will give us the desired solution. Here, the model solver will decide which (100 − p)% of the evacuees are the outliers and then direct them through the 'bypass' node. As the timestep of this bypass node is zero, the cost for the outlier evacuees will not be added to the final objective. However, as they still have to reach the safe nodes through the road network, they will contribute to the congestion on the road.

B. Proof of Hardness
We point out that CT-DCFP is equivalent to the Confluent Quickest Flow Problem in [18], which the authors have shown the above inapproximability result. In addition, they have also provided a bicriteria hardness result where it is not possible to approximate the completion time within some constant factor even when satisfying only some constant fraction of the demand. For O-DCFP, note that it is a generalization of Problem A-DCFP where p = 1. Therefore, any hardness result for A-DCFP also holds for O-DCFP. Thus, we focus on the following two theorems: The approximation hardness proof of A-DCFP is similar to the one in [18]. The main difference is in its analysis since the objective in the two problems differ. For brevity and completeness, we outline the reductions but omit certain proofs and ask the readers to refer to [18] for more details.
The approximation hardness result relies on the NPhardness of the capacitated version Two Distinct Path Problem.
Problem 5 (Two Distinct Path Problem). Let G be a graph with two sources x 1 , x 2 and two sinks y 1 , y 2 . Every edge is labelled with either 1 or 2. Determine if there exists two edgedisjoint paths P 1 , P 2 such that P i is a path from x i to y i for i = 1, 2 and P 2 only uses edges with label 2 (P 1 can use any edge).
The above problem is known to be NP-hard [30]. Other variations of the problem such as uncapacitated, undirected/directed, edge/node-disjoint paths are also known to be hard (see e.g. [31], [32] and [33]).
Proof of Theorem 2. Given an instance I of the Two Disjoint Path Problem, consider constructing the following graph G where we attach safe node t to y 1 , y 2 with an edge of capacity 1, 2 respectively. For i = 1, 2, we also add a source x i and attach it to x i with an edge of capacity i. Every edge with label i also has capacity i. Sources s − i has M * i evacuees for some large M and each edge has a travel time of 1. The upperbound completion time is set to be M 2 n.
If there exists two disjoint paths in I, then a valid schedule simply sends i evacuees at every time step, where the last group of people leaves their sources at time M . Since each path has length at most n = |V (G)|, the i evacuee that left source s i at time k is guaranteed to arrive by time k + n. Then, the total evacuation time is at most M k=1 3(k + n) = 3M 2 /2 + 3M/2 + 3M n, resulting in an average evacuation time of roughly M/2.
If I does not have two disjoint paths, then consider the two paths P 1 , P 2 in a solution to A-DCFP in G. If P 1 , P 2 intersects before t, since it is a confluent flow, the edge following their node of intersection is a single-edge cut that separates the sources from the sink. If the two paths only intersects at t, since we are in a NO-instance of I, P 2 must have used an edge of capacity 1. Then, deleting that edge along with s 1 x 1 also separates the sources from the sink. Note that in both cases, the cut has capacity at most 2. Then, at every time step, at most 2 evacuees can cross the cut. Thus, for all 3M evacuees to cross the cut, it takes at least 3M/2 time steps. Due to this bottleneck, it follows that the smallest total evacuation time is at least 3M/2 k=1 k ≥ 9M 2 /4, giving an average of at least 3M/4.
Since the two instances has a gap of 3/2− where depends on the choice of M , our result follows.
The proof of Theorem 3 follows a similar structure. We use the same setup as the proof of log-hardness of Quickest Flow Time in [34] (Theorem 7) with an arbitrarily large upperbound on completion time. Refer to [34] for more details. Note that the resulting graph contains N sources, one safe node and a total of M 2 log n evacuees. We can similarly show that a YESinstance of the Two-Disjoint Path problem leads to an average evacuation time of at most M 2 /2. Meanwhile, if it is a NOinstance, then by Lemma 23 in [34], there is cut of capacity at most 2, creating a bottleneck. Using similar analysis as before, one can show that the average evacuation time is at least M 2 log N/4. Then, theorem 3 follows immediately.
Proof of Theorem 1. The proofs for all three problems are similar and hence we only focus on A-DCFP here. Given an instance of kNDP where G is a subgraph of a grid and k = O(n), we first swap the location of the sources and sinks such that all sinks lie on the outer boundary. One can easily check that subdividing any edge and adding a leaf to any vertex of degree less than 4 still ensures that G remains a subgraph of a grid. Then, we claim that without loss of generaltiy, we can assume that all sources and sinks are degree-1 vertices. This can be accomplished by first subdivide every edge and shift the sources/sinks to an edjacent newly added vertex. Then, add a single edge to it and shift the source/sink to the new pendant vertex. This ensures that every source and sink is incident to only one edge. For any edge incident to a source s i or a sink t i , assign it a capacity of i; every other edge has a capacity of k. Each source s i has demand M * i where M ≥ n 3 . This means in total, there are M * k 2 /2 evacuees.
In a YES-instance, every source follows its designated path. At every timestep, each source s i sends i people, ensuring a total of k 2 /2 people leaves the sources every timestep. Since each disjoint path has length at most n, anyone leaving at time t is guaranteed to arrive by time t + n. This implies that the last group, leaving at time M also arrives by time M + n, achieving an average arrival time of at most 1 M * k 2 /2 (k 2 /2) * M t=1 (t + n) ≤ M/4 + n. In a NO-instance, in a solution to A-DCFP, consider the cut formed by the edges incident to the sinks. We claim that the amount of flow through this cut at any point in time is at most k 2 /2 − 1. Assume for the sake of contradiction that there exists some point in time where the flow across the cut is at least k 2 /2. Since the total capacity of the cut is k 2 /2, then every edge is used at full capacity at that point in time; in particular, every sink is used in the final routing. Since the routes are confluent, every source is matched to a distinct route. Since we are in a NO-instance, there exists i, j such that source s i is routed to sink t j and i < j. Since source s i can send at most i flow at any point in time, sink t j can never receive its full capacity of j, a contradiction.
Then, at any point in time, at most k 2 − 1 people can arrive at the destinations. This implies we need at least T = M * k 2 /(2(k 2 − 1)) rounds and thus a total evacuation time of at least (T 2 /2)(k 2 − 1). This implies an average evacuation time of at least M k 2 /(4(k 2 − 1)) = M/4 + M/(4(k 2 − 1)). By our choice of M , M/(4(k 2 − 1)) ≥ n, causing a gap between the YES and NO instances, proving our theorem.