Investigation of ZnO doping on LaFeO₃/Fe₂O₃ prepared from yarosite mineral extraction for ethanol gas sensor applications

1.   Introduction

The theory of mean field games, (MFGs for short), aims at studying deterministic or stochastic differential games (Nash equilibria) as the number of agents tends to infinity. It supposes that the rational agents are indistinguishable and individually have a negligible influence on the game, and that each individual strategy is influenced by some averages of quantities depending on the states (or the controls as in the present work) of the other agents. MFGs have been introduced in the pioneering works of J. M. Lasry and P. L. Lions ^[17,18,19]. Independently and at approximately the same time, the notion of mean field games arose in the engineering literature, see the works of M. Y. Huang, P. E. Caines and R. Malhamé ^[14,15].

The present work deals with numerical approximations of mean field games in which the agents interact through both their states and controls; it follows a more theoretical work by the second author, ^[16], which is devoted to the mathematical analysis of the related systems of nonlocal partial differential equations. There is not much literature on MFGs in which the agents also interact through their controls, see ^{[7,8,10,12,13,16]}. To stress the fact that the latter situation is considered, we will sometimes use the terminology mean field games of control and the acronym MFGC.

The MFGC that is considered in the present work is described by the following system of nonlocal partial differential equations:

Here, the state space $\Omega$ is a bounded domain of $\mathbb{R}^d$ with a piecewise smooth boundary, while the controls are vectors of $\mathbb{R}^d$. The time horizon is $T > 0$, and the parameter $\nu > 0$ is linked to the level of noise in the trajectories (the dynamics of a given agent is described by $dX_t = \alpha_t dt +\sqrt{2\nu} dB_t$ where $(B_t)$ is a standard $d$-dimensional Brownian motion and $\alpha_t\in \mathbb{R}^d$ is the control at time $t$). A special feature of the present model is that the third argument of the Hamiltonian $H$ is a probability measure $\mu(t)\in \mathcal{P} \left(\Omega\times \mathbb{R}^d \right)$, which stands for the joint probability of the states and optimal controls of the agents at time $t$. In (1), $u: [0, T]\times \Omega \to \mathbb{R}$ and $m : [0, T]\times \Omega \to \mathbb{R}_+$ respectively stand for the value function of a representative agent and the distribution of states. The first, respectively second line in (1) is the Hamilton-Jacobi-Bellman equation (HJB for short) which leads to the optimal control of a representative agent, respectively the Fokker-Planck-Kolmogorov equation (FP for short) which describes the transport-diffusion of the distribution of states by the optimal control law. The HJB equation is a backward w.r.t. time parabolic equation and is supplemented with a condition at $t = T$ which involves the terminal cost function $\phi: \Omega\to \mathbb{R}$, whereas the FP equation is forward w.r.t. time and is supplemented with an initial condition, which translates the fact that the initial distribution of states is known. The HJB equation also involves the so-called coupling function $f:\Omega\times \mathbb{R}_+\to \mathbb{R}$, or, in other words, $f(X_t, m(t, X(t)))$ is part of the running cost of a representative agent. We shall discuss later the boundary conditions on $(0, T)\times \partial \Omega$ associated with the first two equations in (1). The third equation in (1) gives the connection between $\mu(t)$ and $m(t)$, and shows in particular that $\mu(t)$ is supported by a $d$-dimensional geometrical object, namely the graph of the feedback law: $\Omega\to \mathbb{R}^d$, $x\mapsto - H_p (x, \nabla u (t, x), \mu(t))$.

1.1.   A brief discussion on the mathematical analysis of (1)

Recall that the Hamiltonian of the problem is $(x, p, \mu)\mapsto H(x, p, \mu)$, $(x, p, \mu)\in \Omega\times \mathbb{R}^d \times \mathcal{P} \left(\Omega\times \mathbb{R}^d \right)$.

From the viewpoint of mathematical analysis, a priori estimates for (1) are more difficult to obtain than in the case when the agents interact only via the distribution of states $m$. Indeed, in the latter case, if for example the costs $f$ and $\phi$ are uniformly bounded, then a priori estimates on $\left\| {u} \right\|$ stem from the maximum principle for second-order parabolic equations. By contrast, since the Hamiltonian in (1) depends non-locally on $\nabla_x u$, the maximum principle applied to the HJB equation only permits to bound $\left\| {u} \right\|$ by a quantity which depends (quadratically under standard assumptions on $H$) on $\left\| {\nabla_xu} \right\|$, and this information may be useless without additional arguments.

If the agents interact only through the distribution of states and if the Hamiltonian depends separately on $p$ and $m$, a natural assumption is that the latter is monotone with respect to $m$, see ^[17,18,19]; it implies existence and uniqueness of solutions, see in ^[19,20]. Such an assumption is quite sensible in many situations, since it models the aversion of the agents to highly crowded regions of the state space. It is possible to extend these arguments to MFGCs, see ^[10] for a probabilistic point of view and ^[8] for a PDE point of view, and the monotonicity assumption then means that the agents favor controls that are opposite to the main stream. In ^[16] and in the present work, we prefer to avoid such an assumption, because it is generally not satisfied, at least in models of crowd motion: indeed, in models of traffic or pedestrian flows, a generic agent would rather try to adjust her speed (control) to the average speed in a neighborhood of her position.

The third equation in (1) can be seen as a fixed point problem for $\mu$ given $u$ and $m$, which turns to be well-posed under the Lasry-Lions monotonicity assumption adapted to MFGC, provided that $u$ and $m$ are smooth enough. We shall replace this assumption by a new structural condition which has been introduced in ^[16], namely that $H_p$ depends linearly on the variable $\mu$ and is a contraction with respect to $\mu$ (using a suitable distance on probability measures). In the context of crowd motion, this structural condition is satisfied if the representative agent targets controls that are proportional to an average of the controls chosen by the other agents nearby, with a positive proportionality coefficient smaller than one. Were this coefficient equal to or larger than one, it would be easy to cook up examples in which there is no solution to (1) or even to the $N$-agent game, see Remark 4.3 below.

In ^[16], the focus is put on the existence of solutions rather than on their uniqueness. Indeed, without the monotonicity condition, uniqueness is unlikely in general if $T$ is not small. Consider for example a game in which the function $\phi$ has two perfectly symmetrical minima (the targets), $f = 0$ and where the initial distribution of states has the same symmetry. If $H$ depends on $|\nabla u(t, x)|$ only (no interaction through the controls), then a representative agent will simply travel to the minimum which is closest to her initial position. On the contrary, if the representative agent favors a control close to the average one, then there are at least two symmetrical solutions in which the whole distribution moves toward one of the two minima.

