# Data-driven model predictive control of battery storage units

Johannes B. Lipka1 and Christian A. Hans2 1Johannes B. Lipka is with Siemens
AG, Germany<EMAIL_ADDRESS>A. Hans is with the
University of Kassel, Germany<EMAIL_ADDRESS>
###### Abstract

In many state-of-the-art control approaches for power systems with storage
units, an explicit model of the storage dynamics is required. With growing
numbers of storage units, identifying these dynamics can be cumbersome. This
paper employs recent data-driven control approaches that do not require an
explicit identification step. Instead, they use measured input/output data in
control formulations. In detail, we propose an economic data-driven model
predictive control (MPC) scheme to operate a small power system with input-
nonlinear battery dynamics. First, a linear data-driven MPC approach that uses
a slack variable to account for plant-model-mismatch is proposed. In a second
step, an input-nonlinear data-driven MPC scheme is deduced. Comparisons with a
reference indicate that the linear data-driven MPC approximates the nonlinear
plant in an acceptable manner. Even better results, however, can be obtained
with the input-nonlinear data-driven MPC scheme which provides increased
prediction accuracy.

##  I Introduction

Battery storage units are integral to future power systems. They exhibit
diverse dynamics based on technology and production. Moreover, their dynamics
may change over time due to ageing. Effective control of future power systems
with storage units necessitates consideration of these dynamics, typically
through a system identification step before control synthesis. In systems with
multiple batteries, each with different dynamics, identifying and updating
their models can be cumbersome. This motivates alternative data-driven control
schemes which do not require an explicit identification step.

Willems et al. [1] provided foundations for modeling of linear time-invariant
(LTI) systems using input/output (I/O) data. In the context of behavioural
systems, they proved that if the input signal of an LTI system is persistently
exciting (see Section II), then the Hankel matrix of sufficiently many I/O
measurements can span the vector space of all possible trajectories of the
system. The authors of [2, 3, 4, 5] extended the results of [1] to LTI state
models and to feedback-linearizable Hammerstein and Wiener systems. In [6, 7,
8] the theory was further extended to unknown LTI systems, ensuring trajectory
tracking and constraint satisfaction under some easy-to-hold assumptions. The
authors of [9] presented a trajectory-based MPC scheme with terminal
constraints and could find conditions for exponential stability in presence of
measurement noise. In [10, 11], these conditions were extended for slow
changing nonlinear systems.

Despite recent advancements in data-driven control theory and their
application to practical settings, the utilization of data-driven approaches
in economic MPC schemes has just recently started (see [12, 13, 14]). To the
authors’ best knowledge, data-driven MPC also has not been applied to plants
with battery storage systems. Furthermore, there exists only a small number of
data-driven control approaches for power systems [6, 12, 15, 16]. Moreover,
the extended fundamental lemma from [2] has been rarely employed in an MPC
context up to know.

In this paper, we aim to further bridge the gap between data-driven control
theory and applications in the power systems domain. Our contributions
include:

  * •

Developing different economic data-driven MPC approaches to control an
islanded power system with nonlinear Hammerstein-type battery dynamics. In
detail, we propose one data-driven MPC approach that considers LTI dynamics
and one which considers input-nonlinearities with a known structure.

  * •

Showcasing how the extended fundamental lemma of [2] can be practically
employed in an MPC formulation to control battery storage systems.

  * •

Quantitatively comparing the prediction capabilities of the different MPC
approaches with a reference controller that has perfect knowledge of the
system’s nonlinearities in case studies. More specifically, we evaluate the
prediction accuracy and associated constraint violations of the proposed MPC
schemes.

The paper is structured as follows. Section II focuses on mathematical
preliminaries. In Section III, the running example of an islanded grid with
input-nonlinear battery dynamics is introduced. The cost function and a
reference MPC are discussed in Section IV. A linear data-driven MPC scheme is
presented in Section V. In Section VI, a nonlinear data-driven MPC scheme that
deals with Hammerstein-type system dynamics of known structure is presented.
Section VII concludes the work.

##  II Mathematical Preliminaries

In what follows, first the notation is discussed. Then, Hankel matrices and
persistence of excitation are defined and Willems’ fundamental lemma for LTI
systems is presented. Finally, Hammerstein systems are introduced and the
extension of Willems’ fundamental lemma for feedback-linearizable systems is
recalled.

###  II-A Notation

Let the set of real numbers, nonnegative real numbers, positive real numbers
and negative real numbers be $\mathbb{R}$, $\mathbb{R}_{\geq 0}$,
$\mathbb{R}_{>0}$ and $\mathbb{R}_{<0}$, respectively. Integers, nonnegative
integers and positive integers are denoted by $\mathbb{Z}$, $\mathbb{N}_{0}$
and $\mathbb{N}$, respectively. Let $\mathbb{X}=\\{x(k)\\}_{k=k_{A}}^{k_{B}}$
denote a sequence, i.e., an enumerated and ordered set
$\\{x(k_{A}),\dots,x(k_{B})\\}$ of elements $x(k)\in\mathbb{R}^{n}$ with
$k_{A}\in\mathbb{N}_{0}$, $k_{B}\in\mathbb{N}_{0}$ and $k_{A}\leq k_{B}$. The
notation $x(k|t)$ is used to refer to a prediction performed at time instant
$t\in\mathbb{N}_{0}$ for prediction step $k\in\mathbb{N}_{0}$, i.e., time
instant $t+k$. $X=[x(k)]_{k=k_{A}}^{k_{B}}$ is shorthand for
$[x(k_{A})^{T}~{}\cdots~{}x(k_{B})^{T}]^{T}$. Moreover, $|a|$ refers to the
absolute value of $a\in\mathbb{R}$.

