Simulation-Assisted Optimization for Large-Scale Evacuation Planning with Congestion-Dependent Delays

First, 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 (Section 5). It is designed based on the Large Neighborhood Search (LNS) framework. In this paper, we work with three objectives: minimizing ($i$) average evacuation time, ($ii$) evacuation completion time, and ($iii$) average evacuation time of ‘non-outlier’ evacuees. We show that all of these three optimization problems are hard to approximate within a logarithmic factor. In MIP-LNS, we model the problems as Mixed Integer Programs (MIP) and then find solutions to these programs using a combination of heuristic search and mathematical optimization. A key technical challenge involves modeling flows over time; this is best achieved using time expanded graphs. But it also leads to substantial increase in the size of the MIP and results in a significant increase in computing resources (time and space). This necessitates the need to combine heuristic search methods developed in the AI literature with MIP techniques developed in the OR literature.

Second, we illustrate how our approach can scale to large problem sizes and can be applied to realistic real-world problems. We choose Harris county in Houston, Texas as our study area and apply MIP-LNS. The county has about $1.5$ million households, spans an area of $1,778$ square miles, and has been affected by several hurricanes in the past. We use real-world road network data and a synthetic population data to construct a realistic problem instance. It is ten times larger (in terms of the size of the time expanded graph, and the number of evacuating vehicles) than the problem instance of our baseline Hafiz Hasan and Van Hentenryck (2021). We show that, within a given time limit, MIP-LNS finds solutions for our problem instance that are on average $13\%$, $20.7\%$ and $58.43\%$ better than the baseline method in terms of average evacuation time, evacuation completion time and optimality guarantee of the solution, respectively (Section 7.2).

Finally, we present MIP-LNS-SIM, where we treat the travel time on each edge as a delay parameter and utilize agent-based simulation to learn the parameter values (Section 6). The MIP model, with the learned parameter values, is then solved to find optimized evacuation plans. Agent-based simulations provide a natural approach to capture the delays one incurs due to congestion – dynamic flow problems cannot capture these delays. Our approach is an example of methods that combine simulation and optimization methods (SO) Gosavi and others (2015); Amaran et al. (2016) considered in OR and has become increasingly popular in AI Van Hentenryck (2013); Kambhampati (2020); Doppa (2021). Through our experiments, we show that MIP-LNS-SIM outperforms MIP-LNS in terms of average evacuation time, evacuation completion time, and average time spent on the road ($10\%,17\%,77\%$ improvement respectively) when delay due to congestion is considered (Section 7.3). In addition, MIP-LNS-SIM has a significantly lower percent error ($6\%$) in estimated evacuation completion time compared to MIP-LNS ($76\%$), demonstrating the efficacy of MIP-LNS-SIM in evacuation planning subject to congestion constraints.

