Continuous Limit of a Crowd Motion and Herding Model: Analysis and Numerical Simulations

In this paper we study the continuum limit of a cellular automaton model used for simulating human crowds with herding behaviour. We derive a system of non-linear partial differential equations resembling the Keller-Segel model for chemotaxis, however with a non-monotone interaction. The latter has interesting consequences on the behaviour of the model solutions, which we highlight in its analysis. In particular we study the possibility of stationary states and the formation of clusters. We also introduce an efficient numerical simulation approach based on appropriate hybrid discontinuous Galerkin, which in particular allow flexible treatment of complicated geometries. Extensive numerical studies also provide a better understanding the strengths and shortcomings of the herding model, in particular we examine trapping effects of crowds behind non-convex obstacles.


Introduction
In recent years, the modelling of the behaviour of (large) human crowds became more and more important.This is mainly due to the fact that an increasing number of people live in cities or other areas with a high number of inhabitants.Of special interest in this context are emergency situations such as evacuations or jamming effects in large events, where herding effects seem to play an important role.The majority of currently used models has a limited capability to simulate the latter, which motivates further research on herding models as the one we shall consider in this paper.In general, models for the behaviour of human crowds are divided into two categories: microscopic and macroscopic models.In microscopic approaches, every individual in the crowd is treated separately (as a "particle").Among the first of these were the so-called social force models, cf.Helbing et al. (2002).Further examples for these kind of models are cellular automata, cf.Fukui & Ishibashi (1999); Muramatsu & Nagatani (2000).In contrast, macroscopic models deal with the whole crowd, often represented by a density function depending on ('continuous') space and time.One example for this class is the model by Hughes, cf. Hughes (2002), which is motivated by considering the crowd as a "thinking fluid".While that one is a classical PDE model (cf. Di Francesco et al. (2010)) there are also models using other mathematical techniques such as optimal transport or mean field games, cf.Maury et al. (2009); Dogbé (2010); Lasry & Lions (2007).For an extensive review on crowd modelling, see cf.Helbing (2001).We finally note that crowd motion models share many features with traffic models, cf.Aw & Rascle (2000).This paper deals with the analysis and numerical simulation of a macroscopic model for the motion of a human crowd, derived by (formally) passing to a continuous limit from a microscopic cellular automata model developed by Kirchner and Schadschneider, cf. Kirchner & Schadschneider (2002).In their approach, the crowd is considered as a group of a finite number of individuals living on a rectangular two-dimensional grid.Giving a discrete time step, the model provides for each individual in a given cell the probability to jump into a neighbouring cell.This probability is determined by several factors.First of all, individuals are not allowed to jump to an occupied cell (size exclusion, cf.Simpson et al. (2009)).Furthermore, there exist two driving forces, respectively, a static field S and a dynamic field D on which the jump-probability depends exponentially.The static field provides the individuals with a sense of their environment, increasing towards locations they want to reach, such as doors.The dynamic field D, being zero at the initial time, models the tendency of people to follow others (commonly referred to as herding, cf.Helbing (2001)).On a given cell, the value of D increases whenever a individual leaves the cell.Furthermore, it is assumed that D relaxes and diffuses.With these ingredients at hand, the probability to jump to a cell i, j is given by Here n i,j is the so-called occupation number being one if a cell is occupied and zero otherwise.Using the positive constants k D and k S , the relative influence of the two fields can be regulated and N is a normalisation factor.Formally passing to the limit of cell length ∆x and discrete time step size ∆t to zero under the scaling ∆t = (∆x) 2 we arrive at the following system of partial differential equations This model can be seen as a non-linear variant of the famous Patlak-Keller-Segel model for chemotatic movement of cells (cf.Keller & Segel (1970);Horstmann (2003Horstmann ( 2004))).From the modelling point of view, the main difference is the fact that in these models the dynamic field corresponds to a so-called chemoattractant.This is a chemical emitted by the cells if e.g. they found a food source and which has the function to lead other cells to this source.In our model, in turn, D only represents information which can be understood as a "virtual chemoattractant", cf.Kirchner & Schadschneider (2002).Furthermore, due to the finite size exclusion in the microscopic approach, this model features a non-linear mobility, proportional to 1−ρ.This prevents the model from showing a blow-up phenomena as known from the original Keller-Segel model for initial data with a mass above a certain threshold (cf.Blanchet et al. (2006)).In this paper we will consider analytical issues such as linear stability of stationary solutions and provide extensive numerical simulations in one and two space dimensions.We shall also provide several results yielding improved understanding of the model by Kirchner and Schadschneider.
• A continuum limit exhibiting connections to (non-linear) chemotaxis and animal herding models (Section 2).
• A discussion of appropriate boundary conditions for modelling crowds in contained environments and a discussion of stationary states, related to the possibility of congestion.
• A discussion of numerical schemes allowing efficient simulations in setups with complicated geometries (Section 4).
• A discussion of particular limitations of the herding model, which occurs in a simulation with non-convex obstacles.Here the local definition of the dynamic fields yields counterintuitive behaviour since the herding only affects local movement of the crowd but not the trend to follow those finding escape routes (Section 5).
2. Derivation from the Microscopic Model Figure 1.Two particles on a one-dimensional cell grid We present the formal derivation of the macroscopic from the microscopic model in one space dimension as the procedure is analogous in higher dimensions.We consider M cells of width ∆x and denote the i-th cell by x i , i = 1, . . ., M .Furthermore, we assume discrete time steps of size ∆x and write t k = k∆t.From Kirchner & Schadschneider (2002) we use the probability for a particle located in either x i−1 or x i+1 to jump into the cell x i (cf. with . Here we denoted by ρ(x, t) the probability to find a particle at position x and at time t.The given function is ξ(x i ) either one or zero at position x i .It is used to model walls and other obstacles.
Using the same notation for D, we can write i.e. the probability to find a particle at x, t + ∆t is given the probability that a particle jumps into this cell minus the probability that a particle leaves the cell if it was already occupied.The value of D increases, whenever a particle leaves an occupied field.To simplify the notation, we introduce h := ∆x, and from now on write n i (t) := n(x i , t), etc.. Taylor expanding (2.1) and taking only the first order terms into account we arrive at Dividing this expression by h 2 , choosing the scaling ∆t = h 2 and passing to the limit h → 0 we obtain the following limiting equation For (2.2), we apply the same procedure.
In the limit ∆t = h 2 → 0, the last term on the r.h.s vanishes and we obtain As proposed in Kirchner & Schadschneider (2002), we also add diffusion and relaxation of D to this equation and thus we finally arrive at with some fixed constants δ, κ ≥ 0.
Remark 2.1.In (Kirchner 2002, Sec. 3.5.2),an alternative definition of the dynamic field is given, namely that the value of the D in a certain cell is increased whenever the cell is occupied, i.e.
We remark that in the continuous limit, this results in a system with a linear coupling in the D equation which reads ) (2.7) This system has already been analysed extensively in the context of chemotaxis, cf.Di Francesco & Rosado (2008).