###  II-B Hankel Matrix

A sequence $\mathbb{X}=\\{x(k)\\}_{k=0}^{N-1}$ of elements
$x(k)\in\mathbb{R}^{n}$ of length $N\in\mathbb{N}$ can induce a Hankel matrix
of order $L\in\mathbb{N}$, $L<N$, of the form

$H_{L}(\mathbb{X}):=\begin{bmatrix}x(0)&x(1)&\cdots&x(N-L)\\\
x(1)&x(2)&\cdots&x(N-L+1)\\\ \vdots&\vdots&\ddots&\vdots\\\
x(L-1)&x(L)&\cdots&x(N-1)\end{bmatrix}.$ (1)

In detail, $H_{L}(\mathbb{X})$ is constructed by cutting the sequence
$\mathbb{X}$ into $N-L+1$ snippets of length $L$. Each column of
$H_{L}(\mathbb{X})$ represents a snippet which is shifted one time step ahead
compared to the column left of it. To proceed with what follows, we recall [2,
Definition 1].

######  Definition 1 (Persistence of Excitation)

A sequence $\mathbb{X}$=$\\{x(k)\\}_{k=0}^{N-1}$ with $x(k)\in\mathbb{R}^{n}$
is persistently exciting of order $L$ if

$\operatorname{rank}(H_{L}(\mathbb{X}))=nL\hskip 2.84544pt.$ (2)

###### Remark 1

If a sequence is persistently exciting of order $L$, then it is also
persistently exciting of order $\tilde{L}\leq L$ [9].

###  II-C Data-Driven Representation of LTI Systems

Let $G$ be a controllable LTI system and

$\displaystyle x(k+1)$ $\displaystyle=Ax(k)+Bu(k),$ (3a) $\displaystyle y(k)$
$\displaystyle=Cx(k)+Du(k),$ (3b)

with $x(0)=x_{0}$ be a minimal realization of $G$. Here,
$x(k)\in\mathbb{R}^{n}$ is the state vector, $y(k)\in\mathbb{R}^{p}$ the
output vector and $u(k)\in\mathbb{R}^{d}$ the input vector. Measuring the
discrete-time input and output signals of (3) for $N$ consecutive times allows
us to form the input and output sequences
$\mathbb{U}^{m}=\\{u^{m}(k)\\}^{N-1}_{k=0}$ and
$\mathbb{Y}^{m}=\\{y^{m}(k)\\}^{N-1}_{k=0}$. If $\mathbb{U}^{m}$ is
persistently exciting of order $L$, then the following theorem from [2,
Theorem 3] which is based on [1] holds.

###### Theorem 1

An input sequence $\\{u(k)\\}^{L-1}_{k=0}$ and its corresponding output
sequence $\\{y(k)\\}^{L-1}_{k=0}$ represent a trajectory of $G$ if and only if
there exists a vector $\alpha\in\mathbb{R}^{N-L+1}$ such that with
$\mathbf{U}=[u(k)]^{L-1}_{k=0}$ and $\mathbf{Y}=[y(k)]^{L-1}_{k=0}$ it holds
that

$\begin{bmatrix}\mathbf{U}\\\
\mathbf{Y}\end{bmatrix}=\begin{bmatrix}H_{L}(\mathbb{U}^{m})\\\
H_{L}(\mathbb{Y}^{m})\end{bmatrix}\alpha.$ (4)

In other words, the space of all possible I/O trajectories of $G$ is spanned
by the matrix $\begin{bmatrix}H_{L}(\mathbb{U}^{m})^{T}\hskip
8.5359ptH_{L}(\mathbb{Y}^{m})^{T}\end{bmatrix}^{T}$. Thus, the I/O behavior of
$G$ can be described by time-shifted windows of measured I/O signals supposed
that $\mathbb{U}^{m}$ is persistently exciting. Note that the persistence of
excitation criterion does not need to be checked for the output sequence
$\mathbb{Y}^{m}$ [1, 2].

When forecasting future trajectories of a system, the initial conditions need
to be included. In a minimal realization of an LTI system $G$ of order
$n\in\mathbb{N}$, the associated input and output sequences $\mathbb{U}$ and
$\mathbb{Y}$ of length larger than or equal to $n$ induce a unique internal
state trajectory [2]. This property can be used to describe the initial state:
Let $\tilde{n}\geq n$, $\tilde{n}\in\mathbb{N}$ represent an upper
approximation of the order of $G$. Then, an input sequence $\mathbb{U}^{m}$
and its corresponding output sequence $\mathbb{Y}^{m}$ of length $\tilde{n}$
(or larger) describe a unique initial state of $G$.

