Conservation laws of mean field games equations

Roman Kozlov

Department of Business and Management Science, Norwegian School of Economics

Helleveien 30, 5045, Bergen, Norway;

E-mail address<EMAIL_ADDRESS>
Abstract

> Mean field games equations are examined for conservation laws. The system of
> mean field games equations consists of two partial differential equations:
> the Hamilton-Jacobi-Bellman equation for the value function and the forward
> Kolmogorov equation for the probability density. For separable Hamiltonians,
> this system has a variational structure, i.e., the equations of the system
> are Euler-Lagrange equations for some Lagrangian functions. Therefore, one
> can use the Noether theorem to derive the conservation laws using
> variational and divergence symmetries. In order to find such symmetries, we
> find symmetries of the PDE system and select variational and divergence
> ones. The paper considers separable, state-independent Hamiltonians in one-
> dimensional state space. It examines the most general form of the mean field
> games system for symmetries and conservation laws and identifies particular
> cases of the system which lead to additional symmetries and conservation
> laws.

Keywords:

mean field games, conservation law, symmetry, Noether theorem

##  1 Introduction

The Partial Differential Equations (PDEs) approach to Mean Field Games (MFG)
theory is a rapidly developing field of research that has its roots in the
seminal work of J.-M. Lasry and P.-L. Lions in 2006 [1, 2, 3]. Another
approach to the development of MFG theory was suggested by M. Huang, R.
Malhamé, and P. Caines [4, 5]. The theory provides a mathematical framework
for understanding the behaviour of large populations of interacting agents,
where the entire population influences each agent’s behaviour. A continuum of
agents represents such a population [6]. Every again knows the distribution of
the other agents and uses this to implement its strategy.

The mean field game theory has found a broad range of applications in various
fields [7]. In economics, MFG theory has been used to study market competition
[8, 9, 10, 11], price formation [12, 13], and financial markets [14, 15, 16,
17, 18]. In particular, MFG theory has been applied to understand the
behaviour of large populations of investors and traders in financial markets
[19]. MFG models have also been used to study traffic flow and congestion in
transportation systems, as well as to model the behaviour of crowds in social
dynamics [20, 21, 22, 23]. MFG theory has been used to understand the dynamics
of learning [24]. Additionally, MFG theory has been applied to the study of
energy systems, such as electricity markets and renewable energy systems [25,
26, 27, 28, 29], and to the optimization of resource allocation in engineering
applications [30, 31]. Mathematical aspects in different frameworks of the MFG
theory are given in several review papers [32, 33] and books [34, 35, 36, 37].

Conservation laws play an essential role in differential equations. They
reflect the intrinsic properties of the equations. They also help to ensure
the consistency and stability of the system. In the present paper, we analyse
the system of mean field games equations for conservation laws. This system
consists of two PDEs: the Hamilton-Jacobi-Bellman (HJB) equation for the value
function of a typical agent and the forward Kolmogorov equation (also known as
the Fokker-Plank equation) for the evolution of the density of the agent
population. It should be noted that the first equation is backward in time,
and the latter is forward in time. There are several approaches to derive
conservation laws [38]. In the paper, we exploit the variational structure of
the MFG equations. For separable Hamiltonians, this system has a variational
structure, i.e., the equations of the system are Euler-Lagrange equations for
some Lagrangian functions. Therefore, one can use symmetries and the Noether
theorem [39] to derive the conservation laws.

Lie point symmetries provide efficient tools to analyse nonlinear differential
equations [40, 41, 42]. The symmetries are related to essential properties of
the models represented by the differential equations:

  * •

Symmetries generate transformations which transfer solutions of the
differential equations into other solutions. In particular, they can be used
to find new solutions from the known ones.

  * •

Symmetries of PDEs can be used to find particular solutions. These solutions
are called invariant solutions.

  * •

Symmetries of variational differential equations can be used for Noether’s
theorem, which provides conservation laws.

For the application of the Noether theorem, one needs variational formulation
of the considered differential equations and variational or divergence
symmetries. Such symmetries are also symmetries of the considered differential
equations. Therefore, the standard approach is to find symmetries of the
differential equations, select those that are also variational or divergence
symmetries, and apply the Noether theorem to obtain conservation laws. These
steps are implemented in the paper. In particular, we start with the most
general case of the differential equations, for which there are a few
symmetries, and continue with identifying particular cases, for which there
are additional symmetries. After that, we select symmetries suitable for the
Noether theorem and provide the conservation law.

The paper is organized as follows. In the next Section we review the mean
field games equations and specify the cases with the variational formulation.
Section 3 describes Lie point symmetries and the Noether theorem. This section
is written for the particular system which is analysed in Sections 4 and 5,
where we present the symmetry analysis and the conservation laws. Section 6
discusses the interpretation of the obtained conservation laws. Final Section
7 provides concluding remarks. Some technical details are extracted into the
Appendices.

##  2 Variational MFG equations

The system of mean field games equations consists of two coupled PDEs: the
Hamilton-Jacobi-Bellman equation and the forward Kolmogorov equation [37]:

$\displaystyle-u_{t}-\varepsilon\Delta u+H(x,\nabla u,m)=0\qquad\mbox{in}\
[0,T]\times\mathbf{R}^{d},$ (2.1a) $\displaystyle m_{t}-\varepsilon\Delta
m-\mbox{div}(mH_{\nabla u}(x,\nabla u,m))=0\qquad\mbox{in}\
[0,T]\times\mathbf{R}^{d}.$ (2.1b)

The scalar functions $u(t,x)$ and $m(t,x)$ stand for the value function and
the probability density. They depend on time $t\in[0,T]$ and spatial
independent variable $x\in\mathbf{R}^{d}$. $H(x,\nabla u,m)$ is the
Hamiltonian associated to the game and $H_{\nabla
u}=(H_{u_{x_{1}}},\cdots,H_{u_{x_{d}}})$. Parameter $\varepsilon\geq 0$
characterizes the noise. It should be noted that the HJB equation (2.1a) is
backward in time and the Kolmogorov equation (2.1b) is forward in time. Thus,
the initial-terminal conditions are given as

$m(0,x)=m_{0}(x),\qquad u(T,x)=G(x,m(T,x))\qquad\mbox{in}\ \mathbf{R}^{d}.$
(2.2)

Since $m(t,x)\geq 0$ is the probability density, we have

$\int mdx=1.$ (2.3)

The solution $(u,m)$ of the coupled system (2.1) represents a Nash equilibria
of the mean field game. We refer to [37] for detailed assumptions one needs
for the existence and uniqueness of the solutions. A sketch of the derivation
of the system is given in Appendix A.

###  2.1 Separable Hamiltonians

A particular case of the system (2.1)

$\displaystyle-u_{t}-\varepsilon\Delta u+H(x,\nabla u)-f(x,m)=0,$ (2.4a)
$\displaystyle m_{t}-\varepsilon\Delta m-\mbox{div}(mH_{\nabla u}(x,\nabla
u))=0$ (2.4b)

corresponds to separate Hamiltonians $H(x,\nabla u,m)=H(x,\nabla u)-f(x,m)$.
This system has a variational formulation. It is provided by the Lagrangian
function

$L=-mu_{t}+\varepsilon\nabla m\cdot\nabla u+mH(x,\nabla u)-F(x,m),$ (2.5)

where

$F(x,m)=\int f(x,m)dm,$ (2.6)

i.e., equations (2.4) are Euler-Lagrange equations

$\displaystyle{\delta L\over\delta m}={\partial L\over\partial
m}-D_{t}{\partial L\over\partial m_{t}}-\sum_{i=1}^{d}D_{x_{i}}{\partial
L\over\partial m_{x_{i}}}=0,$ (2.7a) $\displaystyle{\delta L\over\delta
u}={\partial L\over\partial u}-D_{t}{\partial L\over\partial
u_{t}}-\sum_{i=1}^{d}D_{x_{i}}{\partial L\over\partial u_{x_{i}}}=0$ (2.7b)

for Lagrangian (2.5). Here $D_{t}$ and $D_{x_{i}}$ are total differentiation
operators for time $t$ and spatial variables ${x_{i}}$, respectively.

###  2.2 One-dimensional state space with state-independent Hamiltonians.

In this paper, we examine the simplest case of system (2.4): the case of one-
dimensional state space $d=1$ with state-independent functions $H=H(u_{x})$
and $f=f(m)$. The equations take the form

$\displaystyle F_{1}=-u_{t}-\varepsilon u_{xx}+H(u_{x})-f(m)=0\qquad\mbox{in}\
[0,T]\times\mathbf{R},$ (2.8a) $\displaystyle F_{2}=m_{t}-\varepsilon
m_{xx}-m_{x}H^{\prime}(u_{x})-mH^{\prime\prime}(u_{x})u_{xx}=0\qquad\mbox{in}\
[0,T]\times\mathbf{R}.$ (2.8b)

We consider second-order MFG systems restricting to $\varepsilon>0$.
Assumptions

$H^{\prime\prime}(u_{x}){\not\equiv}0,\qquad f^{\prime}(m){\not\equiv}0;$
(2.9)

guarantee that the equations are genuinely coupled, i.e., none of these
equations can be considered independently from the other.

Equations (2.8) are Euler-Lagrange equations for the Lagrangian function

$L=-mu_{t}+\varepsilon m_{x}u_{x}+mH(u_{x})-F(m)$ (2.10)

with

$F(m)=\int f(m)dm.$ (2.11)

In Section 4, we analyse the system of equations (2.8) for admitted symmetries
and in Section 5 we derive the conservation laws of these equations.

###  2.3 Other variational cases

Between the most general variational system (2.4) and the simplest variational
system (2.8), two particular MFG systems might be worthy of a separate study.
They correspond to

  * •

Multi-dimensional state space $x\in\mathbf{R}^{d}$ with state-independent
Hamiltonians.

For state-independent functions $H=H(\nabla u)$ and $f=f(m)$, we get the
particular case of the system (2.4)

$\displaystyle-u_{t}-\varepsilon\Delta u+H(\nabla u)-f(m)=0\qquad\mbox{in}\
[0,T]\times\mathbf{R}^{d},$ (2.12a) $\displaystyle m_{t}-\varepsilon\Delta
m-\mbox{div}(mH_{\nabla u}(\nabla u))=0\qquad\mbox{in}\
[0,T]\times\mathbf{R}^{d},$ (2.12b)

provided by Lagrangian

$L=-mu_{t}+\varepsilon\nabla m\cdot\nabla u+mH(\nabla u)-F(m)$

with $F(m)$ given by (2.11).

  * •

One-dimensional state space

If the state space is one-dimensional, system (2.4) gets simplified as

$\displaystyle-u_{t}-\varepsilon u_{xx}+H(x,u_{x})-f(x,m)=0\qquad\mbox{in}\
[0,T]\times\mathbf{R},$ (2.13a) $\displaystyle m_{t}-\varepsilon
m_{xx}-(mH_{u_{x}}(x,u_{x}))_{x}=0\qquad\mbox{in}\ [0,T]\times\mathbf{R}.$
(2.13b)

It’s Lagrangian function is

$L=-mu_{t}+\varepsilon m_{x}u_{x}+mH(x,u_{x})-F(x,m).$

##  3 Lie point symmetries and the Noether theorem

This section reviews Lie point symmetries and the Noether theorem [40, 41,
42]. For convenience, the discussion is adopted to the system of PDEs (2.8)
and first-order Lagrangian function (2.10), which are analysed in this paper.

Lie point symmetries of system (2.8) are given by generators of the form

$X=\xi^{t}{\partial\over\partial t}+\xi^{x}{\partial\over\partial
x}+\eta^{u}{\partial\over\partial u}+\eta^{m}{\partial\over\partial m}$ (3.1)

with generator coefficients

$\xi^{t}=\xi^{t}(t,x,m,u),\qquad\xi^{x}=\xi^{x}(t,x,m,u),\qquad\eta^{u}=\eta^{u}(t,x,m,u),\qquad\eta^{m}=\eta^{m}(t,x,m,u).$

For application to differential equations a generator should be prolonged to
all variables involved in the considered differential equations. In our case,
we get the prolonged generator

