# Modularized Bilinear Koopman Operator for Modeling and Predicting Transients
of Microgrids

Xinyuan Jiang, Yan Li, and Daning Huang X. Jiang and Y. Li are with the
Department of Electrical Engineering, The Pennsylvania State University,
University Park, PA 16802, USA (e-mail: yql5925@psu.edu). D. Huang is with the
Department of Aerospace Engineering, The Pennsylvania State University,
University Park, PA 16802, USA (e-mail: daning@psu.edu).

###### Abstract

Modularized Koopman Bilinear Form (M-KBF) is presented to model and predict
the transient dynamics of microgrids in the presence of disturbances. As a
scalable data-driven approach, M-KBF divides the identification and prediction
of the high-dimensional nonlinear system into the individual study of
subsystems; and thus, alleviating the difficulty of intensively handling high
volume data and overcoming the curse of dimensionality. For each subsystem,
Koopman bilinear form is applied to efficiently identify its model by
developing eigenfunctions via the extended dynamic mode decomposition method
with an eigenvalue-based order truncation. Extensive tests show that M-KBF can
provide accurate transient dynamics prediction for the nonlinear microgrids
and verify the plug-and-play modeling and prediction function, which offers a
potent tool for identifying high-dimensional systems. The modularity feature
of M-KBF enables the provision of fast and precise prediction for the
microgrid operation and control, paving the way towards online applications.

###### Index Terms:

Modularized Koopman bilinear form (M-KBF), data-driven modeling, Koopman
eigenanalysis, extended dynamic mode decomposition (EDMD), transient dynamics
prediction, microgrids, distributed energy resources (DERs).

##  I Introduction

Distributed Energy Resources (DERs), such as photovoltaic (PV) and wind power,
is seen as a great opportunity to achieve the target of modernizing the power
systems. Microgrids have been developed to integrate those DERs. Considering
most DERs are integrated into microgrids through power-electronic interfaces,
the system’s inertia is significantly reduced. Consequently, microgrids are
sensitive to disturbances like PV fluctuations and load disturbances; and
thus, it is of primary importance to investigate the transient dynamics of
microgrids subjected to disturbances.

Efforts on the modeling and control of transient dynamics since the 1980s
mainly focused on the mathematical model-based simulation and analysis [1].
Despite that several control approaches have been developed for stabilizing
the system during transients, it is still elusive that _how to analyze and
predict the system’s transient dynamics when an accurate system model is
unavailable._ The rapidly developing machine learning and artificial
intelligence technologies [2, 3, 4] provide a potent means to resolve this
challenge from the data-driven perspective.

Recognizing the fact that microgrid is a quintessentially nonlinear dynamical
system, there are several existing data-driven approaches to identify the
transient dynamics model of a nonlinear system through its operating data.
They can fall into three major categories: (1) linear models as local
linearization of a globally nonlinear system, (2) nonlinear models to directly
capture the global nonlinearity, and (3) linear models in the embedded space
to reproduce the global nonlinearity.

First, the local linear models are the most commonly used system
identification method and develop a high-order linear system with input and
output to approximate the system dynamics near an equilibrium point. From the
perspective of linear system identification [5], existing methods include
Prony analysis [6, 7] and state space methods, e.g., Minimal Realization
algorithm [8], Eigenvalue Realization Algorithm (ERA) [9], Matrix Pencil
method [10], Hankel Total Least Squares (HTLS) [11], subspace identification
[12], Dynamic Mode Decomposition (DMD) [13], etc. Beyond the linear system
perspective, there are nonparametric spectrum estimation such as Welch
periodgram [14] and parametric methods, including Yule-Walker [15], Frequency
Domain Decomposition [16], etc. These aforementioned approaches generate
accurate model for linear systems given sufficient system responses and are
straightforward to perform with guaranteed convergence. However, the linear
system methods do not have the extrapolation ability due to their nature of
local linearization; and thus, it is not suitable to directly apply them to
identify microgrid systems that are typically nonlinear for the entire
operating envelope.

Second, several data-driven methods have also been developed to identify a
nonlinear system to capture the transient dynamics globally over the entire
state space. Representative methods include Taylor series [17], Volterra
series [18], Sparse Identification (SINDy) [19], etc. These methods can
theoretically identify an accurate model if the correct nonlinear terms are
used. However, their implementation usually result in a non-convex
optimization formulation and need a large computational effort; and the
identified models are limited by the high algorithmic complexity in the
prediction stage, which poses a challenge for control applications.

Third, to reconcile the dilemma between local linearization and global
nonlinearity, the global linearization method, represented by the Koopman
operator theory, has been seen as a promising paradigm for the data-driven
modeling and control of nonlinear systems [20, 21, 22, 23, 24] and has been
widely used in the domain of power systems [25, 26, 27]. The Koopman theory
states that a nonlinear dynamical system can be represented using an infinite-
dimensional linear operator, that is characterized by an infinite set of
eigenvalues and eigenfunctions, or eigenpairs, on the so-called embedded
space. Using a converged, truncated set of eigenpairs, the nonlinear dynamics
can be represented by a finite-dimensional linear system on a Koopman
invariant embedded subspace with the eigenfunctions as new coordinates.
Compared to the traditional nonlinear system identification models, the
Koopman-based model evolves linearly in the embedded space, leading to a
simplified representation of the original nonlinear system that is amenable
for efficient control applications.

In the Koopman theoretic framework, the data-driven modeling centers on the
identification of a set of eigenpairs that is sufficiently comprehensive to
represent the nonlinear dynamics linearly in the embedded subspace and
reconstruct it in the original state space. The eigenfunctions are usually
constructed as functionals of the observable functions of the state variables
in the nonlinear dynamics. The existing methods for identifying Koopman
eigenfunctions can generally be divided into two types. The first type is the
Extended Dynamic Mode Decomposition (EDMD) and its variants, where a
predetermined set of observable functions is given [28]. The observable
functions lift the original state in the training data into an embedded space,
on which a linear system identification method such as DMD can be applied to
identify a linear model as well as a Koopman invariant subspace spanned by the
correspondingly approximated Koopman eigenfunctions. The set of observable
functions are required to approximate a Koopman invariant subspace as much as
possible; otherwise, the identified model may be no better than a local linear
model. Moreover, the EDMD-type methods are usually limited by the curse of
dimensionality, as the number of observable functions may grow exponentially
with the number of states, which makes this type of method infeasible to
directly use for identifying high-dimensional nonlinear dynamical systems like
power systems and microgrids. The second type of Koopman methods directly
solves for Koopman eigenfunction basis to build the linear system [29, 30].
This type of methods exploit the connection of the Koopman operators to state-
space geometry either within a basin of attraction or off-attractor. And
hence, they are limited in the case where the state-space geometry is altered
by inputs.

To predict transient dynamics in the Koopman framework, inputs that represent
the disturbances and cause the transient response need to be included in the
linear evolution of eigenfunctions. Previously, the nonlinear system’s inputs
were added to the Koopman model linearly, which resulted in a linear system
with inputs [26]. However, this assumption applied to the function space is
not consistent with how the inputs influence the dynamics of the original
nonlinear system in the state space, which significantly limits the applicable
range of Koopman model. Recently, Koopman model called Koopman Bilinear Form
(KBF) in the context of control-affine systems were proposed in [31, 32],
where the bilinear input has a precise connection to the original state-space
input.

In this paper, a Modularized Koopman Bilinear Form (M-KBF) is developed as a
scalable data-driven modeling approach for efficiently predicting the
transient dynamics of microgrids integrated with power-electronics-interfaced
DERs. The data-driven modeling of a whole microgrid system is divided into the
separate identification of backbone subsystem and DERs dynamics. KBF is used
to develop a precise data-driven model for each DER, since the typical power-
electronics-interfaced DERs have been determined to be control-affine to the
input disturbance from the current or voltage at the DER integration point.
The identified KBF modules for DERs are then integrated into the backbone
subsystem to assemble the entire data-driven model for the whole microgrid.
The novelties of the presented work are summarized below.

  1. 1.

M-KBF unlocks the potential of KBF for high-dimensional nonlinear systems by
overcoming the curse of dimensionality. It develops the data-driven model of
high-dimensional systems by first breaking it into several manageable
subsystems that each is modeled by a KBF model and then connecting the
subsystems via the accurately modeled inputs.

  2. 2.

M-KBF offers modularity in both data-driven modeling and prediction. It does
not need to send high volume data to a central computational resource; and
thus, issues caused by data acquisition like data privacy and poisoning can be
alleviated or fully avoided. Meanwhile, local distributed computational
resources can be exploited for identifying large-scale dynamical systems.

  3. 3.

The modularity design enables M-KBF with the plug-and-play versatility. The
KBF models for subsystems can be built from a relatively small amount of
training data of only the subsystem under consideration. Moreover, the
identified KBF models for subsystems like DERs can be repeatedly used for
other systems, which is advantageous over previous data-driven methods in the
Koopman framework that focus on modeling the whole system.

  4. 4.

M-KBF ensures accurate predictions because each identified KBF models are
independent to the initial choice of the global DQ reference frame, which is
found to be important for nonlinear system identification. It is realized by
introducing random DQ frames to the measurement data set. This property
signifies a major difference to the linear data-driven models.

Besides, considering the fact that the voltage and current usually respond to
input change with no time delay, a modification to the explicit solution
procedure for the connected data-driven model is also presented to reduce its
intrinsic time delay effect on prediction. The prediction obtained from M-KBF
can potentially be applied to provide guidance for system control, such as
performing model predictive control, which is the authors’ next work.

The reminder of the paper is organized as follows. Section II establishes the
modularized Koopman bilinear form. Section III introduces the determination of
eigenfunctions through the EDMD method with order truncation. Section IV
presents the prediction of transient dynamics based on M-KBF. In Section V,
tests on a microgrid system verify the effectiveness and efficiency of the
presented method in modeling and predicting transients. Conclusions are drawn
in Section VI.

##  II Modularized Koopman Bilinear Form

The essential idea of M-KBF is to decouple a large system into several smaller
subsystems, identify the data-driven model for each subsystem separately
through KBF, and then integrate them into the network model with inputs and
outputs of bus voltages or currents. Considering the transient dynamics of
microgrids is dominated by the nonlinear dynamics of DERs while the grid can
be modeled as a set of linear constraints between node currents and voltages,
the data-driven modeling of the whole system can be separated into identifying
the grid’s algebraic constraint and the dynamics of DER subsystems. The main
advantages of the presented M-KBF are that the modeling and prediction
architecture is scalable and the identified model is linear or control-affine
in the lifted space so linear system identification methods can be leveraged
[27].

###  II-A Identification of Hybrid Network Parameter Matrix

