Simulating competitive egress of noncircular pedestrians

R. C. Hidalgo,1,* D. R. Parisi,2,3,† and I. Zuriguel1,‡ 1Departamento de Física, Facultad de Ciencias, Universidad de Navarra, E-31080 Pamplona, Spain 2Instituto Tecnológico de Buenos Aires, Lavarden 389, (C1437FBG) C. A. de Buenos Aires, Argentina 3Consejo Nacional de Investigaciones Científicas y Técnicas, Godoy Cruz 2290 (C1425FQB), C. A. de Buenos Aires, Argentina (Received 13 December 2016; published 25 April 2017)


I. INTRODUCTION
Crowds in competitive behavior trying to cross a geometrical constriction can be observed in many circumstances.As an example, this happened in emergency egress (such as in The Station nightclub, USA [1,2] and Hillsborough Stadium, England), in aggressive shopping (black Friday stampedes [3][4][5][6]), and in everyday situations within transport systems [7].Undoubtedly, understanding and simulating this kind of competitive flow in confined geometries would improve the design of pedestrian facilities, achieving greater safety performance.
In the last two decades, an important number of pedestrian models have been developed and used to study competitive egress.Among these models, we can distinguish two approaches: cellular automata and force-based models.Cellular automata operate in discrete space and time, and the states of cells change following predefined rules.These models were introduced by Kirchner and Schadschneider [8] and have been broadly accepted, primarily because they are very fast in computational terms and they allow simulating a large amount of individuals.The force-based models (also known as "second-order" or Newtonian models) however, seem to provide a more natural and accurate description of the dynamics of the system as they consider continuous space and time allowing one to compute contact forces and pressures.In this sense, the work of Helbing et al. [9] investigating egress dynamics using the social force model represents a turning point.
Despite the large amount of models developed, the lack of experimental studies on competitive egress under controlled conditions was thought to be a drawback that prevents model validation.Very recently, experiments of competitive-human egress have been reported [10,11] confirming the existence of the faster-is-slower effect predicted by the Helbing-Farkas-Vicsek model [9].These new data open a singular opportunity for calibrating or developing models whose outputs approach those of the experiments.In particular, a key observable of these competitive egress drills was the statistics of time lapses between two consecutive outgoing pedestrians (δt).Its distribution reveals a power-law tail, in accordance to other systems of bottleneck flow such as grains, colloids, and sheep [12,13].
This feature can be displayed by a version of the social force model implemented in [13] but for unrealistic values of δt (δt ∈ [∼ 10 2 -10 6 ] s, while the experimental values are in the range, δt 10 1 s).Beyond this work, force-based models showing the power-law tail of the δt distribution are not found in the literature (note that some of these models perhaps are able to reproduce this feature but they have not looked at it).Very recently, a cellular automaton model was introduced that is able to reproduce, at least semiquantitatively, the experimental power-law distributions [14].
In order to fill the existing gap, we propose a force-based model considering noncircular particles, more specifically, spherocylinders [15].Note that the shape of the particles has been shown to play an important role in dense granular flows [16][17][18], and in particular in granular flows through bottlenecks [19][20][21].Thus, it is expected that in competitive pedestrian dynamics, where the "granular" or contact component is important, this parameter is also relevant and its implementation could improve the model performance.
Almost all pedestrian models considered circular pedestrian shape, however, there exist a few antecedents of pedestrian models with noncircular shape.Thompson and Marchant [22] proposed to simulate pedestrians by superposing three aligned circles ["three circle model"(TCM)].This configuration facilitates the calculation of interpersonal distance, which is used for computing the corresponding velocity of the simulated pedestrian.However, the model is not force based, and neither contact nor frictional forces are calculated.The same considerations hold for the simulations of elliptical pedestrians flowing within a crowd presented in [23].A force-based model that considers a repulsive force among ellipses was introduced by Chraibi et al. [24].Indeed, the proposed force was a function of the border-to-border distance of the ellipses along the line connecting both simulated pedestrians.Importantly, overlapping pedestrians was not desirable in this model and in consequence, tangential friction was not included.Despite the undoubtable novelty of these works, we must note that they were implemented in a noncompetitive regime where it is expected that contact and friction are not as relevant as in the competitive scenario.In the only antecedent we are aware of, pedestrians are simulated with shapes different to the disk in competitive conditions in the work of Alonso-Marroquín et al. [25].They describe the shape of a human body [its two-dimensional (2D) transversal section at chest height] using spheropolygons and considered the normal and tangential contact forces between pedestrians.Nevertheless, the proposed scenario was complex and experimental validation was clearly unfeasible.
In this work, we implement a force-based model where the pedestrians projection in the 2D plane is simulated by spherocylinders.Importantly, the simulated scenario is a single room evacuation process in highly competitive conditions where the contacts among pedestrians have been shown to be persistent and where it is expected that torques among the anisotropic bodies become important in determining the global behavior.In Sec.II we will explain the model.Then, in Sec.III we will globally describe the evacuation behavior, and the statistics of exiting times.Finally we will analyze the effect of three different variables of the model as (i) the desired pedestrian speed; (ii) the strength with which pedestrians tend to orient in the direction perpendicular to the desired movement direction; and (iii) the noise that allows pedestrians rotation to escape from a jammed situation.

