Network Simulation for Pedestrian Flows with HyDEFS

The reliable simulation of pedestrian movement is an essential tool for the security aware design and analysis of buildings and infrastructure. We developed HyDEFS, an event-driven dynamic flow simulation software which is designed to simulate pedestrian movement depending on varying routing decisions of the individual users and varying constraints. HyDEFS uses given density depending velocities to model congestions and evaluates flow distributions with respect to average and maximum travel time. This is of particular importance when considering evacuation scenarios. We apply HyDEFS on two small networks and cross validate its results by time-discrete and time-continuous calculations.


Introduction
Modelling and simulation of pedestrian flows are often applied to the design of traffic infrastructures or the evaluation of rescue routes [1][2][3]. The goal of the investigations is to avoid long lasting congestions or to minimize evacuation times of a building. In order to consider individual pedestrian characteristics such as free walking speed, size, reaction times, etc. mostly agent based models are used. These models calculate locations and speed of individual pedestrians and their interaction with the environment rule-based or using coupled ordinary differential equations [4].

Collective Dynamics 5, A24:1-16 (2020) Licensed under
With respect to route selections these models and simulations use simplifying assumptions. For example, that all pedestrians choose the shortest path [5] or the fastest route [6]. When choosing a route along the shortest path, it is implicitly assumed that all pedestrians have complete knowledge of the possible routes. In real systems, however, this knowledge depends on the users of the building. While in schools and office buildings it can be assumed that the users are familiar with the infrastructure, in subway stations or railway stations, however, this is not necessarily the case. The simplification is also significant when assuming that pedestrians choose the fastest route, since knowledge about the occupancy level would only be given if all routes can be viewed by the pedestrian.
For a more realistic modelling of route selection and knowledge about building structures, first approaches exist which, for example, consider the perception of congestion in front of doors for a rerouting [7], introduce different knowledge levels and uncertainties about distances [8,9], modelling the influence of signage [10][11][12] or incorporating knowledge about the building by modelling cognitive maps, landmarks or generalised knowledge [13]. Although this progress has been achieved, in practice analyses are limited to a few explicit route selection scenarios. For time and cost reasons, it is impossible to consider all possible factors.
The discussion in the previous paragraph considers the system from the pedestrian's point of view. However, this process can also be seen as the flow of people on a network of paths (e.g. escape routes from a building).
For the modelling of flows on networks there are mathematically, physically and engineeringly motivated approaches. Since pedestrians could be considered as directed flows of self-driven particles, the approaches to describe flows on networks for Internet or vehicular traffic are similar. In mathematics mostly PDEs are used to describe the dynamics on the edges [14][15][16][17], in physics and engineering microscopic or agent-based models are used [18][19][20][21][22][23][24][25][26]. In order to model flows at junctions and crossings of a network, the mathematical models introduced in [14,15] rely on assumptions for the distribution of the flows. However, empirical data for the distribution are rare [27,28] and could depend on various factors including properties of the network or its infrastructure, like signage or traffic lights. By using agent-based models, individual route decisions can be taken into account, depending on, besides others, knowledge about shortest paths and psychological factors such as keeping away from walls or following others [18]. For networks of escape routes in buildings the simulations in [7,13] show how different assumptions for the strategy of route choices influence the evacuation time. To consider larger networks, discrete formulations of PDE models are also used for pedestrians [19]. In networks for vehicle traffic, the effects of route decisions on vehicle traffic (traffic lights and traffic jams) [20,21] or network properties such as the betweenness [22,23] were investigated. But also here the models are based on assumptions how the decisions are made.
Considering the complexity of the system already small networks for driven particles show complex phenomena such as the improvement of the evacuation time by reducing the flow at doors [24] or the increase of travel times by introducing an additional route option outlined by the Braess paradox [29,30]. Recent work on the Braess paradox reveals a complex structure of the phase diagram and its dependence on the length ratios of the routes [25] or on assumptions about the route choice of the users [26].
The work presented in [25,26] already studied all possible distributions of flows on a network by using a microscopic model and analysing the phase space by a fine discretization. It precisely maps the dynamics of the flows on the network but requires a high computing effort.
The discussion shows that a realistic evaluation of small networks is already difficult and requires that the models and simulations take into account a high variability of the user's route decisions.
To be able to consider all possible distributions of flows and to parametrize properties of the network or its infrastructure we aim to introduce a simpler approach. First a simulation software HyDEFS based on an event-driven dynamics is introduced. Simulation results of a simple scenario are presented and then these are verified against analytic results for a continuous, as well as a discrete, flow model of the same scenario. Finally, a more complex scenario is studied with the new method.
The approach presented in this article is aimed at allowing to investigate such multiparameter variability also in larger networks. This is also the basis for a future analytic description of the system in order to optimize or to increase the robustness of the network with respect to variations of the distribution of the flows.

