# Team-Optimal Solution of Finite Number of Mean-Field Coupled LQG Subsystems

Jalal Arabneydi and Aditya Mahajan This work was supported by the Natural
Sciences and Engineering Research Council of Canada through Grant NSERC-RGPIN
402753-11.Jalal Arabneydi and Aditya Mahajan are with the Department of
Electrical and Computer Engineering, McGill University, Montreal, QC, Canada.
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

A decentralized control system with linear dynamics, quadratic cost, and
Gaussian disturbances is considered. The system consists of a finite number of
subsystems whose dynamics and per-step cost function are coupled through their
mean-field (empirical average). The system has mean-field sharing information
structure, i.e., each controller observes the state of its local subsystem
(either perfectly or with noise) and the mean-field. It is shown that the
optimal control law is unique, linear, and identical across all subsystems.
Moreover, the optimal gains are computed by solving two decoupled Riccati
equations in the full observation model and by solving an additional filter
Riccati equation in the noisy observation model. These Riccati equations do
not depend on the number of subsystems. It is also shown that the optimal
decentralized performance is the same as the optimal centralized performance.
An example, motivated by smart grids, is presented to illustrate the result.

Proceedings of IEEE Conference on Decision and Control, 2015.

##  I Introduction

###  I-A Motivation and literature overview

Team theory investigates multi-agent decision problems in which all agents
share a common objective. Such problems arise in various applications
including cyber-physical systems, networked control systems, surveillance and
sensor networks, communication networks, smart grids, robotics, and
organizational economics.

Not much is known regarding the optimal control of general team problems. Most
of the results in the literature are for specific information structures such
as delayed state sharing [1], delaying observation sharing [2], periodic
sharing [3], belief sharing [4], mean-field sharing [5] and partial-history
sharing [6]. We refer the reader to [7] for an overview.

There are two main challenges in solving team problems:

  1. 1.

Team problems are _conceptually_ difficult due to the decentralized nature of
the information available to the controllers. In particular, controllers need
to cooperate and coordinate to minimize a common cost but they have different
information about the state of the environment. Due to this discrepancy of
information, dynamic programming, which is one of the main solution techniques
for the optimal design of stochastic systems, does not work directly for team
problems.

  2. 2.

Team problems are _computationally_ difficult. Even when a dynamic programming
decomposition is obtained, the solution complexity increases exponentially or
double exponentially with the number of controllers. So, it is difficult to
use the dynamic programming solution for systems with large number of
controllers.

Similar challenges exist in dynamic games; however, they have been
successfully resolved for a class of models known as _mean-field games_ [8, 9,
10, 11, 12, 13, 14]. The salient feature of these models is that the
agents/players/controllers are coupled in their dynamics and per-step cost
functions only through the mean-field (i.e., the empirical average or the
empirical distribution). If the number of players are large, then an
approximately optimal solution for these models can be obtained by solving the
infinite population limit. The infinite population limit is easier to solve
than the finite population model because, when the population is
asymptotically large, the action of a single controller does not affect the
mean-field. Therefore, the optimal solution can be obtained by solving two
coupled equations: a backward dynamic programming equation that determines the
best-response of a controller given the trajectory of the mean-field; and a
forward Fokker-Planck equation that determines the evolution of the mean-field
given the strategies of the individual controllers. A consistent solution of
these equations determine Nash equilibrium strategies. A desirable feature of
the solution is that the resultant control laws and the complexity of the
solution do not depend on the number of controllers.

Motivated by these results, we investigate _team-optimal_ (rather than Nash
equilibrium) solution of mean-field models. In this paper, we focus on linear-
quadratic models. In particular, the system consists of a _finite_ number of
subsystems; the dynamics of the subsystems are coupled through the mean-field
of the states. The per-step cost is also coupled through the mean-field. We
assume that the controllers have _mean-field sharing_ information structure
introduced in [5]; that is, each controller observes it’s local state (either
completely or with noise) and the mean-field of the states of all subsystems.
In some applications such as communication networks, the sharing of the mean-
field happens naturally. In other applications, such as robotic teams, it is
possible to share the mean-field using distributed consensus algorithms.

The main difference between our setup and mean-field games is the following:
(i) We assume a finite number of controllers and obtain team-optimal control
strategies. In mean-field games, one assumes asymptotically large number of
controllers and obtains Nash equilibrium strategies. (ii) We assume that each
controller observes the mean-field. In mean-field games, this assumption is
not made. However, when the number of controllers goes to infinity, both the
information structures are equivalent because the mean-field becomes a
deterministic process.

###  I-B Contributions and salient features of the results

The linear quadratic mean-field model is a decentralized system with non-
classical information structure. Our main result is to show that linear
control strategies are optimal. The model is neither partially nested nor
quadratic invariant. In contrast, almost all positive results for linear
quadratic systems with non-classical information structures are for models
that are either partially nested or quadratic invariant.

We also show that the optimal control laws at each subsystem are identical. In
general, identical control laws are _not_ optimal for systems with
exchangeable subsystems (see [5] for one an example). So, a natural question
is the following: Given a decentralized system with exchangeable subsystems,
under what conditions are optimal control laws identical across subsystems.
This question warrants further investigation.

Our solution, and the solution complexity, do not depend on the number of
subsystems. Irrespective of the number of subsystems, the optimal control law
is given by the solution of _two_ (backward) Riccati equations in the full
observation model, and an additional filter Riccati equation in the noisy
observation model. The parameters of these Riccati equations do not depend on
the number of subsystems. Thus, the optimal control laws can be computed
without any knowledge of the number of subsystems.

Since the optimal control laws do not depend on the number of subsystems, our
results remain valid in the infinite population limit. In the infinite
population limit, the mean-field becomes a deterministic process that can be
pre-computed at all controllers. Thus, mean-field sharing is equivalent to the
completely decentralized information structure (where each controller only
observes its local state) considered in [8, 9]. Thus, our results may be
viewed as an alternative derivation of the infinite population mean-field
games solution.

As an intermediate step of the proof, we show that the decentralized control
laws—and, hence, the decentralized performance—are same as centralized (i.e.,
when all control actions are chosen by a single controller that observes all
the information available to all the decentralized controllers). A natural
follow up question that warrants further investigation is the following: In
general decentralized control problems, when is decentralized performance the
same as the centralized performance? Or more generally, when is the
performance under one information structure the same as the performance under
a finer information structure?

From an implementation point of view, the above feature has an interesting
consequence. If we have the freedom to design the information structure, then
in the linear-quadratic mean-field systems there is no advantage of sharing
anything beyond the mean-field. Since, it is easy to share the mean-field
using distributed consensus algorithms, it is relatively easy to implement the
optimal centralized solution in a distributed manner.

Furthermore, as we already argued, in the infinite population limit, mean-
field sharing is equivalent to the completely decentralized setup (i.e., each
controller only observes its own local state). Since, there is no advantage of
sharing any information beyond the mean-field—which is already computable at
each controller—the optimal completely decentralized solution is also the
optimal centralized solution. An immediate corollary is that the control laws
obtained using the approach of mean-field games are also the optimal
centralized solution.

##  II Problem formulation and the main results

###  II-A Notation

Uppercase letters denote random variables/vectors and lowercase letters denote
their realization. $\mathds{P}(\cdot)$ denotes the probability of an event and
$\mathds{E}[\cdot]$ denotes the expectation of a random variable. $\mathds{R}$
denotes the set of real numbers.

For a sequence of (column) vectors $X$, $Y$, $Z$, …, the notation
$\operatorname{vec}(X,Y,Z,\dots)$ denotes the vector
$[X^{\intercal},Y^{\intercal},Z^{\intercal}]^{\intercal}$. The vector
$\operatorname{vec}(X_{1},\dots,X_{t})$ is also denoted by $X_{1:t}$.

Superscripts index subsystems and subscripts index time. Thus, $X^{i}_{t}$
denotes a variable at subsystem $i$ at time $t$. Bold letters denote
collection of variables at all subsystems. For example, $\mathbf{X}$ denotes
$(X^{1},\dots,X^{n})$. $\langle\mathbf{X}\rangle$ denotes the mean-field of
$\mathbf{X}$, i.e, $\langle\mathbf{X}\rangle=\frac{1}{n}\sum_{i=1}^{n}X^{i}$.

###  II-B Model and problem formulation

Consider a decentralized control system with $n$ homogeneous subsystems that
operate for a finite horizon $T$. A controller is co-located with each
subsystem. Let $X^{i}_{t}\in\mathbb{R}^{d_{x}}$ denote the state of subsystem
$i$ and $U^{i}_{t}\in\mathbb{R}^{d_{u}}$ denote the control action taken by
controller $i$, $i\in\\{1,\ldots,n\\}$, at time $t$. Let
$Z_{t}=\langle\mathbf{X}_{t}\rangle$ denote the mean-field of the states at
time $t$.

####  II-B1 System dynamics

The initial states $\\{X^{i}_{1}\\}_{i=1}^{n}$ are independent Gaussian random
variables with zero mean and covariance $\Sigma_{X}$. They evolve in time as
follows:

$X^{i}_{t+1}=A_{t}X^{i}_{t}+B_{t}U^{i}_{t}+D_{t}Z_{t}+W^{i}_{t},\quad
i\in\\{1,\ldots,n\\},$ (1)