II. THE MODEL
We adopted a hybrid CPU-GPU discrete element algorithm [26,27] which benefits from the highly parallel structure of the NVIDA graphics processing units (GPUs).For this reason GPUs become very effective when processing large blocks of data and have emerged as a serious alternative for parallel computing.Our algorithm has already been successfully used in describing the behavior of spherical and nonspherical inert grains [28][29][30].
The approach consists in numerically simulating the competitive egress of pedestrians through a narrow door, using nonspherical particles and taking into account the correlations between the particle's shape and the driving forces and torques.The proposed model is based on three main ingredients.
First, we use a better approximation for the pedestrian's shape than the traditionally used disk.Indeed, pedestrians are modeled as spherocylinders and thus, they are characterized by their length l and spheroradius r.Therefore, the aspect ratio reads as ζ = (l + 2r)/2r.Pedestrians are considered as self-propelled particles with a driving mechanism which was first proposed for vehicular traffic by Pipes [31] and used for pedestrian simulations by Helbing et al. [9,32].The driving force F Di reads as where m i is the pedestrians' mass, v i is its velocity, and v d = v d e i accounts for the desired velocity.As usual, the normalized vector e i points to the desired target.Here, the target is the closest point of the door aperture, which is defined as a segment whose length is slightly shorter than it [see Fig. 1(b)].Moreover, τ represents the characteristic time it will take for a pedestrian to achieve v d if there were no interactions with other pedestrians.
As a second ingredient we propose that the pedestrians have a desire to advance with the body aligned perpendicularly to the movement direction.Note that this is necessary if we want to mimic the pedestrian alignments observed in evacuation drills as shown in Fig. 2. Otherwise ellipses tend to get aligned with their long axis in the direction of the flow as reported for the case of inert grains [21].In order to implement this preferential alignment we define a unitary vector perpendicular to the particle's longer side n e , and compare its orientation with the acting direction e i of the driving force, F Di [see Fig. 1(b)].Therefore, the desire to turn to the prescribed direction is given by the existence of a driving torque τ D , which acts aligning FIG. 2. Snapshot corresponding to an evacuation drill in highly competitive conditions as reported in [10,11].Spherocylinders are superimposed on the image by hand with an orientation that is, approximately, that of the pedestrians.Clearly, most of pedestrians align perpendicular to the preferred movement direction but, near the door, high shear forces lead to rotations and emergence of disorder.
the vector n e with the direction of the target e i .This driving torque reads as (2) Note that this driving torque has a linear part that depends on the angular deviation of n e from the target ( θ ) and S D is the angular strength.Additionally, it also includes a viscous term, which reflects a dependence on the angular velocity θ through the damping coefficient β.In all cases presented here, we have used β = 4.5 √ S D , which guarantees overdamped conditions.
Finally, the third ingredient in the model is the introduction of a sinusoidal noise R(t), which mimics the pedestrian desire to rotate the body in any direction when it gets trapped in a given position (with a given orientation).R(t) is characterized by the amplitude η and the period T , which is set to T = 1.0 s, in all cases.The introduction of this periodicity, which makes this "noise" of different nature than previously used ones [13,32,33], was inspired by visual observation of evacuation drills where jammed pedestrians seem to perform rotational movements during a time lapse of the order of 1 s.
As stated before, pedestrians move at predetermined desired velocities toward the exit.But on their way out, they cannot occupy the same physical space as other pedestrians or walls.In the model, this is implemented by a contact force akin to the granular force.In our system the granular force F Gi acting on particle i [see Fig. 1(a)] is computed as where F wi is the contact force of the simulated pedestrian with the wall, and the sum accounts for the forces F ij exerted by the N c contacting neighbors j on i.More details about the numerical implementation of the granular interaction among two spherocylinders can be found in the Appendix.
In discrete element modeling, the translational motion of each simulated pedestrian is governed by Newton's second law of motion: and the rotational motion is ruled by the Euler equations: where r ij is a vector pointing from the center of particle i to the contact point with particle j .We have implemented a velocity Verlet numerical algorithm to integrate both the translational and rotational equations of motion of each simulated pedestrian.With all these ingredients, we simulate an evacuation process where 192 persons try to leave a room through a small door of width W = 0.70 m.We consider a square room of 8 m × 8 m.Each pedestrian is modeled as a spherocylinder with the long axis uniformly distributed in the range [0.35-0.5 m], the short axis in [0.24-0.33 m] and a mean mass of 67 kg with a normal distribution truncated at [45-114 kg].These values correspond to the real experimental drills performed by Garcimartín et al. [11].Initially, the pedestrians are homogeneously distributed within the room, facing the door, at the vertices of a 12 × 16 square lattice (16 individuals arranged in the direction perpendicular to the door).Then, they begin to move toward the exit, in such a way that the direction of the desired velocity e i points to the closest point of the door [Fig.1(b)].In order to avoid transient effects [34] related to the system's initial state and the reduction of pedestrians in the room, we simulated steady-state situations.To this end, we have implemented periodic boundary conditions in the direction perpendicular to the door as was previously done to simulate room evacuation [35] (in the direction parallel to the door we also implemented periodic boundary conditions but the room was wide enough to avoid pedestrians going from one side to the other).Then, once a pedestrian passes through the exit, it is reinserted at the back of the room (at 8 m from the door).Following this procedure steady flow conditions are quickly reached and the simulated pedestrians concentrate behind the door forming a dense cluster as shown in Fig. 3.The mean size of the cluster fluctuates around N c ∼ 180 pedestrians.
As is shown below, our model reproduces that the acting direction e i of the driving force F Di is correlated with the particle's shape.Thus, as the angular strength S D is increased the pedestrians get more aligned with their longer side perpendicular to e i , which points to the target.This fact has a significant impact on the particle's flow because it determines the particle's effective size.Moreover, we find that the amplitude of the sinusoidal noise η is also important as it determines the ability of the pedestrians to escape from a jammed configuration.For this reason, we have systematically studied the influence of these two parameters (S D and η) on the flow properties for different competitiveness levels.The latter are simulated by changing the desired velocity of the evacuating pedestrians, but similar effects could be achieved by modifying the relaxation time τ of the driving force [see Eq. ( 1)].