For example, the initial conditions $x(t)=x_{t}$ of a minimal realization of
$G$ can be described by the $\tilde{n}$ most recent I/O measurements at time
$t\geq\tilde{n}$. That is, the sequences
$\mathbb{U}_{t}^{m}=\\{u^{m}(k|t)\\}^{-1}_{k=-\tilde{n}}$ and
$\mathbb{Y}_{t}^{m}=\\{y^{m}(k|t)\\}^{-1}_{k=-\tilde{n}}$ induce the unique
initial state $x(t)$. In conclusion, we can adequately describe the behavior
and initial state of an LTI system using I/O data and an upper approximation
of its order $n$ [2].

###  II-D Data-Driven Representation of Hammerstein Systems

Consider a Hammerstein system of the form

$\displaystyle x(k+1)$ $\displaystyle=Ax(k)+B\psi(u(k)),$ (5a) $\displaystyle
y(k)$ $\displaystyle=Cx(k)+D\psi(u(k)),$ (5b)

with $x(0)=x_{0}$ where $\psi:\mathbb{R}^{d}\rightarrow\mathbb{R}^{\tilde{d}}$
is a nonlinear mapping of input $u$ from $\mathbb{R}^{d}$ to
$\mathbb{R}^{\tilde{d}}$ . If $\tilde{d}=1$, then $\psi$ takes the form

$\psi(u)=\alpha_{1}\psi_{1}(u)+\alpha_{2}\psi_{2}(u)+\ldots+\alpha_{r}\psi_{r}(u)$
(6)

where $\psi_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ are linear and nonlinear
mappings and $\alpha_{i}\in\mathbb{R}$ coefficients which are nonzero for at
least one $i\in\\{1,\dots,r\\}$ [2]. Let $G$ be a Hammerstein system of the
form (5). For simplicity, let $\tilde{d}=1$ such that $\psi$ can be described
by (6). Then, a new auxiliary linear input vector

$v(k)=\begin{bmatrix}\psi_{1}(u(k))&\psi_{2}(u(k))&\cdots&\psi_{r}(u(k))\end{bmatrix}^{T}$
(7)

can be formed. This vector can then be used in following theorem from [2,
Proposition 5].

###### Theorem 2

Suppose that $\mathbb{U}^{m}=\\{u^{m}(k)\\}^{N-1}_{k=0}$ and
$\mathbb{Y}^{m}=\\{y^{m}(k)\\}^{N-1}_{k=0}$ are measured I/O sequences of the
Hammerstein system (5) with order $n$. Moreover, let
$\mathbb{V}^{m}=\\{v^{m}(k)\\}^{N-1}_{k=0}$ be the corresponding linear input
sequence with
$v(k)=[\psi_{1}(u^{m}(k))^{T}~{}\psi_{2}(u^{m}(k))^{T}~{}\cdots~{}\psi_{r}(u^{m}(k))^{T}]^{T}$.
If $\mathbb{V}^{m}$ is persistently exciting of order $L$, then an I/O
sequence $\\{v(k)\\}^{L-1}_{k=0}$, $\\{y(k)\\}^{L-1}_{k=0}$ is a trajectory of
system (5) if and only if there exists a vector $\alpha\in\mathbb{R}^{N-L+1}$
such that

$\begin{bmatrix}\mathbf{V}\\\
\mathbf{Y}\end{bmatrix}=\begin{bmatrix}H_{L}(\mathbb{V}^{m})\\\
H_{L}(\mathbb{Y}^{m})\end{bmatrix}\alpha$ (8)

where $\mathbf{Y}=[y(k)]^{L-1}_{k=0}$ and $\mathbf{V}=[v(k)]^{L-1}_{k=0}$ with
$v(k)$ from (7).

##  III Model of an islanded grid

The following section presents a minimal example of an islanded grid which is
considered throughout this work to illustrate how data-driven models of
battery dynamics can be employed. Note that the presented results are not
restricted to this structure and can be easily extended to more complex ones.
Figure 1 depicts the grid consisting of a load, an energy storage unit, a
conventional and renewable generator with power $w_{d}$, $p_{s}$, $p_{t}$,
$p_{r}$, respectively. Arrows indicate power direction: positive values
describe power provided to the grid, negative values describe power consumed.

Figure 1: Islanded grid composed of storage unit, conventional generator,
renewable generator and a load. Layout from [17].

###  III-A Renewable Generator

The power of the renewable generator can be modeled by

$p^{min}_{r}\leq p_{r}(k)\leq w_{r}(k),$ (9)

where $w_{r}(k)\in\mathbb{R}$ is the weather-dependent available renewable
power at time-step $k\in\mathbb{N}_{0}$ and $p^{min}_{r}\in\mathbb{R}_{\geq
0}$ the minimum power. Note that in (9) the renewable infeed $p_{r}(k)$ is
assumed to be curtailable to values below the available power $w_{r}(k)$.

###  III-B Conventional Unit

We consider a conventional unit whose change in output power $\Delta
p_{t}(k)=p_{t}(k)-p_{t}(k-1)$ is not restricted. The on/off status of the unit
is represented by a binary variable $\delta(k)$, i.e.,

$\delta(k)=\begin{cases}0,&\text{if the conventional unit is disabled},\\\
1,&\text{if the conventional unit is enabled}.\end{cases}$