Going back to existence results, it is now frequent in the MFG literature to obtain energy estimates by testing the HJB equation by $m$, the FP equation by $u$, summing the resulting equations, integrating in the time and state variables then making suitable integrations by parts. If the value function is uniformly bounded from below (which is often the case even if there are interactions through the controls), this results in a relationship between the $L^{\infty} \left([0, T]; L^1 \left(\Omega \right) \right)$-norm of the positive part of $u$ and the $L_m^2 \left([0, T]\times\Omega; \mathbb{R}^d \right)$-norm of $\nabla_xu$; this observation can then be used to obtain additional a priori estimates, which may be combined with the ones obtained from the maximum principle and discussed above, and finally with Bernstein method. This strategy has been implemented in ^[16] and leads to the existence of a solution to (1) under suitable assumptions. By and large, existence of classical solutions was proved in ^[16] for periodic problems in any of the following cases:

● short time horizon

● monotonicity

● small enough parameters (in particular the contraction factor mentioned above, see the parameters $\lambda$ and $\theta$ in (7)–(10) below)

● weak dependency of the average drift on the state variable

Note that in ^[11], existence and uniqueness have been proved with probabilistic arguments in the case where the Hamiltonian depends separately on $p$ and $\mu$.

1.2.   A more detailed description of the considered class of MFGCs

For any $x\in \partial \Omega$, let $n(x)$ be the outward pointing unit normal vector to $\partial \Omega$ at $x$. The dynamics of a representative agent is given by

where $\left(B_t \right)_{t\in \left[0, T \right]}$ is a standard $d$-dimensional Brownian motion, $({\alpha}_t)_{t\in [0, T]}$ is the control, a stochastic process adapted to $(B_t)$ with values in $\mathbb{R}^d$, $L_t$ is the local time of the process $X_t$ on $\partial \Omega$. We assume that $X_{|t = 0}$ is a random variable in $\overline \Omega$, independent of $(B_t)$ for all $t > 0$, whose law is $\mathcal{L} \left(X_{|t = 0} \right) = m_0$.

In what follows, the part of the running cost which models the interactions via the controls will involve an average drift $V\in \mathbb{R} ^d$:

for $(t, x)\in[0, T]\times\Omega$, where $K$ is a nonnegative kernel, $Z(t, x) = \int_{\Omega}K(x, y)dm(t, y)$ is a normalization factor, $\mu(t, \cdot, \cdot)$ is the law of the joint distribution of the states and the controls, and $m(t, \cdot)$ is the law of the distribution of states.

Hereafter, we are going to focus in MFGCs in which the cost to be minimized by a representative agent is

where ${\lambda}, {\theta}$ are real numbers such that $-1 < {\lambda} < 1$ and $0 < {\theta}\leq 1$. This leads us to define the Lagrangian

which is convex with respect to $\alpha$, and its Legendre transform (with respect to $\alpha$):

With this definition of the running cost, the first three lines of (1) can be written as follows:

We can now specify the boundary conditions on $(0, T)\times \partial \Omega$. First, assuming that $u$ is smooth enough, the optimal control of a representative player is given by the feedback law:

The reflection condition at the boundary translates into the Neumann boundary conditions:

Now, from the definition of $m$, $\int_\Omega m(t, x)\psi(x) dx = \mathbb{E}\left(\int_\Omega m_0(x) \psi(X_{t, x}) dx\right)$ for all smooth enough test-function $\psi$ such that $\frac{\partial \psi}{\partial n}(x) = 0$ on $\partial \Omega$. Taking the time derivative of the latter equality, and using (13), we deduce that

Note that (8), (13) and (14) imply that the total mass $\int_\Omega m(t, x)dx$ is conserved.

We have studied or will sometimes consider other kinds of boundary conditions on some parts of the boundary, for example:

● periodic conditions, i.e., we set (7)–(12) in the flat torus $\Omega = \mathbb{R}^d/ \mathbb{Z}^d$

● At a part of the boundary standing for an exit, there is an exit cost. This yields a Dirichlet condition on $u$. The Dirichlet condition on $m$: $m = 0$ at exits, means that the agents stop taking part in the mean field game as soon as they reach the exit. Such boundary conditions will be simulated numerically in paragraph 4.2 below.

● If a part of the boundary stands for an entrance, then it is natural to impose a Dirichlet condition on $u$ to prevent the agents from crossing the entrance outward, and a flux-like condition on $m$ to specify the entering flux of agents, see also paragraph 4.2 below.

1.3.   Organization of the paper

Section 2 is devoted to the description of the finite difference scheme; it is based on a monotone upwind scheme for the HJB equation, and a scheme for the FP equation which is obtained by differentiating the discrete HJB equation and taking the adjoint. It is an extension of the finite difference schemes proposed and studied by the first author and I. Capuzzo-Dolcetta in ^[2,3] for simpler MFGs. Designing an efficient strategy for solving the system of non linear equations arising in the discrete version of MFGCs is a challenging task, in particular because

● the underlying system of PDEs is forward-backward and cannot be solved by simply marching in time

● there is no underlying variational structure

● Equation (9) is non local.

In Section 3, we propose a strategy for solving the system of non linear equations: it is a continuation method in which the viscosity parameter (for instance) is progressively decreased to the desired value, and for each value of the latter parameter, the system of non linear equations is solved by a Newton method with inner iterations of a bi-gradient like algorithm.

Finally, in Section 4, we discuss the results of numerical simulations in two cases. In the first example, we show in particular that there exist multiple solutions and we discuss how the iteration count of the solver is affected by the variation of the parameters in the model. The second case is a model for a crowd of agents crossing a rectangular hall from an entrance to an exit: We consider situations in which queues occur, and we show how the solution depends on the parameters.

2.   Finite difference methods

In order to approximate (7)–(12) numerically, we propose a finite difference method reminiscent of that introduced in ^[2,3] for MFGs without coupling through the controls. The important features of such methods are as follows

● they are based on monotone finite difference schemes for (7). Hence, a comparison principle still holds at the discrete level.

● the special structure of (7) and (8) is preserved at the discrete level, namely that the FP equation can be obtained by differentiating the HJB equation w.r.t. $u$ and by taking the adjoint of the resulting equation. This results in a monotone approximation of (8), which ensures that the discrete version of $m$ remains non negative at all time if $m_0$ is non negative. The discrete FP equation will also preserve mass.

To simplify the discussion, let us focus on the case when $d = 2$ and $\Omega = (0, 1)^2$. More complex domains (even with holes) can be handled by the present method, but an additional effort would be needed if the domain had boundaries not aligned with the axes of the underlying grid that will be introduced soon. We also assume for simplicity that the boundary conditions are of the type (13) and (14) on the whole $\partial \Omega$.