III. RESULTS
As we mentioned above, we examine the evacuation process of 192 persons, which leave the room through a small door, with periodic boundary conditions both in the direction parallel to the door and perpendicular to it.As pedestrians that exit the room are reinserted from the opposite side, a steadystate condition quickly establishes, allowing us to perform a long simulation instead of the traditional procedure of testing several egresses.We perform numerical simulations of T end = 1600 s duration for each set of parameters.During this time, the number of evacuated pedestrians ranges from 2000 to 4000, depending on the simulation parameters.

A. Discharge curves
In Fig. 4(a) we report the number of egressed pedestrians as a function of time (which is known as the discharge curve) for three desired velocities.The driving torque is ruled by S D = 20 N m and η = 1 in units of mgL, where g = 10 m/s 2 and m and L are the largest pedestrian's mass and long axis, respectively.It can be seen that lower evacuation flow rates (lower y values in the graph) are obtained as the desired velocity increases which is compatible with the faster-isslower (FIS) effect predicted in [9] and proved experimentally in [10].Interestingly, very different behavior is observed at the very beginning of the simulation, before the stationary state is reached [inset of Fig. 4(a)].Clearly, when t 8 s the flow rate is higher for greater desired speeds, in agreement with recent experimental findings [11].This behavior has been related to the short transient that takes place when pedestrians first approach the door and found enough free space to reach its desired velocity.After this, agglomeration behind the door occurs and the stationary regime is established.
With the aim of looking to the temporal evolution of the flow rate, we define the instantaneous specific flow (Q I ) which is the flow rate obtained during the passage of 200 consecutive pedestrians: where N tot is the number of pedestrians evacuated in 1600 s, the subindex i corresponds to each pedestrian, and t i is the time at which the pedestrian i escaped from the room.In Fig. 4(b) we display the evolution of Q I over time, confirming that the steady-state conditions are reached and the flow fluctuates around constant values for each desired velocity.In addition, it is evidenced that less competitive conditions (v d = 1 m/s) lead to smoother discharge curves with smaller fluctuation amplitude than in the high competitive scenario (v d = 3 m/s).Moreover, in the latter case, long flow interruptions (clogs) are observed.Evidently, the way in which these flow interruptions appear in a plot such as Fig. 4(b) strongly depends on the window chosen to calculate the instantaneous flow rate: the shorter the window the higher the number of times that Q I takes values near or equal to zero.Therefore, we perform an alternative characterization of these flow interruptions by looking at the time lapses between successive outgoing individuals (δt).This parameter has been already used when examining an experimental system of pedestrian egress [10,11] and it will be analyzed in what follows.