The HyDEFS software
The experiments in this work have been done with the HyDEFS (Hybrid Density-dependent Event-driven Flow Simulation) package developed by the authors. The purpose of this software is the simulation of pedestrian flows under varying routing decisions or side constraints. The software allows the simulation of "groups" of persons moving within a network, where the velocity depends on the density, i.e., on the number of persons that are currently on the same edge.
We consider a network G = (V, E) with n nodes and m (directed) edges e = (v, w), allowing to move from node v to node w. Groups of persons may enter the network at some of the nodes, v ∈ V in ⊂ V , and leave it at some other nodes, v ∈ V out ⊂ V .
In our approach, an event ε = (t arr , q arr , v cur , e in ) represents a group of q arr > 0 persons arriving at time t arr at a node v cur , coming from one of v cur 's incoming edges, e in . In particular, the persons entering the network are the events known at the beginning of the simulation, ε in = (t in , q in , v in , 0), where 0 stands for "not coming from an edge." Throughout the simulation, the currently known event points are kept in an event point schedule E , which is organized as a priority queue w.r.t. the events' time component.
If an event point ε = (t arr , q arr , v cur , e in ) corresponds to a group arriving at a "terminal" node v cur ∈ V out and thus leaving the network, then this fact is recorded for the final statistics on, among others, average and maximum time of leaving. Otherwise, the group is distributed over the edges e out starting at node v cur . To this end, each edge of the network has been assigned a probability prob(e), such that { current load (number of persons) on the edge } starting with E = / 0, insert all (t in , q in , v in , 0) corresponding to groups entering the network record that q arr persons have left the network at time t arr else { distribute the group over the edges starting at node v cur } for all outgoing edges e out = (v cur , v next ) such that prob(e out ) > 0 q out = q arr · prob(e out ) increase load(e out ) by q out t next = t arr + len(e out )/v out , where v out = f e out (ρ(e out )) is the current velocity corresponding to density ρ(e out ) = load(e out )/len(e out ) on the edge insert (t next , q out , v next , e out ) into E and a portion prob(e out )·q arr of the group is sent over edge e out . It progresses with a velocity that depends on the (new) density ρ out on the edge, that is, on the "load" of the edge (the number of persons currently moving on it), load(e out ), and on its length, len(e out ). More precisely, the velocity is given by v out = f e out (ρ out ), where ρ out = load(e out )/len(e out ), and f e out is a function that can be prescribed for each edge; see below.
The core of the simulation is summarized in Fig. 1. Note that the term "group of persons" is to be understood in a wider sense, q taking non-integer values in general. Also, subdividing the groups at each node increases the number of event points, but it allows tracing all possible paths and their probabilities with a single simulation, thus obviating the need for a sampling process over single instances, as done in many item-tracking approaches.
Currently, HyDEFS supports the following three types of velocity functions; cf. Fig. 2.
is defined by two limiting velocities, v max and v min , and two thresholds for the density, ρ 1 and ρ 2 . If the current density on the edge is small, ρ ≤ ρ 1 , then the group can move fast, v = v max , whereas high density (ρ > ρ 2 ) leads to slow movement, v = v min . For ρ 1 < ρ ≤ ρ 2 , the velocity v drops linearly from v max to v min . Note that this includes densityindependent velocities (v max = v min ), as well as step functions (ρ 1 = ρ 2 ).
is a "smooth" version of f 1 with v max and v min denoting the asymptotic limits and ρ 1 and ρ 2 indicating where 10 and 90%, resp., of the difference v max − v min are reached.
• f 3 (v max , ρ 1 , ρ 2 ; ρ) describes a velocity that is inversely proportional to the density, with constant value v max for ρ ≤ ρ 1 and dropping to v max /2 at ρ 2 . Note that a length-1 edge with velocity f 3 (1, c, 2 · c) roughly corresponds to a constriction allowing c persons per second to pass, almost independently from the number of persons waiting before the constriction; cf. bottom right picture in Fig. 2.
HyDEFS also includes components for pre-and postprocessing. In particular, the description of the network may contain parameters that can be varied automatically between subsequent runs of the simulation. This feature is used heavily in the following sections.