2.1.   Notations and definitions

For two positive integers $N_T, N_h$, we set $\Delta t = T/N_T$, the time step, and $h = 1/N_h$, the step related to the state variables. Consider the set of discrete times $\mathcal{T}_{\Delta t} = \left\{ {{t_k} = k\Delta t, k = 0, \ldots, {N_T}} \right\}$ and the grid ${\Omega _h} = \left\{ {{x_{i, j}} = (ih, jh), i, j = 0, \ldots, {N_h}} \right\}$.

The goal will be to approximate $u \left(t_k, x_{i, j} \right)$ and $m \left(t_k, x_{i, j} \right)$ respectively by $u^k_{i, j}$ and $m^k_{i, j}$, for all $k\in \{0, \dots, N_T\}$ and $(i, j)\in \{0, \dots, N_h\}^2$, by solving the discrete version of (7)–(14) arising from the finite difference scheme.

Let us define the discrete time derivative of $y: \mathcal{T}_{\Delta t}\times\Omega_h\rightarrow \mathbb{R}$ by the collection of real numbers:

for $k = 0, \dots, N_T-1$, $i, j = 0, \dots, N_h$.

Given a grid function $y:\Omega_h\rightarrow \mathbb{R}$, we introduce the first order right sided difference operators

and the five point discrete laplace operator:

for all $(i, j)$ such that $0\le i, j\le N_h$. Note that the definition of these operators at boundary nodes of $\Omega_h$ needs to extend the grid function $y$ on a layer of nodes outside $\Omega_h$. This is done by using discrete versions of the Neumann conditions (13)–(14): assume that the boundary condition for $y$ is of the type $\frac{\partial y}{\partial n} = z$, where $z$ is a given continuous function. Let $(z^k_{i, j})$ is a suitable grid function approximating $z$. We then choose the discrete version of the latter equation (first order scheme):

for $i, j = 0, \dots, N_h$ and $k = 0, \dots, N_T$.

Remark 2.1. We have chosen a first order scheme for the boundary condition in order to preserve the monotonicity of the discrete Hamiltonian at boundary nodes, see later. Since the overall scheme is monotone as it was already mentioned above, thus first order, it would be pointless to choose a higher order scheme for the boundary conditions.

Remark 2.2. Note that the previously mentioned discrete operators are not needed at the nodes of $\partial \Omega_h$ at which Dirichlet boundary conditions are imposed.

For a grid function $V:\Omega_h\rightarrow \mathbb{R}$, we define the value of $V$ at half-integer indices by linear interpolation:

In order to define the Godunov scheme for (7), we introduce the map $\Phi$: for two grid functions $y:\Omega_h\to \mathbb{R}$ and $V: \Omega_h\to \mathbb{R}^2$, the grid function $\Phi(y, V): \Omega_h\to \mathbb{R}^4$ is defined by

Finally, let us introduce the function $g: \mathbb{R}^4\times \mathbb{R}^2\to \mathbb{R}$:

2.2.   The scheme

With the ingredients defined in paragraph 2.1, we are ready to propose the discrete version of (7)–(10):

for $i, j = 0, \dots, N_h$, $k = 0, \dots, N_T-1$, where the discrete transport operator $\mathcal{T}$ is defined by

for any $i, j = 0, \dots, N_h$. Note that (20) is not enough in order to fully determine $\mathcal{T}(\Phi(u, V), m)$ at the boundary nodes. We also need the following quantities:

for $i, j = 0, \dots, N_h$, where the last identity in each line comes from (30) below. The discrete version of (11) and (12) is

The discrete version of (13)–(14) is

for $i, j = 0, \dots, N_h$ and $k = 0, \dots, N_T-1$.

Note that (22) is an implicit scheme for (7), (recall that (7) is backward w.r.t. time), whereas (23) is an implicit scheme for (8), (recall that (8) is forward w.r.t. time). This explains why no restriction is made on the time step.

The discrete Hamiltonian introduced in (21) $g: \mathbb{R}^4\times \mathbb{R}^2\to \mathbb{R}$, $(q, V)\mapsto g(q, V)$, has the following properties:

H1 (monotonicity) if $q = (q_1, q_2, q_3, q_4)$, then $g$ is nonincreasing with respect to $q_1$ and $q_3$ and nondecreasing with respect to $q_2$ and $q_4$

H2 (consistency) $g \left(q_1, q_1, q_3, q_3, V \right) = H \left(q_1, q_3, V \right)$

H3 (regularity) $g$ is $C^1$.

Remark 2.3. We aim at using Newton iterations in order to solve the system of nonlinear equations arising from the discrete scheme. For that purpose, we need to linearize the discrete version of the Fokker-Planck equation. Therefore, for another small positive parameter ${\varepsilon}$, we may replace the definition (20) of $\Phi$ by the following:

where the $C^1$ approximation $a\mapsto a_{+, {\varepsilon}}$ of $a\mapsto a_+ = a {\bf{1}}_{a > 0}$ is defined by

and $a_{-, {\varepsilon}} = a_{+, {\varepsilon}}-a$. We may modify (28) and (31) accordingly.

For fixed $u$ and $V$, $m$ can be found by solving at each time step a system of linear equations; the matrix of the system has positive diagonal entries and non positive off diagonal entries. This is enough to prove that the positivity of $m$ is preserved, see ^[2,3]: if $m^0$ is a non negative grid function, then so is $m^k$ for all $k\le N_T$. It is also possible to prove that the mass is conserved:

Lemma 2.1. Assume that $(u, m, V)$ satisfies (22)–(26), (29), and (30), (31).

The total mass of $m$ is preserved, i.e.,

for any $k = 0, \cdots, N_T$.

Proof. From (17), we deduce that

Then we sum (27) for $i, j = 0, \dots, N_h$:

From (20), (28), (30) and (31), the sums of the right hand sides of (36) and (35) vanish if $q = \Phi \left(u^k, V^k \right)$, for any $k = 0, \dots, N_T-1$.

From the observations above, summing (23) on $i, j = 0, \dots, N_h$ yields

then the desired result.

Neither the proof of existence for the discrete problem (22)–(26), (29), and (30), (31) nor the convergence properties of the scheme will be discussed in the present work which is mainly devoted to numerical results. They will be addressed in a forthcoming work.

