Handling congestion in crowd motion modeling

We address here the issue of congestion in the modeling of crowd motion, in the non-smooth framework: contacts between people are not anticipated and avoided, they actually occur, and they are explicitly taken into account in the model. We limit our approach to very basic principles in terms of behavior, to focus on the particular problems raised by the non-smooth character of the models. We consider that individuals tend to move according to a desired, or spontanous, velocity. We account for congestion by assuming that the evolution realizes at each time an instantaneous balance between individual tendencies and global constraints (overlapping is forbidden): the actual velocity is defined as the closest to the desired velocity among all admissible ones, in a least square sense. We develop those principles in the microscopic and macroscopic settings, and we present how the framework of Wasserstein distance between measures allows to recover the sweeping process nature of the problem on the macroscopic level, which makes it possible to obtain existence results in spite of the non-smooth character of the evolution process. Micro and macro approaches are compared, and we investigate the similarities together with deep differences of those two levels of description.


1.
Introduction. Congestion phenomena in population dynamics cover a wide range of mechanisms, which can be classified according to their stiffness : 1. Soft congestion: as the distance between individuals becomes smaller, the behavior of a single person is affected by the presence of others; 2. Hard congestion: actual contacts between individuals occur, and the overall motion is perturbed by the fact that two persons may not occupy the same place at the same time. Soft congestion models can be included in both microscopic and macroscopic models. In the first class there are discrete-space models like cellular automatabased models which constrain pedestrians to be located at squares of a fixed grid [7,11,46,53,54]. We can also mention models based on networks as route choice models [9,8] or queuing models [40,62]. A lot of evacuation softwares rest on such models as for example buildingExodus in [32]. Some microscopic models are space continous as the social force model introduced in [35] and its forerunner which has been proposed by [30]. In [35], people are identified to particles submitted to the laws of Newtonian mechanics.
As for macroscopic models, soft congestion is also usually favored. Some of these models (see for example [18], [13] and [14]) own their origins in vehicular traffic models, and deal with the one-dimensional case. In higher dimension, similarly to the microscopic case, many models rely on social forces. One solution, used e.g. by Bellomo and Dogbe in [4] and [24], or Degond in [22], is to add to the equation satisfied by the velocity a repulsive term that forces people to avoid high density areas. Another possibility is to modify directly the velocity by adding a term that make people deviate from their preferred path as soon as they approach a crowded area. The modeling of the velocity is the key point of these strategies, see for example the work of Hughes [36], [37], Coscia [19], or Piccoli [47], [48].
We aim here at addressing the particular issues pertaining to hard congestion. In this spirit, we shall make very crude assumptions on the social aspects triggering crowd dynamics (see Section 6 for possible improvements on these aspects), and simply consider that, at any time, a spontaneous (or desired) velocity field is given, and that it is purely reptilian: everyone's priority is to achieve his own goals, without accounting for others. Actual behavior shall result from some sort of compromise between individual tendencies and congestion constraints. More precisely, we shall consider that the actual velocity is defined as the projection of the spontaneous one onto the set of feasible velocities (i.e. which do not lead to a violation of the non-overlapping constraint).
Those basic principles can be applied to microscopic and macroscopic descriptions of the crowd. In the microscopic setting, the degrees of freedom are the positions of individuals (identified to rigid disks), and the non-overlapping constraints can be written straightforwardly by prescribing a minimal value for the distance between centers. The problem takes the form of a differential inclusion which fits into the general framework of sweeping processes introduced by Moreau in the 70's (see [44]). This approach has been extended more recently to non convex cases, namely the case of uniformly prox-regular sets in [16], in [55] and later in [17]. The perturbed (i.e. inhomogeneous) problem has been studied in [10,55,27,28]. As we shall see, it gives a natural framework for the microscopic model we propose.
In the macroscopic setting, the crowd is seen as a population density and the non-overlapping constraint consists in prescribing a maximal value for this density. Although both microscopic and macroscopic models express the same type of modeling assumptions, the mathematical structure of the macroscopic model is less obvious. Expressed in the Eulerian framework as a transport equation by the actual velocity field, which is defined as the projection of the spontaneous one onto the set of feasible velocities, it involves a conservative transport equation by a field whose regularity cannot be controlled a priori, so that classical results cannot be used. A first attempt to address those issues was proposed in [43], in the case where the spontaneous velocity is the gradient of a given dissatisfaction function (typically the distance to the exit for the evacuation of a building): the framework of optimal transportation allows to reformulate the problem as a gradient flow in the Wasserstein space of measures (see e.g. [61] for an overview of the theory of optimal transport, and [3] for the notion of gradient flow in this setting). We show here that the sweeping process framework, which does not rely on any assumption on the gradient nature of the forcing term, can be extended in the macroscopic situation.
Stemming from similar principles, both microscopic and macroscopic models present nevertheless deep differences. For instance we shall point out the fact that the notion of maximal density is not properly defined in the microscopic setting. Thus the macroscopic one cannot be expected to be obtained from the microscopic one by any homogenization process.
We shall focus here on the particular issues related to hard congestion. In order to highlight the very problems it raises, we will favor a very basic form of the crowd motion model. In particular we shall disregard in the main part of this paper any social or strategical aspects of human behavior. We are aware of the very crude character of these assumptions, and we gathered in Section 6 other aspects of pedestrian behavior that could be included in our approach. Besides, we hope that the framework presented here, and the issues we raise in attempting to establish links between microscopic and macroscopic settings, may be fruitful in other contexts, e.g. micro-macro issues in granular flows or modeling of chemotactic motion of large populations of cells.
2.1. Microscopic setting. We consider a population of N individuals, identified with rigid disks centered at q 1 , q 2 , . . . , q N , with common radius r > 0 (it can be extended straighfowardly to the so-called polydisperse case, i.e. with different radii). We denote by the generalized spontaneous velocity of the crowd (U i is the velocity which the individual i would like to have in the absence of others). We assume here that the desired velocity of i depends on his location only, and that this dependence is the same for all individuals (see Section 6 for more general behavioral models): where U 0 is a given field. The set of feasible configurations is defined as We disregard here obstacles (like walls, furniture) to alleviate notations, but they can be included in the definition of K straightforwardly. Non-overlapping is preserved by prescribing that It leads to define the set of feasible velocities as The model simply writes where the projection is performed according to the euclidean norm over R 2N .
Saddle-point formulation. Problem (2) is self consistent and drives the evolution process, as will be shown in Section 2. Yet, the saddle-point formulation of the projection problem sheds light on the underlying Darcy-like structure of the problem, and introduces a pressure field which will have a straightforward interpretation in terms of modeling, and which will have a continuous counterpart in the macroscopic setting.
At some given configuration q, let us denote by Λ the set of couples (i, j), i < j, such that i and j are in contact. Constraints on the velocity can be written in a matrix form (we drop the dependence of matrix B upon the configuration) where the rows of B correspond to active constraints : The saddle-point formulation of the projection problem consists in finding (u, p) ∈ supplemented by the complementarity condition Pressure p ij can be seen as the interaction force between individuals i and j.

