Quantitative Comparison Between Crowd Models for Evacuation Planning and Evaluation

Crowd simulation is rapidly becoming a standard tool for evacuation planning and evaluation. However, the many crowd models in the literature are structurally different, and few have been rigorously calibrated against real-world egress data, especially in emergency situations. In this paper we describe a procedure to quantitatively compare different crowd models or between models and real-world data. We simulated three models: (1) the lattice gas model, (2) the social force model, and (3) the RVO2 model, and obtained the distributions of six observables: (1) evacuation time, (2) zoned evacuation time, (3) passage density, (4) total distance traveled, (5) inconvenience, and (6) flow rate. We then used the DISTATIS procedure to compute the compromise matrix of statistical distances between the three models. Projecting the three models onto the first two principal components of the compromise matrix, we find the lattice gas and RVO2 models are similar in terms of the evacuation time, passage density, and flow rates, whereas the social force and RVO2 models are similar in terms of the total distance traveled. Most importantly, we find that the zoned evacuation times of the three models to be very different from each other. Thus we propose to use this variable, if it can be measured, as the key test between different models, and also between models and the real world. Finally, we compared the model flow rates against the flow rate of an emergency evacuation during the May 2008 Sichuan earthquake, and found the social force model agrees best with this real data.


Introduction
Crowd simulation is an area that has been the subject of a significant amount of multidisciplinary work over the last few decades [1][2][3]. Its applications range from simulating crowds for movies [4,5] and games [6,7] to analyzing pedestrian behavior [8][9][10] and preparing for fire evacuations and similar emergencies [11][12][13]. The earliest attempts to simulate crowds generally adopted a macroscopic approach [14,15], where there is no explicit notion of an individual. Later, with increasing computational resources and with availability of observational data on an individual level [16][17][18], modelers were able to develop microscopic approaches [5,8] for application in areas where it was necessary to model and analyze individuals in the crowd. For example, in a simulation of evacuation, knowl-edge of the movement of crowds could reveal methods to improve crowd flow and evacuation speed.
One of the earliest and seminal works in individualbased motion planning was Craig Reynolds's model of coordinated animal motion such as bird flocks and fish schools [5]. Okazaki and Matsushita [19] assigned magnetic poles to goals, agents and obstacles to model movement. Subsequently, Helbing's Social Force Model [8] was developed, which is still one of the most popular models of movement in crowds. More recently, there have been several velocity-based approaches to motion planning, such as the synthetic vision based model [20] and the Reciprocal Velocity Obstacle Model [21][22][23].
The advocation of simulation-based analysis has become increasingly common over the last decade. Some well know applications include analysis of the yearly Mus-arXiv:1401.0366v2 [cs.MA] 10 Jan 2014 lim Hajj [24], or more recently the Love Parade disaster, Germany 2010 [25]. In the case of the latter, models and expertise had been used to guarantee the safety of the event, only for unforeseen circumstances to result in the deaths of 21 individuals. Clearly these real world examples emphasize the critical role that crowd modelling plays in safety preparation and planning, this in turn emphasizes the need for understanding the model dynamics, limitations and similarities. The extent to which these models are capable of accurately predicting the motion of a crowd is therefore critical for planning and safety.
Ideally these models should be validated against realworld data from all scenarios, for all cultures and for all varieties of crowd composition. Unfortunately real world data regarding egress or emergency situations is limited, often incomplete (the initial conditions are hard to know) and certainly not controlled. Even in the rare circumstances where data is available this often describes measurements and phenomena at the macro-scale, e.g., flow, density, average speed, etc. Because these models of individuals exhibit emergent behaviors at the mesoscopic and macroscopic scales, it is in general hard to tell whether the dynamics at the individual level are correct. The best researchers can realistically hope for is to collect microscopic data for a single scenario to calibrate the model. An alternative approach to the data-intensive one is to instead look quantitatively at the fundamental dynamics of the models and understand where these models differ, at the same time identifying common aspects of the models. This will aid in understanding fundamental aspects of all crowd models and offer insight into real-world dynamics of human crowds.
While much of the existing research does offer basic comparison of proposed models with existing models to demonstrate their usefulness, there is no existing quantitative comparison of the differences and similarities of these models to tell which is most accurate. The objective of this paper is to demonstrate a methodology for quantitative comparison of simulation models in a simple egress simulation. For this, we choose three popular models that are structurally very different. Since our focus is on the comparison methodology, we do not worry about specific versions of the models, in particular recent modifications and enhancements that are supposed to produce more realistic crowd behaviors.
The contribution of the paper is then three-fold: firstly, the analysis and comparison of these models provides interesting insight into the consequence of adopting each in particular forms of crowd or egress simulation. Secondly, the systematic approach we describe could be used in future to compare further models and develop a standard method of comparison for crowd simulation. Finally, the paper conclusion identifies a single measurable metric that is most effective in distinguishing the behaviour of the models.
The remainder of this paper is organised as follows. The models that are consider in the comparison are first described in Section 2; Section 3 describes the experiments that were conducted and the methodology for anal-ysis. Significant observations from the simulations are presented and analyzed in Section 4. Finally, Section 5 concludes the paper.