For the existence of solutions of the discrete problem, a fixed point argument should be used. It requires a priori estimates for $u$. Obtaining these estimates is the most technical and difficult part in the analysis in the continuous case, see ^[16]: the key step consists of obtaining an estimate on $\| \nabla_x u\|_{\infty}$ by Bernstein method. This difficulty also exists at the discrete level, and the situation may be worse, because adapting Bernstein estimates for solutions of finite difference schemes is well known to be quite tricky when it is possible.

However, if we know a priori that in the continuous problem, $u$ is Lipschitz continuous with respect to the state variable uniformly in time, then, for $M$ large enough, we may obtain an equivalent problem

● either by modifying (3), i.e the definition of $V$, as follows:

● or by modifying (9) as follows:

where $\Psi_M : \mathbb{R}_+\to \mathbb{R}_+$ is a smooth, non increasing cut-off function such that

Note that the second truncation (37) is more convenient for us, since it is made at the level of the system of PDEs. The analysis of the discrete version of the new system then becomes much easier since the nonlocal nonlinearity becomes bounded uniformly with respect to $u$ and $m$. It is then possible to use the arguments contained in e.g., ^[2,3] and obtain the existence of solutions for the discrete problem. In this setting, it seems also possible to obtain convergence of sequences of solutions of the discrete problem to a weak solution of the continuous problem, for sequences of grid parameters tending to $0$.

2.3.   Solving the discrete version of the Hamilton-Jacobi-Bellman equation

For brevity, we will use the notation $(y_{i, j}^k)$ for a grid function, omitting that the indices $i, j$ take their values in $\{0, \dots, N_h\}$ and that $k$ takes its values in $\{0, \dots, N_T\}$. This paragraph is devoted to solving the system of nonlinear equations satisfied by the grid function $(u_{i, j}^k)$, given the grid functions $(m^k_{i, j})$ and $(V_{i, j}^k)$. Let ${\widetilde f} = ({\widetilde f}^k_{i, j})$ denote the grid function defined by

We may then rewrite (22), (30) and the first identity in (29) in the compact form

Finding $u$ given $m$ and $V$ amounts to solving a discrete version of a nonlinear parabolic equation posed backward in time with Neumann boundary conditions. This is much simpler than solving the complete forward-backward system for $(u, m, V)$, because a backward time-marching procedure can be used. Since, as it was already observed above, the scheme is implicit, each time step consists of solving the discrete version of a nonlinear elliptic partial differential equation which is local because $V$ is given. Starting from the terminal time step $N_T$ for which (29) gives an explicit formula for $u^{N_T}$, the $k$-th step of the backward loop consists of computing $u^k$ by solving (22) and (30) given $u^{k+1}$, $m^{k+1}$ and $V^k$. This is done by means of Newton iterations.

For completeness, let us give a few details on Newton iterations. We introduce the operator $\mathcal{R}_u$,

for $u', v', f':\Omega_h\rightarrow \mathbb{R}$, $V':\Omega_h\rightarrow \mathbb{R}^2$, and $i, j = 0, \dots, N_h$. Note that in (40), $u'$ is extended outside $\Omega_h$ thanks to (30). We aim at approximating the solution of

Starting from an initial guess noted $u^{k, 0}:\Omega_h\rightarrow \mathbb{R}$, the Newton iterations consist of computing by induction a sequence $u^{k, \ell}$ of approximations of the solution to (41): given $u^{k, \ell}$, $u^{k, \ell+1}$ is found by solving the system of linear equations

In the latter equation, the Jacobian matrix $d_u \mathcal{R}_u \left(u^{k, \ell}, u^{k+1}, {\widetilde f}^k, V^k \right)$ is sparse since the PDE is local, and invertible since the scheme is monotone. Note that $d_u \mathcal{R}_u \left(u^{k, \ell}, u^{k+1}, {\widetilde f}^k, V^k \right)$ depends neither on $u^{k+1}$ nor on ${\widetilde f}^k$, so we can write it $d_u \mathcal{R}_u \left(u^{k, \ell}, \cdot, \cdot, V^k \right)$. The system of linear equations is solved by using special algorithms for sparse matrices. In the numerical simulations presented below, we use the C-library UMFPACK, see ^[1], implementing unsymmetric multifrontal method and direct sparse LU factorization.

Note that one may choose the initial guess $u^{k, 0} = u^{k+1}$. The Newton iterations are stopped when the residual $| \mathcal{R}_u(u^{k, \ell}, u^{k+1}, \tilde f^k, V^k)|$ is small enough, say for $\ell = \ell_0$. Then we set $u^{k} = u^{k, \ell_0}$. It is well known that the Newton algorithm for (22) is equivalent to an optimal policy iteration algorithm, that it is convergent for any initial guess, and that the convergence is quadratic.

2.4.   Solving the discrete version of the Fokker-Planck-Kolmogorov equation

This paragraph is devoted to solving the system of linear equations satisfied by the grid function $(m_{i, j}^k)$, given the grid functions $(u^k_{i, j})$ and $(V_{i, j}^k)$. We may then rewrite (23), (31) and the second identity in (29) in the compact form

Finding $m$ given $u$ and $V$ amounts to solving a discrete version of an linear parabolic equation posed forward in time with Neumann boundary conditions. Since the scheme is implicit, each time step consists of solving the discrete version of a linear elliptic partial differential equation which is local. Starting from the terminal time step $0$ for which (29) gives an explicit formula for $m^{0}$, the $k$-th step of the forward loop consists of computing $m^{k+1}$ by solving (23), (31) given $m^{k}$, $u^{k+1}$ and $V^k$. It is easy to see that the matrix of the latter system of linear equations is exactly the transposed of $d_u \mathcal{R}_u \left(u^{k}, \cdot, \cdot, V^k \right)$ which has been introduced in the previous paragraph. Therefore, it is sparse and invertible. In the numerical simulations presented below, we use the C-library UMFPACK again.

3.   Newton algorithms for solving the whole system (22)–(31)

The solution of the system of nonlinear equations (22)–(31) is not easy because

$1).$ the system is the discrete version of a system of forward-backward equations, which precludes a simple time marching algorithm

$2).$ The easiest instances of MFGs can be interpreted as optimal control problems driven by a partial differential equation, which opens the way to algorithms based on the variational structure. In the MFGC considered here, there is no underlying variational principle and these algorithms cannot be used.

Following previous works of the first author, see ^[2], we choose to use a continuation method (for example with respect to the viscosity parameter $\nu$) in which every system of nonlinear equations (given the parameter of the continuation method) is solved by means of Newton iterations. With Newton algorithm, it is important to have a good initial guess of the solution; for that, we take advantage of the continuation method by choosing the initial guess as the solution obtained with the previous value of the parameter. Alternatively, we have sometimes taken the initial guess from the simulation of the same problem on a coarser grid, using interpolation.

