EVACUATION DYNAMICS INFLUENCED BY SPREADING HAZARDOUS MATERIAL

. In this article, an evacuation model describing the egress in case of danger is considered. The underlying evacuation model is based on continuous network ﬂows, while the spread of some gaseous hazardous material relies on an advection-diﬀusion equation. The contribution of this work is twofold. First, we introduce a continuous model coupled to the propagation of hazardous material where special cost functions allow for incorporating the predicted spread into an optimal planning of the egress. Optimality can thereby be understood with respect to two diﬀerent measures: fastest egress and safest evacuation. Since this modeling approach leads to a pde/ode-restricted optimization problem, the continuous model is transferred into a discrete network ﬂow model under some linearity assumptions. Second, it is demonstrated that this reformulation results in an eﬃcient algorithm always leading to the global optimum. A computational case study shows beneﬁts and drawbacks of the models for diﬀerent evacuation scenarios.


1.
Introduction. An evacuation process in an emergency situation can be considered as the task of moving residents quickly out of a danger zone to safety. Evacuation modeling is an important tool for estimating the total evacuation time, identifying dangerous substructures in the considered area, and optimizing evacuation strategies (see [34,35,18,14,38]). Evacuation models are usually classified into microscopic and macroscopic approaches. Microscopic models simulate the behavior of each individual evacuee. Interaction between evacuees and local phenomena such as lane formation of pedestrians and bottlenecks can be identified in these models [35,25]. Macroscopic approaches model evacuees as part of a homogeneous group and do not take into account individual behavior. Different aspects of an evacuation scenario such as infrastructure, evacuees, and safe locations are considered. With this overall knowledge, system optima yielding among others minimal evacuation times or fastest evacuation routes can be computed. Macroscopic models based on optimization techniques can be designed such that their results correspond to a best-possible evacuation, that means, e.g., that the evacuation times computed are provable lower bounds on the true evacuation time [18,15].
In most applications, macroscopic evacuation models are based on linear dynamic network flows. The big advantage of this modeling approach is the availability of fast and well-investigated algorithms. Dynamic network flow models are discrete macroscopic models that are capable of yielding optimal evacuation strategies [5,4]. In dynamic networks, flow units are assumed to take time to traverse an arc [37]. This generalization of static network flows was first considered by Ford and Fulkerson [9] who also proposed a method to compute a maximal dynamic flow between two specified nodes using temporally repeated flows computed in the underlying static network.
For evacuation modeling, the quickest flow problem, where the goal is to find the minimum time to bring all flow units in the network to a predefined sink node, is also of importance. It can be solved in polynomial time by combining maximum dynamic flows with an appropriate binary search over the time horizon [3]. Extensions and variants of dynamic network flow problems have been examined to include phenomena emerging in evacuation processes. Hamacher and Tjandra [17] solved the earliest arrival problem with time-dependent capacities, enabling models with temporally failure of arcs. Flow and load dependent travel times reflect slow-down of evacuees on congested arcs and are considered in [27,28]. Lately, microscopic aspects are included in macroscopic models as well as both, microscopic and macroscopic models, are combined in common framework to obtain more reliable information about the evacuation process [16,26].
From a pde point of view, evacuation aspects can be found in the analysis of pedestrian movements, see [19,22,23] for an overview. The majority of these models is formulated on 1d-or 2d corridors but not on networks. To map pdes on networks, it is necessary to define coupling conditions at nodes ensuring the conservation of mass. For instance, this is done in traffic flow [6,20] or supply chain modeling [12,7]. Both approaches are close to the pedestrian models and therefore many ideas can be copied.
In this work, we are interested how the propagation of some hazardous material (e.g. a gas cloud) influences the evacuation process. This means, the evacuation model under consideration is combined with an equation for the hazardous material in a clever way. To solve the coupled problem, we derive a new framework and apply the aforementioned dynamic network flow algorithms. Further, we investigate the numerical properties of the network flow approach by comparing computing times for the solution to the standard mixed-integer approach and demonstrate the efficiency of network flow algorithms in this context.
The remainder of this article is organized as follows. In Section 2, we start from a pde-based evacuation model which is coupled to the expansion of spreading hazardous material described by an advection-diffusion equation. In Section 3 we introduce two solution procedures for the coupled problem: on the one hand a straightforward discretization leading to a mixed-integer model and, on the other hand, a reformulation as a dynamic network flow problem. The numerical results are presented in Section 4 and contain comparisons of the two solution approaches and case studies as well.
2. Continuous evacuation model. To introduce an evacuation model that is suitable for regional evacuation purposes, we need a mathematical definition of a EVACUATION DYNAMICS 445 network where roads correspond to arcs and cities to nodes. In the beginning of the evacuation process, all evacuees are located in a source node from where they shall get to safety. Therefore, we define the following: Definition 2.1 (Network definition). For a time horizon T let G = (V, A) be a directed graph with node set V and arc set A. The node set is assumed to subsume a source s ∈ V and a sink d ∈ V . For a node i ∈ V , δ − (i) and δ + (i) is the set of all arcs starting and ending at i, respectively. An arc e ∈ A is associated with an interval [a e , b e ] of length L e . The constant times needed to traverse one arc are called travel times τ e and correspondingly the travel velocities are v e = Le τe . Arc capacities are denoted by piecewise constant functions µ e (t, x) and storage capacities at nodes by time-dependent functions M e (t).
Remark 1. Mathematically, the model is restricted to subsume one source and one sink only. Note that the model can easily be extended to multiple sources and multiple sinks using a super source and a super sink, respectively.
For our modeling, we prefer a macroscopic model composed of ordinary and partial differential equations. Similar to the ideas used in traffic flow network models [6,21] and the supply chain approach [12,7], we stick to the following pde/ode-restricted optimization problem describing the evacuation scenario: The objective function (1a) must be chosen in an appropriate manner. Note that the function Φ might depend on the density of evacuees ρ e on arcs e ∈ A and the load of the queue q i in vertices i ∈ V , where the queues buffer an excess of net inflow, cf. [2]. Possible choices of the objective function in the context of evacuation are discussed below and in Subsection 2.1. The advection equation (1b) describes the movement of people on arcs and ensures conservation of mass. We additionally assume initial values ρ e (0, x) = ρ(x) e,0 . Usually, depending on the application, it is necessary to prescribe boundary data in form of an influx for e ∈ δ + (s). In this model, the relation between the flux f e and the density ρ e is simply linear, i.e. f e = v e ρ e . That means, natural bottleneck situations as in traffic flow models cannot occur. But as we will see later on, this linear structure will help to interpret the evacuation model in the sense of dynamic network flow models. Equation (1c) determines the load of the queue q i in node i. The ordinary differential equation is a rate equation measuring the total ingoing e∈δ − (i) f e and outgoing flux ẽ∈δ + (i) fẽ into this node. Here, we have degrees of freedom to distribute the flux. It is the subject of optimization to determine the best routing of evacuees through the network. The initial load of each queue i ∈ V is given by q i,0 representing the number of residents in vertex i to be evacuated, see equation (1d).
Box constraints (1e) and (1f) are bounded by the respective capacity constraints from above.
Remark 2. The problem (1) is inspired by the evacuation model motivated in [18]. Due to its linear nature, it is the simplest approach to model an evacuation process. This means, either the evacuees move forward towards their safety destination as long as capacity constraints permit escape or, if the capacity is reached, the evacuees have to stay in nodes. Therefore, congestions on arcs are avoided. However, a more sophisticated model should allow for bottleneck situations on arcs, cf. traffic flow models [6,21].
From literature about evacuation modeling, we know that for a fixed time horizon T , the goal of the so-called maximum dynamic flow problem is to send the maximum possible number of flow units from source to sink. The side constraints to this problem are conservation of mass (1b), (1c) and capacity restrictions (1e), (1f) with q i,0 = 0 for all i ∈ V . The goal is to optimize the following objective function However, for a given amount of flow, the quickest flow problem seeks for the minimum time horizon to send all flow units from source to sink through the network. This can be modeled using (1b) -(1f) with q i,0 = 0 for all i ∈ V \ {s} and q s,0 > 0 and the following objective function min Φ(ρ, q) = sup t ∈ [ 0, ∞) | f e (ρ e (t, b e )) = 0 for some e ∈ δ − (d) . ( Further, it has to be ensured that all evacuees are indeed sent through the network. Thus, an additional condition is introduced The main focus of this work is to couple the evacuation model (1) to the spread of some hazardous material (think of a gas cloud). We present a two-stage modeling approach where we first compute the propagation of gas by an advection-diffusion equation in 2d and second map the concentration as extra costs onto the network. The resulting evacuation model is supplemented with a new objective function involving concentration-depending costs.
2.1. Evacuation influenced by hazardous material. Since the task of evacuating people from a certain area has been intensely studied in the past (see e.g. [18,13]), we consider an advanced scenario where time-dependent parameters enter the problem due to an external source, which will be called hazardous material. This may be a cloud consisting of gas, radioactive material (e.g. when considering the evacuation of the area around a nuclear power plant after an MCA), or an ash cloud from a volcanic eruption. In doing so, we first introduce the equations and parameters describing the hazardous material. Then we discuss possibilities to map the new information onto the network. Finally, we introduce the coupled evacuation problem.
The hazardous equation shall set up additional parameters for the model (1). The original model relies on a general graph. Yet, in view of the application, vertices are specified by coordinates in Ω ⊂ R 2 . Furthermore, we assume that the concentration of hazardous material ω obeys the following linear advection-diffusion equation: with diffusion constant κ > 0 and wind vector field d = (d 1 , d 2 ) ∈ R 2 . The starting level of the concentration is 100% (life-threatening value) in the area of the source and 0% elsewhere. On the boundary of Ω, we use Dirichlet conditions with ω = 0. Note that this way of modeling does not include a time-dependent source and thus, the concentration is fixed to the amount given at the beginning of the calculation.
Remark 3. Since the focus of this subsection is on the mapping of equation (5) on the network, we prevent a discussion on more realistic propagation models.
Due to the simple dynamics of the hazardous material (associated with standard numerics), we now turn our attention to the mapping of the values ω onto the network. We present two ways to include the concentration level as an additional parameter in the evacuation model (1) which lead to a coupling to the original model via a new objective function Φ.

Remark 4.
An alternative coupling might be to adjust the capacities µ e (t, x), M e (t) or travel velocities v e of arcs. E.g. in exposed areas the number of accidents may increase since people are getting nervous. This potentially leads to smaller capacity constraints and lower velocities.
In the following, we concentrate on a new objective function approach and introduce time-dependent concentration parameters c e (t) and c i (t) defined on arcs e and vertices i, respectively. We map the concentration value ω to the corresponding positions in the network, e.g. ω| e means that the values of ω are restricted to the arc e, and comment on two integral functionals, called static model and dynamic model.
The travel time is included in the computation of the cost, because this displays that evacuees on an arc with large travel time are more exposed to the concentration than evacuees on arcs with smaller travel times. For a vertex i ∈ V the concentration parameter at time t represents the exposure of evacuees to the hazardous material when waiting in the vertex for one time step. Thus, we get Correspondingly, for a vertex i, we get Naturally, the dynamic model is closer to reality than the static model. But obviously its drawback is the need of the full information on the hazardous material in advance. Thus, the dynamic model does not allow for an online use. Numerical experiments comparing the static and the dynamic model can be found in Section 4.
An appropriate possibility for including the costs in the evacuation model is to minimize the level of concentration to which the evacuees are exposed. This can be done with the definition of new cost-dependent objective functions Φ(ρ, q, c). On the one hand the maximal level of concentration to which the evacuees are exposed is minimized by (10) This objective function is sensible, if there is a level of concentration which directly causes the loss of evacuees, because the optimization seeks a solution avoiding arcs with high values of concentration costs. For example, if a hazardous gas is considered, which causes the death of evacuees if its concentration exceeds a certain percentage in the air but does no real harm below this critical percentage, this objective function will lead to an evacuation plan avoiding the highly contaminated arcs but therefore using -possibly more -arcs with concentration parameters slightly below the critical value. On the other hand the total amount of concentration is minimized by This objective function is useful when the hazard is harmful if the accumulated level exceeds a critical level (e.g. low levels of concentration). In summary, both objective functions can be used to extend the evacuation model (1) by solving the optimization problem min(10) or (11) (12a) Generally, this type of optimization problem can be solved using adjoint-based methods due to the continuous nature of the problem. In case of a linear objectives and constraints or only 'mild' nonlinearities, it is possible to reformulate the original problem into a mixed-integer program. The latter will be explained in Section 3.
Remark 5. In many applications, it is useful to consider the combination of both objective functions (10) and (11). This will lead to a so-called multi-objective optimization problem that requires more involved solution methods.
3. Solution methods. In Section 2 an evacuation model using partial and ordinary differential equations as constraints was introduced. The usual way to solve such an optimization problem is to calculate the adjoint equations as proposed in [24]. This is a straightforward computational approach first deriving the first order optimality system and then applying a gradient scheme to compute a (local) optimum. Since the evacuation problem is purely linear and the computation of the optimality system is rather complex we prefer another approach following the ideas in [11,10]. Therein, a mixed-integer program is introduced by discretizing the underlying network model using suitable numerical methods. This approach works for the traffic flow and supply chain model pretty well and can also be applied for the evacuation model (1). As it turns out, the evacuation model has nice features (linear dynamics on arcs, conservation of mass through nodes and box constraints) and, thus, it is even possible to put the model (1) into the framework of a dynamic network flow model, cf. [18]. This is a favorable approach, since dynamic network flows are well-investigated (see e.g. [1,8,18,37]) and make use of very efficient algorithms leading to the global optimum. α k e := f e (ρ e (t, a e )), β k e := f e (ρ e (t, b e )), γ k e := q e (t), ∀e, k.
As abbreviation for the original objective function Φ(ρ, q, c), we use a straightforward discretized version Φ disc (α, β, γ, c) depending on the respective variables. Then, the mixed-integer program (MIP) of the evacuation model (1) reads: Let us explain the equations in (14) step by step. A special choice of the objective function -corresponding to the maximum dynamic flow problem -is given by the maximization of flux reaching the sink Consequentially, the problem (14) reduces to a linear program (LP) since the objective functions and all constraints are linear. For the problem (14), we only focus to the max-flow objective, since we are interested in computing times of the LP compared to other approaches, see Section 4. Note that only depending on the objective function the problem (14) will become a mixed-integer program. E.g. the nonlinear objective function (10) requires binary variables making the problem a mixed-integer one. The transport equation (1b) is reformulated with the help of the Origin Tracking Algorithm [29,30]. This algorithm has been invented to track down the origin of contamination in water supply networks, but it is general in nature and we use it to detect the inflow on each arc. Using a unit time discretization and integer travel times simplifies the algorithm and allows for an exact calculation of the time step of the inflow. Then, equations (14b) and (14c) represent the exact solution which is independent of the CFL-condition. For the right implementation of the coupling condition (1c), we need variables ζ k e collecting all influx at x = a e , see (14d). The original ordinary differential equation for the queue is discretized using the explicit Euler method (14e). The remaining constraints remain unchanged and can be discretized directly. Basically, the MIP (14) is solved using the commercial solver of CPLEX [39] which is based on Branch-and-Bound and Branch-and-Cut techniques. For a numerical case study, we refer to Section 4.  (1) is known and we can write Furthermore, the flow conservation constraint (1c) in the evacuation model (1) can be rewritten as follows. The left-hand side ∂ t q i is approximated via the differential quotient whereas the right-hand side is reformulated applying (15) in combination with τ e = (b e − a e )/v e , cf. Definition 2.1. Then, we end up with following equation for a fixed vertex i: Now, we introduce the notation to translate equations (15), (16) in the sense of a typical dynamic network flow model. Let us assume that we have an integer time horizon T ∈ N, an equidistant time-grid where ∆t = 1 and we denote the time steps by t k with 0 = t 0 and t kmax = T . The variables x ij (t k ) describe the flow on arc e = (i, j) ∈ A at time step t k and y i (t k ) the flow waiting in vertex i ∈ V at time step t k . As before, τ ij describes the travel time on arc e = (i, j) and correspondingly a i (t k ) and b ij (t k ) are the discrete capacity limits. In contrast to the continuous evacuation model (1), the parameters τ ij , a i (t k ), b ij (t k ) are restricted to N. To shorten notation we use X = (x ij (t k )) (i,j)∈A,0≤k≤kmax , Y = (y i (t k )) i∈V,0≤k≤kmax for the whole flow. Furthermore we assume x ij (t k ) = 0 for t k < 0. The number of initial residents in each vertex i is given as y i,0 ≥ 0. The source s ∈ V is the vertex in the danger zone, where all the evacuees are located at the beginning of the evacuation process. Thus, as we consider only a single source problem, y s,0 > 0 is the only non-zero value in (y i,0 ) i∈V . Summarizing, we state: In other words: The model (17) computes a flow from the source s to the sink d, shortly s-d-flow, with respect to some objective function φ(X, Y ). The flow decomposition (see [1]) of the resulting flow can then be interpreted as routes used by evacuees to safety.
Finally, we comment on the objective functions φ(X, Y ) for the dynamic network flow model. As in Section 2, we start with a maximum flow for a given time horizon T For this objective the constraints (17b) -(17e) are used and the calculation is started with an empty network, meaning y i,0 = 0 for all i ∈ V . The objective function for the quickest flow problem, minimizing the time horizon to clear a network, as described in literature (cf. [18,3,8]) takes the form min φ(X, Y ) = min max 0 ≤ t k ≤ ∞ | x e (t k − τ e ) = 0 for some e ∈ δ − (d) .
Again we use the constraints (17b) -(17e), additionally equipped with kmax k=0 (i,d)∈A x id (t k ) = y s,0 ensuring that all people are evacuated (i.e. reach the sink d).
But for our applications, it is more interesting to deal with problems where cost terms for the propagation of hazardous material are included. To deduce objective functions similar to (10) and (11), we consider the unit intervals I k = [ t k , t k+1 ) between the time steps t k . Then, for both the static and the dynamic model, we take the mean value of the concentration parametersc ij (t k ) andc i (t k ) for e = (i, j) ∈ A and i ∈ V , respectively, over each interval I k yielding where c e (t) and c i (t) are the time-dependent concentration parameters introduced in Section 2. Thus we can now reformulate the objective functions (10) and (11) yielding max max (i,j)∈A, 0≤k≤kmax| xij (t k )>0c for the minimization of the maximal hazard to which the evacuees are exposed and for the minimization of the total amount of hazard to which the evacuees are exposed. Solution methods for the model (17) are explained in Section 4.2.
4. Computational results. In this section, we first present techniques to solve optimization problems described in the previous sections. We start with the numerical solution of the advection-diffusion equation (Section 4.1), which is a subproblem of the optimization problem (12) and its discrete analogon. Then we introduce some basic concepts for the solution of the dynamic network flow model (17) and sketch algorithms to solve the dynamic network flow model for the objective functions (18), (19), (22) and (23). Furthermore we present a comparison of running times in Section 4.2.2. We close this section with a numerical example comparing the quickest flow with the solutions obtained by considering the objective function (23).

Numerical solution of the advection-diffusion equation.
For the numerical solution of equation (5) we use a five-point difference stencil ∆ h for the Laplacian as a discretization and discretize the auxiliary part with This yields a stable and second order consistent and convergent method as long as κ > h/2|d i | for i = x, y.
If we would now use the explicit method for the time discretization, we would run into problems with memory requirements and long running times as we would get restricted in the size of the time step because of the parabolic CFL condition ∆t h 2 ≤ 1 2 (25) to ensure stability of the method. So, we use the trapezoidal rule for time discretization and obtain the two-dimensional version of the Crank-Nicolson Method This method is implicit and not restricted in the step size ∆t of the time discretization. But the Crank-Nicolson Method responds to the discontinuities in the initial conditions with oscillations (see text to equation (5)). This problem has been studied, for example, by Pearson [32] and Osterby [31]. They show that the complete annihilation of the oscillations is only possible by using the implicit method for the time discretization, which would reduce the order of accuracy from second to first. But Pearson [32] suggests a method which does not aim at a complete annihilation of the oscillations but at cumulative damping in several steps to an acceptable level. Therefore the first step is subdivided into M equal steps of length ∆t = ∆t/M and the resulting cumulative damping is large enough, to get rid of the oscillations. In our case using M = 15 is enough meaning not to much extra computational effort has to be spent to remove the discontinuities.

Evacuation solvers.
After the distribution of the hazard has been calculated and mapped onto the network (see Section 2.1), the full network flow model with all parameters is at hand and needs to be solved. Solution methods for the several objective functions mentioned before will now be presented where we start with the description of the general concept of the time-expanded network (Section 4.2.1), which is similar to the concept of a time grid for partial differential equations. Then, we show how a maximal flow problem is solved using the time-expanded network (Section 4.  (17), a time component is included in the network. This can be modeled explicitly by considering a copy of the network at every possible time step. We map the dynamic network onto a static time-expanded network (TEN) [1] where a feasible static flow corresponds to a feasible dynamic flow in the dynamic network. For simplicity of notation, we omit the subscript k of the time step t k in the following and simply write T instead of t kmax . For the right time discretization we refer to Section 3.2.
The capacities b T e of the arcs e ∈ A T ∪ W T are chosen as The time-expanded network contains T + 1 copies of every vertex i ∈ V of the original network N which represent this vertex at all time steps t ∈ {0, . . . , T }. The time copies are connected in such a way that flow which enters an arc e ∈ A at time t arrives at the end-vertex of this arc exactly at time t + τ e (arcs in A T ). The arcs in W T are called holdover arcs and model the possibility to wait at a vertex for one time step. The resulting network has no time component and algorithms for static networks can be applied. The static flow conservation constraints in the TEN correspond to the dynamic flow conservation constraints of the problem we expanded. Hence, we can interpret a static flow f s in N T as a flow over time (X, Y ) in the dynamic network N by setting x ij (t) = f s (it,jt+τ e ) and y i (t) = f s (it,it+1) . Since the time-expanded network considers every time step explicitly it is straight forward to include time dependent parameters on the arcs.  The drawback of this general method is the size of the time-expanded network. The number of vertices and arcs of the network scale with T . So, although time expansion is a useful tool, using static network flow algorithms on the time-expanded network will not give polynomial-time algorithms for the dynamic network flow problems. Note that for capacities and travel times that are constant over time there exist also polynomial time algorithms for the maximum dynamic flow problem (18) and the quickest flow problem (19) (see [37,3]).
We now use the concept of the time-expanded network to give algorithms for the solution of the optimization problems mentioned in Section 3.

4.2.2.
An algorithm for the maximum dynamic flow problem. To solve the maximum dynamic flow problem we have to find a dynamic s-d-flow which maximizes the flow arriving at d in a given time horizon T . Since the time at which the flow starts in the source and the time at which it arrives at the sink are not important -as long as it is within the time horizon T -for the objective function, we add a super source S which is connected to every time copy of s with infinite capacity, and a super sink D which is connected to every time copy of the sink d with infinite capacity. So we reduce the problem to a single source single sink maximum flow problem in the time-expanded network which can be solved efficiently [1]. This solution is then mapped to a flow over time which is a maximum dynamic flow.
Comparison of CPU times. As the maximum flow problem is on the one hand a basic problem and on the other hand easy to implement for both the discrete model (17) and the MIP (14) (which reduces to a LP in case of the maximum flow problem), we use it to compare the CPU times of both models. For the analysis of the running time, depending on the size of a graph, we use a network shown in Figure 2 to which we refer as the chessboard. It consists of a square of vertices, in which vertices are connected from bottom to top and from left to right. Furthermore to the lower left corner an arc from the source and to the upper right corner an heading to the sink are added. We choose this kind of a network as it can easily be implemented and scaled in a straightforward way. We start the analysis with a graph of five vertices in a square -not counting the source and the sink. The capacities of inflow and outflow arc are µ 1 = µ 42 = 40, their travel times τ 1 = τ 42 = 1. The capacities and travel times of all inner arcs, i.e. all other arcs, are each chosen randomly between one and four from an equal distribution. The travel times are kept constant over time but the capacities are chosen for each time step up to the time horizon T = 160. As we choose the parameters of the inner arcs randomly and the running time may be affected by that choice we throw the die ten times. In that way we generate ten graphs with n = 27 vertices and m = 42 arcs with different parameters, which we call G 1 5 , G 2 5 , . . . , G 10 5 -the subscript indicating the size of the square and superscript the number of the die roll. We extend these graphs by adding five vertices to the upper and the right end of the square to end up with a square of ten times ten vertices. We only choose the parameters of the inner arcs which were no element of the corresponding G i 5 randomly between one and four, the other parameters are taken from G i 5 . We repeat this procedure until we reach a size of 50 for the chessboard yielding the following chain of sub-and supergraphs The number of vertices n and the number of arcs m for each graph G i g is calculated for all i by Thus, the largest network has n = 2502 vertices and m = 4902 arcs. For all graphs G i g we compute the optimal solutions for both the discrete model using the time-expanded network (TEN) and the PDE-model being simplified as a LP shown in equation (14). The TEN is solved using the Boost-library [36] whereas the LP is solved using the commercial solver of CPLEX [39].
The optimal objective values of the solutions gained by the TEN and LP coincide for each graph G i g , so we only present the mean solution times (CPU times of the solvers) in Table 1 and Figure 3, where we take the mean over the different die rolls Table 1. Average values of solution times for graphs of equal size for TEN (discrete model) and LP (PDE-model) i for each graph size g. But -because the optimal solution is not unique -the optimal solutions for each graph G i g gained by the TEN and LP themselves differ. We see that the LP is not capable to process the graphs of sizes g > 35, which is due to excessive amount of variables and memory for these graphs. For the largest graphs which we can handle using the LP we need more than 1.5 million variables and the input file for CPLEX [39] needs more than 2 GB memory. Moreover we see that the TEN algorithm is much faster (up to 250 times) computed than the solution of the LP. From Figure 3 we also see that the average solution times of both the LP and the TEN grow more slowly for greater networks. We suppose that this is due to the relatively small time horizon and enhanced preprocessing capabilities of CPLEX [39]. Since for the large networks the time horizon is not large enough for all s-d-paths to be used, less flow is shipped through the network and, thus, more variables are zero, which speeds up the computation.

4.2.
3. An algorithm for the quickest flow problem. After having seen how well the concept of the time-expanded network worked, we now show how to use it for the solution of the quickest flow problem. Therefore we add again a super source S and super sink D as in the maximum dynamic flow problem. In this scenario, however, the time at which flow arrives at the sink is important. Hence, we assign turnstile costs [4] to the arcs that lead from any time copy of the sink to the super sink (cf. Figure 4). These costs represent the time at which flow using this arc arrives at the sink. All other arcs get a cost value of zero. Calculating a minimum cost flow with value y s,0 and with respect to the turnstile costs in the TEN then gives a quickest flow, since if flow could arrive earlier in the network this would reduce the costs of the flow in the time-expanded network.
We reduce the problem again to a single source single sink problem by adding super source and super sink. The costs of the arcs which connect these two new vertices to the remaining network are set to zero. Since waiting in the source is also dangerous we have to ensure that all the flow enters the network at time t = 0. This can be done by connecting the super source only to the time copy of the source at t = 0.
Since the sink is considered to be a save vertex this restriction is not necessary for the super sink which is still connected to all time copies of the sink.
To minimize the sum of the costs (23) it is now sufficient to calculate a minimum cost flow with respect to these costs. To minimize the maximum cost value (22) we can iteratively delete the arcs with the highest costs and check if there is still a feasible solution which can send all y s,0 units of flow to the sink.

4.2.5.
Realistic case studies. Finally, we present a numerical example for the coupled model. Therefore we consider the serial computation (see Algorithm 1) of the transportation of the hazardous material, its mapping onto the network and the calculation of the evacuation plan. In this way, the calculation of the dynamics of the hazard can be seen as a preprocessing step for the evacuation. Since we have shown in Section 4.2.2 that the calculation of the dynamic network flow model is much faster than the calculation of the MIP, we only consider the dynamic network flow model and the corresponding objective functions in this example.

Algorithm 1 Solution algorithm
1: Transport of hazardous material as preprocessing 2: Mapping on network 3: Start evacuation solver We start with the description of the network and its parameters. After that, we calculate the spreading of the hazard for different directions and speeds of wind. For each of these settings we then show snapshots of the resulting optimal evacuation process. The source of the hazard is supposed to be located in the middle of the graph (see circle in Figure 5). We consider all vertices at the boundary, which are shown as squares in Figure 5, as safe destinations, whereas the inner vertices of the graph (dots in Figure 5) are considered to be possible evacuation sources. But we restrict ourselves to only one evacuation source. The lines in Figure 5 show  other lines represent arcs in both directions. This is motivated by the fact that it does not make any sense to drive back into the danger zone after having reached safety. As we want to be as realistic as possible the network presented above mimics a region of about 50km 2 in Southern Germany. In order to calculate the capacities and travel times of the arcs as dimensionless numbers we fix the length of one time step θ to be one minute. For each arc we assume a constant average velocity v e depending on the category of the street considered. Depending on the velocity the cars are assumed to hold an appropriate distance d e to the one driving in front. By appropriate we understand approximately the well-known rule of thumb "halftimes speedometer". Furthermore we assume that each car takes three persons to safety (n pc = 3). Thus we may calculate the capacity of each arc e = (i, j) by b ij = v e d e n pc n e,l θ, where n e,l is the number of lanes of the street depending on its category. The travel time τ ij of each arc e = (i, j) is calculated by where L e is the distance between the vertices linked by the arc.
Inserting the values of Table 2 enables us to determine the capacities of the arcs. Since the solver used cannot handle real numbers for the parameters all parameter values -capacities and travel times -are rounded to integers.
After having set up the network we now look at the solutions of the evacuation model which are gained by the time-expanded network algorithm. Since we only want to consider a single source problem, we do not evacuate all vertices. We concentrate on the vertex east of the source of the hazard, which is located nearest to it. Furthermore we give one example for the evacuation of the vertex southwest of the source of the hazard, which is well connected. We consider two directions of wind. On the one hand we consider a wind blowing east cutting the network in a northern and a southern part. On the other hand, we consider wind in southwest direction, because this region is the most densely populated. Additionally, we consider a calm, i.e. only diffusion. The initial area Ω 0 of the cloud is always 4km 2 large and centered at the source of the hazard.
Independently of the direction of the wind we can calculate the quickest flow objective (19) for the evacuation of the vertex east of the source of the hazard (y s,0 = 8787). But from Figure 6 we see that this is -despite of being the fastest solution (cf. Table 3) -not a good solution, especially if the wind is blowing east, since the arcs are mostly covered by the hazardous cloud.  Thus we use the sum objective (23) minimizing the accumulated sum of hazard to which the evacuees are exposed. This is also a reasonable objective function from the medical point of view if an exposure to radiation is considered, since even an exposure for a duration of several hours can be treated as an one-time one [33]. We use the dynamic model (see Section 2.1) to map the cloud onto the network since this is the most realistic one. This can be seen in Figure 7 and in the comparison of the solutions in Figure 8. Furthermore, we observe that using the static model it is possible that the evacuees head in direction of the centre of the hazardous cloud, which is intuitively no good idea. Using the dynamic model such a behaviour of the solution does not occur (see Figure 8 left picture, arc heading north).
We close this section with a comparison of the optimal solutions of the evacuation model with the sum objective for the different wind directions.   We see from Figures 9 and 10 that all arcs leaving the vertex east of the source are usable for the evacuation in the case of wind in southwestern direction and a lull and that the evacuation time is within a 20 % margin of the time of the quickest flow (cf. Table 3). The accumulated concentration for the diffusion case is higher since the center of the hazard does not move away from the vertex, where the evacuees have to wait as the capacities of the streets are not large enough. Thus a lot of people are affected by the hazard before they can be evacuated. For the wind blowing east the evacuation time doubles (cf. Table 3) as at first only the arc heading south is available (see Figure 11 left picture) and only after the cloud has passed the vertex and decayed, the arc heading north can be used. However, the arc heading east can never be used since it leads directly into the hazard.
Outlook. This article presents a first, rather simple, attempt to incorporate the spread of a hazardous material into an egress model for evacuation of regions based on street networks. With the definition of suitable objective functions, the best possible evacuation routes have been found, either to fastest egress or to safest  Table 3. Comparison of evacuation times and accumulated concentration for different wind directions and strengths.
evacuation. Furthermore, we have demonstrated the effectiveness of using dynamic network flows to compute solutions of the original continuous evacuation problem. Future research will focus on extending the current models to multiple sources and to take traffic jams (similar to [6,21]) and selfishness of evacuees into account. This will make the evacuation model nonlinear and requires adapted dynamic flow algorithms. Moreover, the model for the propagation of hazardous material might be more sophisticated. Think of time and space-dependent wind fields or even Navier-Stokes equations. Finally, instead of considering single objective functions, a mulitcriteria model in which two or more objectives are considered simultaneously might allow for adequate compromise solutions and form the basics of a decision support system.