Taking into account the different functions of grid-forming and grid-following
DERs in the microgrid system [33], the network that connects DERs and power
loads can be represented by a hybrid network parameter matrix [34], as given
in (1), which gives the constraint between node currents and voltages.

$\begin{bmatrix}\mathbf{I}_{\mathcal{P}_{1}}\\\
\mathbf{V}_{\mathcal{P}_{2}}\end{bmatrix}=\mathbf{H}\begin{bmatrix}\mathbf{V}_{\mathcal{P}_{1}}\\\
\mathbf{I}_{\mathcal{P}_{2}}\end{bmatrix},$ (1)

where the partition $\mathcal{P}_{1}$ is the set of nodes connecting to the
grid-forming DERs, and the partition $\mathcal{P}_{2}$ is the set of nodes
connecting to the grid-following DERs and/or power loads. Let $\mathcal{N}$ be
the set of microgrid nodes, then
$\mathcal{P}_{1}\cup\mathcal{P}_{2}\subseteq\mathcal{N}$. The goal is to find
a hybrid network parameter matrix $\mathbf{H}$ that satisfies (1).

Specifically, the grid-forming DER inverters take node current as input and
node voltage as output in (1), while the grid-following inverters take node
voltage as input and node current as output. The power loads have the same
input and output definitions as grid-following DERs. Therefore, the network
has complementary input and output ports at the DER inverter nodes, leading to
the hybrid network parameter matrix $\mathbf{H}$.

The hybrid network parameter matrix can be obtained by using two methods,
namely by obtaining the matrix analytically from the admittance matrix if it
is known, or by solving a least-squares problem based on the measurement of
node voltages and currents. Note that both methods need to meet the condition
given in Remark 1 to guarantee that $\mathbf{H}$ exists. To derive
$\mathbf{H}$ from the admittance matrix, there are two steps [34], including
the elimination of nodes with zero current injections from generations or
loads (e.g., Kron reduction [35]) and the switching of a set of node currents
and voltages (for partition $\mathcal{P}_{2}$). To solve for the hybrid
network parameter matrix from measurements of node voltages and currents, it
is typical to formulate (1) into a least-squares problem, which is well posed
given the existence of $\mathbf{H}$ (Remark 1).

######  Remark 1 (Sufficient condition for the existence of $\mathbf{H}$)

In a connected network, assume all the branches are not electromagnetically
coupled and have nonzero admittance. If all lines are passive, i.e., having
strictly positive real part in line admittance, then the hybrid network
parameter exists [34]. The existence follows from the ability to transform the
admittance matrix through the two steps described above. If, in particular, no
shunt elements (shunt capacitors or constant impedance loads) are included in
the network, then the hybrid matrix always exists for a connected network.

###  II-B Bilinear Koopman Surrogate Model for DERs

Koopman Bilinear Form is a global bilinearization of control-affine nonlinear
systems in an operator theoretic view. Each grid-forming controlled DER with
node current inputs or the grid-following DER with node voltage inputs can be
expressed in the following form of control-affine system, with Vf controlled
DER provided as an example in the Appendix.

$\displaystyle\dot{\mathbf{x}}$
$\displaystyle=\mathbf{f}(\mathbf{x})+\sum_{i=1}^{m}\mathbf{g}_{i}(\mathbf{x})u_{i},$
(2a) $\displaystyle\mathbf{y}$ $\displaystyle=\mathbf{C}\mathbf{x},$ (2b)

where $\mathbf{x}\in\mathbb{X}\subseteq\mathbb{R}^{n}$ is the state vector for
one DER, $\mathbf{u}=[u_{1}\ldots u_{m}]^{\mathsf{T}}\in\mathbb{R}^{m}$ is the
input vector, $\mathbf{y}\in\mathbb{R}^{2}$ is the output vector,
$\mathbf{C}\in\mathbb{R}^{2{\times}n}$ is the output matrix,
$\mathbf{f}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{X}\to\mathbb{X}$ and
$\mathbf{g}_{i}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{X}\to\mathbb{X}$
are the flow and control vector fields for the two typical inverter controls.
The KBF for identifying DER dynamics are then introduced from the following
three aspects.

####  II-B1 Observable Function

An observable function $\varphi\mathrel{\mathop{\mathchar
58\relax}}\mathbb{X}\to\mathbb{C}$ is a complex-valued function of the state
vector $\mathbf{x}$. Let $\mathcal{F}$ be the space of all possible observable
functions so that $\varphi\in\mathcal{F}$. Let
$\mathbf{\Phi}_{\mathbf{u}}(t,\mathbf{x}_{0})$ be the flow map of the system
(2) at time $t>0$ starting from an initial condition $\mathbf{x}_{0}$ with
input $\mathbf{u}$. The time-varying observable
$\psi(t,\mathbf{x})\triangleq\left.\varphi(\mathbf{x})\right|_{\mathbf{x}=\mathbf{\Phi}_{\mathbf{u}}(t,\mathbf{x}_{0})}$
of the system (2) is the solution of the following partial differential
equation, which lays the foundation for bilinearization.

$\begin{split}\frac{\partial\psi}{\partial
t}&=L_{\mathbf{f}}\psi+\sum_{i=1}^{m}u_{i}L_{\mathbf{g}_{i}}\psi,\\\
\psi(0,\mathbf{x})&=\varphi(\mathbf{x}),\end{split}$ (3)

where $L_{\mathbf{f}}\triangleq\mathbf{f}\cdot\nabla$,
$L_{\mathbf{g}_{i}}\triangleq\mathbf{g}_{i}\cdot\nabla,\,i=1,\ldots,m$ are the
Lie derivatives [36] with respect to the drift and control vector fields
$\mathbf{f}(\mathbf{x})$ and $\mathbf{g}_{i}(\mathbf{x})$, which are linear
operators on $\mathcal{C}^{1}(\mathbb{X})$.

####  II-B2 Koopman Operator

The Koopman operator is defined for the system (2) with zero inputs,

$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})$ (4)

Assume ${\mathbf{\Phi}}(t,\mathbf{x}_{0})$ is the flow map of (4) for the time
$t>0$, the continuous time Koopman operator is defined as
$\mathcal{K}^{t}\mathrel{\mathop{\mathchar
58\relax}}\mathcal{F}\to\mathcal{F}$ such that,

$(\mathcal{K}^{t}\varphi)(\cdot)=\varphi\circ\mathbf{\Phi}(t,\cdot),$ (5)

where $\circ$ represents function composition. Based on (5), we can see the
Koopman operator is linear, namely,

$\displaystyle\mathcal{K}^{t}(\alpha_{1}\varphi_{1}+\alpha_{2}\varphi_{2})$
$\displaystyle=\alpha_{1}\varphi_{1}\circ\mathbf{\Phi}(t,\cdot)+\alpha_{2}\varphi_{2}\circ\mathbf{\Phi}(t,\cdot)$
(6)
$\displaystyle=\alpha_{1}\mathcal{K}^{t}\varphi_{1}+\alpha_{2}\mathcal{K}^{t}\varphi_{2}.$

Therefore, the Koopman operator can be characterized by its eigenvalues and
eigenfunctions,

$\mathcal{K}^{t}\phi=e^{t\lambda}\phi,$ (7)

where $\lambda\in\mathbb{C}$ is the Koopman eigenvalue and
$\phi\in\mathcal{F}$ is the corresponding eigenfunction.

Through the definition of Lie derivatives, the infinitesimal generator of the
Koopman operator is equal to the Lie derivative of the drift vector field,
i.e., $L_{\mathbf{f}}=\lim_{t\to 0}\frac{\mathcal{K}^{t}-I}{t}$, where $I$ is
the identity operator, so $L_{\mathbf{f}}$ is also referred to as the Koopman
generator. The eigenvalue and eigenfunction relation in (7) can also be
expressed in terms of the Lie derivative,

$L_{\mathbf{f}}\phi=\lambda\phi.$ (8)

From (8) and the definition of $L_{\mathbf{f}}$, for any two eigenpairs
$(\lambda_{1},\phi_{1})$ and $(\lambda_{2},\phi_{2})$, we have

$L_{\mathbf{f}}(\phi_{1}\cdot\phi_{2})=(\lambda_{1}+\lambda_{2})(\phi_{1}\cdot\phi_{2}).$
(9)

Thus, there exist infinitely many eigenfuncions and eigenvalues for the
Koopman operator.

####  II-B3 Bilinearzation

The goal of bilinearzation of (2) is to choose a set of observable functions,

$T(\mathbf{x})=[\varphi_{1}(\mathbf{x}),\varphi_{2}(\mathbf{x}),\,\ldots,\,\varphi_{N}(\mathbf{x})]^{\mathsf{T}},$
(10)

such that their evolution over time is that of a bilinear system (11), which
mirrors (3) with the observable function basis $T(\mathbf{x})$.

$\dot{\mathbf{z}}=\mathbf{A}\mathbf{z}+\sum_{i=1}^{m}u_{i}\mathbf{B}_{i}\mathbf{z},\quad\mathbf{z}(0)=T(\mathbf{x}_{0}),$
(11)

where $\mathbf{A}$ and $\mathbf{B}_{i}$ will be determined in Section III. The
condition for bilinearizability into (11) with the state embeddings
$T(\mathbf{x})$ being the eigenfunctions of $L_{\mathbf{f}}$ is provided in
Theorem 1 from [31]. And furthermore, bilinearizability condition with finite
eigenfunction embeddings [31] is presented in the following Theorem for
completeness.

###### Theorem 1

If a set of Koopman eigenfunctions
$\\{\phi_{1},\phi_{2},\ldots,\phi_{n}\\},\,n\in\mathbb{N}$ of the unactuated
system forms an invariant subspace of $L_{\mathbf{g}_{i}},\,i=1,\ldots,m$,
then the system (2) is bilinearizable with an $n$ dimensional state space.

Given that the condition of Theorem 1 is satisfied, the control-affine DER
subsystem (2) is bilinearizable with a finite number of eigenfunction
embeddings. Next, we will obtain a finite-dimensional bilinear model (11) by
modifying the Extended Dynamic Mode Decomposition (EDMD) algorithm to include
blinear inputs.

##  III Data-Driven Identification of KBF

The data-driven identification of the KBF equivalent model (11) of DER
subsystems is carried out through a least-squares formulation using a
predetermined set of observable functions. The eigenfunction embeddings and
the KBF system are identified simultaneously from the state and input
measurements of the DER subsystems.

The general rule of selecting a dictionary of observable functions is that
these functions should span a rich subspace so that a certain set of
eigenfunctions can be approximated by their projections onto this subspace.
There are several choices of observable functions [37], including monomials,
radial basis functions, Hermite polynomials, and Chebyshev polynomials. And in
this work we choose monomial functions as the observable dictionary. In this
paper, we use EDMD with additional bilinear inputs to identify the KBF model
(11). As will be shown in the subsection A, the KBF model identified in this
way inherently satisfies the bilinearizability condition in Theorem 1.