where $A_{t}$, $B_{t}$, and $D_{t}$ are matrices of appropriate dimensions and
$\\{W^{i}_{t}\\}_{t=1}^{T}$ is an i.i.d. noise process where $W^{i}_{t}$ is
Gaussian with zero mean and covariance $\Sigma_{W}$. We assume that the
primitive random variables $\\{X^{i}_{1},W^{i}_{1},\dots
W^{i}_{T}\\}_{i=1}^{n}$ are independent.

Thus, the dynamics of the subsystems is coupled to the others through the
mean-field.

####  II-B2 Per-step cost

At time $t$, the system incurs a cost that depends on the local state of the
subsystems and the mean-field as follows:

$c_{t}(\mathbf{X}_{t},\mathbf{U}_{t},Z_{t})=\frac{1}{n}\sum_{i=1}^{n}\Big{[}\mathinner{X^{i}_{t}}^{\intercal}Q_{t}{X^{i}_{t}}+\mathinner{U^{i}_{t}}^{\intercal}R_{t}{U^{i}_{t}}\Big{]}+\mathinner{Z_{t}}^{\intercal}P_{t}{Z_{t}}$
(2)

where $Q_{t}$, $P_{t}$, $S_{t}$, and $R_{t}$ are matrices of appropriate
dimension; $Q_{t}$ and $P_{t}$ are symmetric and non-negative definite and
$R_{t}$ is symmetric and positive definite.

###### Remark 1

The results derived in this paper also apply to a more general per-step cost
of the form

$\frac{1}{n}\sum_{i=1}^{n}\Big{[}\mathinner{X^{i}_{t}}^{\intercal}Q_{t}{X^{i}_{t}}+{X^{i}}^{\intercal}S_{t}Z_{t}+\mathinner{U^{i}_{t}}^{\intercal}R_{t}{U^{i}_{t}}\Big{]}+\mathinner{Z_{t}}^{\intercal}P_{t}{Z_{t}}$

because it can be rewritten in the form (2). □

####  II-B3 Observation model and information structure

We consider two observation models that differ in the observation of the local
state $X^{i}_{t}$ at controller $i$. In the first model, called _full
observation model_ , controller $i$ perfectly observes the local state
$X^{i}_{t}$; in the second model, called _noisy observation model_ ,
controller $i$ observes a noisy version $Y^{i}_{t}\in\mathbb{R}^{d_{y}}$ of
the local state $X^{i}_{t}$ that is given by

$Y^{i}_{t}=C_{t}^{x}X^{i}_{t}+C_{t}^{z}Z_{t}+V^{i}_{t},\quad
i\in\\{1,\ldots,n\\},$ (3)

where the observation noise $\\{V^{i}_{t}\\}_{t=1}^{T}$ is i.i.d. Gaussian
process with zero mean and covariance $\Sigma_{V}$. In addition, we assume
that the noise processes are independent across subsystems and also
independent of $\\{X^{i}_{1},W^{i}_{1},\dots W^{i}_{T}\\}_{i=1}^{n}$.

In both models, in addition to the local measurement of the state of its
subsystem, each controller perfectly observes the mean-field $Z_{t}$.
Following [5], we call this model _mean-field sharing_ information structure.
Controllers perfectly recall all the data they observe. Thus, in the full
observation model, controller $i$ chooses control action according to

$U^{i}_{t}=g^{i}_{t}(X^{i}_{1:t},U^{i}_{1:t-1},Z_{1:t}),$ (4)

while in the noisy observation model, controller $i$ chooses control action
according to

$U^{i}_{t}=g^{i}_{t}(Y^{i}_{1:t},U^{i}_{1:t-1},Z_{1:t}).$ (5)

The function $g^{i}_{t}$ is called the _control law_ of controller $i$. The
collection $\mathbf{g}^{i}=(g^{i}_{1},g^{i}_{2},\ldots,g^{i}_{T})$ is called
the _control strategy_ of controller $i$. The collection
$\mathbf{g}=(\mathbf{g}^{1},\ldots,\mathbf{g}^{n})$ of control strategies of
all controllers is called the _control strategy of the system_.

####  II-B4 The optimization problem

We are interested in the following optimization problem.

###### Problem 1

In the model described above, find a joint strategy
$\mathbf{g}=(\mathbf{g}^{1},\ldots,\mathbf{g}^{n})$ that minimizes the
following cost:

$J(\mathbf{g})=\mathds{E}^{\mathbf{g}}\bigg{[}\sum_{t=1}^{T}c_{t}(\mathbf{X}_{t},\mathbf{U}_{t},Z_{t})\bigg{]},$
(6)

where the expectation is with respect to the measure induced on all system
variables by the choice of the strategy $\mathbf{g}$.

###  II-C The main results

####  II-C1 Full observation model

For the full observation model, our main result is the following:

###### Theorem 1

The optimal control strategy of Problem 1 is unique and given by

${U^{i}_{t}}=K^{x}_{t}X^{i}_{t}+(K^{z}_{t}-K^{x}_{t})Z_{t},\quad
i\in\\{1,\ldots,n\\}.$ (7)

The optimal gains are given as follows. For the terminal time-step
$K^{x}_{T}=K^{z}_{T}=0_{d_{u}\times d_{x}}$ and for $t\leq T-1$,

$K^{x}_{t}=-\left(B_{t}^{\intercal}M^{x}_{t+1}B_{t}+R_{t}\right)^{-1}B_{t}^{\intercal}M^{x}_{t+1}A_{t},$
(8)

and

$K^{z}_{t}=-\left(B_{t}^{\intercal}M^{z}_{t+1}B_{t}+R_{t}\right)^{-1}B_{t}^{\intercal}M^{z}_{t+1}\bar{A}_{t}$
(9)

where $\bar{A}_{t}\coloneqq A_{t}+D_{t}$ and $\\{M^{x}_{t}\\}_{t=1}^{T}$ and
$\\{M^{z}_{t}\\}_{t=1}^{T}$ are solutions of the following two _decoupled_
Riccati equations:

$M^{x}_{T}=Q_{T},\quad M^{z}_{t}=Q_{T}+P_{T},$

and for $t=T-1,\ldots,1$,

$\displaystyle M^{x}_{t}$
$\displaystyle=-A_{t}^{\intercal}M^{x}_{t+1}B_{t}\left(B_{t}^{\intercal}M^{x}_{t+1}B_{t}+R_{t}\right)^{-1}B_{t}^{\intercal}M^{x}_{t+1}A_{t}$
$\displaystyle\quad+A_{t}^{\intercal}M^{x}_{t+1}A_{t}+Q_{t},$ (10)

and

$\displaystyle M^{z}_{t}$
$\displaystyle=-{\bar{A}_{t}}^{\intercal}M^{z}_{t+1}B_{t}\left(B_{t}^{\intercal}M^{z}_{t+1}B_{t}+R_{t}\right)^{-1}B_{t}^{\intercal}M^{z}_{t+1}\bar{A}_{t}$
$\displaystyle\quad+\bar{A}_{t}^{\intercal}M^{z}_{t+1}\bar{A}_{t}+Q_{t}+P_{t}.$
(11)

□

The above result may also be derived under the weaker assumption that the
initial state $\mathbf{X}_{1}=(X^{1}_{1},\dots,X^{n}_{1})$ as well as the
noise $\mathbf{W}_{t}=(W^{1}_{t},\dots,W^{n}_{t})$ are correlated across the
subsystems.

####  II-C2 Noisy observation model

For the noisy observation model, our main result is the following:

###### Theorem 2

For the noisy observation model, the optimal control strategy for Problem 1 is
unique and is given by

$U^{i}_{t}=K^{x}_{t}\hat{X}^{i}_{t}+(K^{z}_{t}-K^{x}_{t})Z_{t},\quad
i\in\\{1,\ldots,n\\},$ (12)

where optimal gains $K^{x}_{t}$ and $K^{z}_{t}$ are the same as in Theorem 1
and $\hat{X}^{i}_{t}=\mathds{E}[X^{i}_{t}\mid
Y^{i}_{1:t},Z_{1:t},U^{i}_{1:t-1}]$ which is generated by the following Kalman
filtering equations: $\hat{X}^{i}_{1}=0$ and for $t>1$,

$\hat{X}^{i}_{t+1}=A_{t}\hat{X}^{i}_{t}+B_{t}U^{i}_{t}+K_{t}^{F}(Y^{i}_{t}-C^{x}_{t}\hat{X}^{i}_{t}-C_{t}^{z}Z_{t})$
(13)

where the Kalman filter gains are given by

$K^{F}_{t}=A_{t}S_{t}{C^{x}_{t}}^{\intercal}(C^{x}_{t}S_{t}{C^{x}_{t}}^{\intercal}+\Sigma_{V})^{-1},$
(14)

and the state estimation error covariances satisfy the (filter) Riccati
equation: $S_{1}=\Sigma_{X}$ and for $t>1$,

$\displaystyle S_{t+1}$
$\displaystyle=-A_{t}S_{t}{C^{x}_{t}}^{\intercal}(C^{x}_{t}S_{t}{C^{x}_{t}}^{\intercal}+\Sigma_{V})^{-1}C^{x}_{t}S_{t}A_{t}^{\intercal}$
$\displaystyle\quad+A_{t}S_{t}A_{t}^{\intercal}+\Sigma_{W}.$ (15)