A simple scenario
We consider a very simple evacuation scenario in a continuous setting. We assume that, over a time interval the door, the exit point is reached by another corridor (10 m from door 1, 20 m from door 2); cf. the sketch in Fig. 3. We assume that, except for the two doors, the corridors are wide enough and there are no obstacles on the way, so that the persons can move with a constant velocity of 1 m/s. For the doors we assume a restricted capacity, allowing c 1 and c 2 persons to pass per second, respectively. This is modelled by the network also shown in Fig. 3, with V in = {1} and V out = {6}, and e 3 and e 5 representing the doors. According to the above assumptions, the remaining edges e 1 , e 2 , e 4 , and e 6 have lengths 10, 20, 10, and 20, resp., and they all have the same (constant) velocity function f 1 (1, 1, 1, 1 ; ρ) ≡ 1. By contrast, e 3 and e 5 have length 1 and the velocity functions f 3 (1, c k , 2 · c k ; ρ). The constant incoming flow is replaced with discrete events, representing groups of size q in arriving at node 1. In this model, the flow can split only at node 1, and we assume that a portion of prob(e 1 ) ∈ [0, 1] follows the outgoing edge e 1 (heading for door 1), leaving prob(e 2 ) = 1 − prob(e 1 ) for the way to door 2.
In the following we report results for this scenario. We first discuss numerical experiments with the HyDEFS software in Sec. 3.1. For the simple scenario these results also can be derived analytically. In Sec. 3.2 we do this with the continuous model, assuming a constant flow entering the network. In Sec. 3.3 a discrete, group-based model is used instead. We observe that the continuous and the discrete model behave rather similarly for this scenario and that the HyDEFS results agree well with both.

Results from simulations with HyDEFS
In all experiments shown in this work, we use T = 10 s, N = 100, q in = 0.1, i.e., the constant flow is replaced with n = 1000 groups entering the network at times 0, 0.01, . . . , 9.99 s, and we fix c 1 = 1.
Since the number of persons arriving per second at node 1 substantially exceeds the capacity of each door, it is clear that both doors should be used even if the way via door 2 is much longer, because the waiting time at door 1 is thus reduced. The simulation results confirm this expectation, indicating that roughly 40% of the persons should use door 2 in order to minimize the maximum and average arrival time at node 6, if both doors have the same capacity c 2 = c 1 = 1; cf. left picture in Fig. 4. We also expect that more persons should be sent to door 2 if the latter's capacity is larger, and again this is confirmed by a simulation with c 2 = 3; cf. right picture in Fig. 4. Indeed, increasing the capacity of door 2 reduces the optimal maximum and average arrival times, as one might hope. The simulation can also be used to optimize the design of the network. Figure 5 shows the maximum and average arrival times if the capacity of door 2 is varied as well. Note that the plots from Fig. 4 correspond to the horizontal lines c 2 = 1 and c 2 = 3 in these images. Going one step further, we can determine how much capacity we must provide for door 2 in order to achieve a desired reduction in the maximum or average arrival time; cf.  Best maximum and average arrival times that can be achieved with a given capacity for door 2, assuming that prob(e 1 ) is chosen optimally in each case.

