Numerical investigation of agent controlled pedestrian dynamics using a structure preserving finite volume scheme

We provide a numerical realisation of an optimal control problem for pedestrian motion with agents that was analysed in Herzog, Pietschmann, Winkler:"Optimal Control of Hughes' Model for Pedestrian Flow via Local Attraction.", arXiv 2011.03580, 2020. The model consists of a regularized variant of Hughes' model for pedestrian dynamics coupled to ordinary differential equations that describe the motion of agents which are able to influence the crowd via attractive forces. We devise a finite volume scheme that preserves the box constraints that are inherent in the model and discuss some of its properties. We apply our scheme to an objective functional tailored to the case of an evacuation scenario. Finally, numerical simulations for several practically relevant geometries are performed.


I
With more and more people living in highly populated areas, the modelling, simulation and control of (large) pedestrian crowds is an important eld of research.In this work, we study the optimal control problem for a regularized version of Hughes' model for pedestrian motion, Hughes, .In our approach, the (continuous) crowd can be controlled by a xed small number of agents that can attract people in their vicinity.In terms of the model, this corresponds to an additional potential term centred at the agents positions.In a previous work, Herzog, Pietschmann, Winkler, , we already studied the well-posedness and optimality conditions of this problem while here, we focus on a numerical http://www.tu-chemnitz.de/mathematik/partdg /peop e/stoetzner/, ORCID ---, max.winkler@mathematik.tu-chemnitz.de,https://www.tu-chemnitz.de/mathematik/partdg /peop e/wink er/, ORCID ---).
implementation of the control problem and extensive numerical examples.In particular, we provide and analyse a nite volume scheme that preserves the box constraints inherent in our problem.
To introduce the model, we x Ω ⊂ R to be a bounded domain with  -boundary Ω.Furthermore  > is an arbitrary time horizon and   := ( , ) × Ω as well as Σ  = ( , ) × Ω denote the space-time cylinder and its lateral boundary, respectively.The boundary is decomposed into two parts: Ω  representing the exits and Ω  the part where the domain is constrained by walls.For theoretical purposes (regularity of solutions) we assume Ω  ∩ Ω  = ∅ meaning that both boundary parts are separated from each other, see Figure . .In a similar way we de ne Σ  = ( , ) × Ω  and Σ  = ( , ) × Ω  .
The unknown variables in our system of equations are the density of the crowd  :   → R + , a potential specifying the current time to escape  :   → R. In addition, there are  agents which may in uence the motion of the crowd via attractive forces.Their positions are denoted by   : ( , ) → R ,  = , . . ., .In addition, each agent is able to regulate the strength by which it acts on the crowd.This is encoded in the intensities   : ( , ) → R + ,  = , . . ., .Both the agent trajectories and interaction strength are summarized in a vector  = ( , . . .,   ) and  = ( , . . .,   ) , respectively.
The mathematical equations describing the movement of a pedestrian crowd in uenced by agents then read as follows.For given agent movement directions  = ( , . . .,   ) with   ∈  ∞ ( , ; R ) the unknowns , ,  are related to each other by means of Let us brie y discuss the meaning of the respective terms: Equation ( .a) states that pedestrians are transported according to the velocity eld , see ( . ) below, while also performing (little) random motion encoded by the Laplacian of .The second equation ( .b) is a modi ed and regularized Eikonal equation whose solution is the distance to the closest exit, mitigating areas of high density via the term on the right-hand side.Here, the additional di usion accounts for the fact that pedestrians do not know their environment exactly.Then, ( .c) governs the motion of the agents, whose speed is also in uenced by the surrounding pedestrian density.The function  : [ , ] → [ , ] is a density-velocity rule, chosen in such a way that  () determines the maximum velocity an individual can move if the density in its current position is .We choose  to be monotonically decreasing meaning that higher densities lead to slower movements.The velocity eld  will re ect the fact that pedestrians are, on the one hand, trying the minimize their exit time which amounts to a drift term in the direction of ∇ and on the other hand, they are attracted by the agents which is realized by additional attractive potentials whose center depends on the agents' positions .This results in a velocity which is the sum of two terms.Furthermore, to account for the e ect that the velocity will deteriorate in regions of high density, it will be modi ed by an additional multiplicative factor  ().As in the equations for the motion of the agents,  is monotonically decreasing and becomes zero at a given maximal density.
The boundary conditions ( . ) allow for an out ow with velocity  on parts of the boundary (Σ  ) while no-ux conditions on the remaining parts are to be interpreted as walls (Σ  ).A detailed description of the involved non-linearities will be given in the next section, but we also refer to Herzog, Pietschmann, Winkler, for more details on the model and the regularizing terms.
Analytical properties of the unregularized Hughes' model introduced in Hughes, (i.e. =  =  = in ( .)), without control, are di cult because of the low regularity of ∇ that appears on a set depending on the solution  of the rst equation, but see Amadori, Goatin, Rosini, ; Amadori, Di Francesco, ; El-Khatib, Goatin, Rosini, .Thus, regularized variants have been considered, see Di Francesco, Markowich, et al., for an instance where  = but  ,  ≠ .In fact, the result there is obtained as a vanishing viscosity limit  → .There is also a number of extensions and variants of the model, aiming to understand additional properties, make it more realistic, or consider di erent settings like graphs, see Burger, Di Francesco, et al., ; Carrillo, Martin, Wolfram, ; Carlini et al., ; Di Francesco, Fagioli, et al., ; Colombo, Gokieli, Rosini, .Control of systems by means of a small number of agents has received lots of interest recently, both on a discrete level (i.e. one considers a large system of ODEs for the motion of individuals coupled to a small number of equations for the agents), see Caponigro et al., ; Himakalasa, Wongkaew, , but also for coupled PDE-ODE systems, Albi, Bongini, et al., ; Albi, Fornasier, Kalise, .Let us emphasize that whenever PDEs are coupled to ODEs in such a fashion that the solution of the PDE needs to evaluated at the solution of the ODE (as in ( .c)), regularity is needed.While in our case, this is obtained by the additional di usion in ( .a), when hyperbolic models for the transport of pedestrian are considered, an additional regularization in the ODE is needed, see e.g.Borsche, Colombo, et al., ; Borsche, Klar, et al., ; Borsche, Meurer, .Finally, let us mention that the ODE-ODE and PDE-ODE perspective are closely related by means of so-called mean eld limits when the number of agents tends to in nity, see Burger, Pinnau, et al., ; Pinnau, Totzeck, and also the recent overview on control of crowds, Banda, Herty, Trimborn, .For the numerical discretization of ( .a), we employ, as we think of the parameter  being small, a nite volume scheme for the spatial discretization, which may also be interpreted as a discontinuous Galerkin scheme.In combination with the Lax-Friedrichs numerical uxes the scheme is stable and preserves the bounds ≤  ≤ inherent in our model.Such structure-preserving discretizations of PDEs gained much attention, e.g., in the context of chemotaxis problems Epshteyn, Kurganov, ; Li, Shu, Yang, ; Strehl et al., ; Guo, Li, Yang, ; Ibrahim, Saad, ; Filbet, .The previously mentioned articles di er also in the choice of the time-stepping scheme.For the treatment of ( . ) we will use an implicit-explicit nite di erence scheme, whereas the di usion-related terms are established implicitly and the convection-related terms explicitly.The equation ( .b) is discretized with standard linear nite elements and for ( .c) we employ a backward Euler scheme.
The paper is organized as follows: In Section we give a precise de nition of our problem and recall the analytical results from Herzog, Pietschmann, Winkler, .In Section we provide a numerical discretization scheme in space and time and analyse some of its properties, in particular, we show that it preserves physical bounds of the density of pedestrians.Corresponding optimization algorithms are discussed in Section .and Section nally provides the results of our numerical experiments.