□

##  III Proof of the main results

The main idea of the proof is as follows. We construct an auxiliary system
whose state, control actions, and per-step cost are equivalent to
$\mathbf{X}_{t}$, $\mathbf{U}_{t}$, and $c_{t}(\cdot)$, respectively (modulo a
change of variables to be described later). However, this auxiliary system is
centrally controlled by a single control that has access to all the
information available to the $n$ decentralized controllers in the original
system. We show that the optimal centralized solution of this auxiliary system
can be implemented in the original decentralized system, and is therefore also
optimal for the decentralized system.

###  III-A The auxiliary system

Define $\bar{X}^{i}_{t}=X^{i}_{t}-Z_{t}$ and
$\bar{U}^{i}_{t}=U^{i}_{t}-U^{z}_{t}$, where
$U^{z}_{t}=\langle\mathbf{U}_{t}\rangle$. The auxiliary system is a
centralized system with state
$\tilde{X}_{t}=\operatorname{vec}(\bar{X}^{1}_{t},\dots,\bar{X}^{n}_{t},Z_{t})$
and actions
$\tilde{U}_{t}=\operatorname{vec}(\bar{U}^{1}_{t},\dots,\bar{U}^{n}_{t},U^{z}_{t})$.
Note that $\tilde{X}_{t}$ is equivalent to $\mathbf{X}_{t}$ and
$\tilde{U}_{t}$ is equivalent to $\mathbf{U}_{t}$.

The dynamics are the same as the model in Sec. II-B. In particular,

$\bar{X}^{i}_{t+1}=A_{t}\bar{X}^{i}_{t}+B_{t}\bar{U}^{i}_{t}+W^{i}_{t}-\langle\mathbf{W}_{t}\rangle$
(16)

and

$Z_{t+1}=(A_{t}+D_{t})Z_{t}+B_{t}U^{z}_{t}+\langle\mathbf{W}_{t}\rangle$ (17)

The per-step cost of the auxiliary model is given by
$c_{t}(\mathbf{X}_{t},\mathbf{U}_{t},Z_{t})$ defined in (2).

As in Sec. II-B, we consider two observation models for the auxiliary system:
full observation and noisy observation. In both cases, there is a _single
centralized controller_ that chooses $\tilde{U}_{t}$ based on the
observations.

In the full observation model, the centralized controller observes
$\tilde{X}_{t}$ and chooses $\tilde{U}_{t}$ according to

$\tilde{U}_{t}=\bar{g}_{t}(\tilde{X}_{1:t},\tilde{U}_{1:t-1}).$

In the noisy observation model, the centralized controller observes
$\mathbf{Y}_{t}=\operatorname{vec}(Y^{1}_{t},\dots,Y^{N}_{t})$ where
$Y^{i}_{t}$ is given by (3) and chooses $\tilde{U}_{t}$ according to

$\tilde{U}_{t}=\bar{g}_{t}(\mathbf{Y}_{1:t},\tilde{U}_{1:t-1},Z_{1:t}).$

In both models, we are interested in the following optimization problem:

###### Problem 2

In the model described above, find a strategy
$\bar{\mathbf{g}}=(\bar{g}_{1},\dots,\bar{g}_{T})$ that minimizes the
following cost:

$\bar{J}(\bar{\mathbf{g}})=\mathds{E}^{\bar{\mathbf{g}}}\Big{[}\sum_{t=1}^{T}c_{t}(\mathbf{X}_{t},\mathbf{U}_{t},Z_{t})\Big{]},$
(18)

where the expectation is with respect to the measure induced on all system
variables by the choice of strategy $\bar{\mathbf{g}}$.

Let $J^{*}$ and $\bar{J}^{*}$ denote the optimal cost for Problem 1 and
Problem 2 respectively. Since the per-step cost is the same in both cases, but
Problem 2 is centralized, we have that $J^{*}\geq\bar{J}^{*}.$ We identify the
optimal control laws for the auxiliary system and show that these laws can be
implemented in, and therefore are optimal for, the original decentralized
system.

A critical step in the proof is to rewrite the cost
$c_{t}(\mathbf{X}_{t},\mathbf{U}_{t},Z_{t})$ in terms of $\tilde{X}_{t}$ and
$\tilde{U}_{t}$. For that matter, we need the following key result.

###### Lemma 1

For any $\mathbf{x}=\operatorname{vec}(x^{1},\dots,x^{n})$ and
$z=\langle\mathbf{x}\rangle$, let $\bar{x}^{i}=x^{i}-z$,
$i\in\\{1,\dots,n\\}$. Then, for any matrix $Q$ of appropriate dimension,

$\frac{1}{n}\sum_{i=1}^{n}\mathinner{x^{i}}^{\intercal}Q{x^{i}}=\frac{1}{n}\sum_{i=1}^{n}\mathinner{\bar{x}^{i}}^{\intercal}Q{\bar{x}^{i}}+\mathinner{z}^{\intercal}Q{z}$
(19)

□

###### Proof

The result follows from elementary algebra and the observation that
$\sum_{i=1}^{n}\bar{x}^{i}=0$. ■

An immediate consequence of Lemma 1 is the following:

###### Corollary 1

For any time $t$,
$c_{t}(\mathbf{X}_{t},\mathbf{U}_{t},Z_{t})=\bar{c}_{t}(\tilde{X}_{t},\tilde{U}_{t})$,
where

$\bar{c}_{t}(\tilde{X}_{t},\tilde{U}_{t})=\frac{1}{n}\sum_{i=1}^{n}\Big{[}\mathinner{\bar{X}^{i}_{t}}^{\intercal}Q_{t}{\bar{X}^{i}_{t}}+\mathinner{\bar{U}^{i}_{t}}^{\intercal}R_{t}{\bar{U}^{i}_{t}}\Big{]}\\\
+\mathinner{Z_{t}}^{\intercal}(Q_{t}+P_{t}){Z_{t}}+\mathinner{U^{z}_{t}}^{\intercal}R_{t}{U^{z}_{t}}.$
(20)

□

Note that the auxiliary model has linear dynamics and in Corollary 1 we have
shown that the cost is quadratic in the state and the control actions. Thus,
the optimal control actions are linear in the state and the corresponding
optimal gains can be obtained by solving an appropriate Riccati equation.
However, the state $\tilde{X}_{t}$ of the auxiliary system belongs to
$\mathds{R}^{(n+1)d_{x}}$, thus, a naive attempt to obtain an optimal solution
will involve solving for $(n+1)d_{x}\times(n+1)d_{x}$ dimensional Riccati
equations. We present an alternative approach in the next section that
involves solving two $d_{x}\times d_{x}$ dimensional Riccati equations
(independent of $n$).

###  III-B Full observation model

The auxiliary system is a stochastic linear quadratic system. From the
certainty equivalence principle [15], we know that the optimal control law is
unique and identical to the control law in the corresponding deterministic
problem, whose dynamics are given by

$\displaystyle\bar{X}^{i}_{t+1}$
$\displaystyle=A_{t}\bar{X}^{i}_{t}+B_{t}\bar{U}^{i}_{t},\qquad
i\in\\{1,\ldots,n\\},$ (21) $\displaystyle Z_{t+1}$
$\displaystyle=(A_{t}+D_{t})Z_{t}+B_{t}U^{z}_{t}.$ (22)

and the per-step cost is $\bar{c}(\tilde{X}_{t},\tilde{U}_{t})$ given by (20).

Note that this system consists on $(n+1)$ components: $n$ components with
state $\bar{X}^{i}_{t}$ and control $\bar{U}^{i}_{t}$, $i\in\\{1,\dots,n\\}$,
and one component with state $Z_{t}$ and control $U^{z}_{t}$. These components
have decoupled dynamics and decoupled cost; and $n$ of these are identical.
Thus, the optimal control law of each component may be identified separately.
Therefore, we have the following:

###### Theorem 3

The optimal control strategy of auxiliary model is unique and given by

$\bar{U}^{i}_{t}=K^{x}_{t}\bar{X}^{i}_{t},\quad U^{z}_{t}=K^{z}_{t}Z_{t},\quad
i\in\\{1,\ldots,N\\},$ (23)

where the gains $K^{x}_{t}$ and $K^{z}_{t}$ are given as in Theorem 1. □

To complete the proof of Theorem 1, note that

$U^{i}_{t}=\bar{U}^{i}_{t}+U^{z}_{t}=K^{x}(X^{i}_{t}-Z_{t})+K^{z}Z_{t}.$

Thus, the control laws specified in Theorem 1 are the optimal _centralized_
control laws, and, a fortrori, the optimal decentralized control laws.

###  III-C Noisy observation model

Figure 1: An example of mean-field control of space heaters. The light lines
show the temperature of individual heaters and the thick line depicts the
mean-field of the trajectories. Note that the individual trajectories adjust
around their initial values so that the mean-field tracks the reference
trajectory $z_{ref}=25$. The tracking is not perfect because of the friction
between tracking the initial states locally and tracking the mean-field
globally.

Define