Basic Properties of the Model
In this section we shall discuss properties of the system derived above.First we note that with an appropriate scaling of t and D, we obtain diffusion coefficients equal to one in both equations and can also remove the factor 2 3 in front of the non-linear coupling in the equation for D. However, for the subsequent analysis, we introduce a new parameter κ ∈ R + and arrive at where n denotes the outward normal on ∂Ω.For realistic boundary conditions on ρ, one would divide the boundary into two parts, namely doors and walls On the wall segments, we prescribe no-flux boundary conditions where the flux is defined by Note that due to the Neumann conditions on D, if ∇S • n = 0, this is equivalent to homogeneous Neumann conditions on ρ (assuming no vacuum and no saturation on ∂Ω w ).At the door we want people to leave the domain with a speed proportional to their density.Therefore, we introduce an outflow velocity v 0 of Euclidean norm 1 pointing outside and prescribe: A typical choice for v 0 would be k S ∇S.Furthermore, we supplement the system with the initial conditions Under appropriate conditions on u I , existence and uniqueness of a weak solution for this system is well-known.For the sake of simplicity we will only state the result with the simplification of homogeneous Dirichlet boundary conditions.We define Then we have the following theorem Theorem 3.1 ((Wrzosek 2004, Thm.2.1)).Assume p > n, p ≥ 2. If furthermore (u I , D I ) ∈ V then for arbitrary t > 0 there exist unique global solutions (u, D) of (3.1), (3.2) with (u(t), D(t)) ∈ V.
We remark that for the realistic boundary conditions stated above, we expect less regularity of the solutions.In particular, the flux may become discontinuous at the edges of doors.The regularity of the flux will depend on the geometry and regularity of the boundary, cf.Markowich (1986); Grisvard (1985).

