Noname manuscript No. (will be inserted by the editor) A Semi-Lagrangian scheme for a modified version of the Hughes model for pedestrian flow Elisabetta Carlini · Adriano Festa · Francisco J. Silva · Marie-Therese Wolfram 6 Received:date/Accepted:date 1 0 2 Abstract In this paper we present a Semi-Lagrangian scheme for a regu- larized version of the Hughes model for pedestrian flow. Hughes originally b e proposed a coupled nonlinear PDE system describing the evolution of a large F pedestriangrouptryingtoexitadomainasfastaspossible.Theoriginalmodel correspondstoasystemofaconservationlawforthepedestriandensityandan 5 1 Eikonal equation to determine the weighted distance to the exit. We consider this model in presence of small diffusion and discuss the numerical analysis of ] the proposed Semi-Lagrangian scheme. Furthermore we illustrate the effect of A small diffusion on the exit time with various numerical experiments. N . Keywords Crowd motion · mean field models · Semi-Lagrangian schemes h t Mathematics Subject Classification (2000) 35Q91 · 65N75 · 60J20 a m [ 1 Introduction 2 v In the last decades crowd dynamics has attracted the attention of many re- 4 searchersinthescientificcommunity.Startingfromthefieldofappliedphysics 2 3 7 E.Carlini 0 DipartimentodiMatematica“G.Castelnuovo”,SapienzaUniversit`adiRoma, E-mail:[email protected] . 1 A.Festa 0 RICAM–JohannRadonInstituteforComputationalandAppliedMathematics,Austrian 6 AcademyofSciences(O¨AW),E-mail:[email protected] 1 : F.J.Silva v XLIM-DMIUMRCNRS7252Facult´edesSciencesetTechniques,Universit´edeLimoges, Xi E-mail:[email protected]. r M-T.Wolfram a Mathematics Institute, University of Warwick, Coventry CV4 7AL and RICAM – Johann RadonInstituteforComputationalandAppliedMathematics,AustrianAcademyofSciences (O¨AW),E-mail:[email protected] 2 ElisabettaCarlinietal. andtransportationresearch,themotionofpedestriancrowdsraisedmoreand more interest in the applied mathematics community. Mathematical models range from the microscopic level, where the individual dynamics are described separately, to the mesocopic and macroscopic level, where the distribution with respect to their velocity and/or position in space is considered. Microscopic models are either force-based, such as the social force model pro- posedbyHelbingandco-workers[30]orlatticebasedlikethecellularautomata modelsproposedin[13,9].Onthemacroscopicleveltheevolutionofthepedes- trian density is usually described by a conservation law, see for example [32, 20,41,21,25]. In these models the velocity field may depend on the current local density, a given external potential and physical constraints due to walls and/orbarriers.Recentlymeanfieldgames,cf.[31,35],havebeenproposedto model the evolution of large pedestrian crowds, see [34]. These models can be derived from stochastic optimal control problems for multi-agent systems as thenumberofindividualstendstoinfinity.Foradetailedoverviewondifferent modeling approaches in pedestrian dynamics we refer to [6,23]. In2002R.Hughesproposedamacropscopicmodelforpedestriandynamicsin [32], which is based on a continuity equation (describing the evolution of the crowddensity)andanEikonalequation(givingtheshortestweighteddistance to an exit). It is given by ∂ m(x,t)−div(m(x,t)f2(m(x,t))∇u(x,t))=0,  t 1 (1) |∇u(x,t)|= ,  f(m(x,t)) where x∈Ω denotes the position in space, t∈(0,T], T ∈R the time and ∇ + thegradientwithrespecttothespacevariablex.Thefunctionmcorresponds to the pedestrian density and u the weighted shortest distance to a target, for example an exit. Hughes proposed different functions penalizing regions of high density, the simplest choice being f(m)=1−m where 1 corresponds to the maximum scaled pedestrian density. In this work, we will assume that f is a general smooth function. System (1) is a highly nonlinear coupled system of partial differential equa- tions.Fewanalyticresultsareavailable,allofthemrestrictedtospatialdimen- sion one. The main difficulty comes from the low regularity of the potential u(x,t),whichisonlyLipchitz-continuous.Forexistenceanduniquenessresults of a regularized problem in 1D and the corresponding Riemann problem we refer to [26,2,3]. In this work we consider a modified version of (1), which served as the basis for the 1D analysis presented by Di Francesco et al. in [26]. It corresponds to  ∂ m(x,t)−ε∆m(x,t)−div(m(x,t)f2(m(x,t))∇u(x,t))=0,  t 1 (2) −ε∆u(x,t)+ 1|∇u(x,t)|2 = .  2 2f2(m(x,t))+δ in Ω×(0,T). Theregularizationparameterδ >0preventstheblow-upofthecostwhenap- ASLschemefortheHughesmodelforpedestrianflow 3 proachingthemaximumdensityone.Thediffusivetermsallowtousestandard analytical techniques from nonlinear PDE theory, see [26]. Diffusive phenom- ena have been observed and studied in pedestrian dynamics [47,36], giving an additional justification of the modification considered. System (2) has to be supplemented with suitable boundary and initial condi- tions. We consider an initial density m of the agents satisfying that m ≥0, 0 0 m ∈L∞(Ω) and the support of m is a subset of Ω. Note that rescaling the 0 0 density m , and possibly modifying the function f in the equation, we can 0 (cid:82) assume that m (x)dx = 1. This normalization is useful in order to pro- Ω 0 vide a probabilistic interpretation of the Fokker-Planck (FP) equation in (2). Possible boundary conditions for the pedestrian density m at the exit are: – a given fixed outflow, corresponding to Neumann boundary condition, – anoutfluxwhichdependsonthepedestriandensity,henceaRobinbound- ary condition, – or a prescribed pedestrian density, giving a Dirichlet boundary condition. Let T denote the common target/goal of the crowd, which is a subset of the boundary i.e. T ⊂∂Ω. We set the pedestrian density to m=0 at the target, hence individuals immediately leave the domain. On the rest of the boundary we impose homogeneous Neumann boundary conditions, i.e. individuals can not penetrate the walls. For the Eikonal equation we set u = 0 at the target and a suitable Dirichlet boundary condition on the rest of the boundary. The above conditions can be summarized as follows: m(x,0)=m (t), on Ω×{0}, m(x,t)=0,0 on T ×(0,T), u(x,t)=0, on T ×(0,T), (3) (uε(∇x,mt)+=fg2((xm))∇um)(x,t)·nˆ(x)=0, oonn ∂∂ΩΩ\\TT××((00,,TT)),, wherenˆ denotestheouternormalvectortotheboundary,whichisassumedto be smooth. Since the theoretical analysis of (2)-(3) has been done in [26] in 1DwithhomogeneousDirichletboundaryconditions,ratherthantacklingthe theoretical analysis of (2)-(3), in this work we focus on the efficient numerical discretization as we detail below. Semi-Lagrangian(SL)schemeshavebeensuccessfullyusedtodiscretizeHamil- ton-Jacobi-Bellman(HJB)equations,see[27]andthereferencestherein.They are based on approximating the characteristics of the problem. A SL scheme has been presented in [18] to deal with linear FP equations. In this work, we use a SL scheme to numerically solve the stationary HJB equation in (2). We propose an extension of the scheme in [18] in order to deal with nonlinear FP equations posed on a bounded domain. One of the main advantages of SL schemes is that they are explicit and allow large time steps. This is of special relevance since we are interested in the behavior of the solutions for arbitrary values of the horizon T which can be large (for example, if we are interested in the evacuation time). Moreover the 4 ElisabettaCarlinietal. SL discretization allows us to run stable simulations for small regularization parameters, closer in the spirit to the original hyperbolic system proposed by Hughes. This paper is structured as follows: in Section 2 we introduce the necessary preliminaries, including the trajectiorial interpretation of both equations, to presentandstudytheSLdiscretizationsinSection3.InSection4weillustrate the influence of the diffusivity on different performance parameters, such as the evacuation time of the crowd or the formation of congestions. 2 Preliminaries In this section we recall the stochastic optimal control interpretation of the HJB as well as the probabilistic interpretation of solutions of FP equations and introduce some notations used throughout this paper. Let Ω ⊂Rd denote a bounded domain with a smooth boundary ∂Ω. Assume that the common target of the crowd is on part of the boundary ∂Ω, hence T ⊂∂Ω. Let us consider a probability space (Ω,F,F,P) (where F is a σ-algebra, P is a probability measure on F, F := (F ) is a filtration in (Ω,F), i.e. s s≥0 F ⊆ F for all s ≥ 0 and F ⊆ F for all 0 ≤ s ≤ s ). We assume that F s s1 s2 1 2 satisfies the usual hypothesis (see e.g. [42]). We denote by E the expectation operator in this probability space. TrajectorialinterpretationoftheHJBequation. Itiswellknownthattheclas- sical solution u of the first equation of (2) can be represented as the value function of an associated stochastic optimal control problem, which we recall now. Given a process α adapted to F (i.e. α(s) is F -measurable for all s) and s satisfyingthatE(cid:0)(cid:82)s|α(r)|2dr(cid:1)<∞foralls≥0(wesaythatαisadmissible), 0 and x∈Ω, we define √ (cid:82)s y (s)=x+ α(r)dr+ 2εW(s) for all s>0, x,α 0 (4) and τ :=inf{s>0; y (s)∈∂Ω}, x,α x,α where W is a d-dimensional Brownian motion adapted to F. Note that the time τ , which corresponds to the first time the trajectory y leaves the x,α x,α domain Ω, is a stopping time for the filtration F (i.e. {τ ≤ s} ∈ F for all x,α s s). Let us fix t ∈ [0,T]. Classical results in stochastic control theory (see e.g. [28]) imply that, if m(·,t) is regular enough, then u(x,t)=inf (cid:26)E(cid:18)(cid:90) τx,α(cid:2)1|α(s)|2+(2f2(m(y (s),t))+δ)−1(cid:3)ds α 2 x,α 0 (5) (cid:17)(cid:27) +g(y (τ )) , x,α x,α ASLschemefortheHughesmodelforpedestrianflow 5 and the optimal feedback law is given by α∗(x,t) = −∇u(x,t) for all s ≥ 0. The function g is supposed to be strictly positive and taking sufficiently large values on ∂Ω\T to incite that agents move towards the target T. Thedependenceonthetimevariablet,seenasaparameterin(5),meritssome additional comments. Indeed, the dependence of u on t is due exclusively to the local density m(x,t) on the right-hand-side of the HJB equation. This implies that the trajectories y (·) in (4) are fictive in the sense that in the x,α optimization process agents take into account the current pedestrian distri- bution m(·,t) only. This a fundamental difference to mean field game models (see [35,34,12]) and mean field type control problems (see [7,19]), in which individuals anticipate the future dynamics of the crowd. Trajectorial interpretation of the nonlinear FP equation. The trajectorial in- terpretation of the nonlinear FP equation is provided through stochastic dif- ferentialequationsofMcKean-Vlasovtype(ormeanfieldtype),see[37,38,39, 46].Moreprecisely,letusconsidertheStochasticDifferentialEquation(SDE) √ dX(t)=b(X(t),µ(X(t),t),t)dt+ 2εdW(t), for all t≥0, (6) X(0)=X0, where b : Rd ×R×R → Rd is a regular vector-valued function, X0 is a + randomvectorinRd,independentoftheBrownianmotionW(·),withdensity m ,andµ(·,t)isthedensityofX(t).Itcanbeshown(see[33])that(6)admits 0 a unique solution and that µ is the unique classical solution of the nonlinear FP equation ∂µ−ε∆µ+div(b(x,µ,t)µ)=0 in Rd×[0,∞[, (7) µ(·,0)=m (·) in Rd. 0 Therefore, if we set b(x,m,t):=−∇u(x,t)f2(m(x,t)) (8) and working on Rd instead of Ω, equation (6) provides a formal probabilistic interpretation of the second equation in (2) with m(·,t) being the density of X(t). Let us point out that the interpretation is a priori only heuristic since u depends implicitly on m. Therefore the definition of b in (8) does not actually fit the framework of [33], where the dependence on the density is explicit. The probabilistic interpretation sketched above is the basis of our SL scheme to solve (2), presented in the next section. To include boundary conditions in the FP equation in (2) we reflect the discrete trajectories at ∂Ω \T and truncate them at T, see [11,29]. Finally, note that in contrast to Mean Field Games, the model considered in this work does not impose dual boundary conditions for the HJB and the FP equation. 6 ElisabettaCarlinietal. 3 The numerical scheme InthissectionweproposeaSLschemetoapproximatethesolutionof (2).The crucialpointisthediscretizationofthenonlinearFPequation,whichisbased on the fact that its solution is a measurable selection of the time-marginal densities of the diffusion defined by (6) (see [33]). We will first propose a SL scheme for a general nonlinear FP equation with smooth coefficients and a given velocity field depending explicitely on the density of the underlying stochasticprocess.Wewillprovethatourschemeisconsistentinanappropri- ate sense. The main feature of the scheme, which can be seen as an extension to the nonlinear case of the scheme proposed in [18,16], is that it is explicit and, at the same time, allows large time steps. This is not the case for e.g. explicit finite-difference schemes where the consistency property is achieved under the classical parabolic CFL condition. Inthecaseofsystem(2)thevelocityfieldinthenonlinearFPequationdepends implicitly on the density m through the solution u of the HJB. Therefore, in order to find an approximation of the velocity field we must solve the station- ary HJB equation at each time step. This is done in Section 3.2, where an adaptation of the fully-discrete scheme proposed in [15], taking into account theDirichletboundaryconditionispresented.Finally,inSection3.3wemerge both schemes to provide the fully-discrete scheme for (2). Letusbeginbyintroducingsomestandardnotation.Forsimplicity,wesuppose that Ω = (0,L)d. Even if this set Ω (and also the domains considered in the numerical simulations) has not a smooth boundary, we prefer to work on a square domain in order to simplify the scheme. Given a time step ∆t > 0 and a space discretization parameter ∆x > 0, let M ∈ N and N ∈ N be such that M∆x = L and N∆t = T. Let us set (x ,t ) := (i∆x,k∆t), where i k i ∈ {0,...,M}d and k = 0,...,N. For a given A ⊆ Ω we set G (A) := {i ∈ ∆x {0,...,M}d : x ∈ A} and call B(G (A)) and B(G (A)) the spaces i ∆x ∆x,∆t of grid functions defined on {x : i ∈ G (A)} and {(x ,t ),i ∈ G (A),k = i ∆x i k ∆x 0,...,N} respectively. GivenastandarduniformtriangulationofΩwithverticesbelongingtoG (Ω), ∆x we denote by {β ; i ∈ G (Ω)} the set of P -basis functions associated to i ∆x 1 this triangulation. We recall that β are continuous functions, affine on each i simplexandβ (x )=δ forallj ∈G (Ω)(whereδ denotestheKronecker i j ij ∆x i,j symbol). Moreover, the functions β have compact support and satisfy that i (cid:80) 0≤β ≤1 and β (x)=1 for all x∈Ω. We consider the following i i∈G∆x(Ω) i linear interpolation operator on Ω (cid:88) I[u](·):= u(x )β (·) for u∈B(G (Ω)). (9) i i ∆x i∈G∆x(Ω) ASLschemefortheHughesmodelforpedestrianflow 7 3.1 A Semi-Lagrangian scheme for a nonlinear Fokker-Planck equation In this section we propose a SL scheme to numerically solve the following nonlinear FP equation (cid:40)∂ m−ε∆m+div(mb(x,m,t))=0 in Rd×(0,T), t (10) m(·,0)=m (·) in Rd, 0 where b : Ω ×R×[0,T] → Rd is a given smooth vector field, depending on m. By an abuse of notation we denote by m the smooth initial datum, now 0 defined on Rd with compact support. In order to formally derive the scheme, we multiply the first equation in (10) by a smooth test function φ with compact support and integrate by parts to get: (cid:90) (cid:90) φ(x)m(x,t )dx= φ(x)m(x,t )dx (11) k+1 k Rd Rd (cid:90) tk+1(cid:90) + [b(x,m(x,t),t)·∇φ(x)+ε∆φ(x)]m(x,t)dxdt. tk Rd We first approximate (11) as (cid:90) φ(x)m(x,t )dx= k+1 Rd (cid:90) [φ(x)+∆tb(x,m(x,t ),t )·∇φ(x)+∆tε∆φ(x)]m(x,t )dx. k k k Rd Using a Taylor expansion we obtain (cid:90) φ(x)m(x,t )dx= k+1 Rd 1 (cid:88)d (cid:90) √ [φ(x+∆tb(x,m(x,t ),t )+ 2dε∆te )]m(x,t )dx+ 2d k k (cid:96) k Rd (cid:96)=1 1 (cid:88)d (cid:90) √ [φ(x+∆tb(x,m(x,t ),t )− 2dε∆te )]m(x,t )dx, 2d k k (cid:96) k Rd (cid:96)=1 where e denotes the (cid:96)-th canonical vector in Rd. (cid:96) We define E =[x1− 1∆x,x1+ 1∆x]×...×[xd− 1∆x,xd+ 1∆x], i i 2 i 2 i 2 i 2 (12) m := 1 (cid:82) m(x,t )dx. i,k (∆x)d Ei k 8 ElisabettaCarlinietal. Approximatingtheintegralsoftheform(cid:82) c(x)m(x,t )dxbysums(∆x)dc(x )m , Ej k(cid:48) j j,k(cid:48) where c is a smooth function, j ∈Zd and k(cid:48) =0,...,N, we get (cid:88) φ(x )m = (13) j j,k+1 j∈Zd d d 1 (cid:88) (cid:88) 1 (cid:88) (cid:88) φ(Φ(cid:96),+[m(x ,t )])m + φ(Φ(cid:96),−[m(x ,t )])m , 2d j,k j k j,k 2d j,k j k j,k (cid:96)=1j∈Zd (cid:96)=1j∈Zd where, for µ∈R, j ∈Zd, k =0,...,N −1 and (cid:96)=1,...,d, we have defined √ Φ(cid:96),±[µ]:=x +∆tb(x ,µ,t )± 2dε∆te . (14) j,k j j k (cid:96) Given i∈Zd setting φ=β in (13), we have i d 1 (cid:88) (cid:88)(cid:16) (cid:17) m = β (Φ(cid:96),+[m(x ,t )])+β (Φ(cid:96),−[m(x ,t )]) m . (15) i,k+1 2d i j,k j k i j,k j k j,k j∈Zd(cid:96)=1 Finally, since m (cid:39) m(x ,t ), setting m = (m ) , (15) gives the fol- i,k i k k i,k i∈Zd lowing explicit scheme for m : i,k m =G(m ,i,k) ∀k =0,...,N −1, i∈Zd, i,k+1 k (16) mi,0 = (cid:82)Ei(m∆0x()xd)dx ∀i∈Zd, in which the nonlinear operator G is defined by d 1 (cid:88) (cid:88)(cid:16) (cid:16) (cid:17) (cid:16) (cid:17)(cid:17) G(w,i,k):= β Φ(cid:96),+[w ] +β Φ(cid:96),−[w ] w , (17) 2d i j,k j i j,k j j j∈Zd(cid:96)=1 for every w ∈B(Zd). Because of the explicit in time discretization the scheme is well-defined. Given the solution m of (16), we associate the function i,k m :Rd×[0,T]→R defined as: ∆x,∆t m (x,t):=m if x∈E and t∈[t ,t [, i∈Zd, k =0,...,N. (18) ∆x,∆t i,k i k k+1 Note that the scheme is conservative by definition, i.e. (cid:90) (cid:90) (cid:88) m (x,t )dx=(∆x)d m = m (x)dx for all k =1,...,N. ∆x,∆t k i,k 0 Rd i∈Zd Rd We extend (17) to B(Zd)×Rd×[0,T] by defining G (v,x,t):=G(v,i,k) if x∈E ∆x,∆t i and t∈[t ,t [, i∈Zd, k =0,...,N −1. k k+1 Followingsimilarcomputationsasinthederivationofthescheme,wecanprove that (16) is consistent. The consistency result in the following Proposition is called weak in order to underline consistency to the weak formulation of (10). ASLschemefortheHughesmodelforpedestrianflow 9 Proposition 1 (Weak consistency) Assume that m : Rd ×[0,T] → R + satisfies: (cid:82) • m(x,t)dx is uniformly bounded in [0,T]. Rd • For all t∈[0,T], m(·,t)∈C2(Rd) and for all x∈Rd, m(x,·) is Lipschitz with a Lipschitz constant independent of x. Set m and m as in (12) and (18). Then, assuming that b is Lipschitz, i,k ∆x,∆t for every φ∈C∞(cid:0)Rd(cid:1) and k =0,...,N we obtain 0 (cid:90) (cid:90) φ(x)m (x,t )dx= φ(x)m(x,t )dx+O(∆x), (19) ∆x,∆t k k Rd Rd and for k =0,...,N −1 (cid:90) φ(x)G (m ,x,t )dx ∆x,∆t k k Rd (cid:90) (cid:90) tk+1(cid:90) = φ(x)m(x,t )dx+ b(x,m(t,x),t)·∇φ(x)m(x,t)dxdt k Rd tk Rd (cid:90) tk+1(cid:90) + ε∆φ(x)m(x,t)dxdt+O(∆x+(∆t)2). tk Rd (20) Inparticular,ifmisdifferentiablew.r.t.tothetimevariableandif(∆x ,∆t ) n n is a sequence of space and time steps such that (∆x ,∆t )→0 and ∆x /∆t →0 n n n n as n→∞, then 1 (cid:90) lim φ(x)[m (x,t )−G (m ,x,t )]dx n→∞∆tn Rd ∆xn,∆tn kn+1 ∆xn,∆tn kn kn (21) (cid:90) = φ(x)[∂ m(x,t)−ε∆m(x,t)+div(b(x,m(x,t),t)m(x,t))]dx, t Rd for kn such that t →t. kn Proof Let C =supp(φ), which is a compact set. By definition (cid:90) (cid:90) (cid:88) φ(x)m (x,t )dx= m φ(x)dx ∆x,∆t k i,k Rd i∈G∆x(C) Ei (cid:90) (cid:88) = m(x ,t )φ(x)dx+O((∆x)2) i k i∈G∆x(C) Ei (cid:90) (cid:88) = m(x,t )φ(x)dx+O(∆x+(∆x)2) k i∈G∆x(C) Ei (cid:90) = m(x,t )φ(x)dx+O(∆x) k Rd 10 ElisabettaCarlinietal. where, we have used that 1 (cid:90) m = m(x,t )dx=m(x ,t )+O(∆x2), i,k (∆x)d k i k Ei whichholdstruebyaTaylorexpansion,sincemisregular.Ontheotherhand, (cid:90) (cid:90) (cid:88) φ(x)G (m ,x,t )dx= G(m ,i,k) φ(x)dx. ∆x,∆t k k k Rd i∈G∆x(C) Ei Now, (cid:90) (cid:88) 1 (cid:88)d (cid:88) (cid:90) G(m ,i,k) φ(x)dx= β (Φ(cid:96),s[m ])m φ(x)dx, k 2d i j,k j,k j,k Ei j∈Zd l=1s∈{+,−} Ei (cid:88) 1 (cid:88)d (cid:88) (cid:90) = β (Φ(cid:96),s[m ]) m(x,t )dxφ(x )+O(∆x), 2d i j,k j,k k i j∈Zd l=1s∈{+,−} Ej where we have used that 1 (cid:90) φ(x)dx=φ(x )+O(∆x). (∆x)d i Ei Therefore, since (cid:88) β (Φ(cid:96),s[m ])φ(x )=I[φ](Φ(cid:96),s[m ])=φ(Φ(cid:96),s[m ])+O((∆x)2), i j,k j,k i j,k j,k j,k j,k i∈G∆x(C) interchanging the sums w.r.t. i and j, we get (cid:90) φ(x)G (m ,x,t )dx ∆x,∆t k k Rd (cid:88) (cid:90) 1 (cid:88)d (cid:88) = m(x,t )dx φ(Φ(cid:96),s[m ])+O((∆x)2+∆x). (22) k 2d j,k j,k j∈Zd Ej l=1s∈{+,−} Note that for x∈E j d 1 (cid:88) (cid:88) φ(Φ(cid:96),s[m ]) 2d j,k j,k l=1s∈{+,−} =φ(x )+∆tb(x ,m ,t )·∇φ(x )+ε∆t∆φ(x )+O((∆t)2), j j j,k k j j (cid:90) tk+1 =φ(x)+ [b(x,m(x,t),t)·∇φ(x)+ε∆φ(x)]dt tk +O(∆x+∆x∆t+(∆t)2). The equality in (20) follows easily from the relation above and (22). Finally, relation(21)followsdirectlyfrom(19)andanintegrationbypartsinthespace variable. (cid:4)