###  III-A Approximated Eigenfunctions

Given the conditions in Theorem 1, the original system (2) can be identified
by the surrogate KBF model (11) in the coordinates of a finite set of
eigenfunctions of $L_{\mathbf{f}}$. One option to approximate such
eigenfunctions of $L_{\mathbf{f}}$ is by applying EDMD on the sampled
trajectories of the system with zero inputs. Since the inputs to DER
subsystems include node current or voltage that cannot be held constant during
system transients, it is not feasible to independently identify the
eigenfunctions using EDMD. Instead, the invariant eigenfunction embeddings and
the bilinear system can be identified together using a predetermined
dictionary of observable functions.

The requirement on the observable functions is that their span approximates an
invariant eigenspace of $L_{\mathbf{f}}$ w.r.t (3) that is defined in Theorem
1. Then given the KBF system (11) and the observable functions
$T(\mathbf{x})$, for each left eigenvector $\mathbf{w}_{i}$ and eigenvalue
$\lambda_{i}$ of $\mathbf{A}$,
$\hat{\phi}(\mathbf{x})=\mathbf{w}_{i}^{*}T(\mathbf{x})$ is proved to be an
eigenfunction of $L_{\mathbf{f}}$ with eigenvalue $\lambda_{i}$ [38]. Assume
that all eigenvalues of $\mathbf{A}$ are distinct. Denote the left eigenvector
matrix $\mathbf{W}$ and eigenvalue matrix $\mathbf{D}$ of $\mathbf{A}$ such
that $\mathbf{W}^{*}\mathbf{A}=\mathbf{D}\mathbf{W}^{*}$. Then the
approximated eigenfunction embeddings
$\tilde{\mathbf{z}}=\mathbf{W}^{*}\mathbf{z}=\mathbf{W}^{*}T(\mathbf{x})$ is a
linear transformation of the predetermined observable function embeddings,
which gives us the following KBF system with the eigenfunction embedded state,

$\dot{\tilde{\mathbf{z}}}=\mathbf{D}\tilde{\mathbf{z}}+\sum_{i=1}^{m}u_{i}\tilde{\mathbf{B}}_{i}\tilde{\mathbf{z}},\quad\tilde{\mathbf{z}}(0)=\mathbf{W}^{*}T(\mathbf{x}_{0}),$
(12)

where
$\tilde{\mathbf{B}}_{i}=\mathbf{W}^{*}\mathbf{B}_{i}(\mathbf{W}^{*})^{-1}$.

###  III-B Discretization of KBF for Identification

Since the measurement data are sampled data, the continuous-time KBF system
(11) needs to be discretized to provide the discrete-time template for
identification. Again, we assume that the span of the observable functions
$T(\mathbf{x})$ is a subspace of invariant eigenspace of $L_{\mathbf{f}}$.
Assume the sampling period is $\Delta t$ and that the inputs $u_{i}(t)$ are
constant within each time step (zero-order hold). The KBF system is linear
time-invariant with the solution (13).

$\mathbf{z}(t+\Delta
t)=\exp\Bigg{[}\Bigg{(}\mathbf{A}+\sum_{i=1}^{m}u_{i}(t+\Delta
t)\mathbf{B}_{i}\Bigg{)}\Delta t\Bigg{]}\mathbf{z}(t)$ (13)

Then, by expanding (13) into Taylor series and taking the first order
approximation on $\Delta t$, we can get the discretized KBF model in explicit
(14a) and implicit (14b) forms:

$\displaystyle\mathbf{z}_{k+1}$
$\displaystyle=\mathbf{A}^{d}\,\mathbf{z}_{k}+\sum_{i=1}^{m}u_{i,k+1}\mathbf{B}_{i}^{d}\,\mathbf{z}_{k}+\mathcal{O}(\Delta
t^{2}),$ (14a) $\displaystyle\mathbf{z}_{k+1}$
$\displaystyle=\mathbf{A}^{d}\,\mathbf{z}_{k}+\sum_{i=1}^{m}u_{i,k+1}\mathbf{B}_{i}^{d}\,\mathbf{z}_{k+1}+\mathcal{O}(\Delta
t^{2}),$ (14b)

where $\mathbf{A}^{d}=\exp(\mathbf{A}\Delta
t),\,\mathbf{B}_{i}^{d}=\mathbf{B}_{i}\Delta t$.

###  III-C Least-Squares Formulation

Suppose that a suitable dictionary of observable functions
$T(\mathbf{x})\in\mathbb{R}^{q}$ is selected, we can identify the discretized
KBF models (14) using least-squares formulation similar to EDMD [38], which
gives us the approximated eigenfunction coordinates and the KBF (11). In the
following, we assume that the explicit model (14a) is used, while the implicit
model (14b) can be obtained in a similar way.

Assume that the measurement data is collected from a single sampled trajectory
which includes many transient responses of the DER under investigation. Denote
the state and input vector data as $\mathbf{x}_{k},\,\mathbf{u}_{k}$ for
$k=1,\,\ldots\,,\,N$. We will organize the measurement data into the following
matrices:

$\displaystyle\mathbf{X}_{1}$
$\displaystyle=\Big{[}T(\mathbf{x}_{1}),\,T(\mathbf{x}_{2}),\,\ldots\,,\,T(\mathbf{x}_{N-1})\Big{]},$
(15a) $\displaystyle\mathbf{X}_{2}$
$\displaystyle=\Big{[}T(\mathbf{x}_{2}),\,T(\mathbf{x}_{3}),\,\ldots\,,\,T(\mathbf{x}_{N})\Big{]},$
(15b) $\displaystyle\mathbf{\Gamma}_{i}$
$\displaystyle=\Big{[}u_{i,2}T(\mathbf{x}_{1}),\,u_{i,3}T(\mathbf{x}_{2}),\,\ldots\,,\,u_{i,N}T(\mathbf{x}_{N-1})\Big{]}.$
(15c)

Then, a least-squares problem based on the discrete-time model (14a) can be
formulated in (16).

$\displaystyle\;\min_{\mathbf{G}}\left\|\mathbf{X}_{2}-\mathbf{G}\mathbf{\Omega}\right\|_{F}^{2},\,\text{where}$
(16a) $\displaystyle\mathbf{G}$
$\displaystyle=\Big{[}\mathbf{A}^{d},\,\mathbf{B}_{1}^{d},\,\ldots\,,\,\mathbf{B}_{m}^{d}\Big{]},$
(16b) $\displaystyle\mathbf{\Omega}$
$\displaystyle=\Big{[}\mathbf{X}_{1}^{\mathsf{T}},\,\mathbf{\Gamma}_{1}^{\mathsf{T}},\,\ldots\,,\,\mathbf{\Gamma}_{m}^{\mathsf{T}}\Big{]}^{\mathsf{T}}.$
(16c)

The solution to (16) is given by
$\hat{\mathbf{G}}=\mathbf{X}_{2}\mathbf{\Omega}^{+}$. The pseudoinverse
$\mathbf{\Omega}^{+}$ can be calculated from the singular value decomposition
(SVD) of $\mathbf{\Omega}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{*}$, where
$\mathbf{\Sigma}$ is a diagonal matrix whose diagonal elements are the nonzero
singular values in descending order, and $\mathbf{U}$ and $\mathbf{V}$ contain
the corresponding left and right singular vectors. The pseudoinverse is

$\mathbf{\Omega}^{+}=\mathbf{V}\mathbf{\Sigma}^{-1}\mathbf{U}^{*}.$ (17)

###  III-D Singular Value Truncation

Since the assumed discretization of KBF is only first-order accurate w.r.t.
the time step $\Delta t$ and the selected observable functions are usually not
guaranteed to represent the eigenspace, using the exact pseudoinverse in the
solution of (16) does not necessarily lead to optimal predictive accuracy.
Therefore the solution of (16) is modified through the technique of singular
value truncation.

Assume that the number of singular values of $\mathbf{\Omega}$ is reduced to
$r\leq q$. We get the truncated singular value matrix
$\tilde{\mathbf{\Sigma}}\in\mathbb{R}^{r{\times}r}$ and singular vector
matrices
$\tilde{\mathbf{U}}\in\mathbb{R}^{mq{\times}r},\,\tilde{\mathbf{V}}\in\mathbb{R}^{N{\times}r}$.
The regularized solution with singular value truncation is,

$\hat{\mathbf{G}}=\mathbf{X}_{2}\tilde{\mathbf{\Omega}}^{+}=\mathbf{X}_{2}\tilde{\mathbf{V}}\tilde{\mathbf{\Sigma}}^{-1}\tilde{\mathbf{U}}^{*}.$
(18)

From the composition of the parameter matrix $\mathbf{G}$, the singular vector
matrix
$\tilde{\mathbf{U}}=\big{[}\tilde{\mathbf{U}}_{0}^{\mathsf{T}},\,\tilde{\mathbf{U}}_{1}^{\mathsf{T}},\,\ldots\,,\,\tilde{\mathbf{U}}_{m}^{\mathsf{T}}\big{]}^{\mathsf{T}}$
is split such that $\tilde{\mathbf{U}}_{i}\in\mathbb{R}^{q{\times}r}$ to get
the discretized KBF system matrices,

$\hat{\mathbf{A}}^{d}=\mathbf{X}_{2}\tilde{\mathbf{V}}\tilde{\mathbf{\Sigma}}^{-1}\tilde{\mathbf{U}}_{0}^{*},\quad\hat{\mathbf{B}}_{i}^{d}=\mathbf{X}_{2}\tilde{\mathbf{V}}\tilde{\mathbf{\Sigma}}^{-1}\tilde{\mathbf{U}}_{i}^{*}.$
(19)

Assume that $\mathbf{P}=\mathbf{X}_{2}\tilde{\mathbf{V}}$ has full column rank
$r$, so that the data matrices are linearly consistent, i.e.,
$\mathrm{Ker}(\mathbf{X}_{2})\subseteq\mathrm{Ker}(\tilde{\mathbf{\Omega}})$
[39], then we can reduce the dimension of the identified discrete KBF model by
taking the state transform $\bar{\mathbf{z}}=\mathbf{P}^{+}\mathbf{z}$, which
gives us the reduced KBF model,