If $\delta(k)=0$, then $p_{t}(k)$ is zero. Otherwise, $p_{t}(k)$ is bounded by
the minimum power $p^{min}_{t}$$\in$$\mathbb{R}_{\geq 0}$ and the maximum
power $p^{max}_{t}$$\in$$\mathbb{R}_{>0}$. This can be considered using the
constraint

$p^{min}_{t}\delta(k)\leq p_{t}(k)\leq p^{max}_{t}\delta(k).$ (10)

###  III-C Energy Storage Unit

Given the availability of economical battery storage units, we focus our
modeling efforts on these types of storages. In detail, we consider storage
units whose dynamics can be modeled along the lines of [17, 18, 19]. Storage
units have physical limitations: the power output is constrained to

$\displaystyle p^{min}_{s}\leq p_{s}(k)\leq p^{max}_{s}$ (11)

with $p^{min}_{s}\in\mathbb{R}_{<0}$ and $p^{max}_{s}\in\mathbb{R}_{>0}$.
Recall that if $p_{s}(k)<0$, the battery charges and if $p_{s}(k)>0$, it
discharges. Let $x(k)$ denote the stored energy, which is bounded by

$x^{min}\leq x(k)\leq x^{max}$ (12)

with $x^{min}\in\mathbb{R}_{\geq 0}$ and $x^{max}\in\mathbb{R}_{>0}$. We
consider nonlinear storage dynamics. In detail, we use a quadratic loss-term
$p_{s}(k)^{2}$ in the dynamics

$\displaystyle x(k+1)$ $\displaystyle=Ax(k)+B_{l}p_{s}(k)+B_{q}p_{s}(k)^{2},$
(13a) $\displaystyle y(k)$ $\displaystyle=x(k),$ (13b)

to capture conversion losses (AC to DC and vice versa) when charging or
discharging [20, 21]. Matrix $A\in\mathbb{R}_{>0}$ is typically used to model
the self discharge. In realistic settings we have that $0\ll A<1$. Moreover,
the input matrices $B_{l}\in\mathbb{R}_{<0}$ and $B_{q}\in\mathbb{R}_{<0}$ are
negative. Note that the dynamics are state-affine but not input-affine.

###  III-D Power Equilibrium

In an islanded grid, the provided power must equal the consumed power. This
can be modeled via

$p_{t}(k)+p_{s}(k)+p_{r}(k)+w_{d}(k)=0.$ (14)

##  IV Reference MPC

In order to compare and evaluate our novel data-driven MPC approaches, we
employ a reference MPC that has perfect knowledge of the battery dynamics. The
structure of the reference MPC is inspired by [17].

In what follows, we will first define a cost function. This function is then
used to formulate an MPC problem which is employed in closed-loop simulations.

###  IV-A Cost Function

At time-step $t$ and prediction instant $k$, a stage cost

$\ell(k|t)=c_{0}(p_{t}(k|t)-p_{r}(k|t))+c_{1}|\Delta\delta(k|t)|+c_{2}\delta(k|t)$
(15)

with weights $c_{1},c_{2},c_{3}\in\mathbb{R}_{>0}$ is considered. The term
$p_{t}(k|t)-p_{r}(k|t)$ rewards renewable infeed and penalizes the use of the
conventional unit. The term $c_{1}|\Delta\delta(k|t)|$ with
$|\Delta\delta(k|t)|=|\delta(k|t)-\delta(k-1|t)|$ penalizes the on and off
switching of the conventional unit. Finally, the conventional unit comes with
the power-independent cost $c_{2}\delta(k|t)$.

###  IV-B Problem Statement

Let $X=[x(k|t)]^{L}_{k=0}$ and $P=[p(k|t)]^{L-1}_{k=0}$ with
$p(k|t)=\begin{bmatrix}p_{t}(k|t)~{}p_{s}(k|t)~{}p_{r}(k|t)\end{bmatrix}^{T}$,
where $L\in\mathbb{N}$ is the prediction horizon. We can then formulate the
following optimization problem with certain inputs $w_{r}(k)$ and $w_{d}(k)$.

######  Problem 1 (Reference MPC)

$\min_{P,X}\textstyle\sum\limits_{k=0}^{L-1}\ell(k|t)\cdot\gamma^{k}$ subject
to (9)–(11), (14) for all $k\in[0,L-1]$, (12) for all $k\in[1,L]$ as well as
the dynamics $x(k+1|t)=Ax(k|t)+B_{l}p_{s}(k|t)+B_{q}p_{s}(k|t)^{2}$ (16a) for
all $k\in[0,L-1]$ with the initial conditions $x(0|t)=x^{m}(t),$ (16b)
$\delta(-1|t)=\delta^{m}(t).$ (16c)

Here, $\delta(-1|t)$ is the current on/off condition of the conventional unit.
Multiplying the predicted cost at time-step $k$ with $\gamma^{k}$,
$\gamma\in(0,1)$, allows to put less emphasis on decisions in the far future.
This can be helpful in receding horizon control. Note that Problem 1 is a
nonconvex, mixed-integer problem with quadratic constraints in the form of the
previously defined dynamics (16a) which can be solved by commercial software
such as Gurobi [22].

At each time step $t$, Problem 1 is solved. From the optimal input trajectory
$[p^{*}(k|t)]^{L-1}_{k=0}$, the first value ${p^{*}(0|t)}$ is applied to the
system and the system is operated for one time step. Then, measurements for
(16b) and (16c) are updated and Problem 1 is solved again in a receding
horizon manner.