T
Let us motivate the remaining quantities arising in the system of equations.The function  : [ ,  max ] → R is a density-velocity relation and is assumed to be  ,∞ (R) ∩   (R) with  ( ) = and  ( max ) = with  max denoting the maximal density.A usual choice is Obviously, a higher density leads to a lower velocity.Throughout the article we set  max = .The movement direction of the crowd described by the function − (, , , ) is modelled as follows.The primary interest of the crowd is to move either towards the closest emergency exit, this is the direction −∇.This is mitigated by the attraction of close by agents which is the direction −∇  , where is an agent potential de ned by with certain parameters ,   > , realizing a repulsion in the near and an attraction in the far range of the agents.This is useful to avoid a high density very close to the respective agent.We refer to Carillo, Huang, Martin, for a more detailed discussion on potentials in the context of ocking problems.
These considerations yield a velocity eld de ned as follows where ℎ is a smoothed projection into the unit ball in R and the factor  () again links the allowed movement speed to the current density, this is, | (, , , )| ≤   () in   .
We brie y introduce the function spaces used in the sequel, see also Denk, Hieber, Prüss, .For a domain Ω ⊂ R we denote by , (Ω),  ∈ N ,  ∈ [ , ∞], the usual Sobolev spaces and by − /, (Γ) for  ≥ the corresponding trace spaces which may be equipped with the Sobolev-Slobodetskij norm.
as well as For the application we have in mind the following spaces are of interest which are equipped with the natural norms .
Spaces with non-integral  and  are de ned as interpolation spaces.
In a previous work, Herzog, Pietschmann, Winkler, , a global (in time) well-posedness and regularity result for ( .)-( .) was established.Furthermore, optimality conditions for related optimal control problems where this system occurs as a constraint were derived.For convenience of the reader we brie y summarize the most important results needed in the present article.
The previous result allows to de ne an operator, the so-called control-to-state operator, with control and state spaces . Furthermore we de ne the set of admissible controls ) and all  = , . . .,  }.
--cbna page of The optimal control problem we study in this article reads Minimize J (; ) ( .c) The objective functional J aims at a fast evacuation of the crowd.By the factor    higher densities at a later time are penalized more.We observe the density in a subregion   =  × Ω where Ω ⊂ Ω is a subregion which the pedestrians must leave.We use the temporal  -norm of the agent movement directions and the intensities as a regularization to guarantee the smoothness required in Theorem . .The regularization parameters  ,  > are arbitrary but positive.
The fourth term in the objective is a barrier used to avoid that the agents walk through walls.The barrier function  ∈   (Ω) is the weak solution of the singularly perturbed problem Here we choose  > to be xed but small.
The control constraint (, ) ∈ Q ad guarantees that the agents do not move faster than the density in their current position allows and that the intensity is bounded by reasonable values.
We have the following well-posedness result and necessary optimality condition.
Furthermore, each local minimizer (, ) ∈ Y × Q ad ,  = (, , ),  = (, ), of ( . ) ful lls for all directions in the tangential cone at (, ), namely  = (, ) ∈ T Q ad (, ), for  = , . . ., , together with the boundary conditions ( . ) and homogeneous initial conditions The proof of the theorem above is very close to those of Theorem .and Theorem .Herzog, Pietschmann, Winkler, .The main di erence is that the model in Herzog, Pietschmann, Winkler, only allows to control the velocity  of the agents, yet not their strength .As for the existence proof, this does not impose any additional di culty due to the uniform  ∞ -boundedness of  as a consequence of the embedding  ↩→  ∞ in one spatial dimension.For the di erentiability result, one has to add the derivatives with respect to , yielding an additional term in ( .a) that, however, can be estimated similarly to the remaining terms.