Macroscopic setting.
In the macroscopic setting the crowd is represented by a density ρ, which we shall assume to be supported within some domain (the room) Ω. In order to alleviate the notations, we shall present the model in the case of a closed room (see the end of Section 3.2 for some details on the way an exit door can be accounted for in this framework). If we set at 1 the saturation value, the set of feasible densities is The density is advected by the actual velocity u (macroscopic counterpart of the microscopic one dq/dt) ∂ t ρ + ∇ · (ρu) = 0, and u is defined as where C ρ is the set of feasible velocities. It can be defined unformally as the set of all those velocities which have a non-negative divergence on the saturated zone [ρ = 1] (where the density may not increase). More precisely, it is defined by duality as with H 1 + (Ω) = {q ∈ H 1 (Ω) , q ≥ 0 a.e. in Ω}.
Saddle-point formulation. The dual expression of C ρ induces a natural saddlepoint formulation for the projection problem, based on a pressure field p: HANDLING CONGESTION IN CROWD MOTION MODELING   5 with the complementarity condition: − Ω u · ∇p = 0 .
3.1. Microscopic model : generalized sweeping processes. We address here the question of existence and uniqueness of solutions to where t → q(t) ∈ R 2N describes the motion of the crowd, C q is the cone of feasible velocities defined by (1), and U(q) is the desired velocity. Before stating the main results of this section, let us describe the general framework into which it fits, namely the sweeping processes, introduced by Moreau in the late 70's (see [44]). The original setting was the following: consider V a Hilbert space and t −→ K(t) ⊂ V a path of convex closed sets in V , with some regularity in time (e.g. K(t) is continuous for the Hausdorff distance). He was interested in the evolution of a point q(t) subjected to remain in K(t). Assuming that the evolution process tends to minimize the norm ofq, one ends up with the following process dq dt where N K (q) is the outward normal cone to K at q: This normal cone has been introduced in [15] and is called more precisely the proximal normal cone. Note that, as K is convex, this normal cone identifies to ∂I K , the subdifferential of the indicatrix function of K(t). The following discrete process (so called catching up algorithm), which is used by Moreau to establish well-posedness of the problem, gives a clear idea of the evolution mechanism. Consider a time step τ > 0, an initial condition q 0 ∈ K(0), successive positions are built according to q n = P K(t n ) (q n−1 ).
As K(t n ) is closed and convex, the projection is well-defined and contractant, which allows to establish convergence results for the sequence of piecewise constant solutions q τ . It is clear that the recursive process extends straightforwardly to the case where the projection onto K is properly defined (i.e. single-valued) into its neighborhood. The finite-dimensional sets satisfying this property were introduced by Federer in [29] under the name of positively reached sets. Then, they were called p-convex sets by Canino in [12] and later proximally smooth sets by Clarke, Stern and Wolenski in [15]. The final name "uniformly prox-regular set" will be given by Rockafellar et al. in [49,50] (see Definition 3.1 below). The crowd motion model differs slightly from the sweeping process, as the feasible set is fixed whereas the point q (which represents the whole crowd) tends to evolve according to some given velocity U. The same principles can be extended straightforwardly, in particular: assuming U is not too large, and K is uniformly prox-regular, then for a sufficiently small time step τ , the catching up algorithm can be used to build discrete solutions. In the same manner, one may write the evolution process as dq dt − N K (q(t)) ∋ U(q(t)).
Note that this new formulation is a direct consequence of (2): as N K (q) and C q are mutually polar, the identity operator can be written as I = P NK (q) + P Cq (see [45]). Equivalence between both formulations is not obvious, as an identity has been replaced by an inclusion, yet we shall see that it holds true under some conditions. The catching up algorithm reads as follows Step 1 (prediction):q n+1 = q n + τ U(q n ); Step 2 (correction): q n+1 = P K q n+1 . Let us now express in a rigorous manner that this is indeed an algorithm for τ sufficiently small, and that this strategy allows to build a solution to our problem.
Let us first give a proper definition of prox-regularity (see Theorem 1.3 in [50]).
is the proximal normal cone to K at q (defined by (7)). It is said to be uniformly prox-regular with a constant η if it is η prox-regular at every point of its boundary.
Proof. First it can be shown that K ij = {q ∈ R 2N , D ij ≥ 0} is uniformly proxregular. In this case, tools of differential geometry can be used to compute the constant of prox-regularity which corresponds to the smallest radius of curvature (i.e. the largest eigenvalue of Weingarten operator). After calculation, we show that K ij is uniformly prox-regular with constant r √ 2. One may wonder whether the intersection of such sets (which is the case for K) is uniformly prox-regular with a constant depending only on the constants of prox-regularity of the smooth sets. From a general point of view, this is wrong as illustrated in Figure 1. Indeed, we have plotted in solid line the boundary of a set S which is the intersection of two identical disks' complements. This set is uniformly prox-regular but its constant of prox-regularity (equal to the radius of the disk plotted in dashed line) tends to zero when the disks' centers move away from each other. In this situation, the scalar product between the normal vectors n 1 and n 2 (see Fig. 2) tends to -1, which suggests that the prox-regularity constant also depends on the angle between normal vectors.  Consequently the proof of the uniform prox-regularity of K rests on a good estimate of the angles between gradients of active constraints. The existence of such a constant η relies on the positive linearly independence of gradients of active constraints. More precisely it can be shown that there exists γ > 1 such that for all q ∈ K, It can be deduced from this inequality that K is η-prox-regular with η = r √ 2/γ (a detailed proof of this result can be found in [42,58]). Furthermore, the given upper bound is obtained by considering a configuration q of aligned disks and by computing the constant of local prox-regularity at q (see Proposition 3.18 in [60] for more details).
Note that the prox-regularity coefficient of K degenerates (i.e. goes to 0) as the size of the disks goes to zero, for a constant total mass. Furthermore the proxregularity coefficient of K depends on N with N 3/2 for N large enough (see [52]). Now the following result can be established: Let U be Lipschitz and bounded, for all T > 0 and all q 0 ∈ K, there exists a unique absolutely continuous solution t −→ q(t) to Moreover, this solution satisfies the differential equation The well-posedness of the differential inclusion can be proved thanks to results in [27,28]. Then Proposition 3.3 in [5] claims that the solution satisfies the following differential equation which is equivalent to (2) since the cones N K (q) and C q are mutually polar.
3.2. Macroscopic model. We are concerned here with the question of existence of solution to the macroscopic model: given a desired velocity field x → U(x) defined in Ω, given an initial density ρ 0 , find ρ(x, t) such that where the actual velocity field u verifies and C ρ is the set of feasible velocities defined by (5). Eq. 10 is meant in a weak sense: The overall problem can be written as a non-local and nonlinear transport equation ∂ t ρ + ∇ · (ρu(ρ)) = 0, where the notation u(ρ) expresses a dependence of u upon the whole density field ρ, through the projection of U onto C ρ . The nonsmooth character of this projection and the fact that the regularity of u is not controlled (the projection is performed in the L 2 sense), rule out the possibility to use standard tools to define solutions to this equation. One may wonder whether the catching-up approach, which proved successfull in the microscopic setting, can be followed in the present context. The prediction step is straightforward: one transports ρ according to the desired velocity field U during a time step τ > 0. As for the correction (or projection) step, it appears immediately that the standard eulerian way to measure distances between densities (by estimating a given norm of their difference) is not suitable. In the microscopic setting, the catching up approach proved successful because the difference between two configurations identifies with a displacement field, and this is due to the very Lagrangian description of individuals. Recovering this feature calls for a new way to measure differences between densities, more respectful of the Lagrangian description of the motion, and this way is the Wasserstein setting.
Optimal transportation and Wasserstein distance. We assume here that Ω is convex. Let t : Ω −→ Ω be a measurable mapping, and µ ∈ P(Ω) a probability measure. We say that t pushes forward µ onto ν (written for any measurable set A. For given measures µ 0 and µ 1 , we denote by Π(µ 0 , µ 1 ) the set of transport maps between µ 0 and µ 1 . The quadratic Wasserstein distance Notice that this is a sloppy definition that holds for atomless measures, which is anyway the case we are interested in; for general measures one should pass through the so-called transport plans, which we do not want to introduce here, and we address the interested reader to [3,61]. However, this quantity W 2 can be proven to be a distance on the space P(Ω) and it makes the space of probability measures a geodesic space, i.e. every pair of points is linked by a curve (which is in this case a curve of measures), such that the length of this curve exactly equals the distance between the points. In case in the minimization defining W 2 there is an optimal transport map t (which is the case if µ 0 is absolutely continuous with respect to the Lebesgue measure L d ), then this geodesic curve is known explicitly and it is given In the spirit of the structure that we are imposing on P, we define the outward normal cone to a subset K ⊂ P (we actually apply the definition of the strong Fréchet subdifferential in [3] for all t such that t # ρ ∈ K. It makes it possible to express the problem like we did in the microscopic setting (Eq. (8)): and u transports ρ according to (10).
We may now write the catching up algorithm (which we consider for the time being as a theoretical tool): given an initial density ρ 0 and a time step τ > 0, build ρ 1 , . . . , ρ n according to where the projection is performed in the Wasserstein sense and the set K is the constrained set introduced in (4). As for the microscopic model, we must check that both steps are well defined (possibly under some restriction on the time step τ ), and that the obtained discrete solutions converge to a limit which can be identified as a solution to problem (10)(11). We shall assume here two properties of U. First we require that it tends to keep people within the room Ω, i.e. that for τ sufficiently small (id + τ U) maps Ω onto Ω. Obviously, in the case of an emergency evacuation, with an open exit door, we will of course alleviate this assumption by allowing U to cross the exit line (we will discuss later on how to "catch up" the possible part of the mass that exits Ω because of this). Secondly, we need it regular enough (Lipschitz continuous is sufficient), so that, still for τ sufficiently small (id + τ U) preserves the absolute continuity of the measure, i.e. ρ < < L d implies (id + τ U) # ρ < < L d . This assumption being made, the prediction step is well defined for any measurable field U.
The key point is the projection step, which differs significantly from the microscopic situation. It is a variational problem in the space of measures, which consists in minimizing the distance to a fixed measure among elements of K. As for existence, standard compactness arguments may be applied, but the question of uniqueness is more delicate.
First of all, let us rule out the possibility to establish uniform prox-regularity in the spirit of Definition 3.1. Consider the one-dimensional situation. A Dirac mass at zero projects onto the characteristic function of (−1/2, 1/2). More generally, a combination of Dirac masses µ = α n δ xn , α n = 1, (we assume that the distance between supports of any two of them is always larger than half the total mass carried by the couple, so that they do not interact) projects onto the sum of characteristic function of the flattened Dirac masses: The Wasserstein cost can be computed straightforwardly, it is On the other way around, consider a density ρ ∈ K (i.e. ρ(x) ≤ 1 almost everywhere) and ω the largest open set such that ρ(x) = 1 a.e. in ω. This open set ω is a countable union of open intervals of lengths (α n ) 1≤n<N (with possibly N = +∞). Among all densities which projects onto K at ρ, the farthest is the combination of Dirac masses with weights α n , supported at the middles of intervals, and it is at distance ( α 3 n /12) 1/2 . As a consequence, η prox-regularity can be expected at some points, but η is not bounded away from 0. Note that some measures are the projection of themselves only, even if they saturate the constraint at every point of their support: consider e.g. a dense open set of total measure 1 in (−1, 1), and define ρ as the characteristic function of its complement. As the saturated zone contains no interval, it cannot be the projection of a measure which violates the constraint.
Note that this non prox-regularity can be seen as a natural consequence of the fact that the property in the microscopic setting degenerates as the granularity gets finer (i.e. when the radius goes to 0).
On the other hand, the set K enjoys some kind of convexity property: it is geodesically convex (see [3]) : Definition 3.3. K ∈ P(Ω) is said to be geodesically convex if for any µ 0 , µ 1 ∈ K, the geodesic curve µ t joining µ 0 and µ 1 belongs to K for all t ∈ [0, 1].
This property might suggest that uniqueness of the projection can be obtained by standard arguments: in a so-called Aleksandrov non-positively curved (NPC) length space (see [1]), considering 2 minimizers and the measure halfway along the geodesic, convexity properties of the distance function µ → W 2 (ρ, µ) 2 lead to a contradiction.
Unfortunately, as soon as d ≥ 2, P(Ω) is not NPC and the previous distance function even turns out to have concavity properties.
Yet, despite the previous assertions (non prox-regularity and the concavity of the squared distance), other kind of interpolation curves (called generalized geodesics) can be used to assert well-posedness of the projection problem, without any restriction on the distance from K. Definition 3.4. (Generalized geodesics, see [3], p. 207) A generalised geodesic joining µ 0 to µ 1 with base ρ is a curve defined by where r 0 (resp. r 1 ) is an optimal transport map from ρ to µ 0 (resp. µ 1 ).
Notice that when r 0 = id this gives again the expression of a geodesic in P according to the distance W 2 .
This generalized notion can be used to establish the following proposition Proposition 2. For any ρ ∈ P(Ω), the Wasserstein distance to the set K of admissible densities defined in (4) is attained at a unique point P K ρ. The projection operator P K is continuous.
Now that the catching-up algorithm is defined properly, the obtained sequence of discrete solutions allows to obtain an existence result: Theorem 3.5. Let ρ τ be the piecewise constant interpolation of the discrete densities given by the catching-up algorithm (13). If U is a C 1 velocity field, and ρ 0 ∈ K, then ρ τ converges as τ tends to 0 to a solution of the macroscopic problem Proof. Let us first describe the discrete quantities we obtain thanks to the catching up algorithm. We denote by r n+1 the optimal transport between the projected density ρ n+1 = P Kρ n+1 andρ n+1 , and we write t n+1 = (id + τ U) −1 (see figure 3). These transport maps allow to define a discrete velocity as follows: We then define two different interpolations of these quantities. First, the piecewise constant interpolation is given by : if t ∈]nτ, (n + 1)τ ].
Using [43] for the projection part, it is possible to prove that v τ satisfies the following discrete decomposition U = v τ + ∇p τ + ε τ , where ε τ converges uniformly to 0 as τ tends to 0. We also define a continuous interpolation of (ρ n ) n as follows It is possible to prove a priori estimates on these interpolated curves, and therefore prove that they both converge to the same limit. The continuous interpolation of (ρ n ) n gives that the limit density satisfies a transport equation, and the discrete decomposition of the piecewise constant interpolation proves that the limit velocity is indeed the projection of the desired velocity onto C ρ .
We finish this section by considering the case where, due to modeling reasons, we allow the vector field U to let the mass exit through a part of the boundary. This means that id + τ U is no longer supposed to map Ω into Ω but the segment connecting x to x + τ U(x) is allowed to cross ∂Ω in a prescribed subset Γ ⊂ ∂Ω which stands for the exit. In such a case the transport step of the algorithm does not need to be changed, but we have to face, during the projection step, a densitỹ ρ n+1 which is no longer supported in Ω. To bring back the mass to the original domain we want to consider a modified set K, that we call K Γ : The reason for this choice is the following: we consider that as soon as a particle reaches Γ, instead of following its movement after Γ, we leave it on it. This is done for simplicity, but it only means that we are no longer concerned with what happens to the particles that have reached Γ, not that they are really blocked on the exit. Obviously, we needed to withdraw the density constraint on Γ, so as to let particles stay on it, and also to represent the fact that Γ stands actually for everything that happens at the door and beyond.
Once the set K Γ is introduced, the projection step is done by defining which means that we take a measure whose support may go beyond Ω and we project it onto the measures over Ω satisfying some extra density constraint. We stress that the same kind of argument (projecting on a set of measures concentrated on Ω) could also be used, in the case with no exit, to withdraw the (yet, natural) assumption that id + τ U maps Ω into Ω.
The mathematical problem with this modified K Γ is much trickier. One of the difficulties, that prevent the usual theory to be applied, is the fact that the set K Γ loses some of the properties that K had previously (in particular it is no more geodesically convex: the geodesic -for the W 2 distance -between two points of K Γ could go out of K Γ ). In particular we have no clue about the uniqueness of the projection, but the algorithm works in the same way if we accept to take any minimizer of the distance.
3.3. Gradient flow setting. In case the spontaneous velocity is the gradient of some dissatisfaction function (e.g. distance to the exit in case of an emergency evacuation), the microscopic model can be put into a gradient flow form: where ∂ϕ is the Fréchet subdifferential of ϕ defined by The function Ψ that we consider is the total dissatisfaction, i.e. Ψ(q) = i D(q i ), where D(x) may be, as we said, the distance to the exit. Indeed, as K is uniformly prox-regular, the proximal normal cone identifies with the Fréchet subdifferential of I K hence ∂ϕ = ∇Ψ + N K (q) (see [51] for more details).
In the macroscopic case, under the same assumption that U is a gradient −∇D, then one defines the global dissatisfaction function in a continuous setting: and the overall process can be seen as a gradient flow in the Wasserstein space, i.e. ρ is advected by u, where ϕ = Ψ + I K , and ∂ϕ is defined in the following sense (see [3]): Definition 3.6. The strong Fréchet subdifferential of a function ϕ at ρ is the set of fields u such that for all transport maps t, the following inequality holds true Notice that the definition we gave before of N K is nothing but a particular case of this one, when ϕ = I K .
The main advantage of this gradient-flow formulation is the fact that it allows to handle vector fields U less regular than what we need in the general case (where it is natural to require U to be Lipschitz continuous, i.e. D ∈ C 1,1 ). For the whole theory on gradient flows, first in a finite-dimensional setting, then in Hilbert spaces, and finally in metric spaces and particularly in P(Ω), we refer again to [3].
In particular, the catching up algorithms made by the coupling of a prediction and a correction step, may be replaced in the case of a gradient structure by a single-step procedure, called proximal algorithm, where This algorithm has also been formalized in a metric setting, where it is known as minimizing movements (see [21,2]): it has first been applied to the case of P(Ω) with the Wasserstein metric by Jordan, Kinderlehrer and Otto in [38]. The details about a gradient flow approach to the macroscopic framework of crowd motion are contained in [43], where both the case with no exit (i.e. with the constraint ρ ∈ K) and with exit (ρ ∈ K Γ ) are dealt with, the second one being much trickier than the first.