$X=\xi^{t}{\partial\over\partial t}+\xi^{x}{\partial\over\partial
x}+\eta^{u}{\partial\over\partial u}+\eta^{m}{\partial\over\partial m}\\\
+\zeta_{t}^{u}{\partial\over\partial
u_{t}}+\zeta_{t}^{m}{\partial\over\partial
m_{t}}+\zeta_{x}^{u}{\partial\over\partial
u_{x}}+\zeta_{x}^{m}{\partial\over\partial
m_{x}}+\zeta_{xx}^{u}{\partial\over\partial
u_{xx}}+\zeta_{xx}^{m}{\partial\over\partial m_{xx}}.$ (3.2)

The coefficients for the derivatives are computed according to the standard
prolongation formulas [40, 41, 42]

$\zeta_{t}^{u}=D_{t}(\eta^{u})-u_{t}D_{t}(\xi^{t})-u_{x}D_{t}(\xi^{x}),\qquad\zeta_{t}^{m}=D_{t}(\eta^{m})-m_{t}D_{t}(\xi^{t})-m_{x}D_{t}(\xi^{x}),$
$\zeta_{x}^{u}=D_{x}(\eta^{u})-u_{t}D_{x}(\xi^{t})-u_{x}D_{x}(\xi^{x}),\qquad\zeta_{x}^{m}=D_{x}(\eta^{m})-m_{t}D_{x}(\xi^{t})-m_{x}D_{x}(\xi^{x}),$
$\zeta_{xx}^{u}=D_{x}(\zeta_{x}^{u})-u_{tx}D_{x}(\xi^{t})-u_{xx}D_{x}(\xi^{x}),\qquad\zeta_{xx}^{m}=D_{x}(\zeta_{x}^{u})-m_{tx}D_{x}(\xi^{t})-m_{xx}D_{x}(\xi^{x}),$

where $D_{t}$ and $D_{x}$ are total differentiation operators for independent
variables $t$ and $x$, respectively, i.e.

$D_{t}={\partial\over\partial t}+u_{t}{\partial\over\partial
u}+m_{t}{\partial\over\partial m}+\cdots,\qquad D_{x}={\partial\over\partial
x}+u_{x}{\partial\over\partial u}+m_{x}{\partial\over\partial m}+\cdots$

To find Lie point symmetries admitted by differential equations, one applies
the infinitesimal invariance criterion [40, 41, 42]. It requires that the
action of the prolonged generator on the considered differential equations
provides identities on the solutions of these equations. For system (2.8), the
infinitesimal invariance criterion takes the form

$\displaystyle E_{1}=\left.X(F_{1})\right|_{F_{1}=0,F_{2}=0}=0,$ (3.3a)
$\displaystyle E_{2}=\left.X(F_{2})\right|_{F_{1}=0,F_{2}=0}=0.$ (3.3b)

These equations are PDEs for the coefficients of the generator (3.1). Solving
them, one finds the admitted generators.

Lie point symmetries of the differential equations are needed for the Noether
theorem.

###### Theorem 3.1

(E. Noether [39]) Let Lagrangian function (2.10) satisfy the equation

$XL+L(D_{t}(\xi^{t})+D_{x}(\xi^{x}))=0$ (3.4)

for generator $X$ given by (3.1), then the generator $X$ is a symmetry of the
system of Euler-Lagrange equations (2.8) and there holds the conservation law

$D_{t}[T^{t}]+D_{x}[T^{x}]=0,$ (3.5)

where

$\displaystyle T^{t}=\xi^{t}L+(\eta^{u}-\xi^{t}u_{t}-\xi^{x}u_{x}){\partial
L\over\partial u_{t}}+(\eta^{m}-\xi^{t}m_{t}-\xi^{x}m_{x}){\partial
L\over\partial m_{t}},$ (3.6a) $\displaystyle
T^{x}=\xi^{x}L+(\eta^{u}-\xi^{t}u_{t}-\xi^{x}u_{x}){\partial L\over\partial
u_{x}}+(\eta^{m}-\xi^{t}m_{t}-\xi^{x}m_{x}){\partial L\over\partial m_{x}}.$
(3.6b)

There is a natural extension of the Noether theorem.

###### Remark 3.2

([43]) If instead of relation (3.4), there holds

$XL+L(D_{t}(\xi^{t})+D_{x}(\xi^{x}))=D_{t}(V^{t})+D_{x}(V^{x})$ (3.7)

with some functions $V^{t}(t,x,u,m)$ and $V^{x}(t,x,u,m)$, then the generator
$X$ is a symmetry of the system of Euler-Lagrange equations (2.8) and there
holds the conservation law

$D_{t}[T^{t}-V^{t}]+D_{x}[T^{x}-V^{x}]=0$ (3.8)

with $T^{t}$ and $T^{x}$ given by (3.6).

Symmetries of the variational equations which satisfy (3.4) are called
variational symmetries and symmetries with property (3.7) are called
divergence symmetries.

##  4 Symmetry analysis of system (2.8)

This section provides analysis of system (2.8),(2.9) for admitted symmetries.
Lie point symmetries of system (2.8) are provided by the infinitesimal
invariance criterion (3.3). In detail, these equations take the form

$E_{1}=\left[-\zeta_{t}^{u}-\varepsilon\zeta_{xx}^{u}+H^{\prime}(u_{x})\zeta_{x}^{u}-f^{\prime}(m)\eta^{m}\right]_{F_{1}=0,F_{2}=0}=0,$
$E_{2}=\left[\zeta_{t}^{m}-\varepsilon\zeta_{xx}^{m}-\zeta_{x}^{m}H^{\prime}(u_{x})-m_{x}H^{\prime\prime}(u_{x})\zeta_{x}^{u}\right.\\\
\left.-\eta^{m}H^{\prime\prime}(u_{x})u_{xx}-mH^{\prime\prime\prime}(u_{x})\zeta_{x}^{u}u_{xx}-mH^{\prime\prime}(u_{x})\zeta_{xx}^{u}\right]_{F_{1}=0,F_{2}=0}=0.$

These equations are used to find symmetries admitted for arbitrary $H(u_{x})$
and $f(m)$. They are also used to identify special cases of $H(u_{x})$ and
$f(m)$ which lead to extensions of the symmetries admitted for arbitrary
$H(u_{x})$ and $f(m)$. The preliminary analysis of these equations aimed at
identifying the special cases of the Hamiltonian function is given in Appendix
B. It results in four special cases of $H(u_{x})$ which should be analysed in
addition to the general case of arbitrary $H(u_{x})$. Thus, we have the
following cases to examine

  * •

General case: arbitrary $H(u_{x})$

  * •

Special cases:

    * –

power

$H(u_{x})=h(u_{x}+q)^{p}+h_{2}u_{x}^{2}+h_{1}u_{x}+h_{0},\qquad h\neq 0,\qquad
p\neq 0,1,2,3;$ (4.1)

    * –

exponential

$H(u_{x})=he^{ku_{x}}+h_{2}u_{x}^{2}+h_{1}u_{x}+h_{0},\qquad h\neq 0,\qquad
k\neq 0;$ (4.2)

    * –

cubic (particular power)

$H(u_{x})=hu_{x}^{3}+h_{2}u_{x}^{2}+h_{1}u_{x}+h_{0},\qquad h\neq 0;$ (4.3)

    * –

quadratic (particular power)

$H(u_{x})=hu_{x}^{2}+h_{1}u_{x}+h_{0},\qquad h\neq 0.$ (4.4)

Here $h$, $h_{2}$, $h_{1}$, $h_{0}$, $p$, $q$ and $k$ are constants.

Before we consider these cases, it is appropriate to simplify them by
transformations which do not change the form of equations (2.8). Such
transformations are called equivalence transformations. We apply the following
steps.

  1. 1.

Transformations of the form

$u\rightarrow\bar{u}=u+Ax,\qquad A=\mbox{const}$ (4.5)

are used to remove parameter $q$ in Hamiltonian (4.1), to scale coefficient
$h$ in Hamiltonian (4.2) and to remove the term $h_{2}u_{x}^{2}$ in
Hamiltonian (4.3).

  2. 2.

Transformations

$u\rightarrow\bar{u}=u+Bt,\qquad B=\mbox{const}$ (4.6)

are used to remove constant terms $h_{0}$ in all Hamiltonians (4.1)-(4.4) as
well as in functions $f(m)$ (if $f(m)$ contains a constant term).

  3. 3.

The transformation

$x\rightarrow\bar{x}=x+h_{1}t$ (4.7)

removes the terms $h_{1}u_{x}$ in all Hamiltonians (4.1)-(4.4).

  4. 4.

Scaling of independent and dependent variables

$t\rightarrow\bar{t}=C_{1}t,\qquad x\rightarrow\bar{x}=C_{2}x,\qquad
u\rightarrow\bar{u}=C_{3}u,\qquad m\rightarrow\bar{m}=C_{4}m,$ (4.8)

where $C_{1}$, $C_{2}$, $C_{3}$ and $C_{4}$ are constants such that
$C_{1}C_{2}C_{3}C_{4}\neq 0$, are used to make a suitable choice of the
leading term coefficient $h$ in Hamiltonians (4.1), (4.3) and (4.4). Scaling
transformations can be chosen so that they satisfy

$C_{2}C_{4}=1.$

In addition to the form of equations (2.8) such scaling transformations
preserve normalization (2.3). Function $f(m)$ is not specified, and it can be
changed by the scaling transformations.

Implementing these steps, we arrive at the following simplified Hamiltonians
for the special cases:

  * •

power

$H(u_{x})={1\over p}u_{x}^{p}+h_{2}u_{x}^{2},\qquad p\neq 1,2,3;$ (4.9)

  * •

exponential

$H(u_{x})={1\over k}e^{ku_{x}}+h_{2}u_{x}^{2};$ (4.10)

  * •

cubic

$H(u_{x})={1\over 3}u_{x}^{3};$ (4.11)

  * •

quadratic

$H(u_{x})={1\over 2}u_{x}^{2}.$ (4.12)

Here $h_{2}$ is an arbitrary constant.

Detailed analysis of determining equations (3.3) for the general and special
cases of the Hamiltonian provides the following result.

###### Theorem 4.1

Up to equivalence transformations (4.5)–(4.8), system (2.8),(2.9). has the
following cases of admitted symmetries:

For arbitrary $H(u_{x})$ and arbitrary $f(m)$, the system admits the three
symmetries

$X_{1}={\partial\over\partial t},\qquad X_{2}={\partial\over\partial x},\qquad
X_{3}={\partial\over\partial u}.$ (4.13)

If $f(m)=\alpha\ln m$, there exists the additional symmetry

$X_{f}=\alpha t{\partial\over\partial u}-m{\partial\over\partial m}.$ (4.14)

There are particular cases of $H(u_{x})$ which lead to the extensions of the
admitted symmetry group.

a) For power $H(u_{x})={1\over p}u_{x}^{p}$ (Hamiltonian (4.9) with $h_{2}=0$)
and power $f(m)=\alpha m^{\gamma}$, there is the additional generator

$X^{a}_{4}=2(p-1)t{\partial\over\partial t}+(p-1)x{\partial\over\partial
x}+(p-2)u{\partial\over\partial u}-{p\over\gamma}m{\partial\over\partial m}.$
(4.15)

b) For exponential $H(u_{x})={1\over k}e^{ku_{x}}$ (Hamiltonian (4.10) with
$h_{2}=0$) and power $f(m)=\alpha m^{\gamma}$, there exists the additional
symmetry

$X^{b}_{4}=2t{\partial\over\partial t}+x{\partial\over\partial
x}+\left(u-{x\over k}\right){\partial\over\partial
u}-{1\over\gamma}m{\partial\over\partial m}.$ (4.16)

c) For quadratic $H(u_{x})={1\over 2}u_{x}^{2}$ and arbitrary $f(m)$ there is
the generator

$X^{c}=t{\partial\over\partial x}-x{\partial\over\partial u}.$ (4.17)

If $f(m)=\alpha m^{2}$, then there is the particular case of symmetry (4.15),
corresponding the $p=2$, namely

$X^{c}_{4}=\left.X^{a}_{4}\right|_{p=2,\gamma=2}=2t{\partial\over\partial
t}+x{\partial\over\partial x}-m{\partial\over\partial m},$ (4.18)

and one more symmetry

$X^{c}_{5}=t^{2}{\partial\over\partial t}+tx{\partial\over\partial
x}+\left(-{x^{2}\over 2}+\varepsilon t\right){\partial\over\partial
u}-tm{\partial\over\partial m}.$ (4.19)