(c) Stationary Solutions
In this section we examine the stationary system, i.e.
Here we assume no flux boundary conditions on ρ and D. For the first equation, it suffices to find ρ ∞ such that the total flux is zero, i.e.
Integrating this equation we obtain and thus (with (3.8) The constant K is implicitly determined by the total mass To obtain the stationary solutions D ∞ to the second equation, we have to solve the following non-linear elliptic mean-field problem where the right hand side is bounded between zero and minus one.We point out that uniqueness of a solution can be obtained by Banach's fixed point theorem if either δ is large or k D is small enough.
Remark 3.2 (Constant Stationary Solutions).We consider the special case of constant stationary solutions.We note that ρ ∞ = const, D ∞ = const solves (3.6), if S = 0. From (3.9) we find Assuming δ > 0, equation (3.7) reduces to the following simple algebraic expression Thus as long as δ > 0, there always exists constant stationary solutions (ρ ∞ , D ∞ ), uniquely determined by the total mass M .

(i) Linear Stability
Constant stationary solutions correspond to situations without congestion and it is thus important to understand if they prevail under perturbation.We shall now derive stability properties of the stationary solutions considered above.For simplicity, we set k S = k D = 1 and S = 0 in the following and examine the perturbation of the constant equilibria ((u, v) ∈ V) Using this ansatz in (3.1), (3.2), we arrive at and Dividing by and letting → 0 we obtain the following linearisation (3.12) We first analyse this system in one space dimension using the Fourier series method: The ansatz leads to (for all n ∈ Z) The corresponding eigenvalues of this 2 × 2 ODE system are Consider the larger one only, then the condition for λ 2 to be negative is With Dirichlet b.c. for the perturbations, the first mode is the relevant one, and the condition becomes We point out that, at least for simple geometries of Ω, this procedure also works in two space dimension.For example for the case of a quadratical domain Ω, the Fourier series is given by and instead of (3.13) we obtain which however for the first mode still leads to (3.14).In Sec, a, we shall present numerical simulations confirming this result in both one and two space dimensions.
(ii) Plateau Solutions In the case of a constant stationary state being linearly unstable it is natural to look for other stationary solutions or at least meta-stable ones in the dynamics.In particular a large coupling constant k D and a small diffusion coefficient κ are promising parameters to vary, as suggested by the linear stability condition.Rescaling the coupling coefficient to 1 would correspond to a small diffusion coefficient in the equation for ρ, and the case of small diffusion coefficients is indeed the one investigated in the literature:  2010)).This is corresponding to the well-studied limit of small Debye length in the semiconductor drift-diffusion equations (cf.Markowich et al. (1990)).
Here we are facing a mixed situation compared to the two types of models above, which are either based on attractive of repulsive interactions.In our case interactions are attractive for small densities, but repulsive for large densities, which can be understood easily in a quasi-stationary model with zero diffusion in D, i.e., which, inserted in (3.1), yields the forward-backward diffusion equation Note that the equation is always forward parabolic for ρ ≥ 1 2 and around ρ = 0, however for γ := k D δ large there exist intermediate densities such that the diffusion coefficient becomes negative.
From the modelling viewpoint we expect congestion effects to render the constant states linearly unstable.A congestion should be related to a plateau of higher density.We will thus Article submitted to Royal Society investigate asymptotics of plateau-like stationary solutions for small parameter = or (3.17) with c K = 1/K.In order to gain understanding of the asymptotics it is convenient to assume c is given (we then obtain solutions for different masses by varying c), which we shall do in the following.Setting directly = 0 we obtain D ∞ = ρ ∞ (1 − ρ ∞ ) and thus the fixed-point equation (3.18) This equation has at least one fixed point since 0 For many values of c K and γ, this equation has a unique fixed point, thus we only expect the regular asymptotics of constant stationary solutions.The situation is more interesting if there are multiple fixed points, namely two stable and one unstable one, which is indeed obtained for c K = 91 and γ = 20, see Fig. 2 for an illustration.The value of c K is obtained using a piecewise constant approximation of the numerical solution (3.17) and work with the resulting equation for D We define Γ := ∂ Ω \ ∂Ω and denote by d the signed distance to Γ(t), where we choose d to be negative outside and positive inside the plateau.In order to resolve the local change in normal direction we define a local variable z := d(x) and expand the solutions in the form The leading order equation in the interfacial layer is given by This determines the exact form of the transition between the boundary values on the plateaus.
The next order is given by with homogeneous boundary values.In one spatial dimension (due to the absence of curvature effects, ∆d = 0) we have which has a trivial solution g 1 = 0.In higher spatial dimension the mean curvature (equal to ∆d) in the first-order determines the shape of the interface.In order to extract the appropriate equation, we use the non-trivial solution ψ 1 (z) of the homogeneous problem where 0 < λ < 1 is the smallest eigenvalue of the operator.Multiplying the first-order equation by ψ and integrating with respect to z we obtain Article submitted to Royal Society with the coefficients The simplest solution we expect only depends on z and not explicitly on x.Thus, we have This means in 2 spatial dimensions that the interface is a circle or a part of a circle if cut by ∂Ω.
As we shall see such solutions are indeed the stationary states we see in numerical solutions.