$\displaystyle\mathring{X}^{i}_{t}$
$\displaystyle=\mathds{E}[\bar{X}^{i}_{t}\mid\mathbf{Y}_{1:t},\tilde{U}_{1:t-1},Z_{1:t}],$
$\displaystyle\breve{X}^{i}_{t}$
$\displaystyle=\mathds{E}[X^{i}_{t}\mid\mathbf{Y}_{1:t},\tilde{U}_{1:t-1},Z_{1:t}].$

Since $X^{i}_{t}=\bar{X}^{i}_{t}+Z_{t}$, we have that
$\breve{X}^{i}_{t}=\mathring{X}^{i}_{t}+Z_{t}$.

From the Separation Theorem for linear quadratic Gaussian systems [15], we
know that the optimal (centralized) control law for the auxiliary system is
given by

$\bar{U}^{i}_{t}=K^{x}_{t}\mathring{X}^{i}_{t},\quad
U^{z}_{t}=K^{z}_{t}Z_{t},\quad i\in\\{1,\ldots,N\\}.$

Or equivalently,

$U^{i}_{t}=K^{x}_{t}\breve{X}^{i}_{t}+(K^{z}_{t}-K^{x}_{t})Z_{t},\quad
i\in\\{1,\ldots,n\\}.$

Thus, to prove Theorem 2, we need to show that
$\breve{X}^{i}_{t}=\hat{X}^{i}_{t}$ (where $\hat{X}^{i}_{t}$ is defined in
Theorem 2). This is proved below.

###### Lemma 2

In the auxiliary system with noisy observation,
$\breve{X}^{i}_{t}=\hat{X}^{i}_{t}$. □

###### Proof

Recall that $X^{i}_{t}$ evolves according to (1) and $Y^{i}_{t}$ is given by
(3). Therefore, for any $i,j\in\\{1,\dots,n\\}$ such that $i\neq j$,
$(X^{i}_{1:t},Y^{i}_{1:t})$ is conditionally independent of
$(X^{j}_{1:t},Y^{j}_{1:t})$. Consequently,

$\displaystyle\breve{X}^{i}_{t}$
$\displaystyle=\mathds{E}[X^{i}_{t}\mid\mathbf{Y}_{1:t},\tilde{U}_{1:t-1},Z_{1:t}]$
$\displaystyle=\mathds{E}[X^{i}_{t}\mid
Y^{i}_{1:t},\tilde{U}_{1:t-1},Z_{1:t}]$
$\displaystyle=\mathds{E}[X^{i}_{t}\mid
Y^{i}_{1:t},\mathbf{U}_{1:t-1},Z_{1:t}]$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mathds{E}[X^{i}_{t}\mid
Y^{i}_{1:t},U^{i}_{1:t-1},Z_{1:t}]$ $\displaystyle=\hat{X}^{i}_{t},$

where $(a)$ follows because $X^{i}_{t}$ depends on $\mathbf{U}_{1:t-1}$ only
through $U^{i}_{1:t-1}$.

■

##  IV Numerical Example

We simulate an example of inspired by the _collective target tracking mean-
field model_ of [16]. We consider a population of space heaters. The state
$X^{i}_{t}$ denotes the room temperature at heater $i$ at time $t$. We assume
that the inital mean temperature is $z_{1}$ degrees, with individual
temperatures distributed according to Gaussian distribution with unit variance
around $z_{1}$ degrees. We assume that the temperature dynamics are linearized
around the operating point and are given by

$X^{i}_{t+1}=a(X^{i}_{t}-x_{0})+bU^{i}_{t}+W^{i}_{t}$

where $x_{0}$ is the ambient temperature, $W^{i}_{t}$ is the randomness due to
the environment and $U^{i}_{t}$ is the control action of a local controller.

It is desired that the mean temperature increases to 21 degrees (denoted by
$z_{ref}$) in $T$ time steps (we assume $T=90$ minutes). The per-step cost
function is

$\frac{1}{n}\sum_{i=1}^{n}\Big{[}q(X^{i}_{t}-X^{i}_{1})^{2}+r{U^{i}_{t}}^{2}\Big{]}+p(Z_{t}-z_{ref})^{2}$

The rationale of this cost function is as follows. It is assumed that the
inital temperature is the comfort level of user $i$; so we penalize local
deviations from the initial temperature. The second term corresponds to the
local control energy. The last term is the penalty for the mean-field
deviating from the reference mean-field. The objective is to minimize the
expected total cost over a finite horizon.

The above model is a tracking problem. The results of Theorem 1 extend to
tracking problems in a natural manner. We solve the Riccati equations for the
following values of the parameters:

$\displaystyle n$ $\displaystyle=30,$ $\displaystyle a$ $\displaystyle=0.8,$
$\displaystyle b$ $\displaystyle=1,$ $\displaystyle q$ $\displaystyle=0.5,$
$\displaystyle p$ $\displaystyle=1,$ $\displaystyle r$ $\displaystyle=1,$
$\displaystyle x_{0}$ $\displaystyle=22,$ $\displaystyle z_{1}$
$\displaystyle=22,$ $\displaystyle z_{ref}$ $\displaystyle=25,$ $\displaystyle
T$ $\displaystyle=90,$ $\displaystyle W^{i}_{t}$
$\displaystyle\sim\mathcal{N}(0,1),$ $\displaystyle X^{i}_{1}$
$\displaystyle\sim\mathcal{N}(22,2).$

The resultant trajectories are shown in Fig. 1.

##  V Conclusion

We consider team-optimal control of finite number of mean-field coupled LQG
subsystems with mean-field sharing information structure. The optimal control
law is unique, identical for all controllers, and linear in the current state
and the current mean-field; the optimal gains are obtained by solving two
decoupled Riccati equations for the full observation case and an additional
filter Riccati equation for the noisy observation case. These Riccati
equations do not depend on the number of subsystems.

To prove the main results, we consider an auxiliary centralized system whose
dynamics and cost are equivalent to the original decentralized system. It is
shown that the optimal centralized control law are implementable with mean-
field sharing. Thus, the optimal decentralized control laws are identical to
the optimal centralized control laws.

Although we only presented the results for the finite horizon optimal
regulation problem, the results extend to infinte horizon and to optimal
tracking problems in a natural manner. Moreover, it is possible to extend the
results to non-homogeneous population consisting of multiple types.

## References

  * [1] M. Aicardi, F. Davoli, and R. Minciardi, “Decentralized optimal control of Markov chains with a common past information set,” _IEEE Trans. Autom. Control_ , vol. 32, no. 11, pp. 1028–1031, 1987.
  * [2] A. Nayyar, A. Mahajan, and D. Teneketzis, “Optimal control strategies in delayed sharing information structures,” _IEEE Trans. Autom. Control_ , vol. 56, no. 7, pp. 1606–1620, Jul. 2011.
  * [3] J. M. Ooi, S. M. Verbout, J. T. Ludwig, and G. W. Wornell, “A separation theorem for periodic sharing information patterns in decentralized control,” _IEEE Trans. Autom. Control_ , vol. 42, no. 11, pp. 1546–1550, Nov. 1997\.
  * [4] S. Yüksel, “Stochastic nestedness and the belief sharing information pattern,” _IEEE Trans. Autom. Control_ , pp. 2773–2786, Dec. 2009.
  * [5] J. Arabneydi and A. Mahajan, “Team optimal control of coupled subsystems with mean field sharing,” in _Proc. 53rd IEEE Conf. Decision and Control_ , Los Angeles, CA, Dec. 2014, pp. 1669–1674.
  * [6] A. Nayyar, A. Mahajan, and D. Teneketzis, “Decentralized stochastic control with partial history sharing: A common information approach,” _IEEE Trans. Autom. Control_ , vol. 58, no. 7, pp. 1644–1658, Jul. 2013.
  * [7] A. Mahajan, N. Martins, M. Rotkowitz, and S. Yüksel, “Information structures in optimal decentralized control,” in _Proc. 51st IEEE Conf. Decision and Control_ , Maui, Hawaii, Dec. 2012, pp. 1291 – 1306.
  * [8] M. Huang, P. Caines, and R. Malhame, “Large population cost-coupled LQG problems with non-uniform agents: individual mass behavior and decentralized epsilon-nash equilibrium,” _IEEE Trans. Autom. Control_ , vol. 52, no. 9, pp. 1560–1571, Sep. 2007.
  * [9] ——, “An invariance principle in large population stochastic dynamic games,” _Journal of Systems Science and Complexity_ , vol. 20, pp. 162–172, 2007.
  * [10] P. Lions and J. Larsy, “Mean-field games,” _Japanese Journal of Mathematics_ , vol. 2, no. 1, pp. 229–260, Mar. 2007.
  * [11] T. Li and J.-F. Zhang, “Asymptotically optimal decentralized control for large population stochastic multiagent systems,” _IEEE Trans. Autom. Control_ , vol. 53, no. 7, pp. 1643–1660, 2008.
  * [12] O. Guéant, J.-M. Lasry, and P.-L. Lions, _Paris-Princeton Lecture on Mathematical Finance_. Springer, 2011, ch. Mean-field games and applications, pp. 205–266.
  * [13] M. Huang, P. E. Caines, and R. P. Malhamé, “Social optima in mean field LQG control: centralized and decentralized strategies,” _IEEE Trans. Autom. Control_ , vol. 57, no. 7, pp. 1736–1751, 2012.
  * [14] R. Elliott, X. Li, and Y.-H. Ni, “Discrete time mean-field stochastic linear-quadratic optimal control problems,” _Automatica_ , vol. 49, no. 11, pp. 3222–3233, 2013.
  * [15] P. E. Caines, _Linear stochastic systems_. John Wiley & Sons, Inc., 1987.
  * [16] A. Kizilkale and R. Malhame, “Mean field based control of power system dispersed energy storage devices for peak load relief,” in _IEEE Conference on Decision and Control (CDC)_ , Dec. 2013, pp. 4971–4976.