Proof. The system of determining equations (3.3) is an overdetermined system
of PDEs for the coefficients of generator $X$, given by (3.1). These equations
are analysed for the case of arbitrary Hamiltonian $H(u_{x})$ and for four
special Hamiltonians (4.9)-(4.12). For the considered cases of the Hamiltonian
function the determining equations split (first for the derivatives of $u$ and
$m$, then for the variables) into equations that specify the symmetry
coefficients.

Particular cases of $f(m)$, namely $f(m)=\alpha\ln m$, $f(m)=\alpha
m^{\gamma}$ and $f(m)=\alpha m^{2}$, are identified from the detailed analysis
of the determining equations as cases that provide additional symmetries.
$\Box$

The admitted symmetries can be divided into the basis symmetry groups, which
is admitted for all functions $H(u_{x})$ and $f(m)$, and the additional
symmetries, which exist for the particular functions $H(u_{x})$ and $f(m)$.
The basis symmetry group consists of three generators
$\\{X_{1},X_{2},X_{3}\\}$, given in (4.13). For convenience, we present the
extensions of the basic symmetry group in Table 1.

$\begin{array}[]{|c|c|c|c|c|}\hline\cr&\mbox{arbitrary $f(m)$}&f(m)=\alpha\ln
m&f(m)=\alpha m^{\gamma},\ \gamma\neq 2&f(m)=\alpha m^{2}\\\
\hline\cr\mbox{arbitrary $H(u_{x})$}&-&X_{f}&-&-\\\ \hline\cr H(u_{x})={1\over
p}u_{x}^{p}&-&X_{f}&X^{a}_{4}&\left.X^{a}_{4}\right|_{\gamma=2}\\\ \hline\cr
H(u_{x})={1\over
k}e^{ku_{x}}&-&X_{f}&X_{4}^{b}&\left.X^{b}_{4}\right|_{\gamma=2}\\\ \hline\cr
H(u_{x})={1\over
2}u_{x}^{2}&X^{c}&X_{f},X^{c}&X^{c},X_{4}&X^{c},X^{c}_{4},X^{c}_{5}\\\
\hline\cr\end{array}$

Table 1. Extensions of the basic symmetry group with generators
$\\{X_{1},X_{2},X_{3}\\}$, which is admitted for arbitrary $H(u_{x})$ and
arbitrary $f(m)$. Classification is performed up to equivalence
transformations (4.5)–(4.8).

###### Remark 4.2

In Theorem 4.1, the case of cubic Hamiltonian (4.11) stands as a particular
case of power Hamiltonian (4.9). The cubic Hamiltonian was identified
separately from the power Hamiltonian in the preliminary analysis of the
determining equations, given in Appendix B. However, it has the same symmetry
properties as the general power Hamiltonian case.

###### Remark 4.3

It is possible to employ scaling transformation to get rid of the parameter
$\alpha$ in the particular cases of $f(m)$ without changing the form of
equations (2.8). However, for some parameter values, such scaling
transformations do not keep normalization relation (2.3). For this reason,
parameter $\alpha$ is kept.

##  5 Conservation laws of system (2.8)

The cases of $H$ and $f$, which should be considered, are specified by Theorem
4.1. It remains to select variational and divergence symmetries and apply the
Noether theorem to derive the conservation laws.

###  5.1 General case: arbitrary $H(u_{x})$

For arbitrary Hamiltonian, two cases of the function $f(m)$ exist to consider.

1) Arbitrary $f(m)$

In this case, we consider

$F(m)=\int f(m)dm$ (5.1)

corresponding to arbitrary $f(m)$.

Equations (2.8) admit symmetries $X_{1}$, $X_{2}$ and $X_{3}$, given in
(4.13). These three symmetries are variational.

  * •

The symmetry $X_{1}$ provides the conservation law

$D_{t}\left[\varepsilon
m_{x}u_{x}+mH(u_{x})-F(m)\right]-D_{x}\left[\varepsilon(u_{t}m_{x}+m_{t}u_{x})+mu_{t}H^{\prime}(u_{x})\right]=0.$
(5.2)

  * •

The symmetry $X_{2}$ leads to the conservation law

$D_{t}\left[mu_{x}\right]-D_{x}\left[mu_{t}+\varepsilon
m_{x}u_{x}+m(u_{x}H^{\prime}(u_{x})-H(u_{x}))+F(m)\right]=0.$ (5.3)

  * •

The symmetry $X_{3}$ gives the conservation law

$-D_{t}\left[m\right]+D_{x}\left[\varepsilon
m_{x}+mH^{\prime}(u_{x})\right]=0.$ (5.4)

Note that this conservation law is the Kolmogorov equation (2.8b) itself. It
reflects the conservation of the probability density $m(t,x)$.

2) Logarithmic $f(m)=\alpha\ln m$

In this case

$F(m)=\alpha(m\ln m-m).$ (5.5)

The generator $X_{f}$ (see (4.14)) is neither variational nor divergence
symmetry. Thus, it does not provide a conservation law.

###  5.2 Special cases of $H(u_{x})$

There are three exceptional cases of Hamiltonian: (general) power, exponential
and quadratic.

a) Power $H(u_{x})={1\over p}u_{x}^{p}$ and power $f(m)=\alpha m^{\gamma}$

In this case

$F(m)={\alpha\over\gamma+1}m^{\gamma+1}.$ (5.6)

Symmetry $X^{a}_{4}$, given by (4.15), is variational, provided that

$\gamma=\gamma_{*}={p\over 2p-3}.$ (5.7)

Thus we consider the symmetry

${Y}^{a}_{4}=\left.X_{4}\right|_{\gamma={\gamma_{*}}}=2(p-1)t{\partial\over\partial
t}+(p-1)x{\partial\over\partial x}+(p-2)u{\partial\over\partial
u}-(2p-3)m{\partial\over\partial m}.$ (5.8)

It provides the conservation laws

$D_{t}\left[2(p-1)t\left(\varepsilon m_{x}u_{x}+{1\over
p}mu_{x}^{p}-{\alpha\over\gamma_{*}+1}m^{\gamma_{*}+1}\right)+(p-1)xmu_{x}-(p-2)mu\right]\\\
-D_{x}\left[2(p-1)t\left(\varepsilon(u_{t}m_{x}+m_{t}u_{x})+mu_{t}u_{x}^{p-1}\right)\phantom{7\over
7}\right.\\\ +(p-1)x\left(mu_{t}+\varepsilon m_{x}u_{x}+{p-1\over
p}mu_{x}^{p}+{\alpha\over\gamma_{*}+1}m^{\gamma_{*}+1}\right)\\\
\left.\phantom{51\over 51}-(p-2)u(\varepsilon
m_{x}+mu_{x}^{p-1})+(2p-3)\varepsilon mu_{x}\right]=0.$ (5.9)

b) Exponential $H(u_{x})={1\over k}e^{ku_{x}}$ and power $f(m)=\alpha
m^{\gamma}$

As in the previous case, $F(m)$ is given by (5.6). The symmetry $X^{b}_{4}$,
given in (4.16), is a divergence symmetry for

$\gamma={1\over 2}$ (5.10)

with

$V^{t}_{4}=0,\qquad V^{x}_{4}=-{\varepsilon\over k}m.$ (5.11)

i.e., for this symmetry, there holds

$XL+L(D_{t}(\xi^{t})+D_{x}(\xi^{x}))=D_{x}\left(-{\varepsilon\over
k}m\right).$

Therefore, we get the symmetry $X^{b}_{4}$ specified as

${Y}^{b}_{4}=\left.X^{b}_{4}\right|_{\gamma={1\over
2}}=2t{\partial\over\partial t}+x{\partial\over\partial x}+\left(u-{x\over
k}\right){\partial\over\partial u}-{2}m{\partial\over\partial m}.$ (5.12)

Using this symmetry, we obtain the conservation law

$D_{t}\left[2t\left(\varepsilon m_{x}u_{x}+{1\over k}me^{ku_{x}}-{2\alpha\over
3}m^{3/2}\right)+xm\left(u_{x}+{1\over k}\right)-mu\right]\\\
-D_{x}\left[2t\left(\varepsilon(u_{t}m_{x}+m_{t}u_{x})+mu_{t}e^{ku_{x}}\right)+x\left(mu_{t}+\varepsilon
m_{x}u_{x}+mu_{x}e^{ku_{x}}+{\varepsilon\over k}m_{x}+{2\alpha\over
3}m^{3/2}\right)\right.\\\ \left.-u(\varepsilon
m_{x}+me^{ku_{x}})+2\varepsilon mu_{x}-{\varepsilon\over k}m\right]=0.$ (5.13)

c) Quadratic $H(u_{x})={1\over 2}u_{x}^{2}$

For the quadratic Hamiltonian, there are two cases of $f(m)$.

1) Arbitrary $f(m)$

For arbitrary $f(m)$, we consider the corresponding function $F(m)$ as given
in (5.1). The symmetry $X^{c}$, given in (4.17), is a divergence symmetry with

$V^{t}=0,\qquad V^{x}=-{\varepsilon}m.$ (5.14)

This symmetry leads to the conservation law

$D_{t}\left[m(tu_{x}+x)\right]\\\ -D_{x}\left[t\left(mu_{t}+\varepsilon
m_{x}u_{x}+{1\over 2}mu_{x}^{2}+F(m)\right)+x(\varepsilon
m_{x}+mu_{x})-{\varepsilon}m\right]=0.$ (5.15)

2) Quadratic $f(m)=\alpha m^{2}$

Now

$F(m)={\alpha\over 3}m^{3}.$ (5.16)

It is a particular case of (5.6). The symmetry

${Y}^{c}_{4}=\left.{Y}^{a}_{4}\right|_{p=2}=2t{\partial\over\partial
t}+x{\partial\over\partial x}-m{\partial\over\partial m},$

i.e., symmetry (5.8) with $p=2$, is variational. It provides a particular case
of conservation law (5.9) corresponding to $\gamma_{*}=p=2$, namely

$D_{t}\left[2t\left(\varepsilon m_{x}u_{x}+{1\over 2}mu_{x}^{2}-{\alpha\over
3}m^{3}\right)+xmu_{x}\right]\\\
-D_{x}\left[2t\left(\varepsilon(u_{t}m_{x}+m_{t}u_{x})+mu_{t}u_{x}\right)+x\left(mu_{t}+\varepsilon
m_{x}u_{x}+{1\over 2}mu_{x}^{2}+{\alpha\over 3}m^{3}\right)+\varepsilon
mu_{x}\right]=0.$ (5.17)

The symmetry $X_{5}^{c}$. given in (4.19), is a divergence one with

$V^{t}=0,\qquad V^{x}=-{\varepsilon}xm.$ (5.18)

It gives the conservation law

$D_{t}\left[t^{2}\left(\varepsilon m_{x}u_{x}+{1\over
2}mu_{x}^{2}-{\alpha\over 3}m^{3}\right)+{tx}mu_{x}+{x^{2}\over
2}m-\varepsilon tm\right]\\\
-D_{x}\left[t^{2}\left(\varepsilon(u_{t}m_{x}+m_{t}u_{x})+mu_{t}u_{x}\right)+tx\left(mu_{t}+\varepsilon
m_{x}u_{x}+{1\over 2}mu_{x}^{2}+{\alpha\over 3}m^{3}\right)\right.\\\
\left.+{x^{2}\over 2}(\varepsilon
m_{x}+mu_{x})-\varepsilon^{2}tm_{x}-{\varepsilon}xm\right]=0.$ (5.19)

The obtained results can be summarized as a theorem.

###### Theorem 5.1

Up to equivalence transformations (4.5)–(4.8), system (2.8),(2.9) has the
following cases of conservation laws:

For arbitrary $H(u_{x})$ and arbitrary $f(m)$, the system has three
conservation laws (5.2), (5.3) and (5.4).

There are particular cases of $H(u_{x})$ which lead to the additional
conservation laws.

a) For power $H(u_{x})={1\over p}u_{x}^{p}$ and power $f(m)=\alpha m^{\gamma}$
with $\gamma=\gamma_{*}={p\over 2p-3}$, there is the additional conservation
law (5.9).

b) For exponential $H(u_{x})={1\over k}e^{ku_{x}}$ and $f(m)=\alpha\sqrt{m}$,
there exists the additional conservation law (5.13).