Numerical Simulations
In this section we shall discuss numerical methods used for computational experiments.Note that we will not give any details on the calculation of the static field S. In the original definition of the model, attractive regions such as doors or escape routes can be modelled as regions where the value of S is large compared to other parts of the domain, cf.Kirchner (2002).In the cases we are concerned with here we only consider the case where S is given by a constant minus the distance to the door, i.e.
The constant S max is typically defined as the maximum distance to the door.This ensures that S assumes is maximum value at and decreases with increasing distance from the door.As the geometry in our examples is very simple, the distance to the door can be obtained by basic geometric considerations.However, for more complex geometries, it might be necessary to use an Euclidean shortest path algorithm such as Dijkstra's algorithm, cf.Dijkstra (1959).Finally we remark that to visualise the results of the two dimensional finite element simulations, we used the freely available tool Visit (cf.Visit (Online)).

Simulations in one space dimensions
In one space dimension, we use a simple semi-implicit finite difference scheme to solve a linearised version of system (3.1),(3.2), namely with x ∈ Ω 1 := [0, 1] and supplemented with an initial datum ρ 0 and Neumann B.C. for both ρ and D. We divide the domain into N equidistant intervals of length ∆h and denote by ρ i (t) := ρ(i∆x, t), D i (t) := D(i∆x, t), etc. the values of the solution at each grid point at time t.We also discretise time in portions of size ∆t and write t k := k∆t.Then our semi-implicit scheme reads where ∂ i x denotes the discrete first derivative at the ith grid point obtained by a central difference quotient, i.e.
To implement the scheme, we used matrix and vector classes from the NgSolve package, cf.Schöberl (1997) and the sparse direct solver Pardiso, cf.Schenk & Gärtner (20042006).

Simulations in two space dimensions
In two space dimensions, we use a hybrid discontinuous Galerkin method as described in Egger (2009).This method has initially been developed for equations of the form ρ t = div(∇ρ + ρv), with some given velocity vector field v.The method already includes upwind stabilisation for the convection turn and we shall apply it to the linearised system Here, ρ and D are asummed to be given functions.In our implementation, at time t n we take ρ = ρ(t n−1 ) and D = D(t n−1 ).Then, we solve the complete linearised system to obtain ρ(t n+1 ) and D(t n+1 ) and repeat this procedure until we reach t = t final .
Remark 4.1.To guarantee the continuity of the normal component of the numerical flux over interior element edges we project j onto the space H 1 (Ω) in every time step.This is necessary to strictly enforce the mass conserving property of the scheme.A similar strategy and more details can be found in Bastian & Rivière (2003).

