# Mean Field Models to Regulate Carbon Emissions in
Electricity Production

René Carmona111Department of Operations Research and Financial Engineering,
Princeton University, Princeton, NJ 08544<EMAIL_ADDRESS><EMAIL_ADDRESS>lauriere@princeton.edu). Gökçe Dayanıklı 11footnotemark:
1 Mathieu Laurière 11footnotemark: 1

###### Abstract

The most serious threat to ecosystems is the global climate change fueled by
the uncontrolled increase in carbon emissions. In this project, we use mean
field control and mean field game models to analyze and inform the decisions
of electricity producers on how much renewable sources of production ought to
be used in the presence of a carbon tax. The trade-off between higher revenues
from production and the negative externality of carbon emissions is quantified
for each producer who needs to balance in real time reliance on reliable but
polluting (fossil fuel) thermal power stations versus investing in and
depending upon clean production from uncertain wind and solar technologies. We
compare the impacts of these decisions in two different scenarios: 1) the
producers are competitive and hopefully reach a _Nash Equilibrium_ ; 2) they
cooperate and reach a _Social Optimum_. We further introduce and analyze the
impact of a regulator in control of the carbon tax policy, and we study the
resulting Stackelberg equilibrium with the field of producers.

Keywords. Mean field games; Mean field control; Carbon emission; Stackelberg
equilibrium

_Acknowledgments._ The authors were partially supported by NSF DMS-1716673,
ARO W911NF-17-1-0578 and AFOSR FA9550-19-1-0291.

##  1 Introduction

Nowadays, it is widely accepted that the most serious threat to ecosystems is
the global warming fueled by the uncontrolled increase in carbon emissions,
and for the last twenty-some years, starting with the Kyoto Protocol in 1997,
international treaties have sprung out in hope to address this negative
externality. The most recent of these treaties is the _Paris Agreement_ with
$196$ signatories aiming at keeping the increase in temperature below
$2^{\circ}$C. According to the Environmental Protection Agency, electricity
production ($27\%$) and transportation ($28\%$) claim the lion share of the
total Greenhouse Gas emissions in the US.
222https://www.epa.gov/ghgemissions/global-greenhouse-gas-emissions-data
Throughout the world, local and federal governments try to disincentivize
reliance on polluting means of production by introducing carbon taxes or cap-
and-trade programs. In the latter case, regulators put a limit on the
allowable quantity of Green House Gas (GHG), any quantity above this limit
having to be covered by emission certificates (allowances) or the payment of a
penalty. In the former case, whether they are levied _upstream_ or
_downstream_ , carbon taxes aim at penalizing the use of fossil fuels for
their carbon content. The interested reader is referred to [14, 15, 13] for a
review of the state of affairs in the early days of the European Union
Emission Trading System, and mathematical treatments of thorough partial
equilibrium models for the comparison of realistic implementations of these
policies in the electricity sector.

As explained earlier, here we also concentrate on the electricity sector and
we propose a model for the analysis of the impact of investments in clean
means of production (e.g. solar and wind). While a model of the electricity
sector should comprise at least three types of agents: electricity producers,
resellers / retailers, and the end-users, we shall concentrate our modeling
effort on the producers. Until the challenging technological problem of
electricity storage is resolved at a larger scale, the demand for this
commodity remains inelastic, and we shall penalize the producers for not
matching the demand. In the following, we shall use the term _renewable_ to
mean electricity produced from wind turbines or solar panels. Alternatively,
we shall use the term _non-renewable_ to mean electricity produced by burning
fossil fuels like coal, crude oil or natural gas. We chose this convention for
convenience, even if this literary license is not completely accurate.

Individual producers control over time their usage of fossil fuels, and hence,
the amount of $CO_{2}$ emissions they are responsible for. They also control
their possible investment in solar or wind production, should they decide to
go that route. Notice that while the decision to use fossil fuels and how
much, changes over time, the investment in solar panels or a wind farm is a
one-time decision made at the beginning of the time period under
consideration. In our model, producing electricity from renewable sources
involves an initial investment and no extra cost over time since the marginal
cost of running these production assets is practically zero (except from
maintenance costs and possible subsidies which we ignore here). While the zero
cost of production is an attractive feature, it comes with the very high risks
due to the difficulties to predict the weather and the uncertainty associated
with the high volatility of these predictions. On the other hand, production
from traditional power plants is more predictable, the costs depending upon
the prices of the fuels and the price put on the $CO_{2}$ emissions by the
regulator. Each producer has to find the right balance between the pros and
the cons of the two major means of production we single out in our stylized
model. The overarching goal is to decarbonize to meet climate goals harnessing
demand-side policies through the establishment of a tax, as well as supply-
side resources including wind and solar production technologies.

Our economic model is based on the premises that the individual producers and
the regulator have only access to aggregate quantities. Basically, they only
have access to the statistical distributions of the productions, emissions,
investments, etc of the individual producers. As a result, we propose two
separate frameworks for the individual producers to optimize the mix of
renewable and nonrenewable production they should include in their portfolios.
We compute and compare the optimal centralized strategies by solving mean
field control problems, and the optimal decentralized strategies by solving
mean field game problems. In both cases, we monitor the effect of a carbon tax
on the optimal and equilibrium decisions. Quantifying the differences between
the two approaches is reminiscent of what is known as the Price of Anarchy
(PoA).

Among the conclusions drawn from the analysis of our model, we confirm that a
carbon tax is an effective incentive for the use of renewables. Also intuitive
is the fact that in the absence of a carbon tax, the overall pollution is
greater when producers compete than when they cooperate. Less obvious is the
fact that cooperating producers can have a lower level of global pollution
than when they compete even if the carbon tax is significant. We also show
that stricter regulations tend to reduce the differences between competitive
and cooperative equilibria. Further, we argue that the best way for the
regulator to enforce demand matching is to incentivize competition over
cooperation among the producers.

Mean Field Game (MFG) models appeared simultaneously and independently in the
original works of [9] and [23]. The thrust of these works was to propose a
paradigm to overcome the challenges of the search for Nash equilibria in large
games by considering models for which the interactions between the players
were of a mean field type, and deriving effective equations in the limit when
the number of players goes to infinity. Models in which a single player plays
a different role from the field of remaining players were introduced and
studied under the name of MFGs with major and minor players. In their
Stackelberg version, they had a significant impact on problems in economic
contract theory. See for example [7], [26], [20], or [16], [19] or [5]. Notice
that in these models, the major player uses a time dependent control, while we
shall assume that regulator uses time independent controls.

Using mean field models for energy applications is very natural. Competition
in the oil industry and the impact of the renewable energy competition was
analyzed in [21] and [17]. The early work [21] was extended with the addition
of a regulator in [1]. In [3], optimal entry and exit times for two types of
agents, electricity producers using either renewable or nonrenewable energy
resources, are analyzed using MFGs. Competition among electricity producers is
analyzed in [18] by using Mean Field Type Game where the mean field
interactions come through conditional expectation of the electricity price and
in [4] by using a model where the interactions enter the electricity spot
price. In [2] and [19], electricity consumers constitute the mean field
population and a single electricity producer plays the role of the principal,
in contrast to our model where we take the electricity producers as the mean
field population and the regulator as the principal.

Mean field models have also been used to model environmental impacts. In [6],
a MFG model is proposed to model climate change negotiations among countries
interacting through a $CO_{2}$ emission permit market. Emission certificate
markets are also studied in [27] and [28], again without the presence of a
regulator.

The paper is structured as follows. In Section 2, we introduce the minor
players’ model and the various equilibrium notions used in the sequel. In
Section 3 (resp. 4), the main theoretical (resp. numerical) results for the
minor players’ model are given. In Section 5, we introduce the regulator and
define the relevant notions of equilibrium. Finally, we provide numerical
results for the combined model with minor players and the regulator in Section
6 and we summarize our findings in a short Section 7.

##  2 Mean Field Model for Electricity Producers

###  2.1 N-Player Model

Although we will focus on mean field limits involving an infinite number of
players, we start with the description of what the $\mathcal{N}$-player
version of the game would be. For symmetry reasons, we assume that the total
electricity demand is split equally between all the agents, and each agent
faces the same demand, say $D_{t}$ at time $t$. The state of producer $i$ is
five-dimensional: instantaneous electricity production
$Q^{i}_{t}\in\mathbb{R}_{+}$, instantaneous irradiance
$S^{i}_{t}\in\mathbb{R}_{+}$, instantaneous emission level
$E^{i}_{t}\in\mathbb{R}_{+}$, cumulative pollution
$P^{i}_{t}\in\mathbb{R}_{+}$, and instantaneous nonrenewable energy production
$\tilde{N}^{i}_{t}\in\mathbb{R}_{+}$. Producer $i$ controls their state by
choosing at time $t=0$ an initial investment $R^{i}_{e}\in\mathbb{R}$ in
renewable production assets (e.g. for the purchase of solar panels), and at
each subsequent time $t$, by choosing the rate of change
$N^{i}_{t}\in\mathbb{R}$ in nonrenewable energy production. Notice that their
first control $N^{i}_{t}$ is time dependent while their second control
$R^{i}_{e}$ is time independent. This will be one of a challenging aspect of
our model.

######  Remark 2.1.

For the sake of definiteness, we use the terminology of solar power
production. However, other types of renewable energy can be modelled in a
similar way. For example, for wind power, $S^{i}_{t}$ would stand for the
instantaneous output of a wind farm and $R_{e}$ would be the corresponding
initial investment.

So with these proviso out of the way, define the time evolution of the state
of producer $i$ as:

$\displaystyle dQ^{i}_{t}$
$\displaystyle=\underbrace{\kappa_{1}N^{i}_{t}dt}_{\text{Term
1}}+\underbrace{\kappa_{2}R_{e}^{i}\left(\alpha\cos(\alpha
t)dt+{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}dS^{i}_{t}}\right)}_{\text{Term
2}}$
$\displaystyle{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color<EMAIL_ADDRESS>$\displaystyle=(\theta-S^{i}_{t})dt+\sigma_{0}d\widecheck{W}^{i}_{t}$
$\displaystyle dE^{i}_{t}=\delta N^{i}_{t}dt+\sigma_{1}dW^{i}_{t}$
$\displaystyle dP^{i}_{t}$ $\displaystyle=E^{i}_{t}dt$ $\displaystyle
d\tilde{N}^{i}_{t}=N^{i}_{t}dt$

The instantaneous electricity production changes depend on the instantaneous
nonrenewable energy usage (given by term 1) and the instantaneous yield from
the renewable energy investment (given by term 2). This second term includes a
seasonality component (sinusoidal term) and a random shock for the variability
of the sun irradiance. The form of the seasonality component was chosen for
the sake of simplicity. It can easily be extended to several harmonics to
include nightly and daily, monthly and yearly effects. In any case, we have
$Q_{t}^{i}=Q^{i}_{0}+\kappa_{1}\tilde{N}^{i}_{t}+\kappa_{2}R^{i}_{e}(\sin(\alpha
t)+S^{i}_{t})$ where $\kappa_{1},\kappa_{2}>0$ are constants that give the
efficiency of the production from nonrenewable and renewable energy,
respectively. The constant $\alpha>0$ gives the period of the seasonality of
the renewable energy.

We model the idiosyncratic noise terms $S^{i}_{t}$ in the renewable
productions as independent stationary processes. For the sake of definiteness,
we assume that they are Ornstein-Uhlenbeck processes with the same mean
$\theta>0$ and volatility $\sigma_{0}>0$, the $\widecheck{W}^{i}$ being
independent Wiener processes.

The dynamics of the instantaneous emissions $E^{i}_{t}$ have two components:
the contribution from the production from nonrenewable energy power plants,
and idiosyncratic random shocks with constant volatility $\sigma_{1}>0$ given
by independent Wiener processes $W^{i}$, also independent of the
$\widecheck{W}^{i}$’s. Despite its simple form, the choice of the constant
$\delta$ could include the effects of some abatement measures such as carbon
capture, sequestration and the use of filters.

Using the notation $\tilde{N}^{i}_{t}$ for the instantaneous nonrenewable
given by $\tilde{N}^{i}_{t}=\tilde{N}^{i}_{0}+\int_{0}^{t}N^{i}_{s}ds$, the
expected cost to producer $i$ over the whole period is:

$C^{\mathcal{N}}(N^{i},R^{i}_{e};\bar{Q})=\mathbb{E}\Big{[}\int_{0}^{T}\Big{[}\underbrace{c_{1}|N^{i}_{t}|^{2}}_{\text{Term
1}}+\underbrace{p_{1}\tilde{N}^{i}_{t}}_{\text{Term
2}}+\underbrace{c_{2}|Q^{i}_{t}-D_{t}|^{2}}_{\text{Term 3}}-\\\
\underbrace{c_{3}\big{(}\rho_{0}-\rho_{1}(D_{t}-\bar{Q}_{t})\big{)}Q^{i}_{t}}_{\text{Term
4}}\Big{]}dt+\underbrace{\tau|P^{i}_{T}|^{2}}_{\text{Term
5}}+\underbrace{p(R_{e}^{i})\Big{]}}_{\text{Term 6}}$ (2.1)

where $\bar{Q}=\sum_{j=1}^{\mathcal{N}}Q^{j}/\mathcal{N}$ and
$p:\mathbb{R}_{+}\mapsto\mathbb{R}_{+}$ is the investment in renewable energy
function.

Term 1 with $c_{1}>0$, is a penalty for attempting to ramp up and down
nonrenewable energy power plants too quickly. Term 2 represents the costs of
the fossil fuels used in nonrenewable power plants. The constant $p_{1}>0$ can
be understood as the average cost of one unit of fossil fuel. In lieu of
storage which is not included in our models because of its scarcity, Term 3
with $c_{3}>0$, imposes a penalty on producers for not matching the demand.
Term 4 represents the revenues from electricity production,
$\big{(}\rho_{0}+\rho_{1}(D_{t}-\bar{Q}_{t})\big{)}$ being the inverse demand
function which is assumed to be linear in excess demand or supply and
$\rho_{0},\rho_{1}>0$ are constants. It captures the fact that the price
increases if there is excess demand, and it decreases if there is excess
supply. We assume that the producers are selling what they produce. This term
introduces the mean field interactions into the model. Term 5 gives the carbon
tax levied by the regulator. We emphasize its role by assuming it is
proportional to the square of the terminal pollution. Term 6 is the initial
investment in renewable electricity production, for example the cost of solar
panels installed at the beginning of the period.

###  2.2 The Mean Field Model

In the mean field we focus on a representative producer interacting with the
field of the other producers, so we drop the superscript $i$ and the dynamics
equations become:

$\displaystyle dQ_{t}$
$\displaystyle=\kappa_{1}N_{t}dt+\kappa_{2}R_{e}\left(\alpha\cos(\alpha
t)dt+{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}(\theta-
S_{t})dt+\sigma_{0}d\widecheck{W}_{t}}\right)$
$\displaystyle{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color<EMAIL_ADDRESS>$\displaystyle=(\theta-S_{t})dt+\sigma_{0}d\widecheck{W}_{t}$ $\displaystyle
dE_{t}=\delta N_{t}dt+\sigma_{1}dW_{t}$ $\displaystyle dP_{t}$
$\displaystyle=E_{t}dt$ $\displaystyle d\tilde{N}_{t}=N_{t}dt$

where $W$ and $\widecheck{W}$ are independent Wiener processes Accordingly,
the total expected cost becomes:

$C(N,R_{e};\bar{Q})=\mathbb{E}\Big{[}\int_{0}^{T}\Big{[}c_{1}|N_{t}|^{2}+p_{1}\tilde{N}_{t}+c_{2}|Q_{t}-D_{t}|^{2}-\\\
c_{3}\big{(}\rho_{0}+\rho_{1}(D_{t}-\bar{Q}_{t})\big{)}Q_{t}\Big{]}dt+\tau|P_{T}|^{2}+p(R_{e})\Big{]}$
(2.2)

where $\bar{Q}_{t}=\mathbb{E}[Q_{t}]$. We shall sometimes use the notation
$\bar{Q}_{t}(N,R_{e})$ to emphasize the fact that the expectation is computed
under the state dynamics controlled by the admissible control $(N,R_{e})$.

###  2.3 Equilibrium Notions

We consider two different models: mean field game (MFG) and mean field control
(MFC). In the mean field game model, producers behave competitively and
minimize their total expected costs (search for their best responses) given
the other players’ decisions. A Nash equilibrium is then characterized as a
fixed point of the best response map so defined. In the sequel, we restrict
our attention to admissible strategies $(N,R_{e})$ such that
$\mathbb{E}[\int_{0}^{T}|N_{t}|^{2}dt]<+\infty$ and $R_{e}\in\mathbb{R}_{+}$.

######  Definition 2.2 (MFG Nash Equilibrium).

An admissible strategy $(\hat{N},\hat{R}_{e})$ is called an MFG equilibrium if
the associated production mean $\bar{Q}(\hat{N},\hat{R}_{e})$ is such that for
any admissible $(N,R_{e})$, we have

$C\Big{(}{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}(N,R_{e})};\bar{Q}(\hat{N},\hat{R}_{e})\Big{)}\geq
C\Big{(}(\hat{N},\hat{R}_{e});\bar{Q}(\hat{N},\hat{R}_{e})\Big{)}.$

In the mean field control case, we assume that the producers cooperate and
leave the choice of the control to a social planner minimizing the total
expected cost as defined in (2.2). In this case, if one player changes their
behavior, every player changes in the same way, and the mean field is
affected. The problem is now an optimal control problem.

######  Definition 2.3 (Social Planner’s MFC Optimum).

An admissible strategy $(N^{*},R^{*}_{e})$ is called an MFC optimum if the
associated production mean $\bar{Q}(\hat{N},\hat{R}_{e})$ is such that for any
admissible $(N,R_{e})$, we have

$C\Big{(}{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}(N,R_{e})};\bar{Q}({\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}N,R_{e}})\Big{)}\geq
C\Big{(}(N^{*},R^{*}_{e});\bar{Q}(N^{*},R^{*}_{e})\Big{)}.$

##  3 Theoretical Results

The presence of both time dependent and time independent controls creates
challenges in the regular probabilistic approach. See Remark 3.5. For this
reason, we first notice that:

$\inf_{(N_{t})_{t},R_{e}}C(N,R_{e};\bar{Q})=\inf_{R_{e}}\inf_{(N_{t})_{t},}C(N,R_{e};\bar{Q}),$

and we assume that $R_{e}$ is fixed in a first analysis. Next, we rewrite the
model in matrix form using $X_{t}:=[Q_{t}\quad S_{t}\quad E_{t}\quad
P_{t}\quad\tilde{N}_{t}]^{\top}$ as $5$-dimensional state process at time $t$.

$\inf_{(N_{t})_{t}}\tilde{C}\Big{(}N;R_{e},\bar{X}\Big{)}=\inf_{(N_{t})_{t}}\mathbb{E}\Bigg{[}\int_{0}^{T}\Bigg{[}\frac{R}{2}|N_{t}|^{2}+H^{\top}_{t}X_{t}+\bar{X}_{t}^{\top}FX_{t}+X_{t}^{\top}GX_{t}+J_{t}\Bigg{]}dt+X^{\top}_{T}S_{T}X_{T}+p(R_{e})\Bigg{]}$
(3.1) $dX_{t}=\Big{(}AX_{t}+B\cdot N_{t}+C_{t}\Big{)}dt+\Sigma
d\widetilde{W}_{t}$

where $R$ and $J_{t}$ are the scalars given by $R=2c_{1}$ and
$J_{t}=c_{2}D_{t}^{2}$ and:

$H_{t}=\begin{bmatrix}-(2c_{2}+c_{3}\rho_{1})D_{t}-c_{3}\rho_{0}\\\ 0\\\ 0\\\
0\\\ p_{1}\end{bmatrix}F=\begin{bmatrix}c_{3}\rho_{1}&0&0&0&0\\\ 0&0&0&0&0\\\
0&0&0&0&0\\\ 0&0&0&0&0\\\
0&0&0&0&0\end{bmatrix}G=\begin{bmatrix}c_{2}&0&0&0&0\\\ 0&0&0&0&0\\\
0&0&0&0&0\\\ 0&0&0&0&0\\\
0&0&0&0&0\end{bmatrix}S_{T}=\begin{bmatrix}0&0&0&0&0\\\ 0&0&0&0&0\\\
0&0&0&0&0\\\ 0&0&0&\tau&0\\\ 0&0&0&0&0\end{bmatrix}$
$A=\begin{bmatrix}0&-\kappa_{2}R_{e}&0&0&0\\\ 0&-1&0&0&0\\\ 0&0&0&0&0\\\
0&0&1&0&0\\\ 0&0&0&0&0\end{bmatrix}B=\begin{bmatrix}\kappa_{1}\\\ 0\\\
\delta\\\ 0\\\ 1\end{bmatrix}C_{t}=\begin{bmatrix}\kappa_{2}R_{e}\Big{(}\alpha
cos(\alpha t)+\theta\Big{)}\\\ \theta\\\ 0\\\ 0\\\
0\end{bmatrix}\Sigma=\begin{bmatrix}\kappa_{2}R_{e}\sigma_{0}&0\\\
\sigma_{0}&0\\\ 0&\sigma_{1}\\\ 0&0\\\
0&0\end{bmatrix}\widetilde{W}_{t}=\begin{bmatrix}\widecheck{W}_{t}\\\
W_{t}\end{bmatrix}.$

Furthermore, we define:

$a=\frac{1}{2}\Sigma\Sigma^{\top}=\frac{1}{2}\begin{bmatrix}(\kappa_{2}R_{e}\sigma_{0})^{2}&\kappa_{2}R_{e}\sigma_{0}^{2}&0&0&0\\\
\kappa_{2}R_{e}\sigma_{0}^{2}&\sigma_{0}^{2}&0&0&0\\\
0&0&\sigma_{1}^{2}&0&0\\\ 0&0&0&0&0\\\ 0&0&0&0&0\end{bmatrix}$

and the value function $u(t,X)$ as:

$u(t,X)=\inf_{(N_{t})_{t}}\mathbb{E}\Bigg{[}\int_{t}^{T}\Big{[}\frac{R}{2}|N_{s}|^{2}+H^{\top}_{s}X_{s}+\bar{X}_{s}^{\top}FX_{s}+X_{s}^{\top}GX_{s}+J_{s}\Big{]}ds+X^{\top}_{T}S_{T}X_{T}+p(R_{e})\Big{|}X_{t}=X\Bigg{]}$
(3.2)

######  Lemma 3.1 (ODE System for the MFG).

Given $R_{e}$, if there exists a function
$t\mapsto(\eta_{t},r_{t},\bar{X}_{t})$ solving the following system of
Ordinary Differential Equations (ODEs):

$\displaystyle\dfrac{d{\eta_{t}}}{dt}-{\eta_{t}}BR^{-1}B^{\top}{\eta_{t}}+A^{\top}{\eta_{t}}+{\eta_{t}}A+2G=0,$
$\displaystyle{\eta_{T}}=2S_{T}$ (3.3a)
$\displaystyle-\dfrac{d{r_{t}}}{dt}=\left(A^{\top}-{\eta_{t}}BR^{-1}B^{\top}\right){r_{t}}+{\eta_{t}}C_{t}+H_{t}+F^{\top}\bar{X}_{t},$
$\displaystyle{r_{T}}=0$ (3.3b)
$\displaystyle\dfrac{d\bar{X}_{t}}{dt}=(A-BR^{-1}B^{\top}{\eta_{t}})\bar{X}_{t}-BR^{-1}B^{\top}{r_{t}}+C_{t},$
$\displaystyle\bar{X}_{0}=\bar{x}_{0}$ (3.3c)