c) For quadratic $H(u_{x})={1\over 2}u_{x}^{2}$ and arbitrary $f(m)$ there is
conservation law (5.15). If $f(m)=\alpha m^{2}$, then there also hold
conservation laws (5.17) and (5.19).

Thus, three conservation laws exist for all $H(u_{x})$ and $f(m)$. The
symmetries which provide additional conservation laws are given in Table 2.

$\begin{array}[]{|c|c|c|c|}\hline\cr&\mbox{arbitrary $f(m)$}&f(m)=\alpha
m^{\gamma},\ \gamma\neq 2&f(m)=\alpha m^{2}\\\ \hline\cr\mbox{arbitrary
$H(u_{x})$}&-&-&-\\\ \hline\cr H(u_{x})={1\over p}u_{x}^{p},\ p\neq
2&-&{Y}^{a}_{4}&-\\\ \hline\cr H(u_{x})={1\over
k}e^{ku_{x}}&-&{Y}^{b}_{4}&-\\\ \hline\cr H(u_{x})={1\over
2}u_{x}^{2}&X^{c}&X^{c}&X^{c},{Y}^{c}_{4},X^{c}_{5}\\\ \hline\cr\end{array}$

Table 2. Variational and divergence symmetries admitted in addition to
symmetries (4.13), which correspond to arbitrary $H(u_{x})$ and arbitrary
$f(m)$. Classification is performed up to equivalence transformations
(4.5)–(4.8).

##  6 Interpretation of conservation laws

Some conservation laws allow intuitive interpretation. Finite group
transformations

$t\rightarrow\bar{t},\qquad x\rightarrow\bar{x},\qquad
u\rightarrow\bar{u},\qquad m\rightarrow\bar{m}$

for generators (3.1) solve the ODE system [40, 41]

${d\bar{t}\over
da}=\xi^{t}(\bar{t},\bar{x},\bar{u},\bar{m}),\qquad{d\bar{x}\over
da}=\xi^{x}(\bar{t},\bar{x},\bar{u},\bar{m}),\qquad{d\bar{u}\over
da}=\eta^{u}(\bar{t},\bar{x},\bar{u},\bar{m}),\qquad{d\bar{m}\over
da}=\eta^{m}(\bar{t},\bar{x},\bar{u},\bar{m})$

with the initial conditions

$\bar{t}|_{a=0}=t,\qquad\bar{x}|_{a=0}=x,\qquad\bar{u}|_{a=0}=u,\qquad\bar{m}|_{a=0}=m.$

Here $a$ is the group parameter. Its value $a=0$ corresponds to the identity
transformation. We remark the solutions can be presented by the exponentiation
[42]

$t\rightarrow\bar{t}=e^{aX}t,\qquad x\rightarrow\bar{x}=e^{aX}x,\qquad
u\rightarrow\bar{u}=e^{aX}u,\qquad m\rightarrow\bar{m}=e^{aX}m.$

Specifically, we consider a particular case of the mean field games described
in Appendix A. Let the dynamics of the agents be described by the SDE

$dX_{s}=\alpha_{s}ds+\sqrt{2\varepsilon}dB_{s},$ (6.1)

where $\alpha_{s}$ is the agent’s control. Brownian motions $B_{s}$ for
different agents are assumed to be independent. Each agent minimizes

$\mbox{E}\left[\int_{0}^{T}\left({\alpha_{s}^{2}\over
2}+f(m(s,X_{s}))\right)ds+G(X_{T},m(T,X_{T}))\ |\ X_{0}=x\right].$ (6.2)

The running costs have two terms. The first term ${\alpha_{s}^{2}\over 2}$
requires to minimise the state changes. The second term $f(m)$ represents the
interaction of the agents. It can be a willingness to avoid others if $f(m)$
is an increasing function or a willingness to be like others if $f(m)$ is a
decreasing function.

We obtain the Hamiltonian

$H(x,p,m)={p^{2}\over 2}-f(m),$

the MFG system

$\displaystyle-u_{t}-\varepsilon u_{xx}+{u_{x}^{2}\over 2}=f(m),$ (6.3a)
$\displaystyle m_{t}-\varepsilon m_{xx}-(mu_{x})_{x}=0.$ (6.3b)

and the optimal feedback control

${\alpha_{s}}=-u_{x}(t,x).$ (6.4)

System (6.3) has four conservation laws. Conservation law (5.3) corresponds to
the symmetry $X_{2}$, which represents the translation of the state variable
$x$

$\bar{t}=t,\qquad\bar{x}=x+a,\qquad\bar{u}=u,\qquad\bar{m}=m.$

Taking (6.4) into account, it is possible to consider this conservation law as
the conservation of the average control value. In mechanics, translations of
the spatial coordinates lead to the conservation of momenta.

Conservation law (5.4) is provided by translations of $u$,

$\bar{t}=t,\qquad\bar{x}=x,\qquad\bar{u}=u+a,\qquad\bar{m}=m,$

given by the symmetry $X_{3}$. It has the same interpretation for the
quadratic Hamiltonian as for the arbitrary Hamiltonian. It represents the
conservation of the probability measure (2.3). It is analogous to the
conservation of mass in mechanical systems.

Conservation law (5.15) is derived with the help of the Galileo-type symmetry
$X^{c}$. This symmetry generates the group transformation

$\bar{t}=t,\qquad\bar{x}=x+at,\qquad\bar{u}=u-ax-{a^{2}\over
2}t,\qquad\bar{m}=m.$

It can be viewed as the law describing the change of the average state value.
In mechanics, the Galileo transformation corresponds to the law of motion of
the centre of mass.

There seems to be no natural interpretation for conservation law (5.2) as well
as for the additional conservation laws (5.17) and (5.19), which exist for
quadratic $f(m)$. These conservation laws correspond to the time translation

$\bar{t}=t+a,\qquad\bar{x}=x,\qquad\bar{u}=u,\qquad\bar{m}=m,$

the scaling

$\bar{t}=e^{2a}t,\qquad\bar{x}=e^{a}x,\qquad\bar{u}=u,\qquad\bar{m}=e^{-a}m,$

and the projective transformation in $(t,x)$ plane

$\bar{t}={t\over 1-at},\qquad\bar{x}={x\over 1-at}$

supplemented with the transformation of the dependent variables

$\bar{u}=u-{ax^{2}\over 2(1-at)}-\varepsilon\ln(1-at),\qquad\bar{m}=(1-at)m.$

The transformations are generated by symmetries $X_{1}$, $Y^{c}_{4}$ and
$X^{c}_{5}$, respectively. In mechanics, the time translation provides the
conservation of energy. Conservation laws corresponding to the scaling
transformations and the projective transformations are known in mechanics,
kinetic theory gases, hydrodynamics, etc. [44]. They have no specific names.

##  7 Conclusion

The paper considers system (2.8),(2.9). This system is the simplest case of
the mean field games system (2.4), which has a variational formulation due to
the separable Hamiltonian. System (2.8) corresponds to state-independent
Hamiltonians in one-dimensional state space. The paper provides Lie point
symmetry analysis of the system. Then, the variational and divergence
symmetries are selected. They are used for the application of the Noether
theorem to derive the conservation laws. The presented approach provides the
symmetries and conservation laws for the general case of the system
(2.8),(2.9) with arbitrary functions $H(u_{x})$ and $f(m)$ and identifies the
particular cases of $H(u_{x})$ and $f(m)$ which lead to the additional
symmetries and conservation laws. Some conservation laws have natural
interpretations.

Conservation laws of the first-order mean field games systems, which
correspond to parameter $\varepsilon=0$ in (2.1), were considered in [45, 46].
To our knowledge, there have been no attempts to investigate second-order mean
field games equations for conservation laws that makes the results of this
paper original. The study given here can also be performed for more general
MFG systems, such as

  * •

general system (2.4) for multi-dimensional state space with state-dependent
$H$ and $f$,

  * •

system (2.12) for multi-dimensional state space with state-independent $H$ and
$f$,

  * •

system (2.13) for one-dimensional state space with state-dependent $H$ and
$f$.

It is also possible to consider stationary versions of these mean field games
equations.

For example, system (2.4) with arbitrary $H(x,\nabla u)$ and $f(x,m)$ admits
only two symmetries

$X_{1}={\partial\over\partial t},\qquad X_{2}={\partial\over\partial u}.$

The generator $X_{1}$ is a variational symmetry and the generator $X_{2}$ is a
divergence symmetry of the Lagrangian function (2.5). They provide the
conservation laws

$D_{t}\left[\varepsilon\nabla m\cdot\nabla u+mH(x,\nabla
u)-F(x,m)\right]-\mbox{div}\left[\varepsilon(u_{t}\nabla m+m_{t}\nabla
u)+mu_{t}H_{\nabla u}(x,\nabla u)\right]=0$

and

$-D_{t}\left[m\right]+\mbox{div}\left[\varepsilon\nabla m+mH_{\nabla
u}(x,\nabla u)\right]=0,$

respectively. These conservation laws generalize conservation laws (5.2) and
(5.4), which hold for system (2.8) with arbitrary $H(u_{x})$ and $f(m)$.

Conservation laws reflect fundamental properties of the differential
equations, that helps to understand and analyse the models. For this reason,
it is crucial to preserve them in numerical simulation [47]. In particular,
conservation of measure (2.3) is always respected by numerical methods for the
MFG equations [48]. Finally, it is worth mentioning that symmetries of PDEs
can be used to find particular solutions [40, 41, 42]. This has not been
explored for the mean field games equations.

## Acknowledgments

The author is grateful to Prof. Vladimir A. Dorodnitsyn for discussions
concerning methods of Lie group analysis.