Numerical solution.
4.1. Microscopic model. Moreau's catching up algorithm suggests a strategy to build approximate solutions to the evolution problem. Yet, projecting a configuration q onto K is not straightforward. The scheme we propose extends ideas introduced in [41] for granular flows. It consists in performing a catching up step with K replaced by some kind of inner local convex approximation. More precisely, considering that the configuration q n at time t n is known, the predicted configuration is obtained byq n+1 = q n + τ U(q n ).
Then q n+1 is obtained by projectingq n+1 onto K q n q n+1 = P K q nq n+1 , where K q stands for In Figure 4, we illustrate the set K ⊂ R 2N , intersection of sets K ij whose boundaries are plotted in solid line. The set K q n is delimited by the dashed line. The theoretical and numerical projections, respectivelyq n+1 := P K (q n + τ U(q n )) and q n+1 are represented (for two examples of U(q n )). Indeed, as K is uniformly prox-regular, the projection onto K of q n + τ U (q n ) is well-defined for τ small enough. The replacement of K by the convex set K q n is convenient because it allows us to use classical numerical methods to compute this projection. Indeed, this projection problem can be reformulated in a saddle-point form which can be solved by Uzawa algorithm : find (q n+1 , λ) ∈ R 2N × (R + ) Note that the Lagrange multipliers are not unique in general. Indeed there exist configurations q n with strictly more than 2N active constraints so that the associated gradients G ij (q n ) are linearly dependent. The matrix appearing in Uzawa algorithm is C = B t B where B is the matrix whose rows are vectors G ij (q). In [42], the authors quantify how the condition number of matrix C varies with the parameter η q (setting a lower bound of the local prox-regularity of K at point q) when C is non-singular. As a consequence we expect that the Uzawa algorithm converges less quickly for configurations with low local prox-regularity. In numerical simulations, we noticed indeed that solving the saddle-point problem requires more iterations in case of a jam.
Numerical analysis. In the proposed scheme (16), the replacement of K by the convex set K q n is computationally convenient but arises many difficulties in the numerical analysis. The convergence result is based on the fact that K q is a good approximation of K near q. Theorem 4.1. Let us denote by q τ the continuous piecewise linear function satisfying q τ (t n ) = q n defined by (16). Let U be Lipschitz and bounded, for all T > 0, the sequence (q τ ) τ uniformly converges to q defined in Theorem 3.2, when τ tends to 0.
Proof. Here we just give a sketch of the proof. We refer the reader to [58] for the details. It can be easily proved that (q τ ) τ is bounded and so a convergent subsequence can be extracted. Thanks to the uniqueness result in Theorem 3.2, it suffices to check that the limit function (which is absolutely continuous) satisfies the differential inclusion. Since q n+1 = P K q n (q n + τ U(q n )), it comes for allq, By dividing by τ , we obtain for allq where u n = (q n+1 −q n )/τ represents the discrete actual velocity. We aim at passing to the limit in the previous inequality and the crucial term is the last one. Thus the convergence result rests on the local continuity of the map (q, q 0 ) −→ d Kq 0 (q) in a neighborhood of the set {(q, q 0 ), q = q 0 }. More precisely it can be checked that for all q ∈ K and allq ∈ B(q, r q ), by using the saddle point form of these projections (17). Thanks to the reverse triangle inequality, the Lagrange multipliers (which are not unique in general) are bounded. Then compactness arguments allow us to obtain the convergence of the projections. We deduce that the limit function q satisfies the following differential inclusion: As for all q ∈ Q, N (K q , q) = N (K, q), the required result is proved. The convergence order of this scheme will be specified in a forthcoming paper [6]. This more accurate result is based on metric qualification conditions between sets K ij q = {q , D ij (q) + G ij · (q − q) ≥ 0}. More precisely there exists a constant θ > 0 such for all q ∈ K and allq ∈ B(q, r q ), . This result rests on the reverse triangle inequality (9).