Joint routing and scheduling over networks requires one to study flows over time, as static flows make the unrealistic assumption that flows travel instantaneously (detailed discussion in the supplementary materials Islam et al. (2023)). For this purpose, researchers have defined dynamic flows (Skutella (2009); Ford and Fulkerson (2015)) and used time expanded graphs to solve dynamic flow problems (Romanski and Van Hentenryck (2016); Hafiz Hasan and Van Hentenryck (2021)). 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 $\mathcal{G}^{x}=(\mathcal{N}^{x}=\mathcal{E}^{x}\cup\mathcal{T}^{x}\cup\mathcal{S}^{x},\mathcal{A}^{x}$). To construct it, we first fix a time horizon $\mathcal{H}$ and discretize the temporal domain into discrete timesteps of equal length. Then we create copies of each node at each timestep within $\mathcal{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+T_{e}})$ for each $t\leq\mathcal{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\in\mathcal{S}$ and each $t\leq\mathcal{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 $\mathcal{H}$ times as large as the road network in terms of number of nodes and edges. – yielding a substantially larger problem representation.

A sample evacuation network and its corresponding TEG with time horizon $\mathcal{H}=3$ are shown in Figure (1(a)-1(b)). 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 1(b). An optimal solution of a-dcfp (and ct-dcfp) for this sample problem instance is to use the routes $0\rightarrow 2\rightarrow A$ from source node $0$ and $1\rightarrow 2\rightarrow A$ from source node $1$, where the evacuee at source node $0$ and $1$ leave at timestep $1$ and $0$ respectively.

Now, we present the Mixed Integer Program (1–8) that we use to represent a class of evacuation planning problems. We can have different objectives in the program (Objective 1), each representing a certain planning problem. We use two types of variables: ($i$) Binary variable $x_{e},\forall e\in\mathcal{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 $\phi_{e_{t}},\forall e_{t}\in\mathcal{A}^{x}$, which denotes the flow of evacuees on edge $e$ at timestep $t$.

Constraint (2) ensures that there is exactly one outgoing edge from each evacuation node. Constraint (3) ensures that at each transit node, there is at most one outgoing edge. Constraint (4) enforces that the total flow coming out of every evacuation node is equal to the number of evacuees at the corresponding node. Constraint (5) ensures flow conservation in the time-expanded graph; here, $\delta^{-}(i)$ and $\delta^{+}(i)$ denote the set of incoming and outgoing edges to/from node $i$, respectively. Constraint (6) allows flow on chosen edges only; it also enforces flow capacity on each edge of the time-expanded graph. Constraint (7) defines $\phi$ as continuous and non-negative variable; constraint (8) defines $x$ as binary variable. The constraint that evacuation completion time needs to be less than the given upper bound is implicit in the model, as we set the time horizon of the TEG to this upper bound.

To solve a-dcfp using model (1–8), we represent the total evacuation time using the variables $x$ and $\phi$ as follows:

Here, $t_{s}(e)$ denotes the timestep of the starting node of edge $e$. Note that, minimizing the average evacuation time and the total evacuation time are equivalent. So, the a-dcfp objective would be: $\min_{x,\phi}\sum_{e\in\delta^{-}{(v_{t})}}\phi_{e}t_{s}(e)$.

We have just provided details on how to formulate a-dcfp as a MIP. Details on ct-dcfp and o-dcfp are provided in the supplementary materials Islam et al. (2023).

We use Harris County in Houston, Texas as our study area. We have used data from HERE maps HERE (2023) to construct its 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 major highways, and function class 5 roads are residential roads.

For our experiments, we consider the nodes which connect and lead from function class 3/4 roads to function class 1/2 roads as the start/source locations of the evacuees. 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 safely. As safe locations, we selected six 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 problem instance are provided in Table 2.

We use a synthetic population Adiga et al. (2015) to find the location of the households. We then assign the nearest exit ramp to each household as their source location. We assume that one vehicle is used per household for evacuation.

We performed all our experiments and subsequent analyses on a high-performance computing cluster, with 128GB RAM 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 described here.

1. 1.

   Gurobi In this experiment, we used Gurobi to directly solve model (1-8) with the a-dcfp objective. Gurobi was not able to find any feasible solution within the time limit. However, Gurobi was able to come up with a lower bound for the objective value. We used the lower bound to calculate the optimality guarantee of the solutions.

2. 2.

   Benders Decomposition Hafiz Hasan and Van Hentenryck (2021) presented Benders Convergent (BC) 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 and used it as our baseline.

3. 3.

   MIP-LNS In our experiments with MIP-LNS, for a-dcfp, we used thirty iterations (i.e. $n=30$ in Algorithm 1 line 1). 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 metrics: average evacuation time, evacuation completion time, and optimality 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. If $z_{opt}$ is unknown, we can use a lower bound of it. Table 3 shows a comparison of our solutions with the baseline in terms of the three metrics.

Let the value of a metric $m$ for the baseline 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.

We also applied MIP-LNS to find solutions of ct-dcfp and o-dcfp for our problem instance. Due to limited space, we provide the results in the supplementary materials Islam et al. (2023). In general, the experiment results show that MIP-LNS can effectively solve the problems with different objectives.

Within MIP-LNS-SIM, we set the parameter $m=10$. We then experimented with two values for the parameter $p_{e}$, which are $5,10$. We refer to these two settings as MIP-LNS-SIM-$5\%$ and MIP-LNS-SIM-$10\%$. We ran MIP-LNS-SIM with the two settings and found two different solutions. We then performed agent-based simulations of the entire evacuation process (i.e. evacuate $100\%$ of the evacuees) using a solution of MIP-LNS and then also using solutions from MIP-LNS-SIM-$5\%$ and MIP-LNS-SIM-$10\%$. We used the QueST simulator with the logistic traffic model here. We now compare the simulation results.

Figure (3(a)) shows the departure rate of the evacuees from their initial locations, in the final solution of the three settings. We observe that, MIP-LNS-SIM-$5\%$, and MIP-LNS-SIM-$10\%$ regulates the departure of evacuees to a significant extent (compared to MIP-LNS). As evacuees leave late in these solutions, we might expect the evacuation completion time to be higher in these solutions compared to MIP-LNS. Surprisingly, we observe in Figure (3(b)) that the evacuation completion time is actually smaller in MIP-LNS-SIM-$5\%$, and MIP-LNS-SIM-$10\%$ compared to MIP-LNS. This implies that although evacuees left early in the MIP-LNS solution, they could not reach safety early due to the resulting congestion on the roads. In the MIP-LNS-SIM-$5\%$ and MIP-LNS-SIM-$10\%$ solutions, evacuees departed from their initial location over a longer period of time. This way there was less congestion on the road and the evacuation was completed early even though many people started late.

Figure 4 verifies our last statement. We see that traffic density on the edges (i.e. number vehicles per lane and per km) is higher in the MIP-LNS solution (than MIP-LNS-SIM-$5\%$, MIP-LNS-SIM-$10\%$) throughout the evacuation time period. The higher traffic density then causes the evacuees to spend more time on the road. In summary, MIP-LNS-SIM is better than MIP-LNS in terms of average evacuation time, evacuation completion time, and average time spent on the road ($10\%,17\%,77\%$ improvement respectively, detailed results in the supplementary materials Islam et al. (2023)). These results indicate that MIP-LNS-SIM is better at evacuation planning than MIP-LNS in terms of reducing congestion on the roads.

Finally, Table 4 shows the estimated and the simulated evacuation completion time for Baseline Hafiz Hasan and Van Hentenryck (2021), MIP-LNS, MIP-LNS-SIM-$5\%$, and MIP-LNS-SIM-$10\%$. For instance, MIP-LNS predicts that when using the routes and schedule provided by its solution, the evacuation will be completed within $5.77$ hours. However, when simulated, it actually took $24.29$ hours to evacuate everyone. We also observe that the percent error in estimation decreases considerably in MIP-LNS-SIM-$5\%$ and MIP-LNS-SIM-$10\%$ where we have used $p_{e}=5$ and $10$ respectively. The lower percent error is earned at a cost of higher algorithm run time (MIP-LNS-SIM-$5\%$: $\sim 6.5$ hours, and MIP-LNS-SIM-$10\%$: $\sim 7.55$ hours).