D
In this section we introduce the numerical scheme used to compute approximate solutions of the forward system ( .)-( .).To this end, we introduce a semi-implicit time-stepping scheme and use a nite volume discretization for the density function  and continuous Lagrange nite elements for the potential function.

. S
For the spatial discretization of the system ( .)-( .) we de ne a family of geometrically conforming triangular meshes {T ℎ } ℎ> .For each  ∈ T ℎ we denote by ℎ  = diam( ) the element diameter and by   the diameter of the largest inscribed ball in  .The mesh parameter is then ℎ = max  ∈ T ℎ ℎ  .The mesh family is assumed to be shape regular meaning that there is a constant  > such that ℎ  /  ≤  for all  ∈ T ℎ and all ℎ > .By F i ℎ we denote the set of interior element edges, by F bd ℎ the boundary edges and write F ℎ := F i ℎ ∪ F bd ℎ .Furthermore, to each edge  ∈ F ℎ we associate the normal vector   which is pointing outward in case of a boundary edge and has a xed orientation in case of an interior edge.
We propose a nite volume scheme for the transport equation.As we use the nite element package FE CS for our implementation, we use a notation which is rather usual for discontinuous Galerkin discretizations, see Di Pietro, Ern, for an overview.The nite-dimensional function spaces are de ned by where P  ( ) denotes the space of polynomials on  of degree not larger than  ∈ N .For a function  : Ω → R, we de ne interface averages and jumps in the following way where  , ∈ T ℎ are chosen in such a way that In order to discretize the system ( . ) we use discontinuous approximations for  and continuous ones for .The unknowns in our semi-discrete scheme are The approximate transport direction is then given by with The semi-discretization of ( .a) then reads and Here, proj  ℎ :  (Ω) →  ℎ is some projection operator, (•, •) Ω is the standard  (Ω)-inner product and the bilinear forms are de ned by The parameter   is de ned by where   is the intersection of the orthogonal edge bisectors of  ∈ T ℎ .The term    ℎ  ℎ with  ℎ =   for some  ∈ T ℎ approximates the di usive ux ∇ ℎ •   over the edge  ⊂  .The bilinear form  establishes the convective ux   •    .As numerical ux function (•) * , we choose the Lax-Friedrichs ux, see Rider, Lowrie, , de ned by The stabilization parameter  ∈ R is speci ed later.For the closely related chemotaxis model such an approach has been sucessfully applied in Li, Shu, Yang, ; Guo, Li, Yang, .Of course, also other ux functions are possible, e.g., the central upwind ux, Epshteyn, Kurganov, .
The Eikonal equation ( .b) is discretized in space using standard linear Lagrange elements which yields In our numerical experiments we used the Newton solver integrated in FE CS.The Jacobian is established by automatic di erentiation.
The ordinary di erential equations for the agent trajectories ( .c) depend on a point evaluation  (,   ()) of a function which is discontinuous in the discrete setting.In particular, this term would not be di erentiable with respect to   ().As a remedy, we use instead of a point evaluation, see also Borsche, Colombo, et  .T For the temproal discretization we cover the time interval [ , ] by an equidistant grid   {  }  = with grid points   :=   and grid size  :=  / .The spatial and temporal discretization parameters are summarized in  = (ℎ, ).Moreover, we de ne the space of grid functions with  an arbitrary linear space.If  is again a function space containing functions  : Ω → R we write  (  ) =  (  , •).The functions  ℎ ,  ℎ ,   ,   and   arising in the semidiscrete equations ( .), ( . ) and ( . ) are approximated by grid functions for  = , . . ., .For brevity we write for  = , . . ., and for the transport vector we use We replace the temporal derivative by a di erence quotient and use a semi-implicit time-stepping scheme, more precisely, the di usion-related terms are evaluated implicitly and the convection-related terms explicitly.This yields the fully-discrete system Note that the system of equations ( .)completely decouples and we can compute each variable after the other, in the following order