###  IV-C Case Study

Consider the grid in Figure 1 with the parameters in Table I. Closed-loop
simulations over four weeks with a time resolution of 30 minutes (1334 total
time steps) were performed using the reference MPC. We considered hypothetical
prescient load and weather forecasts to rule out the influence of forecast
uncertainties. Moreover, the reference MPC has perfect knowledge of the
battery dynamics. For our model, we use prediction horizon $L=10$. Running the
simulation resulted in the system trajectory shown in Figure 2. The upper
figure displays the time series of each unit’s power whereas the corresponding
trajectory of the stored energy is displayed in the lower figure.

TABLE I: Parameter values of the islanded grid. Param. | Value |  | Param. | Value |  | Param. | Value
---|---|---|---|---|---|---|---
$c_{0}$ | 1 |  | $p^{max}_{t}$ | 1 |  | $B_{s}$ | -0.5
$c_{1}$ | 0.3 |  | $p^{min}_{s}$ | -1 |  | $B_{q}$ | -0.05
$c_{2}$ | 0.2 |  | $p^{max}_{s}$ | 1 |  | $x^{min}$ | 0.5
$\gamma$ | 0.9 |  | $p^{min}_{r}$ | 0 |  | $x^{max}$ | 6.5
$p^{min}_{t}$ | 0.3 |  | $A$ | 0.99 |  |  |

Figure 2: System trajectory with reference MPC. (pu: per unit, d: days)

##  V Linear Data-Driven MPC

In this section we control the system from Section III with nonlinear battery
dynamics using a linear data-driven MPC approach. We consider unknown battery
dynamics but assume mild nonlinearities, where the nonlinear terms have a much
smaller influence on the dynamics than the linear ones. This allows us to use
a modified version of (4) and build a linear data-driven MPC that accounts for
plant-model-mismatch.

###  V-A Problem Statement

Let $G$ denote our battery dynamics. We assume that $G$ is a controllable LTI
system and that $\tilde{n}$ is an upper approximation of $G$’s minimal order
$n$, i.e., $\tilde{n}\geq n$. For the moment, let us assume that we have no
plant-model-mismatch, that we can measure the battery’s stored energy, i,e.,
$y(k)=x(k)$, and control its power, $u(k)=p_{s}(k)$. Measuring the input and
its corresponding output $N$ consecutive times allows us to create the
measured sequences $\mathbb{U}^{m}$ and $\mathbb{Y}^{m}$. If $\mathbb{U}^{m}$
is persistently exciting of order $L+\tilde{n}$, then based on Theorem 1, the
I/O sequence $\\{u(k)\\}^{L-1}_{k=-\tilde{n}},\\{y(k)\\}^{L-1}_{k=-\tilde{n}}$
is a trajectory of $G$ if and only if there exists an
$\alpha\in\mathbb{R}^{N-(L+\tilde{n})+1}$ such that

$\begin{bmatrix}\mathbf{U}\\\
\mathbf{Y}\end{bmatrix}=\begin{bmatrix}H_{L+\tilde{n}}(\mathbb{U}^{m})\\\
H_{L+\tilde{n}}(\mathbb{Y}^{m})\end{bmatrix}\alpha$

with $\mathbf{U}=[u(k)]^{L-1}_{k=-\tilde{n}}$ and
$\mathbf{Y}=[y(k)]^{L-1}_{k=-\tilde{n}}$, i.e.,

$\begin{bmatrix}[p_{s}(k|t)]^{L-1}_{k=-\tilde{n}}\\\
[x(k|t)]^{L-1}_{k=-\tilde{n}}\end{bmatrix}=\begin{bmatrix}H_{L+\tilde{n}}(\\{p^{m}_{s}(t)\\}^{N-1}_{t=0})\\\
H_{L+\tilde{n}}(\\{x^{m}_{s}(t)\\}^{N-1}_{t=0})\end{bmatrix}\alpha$ (17)

for $t>N-1$.

Recall Remark 1 and Theorem 1: it is easy to see that if $\mathbb{U}^{m}$ is
persistently exciting of order $L+\tilde{n}$, then it is also persistently
exciting of order $L+n\leq L+\tilde{n}$.

The first $\tilde{n}$ elements are used to describe the initial conditions of
the battery using the equality constraint

$\begin{bmatrix}[p_{s}(k|t)]^{-1}_{k=-\tilde{n}}\\\
[x(k|t)]^{-1}_{k=-\tilde{n}}\end{bmatrix}=\begin{bmatrix}[p^{m}_{s}(k)]^{t-1}_{k=t-\tilde{n}}\\\
[x^{m}_{s}(k)]^{t-1}_{k=t-\tilde{n}}\end{bmatrix}.\\\ $ (18)

The remaining $L$ values are used to obtain a forecast for $L$ steps into the
future.

Until now, (17) and (18) where formulated assuming that $G$ is LTI. However,
$G$ is actually not LTI. In order to account for the mismatch between the
linear model I/O data and the I/O data from the nonlinear plant, we employ a
modified version of (17) in the following MPC scheme.

######  Problem 2 (Linear data-driven MPC)