## References

  * [1] J.-M. Lasry and P.-L. Lions, Mean field games. I - The stationary case — Jeux à champ moyen. I - Le cas stationnaire, Comptes Rendus Mathematique, 343(9), 619–625, 2006
  * [2] J.-M. Lasry and P.-L. Lions, Mean field games. II - Finite horizon and optimal control — Jeux à champ moyen. II - Horizon fini et contrôle optimal, Comptes Rendus Mathematique, 343(10), 679–684, 2006
  * [3] J.-M. Lasry and P.-L. Lions, Mean field games, Japanese Journal of Mathematics, 2(1), 229–260, 2007
  * [4] M. Huang, R. P. Malhamé and P. E. Caines, Large population stochastic dynamic games: closed-loop Mckean-Vlasov systems and the Nash certainty equivalence principle. Commun. Inf. Syst. 6, 221–252, 2006
  * [5] M. Huang, R. P. Malhamé and P. E. Caines, Large-population cost-coupled LQG problems with nonuniform agents: individual-mass behavior and decentralized $\epsilon$-Nash equilibria, IEEE Trans. Autom. Control, 52, 1560–1571, 2007
  * [6] R. J. Aumann, Markets with a continuum of traders. Econometrica, 32 (1-2), 39–50, 1964
  * [7] O. Guéant, J.-M. Lasry and P.-L. Lions, Mean field games and applications, Lecture Notes in Mathematics 2003, 205–266, 2011.
  * [8] P. Chan and R. Sircar, Bertrand and Cournot Mean Field Games, Applied Mathematics and Optimization, 71(3), 533–569, 2015
  * [9] P. J. Graber and C. Mouzouni, Variational Mean Field Games for Market Competition, Springer INdAM Series, 28, 93–114, 2018
  * [10] R. Carmona and G. Dayanlkll, Mean Field Game Model for an Advertising Competition in a Duopoly, International Game Theory Review, 23(4), 2150024, 2021
  * [11] R. Chenavaz, C., Paraschiv and G. Turinici, Dynamic Pricing of New Products in Competitive Markets: A Mean-Field Game Approach, Dynamic Games and Applications, 11(3), 463–490, 2021
  * [12] D. A. Gomes and J. Saúde, A Mean-Field Game Approach to Price Formation, Dynamic Games and Applications, 11(1), 29–53, 2021
  * [13] D. Gomes, J. Gutierrez and R. Ribeiro, A Random Supply Mean Field Game Price Model, SIAM Journal on Financial Mathematics, 14(1), 188–222, 2023
  * [14] R. Carmona, J.-P. Fouque and L.-H. Sun, Mean field games and systemic risk, Communications in Mathematical Sciences, 13(4), 911–933, 2015
  * [15] R. Carmona, F. Delarue and D. Lacker, Mean Field Games of Timing and Models for Bank Runs, Applied Mathematics and Optimization, 76(1), 217–260, 2017
  * [16] G. Fu, P. Graewe, U. Horst and A. Popier, A mean field game of optimal portfolio liquidation, Mathematics of Operations Research, 46(4), 1250–1281, 2021
  * [17] G. Fu and C. Zhou, Mean field portfolio games, Finance and Stochastics, 27(1), 189–231, 2023
  * [18] G. Fu, Mean field portfolio games with consumption, Mathematics and Financial Economics, 17(1), 79–99, 2023
  * [19] A. Lachapelle, J.-M. Lasry, C.-A. Lehalle and P.-L. Lions, Efficiency of the price formation process in presence of high frequency participants: a mean field game analysis Mathematics and Financial Economics, 10(3), 223–262, 2016
  * [20] A. Lachapelle and M.-T. Wolfram, On a mean field game approach modelling congestion and aversion in pedestrian crowds, Transportation Research Part B: Methodological, 45(10), 1572–1589, 2011
  * [21] M. Burger, M. Di Francesco, P. A. Markowich and M.-T. Wolfram, On a mean field game optimal control approach modelling fast exit scenarios in human crowds, Proceedings of the IEEE Conference on Decision and Control, 6760360, 3128–3133, 2013
  * [22] M. Burger, M. Di Francesco, P. A. Markowich and M.-T. Wolfram, Mean field games with nonlinear mobilities in pedestrian dynamics, Discrete and Continuous Dynamical Systems - Series B, 19(5), 1311–1333, 2014
  * [23] Y. Achdou and J.-M. Lasry, Mean field games for modelling crowd motion, Computational Methods in Applied Sciences, 47, 17–42, 2019
  * [24] M. Burger, A. Lorz and M.-T. Wolfram, On a Boltzmann mean field model for knowledge growth, SIAM Journal on Applied Mathematics, 76(5), 1799–1818, 2016
  * [25] P. Chan and R. Sircar, Fracking, renewables, and mean field games, SIAM Review, 59(3), 588–615, 2017
  * [26] C. Alasseur, I. Ben Taher and A. Matoussi, An Extended Mean Field Game for Storage in Smart Grids, Journal of Optimization Theory and Applications, 184(2), 644–670, 2020
  * [27] B. Djehiche, J. Barreiro-Gomez and H. Tembine, Price Dynamics for Electricity in Smart Grid Via Mean-Field-Type Games, Dynamic Games and Applications, 10(4), 798–818, 2020
  * [28] R. Aïd, R. Dumitrescu and P. Tankov, The entry and exit game in the electricity markets: a mean-field game approach, Journal of Dynamics and Games, 8(4), 331–358, 2021
  * [29] R. Carmona, G. Dayanıklı and M. Laurière, Mean Field Models to Regulate Carbon Emissions in Electricity Production, Dynamic Games and Applications, 12(3), 897–928, 2022
  * [30] Y. Achdou, P.-N. Giraud, J.-M. Lasry and P.-L. Lions, A Long-Term Mathematical Model for Mining Industries, Applied Mathematics and Optimization, 74(3), 579–618, 2016
  * [31] Y. Achdou, C. Bertucci, J.-M. Lasry, P.-L. Lions, A. Rostand and J. A. Scheinkman, A class of short-term models for the oil industry that accounts for speculative oil storage, Finance and Stochastics, 26(3), 631–669, 2022
  * [32] D.A. Gomes and J. Saúde, Mean Field Games Models - A Brief Survey, Dynamic Games and Applications, 4(2), 110–154, 2014
  * [33] M. Huang, R. P. Malhamé and P. E. Caines, Mean Field Games, In T. Başar and G. Zaccour (eds) Handbook of Dynamic Game Theory, Springer, 2018
  * [34] R. Carmona and F. Delarue, Probabilistic Theory of Mean Field Games with Applications I: Mean Field FBSDEs, Control, and Games, Probability Theory and Stochastic Modelling 83, Springer, 2018
  * [35] R. Carmona and F. Delarue, Probabilistic Theory of Mean Field Games with Applications II: Mean Field Games with Common Noise and Master Equations Probability Theory and Stochastic Modelling 84, Springer, 2018
  * [36] P. Cardaliaguet, F. Delarue, J.-M. Lasry and P.-L. Lions, The master equation and the convergence problem in mean field games, Annals of Mathematics Studies 201, Princeton University Press, 2019
  * [37] P. Cardaliaguet and A. Porretta, An introduction to mean field game theory, in P. Cardaliaguet and A. Porretta (Eds.) Mean Field Games, Lecture Notes in Mathematics 2281 (Cetraro, Italy 2019), pp. 1– 158, 2020
  * [38] R. Naz, F. M. Mahomed and D. P. Mason, Comparison of different approaches to conservation laws for some partial differential equations in fluid mechanics, Appl. Math. Comput., 205 (1), 212–230, 2008
  * [39] E. Noether, Invariante variations problem. Konigliche Gesellschaft der Wissenschaften zu Gottingen, Nachrichten, Mathematisch-Physikalische Klasse Heft 2, pages 235–257, 1918. English translation: Transport Theory and Statist. Phys., 1(3), 1971, 183–207
  * [40] L. V. Ovsiannikov, Group Analysis of Differential Equations, Nauka, Moscow, 1978. English translation, W. F. Ames, Ed., published by Academic Press, New York, 1982
  * [41] N. H. Ibragimov, Transformation Groups Applied to Mathematical Physics, Reidel, Boston, 1985
  * [42] P. J. Olver, Applications of Lie Groups to Differential Equations, Springer, New York, 1986
  * [43] E. Bessel-Hagen, Über die Erhaltungssätze der Elektrodynamik, Math. Ann., 84(3-4), 258–276, 1921
  * [44] A. V. Bobylev and N. Kh. Ibragimov, Interconnectivity of symmetry properties for equations of dynamics, kinetic theory of gases, and hydrodynamics, Matem. Mod., 1(3), 100–109, 1989
  * [45] D. A. Gomes, L. Nurbekyan and M. Sedjro, One-dimensional forward–forward mean-field games, Applied Mathematics and Optimization, 74(3), 619–642, 2016
  * [46] D. Gomes, L. Nurbekyan and M. Sedjro, Conservation laws arising in the study of forward–forward mean-field games, Springer Proceedings in Mathematics and Statistics, 236, 643–649, 2018
  * [47] E. Hairer, C. Lubich and G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, Springer, 2006
  * [48] Y. Achdou and M. Laurière, Mean field games and applications: Numerical aspects, in P. Cardaliaguet and A. Porretta (Eds.) Mean Field Games, Lecture Notes in Mathematics 2281 (Cetraro, Italy 2019), pp. 249–307, 2020
  * [49] W. H. Fleming and R. W. Rishel, Deterministic and Stochastic Optimal Control, Springer Science & Business Media, Berlin, 2012
  * [50] J. Yong and X. Y. Zhou, Stochastic controls: Hamiltonian systems and HJB equations, Springer Science & Business Media, Berlin, 1999

## Appendices

##  Appendix A A sketch of the derivation of the MFG system

In this Appendix, we sketch the derivation of the system (2.1), which
describes Nash equilibrium in a game with a continuum of agents. The
presentation follows [37]. We refer to this book chapter for details and to
[49, 50] for the optimal control theory.

Each agent exercises control $\alpha$ in the stochastic differential equation
in $\mathbf{R}^{d}$

$dX_{s}=b(X_{s},\alpha_{s},m(s,X_{s}))ds+\sqrt{2\varepsilon}dB_{s},$ (A.1)

where $B_{s}$ is a Brownian motion. Here $m(t,x)$ is the density of the
probability measure in $\mathbf{R}^{d}$, which describes the distribution of
the agents. The agent minimizes

$\mbox{E}\left[\int_{0}^{T}L(X_{s},\alpha_{s},m(s,X_{s}))ds+G(X_{T},m(T,X_{T}))\
|\ X_{0}=x\right],$ (A.2)

where $L(x,\alpha,m)$ and $G(x,m)$ are running and terminal costs. They depend
on the agent’s state $x$ and the distribution of agents $m$. The running cost
$L$ also depends on the agent’s control $\alpha$.

This problem provides the value function

$u(t,x)=\mathop{\mbox{inf}}\limits_{\alpha}\mbox{E}\left[\int_{t}^{T}L(X_{s},\alpha_{s},m(s,X_{s}))ds+G(X_{T},m(T,X_{T}))\
|\ X_{t}=x\right]$ (A.3)

The infimum is taken over admissible controls. Function $X_{s}$ satisfies SDE
(A.1) with the initial condition $X_{t}=x$.

Each agent makes an infinitesimally small contribution to the population
dynamics. Therefore, the agent gets a problem of individual optimization with
given time dependent distribution of the other agents. For the value function,
we get PDE (2.1a) with

$H(x,p,m)=\mathop{\mbox{sup}}\limits_{\alpha}[-b(x,\alpha,m)\cdot
p-L(x,\alpha,m)].$ (A.4)

The optimal control has the feedback form $\alpha^{*}=\alpha^{*}(t,x)$ and

$b(x,\alpha^{*}(t,x),m(t))=-H_{p}(x,\nabla u,m(t)).$ (A.5)

If the agents exercise the optimal behaviour and their noises $B_{s}$ are
independent, the evolution of the agents’ distribution is given by the forward
Kolmogorov equations (2.1b)

##  Appendix B Preliminary analysis of the determining equations

Preliminary analysis identifies cases of the Hamiltonian function $H(u_{x})$,
which should be considered.

Lie point symmetries

$X=\xi^{t}{\partial\over\partial t}+\xi^{x}{\partial\over\partial
x}+\eta^{u}{\partial\over\partial u}+\eta^{m}{\partial\over\partial m}$ (B.1)

of the MFG systen

$\displaystyle F_{1}=-u_{t}-\varepsilon u_{xx}+H(u_{x})=f(m),$ (B.2a)
$\displaystyle F_{2}=m_{t}-\varepsilon
m_{xx}-m_{x}H^{\prime}(u_{x})-mH^{\prime\prime}(u_{x})u_{xx}=0$ (B.2b)

are specified by the determining equations

$\displaystyle E_{1}=\left.X(F_{1})\right|_{F_{1}=0,F_{2}=0}=0,$ (B.3a)
$\displaystyle E_{2}=\left.X(F_{2})\right|_{F_{1}=0,F_{2}=0}=0.$ (B.3b)

As discussed in Section 3, the generator $X$ is assumed to be prolonged to all
variables involved in the equations.

Application of the prolonged generator $X$ to the equations $F_{1}=0$ and
$F_{2}=0$ gives the system of two PDEs for the coefficients of the generator
$X$. The generator coefficients depend on the independent and dependent
variables $t$, $x$, $u$ and $m$. The system also depends on the functions
$H(u_{x})$ and $f(m)$ and their derivatives. Overall, the system depends on
the variables
$\\{t,x,u,m,u_{t},m_{t},u_{x},m_{x},u_{tx},m_{tx},u_{xx},m_{xx}\\}$. Since the
action of the prolonged generator $X$ to the MFG equations is considered on
the solutions of the MFG equations, we eliminate two of these variables.
Substituting

$m_{xx}={1\over\varepsilon}(m_{t}-m_{x}H^{\prime}(u_{x})-mH^{\prime\prime}(u_{x})u_{xx})$
(B.4)

and then

$u_{xx}={1\over\varepsilon}(-u_{t}+H(u_{x})-f(m)),$ (B.5)

we eliminate $u_{xx}$ and $m_{xx}$. These substitutions lead to the equations

$\displaystyle
E_{1}(t,x,u,m,u_{t},m_{t},u_{x},m_{x},u_{tx},m_{tx})=\left.X(F_{1})\right|_{(\ref{substitution_m}),(\ref{substitution_u})}=0,$
(B.6a) $\displaystyle
E_{2}(t,x,u,m,u_{t},m_{t},u_{x},m_{x},u_{tx},m_{tx})=\left.X(F_{2})\right|_{(\ref{substitution_m}),(\ref{substitution_u})}=0.$
(B.6b)

In detail, these equations take the form