$\displaystyle\bar{\mathbf{z}}_{k+1}$
$\displaystyle=\bar{\mathbf{A}}^{d}\bar{\mathbf{z}}_{k}+\sum_{i=1}^{m}u_{i,k+1}\bar{\mathbf{B}}_{i}^{d}\bar{\mathbf{z}}_{k},\,\text{where}$
(20a) $\displaystyle\bar{\mathbf{A}}^{d}$
$\displaystyle=\tilde{\mathbf{\Sigma}}^{-1}\tilde{\mathbf{U}}_{0}^{*}\mathbf{X}_{2}\tilde{\mathbf{V}},$
(20b) $\displaystyle\bar{\mathbf{B}}_{i}^{d}$
$\displaystyle=\tilde{\mathbf{\Sigma}}^{-1}\tilde{\mathbf{U}}_{i}^{*}\mathbf{X}_{2}\tilde{\mathbf{V}}.$
(20c)

It can be checked that the reduced system (20) preserves the dynamics of the
unreduced system (19) by observing that $\bar{\mathbf{A}}^{d}$ and
$\hat{\mathbf{A}}^{d}$ share the same nonzero eigenvalues (see [40]) and the
same is true between $\bar{\mathbf{B}}_{i}^{d}$ and
$\hat{\mathbf{B}}_{i}^{d}$.

In summary, the benefits for using singular value truncation include the
following. (i) Reduce KBF model order in (20); (ii) Form linear consistency
between $\mathbf{X}_{2}$ and $\tilde{\mathbf{\Omega}}$ where the reduced
$\tilde{\mathbf{\Omega}}$ is almost the same as $\mathbf{\Omega}$; and (iii)
SVD truncation in effect adds $l_{2}$ regularization to the elements of
$\mathbf{G}$ [41], which in discrete-time makes the solution more stable.

###  III-E Original State Reconstruction

Since the monomial functions used as observable functions $T(\mathbf{x})$
contains the original state variables $\mathbf{x}$, the original state can be
reconstructed linearly by extracting out their entries using a constant matrix
$\mathbf{C}^{\mathbf{x}}$ as in
$\mathbf{x}=\mathbf{C}^{\mathbf{x}}\mathbf{z}$.

For the reduced model (20) with the reduced observable function embedding, the
hypothetical linear reconstruction takes the form
$\mathbf{x}_{recon}=\mathbf{C}^{\mathbf{x}}\mathbf{P}\tilde{\mathbf{z}}=\mathbf{C}^{\mathbf{x}}\mathbf{P}\mathbf{P}^{+}\mathbf{z}$.
However, the loss of information due to the dimension reduction of
$\tilde{\mathbf{z}}$ may prevent the complete linear reconstruction of the
original states. The information loss turns out to be the case for the Vf
model in our example but not for the PQ model, and therefore the system (19)
is used for the Vf model in order to preserve linear reconstruction.

##  IV Prediction of the Whole Dynamical System

Based on the data-driven modules for the network, the individual models of
grid-forming DERs, grid-following DERs, and power loads can be connected to
form a data-driven model for the whole microgrid system, for the prediction of
its transients. The prediction loop that connects each subsystems is provided
in Algorithm 1, where the Vf and PQ controls are used for grid-forming and
grid-following DER controls, respectively, and constant power loads are used
as an example. Nonetheless, it can be easily modified to ZIP loads or include
more comprehensive load dynamics.

Note that in order to identify an accurate KBF from data, we find that it is
important to use the explicit form (14a) for the PQ controlled DER subsystem
and implicit form (14b) for the Vf controlled DER subsystem. The main
difference between these two DERs lies in that the Vf control has node voltage
as an output while the PQ control has node voltage as an input. In particular,
the dynamics of both types of DERs depend on their PLL state to establish the
controller reference frame, which always depends on the node voltage in the
same time step. Since the node voltage is an internal state in the Vf model
that needs to be solved for, it follows that an implicit form need to be used.
For the PQ model, the explicit form is due to the more reactive nature of the
grid-following control where the node voltage is a direct input to the model.

To simplify notation, we partition all the nodes in the system into Vf node
$\mathcal{P}_{1}$, PQ nodes with no local load $\mathcal{P}_{2a}$, PQ nodes
with local load $\mathcal{P}_{2b}$, and load-only nodes $\mathcal{P}_{2c}$,
after Kron reduction. Assume that the following information is given for
prediction: the initial state of the Vf inverter subsystem
$\mathbf{x}^{p}_{0},\,p\in\mathcal{P}_{1}$, initial states of the PQ inverter
subsystems $\mathbf{x}^{p}_{0},p\in\mathcal{P}_{2a}\cup\mathcal{P}_{2b}$, the
initial mixed node current and voltage vector $\mathbf{u}_{0}$ (input to the
map given by the hybrid matrix $\mathbf{y}=\mathbf{H}\mathbf{u}$ as in (1)),
power loads $S_{k}$, and PQ inverter references $R_{k}$ for
$k=1,\,\ldots\,,\,PredLength$.

The function shiftPhase transforms the input and output variables of each DER
model between the global and local reference frames. This point is further
explained in Section V-D.

1

Input : KBF models $\pi_{1}$ (implicit) for Vf inverter and $\pi_{2}$
(explicit) for PQ inverter, and hybrid network parameter matrix $\mathbf{H}$
for the grid

Output : Node voltage and current in $\mathbf{y}_{k}$ and $\mathbf{u}_{k}$

2
$\mathbf{x}^{p}_{0}=\textnormal{{shiftPhase}}(\mathbf{x}^{p}_{0},-x_{2,0}^{p})$;

3
$\mathbf{z}^{p}_{0}=\mathbf{P}_{\pi_{1}}^{+}T_{\pi_{1}}(\mathbf{x}^{p}_{0}),\,p\in\mathcal{P}_{1}$;

4
$\mathbf{z}^{p}_{0}=\mathbf{P}_{\pi_{2}}^{+}T_{\pi_{2}}(\mathbf{x}^{p}_{0}),\,p\in\mathcal{P}_{2a}\cup\mathcal{P}_{2b}$;

5 for  _$k\leftarrow 0$ to $PredLength$_ do

6 $\mathbf{y}_{k}=\mathbf{H}\,\mathbf{u}_{k}$;

7 $\mathbf{u}_{k+1}=\mathbf{u}_{k}$;

8 for  _$p\in\mathcal{P}_{2b}\cup\mathcal{P}_{2c}$_ do // Load change

9
$\mathbf{u}^{p}_{k+1}\mathrel{+}=\overline{(S^{p}_{k+1}-S^{p}_{k})/\mathbf{y}^{p}_{k}}$;

10

11 end for

12 $\mathbf{y}_{k+1}=\mathbf{H}\,\mathbf{u}_{k+1}$;

13 for  _$p\in\mathcal{P}_{2b}\cup\mathcal{P}_{2c}$_ do // Redo load current

14 $\mathbf{u}^{p}_{k+1}=\overline{S^{p}_{k+1}/\mathbf{y}^{p}_{k+1}}$;

15

16 end for

17 for  _$p\in\mathcal{P}_{2a}\cup\mathcal{P}_{2b}$_ do // PQ inverter

18
$\tilde{\mathbf{y}}^{p}_{k+1}=\textnormal{{shiftPhase}}(\mathbf{y}^{p}_{k+1},-x_{2,0}^{p})$;

19
$\mathbf{z}^{p}_{k+1}=\textnormal{{PQModel}}(\mathbf{z}^{p}_{k},\tilde{\mathbf{y}}^{p}_{k+1},R^{p}_{k+1})$;

20
$\mathbf{x}^{p}_{k+1}=\mathbf{C}^{\mathbf{x}}_{\pi_{1}}\mathbf{P}_{\pi_{1}}\mathbf{z}^{p}_{k+1}$;

21 if  _$p\in\mathcal{P}_{2b}$_ then

22
$\mathbf{u}^{p}_{k+1}\mathrel{+}=\textnormal{{shiftPhase}}(\mathbf{C}\mathbf{x}^{p}_{k+1},x_{2,0}^{p})$;

23

24 else

25
$\mathbf{u}^{p}_{k+1}=\textnormal{{shiftPhase}}(\mathbf{C}\mathbf{x}^{p}_{k+1},x_{2,0}^{p})$;

26

27 end if

28

29 end for

30 $\mathbf{y}_{k+1}=\mathbf{H}\,\mathbf{u}_{k+1}$;

31 for  _$p\in\mathcal{P}_{1}$_ do // Vf inverter

32
$\tilde{\mathbf{y}}^{p}_{k+1}=\textnormal{{shiftPhase}}(\mathbf{y}^{p}_{k+1},-x_{2,0}^{p})$;

33
$\mathbf{z}^{p}_{k+1}=\textnormal{{VfModel}}(\mathbf{z}^{p}_{k},\tilde{\mathbf{y}}^{p}_{k+1})$;

34
$\mathbf{x}^{p}_{k+1}=\mathbf{C}^{\mathbf{x}}_{\pi_{2}}\mathbf{P}_{\pi_{2}}\mathbf{z}^{p}_{k+1}$;

35
$\mathbf{u}^{p}_{k+1}=\textnormal{{shiftPhase}}(\mathbf{C}\mathbf{x}^{p}_{k+1},x_{2,0}^{p})$;

36

37 end for

38

39 end for

Algorithm 1 M-KBF Enabled Prediction Algorithm

##  V Numerical Examples

Figure 1: A typical microgrid test system

A typical microgrid system shown in Fig. 1 is used to test and verify the
effectiveness of the M-KBF method for modeling and predicting the nonlinear
transient dynamics of the microgrid. The test system includes 35 buses, one
Vf-controlled DER, ten PQ-controlled DERs, and eight constant power loads. For
the details of the system, the reader is referred to [42]. For one thing,
synthetic data is prepared by modeling the test system via a set of
Differential Algebraic Equations (DAEs) and simulating it through the
Numerical Differenciation Formula (NDF) method with a fixed step size of
$10^{-3}$ second. For another, simulations of the test system’s transient
responses to disturbances are performed to provide training data for M-KBF to
identify the data-driven model and then predict the system’s transient
dynamics. Each transient response is induced by an approximated step change in
constant power loads and/or power references of PQ-controlled DERs. In
practice, the training data set is obtained from the measurement of the
controller and the RLC filter through the advanced metering infrastructure.

Note that in the test system, the dynamical model for all PQ-controlled DERs
are identical, only with different power outputs. It is designed to
demonstrate that the identified KBF model can be repeatedly used to
efficiently identify the whole system, so as to show that the modularity of
M-KBF leads to its plug-and-play versatility. The M-KBF results are analyzed
from the following four aspects.

###  V-A Data-Driven Modeling for Individual DER via KBF

The synthetic training data for the Vf-controlled DER are the transient
responses of the DER under random changes of power loads and other DERs’
outputs, where the changes in this paper is designed up to $\pm 80\%$ of their
nominal values. The training data for the PQ-controlled DERs are the synthetic
combined transient responses of the DERs 13, 16, 18, and 20, under random
changes of their power outputs that are up to $\pm 20\%$ of their nominal
values. While both the $2^{\text{nd}}$-order and $3^{\text{rd}}$-order
monomials can accurately predict the DER’s transient responses up to $\pm
80\%$ changes of power loads (or DER outputs), Fig. 2 and Fig. 3 show the
individual model predictions in response to up to $\pm 100\%$ input changes.