Monte-Carlo Simulation
In order to check the consistency between our PDE and the original model, we implemented a Monte Carlo scheme as described in Kirchner & Schadschneider (2002).We used a Mersenne twister, cf.Matsumoto & Nishimura (1998) to create the pseudo random numbers needed.The main issue here is to deal with so-called "conflicts", i.e. the case when two particles want to jump into the same cell.In our implementation, we followed the strategy described in Kirchner (2002).The basic idea is the following: A new parameter λ ∈ [0, 1] in introduced.If two or more particles want to jump to the same cell, this new parameter determines their behaviour: With probability λ, none of the particles jumps and the cell remains empty.With probability (1 − λ), one particle is chosen randomly and jumps into the target cell.

(a) Linear Stability
In this section we will use the numerical schemes described above to to verify the impact of our linear stability analysis in one and two space dimensions.In the one dimensional case, we choose the constant stationary states ρ ∞ = 0.25 and D ∞ = 0.75.We chose the mesh size h = 5e − 4 and time steps ∆t = 0.1.Furthermore, we choose the parameter = 1, κ = 0.001 and δ = 0.25.From the linear stability analysis (3.14), we expect the system to become unstable for k D ≥ 2.66.To verify this behaviour numerically, we add at time t = 0.5, the perturbation u = 0.01 sin(πx) to ρ and D. The results in Fig. 4 illustrate the numerical behaviour for different k D .(Note that for k D = 3, the instability develops very slowly and thus can hardly be seen in the plotted figure).The results confirm our analytical calculations, cf.(3.13).Furthermore, we can see that the form of the new non-constant equilibria and the speed at which they are approached depend heavily on the value of k D .We experience the same behaviour in the two dimensional case.Here, the perturbation is u 2D = 0.01 sin(πx) sin(πy).In Fig. 5, we only show the case k D = 3 as an example.Videos for both the one-and two-dimensional case are available at Videos (Online).

(b) Plateau Solutions
Neglecting the time derivative in (3.2), choosing κ = 10 −4 we are in the situation of Section ii.Starting with a constant initial guess ρ = 0.25, D = 0.75, both perturbed by 0.01 sin(πx) sin(πy), we can numerically confirm the emergence of plateau solutions as described in Sec.ii in one as well as in two space dimensions, cf.

Limitations of the Model: Non-convex Obstacles
In this subsection, we present the results of a numerical test designed to examine the herding behaviour of our model.The basic set-up is shown in Fig. 8(a): A group of people is located behind an obstacle that prevents them from seeing the exit of the room (i.e. S = 0 behind the obstacle, cf.8(b).However, due to the diffusion in the model, after a certain time some people will move around the obstacle and reach a position from which is it possible to see the door (S = 0).Naively, one would expect that other people will follow these people and therefore be able find their way to the exit faster and that this effect increases with the value of k D .To verify this, we did simulations with several values of k D (namely 0, 1, 3 and 5) and compared the loss of mass vs. time, cf.Fig. 9.However, the results shown confirm this expectation neither in the discrete nor in the continuous case.In fact, people behind the obstacle moving around due to diffusion are creating a large D-field within the obstacle.The stronger the coupling, the more people are held back within the obstacle.
Fig 2) given by Boundary ConditionsFor the boundary conditions, we prescribe homogeneous Neumann conditions for D, i.e.∇D • n = 0 on ∂ΩArticle submitted to Royal Society

•
In the Keller-Segel model with volume filling (corresponding to our model with ρ(1 − ρ) replaced by ρ), meta-stable plateau solutions (asymptotically piecewise constant taking the values zero and one) are obtained for small diffusion coefficients in the equation for ρ (cf.Dolak & C. (2005); Burger et al. (2008)) • In the Poisson-Nernst-Planck equations with size exclusion (corresponding to our model with ρ(1 − ρ) replaced by −ρ and δ = 0) stationary solutions are obtained in the presence of additional external potentials for small diffusion coefficients in the Poisson equation (the equation for D, cf.Burger et al. (

Figure 5 .ArticleFigure 8 .
Figure 5. Linear stability in 2 D: ρ (left column) and D (right column) at time t = 0, t = 60 and t = 300 Article submitted to Royal Society