# Local Effects of Grid-Forming Converters Providing Frequency Regulation
to Bulk Power Grids

Antonio Zecchino, Francesco Gerini, Yihui Zuo,
Rachid Cherkaoui, Mario Paolone  Distributed Electrical Systems Laboratory
(DESL)
Ecole Polytechnique Fédérale de Lausanne, EPFL
Lausanne, Switzerland
{antonio.zecchino, francesco.gerini, yihui.zuo,
rachid.cherkaoui<EMAIL_ADDRESS>Elena Vagnoni This work is supported
by OSMOSE project. The project has received funding from the European Union’s
Horizon 2020 Research and Innovation Program under Grant 773406. Technology
Platform for Hydraulic Machines
Ecole Polytechnique Fédérale de Lausanne, EPFL
Lausanne, Switzerland
<EMAIL_ADDRESS>
###### Abstract

The progressive displacing of conventional generation in favor of renewable
energy sources requires restoring an adequate capacity of regulating power to
ensure reliable operation of power systems. _Battery Energy Storage Systems_
(BESSs) are considered to be promising assets to restore suitable frequency
regulation capacity levels. BESSs are typically connected to the grid with
power-converters, able to operate in either grid-forming or grid-following
modes. This paper quantitatively assesses the impact on the local distribution
grid of BESSs providing frequency regulation to bulk power systems. Specific
metrics are proposed to compare the performance of grid-forming and grid-
following control. Experimental results are obtained taking advantage of a 720
kVA/500 kWh BESS connected to the 20 kV distribution grid of the EPFL campus.
The quantitative evaluation based on suitably proposed metrics confirms the
superior performance of grid-forming strategy, compared to grid-following one.

##  I Introduction

Integrating renewable energy sources in the European energy supply mix has
long been ascertained as a fundamental step for the decarbonisation of the
electric sector. Moving towards this direction raises multiple challenges for
power systems, including the need for regulating power and flexibility
services provision to mitigate the massive integration of non-dispatchable and
stochastic resources. In this respect, converter-interfaced _Battery Energy
Storage Systems_ (BESSs) are considered as perspective assests for grid
frequency regulation in view of their large ramping rates, high round-trip
efficiency and commercial availability [1, 2, 3]. BESSs interface with the
grid through four-quadrant power converters [4], typically controlled to
provide grid ancillary services ranging from frequency regulation up to energy
management [5]. As known, there are two main approaches to control converter-
interfaced BESSs: _Grid-Forming_ (GFR) and _Grid-Following_ (GFL). Nowadays,
the majority of converter-interfaced resources are controlled as GFL sources,
as this operating mode is considered to be sufficiently efficient [6, 7].
Nonetheless, future grids may require GFR devices to provide more frequency
and voltage regulation, stability and black-start capabilities, all services
that are nowadays offered primarily by synchronous machines [8, 9].
Furthermore, recent studies have proved GFR control strategies to outperform
GFL in terms of frequency regulation performance in low-inertia power grids
[10]. The impact of GFR converters on the dynamics of a reduced-inertia grid
has been investigated in [11], which quantitatively proved the good
performance of GFR units in limiting the frequency deviation and in damping
the frequency oscillations in case of large contingencies. The existing
scientific literature lacks of studies on the assessment of the performance of
GFR units in supporting the frequency containment process of large
interconnected power grids. In fact, studies on the GFR units synchronizing
with AC power grids are mostly limited to either simulation [11, 10, 12] or to
experiments on ideal slack buses with emulated voltage [13, 14]. The latter
are capable of identifying the fundamental characteristics of the GFR
controls, as well as to optimize optimal control parameters [13]. However,
they can not assess the dynamic interaction between power grids and the GFR
units. Simulations-based studies have been investigating the impact of large-
scale GFR or GFL units on frequency containment of low-inertia power grids[10,
12], where the power rating of the controlled element (i.e. the BESS) is
comparable to the one of the power system. Nevertheless, in bulk power grids,
the influence of a single frequency control provider can hardly be detected
through frequency measures at the transmission level. An alternative approach
is to measure the influence of units providing frequency regulation services
on the local distribution grid. In this respect, the scope of this paper is to
experimentally assess the performance of GFR and GFL control for converter-
interfaced BESS by measuring with accurate _Phase Measurement Units_ (PMUs)
the influence of the two control strategies on the local distribution grid.
The experimental campaign is conducted by taking advantage of a grid-connected
720 kVA/500 kWh BESS connected to the 21 kV distribution feeder on the EPFL
campus.

The paper is structured as follows: Section II presents the characteristics
GFR and GFL, Section III analyses the state-of-the-art metrics used to compare
GFR and GFL control and proposes a new metric, to better capture local effects
of the two control strategies. Section IV presents and discusses the
experimental results. Finally, Section V summarizes the results and concludes
the paper.

##  II Converter Controls

As mentioned in Section I, GFR and GFL controls are the two approaches
currently used to control converter-interfaced BESSs [8, 15]. For the sake of
consistency, this section reports the description of GFR and GFL controller
proposed by [8].

GFR units are based on a power converter which controls magnitude and angle of
the voltage at the _Point of Battery Coupling_ (PBC). As a consequence, the
knowledge of the fundamental frequency phasor of the grid voltage at the PBC
is not strictly necessary. Fig. 1a presents a general structure of the GFR
control, where the modulated voltage angle is linked with the converter active
power, while the voltage magnitude is regulated by taking into account of
voltage or/and reactive power reference.

GFL units are based on a power converter which controls the currents injection
to adjust the active and reactive power injection. To this end, the injected
current is controlled with a specific phase displacement with respect to the
grid voltage at the PBC. Consequently, the fundamental frequency phasor of the
grid voltage at the PBC is needed at any time to correctly calculate the
converter reference currents. Fig. 1b shows the classical structure of the GFL
control, where the grid voltage angle $\tilde{\theta}_{g}$ is estimated thanks
to a _Phase Locked-Loop_ (PLL) and used by the Park’s transformation. The
outer control loops can adjust the amplitude and angle of the current
reference currents to inject required amount of active and reactive power.
Grid-supporting GFL converters can be applied to adjust active and reactive
power set points under frequency or voltage deviations

(a) Main scheme of a GFR control.

(b) Main scheme of a GFL control.

Figure 1: Main schemes of GFR and GFL controls [8].

##  III Quantification of Frequency Regulation

###  III-A State-of-the-Art Metrics

To evaluate the impact of the frequency regulation provided by BESSs to low-
inertia grids, [10] proposes the hereby listed metrics:

  * •

Frequency _Probability Density Function_ (PDF), where the standard deviation
of the grid-frequency statistical distribution is observed to assess the
capability of the control strategy in limiting the frequency deviations from
the nominal value.

  * •

_Integral Frequency Deviation_ (IFD):

$\displaystyle\text{IFD}=\sum_{i}^{M}\sum_{k=1}^{N}{|f_{k,i}-f_{0}|}$ (1)

where $M$ is the number of measurement units, $N$ is the total sampling number
of frequency measurements for each load, $f_{k,i}$ is the frequency measured
by the measurement unit $i$ at sample $k$ and $f_{0}$ is the nominal
frequency, i.e. 50 Hz.

Both frequency PDF and IFD establish a relation between grid-frequency
measurements and control actions of the converter-interfaced BESS.
Nevertheless, it is not possible to measure the influence of a small frequency
control provider to the bulk power grid, as the frequency can be considered
independent from the single action of a small power converter. Therefore,
frequency PDF and IFD can not be considered as reliable metrics.

###  III-B New metrics

To characterize the local effect of a converter-interfaced BESS providing
frequency regulation and to identify the difference between the GFR and GRL
controls, we proposed new metrics based on PMUs measurements performed
locally, i.e., on the distribution feeder where the BESS is installed. The
frequency measures at the _Point of Common Coupling_ (PCC) of the distribution
feeder, as well as the phase angle differences between the voltage controlled
by the converter unit and the PCC voltage are taken into consideration to
derive the following metrics:

  * •

_Relative Rate-of-Change-of-Frequency_ (rRoCoF):

$\displaystyle rRoCoF=\left|\frac{\Delta f_{PCC}/\Delta t}{\Delta P}\right|$
(2)

where $\Delta f_{PCC}$ is the difference between one grid frequency sample and
the next (once-differentiated value) at the bus where the BESS is connected
to, $\Delta P$ is the once-differentiated BESS active power, and $\Delta t$ is
the sampling interval.

  * •

_Relative Phase Angle Difference Deviation_ (rPADD):

$rPADD_{k}=\left|\frac{\Delta\theta_{k}-\Delta\theta_{0}}{\Delta
P_{k}}\right|$ (3)

quantifying the change in the phase-to-neutral voltage angle difference
$\Delta\theta_{k}$, measured by two PMUs installed on nodes at different
voltage level of the local feeder, versus the case with null delivered active
power $\Delta\theta_{0}$.

The relative RoCoF metric measures the RoCoF regulation at PCC. Since the
RoCoF is weighted by the delivered active power of the BESS, it can be used to
compare the local performance of GFR versus GFL converter in large
interconnected electrical systems. On the other side, for the rPADD
computation, the benefit of using metrics based on voltage phase angle is two-
fold:  (i) changes in the phase angle estimation are directly linked with
frequency deviations from nominal frequency [16] and  (ii) it can indicate the
ability of the connected unit in sustaining local voltage by controlling its
angle.  The metric is normalized by the delivered active power of the
converter-interfaced unit to ensure its independence from the grid condition
and the power rating of the converter-interfaced unit.

##  IV Experimental campaign

###  IV-A Experimental Setup

####  IV-A1 Feeder Characteristics

for the experimental campaign, a 20 kV distribution network equipped with a
BESS is considered. The experimental configuration shown in Fig. 2 consists in
a group of buildings of the EPFL campus characterised by a 140 kW base load,
hosting 95 kWp root-top PV installation and a grid-connected 720 kVA/500 kWh
Lithium Titanate BESS. The BESS bidirectional real power flow is denoted by
$P$, and $G$ is the composite power flow seen at the PCC. The aggregated
building demand is denoted by $L$ and, by neglecting grid losses [17] , it is
estimated as $L=P-G$. For both buildings and BESS, positive values denote
injection towards the PCC, while negative denote consumption. The measuring
systems is composed by three PMUs, i.e. PMU 2, PMU 1 and PMU 0, installed
respectively on 0.3 kV, 20 kV and 50 kV side of the feeder. Two transformers
connect the different voltage levels of the feeder. Their nominal voltage and
power, as well as their short circuit voltage expressed in percent are shown
in Fig. 2.

Figure 2: Experimental setup of the EPFL campus MV grid hosting the utility-
scale 720kVA/500 kWh BESS.

####  IV-A2 Controller Characteristics

the BESS converter is capable of synchronizing with the grid in both GFL and
GFR control modes. In GFR mode the converter controls the voltage angle to
provide frequency regulation by means of an internal control loop directly
establishing a droop relation between the grid frequency and the active power
to be delivered. In GFL mode, the active power set-point for frequency
regulation is computed in an external loop, as a function of the grid
frequency deviation (from 50 Hz). In particular, the maximum feasible droop,
i.e. 1.44 MW/Hz, is selected for the experiments and implemented in both GFR
and GFL modes to maximise the BESS response.

####  IV-A3 Requirement for PMUs

for the computation of the Relative Integrated Phase Angle Difference
Deviation (rPADD), high-accuracy PMUs are required. Knowing the feeder
topology, the short-circuit voltage of the two transformers and assuming the
line as lossless, it is possible to estimate the angle displacement per unit
of delivered power by the BESS as 0.006 degree/kW. Considering a droop of 1.44
MW/Hz and the need of estimating BESS reactions for frequency deviations
greater than the typical frequency deadband [18], i.e. 10 mHz [18], the
measuring infrastructure must be able to detect angle differences of 0.9
degree. To robustly measure those phase angle differences, the PMU-based
distributed sensing with infrastructure deployed in the EPFL campus have been
utilized [17], thanks to its accuracy in terms of 1 standard deviation equal
to $0.001$ degrees (i.e. 18 $\mu$rad), 90 times smaller than the minimum angle
required for the rPADD computation.

####  IV-A4 Grid condition

as the experimental validation relies on real operating conditions, it is
subjected to the uncontrollable variation of the power system frequency of the
Continental Europe Synchronous Region, to which the medium-voltage grid of the
EPFL campus is connected. For this reason, the experimental campaign has taken
place during the transition of the hour, which typically present fast
frequency ramps due to mechanisms related to electricity markets [19]. Such
frequency ramps have been used to activate BESS to perform considerable active
power exchanges $P$, following the 1.44 MW/Hz droop implemented both in the
GFR and GFL controls. Furthermore, the proposed tests have been carried out
when the base consumption of the local medium voltage $L$ grid presents a
mostly constant profile such that its dynamics are not influencing the local
measures on which the metric are based, i.e., during evening and night hours.
This allows for assuming $\Delta\theta_{0}$ to be constant during the
experiment, enabling a better identification of the impact of the BESS
provided power in the local phase angle displacements.

###  IV-B Experimental Results

Figure 3: Experimental test results using the EPFL BESS converter in GFR mode
with 1.44 MW/Hz droop. Figure 4: Experimental test results using the EPFL BESS
converter in GFL mode with 1.44 MW/Hz droop

The experimental assessment of the BESS response are shown in Fig. 3. At the
top, the frequency measured by the PMUs and the BESS active power response are
provided. The frequency ramp at the hour transition (21:00 CET) are noticeably
reflecting on the BESS power ramp. The angle displacement is computed as the
angle difference between PMU 1 (21 kV) and PMU 2 (0.3 kV):

$\Delta\theta_{k}=\theta_{k,PMU1}-\theta_{k,PMU2}$ (4)

Before the activation of the BESS control at 20:50 CET, the phase angle
difference $\Delta\theta_{0}$ is computed as:

$\Delta\theta_{0}=\theta_{0,PMU1}-\theta_{0,PMU2}$ (5)

The values of $\Delta\theta_{0}$ and $\Delta\theta_{k}$ can be inserted in Eq.
3 for the rPADD computation. For the sake of comparison, a similar
experimental assessment is conducted in GFL mode and the results are shown in
Fig. 4.

To quantitatively assess and compare the two cases with respect to the impact
on the local grid frequency, the two relative metrics presented in Section
III, namely rRoCoF and rPADD, are considered. rRoCoF is computed over a time
window of 60 ms. For the sake of clarity and simplicity of the interpretation,
results of the rRoCoF the rPADD computation are reported in terms of
_Cumulative Distribution Function_ (CDF). Fig. 5 shows the rRoCoF CFD for
multiple tests, expressed in [Hz/s/kW]. It is important to observe the
vertical trend of the red line, representing the GFR case, implying higher
probabilities of having smaller values of rRoCoF. This means that, for a
higher number of instances for which the rRoCoF has been computed, it is more
likely to have smaller values when controlling the BESS in GFR mode than in
GFL mode. This demonstrates that the GFR control can guarantee larger
containment of the rate of change of the frequency, compared to the more
traditional GFL strategy. These results are in line with what is presented in
[10], relatively to the simulation activities on the IEEE 39-bus benchmark
grid. Fig. 6 shows the CDF probability of the rPADD of having a given value
expressed in [deg/kW]. It is clear that larger values of the index mean larger
impact on the local phase angle displacement of the voltages measured at the
two different voltage levels. It can be noticed that in this case larger
values are found for the GFR case, demonstrating that the implemented GFR
control can guarantee larger impact on the local phase angle difference,
compared to the more traditional GFL strategy.

Figure 5: Cumulative distribution function of the rRoCoF [Hz/s/W] for the
experimental tests. Figure 6: Cumulative distribution function of the rPADD
[deg/kW] for the experimental tests.

##  V Conclusion

A novel metric to quantitatively assess the local impact of BESSs providing
frequency regulation to bulk power grid has been proposed. Experimental tests
confirmed the outcome of the analysis carried out by means of simulation
activities on the low-inertia configuration of the IEEE 39-bus emulated grid
presented in [10], validating the positive effects of the grid-forming
strategy in the control of the local frequency. The upgrade of the PMU-based
sensing infrastructure, enabled the assessment of impacts on the local
frequency with the BESS control running with maximum possible frequency droop.
The proposed metrics (rROCOF and rPADD) can potentially be of use for the
quantification of frequency reserve provision to distribution systems.
Moreover, both metrics confirmed the superior performance of the grid-forming
strategy, characterised by a more contained rROCOF and a larger impact in the
local phase angle compared to grid-following mode.

## References

  * [1] B. Dunn, H. Kamath, and J.-M. Tarascon, “Electrical Energy Storage for the Grid: A Battery of Choices,” _Science_ , vol. 334, no. 6058, pp. 928–935, Nov. 2011.
  * [2] F. Sossan, E. Namor, R. Cherkaoui, and M. Paolone, “Achieving the Dispatchability of Distribution Feeders Through Prosumers Data Driven Forecasting and Model Predictive Control of Electrochemical Storage,” _IEEE Transactions on Sustainable Energy_ , Oct. 2016.
  * [3] P. Mercier, R. Cherkaoui, and A. Oudalov, “Optimizing a Battery Energy Storage System for Frequency Control Application in an Isolated Power System,” _IEEE Transactions on Power Systems_ , vol. 24, no. 3, pp. 1469–1477, Aug. 2009.
  * [4] G. Wang, G. Konstantinou, C. D. Townsend, J. Pou, S. Vazquez, G. D. Demetriades, and V. G. Agelidis, “A Review of Power Electronics for Grid Connection of Utility-Scale Battery Energy Storage Systems,” _IEEE Transactions on Sustainable Energy_ , vol. 7, no. 4, pp. 1778–1790, Oct. 2016.
  * [5] E. Namor and R. Cherkaoui, “Control of Battery Storage Systems for the Simultaneous Provision of Multiple Services,” _IEEE TRANSACTIONS ON SMART GRID_ , vol. 10, no. 3, p. 10, 2019.
  * [6] J. Beerten, O. Gomis-Bellmunt, X. Guillaud, J. Rimez, A. van der Meer, and D. Van Hertem, “Modeling and control of HVDC grids: A key challenge for the future power system,” in _2014 Power Systems Computation Conference_. Wrocław, Poland: IEEE, Aug. 2014, pp. 1–21.
  * [7] J. A. Baroudi, V. Dinavahi, and A. M. Knight, “A review of power converter topologies for wind generators,” _Renewable Energy_ , vol. 32, no. 14, pp. 2369–2385, Nov. 2007.
  * [8] M. Paolone, T. Gaunt, X. Guillaud, M. Liserre, S. Meliopoulos, A. Monti, T. Van Cutsem, V. Vittal, and C. Vournas, “Fundamentals of power systems modelling in the presence of converter-interfaced generation,” _Electric Power Systems Research_ , vol. 189, p. 106811, Dec. 2020.
  * [9] F. Milano, F. Dorfler, G. Hug, D. J. Hill, and G. Verbic, “Foundations and Challenges of Low-Inertia Systems (Invited Paper),” in _2018 Power Systems Computation Conference (PSCC)_. Dublin: IEEE, Jun. 2018, pp. 1–25.
  * [10] Y. Zuo, Z. Yuan, F. Sossan, A. Zecchino, R. Cherkaoui, and M. Paolone, “Performance assessment of grid-forming and grid-following converter-interfaced battery energy storage systems on frequency regulation in low-inertia power grids,” _Sustainable Energy, Grids and Networks_ , Sep. 2021\.
  * [11] Y. Zuo, M. Paolone, and F. Sossan, “Effect of voltage source converters with electrochemical storage systems on dynamics of reduced-inertia bulk power grids,” _Electric Power Systems Research_ , vol. 189, p. 106766, Dec. 2020.
  * [12] J. Liu, Y. Miura, and T. Ise, “Comparison of Dynamic Characteristics Between Virtual Synchronous Generator and Droop Control in Inverter-Based Distributed Generators,” _IEEE Transactions on Power Electronics_ , vol. 31, no. 5, pp. 3600–3611, May 2016.
  * [13] T. Qoria, E. Rokrok, A. Bruyere, B. Francois, and X. Guillaud, “A PLL-Free Grid-Forming Control With Decoupled Functionalities for High-Power Transmission System Applications,” _IEEE Access_ , vol. 8, pp. 197 363–197 378, 2020.
  * [14] R. Rosso, S. Engelken, and M. Liserre, “Robust Stability Analysis of Synchronverters Operating in Parallel,” _IEEE Transactions on Power Electronics_ , vol. 34, no. 11, pp. 11 309–11 319, Nov. 2019.
  * [15] J. Rocabert, A. Luna, F. Blaabjerg, and P. Rodríguez, “Control of Power Converters in AC Microgrids,” _IEEE Transactions on Power Electronics_ , vol. 27, no. 11, pp. 4734–4749, Nov. 2012.
  * [16] “IEEE Standard for Synchrophasor Measurements for Power Systems,” IEEE, Tech. Rep., 2011.
  * [17] M. Pignati, M. Popovic, S. Barreto, R. Cherkaoui, G. Dario Flores, J.-Y. Le Boudec, M. Mohiuddin, M. Paolone, P. Romano, S. Sarri, T. Tesfay, D.-C. Tomozei, and L. Zanni, “Real-time state estimation of the EPFL-campus medium-voltage grid by using PMUs,” in _2015 IEEE Power & Energy Society Innovative Smart Grid Technologies Conference (ISGT)_. Washington, DC, USA: IEEE, Feb. 2015, pp. 1–5.
  * [18] “COMMISSION REGULATION (EU) 2016/ 631 - of 14 April 2016 - establishing a network code on requirements for grid connection of generators,” p. 68.
  * [19] T. Weibbach, S. Remppis, and H. Lens, “Impact of Current Market Developments in Europe on Deterministic Grid Frequency Deviations and Frequency Restauration Reserve Demand,” in _2018 15th International Conference on the European Energy Market (EEM)_. Lodz: IEEE, Jun. 2018, pp. 1–6.


We may also combinatorially characterize the unique parameters. For a caterpillar tree, let $ij'$ be a cherry. Assume without loss of generality that from the cherry $ij'$ towards the fixed point, the adjacent leaves are given $i_1, i_2, \dots, i_m, i_{m+1}', \dots$. In particular, we here note that $i_1,\dots, i_m$ are all leaves in the same color as $i$ in the cherry, and $i_{m+1}'$ is the first leaf of the opposite color to $i$ along this traversal (if one exists.)

With the assumption above, the unique parameters of the matroid are exactly  $ij',i_1j', \dots, i_mj'$.

Note that
a unique parameter is an entry of the corresponding symmetric tropical rank 2 matrix that is nonzero and unique up to symmetry.
Therefore, it is enough to consider the case that two leaves $a$ and $b'$ in different colors are in the same branch because $M_{ab}=0$ if the leaves $a$ and $b$ in the same color, are in the same branch. Also, $M_{ab}$ is the minimum of the lengths of the paths from the fixed point to the leaf $a$ and $b$. We now consider when the entry of $M$ is uniquely determined.

The cherry $ij'$ is the only way to achieve $M_{ij}$ which is the length of the path from the fixed point to the leaf $i$ (or $j'$), and so it must be a unique parameter. For $1 \leq l \leq m$, consider the leaf $i_l$ on this branch of the caterpillar tree. The entry $M_{i_lj}$ is the length from $i_l$ to the fixed point. Furthermore, this is the only value that appears in the matrix $M$. Therefore $ij', i_1j', \ldots, i_mj'$ are unique parameters.