$\min_{P,X,\alpha,\beta}\textstyle\sum\limits_{k=0}^{L-1}\ell(k|t)\cdot\gamma^{k}+c_{\alpha}||\alpha||^{2}_{2}+c_{\beta}||\beta||^{2}_{2}$
subject to (9)–(11), (14) for all $k\in[0,L-1]$, (12) for all $k\in[1,L]$ with
$\delta(-1|t)=\delta^{m}(t)$ as well as
$\begin{bmatrix}[p_{s}(k|t)]^{L-1}_{k=-\tilde{n}}\\\
[x(k|t)]^{L-1}_{k=-\tilde{n}}+\beta\end{bmatrix}=\begin{bmatrix}H_{L+\tilde{n}}(\\{p^{m}_{s}(t)\\}^{N-1}_{t=0})\\\
H_{L+\tilde{n}}(\\{x^{m}_{s}(t)\\}^{N-1}_{t=0})\end{bmatrix}\alpha$ (19a) for
$t>N-1$ and $\begin{bmatrix}[p_{s}(k|t)]^{-1}_{k=-\tilde{n}}\\\
[x(k|t)]^{-1}_{k=-\tilde{n}}\end{bmatrix}=\begin{bmatrix}[p^{m}_{s}(k)]^{t-1}_{k=t-\tilde{n}}\\\
[x^{m}_{s}(k)]^{t-1}_{k=t-\tilde{n}}\end{bmatrix}.\\\ $ (19b)

The measured output sequence $\\{x^{m}_{s}(k)\\}^{N-1}_{k=0}$ is composed of
one part that reflects the linear dynamics and another part which stems from
the nonlinear term. In a way, this nonlinear part can be understood as an
error that acts on the linear system. Therefore, large values of $\alpha$
increase the error that the Hankel matrices with the data from the nonlinear
system introduce. To overcome this issue, motivated by [7, 9, 11], we added
the cost term $c_{\alpha}||\alpha||^{2}_{2}$ to punish large values of
$\alpha$. Similar to [11], where state-nonlinearities are investigated, we
further added the slack variable $\beta\in\mathbb{R}^{p\cdot(L-1+\tilde{n})}$
to (19a) to account for model mismatch and used the term
$c_{\beta}||\beta||^{2}_{2}$ to keep $\beta$ small.

###### Remark 2

In [9], a similar slack variable is employed to account for noisy measurement
data and [7, 12] added a slack variable to (19b) instead of (19a). The latter
approach, did however not lead to satisfactory results in our case.

###  V-B Case Study

Let us consider the setup from Section IV-C. The input $u(k)=p_{s}(k)$ and
output $y(k)=x_{s}(k)$ of the battery were measured throughout the simulation.
We used $185$ consecutive measurements from the simulation to form the
measured sequences $\mathbb{U}^{m}=\\{u(k)\\}^{184}_{k=0}$ and
$\mathbb{Y}^{m}=\\{y(k)\\}^{184}_{k=0}$. We considered $\tilde{n}=1$ and
created the Hankel matrices $H_{L+\tilde{n}}(\mathbb{U}^{m})$ and
$H_{L+\tilde{n}}(\mathbb{Y}^{m})$. We positively validated persistence of
excitation of order $L+\tilde{n}$ for $\mathbb{U}^{m}$.

We employed Problem 2 with $c_{\alpha}=5$ and $c_{\beta}=10^{4}$ in closed-
loop simulations. The average violation of the state of charge constraints of
the reference MPC is $0$, which was expected, since it has perfect knowledge
of the battery’s dynamics. The average violation with Problem 2 was nonzero
but small. This illustrates the effectiveness of the cost terms associated
with $\alpha$ and slack variable $\beta$. However, it also shows that the
data-driven MPC did not manage to perfectly describe the actual input-
nonlinear battery dynamics. Nonetheless, one great advantage of using the
data-driven MPC is that the dynamics of the system do not need to be
explicitly known in the MPC formulation and that an identification step can be
disregarded.

We also investigated the prediction capabilities. At each step $t$, each MPC
outputs the optimal future input sequence $\\{p^{*}(k|t)\\}^{L-1}_{k=0}$.
Parallel to our actual grid (which is used for closed-loop simulations), we
use a second identical plant model to which we apply all $L$ predicted optimal
values at each step. This allows to compare the MPCs’ calculated trajectories
$[x^{*}(k|t)]^{L}_{k=1}$ with a perfectly accurate model and calculate the
prediction error of the respective MPC for each step. For the simulations
performed in Sections IV and V, the box plots of the errors can be found in
Figure 3. Despite the superiority of the reference MPC, the data-driven MPC
still manages to do a decent job in predicting the battery’s state trajectory.
However, the prediction capabilities of data-driven MPC can be further
improved, as illustrated in the next section.

Figure 3: Box plots of prediction error for reference MPC (Problem 1) and
linear data-driven MPC (Problem 2).

##  VI Hammerstein-type data-driven MPC

In this section, we formulate a data-driven MPC where only the nature of the
battery’s nonlinearity is known. First, we present the MPC problem. Then, we
use it in closed-loop simulations.

###  VI-A Problem Statement