Models
In this section we describe the three individual (or microscopic) models that are compared in this paper. We use the general term agent to refer to individuals within each of the models. The lattice gas model is a probabilistic approach where the future location of an agent is probabilistically determined based on the current configuration of its neighborhood. This implies that the same initial configuration can produce different results during different runs. On the other hand, in the social force model and reciprocal velocity obstacle model, agents find their ways to their destinations via deterministic collision-avoiding calculations. As a result these models produce the same result on multiple runs for a particular initial configuration.

Lattice gas model
A Cellular Automation (CA) model is one in which space and time are discrete. Furthermore, the state space of a CA model is also discrete and finite. In each time step the values of all cells are updated synchronously based on the values of cells in their neighborhoods. Depending on the type of neighborhood (i.e., von Neumann, Moore), and the type of lattice (triangular, square, hexagonal, etc.), the exact number of cells in the neighborhood of a given cell can vary [26]. Lattice gas models are CA models that make use of a discretized version of the Boltzmann transport equation to model motion [27][28][29]. Advances in this modeling area include the Extended Floor Field Model [30] in which agents interact through virtual traces that act like the pheromones in chemotaxis, and the SWARM information model [31] which uses multiple floor fields to model transmission of knowledge between agents. In this study we use the model by Tajima and Nagatani [32]. In contrast to later CA models with additional features like bi-directional movement [33] and vision impairment during evacuation [29], this original model is simple, and adequate for our model comparison purpose.
A square grid is used, and each cell can be occupied by at most one person. The person performs a random walk biased towards the single exit in the room. For example, in Fig. 1(a), the probability that the person takes a step in the +y direction is P y = (1 − D)/3 + De y /(e x + e y ), while the probability that it takes a step in the −x direction is P −x = (1 − D)/3 + De x /(e x + e y ), since the intended direction e has positive projections along the +y and −x directions. Along the +x direction, for which e's projection is negative, the probability that the person takes a step in the +x direction is P x = (1 − D)/3. Here we see that the random walk probabilities contain an unbiased component (1 − D)/n, which is the same for all n permissible directions (in the Tajima-Nagatani model, the −y Intended direction

Intended direction
Intended direction direction is not permitted), as well as a biased component D = 0.7 that favours movement towards the exit.