B. Distribution of passage times
The survival functions of δt [P (T > δt)] for the simulation data presented in Fig. 4(a) are shown in Fig. 5. Importantly, we obtain the power-law decays reported in recent experimental works characterized by the high degree of competitiveness and the strong contact forces among the individuals [10,11].The fittings to these power laws have been performed with the Clauset-Shalizi-Newman method [36] which, apart from the exponent values (reported in the caption of Fig. 5) provides an estimate of the fit accuracy, based on a number of realizations.In all cases, we find a p value > 0.1.In this figure we also find that, as in the experiments, increasing the desired velocity (our parameter for competitiveness) the probability of obtaining longer clogs augments (i.e., FIS effect).Finally, let us stress that apart from this qualitative agreement we obtain good quantitative correspondence as the maximum duration of long time lapses δt are of the same order of magnitude as in the experiments [10].

C. Burst size
A readily accessible variable that can be calculated from the individual passage times (δt) is the size of the bursts (S) defined as the number of individuals that egress between two consecutive clogs.To this end, it is necessary to define a clog (δt c ) as the time interval above which the flow is said to be interrupted.In order to compare with experiments, we choose δt c = 0.7 s observing that our model reproduces the well-known exponential distribution [11] (see Fig. 6).This exponential behavior indicates that the probability of occurrence of a clog is constant during the stationary egress process as pointed out for granular systems [37][38][39] and the flow of sheep through bottlenecks [40].As expected, when normalized by its mean burst size, all the exponentials collapse into a single curve, indicating that the only difference between distinct cases is the average value.

D. Influence of rotational parameters
Once we have proved the ability of the numerical model to reproduce experimental results of room evacuation we can study the effect that the proposed rotational parameters have on the flow rate.To this end we will look at the average specific flow rate (Q s ) during the whole simulation run which can be computed as Q s = N tot T end b , where N tot is the total number of pedestrians that cross the door in the T end = 1600 s simulated time.Note that, despite the already mentioned existence of a transient during the first 8 s of the evacuation, this time is so small in comparison with the total simulation time that its effect is negligible.
In Fig. 7 we report the specific flow rate versus the desired velocity for different values of directional strength S D .First, we observe that the FIS effect is robust and appears no matter what the value of S D is.The influence of increasing S D is reducing the flow rate, in principle because the effective diameter of simulated pedestrians increases when they have more tendency to align parallel to the door.Interestingly, the optimum desired velocity at which Q s is maximum is the same for all the values of S D studied: Complementarily, we can look at the variations in the distribution of δt when changing this rotational parameter.Figure 8 shows that the lower values of the exponent α are found for the higher values of S D .In other words, the duration of clogs increases with the directional strength.Therefore, the reduction of the average flow rate observed in Fig. 7 when S D was increased can be attributed to the increase of the clogging times reported in Fig. 8.
The other relevant parameter related with the rotation of simulated pedestrians is the angular noise η.In Fig. 9 we display the outcomes obtained for the average specific flow rate versus the desired speed when implementing different values of η.Several interesting features can be observed.First, increasing the rotational noise produces a fluidization of the system generating greater flow rates.Second, the FIS effect is robust under the presented variations of angular noise.And third, the position of the v d at which Q s is maximum seems to be sensitive to this parameter.This result contrasts with the lower sensitivity of this optimal desired speed with the directional strength reported in Fig. 7. Finally, let us stress once more that the values of the specific flow rate obtained are comparable to the ones reported in experiments [11] (about Q s ∼ 4 s −1 m −1 ), in particular for η = 1.25.
In summary, modifying the rotational parameters does affect the absolute values of specific flow rate but the relation between them for different competitiveness is preserved, at least for the ranges explored.Indeed, a maximum of Q s is always present, which means that the FIS effect is robust under these changes in the parameter's values.

