Universal behavior in non stationary Mean Field Games

TL;DR

This paper investigates the behavior of non-stationary Mean Field Games, revealing a universal scaling solution analogous to an ergodic state, derived through a mapping to electrostatics, enhancing understanding of complex dynamic systems.

Contribution

It introduces a universal scaling solution for non-stationary Mean Field Games lacking an ergodic state, expanding theoretical understanding of their long-term dynamics.

Findings

Existence of a universal scaling solution in non-stationary MFGs

Mapping of MFG behavior to electrostatic problems

Insights into long-term dynamics without ergodic states

Abstract

Mean Field Games provide a powerful framework to analyze the dynamics of a large number of controlled objects in interaction. Though these models are much simpler than the underlying differential games they describe in some limit, their behavior is still far from being fully understood. When the system is confined, a notion of "ergodic state" has been introduced that characterizes most of the dynamics for long optimization times. Here we consider a class of models without such an ergodic state, and show the existence of a scaling solution that plays similar role. Its universality and scaling behavior can be inferred from a mapping to an electrostatic problem.

Equations

$$c[{\mathbf{a}}^{\mathbf{i}}]({\mathbf{x}}_t^i,t)$$

$$V[m_{\tau }](\mathbf{x})=U_0(\mathbf{x})+g{m}_t(\mathbf{x}).$$

$${\partial }_{t}m-\frac{1}{\mu }\nabla \left[m\nabla u\right]-\frac{{\sigma }^2}{2}\Delta m=0$$

$${\partial }_{t}u-\frac{1}{2\mu }{\left[\nabla u\right]}^2+\frac{{\sigma }^2}{2}\Delta u=V[m].$$