$E_{1}=\left.\left(-\zeta_{t}^{u}-\varepsilon\zeta_{xx}^{u}+H^{\prime}(u_{x})\zeta_{x}^{u}-f^{\prime}(m)+\eta^{m}\right)\right|_{(\ref{substitution_m}),(\ref{substitution_u})}=0,$
$E_{2}=\left(\zeta_{t}^{m}-\varepsilon\zeta_{xx}^{m}-\zeta_{x}^{m}H^{\prime}(u_{x})-m_{x}H^{\prime\prime}(u_{x})\zeta_{x}^{u}\right.\\\
\left.\left.-\eta^{m}H^{\prime\prime}(u_{x})u_{xx}-mH^{\prime\prime\prime}(u_{x})\zeta_{x}^{u}u_{xx}-mH^{\prime\prime}(u_{x})\zeta_{xx}^{u}\right)\right|_{(\ref{substitution_m}),(\ref{substitution_u})}=0.$

Clearly, these equations can be split for $u_{t}$, $m_{t}$, $m_{x}$, $u_{tx}$
and $m_{tx}$ because the generator coefficients, $H(u_{x})$ and $f(m)$ do not
depend on these variables.

We proceed as follows.

  1. 1.

From

${\partial E_{1}\over\partial
u_{tx}}=2\varepsilon(\xi^{t}_{x}+\xi^{t}_{u}u_{x}+\xi^{t}_{m}m_{x})=0,$

where $\xi^{t}=\xi^{t}(t,x,u,m)$, we obtain

$\xi^{t}=A(t),$

where $A(t)$ is an arbitrary function.

  2. 2.

The equation

${\partial^{2}E_{1}\over\partial
m_{x}^{2}}=2\varepsilon(\xi^{x}_{mm}u_{x}-\eta^{u}_{mm})=0$

for $\xi^{x}=\xi^{x}(t,x,u,m)$ and $\eta^{u}=\eta^{u}(t,x,u,m)$ gives

$\xi^{x}=B_{1}(t,x,u)m+B_{0}(t,x,u)$

and

$\eta^{u}=C_{1}(t,x,u)m+C_{0}(t,x,u).$

  3. 3.

The equation

${\partial^{4}E_{2}\over\partial u_{x}^{2}\partial
m_{x}^{2}}=2(B_{1}(t,x,u)u_{x}-C_{1}(t,x,u))H^{\prime\prime\prime\prime}(u_{x})=0$

identifies the special case of the Hamiltonian function

$H(u_{x})=h_{3}u_{x}^{3}+h_{2}u_{x}^{2}+h_{1}u_{x}+h_{0}.$

This case will be treated separately. It is convenient to split this case into
two: the case of the cubic Hamiltonian

$H(u_{x})=hu_{x}^{3}+h_{2}u_{x}^{2}+h_{1}u_{x}+h_{0},\qquad h\neq 0$ (B.7)

and the case of the quadratic Hamiltonian

$H(u_{x})=hu_{x}^{2}+h_{1}u_{x}+h_{0},\qquad h\neq 0.$ (B.8)

Linear Hamiltonians are excluded by requirement (2.9).

We continue with the assumption that the Hamiltonian has a form different from
(B.7) and (B.8). It follows

$B_{1}(t,x,u)=0,\qquad C_{1}(t,x,u)=0.$

  4. 4.

The equation

${\partial^{2}E_{2}\over\partial m_{x}^{2}}=-2\varepsilon\eta^{m}_{mm}=0$

for $\eta^{m}=\eta^{m}(t,x,u,m)$ provides

$\eta^{m}=G_{1}(t,x,u)m+G_{0}(t,x,u).$

  5. 5.

From

${\partial^{2}E_{2}\over\partial u_{x}\partial m_{t}}=2(B_{0}(t,x,u))_{u}=0$

we get

$B_{0}(t,x,u)=\tilde{B}_{0}(t,x).$

After all these simplifying steps, we obtain the symmetry coefficients in the
form

$\xi^{t}=A(t),\qquad\xi^{x}=\tilde{B}_{0}(t,x),\qquad\eta^{u}=C_{0}(t,x,u),\qquad\eta^{m}=G_{1}(t,x,u)m+G_{0}(t,x,u).$

These such coefficients we obtain the equation

${\partial^{3}E_{2}\over\partial u_{x}^{2}\partial
m_{x}}=(F_{1}(t,x,u)u_{x}+F_{2}(t,x,u))H^{\prime\prime\prime\prime}(u_{x})+F_{3}(t,x,u)H^{\prime\prime\prime}(u_{x})=0,$

where functions $F_{1}$, $F_{2}$ and $F_{3}$ depend on the derivatives of the
functions $\tilde{B}_{0}(t,x)$ and $C_{0}(t,x,u)$. This equation is the
classifying equation for the Hamiltonian function. In addition to the cases
(B.7) and (B.8) found earlier, we get two more cases

$H(u_{x})=h(u_{x}+q)^{p}+h_{2}u_{x}^{2}+h_{1}u_{x}+h_{0},\qquad h\neq 0,\qquad
p\neq 0,1,2,3$ (B.9)

and

$H(u_{x})=he^{ku_{x}}+h_{2}u_{x}^{2}+h_{1}u_{x}+h_{0},\qquad h\neq 0,\qquad
k\neq 0.$ (B.10)

Further examination of the determining equations does not provide any other
special cases of the Hamiltonian function. Thus, in addition to the special
cases (B.7), (B.8), (B.9) and (B.10), we have only the general case of the
arbitrary $H(u_{x})$.

This is the end of the preliminary analysis of the determining equations.
After this we proceed with the specified forms of the Hamiltonian function. We
should return to the determining equations (B.3) and solve them for the
different cases of $H(u_{x})$. It provides the admitted symmetries and the
special cases of function $f(m)$.


format=plain,labelformat=simple, labelsep=period

#  Connect Everywhere: Wireless Connectivity in Protected Areas