Figure 2: Prediction of Vf-controlled DER under disturbances in the $100\%$
range Figure 3: Prediction of PQ-controlled DER on node 31 under disturbances
in the $100\%$ range

From the predictions, we can see that observables play an essential role in
identifying the data-driven model and predicting transient dynamics, in
particular when the system is under large disturbances. Detailed comparisons
are:

  * •

The prediction through the KBF model obtained by the $3^{\text{rd}}$-order
monomials remains accurate when the system is under $\pm 100\%$ changes (note
the training data is up to $\pm 80\%$ changes), for both Vf- and PQ-controlled
DERs.

  * •

The prediction based on the $2^{\text{nd}}$-order monomials experiences
equilibrium point drifts for Vf-controlled DER model, which causes the
deviation from the true data.

  * •

The $1^{\text{st}}$-order KBF provides false prediction for the Vf-controlled
DER as shown in Fig. 2; and the prediction trajectory quickly diverges for the
PQ-controlled DERs, which is too large to show in Fig. 3.

###  V-B SVD Order Truncation and Eigenfunctions

In order to obtain an appropriate KBF model for prediction and overcome the
discretization error when arriving at the template model (14), it is necessary
to apply singular value truncation to the data matrix. Here we show that an
optimal truncation order can be inferred by the eigenavlue distribution. Fig.
4 shows the distribution of the approximated Koopman eigenvalues on the Vf
model with monomial observables of up to $3^{\text{rd}}$ order before and
after the optimal truncation order of $155$ is applied. We can see that:

Figure 4: Identified Koopman eigenvalues before (a) and after (b) SVD
truncation (Blue shows the SVD truncation order)

  * •

Before the SVD order truncation is applied to the KBF model solution, there
exists an eigenvalue pattern repetition of the first $155$ eigenvalues, as
illustrated by the black dashed lines in Fig. 4 (a). This is due to the fact
that the sum of two Koopman eigenvalues is another Koopman eigenvalue with the
associated eigenfunction being the product of the two existing eigenfunctions,
as given in (9). As a consequence, the repeated eigenvalues to the left
correspond to higher-order eigenfunctions that are less likely to be
approximated accurately by the monomials of lower orders, which may result in
an unstable KBF model for prediction.

  * •

After the SVD truncation is applied, the identified Koopman eigenvalues
contain no repetitive pattern and the associated eigenfunctions can be well
approximated by the selected low-order monomials so that the resulting KBF
model is not only stable but can accurately represent the dynamics of the DER
subsystem, as has been demonstrated in Fig. 2 and Fig. 3.

We have tested Vf- and PQ-controlled DER models with both the
$2^{\text{nd}}$-order and $3^{\text{rd}}$-order monomials and find that in
each setup there is a narrow range of SVD truncation order that results in the
KBF model having the lowest prediction error. The optimal truncation order is
inferred by the distribution of the identified Koopman eigenvalues obtained
without truncation as in Fig. 4 (a) by retaining only the number of
eigenvalues with no repetitive pattern.

The coefficient matrices of the identified eigenfunction embeddings with the
optimal SVD truncation orders for the Vf and PQ models are shown in Fig. 5 and
Fig. 6, where the complex coefficients are projected onto the real domain by
Koopman canonical transform [43]. The vertical axes in both figures correspond
to the eigenfunctions in the descending order of the real parts of their
eigenvalues. The yellow rectangle in each figure shows that these
eigenfunctions have large coefficients corresponding to the voltage and
current variables in the DER subsystem, which agrees with the fact that these
variables have fast dynamics in the solution of the DAEs.

Figure 5: Eigenfunction coefficients for Vf model Figure 6: Eigenfunction
coefficients for PQ model

###  V-C Prediction based on M-KBF and Error Analysis

Based on Algorithm 1, we can integrate the identified KBF modules to perform
the entire system’s transient prediction. Fig. 7 and Fig. 8 show the voltage
and current predictions of the connected M-KBF model under $\pm 100\%$ range
disturbances. From the predictions, we can see that:

  * •

The Vf model trained on $\pm 80\%$ range data is capable of predicting $\pm
100\%$ range disturbances in the network. The PQ model trained on $\pm 20\%$
range data from DERs 13, 16, 18, and 20, is capable of predicting every PQ-
controlled DER in $\pm 100\%$ range disturbance situations, as shown in Fig. 7
and Fig. 8.

  * •

The predictions demonstrate the M-KBF’s plug-and-play function, since the
identified KBF model for PQ-controlled DERs is repeatedly used in the system.

  * •

The prediction error for each DER model is smaller in the connected system
than that in the individual test because each DER model’s input in the
connected system is a dynamic feedback from the rest of the system.

The average prediction error of the M-KBF system with $2^{\text{nd}}$-order
and $3^{\text{rd}}$-order monomials when the sytem is under different ranges
of disturbances are provided in Fig. 9. Overall, it shows the error increases
as the disturbance increases. When the Vf model adopts the
$3^{\text{rd}}$-order implicit prediction, we have the smallest errors
compared to the $2^{\text{nd}}$-order implicit and $3^{\text{rd}}$-order
explicit predictions. However, it needs a high computational effort, as
summarized in Table I, which is tested on a 2.9GHz PC. Note that Table I also
shows that M-KBF enables faster than real-time predictions, which provides
space for performing predictive control to the dynamical system.

Table I: Average computational time for $20s$ prediction of the whole system with different Vf models $2^{\text{nd}}$-order(explicit) |  $2^{\text{nd}}$-order(implicit) |  $3^{\text{rd}}$-order(explicit) |  $3^{\text{rd}}$-order(implicit)
---|---|---|---
$5.0064$s |  $6.2007$s |  $7.6465$s |  $20.6726$s
Figure 7: Voltage predictions for the whole system with $\pm 100\%$ range
disturbances Figure 8: Current predictions for the whole system with $\pm
100\%$ range disturbances

###  V-D Independence of KBF to the Global Reference Frame

The global DQ reference frame rotating at the nominal frequency allows one to
represent the instantenous voltages and currents as phasor quantities.
However, the initial phase of the rotating DQ frame is not unique, so that if
all phasors in the initial condition rotate by a same angle corresponding to a
shift in the DQ reference frame, then the system dynamics would be the same in
the new reference frame. This property is seen as the zero eigenvalue in the
linearized system but it is a global nonlinear property, which cannot be
recreated by a linear system model. However, the nonlinear KBF model may still
be dependent on a specific DQ reference frame if only one frame is present in
the training data. To illustrate this issue, Fig. 10 shows the prediction of
the Vf KBF model when the reference frame for the initial condition and the
inputs are shifted by only $0.001$ radian. We can see that the predicted
trajectory deviates from the initial equilibrium point and runs parallel to
the simulated trajectory.

To address this issue, the training data is preprocessed as follows. For the
Vf model, we first shift the transient responses’ reference frame to zero and
then randomly add an angle within $\pm 0.01$ radian. For the PQ model, the
reference frames are only shifted to have zero initial PLL phase angle because
the four PQ DERs are already in different reference frames. Therefore, the
nonlinearity related to the changes between the internal PLL reference frame
and the global one is represented in the KBF model in a limited range of PLL
phases. The prediction remains accurate when the PLL phase is within a $\pm
0.04$ radian inverval for the Vf model, and $\pm 0.05$ radian interval for the
PQ model. These limits could be removed by selecting proper observable
functions that are invariant to different reference frames, which is the
authors’ future work.

Figure 9: Average node voltage prediction errors Figure 10: Counter example of
dependence to DQ reference frame

##  VI Conclusions

The paper contributes a scalable data-driven method, M-KBF, to efficiently
model and predict the transient dynamics of nonlinear microgrids under
disturbances. The data-driven modules for subsystems are developed through
EDMD with eigenvalue-based order truncation and then integrated into a
combined model for microgrid systems. Test results have demonstrated the
effectiveness of M-KBF in providing fast and precise transient predictions.
For future work, M-KBF will be further developed to handle data with noise, to
identify the system by using output-only measurements like bus voltage and
current, and to perform predictive control based on the predictions obtained
from M-KBF.

The DER connection circuit and the Vf-controller (given as an example) are
shown in Fig. 11 and Fig. 12, respectively. Ten differential equations can
then be developed, which are not given here due to the page limit. These
equations form a control affine system, where the state vector is
$\mathbf{x}=\big{[}x_{1},\,x_{2},\,x_{3},\,x_{4},\,x_{5},\,x_{6},\,V_{C}^{D},\,V_{C}^{Q},\,I_{L}^{D},\,I_{L}^{Q}\big{]}^{\mathsf{T}}$,
the input vector is the node current
$\mathbf{u}=\big{[}I_{N}^{D},\,I_{N}^{Q}\big{]}^{\mathsf{T}}$, and the output
to the network is the node voltage that is equal to the state variables
$\mathbf{y}=\big{[}V_{C}^{D},\,V_{C}^{Q}\big{]}^{\mathsf{T}}$.

The input functions of the control affine system, $\mathbf{g}_{1}(\mathbf{x})$
and $\mathbf{g}_{2}(\mathbf{x})$, associated with each input are given by,

$\displaystyle\mathbf{g}_{1}$
$\displaystyle=\Big{[}\,0,\,0,\,0,\,0,\,-\cos{x_{2}},\,\sin{x_{2}},\,\frac{1}{C},\,0,\,-\frac{K_{p}^{Ireg}}{L},\,0\,\Big{]}^{\mathsf{T}},$
(21) $\displaystyle\mathbf{g}_{2}$
$\displaystyle=\Big{[}\,0,\,0,\,0,\,0,\,-\sin{x_{2}},\,-\cos{x_{2}},\,0,\,\frac{1}{C},\,0,\,-\frac{K_{p}^{Ireg}}{L}\,\Big{]}^{\mathsf{T}}.$
(22)

Figure 11: The connection circuit of DER Figure 12: Vf control