$$m(t,x)=\{\begin{array}{rlr}& \frac{3}{4}\frac{(z(t{)}^2-x^2)}{z(t{)}^3}& \text{if x≤z(t)}\\ & 0& \text{ otherwise}\end{array},$$

$$z(t)\sim t^{2/3}.$$

$$v(x,t)=-\frac{\nabla u}{\mu }-\frac{{\sigma }^2\nabla m}{2m},$$

$${\partial }_{t}m+\nabla (mv)=0.$$

$${\partial }_{t}v+\nabla \left[\frac{{\sigma }^4}{2\sqrt{m}}\Delta \sqrt{m}+\frac{v^2}{2}+\frac{g}{\mu }m\right]=0,$$

$$\{\begin{array}{rl}& {\partial }_{t}m+\nabla (mv)=0\\ & {\partial }_{t}v+\nabla \left[\frac{v^2}{2}+\frac{g}{\mu }m\right]=0\end{array}.$$

$$z(t)=3{\left(\frac{|g|}{4\mu }\right)}^{1/3}{t}^{2/3},$$

$$\{\begin{array}{rl}& {\partial }_t{\lambda }_{+}+\left(\frac{3}{4}{\lambda }_{+}+\frac{1}{4}{\lambda }_{-}\right){\partial }_x{\lambda }_{+}=0\\ & {\partial }_t{\lambda }_{-}+\left(\frac{1}{4}{\lambda }_{+}+\frac{3}{4}{\lambda }_{-}\right){\partial }_x{\lambda }_{-}=0\end{array}.$$

$$\{\begin{array}{rl}& {\partial }_{{\lambda }_{-}}x-{\beta }_{+}{\partial }_{{\lambda }_{-}}t=0\\ & {\partial }_{{\lambda }_{+}}x-{\beta }_{-}{\partial }_{{\lambda }_{+}}t=0\end{array},$$

$$\{\begin{array}{rl}& x-{\beta }_{+}t={\omega }_{+}\\ & x-{\beta }_{-}t={\omega }_{-}\end{array},$$

$${\partial }_{{\lambda }_{\pm }}{\omega }_{\mp }=-({\partial }_{{\lambda }_{\pm }}{\beta }_{\mp })t=-\frac{1}{4}t.$$

$${\partial }_{{\lambda }_{+}{\lambda }_{-}}\chi -\frac{1}{2({\lambda }_{+}-{\lambda }_{-})}({\partial }_{{\lambda }_{+}}\chi -{\partial }_{{\lambda }_{-}}\chi )=0.$$

$${\partial }_{\xi \xi }\chi +{\partial }_{\eta \eta }\chi +\frac{1}{\eta }{\partial }_{\eta }\chi =0,$$

$$\{\begin{array}{rl}& \eta t=-E_{\eta }\\ & 2(x-\xi t)=-E_{\xi }\end{array},$$

$$\chi (\eta ,\xi )\approx \frac{Q_0}{\sqrt{{\eta }^2+{\xi }^2}},$$

$$m(t,x)$$

$$v(t,x)$$

$$Q_0=\frac{1}{4\pi }{\int }_{S_{\tilde{t}}}(\vec{E}\cdot \vec{n})dS$$

$$Q_0$$

Full text

Mean Field Games are a powerful framework introduced about a decade ago by Lasry and Lions Lasry and Lions (2006) as an alternative approach to differential game theory when the number of agents becomes large. Their applications are numerous, ranging from finance Lachapelle et al. (2015); Lasry and Lions (2007) and economy Guéant et al. (2011); Achdou et al. (2014) to engineering sciences Mériauxi et al. (2012); Kizilkale and Malhamé (2016), and wherever one has to deal with optimization issues for many coupled subsystems. On a quite general basis, a Mean Field Game involves a set of $N$ players (or agents) which are characterized by a continuous “state variable” ${\bf X}^{i}\in\mathbb{R}^{d}$, $i=1\dots N$, which, depending on the context, may represent a physical position, the amount of resources owned by a company, the house temperature in a network of controlled heaters, etc.. These state variables evolve on a time interval $[0,T]$ according to some controlled dynamics, which we assume here to be described by a linear, d-dimensional, Langevin equation, $d{\bf X}^{i}_{t}={\bf a}^{i}_{t}dt+\sigma d{\bf W}^{i}_{t}$, where each component of ${\bf W}^{i}$ is an independent white noise of variance one, $\sigma$ is a constant, and the “control parameter” is the velocity ${\bf a}^{i}_{t}$. This control is adjusted in time by the agent $i$ in order to minimize a cost functional which in the simplest case can be assumed of the form

[TABLE]

In (1), $V[m_{t}]({\bf x})$ is a functional of the empirical density $m_{t}({\bf x})=\frac{1}{N}\sum_{j}\delta({\bf x}-{\bf X}^{j}(t))$, through which the agents’ optimization problems are coupled. We shall assume $V[m_{\tau}]({\bf x})$ takes the simple form

[TABLE]

For a very large number of players, like in a mean field theory, the fluctuations of the empirical density are neglected and $m_{t}({\bf x})$ becomes a deterministic quantity governed by a Fokker Planck equation. Furthermore, the optimization problems decouple and the optimal value of the cost (1) (the “value function”) for the agent $i$ becomes a function of its state variable ${\bf x}^{i}_{t}$ at time $t$, solution of an Hamilton-Jacobi Bellman equation Bertsekas (2017). The resulting Mean Field Game is defined as a pair of coupled equations, a forward diffusion for the density $m({\bf x},t)$ and a backward optimization for the value function $u({\bf x},t)$, which, in the simple case we consider here takes the form Ullmo et al. (2019)

[TABLE]

The coupling between the two PDE’s comes from two parts: in the Fokker-Planck equation (3), the optimal velocity appears in the drift term and here is proportional to the gradient of the value function, ${\bf a}=-\frac{1}{\mu}\nabla u$; in the Hamilton-Jacobi Bellman equation (4) the term $V[m]$ reflects the dependence of the cost functional (1) on the density. This structure also induces rather atypical boundary conditions: the (forward) Fokker-Planck equation is associated with an initial condition $m(0,{\bf x})=m_{0}({\bf x})$ specifying the initial distribution of agents, while the terminal cost in (1) imposes a final condition for the value function, $u(T,{\bf x})=c_{T}({\bf x})$. This forward-backward structure together with mixed initial-final boundary conditions leads to new challenges when trying to characterize, either analytically or numerically, solutions to this system of equations.

For a large class of settings, and in particular for repulsive interactions, such a system, when confined, exhibits an ergodic (stationary) state, independant on the boundary conditions, which can be rigorously defined in the limit $T\to\infty$ as an hyperbolic fixed point. The importance of this result, as proven by Cardialaguet et al. Cardaliaguet et al. (2013) is twofold: for finite but long enough optimization time $T$, the game will stay very close to this ergodic (time independent) state except possibly in its initial and final parts. Furthermore, the transient dynamics near $t=0$ and $t=T$ completely decouple one from the other, and the mixed boundary problem simplifies accordingly: they describe how the initial and final boundary conditions match with the ergodic state and are characterized by two possibly different time scales, $\tau_{0}$ and $\tau_{T}$.

If however we consider a system of agents with repulsive interactions in an unbounded state space, any initially localized configuration will expand forever and no ergodic state will ever be reached. The natural question we address here is whether some kind of limiting regime still exists for such systems that would, at some level, play a similar role as the ergodic state in confined systems.

A first indication that this is indeed the case is brought out by numerical simulations. In Fig. 1, we report the numerical evolution of a one-dimensional system with a narrow initial distribution of agents with repulsive interactions, in the absence of external potential ($V[m_{t}](x)=gm(t,x)$, $g<0$). What we observe there is that except for transient times near $t=0$ and near $t=T$, the density assumes an almost perfect “inverted parabola form” at all times

[TABLE]

where the prefactor derives from the normalization, and the time dependent scaling factor $z(t)$ grows in time as a power law,

[TABLE]

These features appear for sufficiently long optimization time $T$, and are essentially independent of the initial and final boundary conditions provided the initial distribution has a bounded extension and the final cost is close to zero everywhere.

Another necessary condition for this behavior to appear is the smallness of the healing length $\nu={\mu\sigma^{2}}/{|g|}$ with respect to other lengthscales, but it will be eventually fulfilled once the density distribution gets sufficiently extended (and at all times thereafter), possibly inducing a time shift in (6).

What those numerical results tell us is that in this configuration, the notion of ergodic state has been replaced by the next best thing, namely a universal scaling solution. The goal of this paper is to understand this puzzling result and in particular the $\frac{2}{3}$ scaling exponent.

We now introduce the few formal transformations which allow to show that this result is rather natural and intuitive, once the problem is cast in its proper language. In turn, we will also gain a better understanding of this regime

Hydrodynamic representation

The main idea underlying these transformations comes from the deep link between Mean Field Games and the nonlinear Schrödinger equation, discussed at lengths in Ullmo et al. (2019), which allows for the use of several techniques developed to study Bose-Einstein condensates. One of these techniques, the Madelung substitution, is particularly well suited to deal with the small $\nu$ regime Pethick and Smith (2008). It consists in defining a velocity field $v(t,x)$

[TABLE]

which drives the evolution of the density $m$ as a simple transport/continuity equation

[TABLE]

The evolution for the velocity field $v$ derives from the HJB equation (4)

[TABLE]

and involves a $O(\sigma^{4})$ term. As in the context of cold atoms, this term can be neglected as long as the characteristic length of the system is large in front of the healing length $\nu={\mu\sigma^{2}}/{|g|}$ Pethick and Smith (2008), leading to what is referred to as the Thomas Fermi approximation. This weak noise limit is also one of the requirements for the appearance of the parabolic behavior described in Figure 1. In this approximation the equation read

[TABLE]

These equations are formally very close to those studied for instance in the field of cold atoms Dalfovo et al. (1999), the main differences being the (negative) sign of $g$, which makes the system elliptic rather than hyperbolic, and the nature of the boundary conditions. Within this approximation, it can easily be verified that the ansatz Eq. (5) with

[TABLE]

is a particular solution of the Thomas-Fermi-like equations Eq. (10). The real mystery is therefore not that such a solution does exist, but to understand why it is universal in the large optimization time limit.

Riemann invariant and hodograph transform

To answer this question, we shall turn to an approach developed in the context of non-linear waves Kamchatnov (2000), which relies on the notions of Riemann invariants and hodograph transform. Riemann’s method can be considered an extension of the method of characteristics. It amounts to finding curves (characteristics) on which some quantities (Riemann invariants) are conserved. Here, one can show that there exists a pair of Riemann invariants, $(\lambda_{+}(x,t),\lambda_{-}(x,t))$, namely $\lambda_{\pm}=v\pm 2i\sqrt{{|g|m}/{\mu}}$, so that (10) reads

[TABLE]

Though characteristics do not exist in this context (they are curves in the complex plane $\mathbb{C}^{2}$), this change of variables still allows us to linearize these equations using an hodograph transformation Kamchatnov (2000). Taking the pair $(\lambda_{+},\lambda_{-})$ as independant variables, we express $x$ and $t$ as functions of them, so that the system (12) transforms into a linear one:

[TABLE]

where $\beta_{\pm}=\left(\frac{3}{4}\lambda_{\pm}+\frac{1}{4}\lambda_{\mp}\right)$. This system can be readily integrated once as

[TABLE]

with $\omega_{\pm}$ solution of

[TABLE]

Thus the functions $\omega_{\pm}$ can be expressed as derivatives of a potential $\omega_{\pm}=\partial_{\lambda_{\pm}}\chi$, where $\chi(\lambda_{+},\lambda_{-})$ is solution of an Euler-Poisson-Darboux equation:

[TABLE]

The main difference with the traditional treatment of NLS that we have closely followed until now is that here the Riemann invariants $(\lambda_{+},\lambda_{-})$ are complex conjugates. In terms of its real and imaginary parts, $\xi=\frac{1}{2}(\lambda_{+}+\lambda_{-})=v$ and $\eta=\frac{1}{2i}(\lambda_{+}-\lambda_{-})=\;2\sqrt{{|g|m}/{\mu}}$, equation (16) becomes

[TABLE]

which is the Laplace equation in cylindrical coordinates (with no angular dependence) with $\eta$ and $\xi$ as radial and axial coordinates, respectively. Equations (14) now read

[TABLE]

with $E_{\eta}$ and $E_{\xi}$ the radial and axial components of the electric field, ${\bf E}=-\nabla\chi$. Note that even if (17) is originally a two-dimensional problem, a clear connection with electrostatics emerges when considering it a three dimensional one with axial symmetry.

Potential representation

Through the hodograph transform we have shown that for any potential $\chi$, solution of the Laplace equation (17), there is a solution to the Thomas Fermi equation (10) provided that the relations (18) between $x$, $t$ and the electric field $\bf E$ hold. The linear Laplace equation (and the related electrostatic problem) is clearly significantly simpler than the original non-linear hydrodynamic equations. The price to pay for that simplification is that taking into account the boundary conditions becomes highly non trivial since the locus of the curves $t(\xi,\eta)=0$ or $t(\xi,\eta)=T$ on which these conditions are expressed actually depend on the particular potential $\chi(\xi,\eta)$ considered.

Since the dynamics we are interested in is related to the spreading of the density of agents, the time evolution is associated with a contraction towards the origin in the plane $(\xi,\eta)$ and the equal time curves, $t(\xi,\eta)=\text{const}$, are nested together, larger times closer to the origin. If we consider $\chi$ as generated by a distribution of charge $\rho(\xi,\eta)$, Eq. (17) implies that there is no charges between the curves $t(\xi,\eta)=0$ and $t(\xi,\eta)=T$ but $\rho(\xi,\eta)$ can be non-zero either near the origin (inside the curve $t(\xi,\eta)=T$) or at large distance (outside the curve $t(\xi,\eta)=0$). If the optimization time is long enough we can assume that there exists a large time inteval $[\tilde{t}_{\rm min},\tilde{t}_{\rm max}]$, $0\ll\tilde{t}_{\rm min}$, $\tilde{t}_{\rm max}\ll T$, so that at any time $t\in[\tilde{t}_{\rm min},\tilde{t}_{\rm max}]$, the details of the distibutions of charges both near the origin and outside the domain (using axial symmetry) can be essentially neglected, keeping only the most slowly decaying contribution. In that case a good approximation of $\chi$ is the potential created by a point charge $Q_{0}$ located at the origin

[TABLE]

with a relation between $Q_{0}$ and the boundary conditions of the problem yet to be determined.

The main result of this paper is that the approximation of the charge distribution by a monopole centered at the origin is precisely the observed universal behavior expressed by Eqs. (5)-(6), and in particular the conditions under which this approximation is valid provide the regime of validity of the scaling form Eqs. (5)-(6). Indeed, inserting Eq. (19) into Eq. (18) and inverting the relation between $(\xi\!=\!v,\eta\!=\!\sqrt{|g|m/\mu})$ and $(x,t)$ readily gives

[TABLE]

Equation (20) identifies to (5) with $z(t)$ given by Eq. (11). In addition, normalization of $m(x,t)$ would require that the value of the charge is fixed, $Q_{0}=-2|g|/\mu$. We have now to show that this is actually the case independently of the boundary condition to prove that Eq. (20) is indeed the numerically observed scaling form Eqs. (5)-(6).

For this purpose, we simply apply Gauss’s law

[TABLE]

on a surface ${S_{\tilde{t}}}$ on which time is constant, $t(\eta,\xi)=\tilde{t}$, for some $t\in[0,T]$. Using $x$ and the dumy angle $\theta$ as coordinates on this surface, the surface element$dS$ and normal vector $\vec{n}$ read, respectively $dS=\eta j(x,\tilde{t})d\theta dx$ and $\vec{n}=(n_{\xi},n_{\eta},n_{\theta})=j(x,\tilde{t})^{-1}\left(\partial_{x}\xi,-\partial_{x}\eta,0\right)$, with $j(x,\tilde{t})\equiv\sqrt{(\partial_{x}\xi)^{2}+(\partial_{x}\eta)^{2}}$. We get

[TABLE]

In the absence of mass flow at infinity, the first, time dependent, term integrates to zero. This appears here as a consequence of the fact that for all $\tilde{t}\in[0,T]$, the total charge enclosed in $S_{\tilde{t}}$ is the same and $Q_{0}$ is by construction a constant of motion. Integrating by part the second term and recalling that $\eta=2\sqrt{{|g|m}/{\mu}}$, Eq.(23) yields $Q_{0}=\frac{2g}{\mu}\int_{\mathbb{R}}mdx$ , which, because of the normalization condition $\int_{\mathbb{R}}mdx=1$, is the required result.

To conclude, we have shown that, using the potential representation of the [one dimensional potential-free] Mean Field Game problem we consider, the remarkable, and a priori quite puzzling universal scaling form which shows up in numerical simulations can be derived in a very natural way by approximating the charge distribution creating the electrostatic potential $\chi$ to a simple monopole. The condition that the observation point is sufficiently far from both the charges near the origin and the one near infinity underlying this approximation is equivalent to considering times far from both $t=0$ and $t=T$, which is of course only possible in the limit of very long optimization time. Furthermore, the electrostatic picture would allow for a more precise quantification of the conditions of validity of the approximation, using estimates on boundary conditions. We have thus been allowed to extend in some sense the notion of ergodic state to a situation where a genuine ergodic state cannot exist.

At a more general level, the potential representation of our Mean Field Games underlies the integrability of the hydrodynamical equations (10). Beyond the simple monopolar approximation Eq. (19), we can construct a complete multipolar expansion for the potential $\chi$, and it can be seen that each “charge” in this expansion relates to a conserved quantity of the dynamics. The mapping between the boundary conditions and these charges being thus equivalent to the mapping between boundary conditions and constants of motion. These considerations emphasize the fact that the deep reason behind the scaling law characterizing the 1d potential free MFG that we consider in this paper is actually their integrable character. A more in depth discussion of this question will appear in a subsequent publication.