It suffices to show that all other cases are not unique parameters. Let $a$, $b'$ be leaves on the same branch with the cherry $ij'$ such that either $b' \neq j'$ or $a$ is not one of the $i, i_2, \ldots, i_m$.
If the leaf $a$ is closer to the fixed point than the leaf $b'$, then $M_{ab} = M_{aj} = M_{ai_{m+1}}$, and so it is not a unique parameter. Now if $b'$ be closer to the fixed point than $a$, then $M_{ab} = M_{ib} = M_{i_2b}$. Either way, $ab'$ cannot be a unique parameter.

We proceed to characterize the circuits of the matroid associated with a caterpillar symbic tree. To do this, we encode each column of the distance parameter matrix $D$ using semisimple graphs, i.e.  graphs allowing a loop on each vertex. In other words, if the $ij'$ column of $D$ corresponds to an edge in a graph connecting vertices $i$ and $j$. As a convention, we assume a loop as a cycle of size $1$. Define a connected two odd cycles as a connected graph that two odd cycles connected by one path whose path length can be possibly zero. For a semisimple graph $G=(V,E)$, we identify its edge set as a ground set of the matroid. For example, if we say $G$ is dependent, then the submatrix of $D$ consists of columns corresponding to edges of $G$ is rank deficient. Furthermore, each set $p_i$ in the partition $\{p_1,\dots, p_{n-1}\}$ may also be considered as a set of edges for semisimple graphs. For a matrix $M$ with $\frac{n(n-1)}{2}$ columns and a semisimple graph $G=(V,E)$ of $n$ vertices, let $M_G$ be the submatrix of $M$ consists of columns corresponding to edges of $G$.