Social force model
The social force model is one of the most popular models for motion planning in crowds [34][35][36]. This model is based on the idea that pedestrians move in response to fictitious attractive or repulsive social forces produced by obstacles and other pedestrians. Over the years, several extensions have been made to the social forces model like the modeling of grouping behavior [34]. In this paper, we use the Helbing-Molnár-Farkas-Vicsek (HMFV) social force model [8], which is tweaked from the original model introduced by Helbing and Molnar in 1995. In the original model proposed by Helbing and Molnar [8], three types of forces act on a given agent i with mass m i , instantaneous position r i (t), and instantaneous velocity v i (t). The first is a restoring force that steers the agent towards the desired velocity v 0 at a rate determined by the characteristic time τ i = 1. Here where V 0 is the preferred speed, e i is a vector that points towards the exit, and (1 − p) is the weight given to this desired velocity. With a weight of p, agent i also adapts to the average velocity v j i in its neighborhood. When p is small, agent i moves more along its intended direction e i (t), whereas if p is large, agent i tends to follow where its neighbors are going. We can therefore tune p from p ≈ 0 (self-directed normal egress) to p ≈ 1 (panic-driven herding during emergency evacuations). In this paper, we used p = 0.2 to simulate an emergency evacuation situation. The second is a repulsive force where that mimics the psychological tendency of agents i and j to move away from each other if they are too close.
Here R ij = R i + R j is the sum of radii of the two agents, d ij = |r i − r j | is the physical distance between the two agents, and n ij = (n 1 ij , n 2 ij ) = (r i − r j )/d ij is the unit vector pointing from agent j to agent i. Additionally, A = 2000 N and B = 0.08 m are the repulsion coefficient and the fall-off length of interacting agents respectively [37]. Helbing and Molnar also found it necessary to introduce two other terms in the interaction force when agents i and j are in contact with each other, i.e. d ij < R ij . This counteracting body compression term kη(R ij −d ij )n ij and sliding friction term κη(r ij − d ij )∆v t ji t ij are crucial for getting realistic behaviors of panicking crowds. Here t ij = (−n 2 ij , n 1 ij ) is the tangential direction and ∆v t ji = (v j − v i ) · t ij is the tangential velocity difference.
Finally, in the third force the first and second terms repels agent i from a wall that it is d iW away from, while the third term (which is negative) is introduced to mimic the observation that people move faster near walls when they are in crowded situations. Here, n iW is the normal vector of the wall, and t iW is the tangent vector of the wall. Also, k = 12, 000 kg/s 2 and κ = 24, 000 kg/ms are respectively the body force constant and the sliding friction force constant used.

Reciprocal velocity obstacles model
The third model we consider is the reciprocal velocity obstacles (RVO) model. In this model, the time to the next collision is calculated based on the relative velocities of N agents. Each agent then changes its velocity {v i } N i=1 to maximize this time to collision. By continuously updating the velocities in this manner, collisions are avoided. This algorithm was first proposed by Fiorini et al. [38], but was first used in multi agent systems in [39]. Since then there have been several modifications and improvements to the algorithm [10,22,40,41], although the underlying idea remained the same. In this paper, we use the RVO2 model introduced by Guy et al. [22], where the collision avoidance computation is based on computational geometry and linear programming.
Given a preferred velocity, RVO2 helps an agent find the velocity closest to the preferred velocity that will enable it to avoid collisions with all other agents. This is done by determining for each neighboring dynamic and static obstacle an Optimal Reciprocal Collision Avoidance line (ORCA line) . Each ORCA line determines the half plane in the velocity plane where the agent's velocity should lie to ensure that no collision occurs with that particular obstacle for the next τ seconds. τ is a parameter that specifies the number of seconds for which the chosen velocity should avoid a collision. From the set of half planes thus obtained, the optimal collision-avoiding velocity can be determined as the velocity within the permissible region for all the half planes that is closest to the preferred velocity. This can be determined efficiently using linear programming. Besides the speed and efficiency of the algorithm, another strong appeal of this model is that we need only set one parameter τ . For the experiments in this paper we use two values of τ ; 0.5 seconds or 10 simulation time steps for avoiding dynamic agents and 0.05 seconds or 1 simulation time step for avoiding static obstacles. The value of τ translates in practical terms to how early a walker tries to avoid a collision with an obstacle.

Methodology
In this section we explain the simulations that were carried out and the analysis done. Java 6 and the MASON simulation framework [42] were used for creating and running the simulations. The simulated scenario consists of agents evacuating from a simple rectangular room with a single exit. The dimensions of the room are shown in Fig. 2. Evacuation of the room was simulated with N = 50, 100, 150, 200, 300, 500 and 1000 agents whose locations were randomly distributed within the room. For each N , we ran 100 simulations for each model, and the positions and velocities of all agents at all time steps were collected. For meaningful comparison between the three models, we assumed that each agent is a perfect circle with radius R = 15 cm [43], that has a preferred speed of 1.3 ms −1 , and a maximum speed of 2.6 ms −1 [44]. For the social force model we also assumed that all agents have the same mass of 60 kg.
Where possible, we have used the same settings for all three models to ensure a fair comparison. The values used for the remaining parameters required for the models was given in Sec. 2 along with the respective models. However, to the best of our knowledge, only the social force model has been calibrated against real world data [45]. The other two models are not calibrated beyond choosing a reasonable and believable model and the values chosen are taken such that they produce empirically believable results.