( . )
. Q In this section we study some basic properties for the solutions of ( .).In particular, it is of interest whether the physical bounds observed for the solution of the continuous system ( .)-( .) are transferred to the discrete setting.
The basis functions of the nite element space  ℎ are denoted by {  }  ∈ T ℎ de ned by   |  ≡   , for all  , ∈ T ℎ .Note that by a slight abuse of notation we use the elements of T ℎ as indices here.
Theorem . .The numerical scheme ( .a) is mass conserving in the following sense.Assuming that  = holds, i. e., there are no-ux boundary conditions present at all boundary parts Ω D and Ω W , the solution   ful lls Proof.The assertion is trivially ful lled for  = as the initial condition is established by  ℎ = proj  ℎ ( ).In matrix-vector notation the assertion is equivalent to ì  ì  + = ì  ì   .This follows from ( . ) after using ì  ì In this expression, the stabilization terms (the ones multiplied with ) cancel out.Furthermore, after sorting terms in the sum by the edges  ∈ F i ℎ and denoting by  , , , the two triangles meeting in  , we obtain the terms The last step follows due to    , = −   , which implies ì   ì   = .
Let  ∈ N be xed and assume that   ℎ () ∈ [ , ] for almost all  ∈ Ω.We show  + ℎ ≥ by con rming that the right-hand side of ( . ) has non-negative entries only.Assuming that    ≥ ,  ∈ T ℎ , we show that the entries of ( −    ) are non-negative as well.The entries on the diagnoal have the form where the rst step follows from the assumption ( . ) implying |  ℎ | ≤ .To estimate further we take the geometric mesh quantities and the shape regularity ℎ  ≤    into account and arrive at where the last step is the CFL condition ( .).The previous two estimates con rm (  , −     , ) ≥ for all  ∈ T ℎ .The remaining entries in the matrix, namely (  , −     , ) with The non-negativity follows again from |  ℎ | ≤ .Combining the previous arguments provides the lower bound  + ℎ ≥ .
To show the upper bound we rearrange the equation system ( . ) in the form We may rewrite the transport term using In ( . ) we reformulate the expression involving   by means of with   ∈ T ℎ \ { },   ∩  =  .With the same arguments like above one can show that the entries of  +    are non-negative and together with − ì   ≥ ,  ì ≥ and the M-matrix property of  +   we arrive at the desired bound − ì  + ≥ .
. T Next, we study a discrete version of the optimal control problem ( .).The control and state variables are grid functions in time and thus, we introduce the discrete  ( , )-inner product for functions   ,   : {  }  = →  , with  some Hilbert space, This induces the norm    ( , ; ),ℎ := (  ,   )  ( , ; ),ℎ .For the discrete control space we obtain The discrete state space is de ned by With these de nitions, the discrete optimal control problem related to ( . ) reads as where is  ℎ ∈  ℎ the nite element approximation of ( . ) with rst-order Lagrange elements.Furthermore,   is the solution operator of ( .).Note that, in order to maintain the di erentiability of the barrier term with respect to    , we use a regularization of the point evaluation of the nonsmooth function  ℎ , compare also ( .c).
We may write the control problem ( . ) in the more compact reduced form To deduce a necessary optimality condition we apply the Lagrange formalism.The Lagrange function coupling the discrete state equation ( . ) reads To shorten the notation we write   := ( , ,  , ,  , ).The adjoint equation system determining these variables for a given control and state is ,ℎ ∈  ℎ,D : for  = , . . ., .Note that this can be interpreted as a coupled system involving a parabolic PDE and an ODE that run backward in time.We use the automatic di erentiation feature in FE CS in our implementation.As the forward system completely decouples in each time step, so does the adjoint system and we can compute step by step: With the adjoint states at hand we can assemble the derivatives of the reduced objective ( . ) and end up with the following optimality condition for ( .): Theorem .(Necessary optimality condition).Let (  ,   ) ∈ Y  × Q ,ad be a local solution of ( .).We write   = ( , , . . .,  , ) ∈ U and   = ( , , . . .,  , ) ∈ C and get the following representation of the gradient of   : This allows an implementation of a projected gradient method which we discuss in the following section.
. O For a solution of the discretized optimal control problem .we propose a projected gradient algorithm.
In this procedure, for a given initial control  ( ) = ( ( )  A formula for the gradient of   has been derived in the previous section already, see ( .).The step length parameter  () > is obtained by an Amijo line search and must ful ll the su cient decrease condition with a decrease parameter  ∈ ( , ) which is usually small (e.g.− ).A reasonable stopping criterion for the projected gradient algorithm is It remains to discuss the realization of the projection operators and we propose a primal dual active set strategy that may also be considered as semismooth Newton method.Note that the operators Π ad are semismooth, see Christof, Wachsmuth, .The evaluation of the projection operator   ad : U  → U ,ad requires to solve the optimization problem ( .a).The unknowns (assuming  = and omitting the agent's index  for a while) are the coe cients of the functions U   ad (  ) =   ì  ∈ R ( + )× for some given U    ì  ∈ R ( + )× .We switch to a matrix-vector notation and de ne as well as the matrix  ∈ R ( + )×( + ) on the left-hand side of the linear system ( . ) inducing the discrete  ( , ), -norm.The Lagrangian for ( . ) reads . The Karush-Kuhn-Tucker system for ( . ) then reads where • is the component-wise multiplication of two vectors.We reformulate the complementarity condition by means of a nonsmooth equation and arrive at the following equivalent form of the KKT system ( . ) This nonlinear system can be solved iteratively by a semismooth Newton method.Given is an initial pair (  ( ) , ì  ( ) ). Successively, one computes the active and inactive set . .,  }, and performs the Newton update This procedure is repeated for  = , , . . .until some termination criterion, e. g.,  ( ì  , ì  , ì ) < tol, is ful lled.

N
This section is devoted to numerical experiments.To establish the discretized system ( . ) the nite element library FE CS was used, complemented by a P implementation of the projected gradient method from Equation ( . ) and the Armijo step size rule from ( .).The computational meshes were created by the mesh generator mshr integrated in FE CS.

. E
In a rst numerical test we solve the problem Equation ( . ) in the domain Ω depicted in Figure .with the following parameters.The initial density  is the sum of Gaussian bells, see also Figure .a.The subdomain Ω where densities are penalized is chosen to cover the region within the walls.Without any controlled agents, most of the people will squeeze through the smaller emergency exits in the south and north while the large exit in the east is rarely used.To improve the evacuation agents were introduced.The initial control (, ) was chosen in such a way that the agent moves straight to the right outside of the room having a constant the intensity.
. E In a second example we consider the domain illustrated in Figure . .The initial density is concentrated near the slit in the wall on the left-hand side.In an uncontrolled evacuation scenario the majority of the people would leave the domain through this slit causing a massive congestion leading to a very slow evacuation of the crowd.The model and algorithm parameters chosen in the current example are as follows: This example shows that the evacuation can be signi cantly improved by using two agents with optimized trajectory and intensity.Interesting is, that the intensity is non-zero only in the time interval  ∈ ( , ).The agents attract the people leading them su ciently far away from the slit in the west and then they stop in uencing the crowd.When being su ciently far away from the slit the people nd the way to the larger exits in the north and south on their own by using the movement direction determined by the potential .

. E
In a third example we consider a square-shaped room with exits in the south, east and north.The exits have di erent width.The model and algorithm parameters are chosen as follows: The initial density is concentrated near the small exit and without a control of the crowd motion most of the people are blocking each other while squeezing through this small exit.Two agents were added in this scenario with the aim attracting the people in such a way that more of them nd the other two exits in the north and east.The computed agent trajectories are quite short.It is interesting to observe that in the time interval  ∈ ( , ) the agents just go to an optimal position su ciently close to the crowd and attract them in the time interval  ∈ ( , ), leading some of the people to the center of the room.At this point the agents drive their intensity to zero meaning that they stop in uencing the crowd.However, when being su ciently far away from the critical exit the people nd the route to the less used exits on their own due to the movement rule determined by the potential .
:= (  ,   ,   ) =   (  ), ( .b) Figure .: Solution of the problem from Section . .The colored background represents the density ; the dots are the agent positions; the black curves are the agent trajectories.
Figure .: Solution of the problem from Section .at various time steps   and intensity of the agents.
Figure .: Solution of the problem from Section .at various time steps   and intensity of the agents.