## References

  * [1] A. R. Sobbouhi and A. Vahedi, “Transient stability prediction of power system; a review on methods, classification and considerations,” _Electric Power Systems Research_ , vol. 190, p. 106853, 2021.
  * [2] X. Li, C. Mishra, S. Chen, Y. Wang, and J. De La Ree, “Determination of parameters of time-delayed embedding algorithm using koopman operator-based model predictive frequency control,” _CSEE Journal of Power and Energy Systems_ , vol. 7, no. 6, pp. 1140–1151, 2021.
  * [3] M. Adibi and J. van der Woude, “A reinforcement learning approach for frequency control of inverted-based microgrids,” _IFAC-PapersOnLine_ , vol. 52, no. 4, pp. 111–116, 2019.
  * [4] S. S. Madani and A. Karimi, “Data-driven lpv controller design for islanded microgrids,” _IFAC-PapersOnLine_ , vol. 54, no. 7, pp. 433–438, 2021.
  * [5] J. Sanchez-Gasca, D. Trudnowski, E. Barocio, J. Hauer, J. Pierre, C. Canizares, H. Huang, T. Rauhala, J. Chow, I. Kamwa, J. Turunen, M. Crow, G. Ledwich, L. Vanfretti, L. Dosiek, R. Martin, V. Vittal, H. Ghasemi, E. Martinez, D. Vowles, M. Gibbard, A. Messina, R. Wies, L. Haarla, B. Pal, and N. Zhou, _Identification of Electromechanical Modes in Power Systems_. United States: IEEE, 2012.
  * [6] J. Hauer, “Application of prony analysis to the determination of modal content and equivalent models for measured power system response,” _IEEE Transactions on Power Systems_ , vol. 6, no. 3, pp. 1062–1068, 1991.
  * [7] D. J. Trudnowski and J. W. Pierre, “Overview of algorithms for estimating swing modes from measured responses,” in _2009 IEEE Power Energy Society General Meeting_ , 2009, pp. 1–8.
  * [8] I. Kamwa, R. Grondin, E. Dickinson, and S. Fortin, “A minimal realization approach to reduced-order modelling and modal analysis for power system response signals,” _IEEE Transactions on Power Systems_ , vol. 8, no. 3, pp. 1020–1029, 1993.
  * [9] J. Sanchez-Gasca and J. Chow, “Performance comparison of three identification methods for the analysis of electromechanical oscillations,” _IEEE Transactions on Power Systems_ , vol. 14, no. 3, pp. 995–1002, 1999.
  * [10] L. L. Grant and M. L. Crow, “Comparison of matrix pencil and prony methods for power system modal analysis of noisy signals,” in _2011 North American Power Symposium_ , 2011, pp. 1–7.
  * [11] G. Liu, J. Quintero, and V. M. Venkatasubramanian, “Oscillation monitoring system based on wide area synchrophasors in power systems,” in _2007 iREP Symposium - Bulk Power System Dynamics and Control - VII. Revitalizing Operational Reliability_ , 2007, pp. 1–13.
  * [12] A. F. El Hamalawy, M. Ammar, H. F. Sindi, M. F. Shaaban, and H. H. Zeineldin, “A subspace identification technique for real time stability assessment of droop based microgrids,” _IEEE Transactions on Power Systems_ , pp. 1–1, 2021.
  * [13] E. Barocio, B. C. Pal, N. F. Thornhill, and A. R. Messina, “A dynamic mode decomposition framework for global power system oscillation analysis,” _IEEE Transactions on Power Systems_ , vol. 30, no. 6, pp. 2902–2912, 2015\.
  * [14] P. Stoica, R. L. Moses _et al._ , _Spectral analysis of signals_. Pearson Prentice Hall Upper Saddle River, NJ, 2005.
  * [15] M. G. Anderson, N. Zhou, J. W. Pierre, and R. W. Wies, “Bootstrap-based confidence interval estimates for electromechanical modes from multiple output analysis of measured ambient data,” _IEEE Transactions on Power Systems_ , vol. 20, no. 2, pp. 943–950, 2005.
  * [16] G. Liu and V. Venkatasubramanian, “Oscillation monitoring from ambient pmu measurements by frequency domain decomposition,” in _2008 IEEE International Symposium on Circuits and Systems_. IEEE, 2008, pp. 2821–2824.
  * [17] H. N. Villegas Pico, B. Mather, and G.-S. Seo, “Model identification of inverter nonlinear control dynamics,” in _2018 IEEE Electronic Power Grid (eGrid)_ , 2018, pp. 1–6.
  * [18] R. Fard, M. Karrari, and O. Malik, “Synchronous generator model identification for control application using volterra series,” _IEEE Transactions on Energy Conversion_ , vol. 20, no. 4, pp. 852–858, 2005.
  * [19] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proceedings of the national academy of sciences_ , vol. 113, no. 15, pp. 3932–3937, 2016.
  * [20] I. Mezić, “Spectral properties of dynamical systems, model reduction and decompositions,” _Nonlinear Dynamics_ , vol. 41, no. 1, pp. 309–325, 2005\.
  * [21] I. Mezić and A. Banaszuk, “Comparison of systems with complex behavior,” _Physica D: Nonlinear Phenomena_ , vol. 197, no. 1-2, pp. 101–133, 2004.
  * [22] I. Mezic, “Koopman operator, geometry, and learning of dynamical systems,” _Notices of the AMS_ , vol. 68, no. 7, pp. 1087–1105, 2021.
  * [23] P. Bevanda, S. Sosnowski, and S. Hirche, “Koopman operator dynamical models: Learning, analysis and control,” _Annual Reviews in Control_ , vol. 52, pp. 197–212, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1367578821000729
  * [24] M. Budišic, R. Mohr, and I. Mezic, “The koopman operator in systems and control: Concepts, methodologies, and applications,” 2020.
  * [25] Y. Susuki and I. Mezic, “Nonlinear koopman modes and coherency identification of coupled swing dynamics,” _IEEE Transactions on Power Systems_ , vol. 26, no. 4, pp. 1894–1904, 2011.
  * [26] M. Korda, Y. Susuki, and I. Mezić, “Power grid transient stabilization using koopman model predictive control,” _IFAC-PapersOnLine_ , vol. 51, no. 28, pp. 297–302, 2018.
  * [27] P. Sharma, V. Ajjarapu, and U. Vaidya, “Data-driven identification of nonlinear power system dynamics using output-only measurements,” _IEEE Transactions on Power Systems_ , 2021.
  * [28] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, “A data–driven approximation of the koopman operator: Extending dynamic mode decomposition,” _Journal of Nonlinear Science_ , vol. 25, no. 6, pp. 1307–1346, 2015.
  * [29] C. Folkestad, D. Pastor, I. Mezic, R. Mohr, M. Fonoberova, and J. Burdick, “Extended dynamic mode decomposition with learned koopman eigenfunctions for prediction and control,” in _2020 american control conference (acc)_. IEEE, 2020, pp. 3906–3913.
  * [30] M. Korda and I. Mezić, “Optimal construction of koopman eigenfunctions for prediction and control,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 12, pp. 5114–5129, 2020.
  * [31] D. Goswami and D. A. Paley, “Bilinearization, reachability, and optimal control of control-affine nonlinear systems: A koopman spectral approach,” _IEEE Transactions on Automatic Control_ , 2021.
  * [32] S. Peitz, S. E. Otto, and C. W. Rowley, “Data-driven model predictive control using interpolated koopman generators,” _SIAM Journal on Applied Dynamical Systems_ , vol. 19, no. 3, pp. 2162–2193, 2020.
  * [33] A. M. Kettner, L. Reyes-Chamorro, J. K. M. Becker, Z. Zou, M. Liserre, and M. Paolone, “Harmonic power-flow study of polyphase grids with converter-interfaced distributed energy resources—part i: Modeling framework and algorithm,” _IEEE Transactions on Smart Grid_ , vol. 13, no. 1, pp. 458–469, 2021.
  * [34] A. M. Kettner and M. Paolone, “On the properties of the power systems nodal admittance matrix,” _IEEE Transactions on Power Systems_ , vol. 33, no. 1, pp. 1130–1131, 2017.
  * [35] F. Dorfler and F. Bullo, “Kron reduction of graphs with applications to electrical networks,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 60, no. 1, pp. 150–163, 2012.
  * [36] A. Mauroy and I. Mezić, “Global stability analysis using the eigenfunctions of the koopman operator,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 11, pp. 3356–3369, 2016.
  * [37] A. L. Bruce, V. M. Zeidan, and D. S. Bernstein, “What is the koopman operator? a simplified treatment for discrete-time systems,” in _2019 American Control Conference (ACC)_. IEEE, 2019, pp. 1912–1917.
  * [38] J. L. Proctor, S. L. Brunton, and J. N. Kutz, “Dynamic mode decomposition with control,” _SIAM Journal on Applied Dynamical Systems_ , vol. 15, no. 1, pp. 142–161, 2016.
  * [39] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. Nathan Kutz, “On dynamic mode decomposition: Theory and applications,” _Journal of Computational Dynamics_ , vol. 1, no. 2, p. 391–421, 2014.
  * [40] D. Butler, “Facts about eigenvalues,” _Textbook, University of California, San diego, US_ , 2015.
  * [41] C. L. Lawson and R. J. Hanson, _Solving least squares problems_. SIAM, 1995.
  * [42] Y. Li, _Cyber-Physical Microgrids_. Springer, 2021.
  * [43] A. Surana, “Koopman operator based observer synthesis for control-affine nonlinear systems,” in _2016 IEEE 55th Conference on Decision and Control (CDC)_. IEEE, 2016, pp. 6492–6499.


# Learning the cost-to-go for mixed-integer nonlinear model predictive control