Measurements
Depending on the context and motivation of the model there are several metrics that are generally used. For CA models that are generally used for studying macroscopic patterns, it is common to measure macro variables like the mean evacuation time [29] and flow rate [32]. For models like social force where the motivation is to study interactions at a more granular level, it is common to make empirical and qualitative observations like the lane effect [8]. Sometimes quantitative metrics like density-dependence of the flow or velocity are also measured [46]. From a graphics perspective, where performance is key, computation time [20], frame rate and run time per frame are generally measured against the number of walkers [23].
In this paper, we performed two main classes of measurements: time-based and distance-based measurements. For time-based measurements, we first measured the evacuation time distributions for different number of agents. To better understand the stages in the evacuation, we also divide the room into six different zones (see Fig. 3) to measure the evacuation time sub-distribution. The outer radii of Zones 1, 2, 3, 4, 5 and 6 are 5 m, 10 m, 15 m, 20 m, 25 m and 30 m respectively. We also measure the flow rate at the exit as a function of time for each model.
For distance-based measurements, we traced the trajectories of all agents to obtain the distribution of total distance travelled going from the initial position to the exit. Here, the total distance travelled D i by agent i is just the sum of its displacements Ti t=0 ||r i (t + 1) − r i (t)||, where T i is the evacuated time of agent i. We also divide the room into 100 × 100 cells, and count the number of agents passing through each cell as the simulation progresses. This is then visualized as a heat map. Finally, we measure the ratio of the total distance travelled D i by agent i to the minimum distance D min it would cover if it evacuated along a straight line (For lattice model, D min is the Manhattan distance between initial position of agent and exit). This in some sense quantifies the 'inconvenience', I i experienced by the agent i during the evacuation.

DISTATIS
After making six different measurements on three models, we discovered that when we look at the flow rates, the lattice gas model is more similar to the RVO2 model. However, if we look at the total distance travelled, the social force model is more similar to the RVO2 model. If we look at the zoned evacuation times, we find that the three models have very different distributions (see Sec. 4.1). Since we do not know a priori which measurements best discriminate the three models (or for that matter, best discriminate models from the real world), we want to be able to compare the three models quantitatively, incorporating information from all six measurements.
To do so, we compute the Jensen-Shannon divergence [47] D between the distributions f kµ (z) and f kν (z) for measurement k of models µ and ν. Here, z is a continuous variable like evacuation time, distance travelled, or inconvenience, is the Shannon information function, and If f kµ (z) and f kν (z) are highly similar, we will get D JS [f kµ , f kν ] ≈ 0, whereas if f kµ (z) and f kν (z) are very different, D JS [f kµ , f kν ] 0, i.e. the Jensen-Shannon divergence qualifies as a distance metric. In this way, we obtain six 3 × 3 distance matrices D k .
Then, we use the DISTATIS method [48] for analyzing multiple distance matrices. This is a generalization of the method of principal component analysis (PCA). STATIS is an acronym for the French expression 'Structuration des Tableauxà Trois Indices de la Statistique', which roughly means 'structuring three way statistical tables'. The difference between a straightforward application of the PCA method and the DISTATIS method is shown in Fig. 4. Instead of six different PCAs for the distance matrices obtained from the six different variables measured, and inevitably end up with conflicting conclusions on which models are more similar, we standardize the six distance matrices, and ask which variables are more similar to each other.
The idea behind DISTATIS is that the data points are similarly clustered in two variables, if their standardized distance matrices in these two variables are similar to each other. Therefore, if we compute the cross correlations between variables, we will discover which variables give more similar outcomes to which other variables through a PCA. Components of the first eigenvector obtained from this PCA tells us how important each variable is. By weighting the distance matrices of each variable by its component in the first eigenvector, we construct a compromise matrix whose matrix elements give us the most reliable similarity/dissimilarity between data points. A final PCA of this compromise matrix then gives the most reliable groups of data points. In Fig. 12 Fig. 4. In many classification problems, we can define different distance metrics for the same set of data points to measure how dissimilar they are to each other. In this figure we show three different distance matrices D1, D2, and D3. If we perform principal component analysis (PCA) on them separately (PCA1, PCA2, PCA3), we are not likely to agree on which data points are more similar to each other. In the DISTATIS method, we first standardize the distance matrices by centering and rescaling, to get S1, S2, and S3. If two data points are indeed more similar to each other than they are with the third, then we expect this similarity structure to be reflected in S1, S2, and S3, i.e. two standardized distance matrices would have very similar patterns of matrix elements. With this in mind, we reshape the matrices S1, S2, and S3 to make them into vectors of distances. We then compute the cross correlations between these vectors, before performing PCA on the correlation matrix C. This PCA tells us which distance matrices agree with each other more, and which less. The result is a consensus between different distance metrics, based on which a single PCA will discover the similarity structure between the data points.