and $s_{0}$ is given by:

$s_{0}=p_{2}R_{e}+\int_{0}^{T}\Big{(}tr(a{\eta_{t}})-\frac{1}{2}{r_{t}}^{T}BR^{-1}B^{\top}{r_{t}}+C_{t}^{\top}{r_{t}}+J_{t}\Big{)}dt,$
(3.4)

then $\hat{N}_{t}(R_{e})=-R^{-1}B^{\top}(\eta_{t}X_{t}+r_{t})$ together with
the constant $R_{e}$ is the MFG equilibrium given $R_{e}$ fixed, and the
expected cost to the representative producer in this equilibrium is:

$\inf_{(N_{t})_{t}}\tilde{C}^{MFG}\Big{(}N;R_{e},\bar{X}\Big{)}=\frac{1}{2}\left(Var(\sqrt{\eta_{0}}X_{0})+\mathbb{E}[\sqrt{\eta_{0}}X_{0}]^{2}\right)+\bar{X}_{0}^{\top}r_{0}+s_{0}.$
(3.5)

######  Theorem 3.2.

Given fixed $R_{e}$, if $T$ is small enough, there exists a unique MFG
equilibrium.

######  Lemma 3.3 (MFC ODE System).

Given $R_{e}$, if there exists a function
$t\mapsto(\eta_{t},r_{t},\bar{X}_{t})$ solving the ODE system (3.3) with the
equation (3.3b) replaced by:

$-\frac{d{r_{t}}}{dt}=\left(A^{\top}-{\eta_{t}}BR^{-1}B^{\top}\right){r_{t}}+{\eta_{t}}C_{t}+H_{t}+F^{\top}{\bar{X}_{t}}+{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}F{\bar{X}_{t}}},\qquad{r_{T}}=0$
(3.6)

and the same $s_{0}$ given by (3.4), then
$N^{*}_{t}(R_{e})=-R^{-1}B^{\top}(\eta_{t}X_{t}+r_{t})$ together with the
constant $R_{e}$ is an optimum for the MFC problem given $R_{e}$, and the
minimal expected cost is

$\inf_{(N_{t})_{t}}\tilde{C}^{MFC}\Big{(}N;R_{e}\Big{)}=\frac{1}{2}\left(Var(\sqrt{\eta_{0}}X_{0})+\mathbb{E}[\sqrt{\eta_{0}}X_{0}]^{2}\right)+\bar{X}_{0}^{\top}r_{0}+s_{0}{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color<EMAIL_ADDRESS>(3.7)

######  Theorem 3.4.

Given fixed $R_{e}$, if $T$ is small enough, there exists a unique MFC
optimum.

######  Remark 3.5.

When $R_{e}$ is treated as a control of the representative producer, both
solutions of the mean field game and the central planner control problem
become more difficult. We illustrate this fact with the examination of what
the probabilistic approach introduced in [10] for MFGs and in [11] for MFC
problems (see also [12] for a presentation in textbook form) would give in
this situation. This approach provides a solution $(\hat{N}_{t},\hat{R}_{e})$
if one can find processes $({X}_{t},{Y}_{t},{Z}_{t})_{t}$ solving the
following Forward-Backward Stochastic Differential Equation (FBSDE):

$\displaystyle dQ_{t}$
$\displaystyle=\kappa_{1}\hat{N}_{t}dt+\kappa_{2}\hat{R}_{e}\left(\alpha\cos(\alpha
t)dt+(\theta-S_{t})dt+\sigma_{0}d\widecheck{W}_{t}\right),\qquad$
$\displaystyle Q_{0}$ $\displaystyle=q_{0}$ $\displaystyle dS_{t}$
$\displaystyle=(\theta-S_{t})dt+\sigma_{0}d\widecheck{W}_{t},\qquad$
$\displaystyle S_{0}$ $\displaystyle=\theta$ $\displaystyle dE_{t}$
$\displaystyle=\delta\hat{N}_{t}dt+\sigma_{1}dW_{t},\qquad$ $\displaystyle
E_{0}$ $\displaystyle=e_{0}$ $\displaystyle dP_{t}$
$\displaystyle=E_{t}dt,\qquad$ $\displaystyle P_{0}$ $\displaystyle=p_{0}$
$\displaystyle d\tilde{N}_{t}$ $\displaystyle=\hat{N}_{t}dt,\qquad$
$\displaystyle\tilde{N}_{0}$ $\displaystyle=\tilde{n}_{0}$ $\displaystyle
dY^{1}_{t}$
$\displaystyle=\Big{(}-{2c_{2}({Q}_{t}-D_{t})}+{c_{3}\big{(}\rho_{0}+\rho_{1}(D_{t}-\bar{Q}_{t})\big{)}}\Big{)}dt+Z_{t}^{1,1}d\widecheck{W}_{t}+Z_{t}^{1,2}dW_{t}\qquad$
$\displaystyle Y_{T}^{1}$ $\displaystyle=0$ $\displaystyle dY^{2}_{t}$
$\displaystyle=\Big{(}\kappa_{2}\widehat{R}_{e}Y_{t}^{1}+Y_{t}^{2}\Big{)}dt+Z_{t}^{2,1}d\widecheck{W}_{t}+Z_{t}^{2,2}dW_{t},\qquad$
$\displaystyle Y_{T}^{2}$ $\displaystyle=0$ $\displaystyle dY^{3}_{t}$
$\displaystyle=Y_{t}^{4}dt+Z_{t}^{3,1}d\widecheck{W}_{t}+Z_{t}^{3,2}dW_{t},$
$\displaystyle Y_{T}^{3}$ $\displaystyle=0$ $\displaystyle dY^{4}_{t}$
$\displaystyle=Z_{t}^{4,1}d\widecheck{W}_{t}+Z_{t}^{4,2}dW_{t},$
$\displaystyle Y_{T}^{4}$ $\displaystyle=2\tau\widehat{P}_{T}$ $\displaystyle
dY^{5}_{t}$
$\displaystyle=-p_{1}dt+Z_{t}^{5,1}d\widecheck{W}_{t}+Z_{t}^{5,2}dW_{t},$
$\displaystyle Y_{T}^{5}$ $\displaystyle=0,$

in which case $(\hat{N},\hat{R}_{e})$ given by:

$\hat{N}_{t}=-\dfrac{Y_{t}^{1}\kappa_{1}+Y_{t}^{3}\delta+Y_{t}^{5}}{2c_{1}}\quad\text{and}\quad\hat{R}_{e}=(p^{\prime})^{-1}\Big{(}-\mathbb{E}\Big{[}\int_{0}^{T}\kappa_{2}Y_{t}^{1}\Big{(}\alpha\cos(\alpha
t)+(\theta{-S_{t}})\Big{)}dt\Big{]}\Big{)}$ (3.8)

is a Nash equilibrium for the MFG. Here, $(p^{\prime})^{-1}$ is the inverse of
the first derivative of the function $p$ (which we assume to exist here),
$X_{t}=[Q_{t},S_{t},E_{t},P_{t},\tilde{N}_{t}]$,
$Y_{t}=[Y^{1}_{t},Y^{2}_{t},Y^{3}_{t},Y^{4}_{t},Y^{5}_{t}]$ and

$Z_{t}=\begin{bmatrix}Z_{t}^{1,1}\quad&Z_{t}^{2,1}\quad&Z_{t}^{3,1}\quad&Z_{t}^{4,1}\quad&Z_{t}^{5,1}\\\
Z_{t}^{1,2}\quad&Z_{t}^{2,2}\quad&Z_{t}^{3,2}\quad&Z_{t}^{4,2}\quad&Z_{t}^{5,2}\\\
\end{bmatrix}^{\top}.$

Moreover, if the equation for $(Y^{1}_{t})_{t}$ is replaced by

$dY^{1}_{t}=\Big{(}-{2c_{2}({Q}_{t}-D_{t})}+{c_{3}\big{(}\rho_{0}+\rho_{1}(D_{t}-{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}2\bar{Q}_{t}})\big{)}}\Big{)}dt+Z_{t}^{1,1}d\widecheck{W}_{t}+Z_{t}^{1,2}dW_{t},$

then $(\hat{N},\hat{R}_{e})$ is an MFC optimum. Even though the set-up of our
models does not seem to fit in the traditional MFGs and MFC problems, the
above claims can still be proved deriving appropriate versions of the
stochastic maximum principle. We do not give the details here. A first reason
is because the strategy of the proof is relatively standard. But the main
reason is that the necessary form (3.8) of the optimal controls, when injected
back into the FBSDE, turns the deceivingly simple looking linear FBSDE into a
path dependent nonlinear FBSDE of the McKean-Vlasov type which we could not
solve analytically.

##  4 Numerical Results

For numerical purposes, given the challenges posed by the solution of the
large FBSDE given above (recall Remark 3.5), we use the analytic approach
based on the ODEs (3.3) and (3.6). We search for the $R_{e}$ and the
corresponding equilibrium $N_{t}$ that minimizes the cost of the minor
players.

As emphasized earlier, the main difference between MFC and MFG is whether the
mean field is affected by the decision of the representative producer (MFC),
or taken to be fixed (MFG). This difference translates into the addition of a
fixed point argument in the MFG case. For pedagogical reasons, we first
discuss the MFC case, then the MFG.

After solving the Riccati equation which is the same in both cases, we solve
the coupled ODE system directly in the MFC case in order to find the mean
field; on the other hand, notice that in the MFG case, the ODEs are decoupled
since the mean field is assumed to be fixed in each iteration of the fixed
point algorithm.

###  4.1 Mean Field Control Algorithm

In order to solve the system of MFC coupled ODEs for $(\bar{X}_{t})_{t}$ and
$r_{t}$ given by equations (3.3c) and (3.6), we discretize the time with
uniform step size $\Delta t$ and solve the following linear equation:

$\begin{bmatrix}\bar{X}\\\ r\end{bmatrix}=M\begin{bmatrix}\bar{X}\\\
r\end{bmatrix}+K$ (4.1)

where $\bar{X}=[\bar{X}_{0},\bar{X}_{\Delta t},\bar{X}_{2\Delta
t},\dots,\bar{X}_{T}]^{\top}$, $r=[r_{0},r_{\Delta t},r_{2\Delta
t},\dots,r_{T}]^{\top}$.

Algorithm 1 Algorithm for Computing the Mean Field Control Cost over
$(N_{t})_{t}$ given $R_{e}$

1:function Optim-MFC-N($R_{e}$)

2: Calculate $(\eta_{t})_{t}$ by solving the Riccati Equation in (3.3a)

3: Solve the coupled $(\bar{X}_{t})_{t}$ and $(r_{t})_{t}$ linear system (4.1)
given $(\eta_{t})_{t}$ and $R_{e}$

4: Calculate $s$ given $R_{e}$, $(r_{t})_{t}$ and $(\eta_{t})_{t}$ using the
equation in (3.4)

5: Calculate the expected cost associated with $R_{e}$,
$\hat{c}:=\inf\tilde{C}^{MFC}(N;R_{e},\bar{X})$ using (3.7)

6: return ($\hat{c},\bar{X}$)

Algorithm 2 Algorithm for Finding a Social Optimum

1:function SocialOpt

2: Search for the optimal $\hat{R}_{e}$ where the optimal cost
$R_{e}\rightarrow c(R_{e})$ and optimal mean field
$R_{e}\rightarrow\bar{X}(R_{e})$ are computed by Optim-MFC-N

3: Let $\hat{c}=c(\hat{R}_{e})$ and $\hat{\bar{X}}=\bar{X}(\hat{R}_{e})$

4: return $(\hat{c},\hat{R}_{e},\hat{\bar{X}})$

###  4.2 Mean Field Game Algorithm