System (13) is a Hammerstein system (see Section II-D) with $B=1$,
$\psi(u)=\alpha_{1}u+\alpha_{2}u^{2}$, $\alpha_{1}=B_{l}$ and
$\alpha_{2}=B_{q}$. To model this system in a data-driven fashion we define
the auxiliary linear input vector.

$v(k)=\begin{bmatrix}\psi_{1}(u(k))\\\ \psi_{2}(u(k))\\\
\end{bmatrix}=\begin{bmatrix}u(k)\\\ u(k)^{2}\\\ \end{bmatrix}.$ (20)

If the nature of the nonlinearity is known, then (based on Theorem 2) we can
express the system behavior by

$\begin{bmatrix}\mathbf{V}\\\
\mathbf{Y}\end{bmatrix}=\begin{bmatrix}H_{L+\tilde{n}}(\mathbb{V}^{m})\\\
H_{L+\tilde{n}}(\mathbb{Y}^{m})\end{bmatrix}\alpha$ (21)

with $\mathbf{V}=[v(k|t)]^{L-1}_{k=-\tilde{n}}$,
$\mathbf{Y}=[y(k|t)]^{L-1}_{k=-\tilde{n}}$,
$\mathbb{V}^{m}=\\{v^{m}(t)\\}^{N-1}_{t=0}$,
$v^{m}(t)=[\psi_{1}(u^{m}(t))~{}\psi_{2}(u^{m}(t))]^{T}$,
$\mathbb{Y}^{m}=\\{y^{m}(t)\\}^{N-1}_{t=0}$, for $t>N-1$ and the initial
conditions by

$\begin{bmatrix}[v(k|t)]^{-1}_{k=-\tilde{n}}\\\
[y(k|t)]^{-1}_{k=-\tilde{n}}\end{bmatrix}=\begin{bmatrix}[v^{m}(k)]^{t-1}_{k=t-\tilde{n}}\\\
[y^{m}(k)]^{t-1}_{k=t-\tilde{n}}\end{bmatrix}.$ (22)

For this approach to work, one has also to consider the relationship between
the two basis functions. In this example, this means that we need to enforce
the equality constraint

$\psi_{2}(u(k|t))=\psi_{1}(u(k|t))^{2}$ (23)

for all $k\in\begin{bmatrix}-\tilde{n},L-1\end{bmatrix}$. If
$v(k|t)=[v_{1}(k|t)~{}v_{2}(k|t)]^{T}$, then another way of formulating (23)
is by

$v_{2}(k|t)=v_{1}(k|t)^{2}.$ (24)

Fortunately, available commercial solvers, e.g., Gurobi, can handle such
quadratic constraints. In what follows, let
$\tilde{p}(k|t)=\begin{bmatrix}p_{t}(k|t)~{}v_{1}(k|t)~{}v_{2}(k|t)~{}p_{r}(k|t)\end{bmatrix}^{T}$
and $\tilde{P}=[\tilde{p}(k|t)]^{L-1}_{k=-\tilde{n}}$. We can then formulate
the following MPC problem.

######  Problem 3 (Hammerstein-type data-driven MPC)

$\min_{\tilde{P},X,\alpha}\textstyle\sum\limits_{k=0}^{L-1}\ell(k|t)\cdot\gamma^{k}$
subject to (9)–(11), (14) for all $k\in[0,L-1]$, (12) for all $k\in[1,L]$, and
the constraints $p_{t}(k|t)+v_{1}(k|t)+p_{r}(k|t)+w_{d}(k|t)=0,$ (25a)
$v_{2}(k|t)=v_{1}(k|t)^{2},$ (25b) for all $k\in[-\tilde{n},L-1]$ with
$\delta(-1|t)=\delta^{m}(t)$, as well as
$\begin{bmatrix}[v(k|t)]^{L-1}_{k=-\tilde{n}}\\\
[y(k|t)]^{L-1}_{k=-\tilde{n}}\end{bmatrix}=\begin{bmatrix}H_{L+\tilde{n}}(\mathbb{V}^{m})\\\
H_{L+\tilde{n}}(\mathbb{Y}^{m})\end{bmatrix}\alpha$ (25c)

$\begin{bmatrix}[v(k|t)]^{-1}_{k=-\tilde{n}}\\\
[y(k|t)]^{-1}_{k=-\tilde{n}}\end{bmatrix}=\begin{bmatrix}[v^{m}(k)]^{t-1}_{k=t-\tilde{n}}\\\
[y^{m}(k)]^{t-1}_{k=t-\tilde{n}}\end{bmatrix}.$ (25d)

###  VI-B Case Study

Analogously to Section V-B, we simulated the grid and used $N=185$ I/O
measurements to create the data-driven MPC. Results indicate that the
controller based on Problem 3 yields identical performance as the one based on
Problem 1. Moreover, constraints are not violated. The absolute prediction
error of both controllers is below the numerical precision of the solver (see
Figure 4) which highlights the ability of Problem 3 to accurately predict the
future state trajectory of a Hammerstein system given that the basis functions
in (20) are known. Problem 1 had perfect knowledge of the plant’s dynamics
which most of the time prerequisites a system identification step. This is not
required in Problem 3, where the measured data is directly used inside the
control loop. Furthermore, the control scheme of Problem 3 can learn changing
battery dynamics (e.g. $A$ and $B_{q}$ may change due to degradation) by
updating its Hankel matrices using more recent measurements.