A circuit of the matroid associated with the matrix $L$ is either an even cycle or a connected two odd cycles.

(idea?) Note that the matrix $L$ corresponds to a linear combination of coordinate matrices (<ref>).
It means that it is equivalent to the zero matrix in the sense that any $3\times 3$-submatrix has all equal terms of tropical determinant. Therefore, linear relations obtained by any term of $3\times 3$ tropical determinant of any $3\times 3$-submatrix of zero matrix are circuits of the matroid associated with $L$.

Since $L$ is the submatrix of the distance parameter matrix $D$ that consists of rows of $D$, a circuit of the matroid associated with $D$ is also a circuit of the matroid associated with $L$. Hence, to characterize circuits of a matroid associated with $D$, it is enough to show which circuits of the matroid associated with $L$ are still a circuit of that of $D$.

To study the properties of circuits of $L$, we consider the right kernel for each circuit of $L$. For an even cycle, since each vertex has degree $2$, a vector $\begin{bmatrix}
1 & -1 & \cdots & 1 & -1
\end{bmatrix}^\top$ is the right kernel. For a connected two odd cycles, we may consider the following construction of the kernel: For edges of the first odd cycle, we assign $1$ to one of edges incident to the path connecting two cycles (or the other odd cycle if there is no such path). We assign $-1$ and $1$ alternatively to edges of the first cycle so that $1$ is assigned to both edges incident to the path. If there is the path connecting two cycles, for each edge on this path, we assign $-2$ and $2$ alternatively. For the second odd cycle, we also assign $1$ and $-1$ to each edge alternatively that we have the same numbers on both edges incident to the path (or the first odd cycle). In this case, if the last edge of the path has $-2$ (or $2$ respectively) then, we assign $1$ (or $-1$ respectively) to both edges incident to this edge. If there is no path connecting two odd cycles, then we assign $-1$ on both edges incident to the first odd cycle. Since each coordinate of the right kernel corresponds to an edge of the graph, this assignment on each edge induces a vector which is the kernel of connected two odd cycles.  example needed?