4.2.
Macroscopic model. The strategy we propose relies again on the catching up philosophy which already made it possible to establish an existence result. The density is first advected by the spontaneous velocity field (with no consideration of the congestion constraint), and then projected (in the Wasserstein sense) onto the set of feasible densities. The time discretization scheme writes as follows: Given a time step τ > 0, an initial configuration ρ 0 , approximate configurations ρ 1 , . . . , ρ n are built recursively according tõ The first step consists in transporting ρ n according to the given transport map id + τ U. We simply transport the center of each cell at velocity U during a time τ , and then distribute the mass of the transported cell on the different cells of the mesh it intersects, as illustrated in Fig. 5. The transported density writes where A(B) represents the area of the set B. The second step is less standard. The scheme we propose to approximate the id + τ U Figure 5. Lagrangian transport of the density projection onto K with respect to the Wasserstein distance is based on the following considerations 1 : Consider a domain ω and a non-zero density µ supported by ω. For ε > 0, 1 ω + εµ / ∈ K, and its projection onto K is identically 1 in ω. Let us denote by εu the displacement field which corresponds to the optimal transport between 1 ω + εµ and P K (1 ω + εµ). One has (id + εu) # (1 ω + εµ) = P K (1 ω + εµ).
As id + εu is an optimal map, u is a gradient : u = −∇p, where p verifies the Poisson problem −∆p = µ.
Since we are considering a violation of the density constraint on a set ω and assuming that, after the projection, the density will be saturated on ω itself, we are only interested in finding the mass that will exit ω when the displacement is given by εu. This means, at a first order approximation, we only need to estimate the flux of density through the boundary ∂ω, i.e. u · n = −∂p/∂n. Now consider the stochastic interpretation of the Poisson equation Considering that µ is a probability measure, we consider a random variable X ∈ ω which follows the law of probability with density µ. We now consider a Brownian motion stemming from X. Its path crosses ∂ω for the first time at Y . The random variable Y ∈ ∂ω is known to follow the probability law with density −∂p/∂n on ∂ω (see [26]). The idea is therefore to redistribute the exceeding mass of the saturated zone in the following way (see Fig. 6): For each saturated cell, a random walk is started, which transports the exceedind mass m = (ρ − 1)|C|, where |C| is the measure of the cell. When this random walk encounters a non-saturated cell, it gets rid of as much mass as it can, and continues as long as the transported mass is not fully distributed. When all the saturated cells have been treated, the obtained density ρ n+1 is admissible. Let us notice that the gradient flow structure when U = −∇D would also allow for another time-discretized algorithm, namely the proximal (or minimizing movement, or Jordan-Kinderlehrer-Otto) one instead of the prediction-correction one. This would lead to a sequence of minimization problems involving the computation of the distance W 2 . The formulation by means of the transport plans, that we did not develop here for the sake of readability, transforms these problems into linear programming problems, that could a priori be solved through a simplex algorithm. Yet, it turns out to be so slow that 2D problems could not be efficiently solved in this way.
Numerical analysis. The numerical analysis of the transport part of the algorithm is standard, but the projection part is more delicate, and we are not able to provide a rigorous error estimate for the stochastic scheme we propose. Let us simply say here that, beyond the unformal justification of the approach given previously, the asymptotic behavior of the so-called Diffusion Limited Aggregation (DLA) process (see [39]), and the link between the distance based on the energy norm between measures (see [31] or [26]) and the Wasserstein distance give some arguments to support the chosen strategy. 5. From micro to macro. Both microscopic and macroscopic model express the same assumptions at their respective levels: a spontaneous velocity is given, and the system evolves according to a velocity which is the closest to the spontaneous one among all feasible velocities, in a least square sense. Yet, the macroscopic model was not built from the microscopic one by a rigorous homogenization process. We describe here some obstacles to such a process, to shed light on the deep differences between both settings, in spite of formal similarities.