In Mean Field Game, since in each iteration it is assumed that the
$(\bar{X}_{t})_{t}$ is fixed, the ODE for $(r_{t})_{t}$ in equation (3.3b) can
be solved directly by using the following linear equation after we discretize
time:

$r=M_{r}r+K_{r}$ (4.2)

where $r=[r_{0},r_{\Delta t},r_{2\Delta t},\dots,r_{T}]^{\top}$. Then with
this $(r_{t})_{t}$, the time discretization of $(\bar{X}_{t})_{t}$ with
dynamics given by the equation (3.3c) can be written as:

$\bar{X}=M_{\bar{X}}\bar{X}+K_{\bar{X}}$ (4.3)

where $\bar{X}=[\bar{X}_{0},\bar{X}_{\Delta t},\bar{X}_{2\Delta
t},\dots,\bar{X}_{T}]^{\top}$. The numerical algorithms to find the Mean Field
Control and Game Equilibria are given in detail in the following sections.

Algorithm 3 Algorithm for Computing Expected Cost over $(N_{t})_{t}$ given
$R_{e}$, $(\bar{X}_{t})_{t}$

1:function Optim-MFG-N($R_{e},(\bar{X}_{t})_{t}$)

2: Calculate $(\eta_{t})_{t}$ by solving Riccati Equation in (3.3a)

3: Solve the linear system for $(r_{t})_{t}$ in (4.2) given
$(\bar{X}_{t})_{t}$, $R_{e}$ and $(\eta_{t})_{t}$

4: Calculate $s$ given $R_{e}$, $(r_{t})_{t}$ and $(\eta_{t})_{t}$ using the
equation in (3.4)

5: Calculate the cost associated with $R_{e}$ and $(\bar{X}_{t})_{t}$,
$\hat{c}:=\inf\tilde{C}^{MFG}(N;R_{e},\bar{X})$ using (3.5)

6: return $\hat{c}$

Algorithm 4 Algorithm for Finding Nash Equilibrium

1:function NashEq

2: Initialize $(\bar{X}^{0}_{t})_{t}$

3: while $||\bar{X}^{k}-\bar{X}^{k-1}||>\epsilon$ do

4: Search for the optimal $\hat{R}_{e}$ given $\bar{X}^{k}$ where the optimal
cost
$(R_{e},\bar{X}^{{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}k}})\rightarrow
c^{{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}k}}(R_{e},\bar{X}^{{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}k}})$
is computed by Optim-MFG-N

5: Let
$\hat{R}_{e}^{{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}k}}=\operatorname*{arg\,min}_{R_{e}}c^{{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}k}}(R_{e},\bar{X}^{{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}k}})$

6: Compute
$(\bar{X}^{{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}{k+1}}}_{t})_{t}$
given
$\hat{R}_{e}^{k},(\bar{X_{t}}^{{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}k}})_{t}$
by solving the linear equation (4.3)

7: Let
$\hat{R}_{e}=\hat{R}^{{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color<EMAIL_ADDRESS>$\hat{c}=c^{{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}k}}(\hat{R}_{e})$
and
$\hat{\bar{X}}=\bar{X}^{{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color<EMAIL_ADDRESS>
8: return $(\hat{c},\hat{R}_{e},\hat{\bar{X}})$

###  4.3 Numerical Experiments

In the numerical experiments reported below, we use the following parameter
values:

$p_{1}=7/\Delta_{t}$ (dollar/time) |  $\kappa_{1}=0.13$ (MWh/$10^{3}$ cu ft) | $\theta=5$
---|---|---
$p_{2}=10^{4}$ |  $\kappa_{2}=0.1$ (MWh/$10^{3}$ dollars) |  $T=20$ years, $\Delta_{t}=10$ days
$c_{1}=10^{-4}$ (dollar/$10^{3}$ cu ft2) | $\alpha=40\pi$ |  $R_{e}=[0,10^{4}]$ (10,000 dollars)
$c_{3}=1$ (dollar/$10^{3}$ cu ft) | $\delta=0.15$ | $D_{t}=2\times 10^{4}-5\times 10^{2}\cos(80\pi\Delta_{t})$
$\rho_{0}=40/\Delta_{t}$ | $\sigma_{0}=0.01$ | $\bar{X}_{0}=[0,\theta,0,0,0]$
$\rho_{1}=0.1/\Delta_{t}$ | $\sigma_{1}=0.01$ | $Var[X_{0}]=[0,0.1,0,0,0]$

Furthermore, we assume that $p(R_{e})=p_{2}R_{e}$ where $p_{2}$ is a positive
constant.

We focus on natural gas as the source of nonrenewable energy. In our numerical
experiments, we ignore the effect of the COVID-19 pandemic, and we run
simulations for $20$ years starting from March 2020. For the cost of solar
power, we use the current assumption that a $1$MW solar farm needs a $1$M$
investment, 333https://news.energysage.com/solar-farms-start-one/and use daily
peak sun hours data to compute daily average production from solar panels. We
assume that on average, peak hours are around 5 hours in the USA, and we infer
that a solar farm built with an initial investment of $\$10,000$ generates on
average $0.5$MWh in $10$ days. We choose $\alpha$ to take into account the
seasonality, since the sun exposure levels are maximum during summer, minimum
during winter. Therefore, we infer that 1 unit $R_{e}$ corresponds to $10,000
investment and on average it generates $\kappa_{2}(\theta\pm 1)=0.1(5\pm
1)$MWh electricity in $10$ days.

According to the data provided by U.S. Energy Information Administration
(EIA), in 2018, $1,365,822$ Million KWh were produced by the electric sector
by using 10,215 billion cubic feet of natural gas.
(444https://www.eia.gov/totalenergy/data/monthly/pdf/sec7.pdfTable 7.2b &
7.3b). Therefore, we assume that $1000$ cubic feet of natural gas produces
$\sim 0.13$MWh. Again, according to the data provided by EIA, typical natural
gas power plants produce $0.91$ pound carbon emission per kWh electricity
generation, 555https://www.eia.gov/tools/faqs/faq.php?id=74&t=11Since $1000$
cubic feet of natural gas produce around $0.13$ MWh, we conclude that around
$122$ pounds of carbon is emitted with the use of natural gas. Moreover, if
coal or other fossil fuels were used, this carbon pollution would have to be
much more than doubled. Taking this into account and converting to metric tone
we end up with $\delta\approx 0.15$.

We assume that the average demand of electricity for each plant is around the
capacity of the plants. According to EIA
(666https://www.eia.gov/electricity/annual/archive/pdf/epa_2018.pdfTable 4.3),
the average daily capacity of a natural gas plant is around $2166.5$ MWh in
2018. Furthermore, the monthly seasonal component is found by using the
monthly residential electricity consumption in 2018 given in 33footnotemark:
3Table 7.6. Therefore, $10$ day demand is taken sinusoidal to show the
seasonality around $20,000$MWh. According to the data provided by EIA in
777https://www.eia.gov/electricity/annual/html/epa_08_04.htmland
888https://www.eia.gov/dnav/ng/ng_pri_sum_a_EPG0_PEU_DMcf_a.htmnonrenewable
energy has 40% of its fuel cost as the operation and maintenance costs on top
of the fuel cost in 2018 and the price of 1000 cu ft natural gas can be
assumed $5. Therefore, we take $p_{1}=\$7$. Finally by using
999https://www.eia.gov/todayinenergy/detail.php?id=37912we see that the
average price of wholesale electricity is around $40 per MWh, therefore we
take $\rho_{0}=40$.

####  4.3.1 Price of Anarchy (PoA) Analysis

From the heat maps in Figure 1 we see that the expected cost of the
representative producer is increasing with the carbon tax $\tau$ and the
penalty $c_{2}$. The second observation is that as expected, for any given
couple $(\tau,c_{2})$, the expected cost is higher in the Nash equilibrium
than for the Social Optimum. We quantify how inefficient the Nash equilibrium
is, and the effect of $\tau$ and $c_{2}$ on this inefficiency, in other words,
we quantify the adverse effect of the non-cooperative behavior of the
producers, by computing the Price of Anarchy (PoA) defined in (4.4) for
different values of $\tau$ and $c_{2}$.

$PoA(\tau,c_{2})=\frac{\inf_{N_{t},R_{e}}C^{MFG}(N_{t},R_{e};\bar{Q},\tau,c_{2})}{\inf_{N_{t},R_{e}}C^{MFC}(N_{t},R_{e};\bar{Q},\tau,c_{2})}$
(4.4)

The results are given in right subfigure in Figure 1. Since for any given
$(\tau,c_{2})$ the expected cost in a MFG equilibrium is higher, PoA is
expected to be greater than $1$ and as it gets higher, the Nash Equilibrium is
getting less efficient. It can be seen that PoA gets smaller as we increase
$\tau$ and $c_{2}$. This means that for higher levels of $\tau$ and $c_{2}$,
the expected costs of to the producers become closer. In other words, the
impact of the social planner diminishes and the advantages of cooperation
lessen as the regulator imposes stricter regulations.

Figure 1: Expected Cost of Minor Players in MFC (left), in MFG (middle) and
the Price of Anarchy given different penalty for not matching the demand and
tax levels

####  4.3.2 Electricity Production Decomposition Analysis

Figure 2: Production decompostion when there is no tax ($\tau$) and the
penalty for not matching the demand ($c_{2}$) is low (left); when there is no
tax and $c_{2}$ is high (middle); when tax and $c_{2}$ are high (right).

Here, we analyze the effect of the penalty $c_{2}$ for not matching the demand
and the carbon tax $\tau$, on the optimal energy production portfolio in both
MFC and MFG models. Figure 2 shows the total production and the decomposition
of this production over a $20$ year period together with a detailed zoom in
behavior between years $1$ and $3$.

The left subfigure in Figure 2, shows that the demand is not matched in the
MFC case. This is because the penalty $c_{2}$ for not matching the demand is
low and the increased revenue from scarce supply is more advantageous. Here,
we see that in the control setting, producers behave as a big monopoly when
not matching the demand is inexpensive. When the penalty is increased the
middle subfigure in Figure 2 shows that producers try to match the demand and
their behaviors in the MFC and MFG cases are similar. In both of these figures
there is no carbon tax, therefore the producers do not have incentives to
invest in renewable energy, and as a result, all the production is exclusively
from the nonrenewable sources. On the right subfigure in Figure 2 when the
carbon tax is increased we see that the producers have an incentive to invest
in renewable energy.

We also analyze the effect of the planning horizon where we compare the cases
in which the producers are planning for the next $2$ years vs. planning for
the next $20$ years. As it can be seen in the left and middle subfigures in
Figure 3, when the planning horizon is short, the fixed costs of renewable
energy outweigh its advantages. Short-sighted producers do not have an
incentive to invest in renewable energy production.

Figure 3: Planning time horizon effect in MFC (left), and in MFG (middle); End
of the 20 years pollution in MFC and MFG (right)

####  4.3.3 Pollution Analysis

The right subfigure in Figure 3, shows that whatever the level of the carbon
tax, the terminal pollution levels are higher when the producers are
competitive (MFG). Further, in the absence of a carbon tax, producers can
decrease the pollution levels further by cooperating and following a social
planner instead of implementing a carbon tax.

##  5 Models with a Regulator

In this section, we describe how the previous models can be extended to
include a major player in charge of choosing the tax level $\tau$, and the
penalty $c_{2}$ for not matching the demand. This major player represents a
_regulator_ , and we shall often speak of minor players when we talk about the
producers. We also modify the “minor player only” model used previously by
offering the producers the option to withdraw their entire production, de
facto walking away from the contract imposed by the regulator. This decision
is made when the expected cost to the producer is higher than a fixed level
above which producing at such a level of loss does not make sense.

In the new model, the regulator does not have a private state per se. It only
has $2$ controls which are the carbon tax level ($\tau\in\mathbb{R}_{+}$) and
the penalty $(c_{2}\in\mathbb{R}_{+})$. Both controls are assumed to be time
independent. This assumption is especially realistic when the period $[0,T]$
is too short for changes in regulation to make sense. The cost function of the
regulator is given as:

$J\Big{(}\tau,c_{2};\bar{P}_{T},\bar{Q}\Big{)}=\underbrace{\alpha_{1}\big{(}\bar{P}_{T}-\bar{P}^{*}_{T}\big{)}_{+}}_{\text{Term
1}}\underbrace{-\alpha_{2}\tau\big{(}\bar{P}_{T}-\bar{P}_{0}\big{)}}_{\text{Term
2}}\underbrace{+\alpha_{3}\big{|}\tau\big{|}^{2}}_{\text{Term
3}}\underbrace{+\int_{0}^{T}\alpha_{4}\big{|}\bar{Q}_{t}-D_{t}\big{|}^{2}dt}_{\text{Term
4}}\underbrace{+\alpha_{5}c_{2}^{2}}_{\text{Term 5}}.$ (5.1)

The first term is the cost for exceeding the pollution target
$\bar{P}_{T}^{*}\geq\bar{P}_{0}$ announced at $t=0$. Since we use the notation
$x_{+}=\max[0,x]$, there is no penalty if the terminal pollution level is
below the target. The constant $\alpha_{1}>0$ quantify the size of the
penalty. The second term is the revenue from the carbon tax. To prevent the
regulator from choosing an abusive high tax to increase their revenue, Term 3
is added to represent a reputation cost ($\alpha_{2},\alpha_{3}>0$). The joint
roles of Term 4 and Term 5 is to insure that matching the demand is not only
incumbent on the producers, but also on the regulator by choosing
$\alpha_{4}>0$. So the regulator also has the role of a system operator who
can be penalized for _system black-outs_ , sharing their costs with the
producers.

###  5.1 Equilibrium Notions

We analyze two types of equilibria in the models with a regulator. In both
cases, we consider that the regulator announces their policy first, and the
producers react accordingly. This is in the realm of Stackelberg games. We
call the first equilibrium Stackelberg MFC equilibrium. In this case, the
regulator assumes that a social planner chooses the controls used by the
electricity producers. The latter behave like one big monopolistic firm.
Therefore, the regulator chooses the tax level, $\tau$ and penalty coefficient
$c_{2}$, assuming that the producers will settle in a MFC optimum. Note that
in this interpretation the regulator and the social planner are two different
entities. We define this equilibrium formally as:

######  Definition 5.1 (Stackelberg MFC equilibrium).

For every $(\tau,c_{2})$, let
$\Big{(}N^{*}(\tau,c_{2}),R^{*}_{e}(\tau,c_{2})\Big{)}$ be the social
planner’s MFC optimum given the tax level $\tau$ and the penalty coefficient
$c_{2}$. In other words, for every $\tau,c_{2}$ and any admissible
$\big{(}N,R_{e}\big{)}$, we have:

$C\Big{(}{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}\big{(}N,R_{e}\big{)}};\bar{X}\big{(}{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}N,R_{e}}\big{)},(\tau,c_{2})\Big{)}\geq
C\Big{(}\big{(}N^{*}(\tau,c_{2}),R^{*}_{e}(\tau,c_{2})\big{)};\bar{X}\big{(}N^{*}(\tau,c_{2}),R^{*}_{e}(\tau,c_{2})\big{)},(\tau,c_{2})\Big{)}$

where we added the notation $\bar{X}\big{(}N,R_{e}\big{)}$ to emphasize the
parameters for which the mean field term $\bar{X}$ is computed. Then the
strategy profile $(\tau^{*},c_{2}^{*})$ is a Stackelberg MFC equilibrium with
a regulator if, for any admissible $(\tau,c_{2})$:

$J\Big{(}{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}(\tau,c_{2})};\bar{X}\big{(}N^{*}({\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}\tau,c_{2}}),R^{*}_{e}({\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}\tau,c_{2}})\big{)}\Big{)}\geq
J\Big{(}(\tau^{*},c_{2}^{*});\bar{X}\big{(}N^{*}(\tau^{*},c^{*}_{4}),R^{*}_{e}(\tau^{*},c_{2}^{*})\big{)}\Big{)}.$