Evacuation time
As Fig. 5 shows, the evacuation time distributions for all models have a uniform sub-distribution in the middle of the evacuation, which signals congestion at the exit, once we have N > 50 agents. However, in the RVO2 model this uniform sub-distribution is flanked by two small peaks. These two small peaks suggest a higher evacuation rate just before congestion sets in and just before congestion dissipates in the RVO2 model. This is not observed in real crowds, and hence is clearly an artifact of the model.
Comparing the three models, we find that the social force model has the longest evacuation times of up to five minutes. More importantly, careful inspection reveals fluid-like crowd dynamics around the exit in the RVO2 and lattice gas models. These can be seen in our videos http://www.youtube.com/watch?v= P64p3nlH_P4 and http://www.youtube.com/watch?v= vJ0Nzi5Bykw. In the RVO2 model, an obvious back flow can also be seen during the sudden gush of agents toward the exit. Due to the strong repulsion forces, we see a solid-like behavior in the social force model when the exit becomes congested (see our video http://www.youtube. com/watch?v=3Ot7m959_yo).
In addition, we divided the room into six zones, as shown in Fig. 3, based on the distance to the exit, to determine how strongly the crowds mix in the three models. The distributions of zoned evacuation times are shown in Fig. 6. These tell us that mixing is strong in the lattice gas and RVO2 models, but weak in the social force model. This mixing dynamics can be seen more clearly from our videos http://www. youtube.com/watch?v=qeoJotgEUxk (lattice gas), http: //www.youtube.com/watch?v=uZpd5LODYZs (RVO2), and http://www.youtube.com/watch?v=wEB6Ya0o_yw (social force). From the movies of the lattice gas and RVO2 simulations, we see that agents from the nearer zones get pushed to the side once there is congestion at the exit. In these two models, agents prefer to keep moving when it is possible for them to do so. Agents in the social force model behave differently: once the exit becomes congested, they will become nearly stationary and wait for their turn to go through the exit.

Passage density
The passage density heat map of agent locations averaged over 100 simulations provides information about the spatial-temporal trace of agents during their evacuation. This is crucial when analyzing the spatial structure of congestions. In Fig. 7(a), we see a channel leading straight through the exit. This channel acts as an attractor, because once an agent gets onto this, it will be forced statistically towards the exit (P t,y = 1.00).
Though not as pronounced as for the lattice gas model, two left-right symmetric channels also formed in the RVO2 model (Fig. 7(c)). These two channels point roughly from the centers of mass of the left and right halves of the room towards the exit. No such artificial channels can be seen in the passage density map of the social force model (Fig. 7(b)).
Comparing the three models, we also see that the RVO2 model (Fig. 7(c)) has the most compact passage density, whereas the social force model (Fig. 7(b)) has the least compact passage density. In particular, in the region just in front of the exit, the RVO2 model gives rise to very high passage density. This high passage density cannot be attributed to agents stepping over a cell once or twice, as in the lattice gas and social force models. In the RVO2 model, this high passage density is produced by repeated visits to the same cell by individual agents that are always moving.

Total distance traveled distribution and inconvenience
To measure the difficulty (inconvenience) for an agent to squeeze through the crowd during evacuation, we measured the total distance traveled and as well as the inconvenience (defined in Eq. 7). As we can see in Fig. 8(b) and Fig. 8(c), the distributions of total distance traveled are very similar for the social force and RVO2 models. The lattice gas model, on the other hand, has very different distributions of total distance traveled ( Fig. 8(a)). This is because space is continuous in the two former models, but discrete in the lattice gas model.
However, when we compare the distributions of inconveniences in Fig. 9, we find the lattice gas model is more qualitatively similar to the RVO2 model, though the typical inconveniences in the lattice gas model are larger than the those in the RVO2 model. For different N , the distribution peaks around the same inconvenience for both models. In contrast, the peak position of the inconvenience distribution of the social force model depends strongly on the number of agents in the simulation.