Figure 4: Box plots of prediction error for reference MPC and Hammerstein-type
data-driven MPC.

##  VII Conclusions

This paper discussed how novel data-driven MPC schemes can be employed to
control power grids with nonlinear Hammerstein-type storage dynamics. First,
the unknown nonlinearity of the battery was considered to be mild, yet the
data-driven MPC managed to control the nonlinear plant with small constraint
violations. Subsequently, the structure of the nonlinearity was considered to
be known. Using an extended version of Willems’ fundamental lemma, we
developed a data-driven MPC that successfully controlled the battery with
identical results as the reference MPC with perfect knowledge of the dynamics.
Both data-driven approaches allowed to skip a system identification step and
to directly formulate MPC schemes from I/O data.

Future work will focus on theoretical considerations for the case of ”mild”
input nonlinearities and on tests with more realistic I/O data of real-world
batteries.

## References

  * [1] J. C. Willems, P. Rapisarda, I. Markovsky, and B. L. De Moor, “A note on persistency of excitation,” _Systems & Control Letters_, vol. 54, no. 4, pp. 325–329, 2005.
  * [2] J. Berberich and F. Allgöwer, “A trajectory-based framework for data-driven system analysis and control,” in _ECC_ , 2020, pp. 1365–1370.
  * [3] C. De Persis and P. Tesi, “Formulas for data-driven control: Stabilization, optimality, and robustness,” _IEEE TAC_ , vol. 65, no. 3, pp. 909–924, 2020\.
  * [4] M. Alsalti, V. G. Lopez, J. Berberich, F. Allgöwer, and M. A. Müller, “Data-driven nonlinear predictive control for feedback linearizable systems,” 2022.
  * [5] M. Alsalti, V. G. Lopez, J. Berberich, F. Allgöwer, and M. A. Múller, “Data-based control of feedback linearizable systems,” _IEEE TAC_ , pp. 1–8, 2023.
  * [6] L. Huang, J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control for grid-connected power converters,” in _IEEE 58th CDC_ , 2019, pp. 8130–8135.
  * [7] J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control: In the shallows of the deepc,” in _18th ECC_ , 2019, pp. 307–312.
  * [8] ——, “Regularized and distributionally robust data-enabled predictive control,” in _58th IEEE CDC_ , 2019, pp. 2696–2701.
  * [9] J. Berberich, J. Köhler, M. A. Müller, and F. Allgöwer, “Data-driven model predictive control with stability and robustness guarantees,” _IEEE TAC_ , vol. 66, no. 4, pp. 1702–1717, 2021.
  * [10] ——, “Linear tracking MPC for nonlinear systems—part I: The model-based case,” _IEEE TAC_ , vol. 67, no. 9, pp. 4390–4405, 2022.
  * [11] ——, “Linear tracking MPC for nonlinear systems—part II: The data-driven case,” _IEEE TAC_ , vol. 67, no. 9, pp. 4406–4421, 2022\.
  * [12] D. Bilgic, A. Koch, G. Pan, and T. Faulwasser, “Toward data-driven predictive control of multi-energy distribution systems,” _Electric Power Systems Research_ , vol. 212, p. 108311, 2022.
  * [13] Y. Xie, J. Berberich, and F. Allgöwer, “Linear data-driven economic mpc with generalized terminal constraint,” 2022.
  * [14] V. Chinde, Y. Lin, and M. J. Ellis, “Data-Enabled Predictive Control for Building HVAC Systems,” _J. Dyn. Syst. Meas. Control._ , vol. 144, no. 8, 05 2022.
  * [15] P. Mahdavipour, C. Wieland, and H. Spliethoff, “Optimal control of combined-cycle power plants: A data-enabled predictive control perspective,” _IFAC-PapersOnLine_ , vol. 55, no. 13, pp. 91–96, 2022, NECSYS 2022.
  * [16] P. Schmitz, A. Engelmann, T. Faulwasser, and K. Worthmann, “Data-driven mpc of descriptor systems: A case study for power networks,” _IFAC-PapersOnLine_ , vol. 55, no. 30, pp. 359–364, MTNS 2022.
  * [17] C. A. Hans, _Operation Control of Islanded Microgrids_ , 1st ed. Shaker Verlag, 2021.
  * [18] K. M. Chandy, S. H. Low, U. Topcu, and H. Xu, “A simple optimal power flow model with energy storage,” in _49th IEEE CDC_ , 2010, pp. 1051–1057.
  * [19] A. Hooshmand, M. H. Poursaeidi, J. Mohammadpour, H. A. Malki, and K. Grigoriads, “Stochastic model predictive control method for microgrid management,” in _IEEE PES ISGT_ , 2012, pp. 1–7.
  * [20] A. Parisio and L. Glielmo, “A mixed integer linear formulation for microgrid economic scheduling,” in _2011 IEEE SmartGridComm_ , 2011, pp. 505–510.
  * [21] M. Gholami and A. Pisano, “Model predictive operation control of islanded microgrids under nonlinear conversion losses of storage units,” _Electricity_ , vol. 3, no. 1, pp. 33–50, 2022.
  * [22] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2022. [Online]. Available: https://www.gurobi.com

  *[MPC]: model predictive control
  *[LTI]: linear time-invariant
  *[I/O]: input/output