The second equilibrium is called Stackelberg MFG Equilibrium. In this one
regulator assumes that electricity producers are competitive and it chooses
$\tau$ and $c_{2}$ levels by assuming that the minor player population is at
Nash Equilibrium. We can define this Equilibrium formally as

######  Definition 5.2 (Stackelberg MFG Equilibrium).

For every $(\tau,c_{2})$, let
$\Big{(}\hat{N}(\tau,c_{2}),\hat{R}_{e}(\tau,c_{2})\Big{)}$ be the producers
MFG Nash equilibrium given the tax level $\tau$ and the demand satisfaction
coefficient $c_{2}$. In other words, for any admissible
$\big{(}N,R_{e}\big{)}$, we have:

$C\Big{(}{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}\big{(}N,R_{e}\big{)}};{\bar{X}}\big{(}\hat{N}(\tau,c_{2}),\hat{R}_{e}(\tau,c_{2})\big{)},(\tau,c_{2})\Big{)}\geq
C\Big{(}\big{(}\hat{N}(\tau,c_{2}),\hat{R}_{e}(\tau,c_{2})\big{)};{\bar{X}}\big{(}\hat{N}(\tau,c_{2}),\hat{R}_{e}(\tau,c_{2})\big{)},(\tau,c_{2})\Big{)}.$

Then the strategy profile $(\hat{\tau},\hat{c}_{2})$ is a Stackelberg MFG
equilibrium with a regulator if, for any admissible $(\tau,c_{2})$, we have:

$J\Big{(}{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}(\tau,c_{2})};{\bar{X}}\big{(}\hat{N}({\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}\tau,c_{2}}),\hat{R}_{e}({\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}\tau,c_{2}})\big{)}\Big{)}\geq
J\Big{(}(\hat{\tau},\hat{c}_{2});{\bar{X}}\big{(}\hat{N}(\hat{\tau},\hat{c}_{2}),\hat{R}_{e}(\hat{\tau},\hat{c}_{2})\big{)}\Big{)}.$

##  6 Numerical Results in the Presence of a Regulator

###  6.1 Algorithms

To implement the walk-away option of the producers, we modify the SocialOpt
and NashEq algorithms. This is done by simply adding an IF condition to these
algorithms that assigns Accept=1 if the cost of the minor player is lower than
the threshold and Accept=0 otherwise. After the algorithms for the producers
are modified and called ”ModifiedSocialOpt” and ”ModifiedNashEq” respectively,
we implement the Stackelberg Equilibrium Algorithm where we assume that if the
producers reject the contract (Accept=0), the regulator cost is equal to
infinity.

Algorithm 5 Algorithm for Finding Regulator’s Cost

1:function RegulatorCost($\tau,c_{2},\bar{X},\text{Accept}$)

2: if Accept = 1 then

3: Compute the regulator’s cost $J=J(\tau,c_{2};\bar{X})$ by using (5.1)

4: else

5: $J=\infty$

6: return $J$

Algorithm 6 Algorithm for Finding Stackelberg Equilibrium with MFC and MFG

1:function StackelbergEq(Type)

2: if Type = MFC then

3: Search for optimal $(\hat{\tau},\hat{c}_{2})$ couple where optimal mean
field $(\tau,c_{2})\rightarrow\bar{X}(\tau,c_{2})$, investment in renewable
$(\tau,c_{2})\rightarrow R_{e}(\tau,c_{2})$ and minor cost
$(\tau,c_{2})\rightarrow c(\tau,c_{2})$ are computed by ModifiedSocialOpt
algorithm and optimal cost of regulator $(\tau,c_{2})\rightarrow
J(\tau,c_{2};\bar{X}(\tau,c_{2}),\text{Accept})$ is found by using
RegulatorCost

4: else if Type = MFG then

5: Search for optimal $(\hat{\tau},\hat{c}_{2})$ couple where optimal mean
field $(\tau,c_{2})\rightarrow\bar{X}(\tau,c_{2})$, investment in renewable
$(\tau,c_{2})\rightarrow R_{e}(\tau,c_{2})$ and minor cost
$(\tau,c_{2})\rightarrow c(\tau,c_{2})$ are computed by ModifiedNashEq
algorithm and optimal cost of regulator $(\tau,c_{2})\rightarrow
J(\tau,c_{2};\bar{X}(\tau,c_{2}),\text{Accept})$ is found by using
RegulatorCost

6: Let
$(\hat{\tau},\hat{c}_{2})=\operatorname*{arg\,min}_{\tau,c_{2}}J(\tau,c_{2};\bar{X}(\tau,c_{2}))$,
$\hat{\bar{X}}=\bar{X}(\hat{\tau},\hat{c}_{2})$,
$\hat{R_{e}}=R_{e}(\hat{\tau},\hat{c}_{2})$,
$\hat{c}=c(\hat{\tau},\hat{c}_{2})$ and
$\hat{J}=J(\hat{\tau},\hat{c}_{2};\hat{\bar{X}})$

7: return $(\hat{\tau},\hat{c}_{2},\hat{\bar{X}},\hat{R}_{e},\hat{c},\hat{J})$

######  Remark 6.1.

In the two Stackelberg equilibria, the numerical algorithms only differ in the
solution of producers’ problem.

###  6.2 Numerical Experiments

####  6.2.1 Contract Acceptance-Rejection by the Producers

First, we analyze how the producers use their option to walk-away under
different tax and penalty for not matching the demand policies. By symmetry,
the producers are statistically identical so they accept or reject the
contract all together. As we can see from Figure 4, the producers reject the
contract when the carbon tax $\tau$, or the penalty $c_{2}$ for not matching
the demand is high, both in the MFC and MFG cases. Notice that there are cases
for which the MFG fixed point algorithm fails to converge when $\tau$ and
$c_{2}$ levels are in the middle values. This is because when the fixed
average production is high, the revenue is low with respect to the costs, and
the representative producer decides not to produce, and the new average
production decreases. In the next step of the iteration, the fixed average
production is low and since the revenue is high, the representative producer
decides to produce.

Figure 4: Contract Acceptance-Rejection Regions in MFC (left) and in MFG
(right)

####  6.2.2 Analysis of Regulator’s Cost

For the experiments of this section, we used the same parameters as for the
producers’ model in the previous section. For the regulator we used:

$\alpha_{1}=1$, $\alpha_{2}=3.3$, $\alpha_{3}=500$, $\alpha_{4}=0.01$,
$\alpha_{5}=0.25$
---
$\tau\in\\{0,10,15,20,25,30,40,50,75,100\\}$
$c_{2}\in\\{50,100,250,500,750,1000,1500,2000,2500,3000,4000,5000\\}$

Figure 5: Regulator Cost where $c_{2}$ (left) or $\tau$ (middle) kept fixed in
both MFC and MFG settings; the effect of the importance given to the demand
matching by the regulator, $\alpha_{4}$, on the decision of the penalty
imposed on minor player for not matching the demand (right)

First, we analyze the regulator’s expected cost for different values of the
carbon tax given a fixed penalty for not matching the demand. Then we switch
the roles of the two controls of the regulator. The left and middle subfigures
in Figure 5 show that the cost of the regulator is convex as a function of the
carbon tax or the penalty, when the other control is fixed. It also shows the
importance given to demand matching by the regulator by tracking the effect of
$\alpha_{4}$ in regulator’s cost. The right subfigure in Figure 5 shows that
when the tax is fixed, the regulator’s minimum cost is attained at higher
$c_{2}$ values when $\alpha_{4}$ is higher. This shows that the regulator
should impose higher penalties for not matching the demand to producers when
the demand matching is more important. This subfigure is for the MFC case but
similar results hold in the MFG case as well.