IV. CONCLUSIONS
In this work we propose a force-based model introducing an anisotropic shape for pedestrians.This modification increases the degrees of freedom and allows a more realistic description of the room evacuation process in highly competitive conditions.In particular, several features of experimental competitive egress are naturally reproduced: (a) the distributions of passage times between successive pedestrians (δt) exhibit power-law tails [10]; (b) the distribution of normalized burst size displays an exponential tendency, also in accordance to the experimental results [11]; (c) the faster-is-slower effect is present; and (d) the figures of the specific flow rate obtained are also compatible with the measured ones [11].
In addition, a thorough variation of the rotational parameters evidences that the faster-is-slower effect is robust under variations of rotational strength and noise showing that the flow rate grows either when increasing the rotational noise amplitude, or when decreasing the rotational strength.Indeed, we argue that both strategies lead to a reduction of the times during which the system remains blocked as evidenced by the changes of the passage times' survival functions.
The presented model and its capacity for reproducing the new experimental data of pedestrian egress in competitive conditions suggests the importance of considering torques among the individuals, being a step forward in counting on validated models for this kind of situation.Obviously, there are some aspects of the model, such as the way of introducing rotational noise, for which robustness should be analyzed with more depth.Indeed, we believe that this work opens a new field of possibilities allowing one to design particular layouts able to improve the safety of people in competitive situations.

APPENDIX: INTERPARTICLE GRANULAR FORCE
To calculate the granular force between two pedestrians i and j , F ij , we use an algorithm based on the concept of spherocylinder [15], which can be spatially defined by two vertices and a spheroradius, r.Thus, the surface of a spherocylinder is delimited by all the points at distance r from an edge between the two vertices.The detection of a contact between two spherocylinders is therefore reduced to find the closest points between two edges.If the distance between these two points is smaller than 2r a contact force would exist.This contact force F ij can be decomposed as F ij = F n n + F t t, where F n is the component normal to the contact plane n and F t acts on the tangential direction t as in Fig. 1.To determine the normal component F n , we use a linear-dashpot interaction force, which is governed by the overlap distance δ between the two spherocylinders.It accounts for dissipation, assuming a velocity-dependent viscous damping.Hence, the total normal force reads as F n = −k n δ − γ n v t rel , where k n is the spring constant in the normal direction, γ n is the damping coefficient in the normal direction, and v n rel is the normal relative velocity between i and j .The tangential force F t also contains an elastic term and a tangential frictional term and it reads as F t = −k t ξ − γ t v t rel , where γ t is the damping coefficient in the tangential direction, and v t rel is the tangential component of the relative velocity of the overlapping pair.ξ (t) represents the elongation of an imaginary spring with elastic constant k t at the contact, and increases as dξ (t)  dt = v t rel provided that there is an overlap between the interacting particles.Finally, F t is truncated as necessary to satisfy the Coulomb constraint , where μ is the microscopic friction coefficient.
Using those parameters we estimated t c = π m r 2k n and set an integration time step t = t c 300 .We simulate a square room of 8 m × 8 m; to mimic rigid walls, we considered the wall as fixed particles with the same normal and tangential stiffness as the particles.

FIG. 1 .
FIG.1.Sketch illustrating (a) the contact interaction between two spherocylinders (see Appendix for more details); and (b) the orientation of the acting direction e i of the driving force.The blue dots indicate the end points of the target segment that pedestrians aim to reach.

FIG. 3 .
FIG. 3. Snapshot illustrating a typical run where the individuals aggregate behind the door forming a cluster of n ≈ 180 pedestrians.Periodic boundary conditions are implemented in both directions.

25 FIG. 4 .
FIG. 4. Specific flow rate.(a) Discharge curves for S D = 20 N m, η = 1 and three different desired velocities as indicated in the legend.The inset displays the transient at the very beginning of the simulations.(b) Evolution of the instantaneous specific flow rate calculated for a sliding window of 200 outgoing pedestrians.

FIG. 5 .
FIG. 5. Survival function of the time lapses δt between the passage of two consecutive pedestrians obtained for an angular noise η = 1 and directional strength S D = 20 N m.Different symbols are used for different desired velocities as indicated in the legend.The solid lines represent the fittings obtained with the Clauset-Shalizi-Newman method [36] and the corresponding power-law exponents are α 1 = 4.6 ± 0.4, α 2 = 3.1 ± 0.2, and α 3 = 2.6 ± 0.1.

FIG. 6 .FIG. 7 .
FIG. 6. Distribution of burst sizes normalized by its mean, for δt c = 0.7 s.The data have been obtained for simulations with η = 1, S D = 20 N m and three different desired velocities v d as indicated in the legend.The continuous line is an exponential fit of the data.

10 S D = 15 S D = 20 FIG. 8 .
FIG. 8. Survival function of passage times for different directional strengths (S D ). Results obtained for simulations performed with η = 1 and v d = 0.5 m/s.The values of the exponent α are 6.6, 5.9, and 4.7, for S D = 10, 15, and 20, respectively.

FIG. 9 .
FIG. 9. Average specific flow rate versus pedestrian desired speed for S D = 20 N m and different values of the angular noise amplitude η as indicated in the legend.