It is also important to implement the Newton algorithm on a "well conditioned" system. Therefore, we shall not directly address (22)–(31), but we shall rather eliminate the unknowns $u$ and $m$ by using the time-marching loops described in paragraphs 2.3 and 2.4 and see (22)–(31) as a fixed point problem for $(\tilde f, V)$.

Before describing the algorithm, we need to provide a numerical method for obtaining the average drift given $u$ and $m$, i.e., for solving (24)-(25) at least approximately.

3.1.   The coupling cost and the average drift

Let us introduce $\mathcal{F}_f$ which maps the grid function $m$ defined on $\mathcal{T}_{\Delta t}\times \Omega_h$ to the grid function ${\widetilde f}$ given in (38). We also define the maps $\mathcal{Z}$ and $\mathcal{V}$ by

for $u', m':\Omega_h\rightarrow \mathbb{R}$, $V':\Omega_h\rightarrow \mathbb{R}^2$ and $i, j = 0, \dots, N_h$.

It can be proved exactly in the same manner as in the continuous case, see [16,Lemma 2.4], that $V'\mapsto \mathcal{V}(u', m', V')$ is a contraction in the maximum-norm for instance, with a contraction factor $| {\lambda}| {\theta}$.

For a positive integer $L$, we define the map $\mathcal{F}_V$ by :

where ${\widetilde V}^{k, 0} = V^k$ and the sequence ${\widetilde V}^{k, \ell}$ is defined by the following induction:

Remark 3.1. 1). With a slight abuse of notation, if $L = \infty$, then ${\widetilde V}^k$ is the fixed point of the map $V'\mapsto \mathcal{V} \left(u^k, m^{k+1}, V' \right)$.

2). The induction (48) corresponds to Jacobi iterations and can be easily parallelised. It is also possible to implement Gauss-Seidel iterations, which are a little more complex to write. They consist of using the components of ${\widetilde V}^{k, \ell}$ as soon as they are obtained (instead of those of ${\widetilde V}^{k, \ell-1}$) in order to compute the components of ${\widetilde V}^{k, \ell}$ which have not been obtained yet. In our implementation, we have in fact used the Gauss-Seidel iterations with a lexicographic ordering of the components of the grid functions.

Remark 3.2. As we shall see in paragraph 4.1.4, the choice of $L$ does not impact the convergence of the overall iterative algorithm. Therefore, a good choice turns out to be $L = 1$.

3.2.   The linearized operators

3.2.1.   Notation

Let $u, m, f: \mathcal{T}_{\Delta h}\times \Omega_h \to \mathbb{R}$ be generic grid functions standing respectively for discrete versions of the value function, the law of the distribution of states, the right hand side of the discrete HJB equation (22). Let $V: \mathcal{T}_{\Delta h}\times \Omega_h\to \mathbb{R}^2$ be a generic grid function standing for the average drift. Let us introduce the operators obtained by differentiation of the maps $\mathcal{F}_u$, $\mathcal{F}_m$ and $\mathcal{F}_V$:

In the three paragraphs below, we explain how these differential operators can be computed.

3.2.2.   Linearized Hamilton-Jacobi-Bellman equation

The variation $du$ of $u = \mathcal{F}_u (f, V)$ induced by variations $df$ and $dV$ of $f$ and $V$ is given by

More explicitly, $du$ is obtained by linearizing (22), (30) and (29). It satisfies

for $i, j = 0, \dots, N_h$ and $k = 0, \dots, N_T$. The first (respectively second) inner product appearing in (51) takes place in $\mathbb{R}^4$ (respectively $\mathbb{R}^2$). The set of equations (51) is supplemented with the terminal condition

and the boundary conditions

Given $df$ and $dV$, the variation $du$ is found by solving (51), (52) and (53). This is done by marching backward in time and at each time step solving a system of linear equations with a sparse and invertible matrix (of the same form as in the Newton iterations described in 2.3). Again, we use the library UMFPACK for that.

3.2.3.   Linearized Kolmogorov-Fokker-Planck equation

The variation $dm$ of $m = \mathcal{F}_m (u, V)$ induced by variations $du$ and $dV$ of $u$ and $V$ is given by

The grid function $dm$ is obtained by linearizing (23), (31) and (29). It satisfies

for $i, j = 0, \dots, N_h$ and $k = 0, \dots, N_T$. The set of equations (55) is supplemented with the initial condition

and the boundary conditions

Here again, (55), (56) and (57) is solved by marching in time and solving a system of linear equations at each time step (using UMFPACK).

3.2.4.   Linearized coupling costs and average drifts

The variation of ${\widetilde f} = \mathcal{F}_f(m)$ induced by a variation $dm$ of $m$ is obviously given by $d {\widetilde f}_{i, j}^k = cdm^k_{i, j}$, hence $C_{f, m}dm = cdm$.

Let us turn to the variation of ${\widetilde V} = \mathcal{F}_V(u, m, V)$ induced by the variations $du$, $dm$ and $dV$ of $u$, $m$ and $V$. Differentiating (46)–(48) leads to

see (48) for the definition of ${\widetilde V}^{k, L}$.

This permits to compute

To sumarize,

3.3.   The algorithm for solving (22)–(31)

We see (22)–(31) as a fixed point problem for the pair $(f, V)$, which we write

where $\mathcal{G}_f, \mathcal{G}_V$ are defined by

Remark 3.3. Given $u$ and $m$, $V\mapsto W = \mathcal{F}_V \left(u, m, V \right)$ is obtained as follows: for each $k$, computing $W^k$ consists of iterating $W^k\times \mathcal{V}^L(u^k, m^{k+1}, W^k)$ $L$ times starting from $W^k = V^k$. Recall that $\mathcal{V}(u^k, m^{k+1}, \cdot)$ is in (44) and (45), and is a contraction with a unique fixed point. Therefore, $\mathcal{F}_V \left(u, m, \cdot \right)$ has also a unique fixed point which does not depend on $L$. Since $\mathcal{F}_u$ and $\mathcal{F}_m$ do not depend on $L$, neither do the the solutions of (58).

The Newton iterations involve the Jacobian $A(f, V)$ of the map $(f, V)\mapsto \Bigl(\mathcal{G}_f(f, V)), \mathcal{G}_V(f, V)\Bigr)$. We set

and the blocks of the Jacobian $A(f, V)$ are defined by:

where $B_{u, f}$, $B_{u, V}$, $B_{m, u}$, $B_{m, V}$, $C_{V, u}$, $C_{V, m}$, $C_{V, V}$ are given by (49) with $u = \mathcal{F}_u(f, V)$ and $m = \mathcal{F}_m(u, V)$. In an equivalent manner, we can write