Finally, Figure 6 gives 3D plots of the regulator cost as a function of their
controls $\tau$ and $c_{2}$. Here, the minimum is attained at
$(\tau^{*},c_{2}^{*})=(50,1000)$ in the MFC case and at
$(\hat{\tau},\hat{c}_{2})=(75,1000)$ in the MFG case. Also, we see that for
any given tax and penalty, the expected cost of the regulator is higher if the
producers are cooperative instead of competitive when the regulator gives more
importance to demand matching. This is because for any given couple
$(\tau,c_{2})$, in the cooperative setting producers are behaving like a big
monopolistic firm and care less about matching the demand than in the
competitive setting in order to maximize their revenues by keeping the prices
higher. When the demand matching is important for the regulator, the regulator
benefits from the competition among the electricity producers even if this
competition creates adverse effect for the producers themselves.

Figure 6: Regulator cost given admissible $c_{2}$ and $\tau$ values in when
the regulator cares about matching the demand in MFC (left) and in MFG
(middle); the difference of the regulator’s cost between MFC and MFG given any
admissible $c_{2}$ and $\tau$ couples (right)

##  7 Conclusion

In this paper, we investigate the behavior of rational electricity producers
in the presence of a carbon tax. We analyze how they manage the trade-off
between reliance on traditional and predictable fossil fuel power production
assets which emit greenhouse gas and hence cost revenues because of the carbon
tax, and the temptation to invest in clean energy production assets which will
not be the source of emissions but which make matching the demand problematic
because of the volatility of their output. We study a large population of
producers in two different models, a first one in which they compete and
hopefully reach a Nash equilibrium, and a second one in which they cooperate
and rely on the solution of a centralized optimization problem. In a second
set of models, we introduce a regulator choosing the level of the carbon tax
in hope to control the overall emissions in the economy, and a penalty to be
imposed on producers failing to meet the demand in hope to avoid power outages
and reputation costs. Our models are based on recent progress in the theory
and the numerical analysis of mean field games and mean field control
problems.

We showed that when the producers cooperate, they are better off by behaving
like a single monopolistic firm. However, if the regulator raises excessively
the penalty to match the demand, they can take advantage of the competitive
behavior of the producers. While our models remain stylized, they open the
door to more complex models, e.g. involving time dependent policies the
regulator could base on the response of the producers. Furthermore, our models
could be the used to include more features of the energy markets such as
storage, and the interactions between neighboring states or countries.

##  Appendix A Proofs of Theorems

###### Proof.

Proof of Lemma 3.1 Following [12, ch. 3], we write the Hamiltonian:

$H(t,N,X,\bar{X},q)=(Ax)^{\top}q+(B\cdot
N)^{\top}q+C_{t}^{\top}q+\frac{R}{2}|N|^{2}+H_{t}^{\top}X+\bar{X}^{\top}FX+X^{\top}GX+J_{t}$
(A.1)

where $q$ is the adjoint process. Therefore, the optimal $N$ to optimize $H$
can be expressed as:

$\hat{N}(q)=-R^{-1}B^{\top}q.$ (A.2)

By plugging $\hat{N}(q)$ in the Hamiltonian (A.1), the optimal Hamiltonian can
be written as:

$\hat{H}(t,X,\bar{X},q)=-\frac{1}{2}q^{\top}BR^{-1}B^{\top}q+X^{\top}A^{\top}q+C_{t}^{\top}q+H_{t}^{\top}X+\bar{X}^{\top}FX+X^{\top}GX+J_{t}.$
(A.3)

Then the Hamilton Jacobi Bellman (HJB) equation can be written as:

$-\frac{\partial u(t,X)}{\partial
t}-tr(aD^{2}u(t,X))=\hat{H}(t,X,\bar{X}_{t},Du(t,X)),\qquad\qquad
u(T,X)=X^{\top}S_{T}X+p_{2}R_{e},$ (A.4)

where $\bar{X}_{t}=\int_{\mathbb{R}^{5}}Xm(t,X)dX$. The Kolmogorov Fokker
Planck (KFP) equation for our MFG problem can be written as:

$\displaystyle\dfrac{\partial m(t,X)}{\partial
t}-tr(aD^{2}m(t,X))+\nabla_{x}\big{(}m(t,X)(AX-
Br^{-1}B^{\top}Du(t,X)+C_{t})\big{)}$ $\displaystyle=0,\quad m(0,X)=m_{0}(X).$
(A.5)

Introduce the following ansatz for the value function:

$u(t,X)=\frac{1}{2}X^{\top}{\eta_{t}}X+X^{\top}{r_{t}}+s_{t}.$ (A.6)

We have $Du(t,X)=\eta_{t}X+r_{t}$ and $D^{2}u(t,x)=\eta_{t}$. By plugging
(A.6) into the HJB equation in (A.4), we obtain that $\eta_{t}$ is the
solution of the following symmetric matrix Riccati equation:

$\frac{d{\eta_{t}}}{dt}-{\eta_{t}}BR^{-1}B^{\top}{\eta_{t}}+A^{\top}{\eta_{t}}+{\eta_{t}}A+2G=0,\qquad\qquad{\eta_{T}}=2S_{T}.$
(A.7)

This Riccati equation has a unique positive symmetric solution, see [22, ch.
14.3]. By plugging (A.6) into the HJB equation (A.4), we also obtain the
differential equation for $r_{t}$ that is coupled with $\bar{X}_{t}$:

$-\frac{d{r_{t}}}{dt}=\left(A^{\top}-{\eta_{t}}BR^{-1}B^{\top}\right){r_{t}}+{\eta_{t}}C_{t}+H_{t}+F^{\top}{\bar{X}_{t}},\qquad\qquad{r_{T}}=0.$
(A.8)

The differential equation for $\bar{X}_{t}$ can be found by plugging the
ansatz in the KFP equation:

$\displaystyle\frac{d\bar{X}_{t}}{dt}$
$\displaystyle=\frac{d}{dt}\int_{\mathbb{R}^{5}}X\cdot
m(t,X)dX=\int_{\mathbb{R}^{5}}X\cdot\frac{\partial m(t,X)}{\partial t}dX$
(A.9)
$\displaystyle=\int_{\mathbb{R}^{5}}X\cdot\Big{(}tr(aD^{2}m(t,X))-\nabla_{x}\big{(}m(t,X)(AX-
BR^{-1}B^{\top}(\eta_{t}X+r_{t})+C_{t})\big{)}\Big{)}dX$
$\displaystyle=-\int_{\mathbb{R}^{5}}X\cdot\frac{\partial m(t,X)}{\partial
X}(AX-BR^{-1}B^{\top}(\eta_{t}X+r_{t})+C_{t})dX-\int_{\mathbb{R}^{5}}X\cdot
m(t,X)(A-BR^{-1}B^{\top}\eta_{t})dX$
$\displaystyle=(A-BR^{-1}B^{\top}\eta_{t})\bar{X}_{t}-BR^{-1}B^{\top}r_{t}+C_{t}.$

Finally from the HJB equation where the ansatz is plugged in we find that:

$\frac{ds_{t}}{dt}=-tr(a\eta_{t})+\frac{1}{2}{r_{t}}^{\top}BR^{-1}B^{\top}{r_{t}}-C_{t}^{\top}{r_{t}}-J_{t},\qquad\qquad
s_{T}=p_{2}R_{e}.$

Therefore we have:

$s_{t}=p_{2}R_{e}+\int_{t}^{T}\Big{(}tr(a{\eta_{s}})-\frac{1}{2}{r_{s}}^{T}BR^{-1}B^{\top}{r_{s}}+C_{s}^{\top}{r_{s}}+J_{s}\Big{)}ds.$
(A.10)

The expected cost of the representative minor player given fixed mean field
and $R_{e}$ can be calculated by using:

$\displaystyle\begin{split}\inf_{(N_{t})_{t}}\tilde{C}^{MFG}\Big{(}N;R_{e},\bar{X}\Big{)}&=\mathbb{E}\left[u(0,X_{0})\right]\\\
&=\mathbb{E}\left[\frac{1}{2}X_{0}^{\top}\eta_{0}X_{0}+X_{0}^{\top}r_{0}+s_{0}\right]\\\
&=\frac{1}{2}\left(Var(\sqrt{\eta_{0}}X_{0})+\mathbb{E}[\sqrt{\eta_{0}}X_{0}]^{2}\right)+\bar{X}_{0}^{\top}r_{0}+s_{0}.\end{split}$

∎

###### Proof.

Proof of Theorem 3.2 For the existence and uniqueness proof, we make use of
Banach Fixed Point theorem. We follow the line of proof used for a stochastic
system in [25, Thm 5.1]. First we fix $(r_{t}^{1})_{t}$ and $(r_{t}^{2})_{t}$,
then corresponding $(\bar{X}^{i}_{t})_{t}$ can be found by solving the
following ODE:

$d\bar{X}_{t}^{i}=[(A-BR^{-1}B^{\top}\eta_{t})\bar{X}_{t}^{i}-BR^{-1}B^{\top}r_{t}^{i}+C_{t}]dt,\qquad\bar{X}_{0}^{i}=\bar{x}_{0},\qquad
i=\\{1,2\\}.$

Further let $\widetilde{\bar{X}}_{t}=\bar{X}_{t}^{1}-\bar{X}_{t}^{2}$ and
$\tilde{r}_{t}=r_{t}^{1}-r_{t}^{2}$, then we have:

$d\widetilde{\bar{X}}_{t}=[(A-BR^{-1}B^{\top}\eta_{t})\widetilde{\bar{X}}_{t}-BR^{-1}B^{\top}\tilde{r}_{t}]dt,\qquad\widetilde{\bar{X}}_{0}=0.$

Now we introduce $(r_{t}^{i^{\prime}})_{t}$ that solves:

$dr^{i^{\prime}}_{t}=[(\eta_{t}BR^{-1}B^{\top}-A^{\top})r^{i}_{t}-\eta_{t}C_{t}-H_{t}-F^{\top}\bar{X}_{t}^{i}]dt,\qquad
r^{i^{\prime}}_{T}=0,\qquad\forall i\in\\{1,2\\}.$

and let $\tilde{r}_{t}^{{}^{\prime}}=r_{t}^{1^{\prime}}-r_{t}^{2^{\prime}}$.
Then we have the following ODE:

$d\tilde{r}^{{}^{\prime}}_{t}=[(\eta_{t}BR^{-1}B^{\top}-A^{\top})\tilde{r}_{t}-F^{\top}\widetilde{\bar{X}}_{t}]dt,\qquad\tilde{r}^{{}^{\prime}}_{T}=0.$

Therefore we have defined a mapping $r\mapsto r^{{}^{\prime}}$. We now show
that it is a contraction mapping to be able to use the Banach Fixed Point
theorem.

Step 1.  Using Itô’s formula, we write the dynamics for
$||\widetilde{\bar{X}}_{t}||^{2}$:

$d||\widetilde{\bar{X}}_{t}||^{2}=2(\widetilde{\bar{X}}_{t})^{\top}d\widetilde{\bar{X}}_{t}=2(\widetilde{\bar{X}}_{t})^{\top}[(A-BR^{-1}B^{\top}\eta_{t})\widetilde{\bar{X}}_{t}-BR^{-1}B^{\top}\tilde{r}_{t}]dt.$

By using these dynamics we can find a bound for
$||\widetilde{\bar{X}}_{t}||^{2}$:

$\displaystyle||\widetilde{\bar{X}}_{t}||^{2}$
$\displaystyle=\int_{0}^{t}2(\widetilde{\bar{X}}_{s})^{\top}[(A-BR^{-1}B^{\top}\eta_{s})\widetilde{\bar{X}}_{s}-BR^{-1}B^{\top}\tilde{r}_{s}]ds$
$\displaystyle\leq\int_{0}^{t}(||2(A-BR^{-1}B^{\top}\eta_{s})||)||\widetilde{\bar{X}}_{s}||^{2}ds+\int_{0}^{t}||BR^{-1}B^{\top}||2<\widetilde{\bar{X}}_{s},\tilde{r}_{s}>ds$
$\displaystyle\leq\int_{0}^{t}(||2(A-BR^{-1}B^{\top}\eta_{s})||)||\widetilde{\bar{X}}_{s}||^{2}ds+\int_{0}^{t}||BR^{-1}B^{\top}||\big{(}||\widetilde{\bar{X}}_{s}||^{2}+||\tilde{r}_{s}||^{2}\big{)}ds$
$\displaystyle\leq\exp\Big{(}\int_{0}^{t}(2||A-BR^{-1}B^{\top}\eta_{s}||+||BR^{-1}B^{\top}||)ds\Big{)}\int_{0}^{t}(||BR^{-1}B^{\top}||)||\tilde{r}_{s}||^{2}ds$
$\displaystyle\leq C^{(1)}\int_{0}^{T}||\tilde{r}_{s}||^{2}ds$ (A.11)

where the third to last inequalities stem from the Gronwall’s inequality, and
we define
$C^{(1)}=\exp\Big{(}T\big{(}2||A||+2(||BR^{-1}B^{\top}||)||\eta||_{T}+||BR^{-1}B^{\top}||\big{)}\Big{)}||BR^{-1}B^{\top}||$
with $||\eta||_{T}:=\sup_{0\leq t\leq T}||\eta_{t}||$.

Step 2.  We write the dynamics for $||\tilde{r}_{t}^{{}^{\prime}}||^{2}$:

$d||\tilde{r}^{{}^{\prime}}_{t}||^{2}=2(\tilde{r}^{{}^{\prime}}_{t})^{\top}d\tilde{r}^{{}^{\prime}}_{t}=2(\tilde{r}^{{}^{\prime}}_{t})^{\top}[(\eta_{t}BR^{-1}B^{\top}-A^{\top})\tilde{r}_{t}-F^{\top}\widetilde{\bar{X}}_{t}]dt.$

Now we find a bound for $||\tilde{r}_{t}^{{}^{\prime}}||^{2}$ as follows by
using Young’s inequality:

$\displaystyle||\tilde{r}^{{}^{\prime}}_{t}||^{2}$
$\displaystyle=\int_{t}^{T}2(\tilde{r}^{{}^{\prime}}_{s})^{\top}[(A^{\top}-\eta_{s}BR^{-1}B^{\top})\tilde{r}_{s}+F^{\top}\widetilde{\bar{X}}_{s}]ds$
$\displaystyle\leq\int_{t}^{T}||A^{\top}-\eta_{s}BR^{-1}B^{\top}||2<\tilde{r}^{{}^{\prime}}_{s},\tilde{r}_{s}>ds+\int_{t}^{T}||F^{\top}||2<\tilde{r}^{{}^{\prime}}_{s},\widetilde{\bar{X}}_{s}>ds$
$\displaystyle=\int_{t}^{T}(||A^{\top}-\eta_{s}BR^{-1}B^{\top}||+||F^{\top}||)||\tilde{r}^{{}^{\prime}}_{s}||^{2}ds+\int_{t}^{T}(||A^{\top}-\eta_{s}BR^{-1}B^{\top}||)||\tilde{r}_{s}||^{2}ds+\int_{t}^{T}(||F^{\top}||)||\widetilde{\bar{X}}_{s}||^{2}ds.$

Now the expression found in (A) can be plugged in and by using Gronwall’s
inequality:

$\displaystyle||\tilde{r}^{{}^{\prime}}_{t}||^{2}$
$\displaystyle\leq\int_{t}^{T}(||A^{\top}-\eta_{s}BR^{-1}B^{\top}||+|F^{\top}||)||\tilde{r}^{{}^{\prime}}_{s}||^{2}ds+\int_{0}^{T}(||A^{\top}-\eta_{s}BR^{-1}B^{\top}||)||\tilde{r}_{s}||^{2}ds+$
$\displaystyle$\displaystyle \int_{0}^{T} ||F^{\top}|| C^{(1)}
\Big{(}\int_{0}^{T} ||\tilde{r}_{s}||^{2} ds\Big{)} ds$\ignorespaces$
$\displaystyle\leq\int_{t}^{T}(||A^{\top}-\eta_{s}BR^{-1}B^{\top}||+|F^{\top}||)||\tilde{r}^{{}^{\prime}}_{s}||^{2}ds+\int_{0}^{T}(||A^{\top}-\eta_{s}BR^{-1}B^{\top}||)||\tilde{r}_{s}||^{2}ds+$
$\displaystyle$\displaystyle T\Big{[} ||F^{\top}|| C^{(1)} \Big{(}\int_{0}^{T}
||\tilde{r}_{s}||^{2} ds\Big{)} \Big{]}$\ignorespaces$
$\displaystyle\leq\exp\Big{(}{T(||A||+(||BR^{-1}B^{\top}||)||\eta||_{T}+||F^{\top}||)}\Big{)}$
$\displaystyle$\displaystyle\int_{0}^{T}
\Big{[}||A^{\top}-\eta_{s}BR^{-1}B^{\top} || + T||F^{\top}|| C^{(1)}
\Big{]}||\tilde{r}_{s}||^{2} ds$\ignorespaces.$

Now define $||r||_{T}:=\sup_{0\leq t\leq T}||r_{t}||$, then we have:
$||\tilde{r}^{{}^{\prime}}||_{T}^{2}\leq c_{T}||\tilde{r}||_{T}^{2}$, where
$c_{T}$ is

$c_{T}=Te^{{T(||A||+(||BR^{-1}B^{\top}||)||\eta||_{T}+||F^{\top}||)}}\big{(}||A^{\top}||+(||\eta||_{T}+T||F^{\top}||e^{T(2||A||+(2||\eta||_{T}+1)||BR^{-1}B^{\top}||)})||BR^{-1}B^{\top}||\big{)}.$

With small $T$ we have $c_{T}<1$, which concludes the proof. ∎

###### Proof.

Proof of Lemma 3.3 For Mean Field Control problems we have the following HJB
and FP systems. The detailed proof of the derivation can be found in [8, ch.
6].

$\displaystyle-\frac{\partial u(t,X)}{\partial
t}-tr(aD^{2}u(t,X))=\hat{H}(t,X,m(t),Du(t,X))+$
$\displaystyle\int_{\mathbb{R}^{n}}\frac{\partial\hat{H}}{\partial
m}(t,\xi,m(t),Du(t,\xi))(x)m(t,\xi)d\xi,$ $\displaystyle$\displaystyle u(T,X)
= X^{\top} S_{T} X + p_{2} R_{e}$\ignorespaces$ $\displaystyle\dfrac{\partial
m(t,X)}{\partial t}-tr(aD^{2}m(t,X))+\nabla_{x}\big{(}m(t,X)(AX-BR^{-1}$
$\displaystyle B^{\top}Du(t,X)+C_{t})\big{)}=0,\qquad m(0,X)=x_{0}$

where $\frac{\partial\hat{H}}{\partial m}$ denotes the Gâteaux differential of
$\hat{H}$ on $L^{2}(\mathbb{R}^{5})$. As it can be seen, the KFP equation
stays same as in MFC but HJB equation changes. By rewriting (A.3) as:

$\hat{H}(t,X,m,q)=-\frac{1}{2}q^{\top}BR^{-1}B^{\top}q+X^{\top}A^{\top}q+C_{t}^{\top}q+H_{t}^{\top}X+\Big{(}\int_{\mathbb{R}^{n}}\xi
m(\xi)d\xi\Big{)}^{\top}FX+X^{\top}GX+J_{t}.$ (A.12)

We find that

$\displaystyle\int_{\mathbb{R}^{n}}\frac{\partial\hat{H}(t,X,m,q)}{\partial
m}(x)m(\xi)d\xi=\int_{\mathbb{R}^{n}}X^{\top}F\xi
m(\xi)d\xi=X^{\top}F\bar{X}.$ (A.13)

Therefore, the HJB equation becomes:

$-\frac{\partial u}{\partial
t}-tr(aD^{2}u)=\hat{H}(t,X,\bar{X},Du)+X^{\top}F\bar{X}\qquad\qquad
u(X_{T},T)=X_{T}^{\top}S_{T}X_{T}+p_{2}R_{e}.$ (A.14)

We introduce the same ansatz as in (A.6). By plugging this ansatz in the HJB
equation given in (A.14), we end up with the same Riccati equation and
equation for $s_{0}$. Only the differential equation of $r_{t}$ changes as
follows:

$-\frac{d{r_{t}}}{dt}=\left(A^{\top}-{\eta_{t}}BR^{-1}B^{\top}\right){r_{t}}+{\eta_{t}}C_{t}+H_{t}+F^{\top}{\bar{X}_{t}}+{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color@cmyk@fill{0}{0.75}{1}{0.24}F{\bar{X}_{t}}},\qquad\qquad{r_{T}}=0,$

where $\bar{X}_{t}=\int_{\mathbb{R}^{5}}Xm(t,X)dX$. Since the Kolmogorov-
Fokker-Planck equation stayed the same, we have the same expression for the
differential equation of $\bar{X}_{t}$ as in (A.9). We obtain the MFC cost
given any fixed $R_{e}$ by using the ansatz (see e.g. [24] for more details):