Maximal density.
First of all, let us point out that the notion of maximal density is somewhat ambiguous as far as the microscopic model is concerned. Let us consider the case of identical radii (monodisperse situation). The maximal packing density for identical disks is known to be ρ max = π/2 √ 3 ≈ 0.9069 . . . , and corresponds to the triangular lattice (see Fig. 7). Yet the actual density of moving collections of rigid disks is Figure 7. Triangular lattice generally strictly less than this maximal value, which is only attained for this very particular configuration. As an example, in the evacuation situation represented in Fig. 8, the mean density upstream the exit door ranges typically between 0.85 and 0.87. The fact that the actual density is strictly less than the maximal one does not mean that the flow is unconstrained (as the macroscopic setting would suggest). Those considerations call for a clear identification of configurations which saturate the constraint. Such a notion is proposed for the three-dimensional situation in [57], in a Nash equilibrium spirit: We say that a particle (or a set of contacting particles) is jammed if it cannot be translated while fixing the positions of all of the other particles in the system. Note that this property is defined as local jamming in [57]. It is tempting to consider as maximal in some sense any density corresponding to such configurations, for which there are no free disks, so that constraints are activated everywhere. The triangular lattice is clearly jammed, but so is the cartesian lattice (ρ = π/4 ≈ 0.79), and it is possible to build looser jammed configuration (see Fig. 9, with ρ = π √ 3/8 ≈ 0.68).  Constraints on the velocity. Beyond this fuzzy definition of maximal packing, the set of feasible velocities in the microscopic setting is strongly dependent on the local structure and not only on the density, a feature which is lost in the macroscopic approach. Let us first consider the triangular lattice represented in Fig. 7. Having the number of disks going to infinity and the radius going to zero, the density (characteristic function of the solid phase) converges to the uniform density ρ max . Microscopic feasible velocities can be defined in the solid phase, and one may wonder what are the corresponding feasible velocities in the macroscopic limit. Such velocities are clearly constrained in 3 directions. Using self-evident notations e 0 , e π/3 and e 2π/3 to design the principal directions of the lattice, one obtains at the limit some sort of unidirectional expansion constraints in those directions, which can be written e · ∇u · e ≥ 0 , with e = e 0 , e π/3 or e 2π/3 .
It still allows for non-rigid motions at the limit: Fig. 10 represents a feasible velocity field at the microscopic level, whose macroscopic counterpart is Note that, although the microscopic velocity field is tangent to the boundary of K, the corresponding macroscopic field is strictly expansive (∇ · u > 0).