or

Every Newton iteration for solving (58) consists of solving a system of linear equations of the form

In our implementation, this system is solved iteratively by BiCGStab algorithm, see ^[21]. Note that BiCGStab algorithm only requires a function that computes

for any grid functions $\hat f: \mathcal{T}_{\Delta t}\times \Omega_h \to \mathbb{R}$ and $\hat V: \mathcal{T}_{\Delta t}\times \Omega_h \to \mathbb{R}^2$. Ths is done using (67) and combining the methods described in paragraph 3.2. The assembly of the Jacobian matrix is not needed.

4.   Numerical simulations

We are going to discuss the results of numerical simulations in two cases, both related to crowd motion.

4.1.   First example

4.1.1.   Description of the model

The state space is the square $\Omega = (-0.5, 0.5)^2$ and the time horizon is $T = 1$.

We consider the MFGC described by (7)–(14), in which

$1).$ The kernel $k$ in (9) is radial, i.e., of the form $K(x, y) = K_\rho(|x-y|)$. Here $r\mapsto K_\rho(r)$ is a non-increasing $C^1$ function defined on $\mathbb{R}_+$, with

for a positive number $\rho > 0$, which is the radius of the disc in the state space that a reprensative agent uses for computing the average of the controls.

$2).$ The parameter $c$ in (7) is chosen as $c = 10^{-3}$. Recall that the cost $cm(t, x)$ reflects the aversion of a representative player to crowded regions of the state space.

$3).$ At time $t = 0$, the agents are distributed in the top-right and the bottom-left corners of the domain $\Omega$, see the left part of Figure 1. The density is piecewise constant, with values appearing on Figure 1.

$4).$ The terminal cost is also piecewise constant. It takes a small value in the top-left and the bottom-right corners of the domain $\Omega$, see the right part of Figure 1 for the chosen values. This cost attracts the agents to the latter two corners of $\Omega$.

$5).$ The function $f_0$ in (7) is chosen to be proportional to the terminal cost, with a factor $0.1$. This term is linked to the running cost and has the same effect as the terminal cost.

$6).$ Recall that the boundary conditions (13)–(31) rule out the entry or exit of agents. The total mass of $m$ is conserved.

Remark 4.1. Note that the problem is symmetric with respect to the two diagonals of the square domain. The grid of the domain is chosen to have the same symmetry.

Remark 4.2. The two parts of the running cost: $cm(t, X_t)$ and $\frac{\theta}2 \left| {\alpha}- {\lambda} V(t, X_t) \right|^2$ may have competing effects: indeed, the former cost may incitate the agents to spread all over the state space, whereas the latter may result in the agents selecting the same controls therefore staying grouped. Although the constant $c = 10^{-3}$ seems small, it is chosen in such a way that the above-mentioned two costs have the same orders of magnitude for e.g., $\theta = 1$ and $\lambda = 0.9$.

4.1.2.   The results of the simulations: discussion

In the present simulation, we choose the parameters as follows:

On Figure 2, we display snapshots of $m$, $u$ and the optimal feedback law at several times.

Since both the problem and the grid are symmetric with respect to the two diagonals, $m$, $u$ and the feedback have the same symmetry for all times. Let us describe the evolution of the distribution of states, that we can name "gathering-kissing-splitting" referring to the "drafing-kissing-tumbling" phenomenon in fluid mechanics (for the interaction of particles in a fluid).

● The term $cm(t, X_t)$ in the running cost prevents the part of the distribution initially supported in one of the opposite corners (say the bottom-left corner) to travel directly to one of the targets (say the bottom-right corner). On the other hand, the interaction through controls prevents this part of the distribution to split into two equal parts which would travel directly to the two targets, because the agents favor controls close to the local average. Therefore, the part of the distribution initially supported in one of the opposite corners first travels to the center of the domain. This is the gathering phase of the evolution.

● The two parts of the distribution reach the center of the domain, where the local average of the velocity becomes small; the dominating cost then becomes the one which attracts the agents to the bottom-right and top-left corners. Because of the repulsive effect due to the part $cm(t, X_t)$ of the running cost, the part of the distribution initially supported in one of the opposite corners (say the bottom-left corner) and having traveled toward the center of the domains splits into two parts which make each a ninety-degrees turn. This is the kissing phase of the evolution.

● After the kissing phase, the distribution of states splits into parts which travel to the targets.

Finally, we see that the paths followed by the agents is far from being a shortest path to the targets.

Let us compare these results with a simulation of the MFG obtained by cancelling either $\theta$ or $\lambda$ in (7)–(14) while keeping all the other parameters and the grid unchanged. Since the problems remains symmetric with respect to the two diagonals, the obtained results have the same symmetry. On Figure 3, we display snapshots of $m$, $u$ and the optimal feedback law at several times. We see that the evolution is quite different from that displayed on Figure 2, since the two parts of the distribution initially supported in two opposite corners of the domain symmetrically split into two parts each, which travels directly to the targets. The initial splitting phase is due to the coupling cost $cm(t, X_t)$ which favors the spread of the distribution. The path followed by the agents is close to a shortest path to the targets.

4.1.3.   Non-uniqueness of solutions

The solution of the MFGC displayed on Figure 2 is likely to be unstable because the paths followed by the agents are significantly longer than the shortest paths to the targets. We expect that there are other solutions. We are going to see that this is indeed the case. For that purpose, we are going to introduce a vanishing perturbation of the initial distribution which breaks the symmetry of the problem. This will lead to different solutions to (7)–(14). An example of perturbation is displayed on Figure 4. It consists of adding very little mass at the top and the bottom of the domain; the perturbed distribution is no longer symmetric with respect to the diagonals. We expect that the agents initially distributed near the bottom target will go to the right, and that all the agents initially distributed at the bottom of the domain will follow them, because of the interaction through controls. Similarly, we expect that all the agents initially located at the top of the domain will go to the left.

In our simulations, we use a continuation method, i.e., we consider perturbations corresponding to a decreasing sequence of nonnegative parameters $(\pi_n)_{n\in 0, \dots, N}$. The last value $\pi_N = 0$ corresponds to the distribution displayed on the left of Figure 1. For each new value $\pi_{n+1}$ of the parameter, we initialize the Newton iterations described in Section 3 by the solution corresponding to the preceding value $\pi_{n}$. For positive values of $\pi_n$, the simulated solution is not symmetric with respect to the diagonals, and this property is preserved for the last value $\pi_{N} = 0$.