$\displaystyle\begin{split}\inf_{(N_{t})_{t}}\tilde{C}^{MFC}\Big{(}N;R_{e}\Big{)}&=\mathbb{E}\left[u(0,X_{0}){\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color<EMAIL_ADDRESS>&=\mathbb{E}\left[\frac{1}{2}X_{0}^{\top}\eta_{0}X+X_{0}^{\top}r_{0}+s_{0}\right]{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color<EMAIL_ADDRESS>&=\frac{1}{2}\left(Var(\sqrt{\eta_{0}}X_{0})+\mathbb{E}[\sqrt{\eta_{0}}X_{0}]^{2}\right)+\bar{X}_{0}^{\top}r_{0}+s_{0}{\color[rgb]{0.76,0.01,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.76,0.01,0}\pgfsys@color@cmyk@stroke{0}{0.75}{1}{0.24}\pgfsys@color<EMAIL_ADDRESS>
∎

###### Proof.

Proof of Theorem 3.4 For the sake of space, we omit the proof of existence and
uniqueness, which follows the same steps as in the proof of theorem 3.2. ∎

###### Proof.

Proof of Remark 3.5 Assume that the strategy couple $(\hat{N},\hat{R}_{e})$ is
optimal in the mean field game and the corresponding mean field is given as
$\bar{Q}=\bar{Q}(\hat{N},\hat{R}_{e})$. Now, assume that the representative
player deviates from the optimal strategy and uses
$(\widehat{N}+\epsilon\widecheck{N},\widehat{R}_{e}+\epsilon\widecheck{R}_{e})$.
Then:

$\displaystyle\dfrac{dC(\widehat{N}+\epsilon\widecheck{N},\widehat{R}_{e}+\epsilon\widecheck{R}_{e};\bar{Q})}{d\epsilon}|_{\epsilon=0}=\mathbb{E}\Big{[}\int_{0}^{T}\Big{[}$
$\displaystyle
2{c_{1}\widehat{N}_{t}\widecheck{N}_{t}}+{p_{1}\widecheck{\tilde{N}}_{t}}+{2c_{2}(\widehat{Q}_{t}-D_{t})\widecheck{Q}_{t}}$
(A.15)
$\displaystyle-{c_{3}\big{(}\rho_{0}+\rho_{1}(D_{t}-\bar{Q}_{t})\big{)}\widecheck{Q}_{t}}\Big{]}dt+{2\tau\widehat{P}_{T}\widecheck{P}_{T}}+{p^{\prime}(\widehat{R}_{e})}\widecheck{R}_{e}\Big{]}$

Furthermore we have the following dynamics:

$\displaystyle d\widecheck{Q}_{t}$
$\displaystyle=\kappa_{1}\widecheck{N}_{t}dt+\kappa_{2}\widecheck{R}_{e}\left(\alpha\cos(\alpha
t)dt+(\theta-S_{t})dt+\sigma_{0}d\widecheck{W}_{t}\right)\qquad$
$\displaystyle d\widecheck{E}_{t}=\delta\widecheck{N}_{t}dt$ $\displaystyle
d\widecheck{P}_{t}$ $\displaystyle=\widecheck{E}_{t}dt$ $\displaystyle
d\widecheck{\tilde{N}}_{t}=\widecheck{N}_{t}dt,$

with initial conditions:
$\widecheck{Q}_{0}=\widecheck{E}_{0}=\widecheck{P}_{0}=\widecheck{\tilde{N}}_{0}=0$.
We can introduce the adjoint variables with the following dynamics:

$\displaystyle dY^{1}_{t}$
$\displaystyle=\Big{(}-{2c_{2}(\widehat{Q}_{t}-D_{t})}+{c_{3}\big{(}\rho_{0}+\rho_{1}(D_{t}-\bar{Q}_{t})\big{)}}\Big{)}dt+Z_{t}^{1,1}d\widetilde{W}_{t}+Z_{t}^{1,2}dW_{t}\qquad\qquad$
$\displaystyle Y_{T}^{1}=0$ $\displaystyle dY^{2}_{t}$
$\displaystyle=\Big{(}\kappa_{2}\widehat{R}_{e}Y^{1}_{t}+Y^{2}_{t}\Big{)}dt+Z_{t}^{2,1}d\widetilde{W}_{t}+Z_{t}^{2,2}dW_{t}$
$\displaystyle Y_{T}^{2}=0$ $\displaystyle dY^{3}_{t}$
$\displaystyle=-Y_{t}^{4}dt+Z_{t}^{3,1}d\widetilde{W}_{t}+Z_{t}^{3,2}dW_{t}$
$\displaystyle Y_{T}^{3}=0$ $\displaystyle dY^{4}_{t}$
$\displaystyle=Z_{t}^{4,1}d\widetilde{W}_{t}+Z_{t}^{4,2}dW_{t}$ $\displaystyle
Y_{T}^{4}=2\tau\widehat{P}_{T}$ $\displaystyle dY^{5}_{t}$
$\displaystyle=-p_{1}dt+Z_{t}^{5,1}d\widetilde{W}_{t}+Z_{t}^{5,2}dW_{t}$
$\displaystyle Y_{T}^{5}=0.$

Now plug these dynamics in the perturbed cost function (A.15) and apply
integration by parts:

$\displaystyle\dfrac{dC(\widehat{N}+\epsilon\widecheck{N},\widehat{R}_{e}+\epsilon\widecheck{R}_{e};\bar{Q})}{d\epsilon}|_{\epsilon=0}$
$\displaystyle=\mathbb{E}\Big{[}\int_{0}^{T}\Big{[}2{c_{1}\widehat{N}_{t}\widecheck{N}_{t}}+{\big{(}-dY_{t}^{5}+Z_{t}^{5,1}d\widetilde{W}_{t}+Z_{t}^{5,2}dW_{t}\big{)}\widecheck{\tilde{N}}_{t}}+\big{(}-dY_{t}^{1}+Z_{t}^{1,1}d\widetilde{W}_{t}+Z_{t}^{1,2}dW_{t}\big{)}\widecheck{Q}_{t}\Big{]}dt$
$\displaystyle$\displaystyle+ {Y^{4}_{T}\widecheck{P}_{T}} +
{p^{\prime}(\widehat{R}_{e})}\widecheck{R}_{e}\Big{]}$\ignorespaces$
$\displaystyle=\mathbb{E}\Big{[}\int_{0}^{T}\widecheck{N}_{t}\Big{(}2{c_{1}\widehat{N}_{t}}+Y_{t}^{5}+Y_{t}^{1}\kappa_{1}+Y_{t}^{3}\delta\Big{)}dt\Big{]}+\mathbb{E}\Big{[}\int_{0}^{T}\left(\kappa_{2}Y_{t}^{1}\Big{(}\alpha\cos(\alpha
t)+(\theta{-S_{t}})\Big{)}dt+{p^{\prime}(\widehat{R}_{e})}\right)\widecheck{R}_{e}\Big{]}.$

By optimality, the above expression should be equal to 0 for any given
$\widecheck{N}_{t}$ and $\widecheck{R}_{e}$; therefore:

$\widehat{N}_{t}=-\frac{Y_{t}^{1}\kappa_{1}+Y_{t}^{3}\delta+Y_{t}^{5}}{2c_{1}}\quad\text{and}\quad\hat{R}_{e}=(p^{\prime})^{-1}\Big{(}-\mathbb{E}\Big{[}\int_{0}^{T}\kappa_{2}Y_{t}^{1}\Big{(}\alpha\cos(\alpha
t)+(\theta{-S_{t}})\Big{)}dt\Big{]}\Big{)}.$ (A.16)

For the proof of the FBSDE system in the mean field control setting, assume
that strategy couple $(\widehat{N},\widehat{R}_{e})$ is optimal in the mean
field control problem and the corresponding mean field is given as
$\widehat{\bar{Q}}$. Now assume that the representative player deviates from
optimal strategy and uses
$(\widehat{N}+\epsilon\widecheck{N},\widehat{R}_{e}+\epsilon\widecheck{R}_{e})$.
Since in the mean field control version every player is going to deviate, the
mean field term also changes to
$\widehat{\bar{Q}}+\epsilon\widecheck{\bar{Q}}$. Then:

$\displaystyle\dfrac{dC({\widehat{N}}+\epsilon\widecheck{N},\widehat{R}_{e}+\epsilon\widecheck{R}_{e};\bar{Q})}{d\epsilon}|_{\epsilon=0}$
(A.17)
$\displaystyle\qquad=\mathbb{E}\Big{[}\int_{0}^{T}\Big{[}2{c_{1}\widehat{N}_{t}\widecheck{N}_{t}}+{p_{1}\widecheck{\tilde{N}}_{t}}+{2c_{2}(\widehat{Q}_{t}-D_{t})\widecheck{Q}_{t}}-{c_{3}\big{(}\rho_{0}+\rho_{1}(D_{t}-2{\widehat{\bar{Q}}}_{t})\big{)}\widecheck{Q}_{t}}\Big{]}dt+{2\tau\widehat{P}_{T}\widecheck{P}_{T}}+{p^{\prime}(\widehat{R}_{e})}\widecheck{R}_{e}\Big{]}.$

We introduce the same adjoint variables as for the MFG FBSDE except that for
$Y_{t}^{1}$ we take:

$dY^{1}_{t}=\Big{(}-{2c_{2}(\widehat{Q}_{t}-D_{t})}+{c_{3}\big{(}\rho_{0}+\rho_{1}(D_{t}-2\widehat{\bar{Q}}_{t})\big{)}}\Big{)}dt+Z_{t}^{1,1}d\widecheck{W}_{t}+Z_{t}^{1,2}dW_{t}\qquad\quad
Y_{T}^{1}=0.$

By plugging the adjoint variable dynamics in the perturbed cost
(LABEL:eq:perturbedcostmfc) and applying integration by parts we end up with
the same optimality conditions as in (A.16). ∎

## References

  * [1] Y. Achdou, C. Bertucci, J.-M. Lasry, P.-L. Lions, A. Rostand, and J. Scheinkman, A class of short-term models for the oil industry addressing speculative storage, arXiv preprint arXiv:2003.11790, (2020).
  * [2] R. Aïd, M. Basei, and H. Pham, A mckean–vlasov approach to distributed electricity generation development, Mathematical Methods of Operations Research, 91 (2020), pp. 269–310.
  * [3] R. Aid, R. Dumitrescu, and P. Tankov, The entry and exit game in the electricity markets: a mean-field game approach, arXiv: Optimization and Control, (2020).
  * [4] 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 (2020), pp. 644–670.
  * [5] A. Aurell, R. Carmona, G. Dayanikli, and M. Lauriere, Optimal incentives to mitigate epidemics: a Stackelberg mean field game approach, 2020\.
  * [6] O. Bahn, A. Haurie, and R. Malhame, Limit Game Models for Climate Change Negotiations, 07 2017, pp. 27–47.
  * [7] A. Bensoussan, M. H. M. Chau, and S. C. P. Yam, Mean field games with a dominating player, Applied Mathematics & Optimization, 74 (2016), pp. 91–128.
  * [8] A. Bensoussan, J. Frehse, and P. Yam, Mean field games and mean field type control theory, 2013.
  * [9] P. E. Caines, M. Huang, and R. P. Malhamé, Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the nash certainty equivalence principle, Communications in Information and Systems, 6 (2006), p. 221–252.
  * [10] R. Carmona and F. Delarue, Probabilistic analysis of mean field games, SIAM Journal on Control and Optimization, 51 (2013), pp. 2705 – 2734\.
  * [11] , Forward-backward stochastic differential equations and controlled McKean-Vlasov dynamics, Annals of Probability, 43 (2015), pp. 2647–2700.
  * [12] R. Carmona and F. Delarue, FBSDEs and the Solution of MFGs without common noise, Springer International Publishing, Cham, 2018, pp. 215–345.
  * [13] R. Carmona and M. Fehr, The clean development mechanism and CER price formation in the carbon emissions markets, (2010), pp. 341 – 384.
  * [14] R. Carmona, M. Fehr, and J. Hinz, Optimal stochastic control and carbon price formation, SIAM Journal on Control and Optimization, 48 (2009), pp. 2168–2190.
  * [15] R. Carmona, M. Fehr, J. Hinz, and A. Porchet, Market design for emissions markets trading schemes, SIAM Review, 52 (2010), pp. 403–452.
  * [16] R. Carmona and P. Wang, Finite-state contract theory with a principal and a field of agents, Management Science, (2020), p. First on line.
  * [17] P. Chan and R. Sircar, Fracking, renewables, and mean field games, SIAM Review, 59 (2017), pp. 588–615.
  * [18] B. Djehiche, J. Barreiro-Gomez, and H. Tembine, Electricity price dynamics in the smart grid : A mean-field-type game perspective, 2018.
  * [19] R. Élie, E. Hubert, T. Mastrolia, and D. Possamaï, Mean–field moral hazard for optimal energy demand response management, Mathematical Finance, 31 (2021), pp. 399–473.
  * [20] R. Élie, T. Mastrolia, and D. Possamaï, A tale of a principal and many, many agents, Mathematics of Operations Research, 44 (2019), pp. 440–467.
  * [21] O. Guéant, J.-M. Lasry, and P. L. Lions, Mean Field Games and oil production, in The Economics of Sustainable Development, Economica, ed., 2010.
  * [22] V. Kučera, Riccati equations and its solutions, in Control System Advanced Methods, W. S. Levine, ed., CRC Press, Inc, 2nd ed., 2009, ch. 14.
  * [23] J.-M. Lasry and P.-L. Lions, Mean field games, Japanese Journal of Mathematics, 2 (2007), p. 229–260.
  * [24] M. Laurière, On numerical methods for mean field games and mean field type control.  Proc. AMS Short Course, 2020.
  * [25] J. Ma and J. Yong, Forward-Backward Stochastic Differential Equations and their Applications, Springer Berlin Heidelberg, Berlin, Heidelberg, 2007.
  * [26] R. Salhab, R. Malhamé, and J. L. Ny, A dynamic collective choice model with an advertiser, 2016 IEEE 55th Conference on Decision and Control (CDC), (2016), pp. 6098–6104.
  * [27] A. Shrivats, D. Firoozi, and S. Jaimungal, Optimal generation and trading in solar renewable energy certificate (srec) markets, Applied Mathematical Finance, 27 (2020), pp. 131–99.
  * [28] S. Zhang, X. Wang, and A. Shanain, Modeling and computation of mean field equilibria in producers’ game with emission permits trading, Communications in Nonlinear Science and Numerical Simulation, 37 (2016), pp. 238–248.