Figure 10. Example of feasible velocity
Jams.
The considerations which we presented above have a crucial consequence on the behavior of both approaches: the microscopic model has the ability to reproduce blocked jams (which are observed in practice), whereas the macroscopic one does not. The latter property is a consequence of the maximum principle. Let us consider the situation represented in Fig. 11, with a saturated zone (ρ ≡ 1) upstream the exit door. The desired velocity U is assumed to point to the exit door, so that ∇ · U < 0. The actual velocity is U − ∇p, where p solves a Poisson problem in the saturated zone, with homogeneous Neumann condition on the walls, and homogeneous Dirichlet boundary conditions at the exit and on the interface with the non-saturated zone. By virtue of the maximum principle, p is nonnegative over the saturated zone, so that the velocity correction −∂p/∂n is non-negative on the door: people exit quicker than they would if there were no congestion (they are always pushed out by other people behind them). As a consequence, if the desired velocity field tends to move people out of the room, the evacuation process will never stop before the room is empty. The microscopic situation is quite different. The pressures are still nonnegative by definition of the saddle-point formulation (3), but they are likely to act in the "wrong" direction for the individuals which are the closest to the door, as soon as they form an arch (see Figure 12).

Figure 12. Microscopic evacuation (arches)
Mathematical issues (convergence or non-convergence). The precise way for setting convergence questions when passing from the microscopic to the macroscopic model needs to associate a density to every microscopic configuration. An easy way to do that is to associate to every non-overlapping disks configuration the uniform measure (with unit density) on the union of those disks. In this way we obtain a density ρ obviously satisfying ρ ≤ 1 but we also know that, when we let the size of the disks go to zero (with the number of disks increasing consequently to infinity) we cannot obtain any density ρ ≤ 1 in the macroscopic limit. In particular in dimension 2 it is impossible to go beyond 0.91, since we know that the maximal density is realized by the triangular lattice.
This suggests that it is better to rescale the approximating ρ by the maximal density, which is 1 in dimension one, π/2 √ 3 in dimension two. . . Now the question is the following: take a sequence of initial data ρ 0 N , i.e. the densities associated (in the way we described above) to a sequence of microscopic configurations with N non overlapping disks of radius r. Let N → ∞ and N ≈ r −d and consider the limit density ρ 0 . Can we say that, for fixed velocity field U , the constrained evolution stemming in the microscopic model from ρ 0 N converges to that obtained in the macroscopic one from ρ 0 ?
The answer can be expected to be positive if d = 1, but the situation is much more complicated for d ≥ 2. Actually, the unidimensional case gives no ambiguity between jammed and maximal density configuration, which is not the case in higher dimension.
For instance one can consider a sequence of jammed cartesian lattices in a twodimensional domain Ω (a square, for instance), with a concentrating vector field U (for instance U(x) = −x). In this situation, the configuration is completely blocked, and the evolution gives, for any time t > 0, the constant density ρ 0 N . On the other hand, the limit as N → ∞ is a constant density, but it does not activate the constraints. This is due to the fact that the density realized by the cartesian grid is strictly smaller than that of the triangular lattice, which was taken as a reference. Hence, in the evolution, ρ t would differ from ρ 0 and move towards a more concentrated configuration.
But this (jammed configuration which are not of maximal density) is not the only source of problems in the micro-macro limit. Actually, one can consider a similar example with a sequence of triangular grids. In such a case, both ρ 0 N and the limit ρ 0 activate and saturate the constraints. This induces some constraints on the admissible velocities but, as we saw in the paragraph devoted to these constraints, they act differently in the microscopic and in the macroscopic models. In particular, the limits of the admissible velocities for the microscopic problem are those vector fields which are "unilaterally incompressible" in the three main direction of the lattice, which is not at all the same as imposing only a positive divergence. The actual limit constrained space at the macroscopic level necessitates obviously to account for the microscopic structure of the contact network. Beyond packing fraction, different types of order metrics have been introduced in [25] to give a better description of the local structure (see also [57] for a recent review on this matter), but the use of those concepts for evolution problem is still widely open.
This shows some deep differences, for d ≥ 2, between the limit of the microscopic case and the macroscopic one. In the case of a gradient vector field U we could express this issue, thanks to the gradient-flow formulation, in terms of the associated dissatisfaction functionals. Some general results give conditions to translate a variation form of convergence on these functionals (called Γ−convergence) into a convergence of the associated flows, but we do not want to insight this question in details because of its complexity; it would go out of the scope of this paper and it is moreover matter of an ongoing work.
6.1. Strategies. We specify here how the model we presented can be improved in order to account for more elaborate individual behaviors than the purely reptilian one on which we built our approach.
Microscopic setting. We aim here at integrating the fact that individuals are likely to elaborate complex strategies to escape a building. For example, in case of congestion, they may decelerate or try to avoid the jam, instead of keeping pushing inefficiently. The velocity of a person becomes then dependent upon the position of people he can see in front of him. Such a strategy can be defined as follows: We define the set N i (see Fig. 13(a)) containing persons who are near and visible to the individual i: where d i = U 0 (q i )/|U 0 (q i )| and e ij = (q j − q i )/|q j − q i |. The constant α is taken equal to the half angle of view ( i.e. α ≃ 60 • ). Two choices are possible for the individual i if his neighbours (belonging to N i ) walk slower than him.
First, he can decelerate instead of going through the crowd. In this case, his speed s n i (i.e. his velocity's norm) at time t n is made dependent upon his neighbours' behavior at time t n−1 . More precisely, it is computed as a barycenter of his neighbours' speeds, weighted by their relative positions.
Otherwise, he can also be in a hurry and prefer changing his way instead of slowing down. In that case, if there exists another clear way through the set N i , he follows it while keeping his desired speed, or else he will go round the group N i .  Among all directions allowing him to do so, he chooses d new i the closest to the one he wanted (see Fig. 13(b)).
Numerical simulations with these strategies are presented in [59]. If every pedestrian prefers decelerating, there is typically no contact between people and rarefaction waves are observed. On the contrary, the avoiding strategy leads to situation where the constraints are activated.
Macroscopic social forces. Integrating strategies in the macroscopic setting is more delicate, as it is less natural to follow the motion of a single individual in the crowd. However, it is possible to model the influence of local density on the behavior of the crowd by changing the desired velocity. More precisely, it is natural to assume that when people arrive upstream a crowded area, they tend to decrease their desired velocity. In the spirit of [36,23], a two dimensional generalization of road traffic models can be integrated in the model: at time t, people located at x ∈ Ω have a desired velocity U ρ (x) with the same direction as the initial desired velocity U(x) and a norm that decreases when ρ increases: where α(ρ) goes monotonically from 1 to 0 as ρ goes from 0 to 1 (saturation value). The evolution of the crowd then obeys the same principles as before, namely: It can be proved that the density transported by U ρ stays in K: the projection step in the evolution equation becomes useless. Yet, if α(1) > 0 (people continue to push even at saturation density), the congestion constraint is likely to be activated, and the approach we propose in this paper can be adapted to this situation. Note that, from a modeling standpoint, it is natural to use downstream information to determine the desired velocity (people adapt their velocity according to what they see, i.e. to the density in front of them, which is downstream their desired direction), whereas, once the velocity is determined, the transport can be performed using an standard upwind scheme. structure: it is sufficient to alleviate the assumption that the desired velocity depends on the position only, but can be defined differently depending of individuals.
Note that, in case U 1 and U 2 are both gradients of dissatisfaction functions (one for each population), the system presents a gradient-flow structure in the Wasserstein setting (see [20] for a similar problem in the context of cell migration).
6.3. Nash equilibrium approach. We give here some hints on another possibility to define the actual velocity, in a Nash equilibrium spirit. Note that this approach is likely to change the mathematical nature of the evolution problem. In particular, the instantaneous velocity is no longer defined in a unique way, but as any (among infinitely many) equilibrium in the following sense: The actual generalized velocity u = (u 1 , u 2 , . . . , u N ) is such that, for any i, u considered as a function of u i only minimizes the distance to U among feasible velocities, i.e.
u ∈ argmin v∈C i q (u) where C i q (u) is the set C i q (u) = v ∈ R 2 , (u 1 , . . . , u i−1 , v, u i+1 , . . . , u N ) ∈ C q . As soon as two disks are in contact and push against each other, there are infinitely many solutions to this problem, as illustrated by the situation of Fig. 14: the desired velocity of person 1 is 1 (in the horizontal direction), whereas person 2 tends to stay still. In the previous approach (ℓ 2 projection onto the cone of feasible velocity), the actual velocity is 1/2 for both. In the present approach, any diagonal couple (α, α), with α ∈ [0, 1], realizes an instantaneous Nash equilibrium. 1 2 Figure 14. Two disks, Nash approach.
In the macroscopic framework, the problem can be described as follows: For any ω ⊂ R 2 , one prescribes that u |ω ∈ argmin v∈C ω ρ (u) with C ω ρ (u) = v ∈ L 2 (Ω) , v ⊕ u |R 2 \ω ∈ C ρ As indicated previously, this approach changes the mathematical nature of the evolution problem, which calls for further developments. Let us simply say that the solutions we built in the previous sections are particular solutions to this new problem, among infinitely many others.

Numerical experiments.
Micro-macro comparison. The first problem of the comparison between these models lies in the initial configuration: given an initial configuration for the disks in the microscopic setting, it is difficult to choose the macroscopic density that best fits this configuration. As long as uniform densities are considered, we adopt the following method: we first estimate the mean density in the initial microscopic configuration, and then normalize it with the maximal density of disks we encounter throughout the microscopic evolution. In Fig. 15, we present an example of such a computation. Figure 15. Calculation of the macroscopic density: the initial density 0.53 (on the left) is normalized by the maximal density encountered 0.81 (on the right) to obtain a macroscopic density of 0.65.
We then let evolve both systems and compare the configurations at equivalent time steps. The formation of the saturated zones is quite the same in microscopic and macroscopic settings (see Fig. 16). However, as pointed out in Section 5, the behavior of the two models at the exit is very different, in particular the evacuation is faster in the macroscopic model (see Fig. 17).  Pressure field. In the microscopic model, the pressure exerted between people appears naturally as a Lagrange multiplier, and is computed using Usawa algorithm. The macroscopic pressure field does not appear explicitly in the numerical scheme, but it can be recovered by estimating during the projection step a discrete equivalent of the odometer function (see Section 4.2): on each cell, we define the pressure as the total mass emitted by this cell during the stochastic projection. In Fig. 18, we represented the microscopic pressure and the isovalues of the macroscopic pressure field upstream an exit door. In both settings, the pressure field is maximal in the middle of the saturated zone, and decreases at the entrance and the exit of this zone.
Realistic geometries. Let us finally present some examples which correspond to more realistic evacuations situations. The microscopic setting, as suggested in Section 5, describes more properly situations where people leave through narrow exits. Fig. 19 shows the evacuation of the first floor of a large building (Departement of Mathematics at Orsay).
The macroscopic setting, on the other hand, would best fit situations where many people have to evacuate a large domain. We present in Fig. 20 the evacuation of the Stade de France. Even if the benches of the stadium are quite narrow passages, the macroscopic model gives results that correspond to real evacuation configurations.