We propose properties of circuits of the matrix $L$ using their right kernel vectors.

Let $G=(V,E)$ be either an even cycle or a connected two odd cycles. Suppose that $G$ satisfies that $|E\cap p_i|=1$ for some $i$. Then, $G$ is independent to the matroid associated with a caterpillar symbic tree.

We consider matrices $D_G$ and $L_G$.
It is enough to show that the right kernel of $L_G$ we discussed is no more a kernel of $D_G$. Note that $G$ is a circuit of the linear matroid associated with $L$. Because $|E\cap p_i|=1$, there is a row of $D_G$ with all zero entries except only one entry of $1$. Then, the right kernel of $L_G$ is not a kernel of $D_G$.

Let $G=(V,E)$ be either an even cycle or a connected two odd cycles. If $G$ is a circuit of the matroid associated with a caterpillar symbic tree and $|E\cap p_i|$ is nonzero for some $i$, then $|E\cap p_i|=2,3$ or $4$.

Note that $|E\cap p_i|$ being nonzero means that there is a vertex of $G$ whose degree is at least $|E\cap p_i|$. Since $G$ is either an even cycle or a connected two odd cycles, $|E\cap p_i|=2$ is possible clearly. If we consider the case that $G$ is a connected two odd cycles, then there is a vertex whose degree is at least $3$ (i.e. the vertex connecting an odd cycle and the path). Hence, $|E\cap p_i|=3$ is possible. Finally, if $G$ is a connected two odd cycle with no path connecting two cycles then the vertex connecting two cycles must have degree $4$. It gives the case of $|E\cap p_i|=4$.

On the other hand $|E\cap p_i|=1$ is impossible due to Lemma <ref>. Also, since there cannot be a vertex of degree at least $5$, we have $|E\cap p_i|<5$.

Let $G=(V,E)$ be either an even cycle or a connected two odd cycles. If $G$ is a circuit of the matroid associated with a caterpillar tree and $|E\cap p_i|>0$ for some $i$, then there is no edge of $G$ from the main block of $p_i$.

Suppose that $G$ is a circuit but there is $p_i$ such that $|E\cap p_i|>0$ and $G$ has an edge from the main block of $p_i$. Then, since $G$ is a connected graph, we may assume that edges of $G$ from $p_i$ are incident to the edge in the main block of $p_i$. Since the vertex incidents to these edge has degree at least $3$, $G$ can't be an even cycle.

Assuming that $G$ is a connected two odd cycles,
we claim that $|E\cap p_i|=2$. If $|E\cap p_i|=4$, then since there is an edge from the main block of $p_i$, there is a vertex of degree at least $5$ which is a contradiction. Suppose that $|E\cap p_i|=3$. Then, due to the existence of the edge from the main block of $p_i$, we have a vertex of degree $4$, and hence $G$ is a connected two odd cycles without the connecting path. Therefore, at least one edge from $E\cap p_i$ must belong to an odd cycle that contains the edge from the main block of $p_i$. It means that there is another set $p_j$ sharing the same main block as $p_i$, and an edge in $E\cap p_j$ that is an edge from the same odd cycle. Since $|E\cap p_j|\geq 2$ due to Lemma <ref>, a vertex incident to edges in $E\cap p_j$ has a degree at least $3$ which is a contradiction because it is a vertex in a cycle.
Therefore, $|E\cap p_i|=2$.

We further claim that edges in $E\cap p_i$ only belong to a single odd cycle. In other words, an edge in the main block of $p_i$ must either be a part of the path or belong to the other cycle. Assuming that this is not true,
consider an odd cycle that contains the edge from the main block of $p_i$ and an edge in $E\cap p_i$ that is contained in the same odd cycle. Then, similarly to the argument that disproves $|E\cap p_i|=3$, there exists another set $p_j$ sharing the same main block as $p_i$. Additionally, there is a vertex incident to edges from $E\cap p_j$ with a degree of at least $3$, which leads to a contradiction.

Now, consider the kernel vector $v$ of the submatrix $L_G$. If we consider the coordinates of the kernel corresponding to edges in $E\cap p_i$, then they must be equal to $1$ since these two edges are incident to the path or the other odd cycle. Considering the submatrix $A_G$, there is a row whose entries are all zeros except for only two entries of $1$ corresponding to $E\cap p_i$. Then, $v$ is not a kernel of this row. Therefore, $G$ is independent.

For a given caterpillar symbic tree, a semisimple graph $G$ is a circuit of as matroid associated with the tree if and only if $G$ is either an even cycle or a connected two odd cycles satisfying that

* edges of $G$ doesn't have edges corresponding to coloops,

* if $|E\cap p_i|>0$ for some $i$, then $|E\cap p_i|>1$, and

* if $|E\cap p_i|>0$ for some $i$, then there is no edge of $G$ corresponding to columns of the main block of $p_i$.

We first prove the `only if' part. Let $G$ be a semisimple graph that is either an even cycle or a connected two odd cycles. The first condition is clear since a circuit cannot contain a coloop. The second condition follows from Lemma <ref>, and the third condition is obtained by Lemma <ref>.

Conversely, we show the `if' part.

For a semisimple graph $G$ that is either an even cycle or a connected two odd cycles, we show that the kernels we discussed for the submatrix $L_G$ is still a kernel of the submatrix $A_G$ when three conditions are satisfied. To show this, it is enough to consider rows of $N$ and columns from the partition $\{p_1,\dots, p_{n-1}\}$ because they are only nonzero entries occur. For an even cycle, since all vertices have degree $2$, we have that if $|E\cap p_i|>0$, then $|E\cap p_i|=2$. Since edges in $E\cap p_i$ are incident to each other, their corresponding coordinates of the kernel have the different signs $1$ and $-1$. Therefore, it is the kernel of the submatrix of $N$ corresponding to the even cycle. For a connected two odd cycles, due to Lemma <ref>, there are three different cases: $|E\cap p_i|=2,3$ and $4$. For a fixed $p_i$, suppose that $|E\cap p_i|=2$. Since there is no edge in the main block of $p_i$, the path connecting two odd cycles (or the other odd cycle if there is no such path) is not both incident to these two edges. It means that these two edges have coordinates of the kernel with different signs. Hence, rows of $N$ also has the same kernel. Assume that $|E\cap p_i|=3$. It means that these three edges are incident to a vertex of degree $3$. It means that two edges should be edges from one odd cycle and the other edge is the edge from the path. Therefore, their corresponding coordinates of the kernel should be $\begin{bmatrix}
1 & 1 & -2
\end{bmatrix}$. Hence, the corresponding row of $N$ has the same kernel. Finally, consider the case that $G$ has $4$ edges from $l_i$. This is only possible when two odd cycles are directly connected without the path. Then, their corresponding coordinates of the vector must be $\begin{bmatrix}
1 & 1 & -1 & -1
\end{bmatrix}$ and its corresponding row of $N$ must have the same kernel.

Since Theorem <ref> presents a characterization of the circuits of the matroid associated with a caterpillar tree, we may obtain independent sets of this matroid. We note that these independent sets correspond to choices of entries that can be independently chosen to complete a given matrix to a rank $2$ matrix. As this way of figuring out independent sets is not deterministic however, we propose a method to softly verify the independence of a given subset $\mathcal{I}$ of the ground set $E$. Note that a set of rank $2$ symmetric matrices can be expressed as a sum of two rank $1$ symmetric matrices. Therefore, a point $(x,y) = (x_1, \dots, x_n, y_1, \dots, y_n) \in \mathbb{C}^{2n}$ can be transformed into a rank $2$ symmetric matrix using the map $(x,y) \mapsto x \cdot x^\top + y \cdot y^\top$. The Jacobian of this map has dimensions $2n \times \frac{n(n-1)}{2}$. For a subset $\mathcal{I}$ of $E$, if a submatrix of the Jacobian, obtained by selecting the columns corresponding to elements in $\mathcal{I}$, has rank $|\mathcal{I}|$ when evaluated at a generic point in $\mathbb{C}^{2n}$, then $\mathcal{I}$ is an independent set. An issue with this method is that it may produce false positives, indicating that a given set is independent when it is not. This is particularly likely to occur when the input is an inexact floating-point number. To mitigate this risk, exact or algebraic numbers (e.g. $\mathbb{Q}$) may be used instead, as long as they are known to be sufficiently generic.

§ TO DO

* Find a nice combinatorial description of tropical rank $2$ (or symmetric Kapranov rank $2$) symmetric matrices, in terms of max/min subdivisions of $2 \Delta_{n-1}$.

* Use this combinatorial description to determine topological properties of those spaces.  Similarly to:
<https://arxiv.org/abs/0711.0944> or

* Use this combinatorial description to determine bases/circuits for the algebraic matroid of the space of tropical rank $2$ (or symmetric Kapranov rank $2$) symmetric matrices.  Kisun has some conjectures (and theorems?).

* Find the number of top-dimensional faces of the complex as per section 4 of