Evacuation rate
We also measured the flow rate of agents at the exit. As shown in Fig. 10(a), the flow rate of the lattice gas model saturates at around 11.8 persons per second due to the finite sizes of the agents and the exit. The social force model has the lowest flow rate of 2.2 people per second ( Fig. 10(b)), because agents slow down dramatically at the end as a result of strong repulsive forces between agents. From Fig. 10(c)), we see the same peaks before and after the congestion period for the RVO2 model, a consequence of the liquid-like dynamics of RVO2 crowds.
In an attempt to compare simulation results with real data we obtained a number of videos showing reallife evacuations [49,50] during the Sichuan Earthquake on May 12, 2008. We analysed two such videos, http: //www.youtube.com/watch?v=-y38ebiAnQw and http: //www.youtube.com/watch?v=e1yapP3z_L4. Because the geometry of the scene and the viewing angles are not known, we could only measure the flow rates on a secondby-second basis from these two videos (Fig. 11), for comparison against what we measured in our simulations. Unlike the computational study, where we could average over 100 simulations to get smoothly varying flow rates, the flow rates measured from the videos are noisy. However, sensible comparisons can still be made. In http://www. youtube.com/watch?v=-y38ebiAnQw, the crowd is shown to evacuate through the check-out counters of a Walmart supermarket. The flow rate was high, but there was no congestion, because of the wide exit. This flow rate cannot be compared against our simulations of evacuation through a narrow exit.
In http://www.youtube.com/watch?v=e1yapP3z_L4, students were evacuating their classroom through a narrow door. From the video, we see that the width of the door allows no more than three school children to exit simultaneously. Therefore, the exit dimensions seem to correspond to the exit dimensions use in our simulations. From Fig. 11, we see that the average flow rate in the congestion phase is about 3 persons per second. This is far below the congestion flow rate in the lattice gas model, and very close to the congestion flow rate of the social force model. Based on this observable, we find that the social force model agrees best with the real data. This small test also tells us that the three models can be distinguished through their congestion flow rates, if other observables could not be measured.

DISTATIS comparison
Finally, we make use of the DISTATIS method explained in Sec. 3.2 to compare the three models with all six different measurements. Fig. 12 shows the different measurements of the three models projected onto their first two principal components. The barycenters for the three models are also shown. From Fig. 12 we see that based on the total distance traveled (green), the social force and RVO2 models have similar distributions, which are both very different from that of the lattice gas model. On the other hand, the lattice gas and RVO2 models have similar distributions that are significantly different that of the social force model.
More importantly, we find the clustering of some observables: evacuation time, passage density, and flow rate all tells us that the lattice gas model is similar to the RVO2 model, and different from the social force model. We also find the inconvenience and zoned evacuation time provide some indication that the three models are qualitatively different. In fact, the DISTATIS method indicates that    the zoned evacuation time is one key metric that clearly distinguishes the three models. Therefore, any real world experiments should attempt to measure this. This observable maximally discriminates between the three models and therefore it is reasonable to expect it will also discriminate between real-world data and the models.

Conclusion
To conclude, we simulated emergency evacuation from a rectangular room with a single narrow exit using three models: (1) lattice gas, (2) social force, and (3) RVO2, with different starting number of agents in the room. From these simulations, we measured six observables: (1) evacuation time, (2) zoned evacuation time, (3) passage density, (4) total distance traveled, (5) inconvenience and (6) flow rate. Besides qualitative comparison of the various distributions obtained, we also compared the three models quantitatively using the DISTATIS method. Comparing the evacuation time, passage density, and flow rate distributions, we find that the lattice gas and RVO2 models are similar, both of which are very different from the social force model. On the other hand, if we compare the distribution of total distance traveled, the social force and RVO2 models are more similar to each other, and very different from the lattice gas model. We compared the simulated congestion flow rates of the three models to congestion flow rates of school children evacuating from their classroom during the May 2008 Sichuan Earthquake, and find based on this observable that the social force model agrees best with the egress data. Finally, from an analysis of the DISTATIS we have identified the zoned evacuation time as the one observable metric that can best discriminate between these models, and also between models and real-world data