Evacuation time estimate for a total pedestrian evacuation using queuing network model and volunteered geographic information

Estimating city evacuation time is a non-trivial problem due to the interaction between thousands of individual agents, giving rise to various collective phenomena, such as bottleneck formation, intermittent flow and stop-and-go waves. We present a mean field approach to draw relationships between road network spatial attributes, number of evacuees and resultant evacuation time estimate (ETE). We divide $50$ medium sized UK cities into a total of $697$ catchment areas which we define as an area where all agents share the same nearest exit node. In these catchment areas, 90% of agents are within $5.4$ km of their designated exit node. We establish a characteristic flow rate from catchment area attributes (population, distance to exit node and exit node width) and a mean flow rate in free-flow regime by simulating total evacuations using an agent based `queuing network' model. We use these variables to determine a relationship between catchment area attributes and resultant ETE. This relationship could enable emergency planners to make rapid appraisal of evacuation strategies and help support decisions in the run up to a crisis.


I. INTRODUCTION
Interaction between individual agents, city topology, disaster type, evacuation mode, information propagation patterns and stochastic variables can all influence the temporal extent of a city-wide evacuation.Additionally, growing urban populations [1] amplifies the impact of extreme events [2].There is a need to examine factors affecting evacuation time in relation to the latest understanding of crowd dynamics and evacuation behaviour.
Evacuation time estimate (ETE) analysis (a) tells emergency planners if an evacuation plan can reduce hazard exposure time (b) measures effect of uncontrollable events such as adverse weather and (c) assesses whether traffic management actions help reduce it [3].A study of flood evacuation in Netherlands identifies a need for alternative evacuation strategies for coastal areas after it found that it was not feasible to evacuate preventively within a 48 hour warning window [4].EMBLEM2, an empirical study, categorises research findings about evacuees' behaviour in hurricanes into 4 evacuation route system parameters, 16 behavioural parameters and 5 evacuation scope/timing parameters to calculate ETE [5].A sensitivity analysis of radiological emergency micro traffic simulation finds that ETE is sensitive to traffic factors (interaction with pedestrians, intersection traversing time, car ownership, etc.) and route choice mechanisms (shortest path and myopic behaviour) [6].NETVACl, a macro traffic simulation finds that ETE for areas surrounding nuclear power plant sites are sensitive to road network topology, intersection design and control, and a wide array of evacuation management strategies [7].Another study produces ETE for 10 miles radius of 52 nuclear power stations taking consideration of factors such as population density, weather conditions, warning time, response time and confirmation time [8].
Some models take a dynamic network flow approach to minimise evacuation time [9,10] while others use social force based models like EPES to establish optimal earthquake evacuation behaviour [11].The 'Last-Mile' project uses a 'queuing network' model to obtain an optimal evacuation plan for the Indonesian city of Padang using time-dependent network attributes to imitate conditions of a tsunami [12].The underlying flow model simulates traffic taking only free speed, bottleneck capacities and space constraints into account.This approach is preceded by an early evacuation plan optimisation study for Yokosuda city in Japan which uses a combination of the shortest path algorithm and minimal cost flow approach accounting for the capacity limit of each place of refuge [13].
Evacuees' behaviour plays an import role during evacuations.Combination of individual traits and basic social psychological processes such as (a) risk perception, (b) social influence and (c) access to resources predict evacuation behaviour while some population subgroups choose not to evacuate depending on the severity of storm, territoriality, etc. according to a study conducted after Hurricanes Hugo and Andrew [14].Subjective perception of how bad the storm is going to be and the severity of damage also seem to play an important role in evacuation likelihood following a warning [15].The effect of compliance behaviour on ETE has been studied using the EVAQ evacuation model and a case study of the Rotterdam metropolitan area in Netherlands [16].
Crowd dynamics is an important feature in large cities and understanding it is a crucial component of emergency evacuation modelling where Agent Based Mod-elling (ABM) is increasingly being used for large scale simulations to account for many interacting entities [17].The transition between low and high density phases are common in social systems like cities [18].Keeping a constant lower limit on the net-time headway has been shown as one of the key mechanisms behind emergent crowd dynamics [19].Observed collective phenomena in pedestrian crowds include lane formation in corridors and oscillations at bottlenecks in normal situations as well as different kinds of blocked states produced in panic situations [20].Video recordings of the crowd disaster in Mina/Makkah during the Hajj on January 12, 2006 reveal two subsequent, sudden transitions from laminar to stopand-go [21] and turbulent flows [22].The transition to turbulent flow is responsible for sudden eruptions of pressure release comparable to earthquakes, which cause sudden displacements and the falling and trampling of people [23].However, from a macroscopic viewpoint, pedestrian behaviour can be assimilated into a relationship between walking speed v and local density k, variables familiar to the transport research community [24].
Review of existing work highlights a gap in understanding which relates ETE to interaction between city population and their topological make-up.Topologies can vary between parts of cities, one city to another and one region of the planet to another, all growing in complexity at the same time.A 'queuing network' ABM which incorporates pedestrian behaviour and network topology has the potential to define a direct relationship between city topological attributes and their ETE.

II. METHODOLOGY
We will now describe a model used for deriving the necessary quantities required for our analysis.We make the following assumptions across the model: • Evacuation type is a total evacuation scenario, for which exit nodes lie at intersection between major roads and city administrative boundary [25].
• Evacuation mode is by walking only.
• Route to exit node is calculated using Dijkstra's shortest path algorithm [26] (no dynamic routing due to bottlenecks congestion).
We incorporate Weidmann's fundamental diagram to describe pedestrian behaviour [24] shown in FIG. 1 into the model.Eq. ( 1) describes the relationship between density k and velocity v.When k = 0 ped/m 2 , the freeflow velocity v f is 1.34 m/s.At maximum density when The relationship between density k and flow rate Q follows as Eq. ( 2).We can differentiate this equation to  [24] where the relationship between density k and flow ).In this equation, the optimum density kopt = 1.75 ped/m 2 and the corresponding maximum flow Qmax = 1.22 ped/ms.derive the optimum density k opt = 1.75 ped/m 2 when dQ/dk = 0.The corresponding maximum flow rate Q max = 1.22 ped/ms.
As pointed out in the introduction, the fundamental diagram is not able to describe system dynamics far from equilibrium (i.e.high density crowds).Therefore, to ensure that agent movement can occur at high densities, we implement a network link density cap at k cap = 5 ped/m 2 such that minimum velocity is never less than v min = 0.04 m/s.This ensures a minimum flow of Q min = 0.05 ped/(ms) [19].Without this cap, as k approaches k max , v tends towards 0 m/s, which means that the simulations would run indefinitely.
Under these assumptions, we select 50 cities similar in area to City of Bristol (235.82 ± 25% km 2 ).We use network topology approximated from OpenStreetMap (OSM) [27].OSM is a source of Volunteered Geographic Information (VGI) [28] growing in both contributor base and data quality [29][30][31].We further divide these cities into 697 catchment areas (CA), which we define as network components that emerge as agents are assigned to an exit node nearest to their initial position calculated using Dijkstra's shortest path algorithm [26].FIG.2a illustrates CA formation for City of Bristol and FIG.2b shows how the distribution of initial agent distance to their exit node D varies between different CAs.

A. Characteristic Variables
Characteristic variables independent of dynamic agent interaction informs part of our analysis.The first of these is the characteristic flow rate Q c described by Eq. ( 3).It is defined as free-flow time averaged flow whereby we assume infinite link capacity.To illustrate the point, we mark the position of Q c for an example CA in FIG.4a.
We calculate Q c using CA population N , exit node width W and free-flow catchment area traversal time for 90% of all CA agents T 90% f , which is also our second characteristic variable.We estimate N from GRUMPv1 year 2000 population dataset [32] uniformly scaled up by a factor of 9.37% in order to account for the rise in UK population between the years 2000 and 2015 [33].It has a granularity of 1 km 2 .FIG. 3a shows how N is distributed in log scale.Values range from 10 2 to 10 5 , peaking at 10 3.82 ≈ 6, 628 agents.For exit nodes tagged 'motorway' on OSM, we assume W = 7.5 m and for those tagged 'trunk' or 'primary', W = 5 m [25].
where D 90% is the distance to exit node for 90% of all agents and free-flow velocity v f = 1.34 m/s.T 90% f is also marked in FIG.4a.If we ignore all congestion and bottleneck effects, it provides a lower bound estimate of evacuation time for 90% of all CA agents.We use D 90% because it approximates the size of a CA as a scalar without the weight of the last decile skewing the result.FIG.3b shows how D 90% is distributed across all CA in log scale with values ranging between 10 1 to 10 5 metres.It peaks at 10 3.74 ≈ 5, 435 metres, a distance belt within which 90% of all CA agents are situated.

B. Simulated Variables
Simulated variables are obtained by studying dynamic interaction of agents under the following assumptions: • Agents immediately 'walk' to the nearest exit on a signal to evacuate (i.e.pre-movement time is zero), • Agents act independently (complex social behaviours such as family regrouping, co-operation, etc. are not taken into account).
Our first simulated variable Q f is defined as simulated exit node flow rate Q averaged within the free-flow regime (T < T 90% f ).The area under the flow curve for each CA is proportional to the total number of agents.Larger this area before flow transitions to congested phase, bigger the Q f value, precipitating a shorter congestion.Hence, the overall ETE is proportional to Q f .We show the position of Q f for an example CA in FIG.4b.The flow curve Q it is derived from is calculated using Eq. ( 2) where the density parameter k = N/(W L), N is the number of agents arriving at the exit node per time-step and W L is the area of the exit link.FIG.4b also marks the position of T 90% , defined as the time at which 90% of all CA agents arrive at the exit node.Unlike T 90% f , T 90% takes agent interaction and emergent bottlenecks into account.While bottlenecks may be interspersed throughout a CA as shown by the example in FIG. 5

III. LINKING CHARACTERISTIC AND SIMULATED VARIABLES
In this section, we establish the link between characteristic variable (Q c , T 90% f ) and simulated variables (Q f , T 90% ).We aggregate the simulated exit node flow Q observed through absolute simulation time T .Then, we level the basis for comparison between CAs by normalising flow as Q/Q c and time as T /T 90% f .We substitute Q and T for Q f and T 90% and define Q f in relation to Q c .We also define ratio Using these relationships, we derive a general description of T 90% using characteristic variables Q c and T 90% f .Aggregating simulated flow at exit node Q across all CA over absolute simulation time T produces FIG.6a.Looking at 0 < T < 20000 band, we observe that the aggregate flows peak around Q ≈ 0.15 ped/(ms) within a wide 68% confidence interval early on in the simulation which gradually tapers.While the peak signals the transition from free-flow (T ≤ T 90% f ) to congested (T > T 90% f ) regime, the exact point of transition is not clear in this representation.We also observe that as the sample size decreases with elapsing T , there is an increase in fluctuation of aggregate Q.We normalise T by 6 which implies that in general, Q c over-predicts the simulated flow.The flattening of the curve beyond the peak at T /T 90% f ≥ 1 indicates the congested flow regime which carries on up to a maximum of T /T 90% f ≈ 72.This is a significant gap between free-flow and simulated time but only applies to a small number of CA.
For the following part, we randomly divide our 697 CA into two datasets, the first half (the 'training' dataset) to train our model with containing 347 CA and the second half (the 'testing' dataset) to test our model containing 348 CA.
Using the 'training' dataset, we attempt to understand how Q c over-predicts simulated flow Q f .FIG. 7a shows this relationship.The upper bound appears to be defined by It is better defined by a power-law fit (r 2 = 0.79) described by Eq. ( 4) where θ = 0.73 and γ = 1.12.
We look for an equation to estimate the ETE, i.e.T 90% by analysing the relationship between ratios Q f /Q c and T 90% /T 90% f representing Q and T as Q f and T 90% respectively.Q f /Q c estimates the peak of the mean curve in FIG.6b.For T 90% /T 90% f 1, delays due to agent interaction is proportionately greater and as such, T 90% /T 90% f = 1 is the best possible desired outcome.We use the 'training' dataset to derive the relationship seen in FIG.7b between Q f /Q c and T 90% /T 90% f with axes.There is a strong correlation (r 2 = 0.81) between the data points.For values of . The relationship between the two ratios is well described by the power-law of Eq. ( 5) with best fit parameter values α = −1.44 and β = 0.92.
In order to obtain at least the first order estimate of ETE for a new CA without running an ABM simulation,  .The best fit power-law equation is we can equate T 90% solely in terms of characteristic variables T 90% f and Q c .We do this by substituting Eq. ( 4) into Eq.( 5) to obtain Eq. ( 6) where φ = α(θ − 1) = 0.38 and ω = βγ α = 0.78.
We verify Eq. ( 6) using the 'testing' dataset.We use Q c and T 90% f parameters alone to calculate T 90% for each CA in this dataset and compare them against their simulated counterpart.According to FIG. 8, there is a good agreement between the values (r 2 = 0.73) where T 90% simulated = η(T 90% calculated ) ζ .The calculated values under-estimate the simulated values (exponent ζ = 1.16, coefficient η = 0.44) for higher values of T 90% , as shown by the deviation of the trend from mirror diagonal line (exponent ζ = 1.00, coefficient η = 1.00).

IV. CONCLUSIONS
In conclusion, by exploring the underlying relationship between simulated ETE and 697 CA attributes from 50 UK cities, we present a method for calculating ETE, all using CA attributes: population, size and exit node width alone.This method more reliably estimates the ETE when characteristic flow rate is similar to mean of simulated exit node flow rate in the free-flow regime.There are discrepancies which exist between calculated and simulated ETE because statistical analyses do not fully capture the unique attributes of each CA.Hence, ABMs are better placed to deal with problems with many interacting entities.In our future work, we want to search for topological attributes that uniquely describe each CA which explain these discrepancies.

FIG. 1 .
FIG. 1.(a) Pedestrian density-velocity diagram [24] where the relationship between density k and velocity v = v f (1.0 − e −1.913( 1.0 k − 1.0 kmax ) ).In this equation, the free-flow velocity v f = 1.34 m/s.At maximum density kmax = 5.4 ped/m 2 , v = 0 m/s.(b) Pedestrian density-flow diagram [24] where the relationship between density k and flow Q(k) = kv(k) = kv f (1.0 − e FIG. 2. (Colour Online) (a) Catchment areas (CA) are obtained by allocating agents to the exit node nearest to their initial position.For this, we use Dijkstra's shortest path algorithm [26].Each colour in the figure represents one of the 15 City of Bristol CA.(b) Example showing distribution of agent distance to exit node D for 3 City of Bristol CAs denoted by red (CA01), green (CA02) and blue (CA03) histograms.

FIG. 3 .
FIG. 3. (a) Histogram of catchment area (CA) populations N .(b) Histogram of agent distances to exits for 90% of all CA agents D 90% .

FIG. 4 .
FIG.4.An example to show the position of free-flow time for 90% of all catchment area (CA) agents T 90% f and simulated time for 90% of all CA agents T 90% using a City of Bristol CA (CA02).It also marks the position of maximum flow rate Qmax, characteristic free-flow time averaged flow rate Qc, mean of simulated exit node flow rate in free-flow regime Q f .
where observed local density k varies through distance from exit node D and elapsed time T , it is ultimately the exit node flow rate Q that influences the overall ETE.FIG. 5 also illustrates how velocity drops where density is high for a randomly picked agent trajectory.

2 ]FIG. 5 .
FIG. 5. Density k at distance D away from exit node at time T where the density ranges between 0 ≤ k ≤ 5 ped/m 2 for a City of Bristol CA (CA02).Trajectory of a randomly picked agent is shown to illustrate how the agent velocity is reduced where the link density is high.

FIG. 6 .
FIG. 6.(a) Aggregate simulated flow at exit node Q averaged across all CA over absolute simulation time T where the grey region signifies 68% confidence interval which becomes narrower with decreasing amount of aggregate data sample.(b) Aggregate exit node flow rate normalised by characteristic flow rate Q/Qc over time normalised by free-flow time for 90% of CA agents T /T 90% f aggregated from all CA showing the 68% confidence interval.Q/Qc peaks within T /T 90% f < 1 and mean Q/Qc < 1.
FIG. 7. (a) Relationship between characteristic flow Qc and mean of simulated exit node flow rate in free-flow regime Q f .Each CA is represented by a data point.Q f = Qc is the upper bound.Q f = γ(Qc) θ describes the power-law best fit.(b) Relationship between the ratio of mean simulated exit node flow rate in free-flow regime to characteristic flow rate Q f /Qc to the ratio of simulated time to free-flow time for 90% of all CA agents T 90% /T 90% f