Thembelihle Dlamini, Mengistu Abera Mulatu  Department of Electrical and
Electronic Engineering
University of Eswatini
Manzini, Eswatini
{tldlamini<EMAIL_ADDRESS>Sifiso Vilakati  Department of Statistics and
Demography
University of Eswatini
Manzini, Eswatini
<EMAIL_ADDRESS>
###### Abstract

Mobile connectivity has become more important, especially for visitors to
parks and protected areas. However, governments’ policies prohibit cable
wiring in these areas in order to preserve the beauty and protect the
historical interest of the landscape. As a result, in order to provide
emergency and infotainment services, mobile connectivity is of primary
importance in areas such as national parks and historical sites. Through
observed practices from other countries, mobile network operators and other
licensed service providers can cooperate with the administrators of such
protected areas in order to provide mobile connectivity without disturbing the
environment. However, the most pervasive problem is the high energy
consumption of the wireless systems and it becomes expensive to power them
using the electricity grid. One attractive solution is to make use of green
energy to power the communication systems and then share the base station (BS)
infrastructure that is co-located with the edge server, an entity responsible
for computing the offloaded delay-sensitive workloads. To alleviate this
problem, this paper offers a resource management solution that seeks to
minimize the energy consumption per communication site through (i) BS
infrastructure and resource sharing, and (ii) assisted (peer-to-peer)
computation offloading for the energy-constrained communication sites, within
a protected area. Using this resource management strategy guarantees a Quality
of Service (QoS). The performance evaluation conducted through simulations
validates our analysis as the prediction variations observed shows greater
accuracy between the harvested energy and the traffic load. Towards energy
savings, the proposed algorithm achieves a $52\%$ energy savings when compared
with the $46\%$ obtained by our benchmark algorithm. The energy savings that
can be achieved decreases as the QoS is prioritized, within each communication
site, and when the number of active computing resources increases.

###### Index Terms:

Protected areas, edge computing, lookahead, green energy.

##  I Introduction

Communication has become the most critical aspect of human life. Hence, mobile
connectivity is crucial in offering the required communication for humans
wherever they are. Multi-connectivity solutions have been suggested for
wireless users in different kinds of environments, however, due to certain
governments’ policies, the same cannot be said for protected areas [1]. Most
governements do not approve cable wiring in such protected areas in order to
preserve their beauty. Thus, in order to guarantee mobile connectivity in
protected areas, i.e., the provision of emergency services and the deployment
of energy self-sufficient communication sites, future Mobile Networks are
expected to leverage the integration of Multi-access Edge Computing (MEC) and
Energy Harvesting (EH) BS, i.e., the Base Stations are empowered with
computing capabilities and also equipped with EH equipments [2]. The use of EH
systems motivates the need for energizing edge systems with green energy in
order to extend network coverage to protected areas, and also minimizing the
carbon footprint [3]. The renewable-powered BSs are more suitable in national
parks and historical sites, when considering environmental impact, as most
governments do not approve cable wiring in order to preserve the landscape
[1].
Despite of the potential presented by the integration of MEC and EH BS, the
challenges of both resource provisioning and energy consumption still come up
under this new paradigm of softwarized MNs. In addition, the variation of the
haversted energy per communication site brings about the notion of assisted
computation offloading, i.e., one BS accept the offloaded delay-sensitive
workload and then cooperate with neighboring BSs that have enough green energy
and available computing resources (Virtual Machines or containers) in order
for them to compute the task and return the result to it. The task offloading
technique helps to avoid large computation latency at overloaded BSs or
energy-deficient BSs that might degrade the Quality of Service (QoS) to end
users [4]. Therefore, in order to address the connectivity issue, this paper
formulates an optimization problem to manage a shared BS system network (the
BS infrastructure and its co-located computing platform (MEC server)) that is
deployed within a national park. The primary objective of this paper is to
minimize energy consumption through infrastructure sharing and resource
management procedures, and assisted computation offloading for the energy-
constrained BS sites. Here, the offloaded workloads are jointly allocated over
the BS sites and the computing platform resources, in an energy-efficient
manner with a guarantee of QoS. The formulated problem is solved by exploiting
predictions of the BS traffic load, harvested energy, and computation load via
a Limited Lookahead Control (LLC) approach, obtaining energy savings higher
than $50\%$ with respect to myopic allocation strategies, i.e., with respect
to communication site management schemes that do not have lookahead
capabilities.
Related work: Various research contributions for the management of renewable-
powered BS networks that are empowered with computing capabilities exist in
literature. The authors of [5] incorporate renewables into MEC system and
propose an efficient reinforcement learning-based resource management
algorithm for handling dynamic workload offloading. Then, in [2], the authors
propose a controller-based network architecture for managing EH BSs empowered
with computation capabilities where on/off switching strategies allow BSs and
VMs to be dynamically switched on/off towards energy savings, over a limited
prediction horizon. In [6], the green-based load balancing technique is
employed for optimizing MEC performance by exploiting the spatial diversity of
the available green energy to reshape the network load among the BSs. Here,
container-based virtualization was considered. Towards infrastructure sharing,
the work of [7] employs an infrastructure sharing algorithm towards energy
savings. In their work, a game-theoretic framework was proposed in order to
enable the MN operators to individually estimate the switching-off
probabilities that reduce their expected financial cost. Then, in [8], the
problem of infrastructure sharing among MN operators is presented as a
multiple-seller single-buyer business. In their contribution, each BS is
utilized by subscribers from other operators and the owner of the BS is
considered as a seller of the BS infrastructure while the owners of the
subscribers utilizing the BS are considered as buyers. In the presence of
multiple seller MN operators, it is assumed that they compete with each other
to sell their network infrastructure resources to potential buyers. However,
while performing task offloading, care must be taken in order to avoid large
computation delays in overloaded BSs. To avoid these large computation
latencies at overloaded BSs, assisted computation must be employed within the
MEC paradigm. The works of [4] develops a novel online BS peer offloading
framework. Here, a Lyapunov technique is used in order to maximize the long-
term system performance while keeping the energy consumption of BSs below
individual long-term constraints. In this regard, it is noted that the
aforementioned works are not fully offgrid similar to the national park
environment that requires no power lines. Also, the joint consideration of
infrastructure and resource sharing, and assisted computation offloading for
energy-constrained communication sites, within the MEC paradigm has not been
given enough research attention.

Figure 1: The BS system within the protected area where energy is harvested
from the environment.

In order to achieve the objective of this paper, the remainder of this article
is organized as follows: In Section II, the system model is described,
formulating the optimization problem and solving it in an online fashion via
an LLC-based approach. In Section III, some selected results are shown,
quantifying the energy savings and delay cost that can be achieved with our
framework. Finally, in Section IV some concluding remarks are given.

##  II System Model

We consider $n\in\mathcal{N}$ BSs, indexed by $\mathcal{N}=\\{1,\dots,N\\}$,
deployed in a protected area (national park in our case), each endowed with
computing capabilities (BS is co-located with a MEC server that runs $D$
containers), as illustrated in Fig. 1 above. Hence, User Equipments can
offload their computation jobs to corresponding serving BS via wireless
communications for processing. In order to preserve the landscape, the
deployed BSs have greater level of camouflage. Each communication site is
mainly powered by renewable energy harvested from wind and solar radiation,
and it is equipped with an Energy Buffer (EB) (with a maximum capacity of
$B_{n}^{\rm max}$) for energy storage. There is an Energy Manager (EM) which
is responsible for selecting the appropriate energy source to fulfill the EB,
and also for monitoring the energy level of the EB. In addition, the MEC
server house a virtualized Access Control Router (ACR) application which acts
as an access gateway, responsible for accepting and routing the workloads that
are either accepted locally or offloaded to neighboring BSs. For remote clouds
or Internet access, a microwave backhaul or a multi-hop wireless backhaul
relaying (e.g., integrated access and backhaul (IAB)) is used. Moreover, a
discrete-time model is considered, whereby the time is discretized as
$t=1,2,\dots$ time slots of a fixed duration $\tau=$30\text{\,}\mathrm{min}$$.
Communication and computing energy cost: it is worth noting that the total
energy consumed [$\text{\,}\mathrm{J}$] per communication site
$\theta_{n}^{\rm site}(t)$ consist of the BS communications, denoted by
$\theta_{n}^{\rm bs}(t)$, and computing platform processes, related to
computing, caching, and communication, which is denoted by $\theta_{n}^{\rm
mec}(t)$. Thus, the energy consumption model at time slot $t$ is formulated as
follows, inspired by [6]:

$\theta_{n}^{\rm site}(t)=\theta_{n}^{\rm bs}(t)+\theta_{n}^{\rm mec}(t).$ (1)

$0$$0.2$$0.4$$0.6$$0.8$$1$$5$$10$$15$$20$$25$$30$$35$$40$$45$

Normalized traffic load

Time [h]Cluster 1Cluster 2Cluster 3Cluster 4

Figure 2: Normalized BS traffic loads behavior represented as clusters. The
data from [9] has been split into four representative clusters.

Energy levels and job flows: within each BS site there is an intelligent
electro-mechanical switch (I-SW) that aggregates the energy sources to fulfill
the EB level. In this work, the EBs work as either an energy sources or sinks,
depending on traffic load $L_{n}(t)$ (see Fig. 2 for traffic load patterns)
which consists of delay-sensitive workloads $l_{n}(t)$. At time slot $t$, the
amount of energy required for computation and communication process from the
EBs, denoted by $b_{n}(t)\geq\theta_{n}^{\rm site}(t)$, is made up of the
harvested energy (see Fig. 3 above for energy profiles) as follows:

$h_{n}(t)=h_{n}^{c}(t)+h_{n}^{o}(t)\,,$ (2)

where $h_{n}^{c}(t)$ is the portion of harvested energy used for charging the
battery and $h_{n}^{o}(t)$ is the portion of harvested energy immediately used
by the site to support its operations. It is worth noting that the actual
amount of harvested energy is bounded by the maximum amount of energy that
could be harvested from the environment, denoted by $H_{n}^{\rm max}$, at a
certain time slot. Thus, the following hard constraint must hold:

$h_{n}^{c}(t)+h_{n}^{o}(t)\leq H_{n}^{\rm max}\,,\forall n,\forall t.$ (3)

The available EB level, at time slot $t$, located at the BS site (BS $n$)
evolves according to the following dynamics:

$b_{n}(t)=\mu_{n}(b_{n}(t-1)-\theta_{n}^{\rm
site}(t))+\alpha_{n}(h_{n}^{c}(t))\,,$ (4)

where $\mu_{n}\in(0,1]$ is a parameter accounting for the self-discharging
behavior of the battery, and $\alpha_{n}\in(0,1]$ accounts for the losses in
the charging process.

Figure 3: Harvested solar and wind traces [10] used in our simulations.

In the interest of network optimization through infrastructure sharing and
resource allocation, over a given horizon, the BSs can perform assisted
computation by offloading in full or part the computational tasks, accepted
from the connected UEs, to neigboring BSs with sufficient green energy. We
define ${w}_{nm}(t)$ (jobs/slot) as the job flow that BS site
$n\in\mathcal{N}$ offloads towards BS site $m$ (one among its neighboring BSs)
in time slot $t$. Note that ${w}_{nn}(t)$ represents the portion of job flow
that is directly collected by MEC server $n$ that is co-located with the BS.
At this regard, we also note that the energy required by the BS site $n$ in
time slot $t$, obtained from the harvesting systems or EB, is linear and
directly proportional to the job flow as follows:

$b_{n}(t)=\eta\left(w_{nn}+\sum_{m\in\mathcal{N}\setminus\\{n\\}}w_{mn}\right)+\kappa\left(\sum_{m\in\mathcal{N}\setminus\\{n\\}}w_{nm}\right).$
(5)

Note that $\eta>0$ while $\kappa<0$ (with $|\eta|>|\kappa|$). Since the
admitted computation workloads (jobs) are first buffered at the input buffer
at time $t$ and the results are accumulated over a fixed period of time to
form a batch at the output buffer. From this, we note that at every time slot
the number of jobs exiting a node cannot exceed the number of jobs entering
the same node. Thus, the following condition holds:

$\sum_{m\in\mathcal{N}\setminus\\{n\\}}(w_{nm}(t)-w_{mn}(t))\leq
w_{nn}(t)\,,\forall t,\forall n.$ (6)

Computing resources and load distribution: the amount of the accepted workload
$w_{nn}(t)$ is computed using containers as computing resources within the MEC
server. For a fair provisioning of the computing resources, $D(t)$ needs to be
obtained first, and then the workload per container $\lambda_{d}(t)$ is
determined. Firstly, each container can only compute an amount of up to
$\lambda_{\max}$ (considering that virtualization technologies specify the
minimum and maximum amount of load that can be allocated per container) and to
meet the latency requirements, $D(t)$ is obtained as:
$D(t)=\big{\lceil}(w_{nn}(t)/\lambda_{\max})\big{\rceil}$, where
$\big{\lceil}\cdot\big{\rceil}$ returns the nearest upper integer. Secondly,
to distribute the workload among the $D(t)$ containers, a heuristic process
splits the computational workload $\lambda_{d}(t)=\lambda_{\max}$ to the first
$D(t)-1$ containers, and the remaining workload
$\lambda_{d}(t)=w_{nn}(t)-(D(t)-1)\lambda_{\max}$ to the last one.

###  II-A Optimization

Our objective is to improve the overall energy savings of the communication
system through infrastructure sharing principles and resource management
procedures, with a QoS guarantee. The power saving modes within each BS site
is achieved by enabling assisted computation and also performing the following
within the MEC server: autoscaling of containers, contents caching and tuning
of the transmission drivers. To achieve our objective, two cost functions are
defined, one captures the communication system energy consumption and the
other, handles the QoS. This is defined as follows: F1) $\theta_{n}^{\rm
site}(t)$, weighs the energy consumption due to transmission in the BSs and
the computing-plus-communication activities in the MEC server. Then, F2) a
quadratic term $(w_{nn}(t)-w_{nm}(t))^{2}$, which accounts for the QoS. At
this regard, it is worth noting that F1 tends to push the system towards self-
sustainability solutions and F2 favors solutions where the delay sensitive
load is either entirely admitted in the co-located MEC server by the router
application or offloaded to a neighboring BS that is not energy-constrained. A
weight $\Gamma=[0,1]$ is utilized to balance the two objectives F1 and F2. The
corresponding (weighted) cost function is defined as:

$\displaystyle J_{n}(t)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}\overline{\Gamma}\,\theta_{n}^{\rm
site}(t)+\Gamma\,(w_{nn}(t)-w_{nm}(t))^{2}\,,$ (7)

where $\overline{\Gamma}\stackrel{{\scriptstyle\Delta}}{{=}}1-\Gamma$. Hence,
starting from $t=1$ (i.e., $t=1,2,\dots,T$) as the current time slot and the
finite horizon $T$, the following optimization problem is formulated as:

P1 $\displaystyle:$
$\displaystyle\min_{h_{n}^{c}(t),h_{n}^{o}(t),D(t),w_{nm}(t)}\sum_{t=1}^{T}J_{n}(t)$
$\displaystyle\mbox{subject to:}\quad\text{Eqn.}(2),(3),(4),(5),(6).$

To solve P1 in (II-A), the LLC principles [11], Geometric Programming (GP)
technique [12], and heuristics, is used towards obtaining the feasible system
control inputs.

###  II-B Limited lookahead control approach

To solve P1 (Eq. (II-A)), we are required to have a complete knowledge of the
traffic load, harvested energy arrivals, and the workloads accepted by the MEC
server. However, such knowledge is not possible to achieve, in fact, we only
have knowledge of the exogenous processes in the current, $t$, and past time
slots. Thus, we adopt a pragmatic LLC approach to solve the problem in an
online fashion. Next, we define the system state vector as $q_{n}(t)$ and the
input vector as $\varphi_{n}(t)$. The system behavior is described by the
discrete-time state-space equation, adopting the LLC principles [13], as:

$\ q_{n}(t+1)=\phi(q_{n}(t),\varphi_{n}(t))\,,$ (9)

where $\phi(\cdot)$ is a behavioral model that captures the relationship
between $(q_{n}(t),\varphi_{n}(t))$, and the next state $q_{n}(t+1)$. Note
that this relationship accounts for the amount of energy drained, that
harvested $h_{n}(t)$, which together lead to the next buffer level
$b_{n}(t+1)$. We note that state $q_{n}(t)$ and $\varphi_{n}(t)$ are
respectively measured and applied at the beginning of the time slot $t$,
whereas the offered load $L_{n}(t)$ and the harvested energy $h_{n}(t)$ are
accumulated during the time slot and their value becomes known only by the end
of it. This means that, being at the beginning of time slot $t$, the system
state at the next time slot $t+1$ can only be estimated, which is formally
written as:

$\hat{q}_{n}(t+1)=\phi(q_{n}(t),\varphi_{n}(t))\,.$ (10)

For these estimations, we use the forecast values of load $\hat{L}_{n}(t)$ and
harvested energy $\hat{h}_{n}(t)$, from the Long Short-Term Memory (LSTM)
forecasting module (see Section II-C for additional details on forecasting).
Specifically, for each time slot $t$, problem (II-A) is solved obtaining
control actions for the prediction horizon $T$. The control action that is
applied at time $t$ is $\varphi_{n}^{*}(t)$, which is the first one in the
retrieved control sequence.

###  II-C Prediction of exogenous processes

In this section, we describe how the traces used in the simulations were
collected and how their predictions were performed.
Traffic load and Harvested energy: in this work, the amount of harvested
energy $h_{n}(t)$ in time slot $t$ is obtained from open-source solar and wind
traces within farm located belgium [10]. As for the traffic load, open source
MN datasets obtained from the Telecom Italia Mobile (TIM) network (availed
through the Big Data Challenge [14]) are used to emulate the traffic and
computational load. Two exogenous processes are considered in this work: the
harvested energy $h_{n}(t)$ and the BS traffic loads $L_{n}(t)$. In order to
generate the predictions ($\hat{h}_{n}(t),\hat{L}_{n}(t)$), the LSTM neural
networks [15] were adopted. Thus, the LSTM-based predictor has been trained to
give an output of the the forecasts for the required number of future time
slots $T$. The trained LSTM network consists of an input layer, a single
hidden layer consisting of $40$ neurons, for $80$ epochs, for a batch size of
$4$; and an output layer. For training and testing purposes, the dataset was
split as $70\%$ for training and $30\%$ for testing. As for the performance
measure of the model, the Root Mean Square Error (RMSE) is used.
Workload (Jobs) flows: in order to understand the daily traffic load patterns,
the clustering algorithm X-means [16] has been applied to classify the load
profiles into several categories. Here, each BS $n$ is assumed to have a
related load profile ${L}_{n}(t)$ which is picked at random as one of the four
clusters in Fig. 2. In addition, it is assumed that ${L}_{n}(t)$ consists of
$80\%$ delay sensitive workloads $l_{n}(t)$ and the remainder is delay
tolerant workloads. Since the amount of the workload to be computed locally,
$w_{nn}(t)$, is determined by considering the amount of energy required by BS
$n$ at time slot $t$, $b_{n}(t)$, and the next time slot energy to be
accumulated through harvesting process, denoted by $b_{n}(t+1)$. Note that the
harvested energy is predicted. Thus,
$w_{nn}(t)=\frac{b_{n}(t)}{b_{n}(t+1)}l_{n}(t)$.