On Figure 5, we display three different solutions (the distribution of states at different times) obtained with three different sequences of vanishing perturbations of the initial distribution displayed on Figure 1: the solution displayed on the left corresponds to the sequence of perturbations displayed on Figure 4; applying to the pertubation a symmetry with respect to a diagonal, we obtain the solution displayed on the center; the solution displayed on the right is obtained by another kind of perturbation symmetric with respect to one diagonal only, located at the top and at the left of the domain: for this solution, we again notice gathering and kissing phases, and that the paths followed by the agents deviate from the shortest ones.

4.1.4.   Behaviour of the algorithm

In this paragraph, we investigate how the behaviour of the algorithm described in Section 3 is affected by the variations of the parameters in the model. Recall that the algorithm is based on Newton iterations with a BiCGStab inner loop at each step.

We start by the grid size and the viscosity parameter $\nu$: in the experiments reported below, the stopping criterion for the outer Newton iteration is that the normalized Euclidean norm of the residual is smaller than $10^{-8}$. The stopping criterion for the inner BiCGStab iterations is that the ratio of the norms of the current and the initial residuals is smaller than $10^{-7}$ (hence, it involves a relative error). The parameters of the model are set to

Table 1 displays the iterations counts for the outer Newton and the inner BiCGStab loops for different viscosities and grid sizes. We notice that the number of iterations is not very sensitive to the grid size, as expected because the map $\mathcal{G}-I_d$, with $\mathcal{G}$ defined in (59) and (60), is the discrete version of a map which has compactness properties. The iteration counts increase as $\nu$ is decreased.

Table 2 displays the iteration counts for different choices of ${\lambda}$ and ${\theta}$, and fixed values the other parameters:

with a $101\times 101$ nodes grid and $101$ time steps.

In fact, the numbers in Table 2 are obtained by running a continuation method in $\theta$ and $\lambda$: for each cell of the table, the solution corresponding either to the left or the upper neighboring cell is used as an initial guess for the Newton iterations. If a cell has two such neighbors, the choice of the initial guess is made as follows: in the top-right triangular part of the table, strictly above the diagonal, we choose the initial guess corresponding to the left neighboring cell. In the bottom-left triangular part of the table, including the diagonal, we choose the initial guess corresponding to the cell immediately above.

We see that changing ${\lambda}$ and ${\theta}$ mostly impacts the number of iterations of the inner loop. A reason for that is that $V$ is obtained by solving a linear fixed point problem whose contraction factor is ${\lambda} {\theta}$. Hence, it is sensible that the number of iterations necessary to solve the systems of linear equations arising in the Newton steps increases as ${\lambda} {\theta}$ tends to $1$.

Recall that we use a continuation method, consisting of decreasing progressively $\nu$ until it reaches the desired value; in Figure 6 (respectively Figure 7), we plot the average number of BiCGStab iterations versus $\nu$, (respectively the number of Newton iterations versus $\nu$) for different values of ${\lambda}$, ${\theta}$ and $c$. The stopping criteria have been described above. The grid contains $101\times 101$ nodes and there are $101$ time steps. We recover the information contained in Tables 1 and 2, i.e., that $\nu$ mostly impacts the number of iterations in the inner loop. We also notice that when $\nu$ is large, the number of iterations seems to depend more on ${\lambda}$ and $\theta$ than on $c$ and that it becomes highly sensitive to $c$ when $\nu$ is small, but we do not really know how to explain this. The number of iterations in the outer Newton loop seems much less sensitive than the number of inner BiCGStab iterations. Looking at the overall number of BiCGStab iterations, they are not very different for the choices $({\lambda}, {\theta}, c) = (0.9, 1, 0)$ and $({\lambda}, {\theta}, c) = (0, 0, 0.1)$, since the latter case needs more Newton iterations but less BiCGStab iterations per Newton steps.

Finally, we investigate the sensitivity of the algorithm to $L$, the number of fixed point iterations yielding a proxy of $V$, see paragraph 3.1. Table 3 contains the iteration counts of the inner and outer loops for different values of $\nu$ and $L$, and two choices of $(\lambda, \theta)$ (all the other parameters are fixed $c = 10^{-3}$, $\rho = 0.2$). We notice that the choice of $L$ seems to have little impact on the iteration count. For these reasons, $L = 1$ appears to be a good choice.

4.2.   Second example

As a second example, we consider a model for a crowd crossing a hall, with an entrance at the left and an exit at the right. Roughly speaking, a generic agent will interact with those located in a cone ahead of her.

4.2.1.   Description of the model

The state space is the rectangle $\Omega = (-1, 1)\times (-0.1, 0.1)$ and the time horizon is $T = 8$. Let $\ell$ denote the length of the bottom and upper sides of $\partial \Omega$: $\ell = 2$.

The boundary of $\Omega$ is split into three parts $\Gamma_{ND}$, $\Gamma_{NN}$ and $\Gamma_{DD}$ which respectively stand for an entrance, some walls and an exit:

and we consider the following boundary conditions:

where $H$ is defined below.

The system of PDEs is still (7)–(10), with the following new features:

● compared to (5), the part ot the running cost accounting for the interaction through controls is multiplied by a positive normalization factor $a$:

and the associated Hamiltonian by the Legendre's transform is

We choose the normalization factor $a$ as follows:

where ${\widetilde a}$ is a positive constant (independent of $\lambda$ or $\theta$). As we shall see later, the normalization allows us to compare the solutions obtained with different values of ${\lambda}$ and ${\theta}$. We shall also see that this normalization is specific to the example under consideration and may not be relevant in other situations.

● The kernel $K$ is given by

where $\omega$ is the angle made by the vector $y-x$ with the vector $(1, 0)$ and

$1). \; k_\rho$ is defined as in the first example

$2). \; \kappa$ is a nonincreasing $C^1$ function defined such that

for a given angle $\omega_0\in (0, \frac \pi 2)$.

● The parameter $c$ in (7) is chosen to be zero

● We shall consider various functions $f_0$ modeling the cost for staying at a given point in the domain.

The initial condition is $m_0(x) = 10^{-4}$ for all $x\in \Omega$, and the terminal cost is $0$.

4.2.2.   The case when $f$ is constant: comparisons with a one-dimensional problem

In this paragraph, we choose $f_0(x) = F$ for all $x\in \Omega$.

Simplification: a one-shot one-dimensional game We approximate the MFGC by a one-shot one-dimensional game. The state space is the interval $[0, \ell]$. The points $x = 0$ and $x = \ell$ respectively stand for an entrance and an exit. When she enters the domain, a representative agent chooses her drift $\alpha\in \mathbb{R}_+$ once and for all, and her dynamics is given by $x'(t) = {\alpha}$ for $t = \ell/\alpha$. The distribution of drifts is a probability measure $\mu$ on $\mathbb{R}_+$ and the average drift is $V = \int_{\alpha\in \mathbb{R}_+} \alpha d\mu(\alpha)$. There is no entry or exit costs.