Analysis for the continuous model
Due to the simplicity of the scenario, the above results can also be obtained analytically.
In the continuous setting, with N persons entering the network during the time interval where the three summands represent, respectively, the start of motion, the time for unhindered motion with maximum speed v ≡ 1 along path k (i.e., len 1 = 10 + 1 + 10 = 21, len 2 = 20 + 1 + 20 = 41), and a potential delay at door k. If the door's capacity can accomodate the arriving flow, then there is no delay: delay k (t) = 0. Otherwise, t · p k N T units of flow have arrived at the door "before me," and t · c k of these have already passed, such that "I" am delayed by Taking both cases together, we obtain t arr k (t) = len k + t max and in particular the latest arrival from path k is at time t arr k (t) = t arr k (T ) = len k + max yielding the overall latest arrival time In our setting (N = 100, T = 10, len 1 = 21, len 2 = 41, c 1 = 1, c 2 ∈ [0.5, 5]), the second term len 1 + T = 31 cannot contribute to the max because the fourth one, len 2 + T = 51, is larger. In fact, the latter can also be dropped, because if it is larger than the first one, len 2 + T > len 1 + p 1 N c 1 , then p 1 < (len 2 − len 1 + T ) c 1 N = 0.3, and in this case len 2 + (1−p 1 )N c 2 > len 2 + 0.7N 5 = 55, i.e., the third term is even larger. As the first (third) term is strictly monotonically increasing (decreasing) w.r.t. p 1 , t max is minimized if In particular, with N = 100, len 2 − len 1 = 20 and c 1 = 1 fixed, for c 2 = 1 the optimum distribution is given by p * 1 = 0.6 with t max = 81, whereas for c 2 = 3 the optimum value t max = 61 is achieved with p * 1 = 0.4. Both agree very well with the simulation results in Fig. 4, as does the dependency "optimum p 1 vs. c 2 " in comparison with the left picture in Fig. 5.
For the average arrival time we use Eq. 1 to obtain In particular, the seemingly quadratic functions in Fig. 4 are only piecewise quadratic on the three intervals 0 ≤ p 1 ≤ c 1 N/T = 0.1, 0.1 ≤ p 1 ≤ 1 − c 2 N/T = 1 − 0.1c 2 , and 1 − 0.1c 2 ≤ p 1 ≤ 1. Again, the minimum of the average arrival time is taken when there is delay at both doors, i.e., t avg = p 1 len 1 + T 2 and setting dt avg /dp 1 = 0 to obtain the minimizer leads again to Eq. 2, i.e., the optimizers for the maximum and average arrival time coincide. These results also agree very well with those from the simulation.

Analysis for the discrete model
The following analysis shows that a discrete, group-based approach can match the behavior of the continuous model very closely. We consider the case c 1 = c 2 = 1 (cf. left picture in Fig. 4). Let again p 1 = prob(e 1 ), len 1 = 21, len 2 = 41, and for a given value of p 1 ∈ [0, 1], let t arr k (i) denote the arrival time at node 6 for that part of the ith group that takes the path via edge e k , k = 1, 2, and i = 1, . . . , n with n = 1000. Then, with similar arguments as in Sec. 3.2 for the delay at the doors, For any (sufficiently large) fixed number of groups, n, we can now evaluate the maximum arrival time t max and the average arrival time t avg at node 6. In the first case, we obtain t max = max{t arr 1 (n),t arr 2 (n)}. The minimum of t max =: t max (p * 1 ) w.r.t. p 1 is attained when t arr 1 (n) = t arr 2 (n), i.e., when For n = 1000 this implies that the best possible distribution is obtained for p * 1 ≈ 0.6001, achieving the smallest possible maximum arrival time of t max (p * 1 ) = 80.95. Assuming again a (sufficiently large) fixed number of groups, n, and setting s n = ∑ n i=1 i−1 10 = n(n−1) 20 , we can evaluate the average arrival time t avg = t avg (p 1 ) at node 6 for 0.1 ≤ p 1 ≤ 0.9 as = 1 n 2 s n p 2 1 − (20 n + 2 s n ) p 1 + len 2 n + s n .
Differentiating w.r.t. p 1 yields dt avg (p 1 ) dp 1 = 1 n (4 s n p − 20 n − 2 s n ) Thus, for n = 1000 we obtain p * 1 ≈ 0.6001 and t avg (p * 1 ) ≈ 54.97. Note that the calculations above are only meaningful if n > 200, i.e., if N > 20. Moreover, we can observe that the optimal distributions p * 1 for minimizing t max and t avg are equal for all reasonable values of n.