###  II-D The online algorithm

The online algorithm proceeds as follows, in a distributed manner: Starting
from the initial state, the online algorithm constructs, in a breadth-first
fashion, a tree comprising all possible future states up to the prediction
depth $T$. Then, a search set $\mathcal{F}$ consisting of the current system
state is initialized, and it is accumulated as the algorithm traverse through
the tree, accounting for predictions, accumulated workloads at the output
buffer, past outputs and controls. The set of states reached at every
prediction depth $t+p$ is referred to as $\mathcal{F}(t+p)$. Given $q(t)$, we
first estimate the traffic load $\hat{L}_{n}(t+p)$, delay-sensitive jobs
$\hat{l}_{n}(t+p)$, locally acceptable computational load $\hat{w}_{nn}(t+p)$,
harvested energy $\hat{h}_{n}(t+p)$, any expected acceptable computational
jobs from neighbors $\hat{w}_{nm}(t+p)$, and generate the next set of
reachable control actions by applying the input workload and energy harvested.
The energy cost function corresponding to each generated state $q_{n}(t+p)$ is
then computed. Once the prediction horizon is explored, a sequence of
reachable states yielding minimum energy consumption is obtained. The control
action $\varphi_{n}^{*}(t)$ corresponding to $\hat{q}_{n}(t+p)$ (the first
state in this sequence) is provided as input to the system while the rest are
discarded. The process is repeated at the beginning of each time slot $t$.

Table I: Simulation parameters. Parameters depending on $n$ apply to all nodes $n\in\mathcal{N}$. Parameter | Value
---|---
$N$ | $20$
$D$ | $20$
$B_{n}^{\rm max}$ | $100\text{\,}\mathrm{kJ}$
$\mu_{n}$ | 0.9999
$\alpha_{n}$ | 0.900
$\eta$ |  $0.105\text{\,}\mathrm{J}$/(job/slot)
$\kappa$ | -$0.035\text{\,}\mathrm{J}$/(job/slot)
$\lambda_{\rm max}$ | 10 MB

##  III Performance Evaluations

In this section, some selected numerical results for the scenario of Section
II are shown. The parameters that were used in the simulations are listed in
Table I above. Our time slot duration $\tau$ is set to
$30\text{\,}\mathrm{min}$ and the time horizon is set to $T=3$ time slots. For
simulation, Python is used as the programming language.

$0$$0.2$$0.4$$0.6$$0.8$$1$$5$$10$$15$$20$$25$$30$$35$$40$$45$

Normalized values

Time [h]Solar 3 (real)Solar 3 (pred)Cluster 2 (real)Cluster 2 (pred)Wind 3
(real)Wind 3 (pred)

Figure 4: Forecast mean value for $L_{n}(t)$ and $h_{n}(t)$.

Forecasting: In Fig. 4 above, we show the real and forecasted values for the
mobile traffic load (Cluster 2) and harvested energy (Solar 3 and Wind 3) over
the time horizon. We show the one-step predictive mean value at each step of
the online forecasting routine. The obtained average RMSE for the traffic load
and harvested energy processes, both normalized in [0,1] for
$T\in\\{1,2,3\\}$, are $L_{n}(t)=\\{0.013,0.015,0.020\\}$, for solar
$h_{n}(t)=\\{0.030,0.035,0.037\\}$, and for wind
$h_{n}(t)=\\{0.032,0.036,0.039\\}$. Note that the forecasting for $L_{n}(t)$
are more accurate than those of $h_{n}(t)$. This is confirmed by comparing the
average RMSE. The measured accuracy is deemed good enough for the proposed
online optimization.

$0$$10$$20$$30$$40$$50$$60$123456789101112131415161718192021222324

Mean energy savings $[\%]$

Time [h]LLCOPEN

Figure 5: Mean energy savings.

Energy consumption: Our online algorithm (LLC) is benchmarked with another
one, named OPEN (Online BS PEer offloadiNg) from [4]. To address the energy
consumption problem, OPEN make use of the Lyapunov technique and solve the
problem by using only the current information (no forecasting in this case).
Fig. 5 above shows the energy savings obtained for one BS with respect to the
case where no energy management is performed; i.e., the network is dimensioned
for maximum expected capacity. The LLC achieves a energy savings of $52\%$ and
OPEN achieves $46\%$. The hourly obtained energy savings corresponds to the
behavior of national park visitor. In the early morning hours, there are few
visitors resulting into higher energy savings
($$6\text{\,}\mathrm{h}$-$8\text{\,}\mathrm{h}$$). During the midday (busy
period), the number of visitors arriving increases resulting into reduced
amount of energy savings ($$9\text{\,}\mathrm{h}$-$17\text{\,}\mathrm{h}$$).
At this time period, assisted computation offloading is more prevalent. After
$17\text{\,}\mathrm{h}$ there is less activity or no activity thus resulting
into an increase in energy savings.
Quality of Service (QoS): Then, Fig. 6 above shows the average energy savings
with respect to $\Gamma$ when the BSs ($N=20$) are running the LLC algorithm
in a distributed manner. Here, a trade-off is obtained by either computing all
or partial the offload jobs, resulting into either an increase or decrease of
energy savings over the protect area. Again, here the energy savings are
obtained with respect to the case where all the BSs are dimensioned for
maximum capacity. As expected, there is a drop in the energy savings achieved
as the value of $\Gamma$ increases, as QoS is prioritized. It can be observed
that LLC achieves a value of $50\%$ or above when $\Gamma=[0,0.5]$ and OPEN
achieves energy savings of above $50\%$ when $\Gamma=[0,0.2]$. From this we
can observe that LLC outperforms OPEN, thanks to assisted computation
offloading and the foresighted optimization which helps to guarantee low
latency services.

$30$$35$$40$$45$$50$$55$$60$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$

Mean Energy Savings $[\%]$

Weight, $\Gamma$LLCOPEN

Figure 6: Overall energy savings versus weighing parameter.

MEC server energy consumption: Finally, figure 7 above shows the energy
drained within the computing platform when the number of active containers is
varied. Here, we observe that the lower the number of active containers the
lower the energy drained within the server. This correspond to cases where the
computational workload $w_{nn}(t)$ is low (instances of low CPU utilization)
or some of the workload have been offloaded to neighboring peers (i.e.,
assisted computing). From the results, it is observed that LLC consumed less
energy when compared with OPEN.

##  IV Conclusions

The challenge of providing connectivity to protected areas will be one of the
pillars for future mobile networks. To address the connectivity issue, this
paper formulates an optimization problem to manage a shared base station
network (the base station infrastructure and its co-located computing platform
that is deployed within a national park. An online algorithm based on
infrastructure sharing, assisted computation offloading, forecasting,
lookahead principles and heuristics, is proposed with the objective of saving
energy within the base stations deployed in a protected area. Numerical
results, obtained with real-world energy and traffic load traces, demonstrate
that the proposed algorithm achieves mean energy savings of $52\%$ when
compared with the $46\%$ obtained by our benchmark algorithm. The online
algorithm achieves energy savings of about $50\%$ or above when
$\Gamma=[0,0.5]$ and the benchmark achieves energy savings of above $50\%$
when $\Gamma=[0,0.2]$. The energy saving results are obtained with respect to
the case where no energy management techniques are applied in the base
stations. The obtained results show that there exists a trade-off between
energy saving and the QoS. In addition, it is observed that the lower the
number of active containers the lower the energy drained within the computing
platform, and the increase in the number of active containers results in an
increase in the energy being consumed. In this case, assisted computation
helps to reduce the energy drained as some of the load is computed by
neighboring base stations.

$10$$20$$30$$40$$50$$60$$2$$4$$6$$8$$10$$12$$14$$16$$18$$20$

$\theta_{n}^{mec}$[J]

Containers, $D(t)$LLCOPEN

Figure 7: Energy consumption within the computing platform.

## References

  * [1] “BEREC and RSPG joint report on Facilitating mobile connectivity in “challenge areas”,” BEREC, Rīga, Latvia, Tech. Rep., Dec 2017.
  * [2] T. Dlamini, Á. F. Gambín, D. Munaretto, and M. Rossi, “Online Supervisory Control and Resource Management for Energy Harvesting BS Sites Empowered with Computation Capabilities,” _Wireless Communications and Mobile Computing_ , 2019.
  * [3] D. Thembelihle, M. Rossi, and D. Munaretto, “Softwarization of Mobile Network Functions towards Agile and Energy Efficient 5G Architectures: A Survey,” _Wireless Communications and Mobile Computing_ , 2017.
  * [4] Chen, Lixing and Zhou, Sheng and Xu, Jie, “Computation peer offloading for energy-constrained mobile edge computing in small-cell networks,” _IEEE/ACM Transactions on Networking_ , vol. 26, no. 4, pp. 1619–1632, 2018\.
  * [5] X. Jie and R. Shaolei, “Online Learning for Offloading and Autoscaling in Renewable-Powered Mobile Edge Computing,” in _IEEE Global Communications Conference (GLOBECOM)_ , Washington, USA, Dec. 2012.
  * [6] T. Dlamini and S. Vilakati, “LSTM-based Traffic Load Balancing and Resource Allocation for an Edge System,” _Wireless Communications and Mobile Computing_ , 2020.
  * [7] Alexandra Bousia and Elli Kartsakli and Angelos Antonopoulos and Luis Alonso and Christos Verikoukis, “Game-Theoretic Infrastructure Sharing in Multioperator Cellular Networks,” _IEEE Transctions on Vehicular Technology_ , vol. 65, no. 5, 2016.
  * [8] Sanguanpuak Tachporn and Guruacharya Sudarshan and Hossain Ekram and Rajatheva Nandana and Latva-aho Matti, “Infrastructure sharing for mobile network operators: analysis of trade-offs and market,” _IEEE Transactions on Mobile Computing_ , vol. 17, no. 12, pp. 2804–2817, 2018.
  * [9] Mobile and Energy datasets. [Online]. Available: https://github.com/lihles/mobile-datasets
  * [10] “Total generation.” [Online]. Available: https://www.elia.be/en/grid-data/power-generation
  * [11] J. P. Hayes, “Self-Optimization in Computer Systems via On-Line Control: Application to Power Management,” in _Proceedings of the First International Conference on Autonomic Computing_ , Washington,USA, May 2004\.
  * [12] W.-C. Ho, L.-P. Tung, T.-S. Chang, and K.-T. Feng, “Enhanced component carrier selection and power allocation in LTE-advanced downlink systems,” in _2013 IEEE Wireless Communications and Networking Conference (WCNC)_ , Shanghai, China, April 2013.
  * [13] S. Abdelwahed, N. Kandasamy, and S. Neema, “Online control for self-management in computing systems,” in _IEEE Real-Time and Embedded Technology and Applications Symposium (RTAS)_ , Ontario, Canada, May 2004.
  * [14] Open Big Data Challenge. [Online]. Available: https://dandelion.eu/datamine/open-big-data/
  * [15] I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_. MIT Press, 2016.
  * [16] D. Pelleg, A. W. Moore _et al._ , “X-means: Extending K-means with efficient estimation of the number of clusters,” in _Proceedings of the Seventeenth International Conference on Machine Learning (ICML)_ , San Francisco, USA, Jun 2000.

  *[MEC]: Multi-access Edge Computing
  *[EH]: Energy Harvesting
  *[BSs]: Base Station
  *[BS]: Base Station
  *[MNs]: Mobile Network
  *[QoS]: Quality of Service
  *[LLC]: Limited Lookahead Control
  *[VMs]: Virtual Machine
  *[EB]: Energy Buffer
  *[GP]: Geometric Programming
  *[LSTM]: Long Short-Term Memory
  *[TIM]: Telecom Italia Mobile
  *[RMSE]: Root Mean Square Error