With the running cost $L(\alpha, V)+F$, where $L$ is given by (74), the mean field Nash equilibrium reads

Given $V$, one checks that the unique solution of the minimization problem in $\alpha$ is

Hence $\mu = \delta_{\alpha ^*(V)}$, and the static mean field Nash equilibrium reads

which yields

This leads us to choose $a = \widetilde a / (1-\theta \lambda^2)$, which yields $\alpha^* = V = \sqrt{\frac {2F } {\widetilde a}}$. Note that the value corresponding to the mean field Nash equilibrium is

Remark 4.3. Note that if $\theta = \lambda = 1$, then (78) implies that $\alpha^* (V) > V$, while an equilibrium corresponds to $V = \alpha^*(V)$. Thus, there is no mean field equilibrium if $\theta = \lambda = 1$.

Remark 4.4 (Comparison with the one-shot mean field type control problem). The Nash equilibrium must be distinguished from the situation in which the drift of the agents is determined in order to minimize the global cost; the latter problem belongs to the class of mean field type control problems (MFTC for short)) studied by A. Bensoussan and his collaborators, ^[6] and R. Carmona and F. Delarue, ^[9].

In this case, the problem reads

Looking for the solution as the Dirac mass $\delta_{\alpha^{**}}$, we see that $\alpha^{**}$ minimizes

thus

The value of the problem is

from Young's inequality $y\leq \frac{y^2}{2z}+\frac{z}2$, with $y = \sqrt{1-\lambda\theta(2-\lambda)}$ and $z = \sqrt{1-\lambda^2\theta}$. The value of the MFTC is lower than the value of the MFG, which is natural since in the former case, the agents are collaborating to minimize the global cost. Moreover, the values of the two problems coincide only if $\lambda \theta = 0$, i.e., when there are no interactions through controls.

Comparison with the one-shot one-dimensional problem We have made numerical simulations of the MFGC described in paragraph 4.2.1 with $f_0(x) = F = 1$ for all $x\in \Omega$, $a = {\widetilde a}(1- {\lambda}^2 {\theta})^{-1}$ and ${\widetilde a} = 2$. We wish to compare the results with the explicit formulas obtained above for the one-shot one-dimensional MFG. If the latter problem is a good approximation of the former one, we should find $\alpha_1$ close to $\sqrt{\frac{2F}{ {\widetilde a}}} = 1$, therefore mostly independent of $\lambda$ and $\theta$.

Note that the approximation by the one-shot one-dimension mean field game is sensible only if the time remaining to the horizon is large enough, i.e., significantly larger than $\ell/\alpha^* = \ell$. If $T-t$ is small, then the optimal strategy for a representative player located on the left part of the domain is not to move.

On Figure 8, we display the norm of the optimal feedback for several values of $\lambda$ and $\theta$. We see that the optimal drift is close to $1$ in agreement with the explicit formula found for the simplified one-shot one-dimensional MFG. We also see that the agents located at the bottom and top of $\Omega$ head toward the center of the domain (i.e., $x_2 = 0$) before reaching the exit. Thus, the second coordinate of $\alpha$ gets large and there are singularities on both sides of the exit, whose amplitude increases as $\lambda\theta$ tends to $1$.

4.2.3.   Queues

Here, we choose $f$ as follows:

hence, an agent pays a cost for staying in $\Omega$ which is higher in the left part of the domain.

The parameters are

The grid has $101 \times 21$ nodes and there are $101$ time steps.

The evolution of the distribution of states and of the optimal feedback is displayed on Figure 9. We compare these results with a simulation of the MFG obtained by cancelling either $\theta$ or $\lambda$ while keeping all the other parameters and the grid unchanged, see Figure 10. On Figure 9, we see that a well distributed queue takes place from the entrance to the exit, by contrast with Figure 10 where we see that the agent rush and accumulate in the right part of the domain. Similarly, the deceleration is much stiffer in the latter case. It is not surprising that the interactions through controls have the effect of smoothing the distribution of states and the optimal feedback law. On the bottom of Figure 10, we see that when $t$ is close to the horizon, the distribution is mainly concentrated near the middle of the domain but slightly on the right: this corresponds to the agents that have reached the zone where $f_0$ has the smaller value, i.e., $1$, but for which reaching the exit before $T$ becomes too costly. There is also a smaller bump near the entrance corresponding to agents that would pay too high a cost to reach the right part of the domain before $T$. These phenomena are clearly a side effects due to the finite horizon. These side effects are also present on Figure 9, but they are attenuated by the interaction through the controls.

Remark 4.5. In models accounting for congestion effects, the cost of motion depends on the density of the distribution of states and gets higher in the crowded regions. We refer to ^[20] for a pioneering discussion of MFG models including congestion, to ^[5] for the analysis of the system of PDEs that arise with these models, and to ^[4] for numerical simulations. Such models also permit to describe queueing phenomena, because the agents located in crowded regions pay a large cost for moving. In ^[4], since the congestion effects are taken into account in a local manner, the queues take place in small regions of the state domain. By contrast, the present nonlocal model accounts for the fact that the agents anticipate low speed regions, making the traffic more fluid and the distribution of state smoother. Although it is quite possible to do, we have not incorporated congestion effects in the present model. We plan to do it in forthcoming works.

4.2.4.   Stationary regime

We now look for a stationary equilibrium. For that, we use an iterative method in order to progressively diminish the effects of the initial and terminal conditions: starting from $(u^0, m^0, V^0)$, the numerical solution of the finite horizon problem described above, we construct a sequence of approximate solutions $(u^\ell, m^\ell, V^\ell)_{\ell \ge 1}$ by the following induction: $(u^{\ell+1}, m^{\ell+1}, V^{\ell+1})$ is the solution of the finite horizon problem with the same system of PDEs in $(0, T)\times \Omega$, the same boundary conditions on $(0, T)\times \partial \Omega$, and the new initial and terminal conditions as follows:

As $\ell$ tends to $+\infty$, we observe that $(u^\ell, m^\ell, V^\ell)$ converge to time-independent functions. At the limit, we obtain a steady solution of

with the boundary conditions on $\partial\Omega$ given by (71)–(73). On Figure 11, we plot the simulated solutions of the stationary problem with different parameters.

Acknowledgments

This research was partially supported by the ANR (Agence Nationale de la Recherche) through MFG project ANR-16-CE40-0015-01).

Conflict of interest

The authors declare no conflict of interest.