A more complex scenario
The situation changes considerably if, before leaving the scene, all persons have to pass another corridor, which is susceptible to congestion; cf. the sketch in Fig. 7. The additional corridor is modelled by a new edge e 7 of length 10 leading to the new exit node 7; cf. Fig. 7, and e 7 is equipped with the velocity function f 3 (1, 3, 5 ; ρ), i.e., it allows movement with maximum velocity v = 1 for densities up to 3, but for ρ = 5 the velocity drops by one half. Thus, the number of persons leaving the corridor per time unit is decreasing for higher densities. To allow for a more localized behavior (the velocity of persons at one end of the edge needs not depend on the density at the opposite end), the edge is subdivided into five sub-edges of length 2, all sharing the above velocity function. These sub-edges will be referred to as e 7A (adjacent to node 6) to e 7E in the following.
The maximum and average arrival times for this scenario are shown in Fig. 8. The simulation indicates that-in contrast to Sec. 3-increasing the capacity of door 2 beyond a certain bound may lead to higher (maximum and average) arrival times. That is, even if the persons can reach node 6 earlier, this leads to a higher density in the final corridor, which in turn will slow down the movement there.
It is noteworthy that in this scenario, the approximate optimizers for the maximum arrival time (c 2 = 2.81, p 1 = 0.218, with t max ≈ 78.72) and for the average arrival time (c 2 = 1.86, p 1 = 0.470, with t avg ≈ 59.87) are different. They correspond to the bottom left and top right, resp., blue regions in the left picture of Fig. 8. Choosing the "wrong" optimizer does, however, not too much harm in this case: At the max-optimizer, we have t avg ≈ 59.91, whereas t max ≈ 79.33 at the avg-optimizer.
More details of the flow can be seen in Fig. 9, which shows the density on selected edges over time for three different configurations from a vertical line through the maxoptimizer, i.e., p 1 = 0.218, with c 2 being smaller (top left), equal to (top right), and larger than (bottom), resp., the optimal value c 2 = 2.81. In all three cases, the flow arrives at edge e 4 after len(e 1 ) + len(e 3 ) = 11 seconds, and since a total of p 1 N = 21.8 persons come through door 1, it continues for 21.8 seconds, and the edge becomes empty after another len(e 4 ) = 10 seconds. A similar behavior can be observed for e 6 , the specific times depending on the respective capacity of door 2. On e 7 , the different behavior shows mainly at the beginning of the edge, e 7A . The other sub-edges feature a a very similar flow, with a delay of len(e 7x ) = 2 seconds from one sub-edge to the next. In the case c 2 = 3 too much flow arrives at e 7A , leading to reduced velocity on this sub-edge. The flow leaving the sub-edge is small enough such that the later sub-edges do not suffer from the same problem.

Conclusions and future work
Routing decisions of users significantly influence the effectiveness of a pedestrian infrastructure. However, the modelling of path selection is due to uncertainty and can even in evacuation scenarios not be completely controlled. HyDEFS simulates pedestrian movement with a macroscopic event-driven model with respect to varying routing decisions and constraints. We consider in this article rather simple scenarios such that the simulation results can be verified analytically. On the small network (in Section 3) the discrete and the time-continuous calculation yield almost the same optimal solutions (minimizing t avg and t max , respectively) as the simulation with HyDEFS. Interestingly, in the first scenario (Section 3) with c 1 = c 2 = 1, the optimal distributions of flow (with respect to average t avg and maximal arrival time t max ) converge to p * = 1 2 , when the number of groups n goes to infinity (with N = n/10 and T = n/100). This effect can be observed analogously in the discrete and the continuous calculation, as well as in the simulation. In the limit n → ∞ this effect relates to the transition of a time dynamic  Density of the edges e 4 , e 6 leading to node 6 and of the sub-edges e 7A−E leading from node 6 to node 7 for the optimal value p 1 = 0.218 and different values of c 2 : c 2 = 2.00 (smaller than optimal, top left), c 2 = 2.81 (optimal, top right), c 2 = 3.00 (larger than optimal, bottom; note the different density range).
flow to a steady-state maximum flow problem, whose optimal solution is determined by the capacities of the paths (in this case the doors) and is independent of the path lengths.
In the more complex scenario (Section 4) the optimal distribution of pedestrians with respect to average and maximal arrival time do no longer coincide. In larger networks we expect that these optimal solutions will differ significantly, which asks for a multiobjective analysis of network designs. Additionally, other objective functions, like e.g., earliest arrival flow [31] and other functional dependencies of velocity and density should be evaluated on larger networks in future research using HyDEFS.
Compared to the analyses of the phase diagram of the network presented in [25,26], the approach introduced in this work represents a simplification. This is associated with a lower computational effort but also with a less accurate representation of the dynamics of the flows on the network. Whether and for which systems this simplification reflects the essential dynamics has to be investigated in future work. In addition, it could be investigated whether a coupling of models with different levels of detail allows to roughly model the main parts of the network and to limit complex calculations with fine models to only a few edges and nodes of the network.