Christopher A. Orrico  W.P.M.H. Heemels  Dinesh Krishnamoorthy  Eindhoven
University of Technology, 5600 MB Eindhoven, The Netherlands (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>d.krishnamoorthy@tue.nl).  DIFFER -
Dutch Institute for Fundamental Energy Research, 5612 AJ Eindhoven, The
Netherlands (e-mail<EMAIL_ADDRESS>
###### Abstract

Application of nonlinear model predictive control (NMPC) to problems with
hybrid dynamical systems, disjoint constraints, or discrete controls often
results in mixed-integer formulations with both continuous and discrete
decision variables. However, solving mixed-integer nonlinear programming
problems (MINLP) in real-time is challenging, which can be a limiting factor
in many applications. To address the computational complexity of solving mixed
integer nonlinear model predictive control problem in real-time, this paper
proposes an approximate mixed integer NMPC formulation based on value function
approximation. Leveraging Bellman’s principle of optimality, the key idea here
is to divide the prediction horizon into two parts, where the optimal value
function of the latter part of the prediction horizon is approximated offline
using expert demonstrations. Doing so allows us to solve the MINMPC problem
with a considerably shorter prediction horizon online, thereby reducing the
online computation cost. The paper uses an inverted pendulum example with
discrete controls to illustrate this approach.

###### keywords:

Mixed integer model predictive control of hybrid systems, inverse optimization

Model predictive control has emerged as a powerful tool to control dynamic
systems in a wide range of application areas, due to its ability to explicitly
handle constraints and optimize desired performance criterion. This is done by
formulating the control problem as a discrete time finite horizon optimal
control problem, which is repeatedly solved online at each sampling instant,
given the current state of the system. the first control input from the
optimal trajectory is implemented in the system and the process is repeated.

As the domain of MPC application expands to new class of problems, this
requires more complex optimization problem formulations. Problem classes that
commonly arises in many application domains include, discrete actuators (where
the controls only take values from a finite set), switching systems (where the
system dynamics are discontinuous and non-smooth), problems with disjoint
constraint sets, avoidance constraints in path planning, or problems with
specified logics etc (Richards and How, 2005). MPC framework can be applied to
such problem classes by formulating the optimal control problem with both
continuous and integer decision variables. The resulting nonlinear MPC
framework with discrete decision variables require the online solution of
mixed integer nonlinear programming (MINLP) problems at each sampling instant.
However, MINLP problems in general are $\mathcal{NP}$-hard, which makes it
challenging for real-time control.

Model predictive control, as the name suggests, uses a model to predict the
effect of the control actions over a finite rolling horizon, and chooses an
optimal control sequence that optimizes a given performance criterion.
Ideally, we want a sufficiently long prediction horizon, since this would help
the controller to consider the effect of the control actions far into the
future, resulting in a better control policy. However, the size of the
optimization problem, and in turn the computational cost, increases with the
length of the prediction horizon. This is only further amplified in the
presence of discrete decision variables, where the number of nodes in the
branch and bound algorithms increases. Using a short prediction horizon
naively on the other hand leads to myopic policies that can result in
performance degradation.

This paper argues that in want for very short prediction horizons all is not
lost, and that value function approximations enable us to leverage Bellaman’s
principle of optimality, such that the online controller can be reduced to
even as small as a one-step-look-ahead controller without jeopardizing the
control performance. The prediction horizon in the original MINMPC problem can
be divided into two parts, where the optimal value function of the tail part
of the prediction horizon is replaced by a cost-to-go function. Since the
optimal cost-to-go function is not known exactly, this can be approximated
offline, and the online controller is solved with the approximate value
function as its cost-to-go. Replacing the optimal value function with
approximate cost-to-go functions are also commonly referred to as approximate
dynamic programming (ADP) Bertsekas (2011).

This paper uses inverse optimization to learn the value function from expert
demonstrations comprising of optimal state-action pairs (Keshavarz et al.,
2011). The expert demonstrations may either be from the original long horizon
MINMPC solved offline, or other approximate MINMPC strategies, or may even be
from human expert demonstrations. Learning the cost function from (near)
optimal state-action pairs is broadly studied under the context of learning
from demonstrations (LfD). Unlike other learning from demonstration approaches
such as direct policy fitting/imitation learning, learning the cost function
via inverse optimization requires considerably much lesser data that need not
cover the entire feasible state space. To this end, the main contribution of
this paper is a mixed integer NMPC formulation based on value function
approximation, where we show that the myopic mixed integer NMPC can be
formulated with very short prediction horizons that results in comparable
performance as a long prediction horizon, at a fraction of the online
computation cost.

#### Related work

Various strategies and heuristics have been proposed in the literature to
address real-time mixed-integer Nonlinear Model Predictive Control (NMPC)
problems. These include solving the relaxed problem by removing the
integrality constraint and obtaining the integer trajectory through rounding
schemes like Sum-Up-Rounding (Sager, 2009). Another approach involves a
penalty term homotopy method, where the integer variable is replaced with the
constraint $z(1-z)\leq\beta,\;z\in[0,1]$ along with the homotopy
$\beta\rightarrow 0^{+}$ (Sager et al., 2006). Nonlinear functions are often
approximated using piecewise linearization to transform the Mixed-Integer
Nonlinear Programming (MINLP) into Mixed-Integer Linear Programming (MILP) or
Mixed-Integer Quadratic Programming (MIQP). However, such approximations may
not always be applicable. Real-time iteration schemes for MINMPC based on
outer convexification and rounding schemes have been extended (De Mauri et
al., 2020), though this becomes challenging if discrete variables appear in
inequality constraints. Recent interest in deep learning has led to direct
policy approximation approaches Karg and Lucia (2018), where the MPC policy is
approximated using deep neural networks trained on large amounts of offline-
generated data. However, this requires extensive data covering the entire
feasible state space and poses challenges for any online updates and changes
in the controller.

##  1 Problem formulation

Consider a mixed integer nonlinear MPC problem

$\displaystyle V^{*}(x(t))=\min\;$
$\displaystyle\sum_{k=0}^{N-1}\ell(x_{k},u_{k},z_{k})+\ell_{N}(x_{N})$ (1a)
s.t. $\displaystyle x_{k+1}=f(x_{k},u_{k},z_{k})$ (1b) $\displaystyle
u_{k}\in\mathcal{U},\quad z_{k}\in\mathcal{Z}$ (1c) $\displaystyle x_{0}=x(t)$
(1d)

where $x\in\mathbb{R}^{n_{x}}$ are the set of states,
$u\in\mathcal{U}\subseteq\mathbb{R}^{n_{u}}$ denotes the set of continuous
controls, $z\in\mathcal{Z}\subseteq\mathbb{Z}^{n_{z}}$ denotes the set of
discrete controls, $N$ denotes the length of the prediction horizon,
$\ell:\mathcal{X}\times\mathcal{U}\times\mathcal{Z}\rightarrow\mathbb{R}$ and
$\ell_{N}:\mathcal{X}\rightarrow\mathbb{R}$ denote the stage and terminal
costs respectively, and
$f:\mathcal{X}\times\mathcal{U}\times\mathcal{Z}\rightarrow\mathbb{R}^{n_{x}}$
denotes the system dynamics. Solving the MINLP (1) at each sampling instant
and implementing the first control elements $[u_{0},z_{0}]^{\mathsf{T}}$
results in the mixed integer MPC policy
$[u(t),z(t)]^{\mathsf{T}}=\pi_{\text{mpc}}(x(t))$.

Problem setting: The MINMPC problem (1) is designed to be optimal, but
computationally intensive and is difficult to solve online within the desired
sampling time.

Objective: Using the MINMPC problem (1) as the “expert” that can be queried
offline, we want to build a simpler control policy that is amenable for online
implementation.

The optimal value function of the OCP starting from $x_{1}$ can be denoted by

$\displaystyle V_{\eta}=\min\;$
$\displaystyle\sum_{k=1}^{N-1}\ell(x_{k},u_{k},z_{k})+\ell_{N}(x_{N})$ (2a)
s.t. $\displaystyle x_{k+1}=f(x_{k},u_{k},z_{k})$ (2b) $\displaystyle
g(x_{k},u_{k},z_{k})\leq 0$ (2c) $\displaystyle\qquad\forall k=1,\dots,N-1$

That is $V_{\eta}$ is the optimal-cost-to-go when starting from state $x_{1}$.
Then by Bellman’s principle of optimality, we can truncate the prediction
horizon and equivalently solve

$\displaystyle\min_{u,z}\;$
$\displaystyle\ell(x(t),u,z)+V_{\eta}(f(x(t),u,z))$ (3a) s.t. $\displaystyle
g(x(t),u,z)\leq 0$ (3b)

which is significantly easier to solve than the original MINMPC problem (1).
By simply reducing the length of the prediction horizon, and consequently, the
number of decision variables, mixed integer problems can be solved more
efficiently.

The challenge with the DP recursion is that calculating the optimal cost-to-go
$V_{\eta}$ from all possible states $x\in\mathcal{X}$ can be intractable and
prohibitively time consuming. This can be addressed using a class of methods
known as approximate dynamic programming (ADP), where one can approximate the
optimal cost-to-go function $V_{\eta}$ in (3) by a convex function
approximation $\mathcal{V}(x):=x^{\mathsf{T}}\mathbf{P}x$, parameterized by
$\mathbf{P}$. The task is then to find $\mathbf{P}$ such that the myopic
policy obtained by solving (3) mimics the full horizon MINMPC policy (1). To
do this, we first generate a set of optimal state-action pairs
$\mathcal{D}:=\\{(x_{i},[u_{i},z_{i}]^{\mathsf{T}})\\}_{i=1}^{M}$ by solving
the full horizon MINMPC offline. This dataset is then used to impute the
parameter $\mathbf{P}$ of the myopic policy, such that the dataset
$\mathcal{D}$ is approximately consistent with the data set $\mathcal{D}$. To
do this, we use inverse optimization.

#### Imputing the cost-to-go function:

Denoting the set of decision variables of the myopic MPC (3) by
$\mathbf{w}:=[u,z]^{\mathsf{T}}$, the KKT condition can be expressed as,

$\displaystyle\nabla_{\mathbf{w}}\ell(x,\mathbf{w})+2f(x,\mathbf{w})^{\mathsf{T}}\mathbf{P}f(x,\mathbf{w})+\lambda^{\mathsf{T}}\nabla_{\mathbf{w}}g(x,\mathbf{w})$
$\displaystyle=0$ $\displaystyle g(x,\mathbf{w})$ $\displaystyle\leq 0$
$\displaystyle\lambda^{\mathsf{T}}g(x,\mathbf{w})$ $\displaystyle=0$
$\displaystyle\lambda$ $\displaystyle\geq 0$

Suppose the optimal state-action data pairs
$\mathcal{D}:=\\{(x_{i},\mathbf{w}_{i})\\}:=\\{(x_{i},[u_{i},z_{i}]^{\mathsf{T}})\\}$
obtained by querying the full MINMPC (1) satisfies the KKT conditions of the
myopic MPC stated above, then the myopic MPC policy is said to be consistent
with the dataset $\mathcal{D}$. In other words, we want to find $\mathbf{P}$,
such the the KKT conditions evaluated at the datapoints
$\\{(x_{i},\mathbf{w}_{i})\\}$ are close to zero.

To this end, the problem of imputing the cost-to-go function can be formulated
as

$\displaystyle\min_{P,\\{\lambda_{i}\\}}$
$\displaystyle~{}\sum{}_{i=1}^{M}\|r_{\textrm{stat}}(x_{i},\mathbf{w}_{i})\|_{2}^{2}+\|r_{\textrm{comp}}(x_{i},\mathbf{w}_{i})\|_{2}^{2}$
(4a) s.t. $\displaystyle P\succeq 0,$ (4b) $\displaystyle\lambda_{i}\succeq
0,\;\;\forall i=1,\dots,M$ (4c)

where KKT residuals for the $i^{th}$ datapoint is defined as

$\displaystyle r_{\textrm{stat}}(x_{i},\mathbf{w}_{i})$
$\displaystyle=\nabla_{\mathbf{w}}\ell(x_{i},\mathbf{w}_{i})+2f(x_{i},\mathbf{w}_{i})^{\mathsf{T}}\mathbf{P}f(x_{i},\mathbf{w}_{i})$
$\displaystyle\qquad+\lambda_{i}^{\mathsf{T}}\nabla_{\mathbf{w}}g(x_{i},\mathbf{w}_{i})$
(5a) $\displaystyle r_{\textrm{comp}}(x_{i},\mathbf{w}_{i})$
$\displaystyle=\,\lambda_{i}^{\mathsf{T}}g(x_{i},\mathbf{w}_{i})$ (5b)

Note that offline problem of imputing the cost function (4) is a semi-definite
program (SDP) which can be solved efficiently. The solution $P$ can then be
used in a myopic MPC controller (3) with a prediction horizon as small as
$N=1$ that approximates the control decisions of the original long prediction
horizon MINMPC controller (1).

##  2 Illustrative example

To demonstrate the efficacy of myopic MPC with IOC-imputed cost-to-go
weighting parameters, we apply this method to the Lokta-Volterra fishing
problem, which sager2006 present as a foundational benchmark problem for
mixed-integer control of hybrid nonlinear systems. This benchmark problem
seeks to bring the oscillations of predator and prey populations close to a
steady-state solution through an optimal fishing allowance. The dynamic system
(6a) describes the evolution of the prey population $x_{1}$ and the predator
population $x_{2}$ as

$\displaystyle\begin{bmatrix}\dot{x_{1}}\\\
\dot{x_{2}}\end{bmatrix}=\begin{bmatrix}x_{1}-x_{1}x_{2}-c_{1}x_{1}u\\\
-x_{2}+x_{1}x_{2}-c_{2}x_{2}u\end{bmatrix},$ (6a) $\displaystyle
u\in\\{0,1\\},$ (6b) $\displaystyle[x_{1},x_{2}]^{\top}\geq 0.$ (6c)

where the growth of both populations are coupled and partly controlled by
fishing rates governed with coefficients $c_{1}$ and $c_{2}$. The system
includes a binary input variable $z_{k}=u$ (6b) representing the decision to
fish ($u=1$) or not to fish ($u=0$) and a state constraint (6c) requiring the
fish populations to be positive.

We first formulate this decision-making problem as a mixed-integer optimal
control problem of the form(1). The controller objective parameters are tuned
to achieve good reference tracking for a prediction horizon length $N=??$ that
approximates the solution of $N=\infty$ with discrete sampling time $t_{s}$.
We construct the controller in Matlab using the CasADi toolbox (Andersson et
al., In Press, 2018) and the Bonmin mixed integer nonlinear programming solver
(Bonami et al., 2008). The resulting solution set size $M=120$ required
$t_{\textrm{CPU}}=2.07\times 10^{3}$ (s) to solve on an Intel(R) CoreTM i5
processor running at 2.40 GHz with 15.7 GB of useable RAM.

Since this is not amenable for online implementation, we aim to impute a
convex cost-to-go function $\mathcal{V}(x):=x^{\mathsf{T}}\mathbf{P}x$, such
that the we can solve a myopic MPC online with $N=1$ and $\mathcal{V}(x)$ as
the cost-to-go function. To compute $\mathbf{P}$, the MINMPC controller then
is run in an offline feedback control simulation with perfect state
observation, for three initial states, resulting in 3 trajectories of data
$\mathcal{D}$.

###  2.1 Imputing the cost-to-go offline

Using the state-action data set $\mathcal{D}$ obtained from the offline MINMPC
simulations, we solve the SDP for $P$ that minimizes $r_{stat}^{(m)}$ and
$r_{comp}^{(m)}$ (4), yielding

$\displaystyle P=\begin{bmatrix}3.07\times 10^{-4}&2.13\times 10^{-5}\\\
2.13\times 10^{-5}&2.06\times 10^{-3}\end{bmatrix},$

and $\|r_{stat}\|_{\infty}=1.29\times 10^{-6}$ and
$\|r_{comp}\|_{\infty}=1.91\times 10^{-6}$. As required in (4), the $P$ matrix
is positive definite. Moreover, the maximum residuals for both the
stationarity and complementary slackness KKT conditions are vanishingly small.
We conclude, therefore, that the imputed cost-to-go parameter $P$ is
consistent with all observed solutions to the Lotka-Volterra fishing MINMPC
problem. Consequently, the $P$ matrix is appropriate as a tuning matrix to
build the myopic MPC controller in (3).

###  2.2 Online Controller Performance

The goal for the MPC controller is to track a state reference without
violating system constraints. Within the context of this work, however, we
introduce an additional key performance goal of the maximum computation time
per controller decision. The controller performance is tested in a discrete
time simulation of the dynamic system, with added plant-model mismatch
(implemented as $10\%$ error in the fishing coefficients $c_{1}$ and $c_{2}$
of (6a)) and synthetic measurement noise (implemented as a zero-mean Gaussian
white noise). The results of the control simulations for the Lotka-Volterra
fishing problem (including the offline expert demonstrations) are given in
Fig. 1.

As shown in Fig. 1, the myopic MPC controller using a $1$ step prediction
horizon and the cost-to-go parameter imputed using IOC closely reproduces the
control actions computed by the full horizon MINMPC controller with slight
deviation. It is worth noting here the MINMPC controller does not drive the
state to the state reference unless $N\gg 1$, meaning that the imputed cost-
to-go parameter is responsible for controller performance.

Figure 1: (a) Prey and (b) predator state population for the three training
solution set trajectories (in light blue) computed offline (without plant-
model mismatch and measurement noise), the trajectory computed with MINMPC
controller (in blue), and the trajectory computed with the myopic MPC
controller (in red). The reference is indicated by a black, dashed line. (c)
The control decisions for both the training sets and the online controllers.
(d) The computation time per controller decision.

While the ability of the myopic MPC controller to achieve reference tracking
performance that approximates that of the long prediction horizon controller
is satisfactory, what is more striking is the reduction in computation time
required to solve the control problem. Fig. 1d shows the controller
computation time at each simulation step. The MINMPC controller was very
computationally expensive, requiring a maximum computation time of
$t_{\textrm{CPU}}=217$ (s) per decision to solve. However, the maximum
computation time of the myopic MPC controller was only $t_{\textrm{CPU}}=54.1$
(ms). By greatly reducing the control problem complexity, the myopic MPC
controller was able to reduce the maximum computation time per controller
decision by more than three orders of magnitude over the MINMPC controller.
Hence, the myopic MPC method takes a controller that is computationally
intractable and makes it solvable in real time. The simulations were performed
on an Intel(R) CoreTM i5 processor running at 2.40 GHz with 15.7 GB of useable
RAM.

##  3 Conclusion

In this work, we have demonstrated that the KKT residual minimization method
for IOC can impute the cost-to-go weighting parameter of a $1$-step prediction
horizon MPC problem from a set of action-state pair solutions computed offline
by a long prediction horizon MINMPC controller. This parameter is shown to be
consistent with the offline solution set through the definition of approximate
optimality. We then show that the myopic MPC controller using this cost-to-go
parameter closely replicates the reference tracking performance of the
original MPC controller. The method that we demonstrate in this paper offers a
powerful tool for simplifying computationally complex MINMPC controllers.

There remain several open questions that this proof-of-concept demonstration
raises. While we impute a parameter for a quadratic cost-to-go term because it
leads to an IOC that is solvable as an SDP, a quadratic term may not be
adequate approximate the cost-to-go of certain MPC problems. Bellman’s
principle of optimality indicates that some function must exist that is
consistent with a given solution set of expert controller demonstrations.
Broadening the definition of the imputed objective term to any convex function
and solving the same IOC problem would broaden the applicability of the method
described here. Additionally, complete characterization of stability and
recursive feasibility properties of the myopic MPC controller is outside of
the scope of this paper. Nonetheless, the authors consider such a
characterization invaluable and plan to explore it with the future work.

Lastly, the question of computational costs of the offline solution set is not
addressed in this work. While the cost of computing offline solutions was
acceptable for the two simple benchmark problems shown here, this may not be
the case for control problems with larger state spaces, not to mention hybrid
systems with combinatorial integer actuators. The authors seek to combine the
imputed objective method shown here with methods to make offline data
generation more efficient. One example of this may be data augmentation method
demonstrated by Krishnamoorthy (2021). Answering the above open research
questions would hopefully yield a generalized approach to reducing controller
complexity with stability and feasibility guarantees from efficiently
generated training data for any MPC problem.

## References

  * Andersson et al. (In Press, 2018) Andersson, J.A.E., Gillis, J., Horn, G., Rawlings, J.B., and Diehl, M. (In Press, 2018).  CasADi – A software framework for nonlinear optimization and optimal control.  _Mathematical Programming Computation_.
  * Bertsekas (2011) Bertsekas, D.P. (2011).  Dynamic programming and optimal control 3rd edition, volume ii.  _Belmont, MA: Athena Scientific_.
  * Bonami et al. (2008) Bonami, P., Biegler, L.T., Conn, A.R., Cornuéjols, G., Grossmann, I.E., Laird, C.D., Lee, J., Lodi, A., Margot, F., Sawaya, N., et al. (2008).  An algorithmic framework for convex mixed integer nonlinear programs.  _Discrete optimization_ , 5(2), 186–204.
  * De Mauri et al. (2020) De Mauri, M., Van Roy, W., Gillis, J., Swevers, J., and Pipeleers, G. (2020).  Real time iterations for mixed-integer model predictive control.  In _2020 European Control Conference (ECC)_ , 699–705. IEEE.
  * Karg and Lucia (2018) Karg, B. and Lucia, S. (2018).  Deep learning-based embedded mixed-integer model predictive control.  In _2018 European Control Conference (ECC)_ , 2075–2080. IEEE.
  * Keshavarz et al. (2011) Keshavarz, A., Wang, Y., and Boyd, S. (2011).  Imputing a convex objective function.  In _2011 IEEE international symposium on intelligent control_ , 613–619. IEEE.
  * Krishnamoorthy (2021) Krishnamoorthy, D. (2021).  A sensitivity-based data augmentation framework for model predictive control policy approximation.  _IEEE Transactions on Automatic Control_ , In-Press.
  * Richards and How (2005) Richards, A. and How, J. (2005).  Mixed-integer programming for control.  In _American Control Conference, 2005. Proceedings of the 2005_ , 2676–2683. IEEE.
  * Sager (2009) Sager, S. (2009).  Reformulations and algorithms for the optimization of switching decisions in nonlinear optimal control.  _Journal of Process Control_ , 19(8), 1238–1247.
  * Sager et al. (2006) Sager, S., Bock, H.G., Diehl, M., Reinelt, G., and Schloder, J.P. (2006).  Numerical methods for optimal control with binary control functions applied to a lotka-volterra type fishing problem.  In _Recent Advances in Optimization_ , 269–289. Springer.