# Learning Robust Model Predictive Control for Voltage Control of Islanded
Microgrid

Sahand Kiani, Hamed Kebriaei∗, , Mohsen Hamzeh, and Ali Salmanpour
∗Corresponding Author (Hamed Kebriaei)S. Kiani, H. Kebriaei, M. Hamzeh, A.
Salmanpour are with the School of Electrical and Computer Engineering, College
of Engineering, University of Tehran, Tehran, Iran. (emails:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper proposes a novel control design for voltage tracking of an islanded
AC microgrid in the presence of nonlinear loads and parametric uncertainties
at the primary level of control. The proposed method is based on the Tube-
Based Robust Model Predictive Control (RMPC), an online optimization-based
method which can handle the constraints and uncertainties as well. The
challenge with this method is the conservativeness imposed by designing the
tube based on the worst-case scenario of the uncertainties. This weakness is
amended in this paper by employing a combination of a learning-based Gaussian
Process (GP) regression and RMPC. The advantage of using GP is that both the
mean and variance of the loads are predicted at each iteration based on the
real data, and the resulted values of mean and the bound of confidence are
utilized to design the tube in RMPC. The theoretical results are also provided
to prove the recursive feasibility and stability of the proposed learning
based RMPC. Finally, the simulation results are carried out on both single and
multiple DG (Distributed Generation) units.

Note for Practitioners– In this paper, we present a new way to control the
voltage in an islanded microgrid to improve Power Quality (PQ). The method we
propose is based on an online optimization technique called Tube-Based Robust
Model Predictive Control. It can handle uncertainties and disturbances that
occur when the microgrid operates independently, ensuring the voltage remains
stable. However, there’s a challenge with this method. It tends to be too
cautious because it assumes the worst-case scenario for uncertainties. To make
the control more efficient, we improve it by combining a learning-based
technique called Gaussian Process regression with RMPC. The advantage of using
GP is that it predicts the uncertainty of the electrical devices based on real
data. We use these predictions to design the control in RMPC more accurately.
We also provide theoretical results to show that our new learning-based
control is reliable and stable. We tested our approach through computer
simulations on different scenarios with one or multiple power sources in the
microgrid. The results show the effectiveness of our control design in
regulating the voltage even with uncertain and nonlinear loads. Overall, this
paper suggests a practical and reliable way to control the voltage in an
independent microgrid using a combination of online optimization and learning
techniques.

###### Index Terms:

Gaussian Process Regression, Islanded Microgrid, Robust Model Predictive
Control, Voltage Control.

##  I Introduction

S everal reasons make the microgrids work in an islanded mode, such as faults
in the main grid, high prices of grid’s power, and supplying remote areas [1].
Improving local reliability, power quality, providing lower investment costs,
and reducing emissions are microgrids’ main features [2]. Voltage regulation
is a vital requirement to maintain the power delivery stable and consistent in
an islanded microgrid. Proper voltage regulation ensures reliable and
efficient operation of the microgrid, especially in remote and isolated areas
where grid connections are not available [3].

Traditional voltage regulation methods, such as PI control, fuzzy logic
control, sliding mode control, and adaptive control have been widely used in
microgrids [4][5]. PI control is a simple and reliable method that provides
good stability, but it has limited performance when dealing with nonlinear and
uncertain systems [6]. Fuzzy logic control is a more flexible method that can
handle nonlinear and uncertain systems but requires more complex modeling and
control design [7]. In recent years, researchers have focused on developing
advanced control strategies, such as MPC, to improve the performance and
stability of islanded microgrids [8][9]. MPC is an optimization-based control
method that can handle both linear and nonlinear systems with constraints,
making it more suitable for complex and uncertain systems. It has been shown
that in terms of steady-state performance, MPC outperforms sliding mode
control [10]. Explicit prediction of future plant behavior and the
computational efficiency of MPC using current systems measurements made this
online method favorable in microgrids [11][12]. In [13], the authors conducted
simulations and experimental tests to compare the performance of the proposed
MPC method with traditional control methods such as PI control and fuzzy logic
control. The results showed that the MPC-based control method outperformed
traditional methods in terms of voltage regulation performance, especially
under conditions of varying loads and disturbances.

To face with the uncertainties and disturbances existing in the model, Tube-
Based RMPC is proposed [14]. The feedback law tights the real state
trajectories inside a bound centered around the nominal system trajectories.
However, this method can be conservative since the control law is obtained
based on the worst-case scenario of uncertainties [15].

To mitigate conservativeness in Tube-Based RMPC, various regression methods
have been proposed to improve the accuracy of the tube approximation [16]. One
such method is Gaussian Process Regression (GPR). GPR is a supervised
learning-based method that has shown promise in estimating and compensating
for disturbances in MPC while reducing conservativeness and improving
computational efficiency [17]. It has also been used in power systems studies
for distribution system voltage control [18]. GPR estimates a bound of
confidence of prediction, in addition to the mean value, which is directly
used in robust controller design [19]. By using GPR to estimate disturbances
in real time, the performance of RMPC is enhanced by reducing the need for
conservative control strategies.

In this paper, we propose the Learning RMPC method that utilizes GPR to
estimate the current of the load as the disturbance, while also considering
the upper bound of parametric uncertainties in the robust controller design.
To the best of the authors’ knowledge, neither Tube-Based RMPC nor Learning
RMPC has been used to control the voltage of an islanded microgrid at the
primary level. The proposed method is initially developed for a single-DG
islanded microgrid, where nonlinear loads and parametric uncertainties can
significantly affect the system’s performance and power quality. Nevertheless,
we also demonstrate its effectiveness for multiple DGs as well. Since the
control system of an islanded microgrid is hierarchical with different layers,
in the case of multiple DG units, power sharing control loop is added.
Moreover, the recursive feasibility and stability of the proposed method have
been verified analytically. The main contributions of this paper are:

  * •

Using Tube-Based RMPC to enhance the performance of voltage regulation of an
islanded microgrid in the presence of parametric uncertainties and
disturbances (including harmonic loads) to meet PQ standards

  * •

Using Gaussian Process as a learning method to calculate the tubes online with
less conservativeness and estimate the loads to be used in shaping the control
law

  * •

Proposing recursive feasibility and stability analysis for the system, in
which the effect of GP is considered

  * •

Demonstrating the effectiveness of the proposed control method at the primary
level while multiple DGs exist in an islanded microgrid

Notation: The Euclidean norm of $x$ vector is denoted by $\|x\|$.
Respectively, The term $x^{T}Qx$ denotes $\|x\|^{2}_{Q}$. The Pontryagin
difference of sets $A\subseteq\mathbb{R}^{n}$ and $B\subseteq\mathbb{R}^{n}$
is denoted by $A\ominus B=\\{a|a+b\in A,\forall b\in B\\}$, and the Minkowski
sum is $A\oplus B=\\{a+b|a\in A,b\in B\\}$. We define the addition of set and
vector as $x\oplus B=\\{x+b|x\in\mathbb{R}^{n},b\in B\\}$. The set
multiplication is defined as if $M\in R^{m\times n}$ then $MA:=\\{Ma|a\in
A\\}$. $I_{n}$ is used to show the identity matrix of size n. The pseudo-
inverse of a $m\times n$ matrix $X$ when $n\leq m$ is
$A^{\dagger}=(A^{T}A)^{-1}A^{T}$. The set $\psi\subseteq\mathbb{R}^{n}$ is
positive invariant if and only if no solution starting inside $\psi$ can leave
$\psi$. i.e. $\forall x(0)\in\psi,\,\,\varphi(k,x(0))\in\psi\,\,\forall k>0$,
for dynamical system ${x(k+1)}=f(x(k))$ with $\varphi(k,x(0))$ as its
solution.

##  II Single-DG In An Islanded Microgrid

A single-DG in an islanded microgrid is considered. Overhauls, faults, high
energy prices, and supplying remote areas are among the reasons to switch a
microgrid from the grid-connected mode to the islanded mode. Various
disturbances might arise when the microgrid is in islanded mode, such as
current of loads. The parametric uncertainties in the inverter’s LC filter are
another source of disturbance. This paper takes into account the effect of
parametric uncertainties and loads as disturbances to the microgrid and
proposes an appropriate control method to compensate for such effects. The
microgrid is defined by three main components in islanded mode: DG, converter,
and loads.

Figure 1: Single-DG unit connected to an unknown load

As shown in Fig. 1, by neglecting the nonlinear dynamics of the DC-Bus and
using Kirchhoff’s voltage and current laws, the dynamical model can be derived
as follows

$\begin{cases}i_{f,abc}=i_{o,abc}+C_{f}\frac{dv_{abc}}{dt}\\\
u_{abc}=L_{f}\frac{di_{f,abc}}{dt}+R_{f}i_{f,abc}+v_{abc}\end{cases}$ (1)

where $i_{f,abc}$ and $i_{o,abc}$ stand for three-phase inductor current and
inverter output current. respectively, $u_{abc}$ and $v_{abc}$ show the
inverter output voltage before LC filter and output terminal voltage.
Moreover, $R_{f}$, $L_{f}$, and $C_{f}$ are the parameters of the filter. To
exploit the advantages of the d-q frame for smoother analysis, Park
transformation is used as follows

$\begin{bmatrix}\dot{V_{d}}\\\ \dot{V_{q}}\\\ \dot{I}_{fd}\\\
\dot{I}_{fq}\end{bmatrix}=\begin{bmatrix}\omega_{0}V_{q}+\frac{1}{C_{f}}I_{fd}-\frac{1}{C_{f}}I_{od}\\\
-\omega_{0}V_{d}+\frac{1}{C_{f}}I_{fq}-\frac{1}{C_{f}}I_{oq}\\\
-\frac{1}{L_{f}}V_{d}-\frac{R_{f}}{L_{f}}I_{fd}+\omega_{0}I_{fq}+\frac{1}{L_{f}}V_{od}\\\
-\frac{1}{L_{f}}V_{q}-\omega_{0}I_{fd}-\frac{R_{f}}{L_{f}}I_{fq}+\frac{1}{L_{f}}V_{oq}\end{bmatrix}$
(2)

Equation (2) forms the small-signal dynamical model of a single-DG unit in d-q
frame. To make the model compatible with MPC, the model is discretized so that
the discrete state-space model can be expressed as follows

$\displaystyle
x(k+1)=Ax(k)+Bu(k)+\underbrace{Eu_{d}(k)}_{w_{1}(k)}+\underbrace{\Delta
Ax(k)+\Delta Bu(k)}_{w_{2}(x,u,k)}$ (3) $\displaystyle y=Cx(k),$

where $x=[V_{d}\;V_{q}\;I_{fd}\;I_{fq}]^{T}\in X\subseteq R^{4}$ denotes the
state vector and $u=[V_{od}\;V_{oq}]^{T}\in U\subseteq
R^{2},u_{d}=[I_{od}\;I_{oq}]$ are control and disturbance input, respectively,
and $y=[V_{d}\;V_{q}]\in R^{2}$ is the output vector. The vectors $w_{1}(k)$
and $w_{2}(x,u,k)$ represent the effect of the current of the loads and
parametric uncertainties in the 3-phase lines, respectively. During the
operation of an inverter, the capacitance, the inductance, and the resistance
of the LC filter vary as time passes. These changes are bounded,
unpredictable, and adversely affect voltage tracking. This effect is captured
by the uncertainty matrices $\Delta A$ and $\Delta B$, obtained by
substituting $C_{f}+\Delta C_{f}$ for $C_{f}$, $L_{f}+\Delta L_{f}$ for
$L_{f}$, and $R_{f}+\Delta R_{f}$ for $R_{f}$ in (2), respectively. The
matrices will be defined as

$\displaystyle A$
$\displaystyle=\begin{bmatrix}0&\omega_{0}&\frac{1}{C_{f}}&0\\\
-\omega_{0}&0&0&\frac{1}{C_{f}}\\\
-\frac{1}{L_{f}}&0&-\frac{R_{f}}{L_{f}}&\omega_{0}\\\
0&-\frac{1}{L_{f}}&-\omega_{0}&-\frac{R_{f}}{L_{f}}\end{bmatrix},B=\begin{bmatrix}0&0\\\
0&0\\\ \frac{1}{L_{f}}&0\\\ 0&\frac{1}{L_{f}}\end{bmatrix}$ $\displaystyle E$
$\displaystyle=\begin{bmatrix}-\frac{1}{C_{f}}&0\\\ 0&-\frac{1}{C_{f}}\\\
0&0\\\ 0&0\end{bmatrix},C=\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ \end{bmatrix}$

This paper follows a novel approach to facing uncertainties and disturbances
stated in (3). $w_{1}(k)$, known as the current of the load, is the measured
disturbance. GP estimator is utilized to learn the mean and the variance of
$w_{1}(k)$. On the other hand, the upper bound for the parametric uncertainty
$w_{2}(x,u,k)$ is predetermined and known, i.e. $\|w_{2}(x,u,k)\|_{\infty}\leq
L_{2}$. It should be noted that the pattern of parametric uncertainties
related to the LC filter depends on the materials used in the inductor and
capacitance, electric flux, and temperature. As a result, in order to estimate
these values, the data should be gathered in laboratory circumstances or from
the designer. Therefore, we only use a worst-case bound for parametric
uncertainty $w_{2}(x,u,k)$. The estimated mean and variance of $w_{1}(k)$ and
the bound of $w_{2}(x,u,k)$ will be used to design Tube-Based RMPC. In what
follows, we introduce the GP estimator of load $w_{1}(k)$ and Tube-Based RMPC
as two main components of the proposed scheme.

##  III Preliminary Results

###  III-A Load’s current prediction using Gaussian Process

Loads are considered the most regular disturbances in islanded microgrids. The
current of the loads can be looked as a time series and estimated at each time
instant. Here, GPR is employed as one of the most useful Bayesian
nonparametric models. By using GPR, not only an estimate from the load current
is attained, but also we can estimate a confidence interval which will be used
as a bound in designing Tube-Based RMPC. Since we estimate this bound using
real data, it would be less conservative compared to the worst-case bound.

Consider the standard linear regression model with Gaussian noise

$\tilde{w}_{1}(k_{i})=w_{1}(k_{i})+\varepsilon(k_{i})$ (4)

where $\widetilde{\omega}_{1}(k_{i})$ denotes measured value,
$\varepsilon(k_{i})$ represents an additive noise. $\omega_{1}(k_{i})$ is the
disturbance of the DG at time instant $k_{i}$. The noise $\varepsilon\sim
N(0,\delta^{2}_{\varepsilon})$ has the property of independent and identically
distributed (i.i.d) Gaussian distribution with zero mean and variance
$\delta^{2}_{\varepsilon}$. We introduce
$\widetilde{\textbf{w}}_{1}(\tau)=[\widetilde{\omega}_{1}(k_{1}),\widetilde{\omega}_{1}(k_{2}),...,\widetilde{\omega}_{1}(k_{n})]$
as measured value vector, $\tau=[k_{1},k_{2},...,k_{n}]$ as the sample time
vector and
${\textbf{w}}_{1}(\tau)=[{\omega}_{1}(k_{1}),{\omega}_{1}(k_{2}),...,{\omega}_{1}(k_{n})]$
as load disturbance vector. The aim is to assign a distribution over
${\textbf{w}}_{1}(\tau)$ given the measured values of
$\widetilde{\textbf{w}}_{1}(\tau)$ using GP [20]. The distribution function of
the measured vector is considered to be drawn from a multivariate Gaussian
distribution as follows [21]

$p(\tilde{\textbf{w}}_{1}(\tau))=\mathcal{N}(\mu(\tau),{C}(\tau,\tau))$ (5)

where $\mu(\tau)=[m(k_{1}),\dots,m(k_{n})]$ is denoted as the vector of mean
functions $m(\tau)=\mathbb{E}[\tilde{w}(\tau)]$. ${C}(\tau,\tau)$ defined as
follows

${C}(\tau,\tau)=\begin{pmatrix}c(k_{1},k_{1})&c(k_{1},k_{2})&\dots&c(k_{1},k_{n})\\\
c(k_{2},k_{1})&c(k_{2},k_{2})&\dots&c(k_{2},k_{n})\\\
\vdots&\vdots&\vdots&\vdots\\\
c(k_{n},k_{1})&c(k_{n},k_{2})&\dots&c(k_{n},k_{n})\end{pmatrix}$ (6)

is known as the covariance matrix whose elements are covariance kernel
functions between two training inputs. The kernel function is the degree of
freedom that plays an essential role in the estimation. As it has been
discussed in [21], choosing the kernel function depends on the knowledge of
data. In this paper, the Radial Basis Function (RBF), also known as Squared
Exponential (SE) kernel, is used as the kernel function
$c(k_{i},k_{j})=h^{2}exp[-(\frac{k_{i}-k_{j}}{\lambda})^{2}]$ in the
covariance matrix of GP. The SE kernel has two hyperparameters. The length
scale $\lambda$ determines the length of the wiggle in the function, and the
output-scale $h$ determines the average distance of the function away from the
mean.

In order to generate a prediction using GPR, we need to compute the posterior
distribution function of $\bf{w_{1}}$ and then predict unseen data by
calculating the predictive posterior [22]. The posterior function consists of
two components; likelihood function and prior distribution. The likelihood
function is defined as follows

$p(\tilde{\textbf{w}}_{1}(\tau)|\textbf{w}_{1}(\tau))\sim\mathcal{N}(\tilde{\textbf{w}}_{1}(\tau)|\textbf{w}_{1}(\tau),\sigma^{2}_{N}I)$
(7)

where the mean of the function is centered on an arbitrary $\textbf{w}_{1}$.
The second component is the prior defined as follows

$p(\textbf{w}_{1}(\tau)|\tau)\sim\mathcal{N}(\textbf{w}_{1}(\tau)|\textbf{0},C(\tau,\tau))$
(8)

It is assumed the mean of the prior is always a zero vector. Nevertheless, as
we see in what follows, GP is able to model a general mean function from the
kernel-based covariance function. From (7) and (8), the posterior over
function is achieved by the likelihood function times the prior over function

$p(\textbf{w}_{1}(\tau)|\tau,\tilde{\textbf{w}}_{1}(\tau))\propto
p(\tilde{\textbf{w}}_{1}(\tau)|\textbf{w}_{1}(\tau))p(\textbf{w}_{1}(\tau)|\tau)$
(9)

Both the likelihood and the prior are Gaussian. As a result, posterior over
function will have a Gaussian distribution as follows

$\displaystyle
p(\textbf{w}_{1}(\tau)|\tau,\tilde{\textbf{w}}_{1}(\tau))\sim\mathcal{N}(\textbf{w}_{1}(\tau)|\bar{\mu},\bar{\Sigma})$
(10)
$\displaystyle\bar{\mu}={C}(\tau,\tau)[{C}(\tau,\tau)+\sigma^{2}_{N}I]^{-1}\tilde{\textbf{w}}_{1}(\tau)$
$\displaystyle\bar{\Sigma}={C}(\tau,\tau)[{C}(\tau,\tau)+\sigma^{2}_{N}I]^{-1}\sigma^{2}_{N}I$

Once the posterior over function is obtained, The predictive posterior can be
computed as follows, which is used to predict the mean and variance of unseen
data.

$\displaystyle p$
$\displaystyle(\tilde{\textbf{w}}_{1}^{\ast}(\tau)|\tau_{\ast},\tau,\tilde{\textbf{w}}_{1}(\tau))=$
(11) $\displaystyle\int
p(\tilde{\textbf{w}}_{1}^{\ast}(\tau)|\tau_{\ast},\textbf{w}_{1}(\tau),\tau)p(\textbf{w}_{1}(\tau)|\tau,\tilde{\textbf{w}}_{1}(\tau))\mathrm{d}\textbf{w}_{1}(\tau)$

where $\tilde{\textbf{w}}_{1}^{\ast}(\tau)$ denotes to a vector of new unseen
measured values with the new sample time matrix $\tau_{\ast}$. The term in the
integral, the posterior over function, is computed in (10). The training data
set is obtained from the measured value of the output current of the inverter,
which is equal to the load’s current while there is a single-DG, at $k_{1}$ to
$k_{n}$. The likelihood of new unseen data $\tilde{w}_{1}(k_{n+1})$ can be
calculated, and it is a Gaussian distribution as well.

###  III-B Tube-Based Robust Model Predictive Control

By neglecting the disturbances and uncertainties of the real system in (3),
the nominal system will be defined as follows

$\tilde{x}(k+1)=A\tilde{x}(k)+B\tilde{u}(k);\hskip
5.0pt\tilde{x}(0)=\tilde{x}_{0},$ (12)

where $\tilde{x}(k)\in\bar{X}$ and $\tilde{u}(k)\in\bar{U}$ are the nominal
trajectory and the nominal input, respectively. The sets $\bar{X}$ and
$\bar{U}$ will be defined in a way so that the robust performance is achieved.
The MPC controller works in this way: At each MPC shot $t$,
$t\in\\{0,\;1,\dots,\;T-N+1\\}$, the control law is computed for the period
$[t,t+1,\dots,t+N-1]$ by solving an optimization problem. Nevertheless, only
the first element of the computed input sequence is applied to the system, and
the state of the system evolves for one step accordingly. This process
continues for the next shots until the control inputs for the whole control
horizon are computed. The control and prediction horizon of MPC are considered
to be $T$ and $N$, respectively. In our problem, the objective is to ensure
that the voltage of the load terminal is tracked in the reference signal
$r(t)$. The reference signal is a setpoint while a single-DG unit is available
in the islanded microgrid. However, by considering multiple DG units, the
reference signal varies as time passes. The optimal nominal trajectories of
the single-DG can be obtained by solving the following nominal MPC controller.

$\displaystyle\min_{\tilde{x}_{0|0},\tilde{u}_{t}}\quad$
$\displaystyle\sum_{k=0}^{N-1}{\|\tilde{x}_{k|t}-r(t)\|^{2}_{Q}+\|\tilde{u}_{k|t}\|^{2}_{R}+\|\tilde{x}_{N|t}-r(t)\|^{2}_{P}}$
(13a) s.t. $\displaystyle\tilde{x}_{k+1|t}=A\tilde{x}_{k|t}+B\tilde{u}_{k|t},$
(13b) $\displaystyle\tilde{x}_{k|t}\in\bar{X},\tilde{u}_{k|t}\in\bar{U},$
(13c) $\displaystyle(\tilde{x}_{N|t}-r(t))\in\tilde{X}_{f}\subseteq\bar{X},$
(13d) $\displaystyle(\tilde{x}_{0|t}-x_{t})\in S_{K}(\infty),$ (13e)

where $\tilde{x}_{k|t}$ and $\tilde{u}_{k|t}$ are known as the predicted state
and control input of the $t^{\text{th}}$ MPC shot [23], which give $k$ step
ahead prediction of the state and input starting from time-step $t$, and
$\tilde{u}_{t}=[\tilde{u}_{0|t},\ldots,\tilde{u}_{N-1|t}]$. The matrices $Q$
and $R$ are semi-positive definite and positive definite. The terminal
weighting matrix $P$, is a positive definite matrix. $\tilde{X}_{f}$
represents the terminal set and is designed to ensure the stability condition
[24]. The set $S_{K}(\infty)$ will guarantee to keep the error bounded. The
calculation method of $S_{K}(\infty)$ is discussed at the end of this section.
Practically, YALMIP Toolbox is utilized [25] to solve this problem, which uses
Quadratic Programming (QP) solver. By solving the MPC problem in (13), the
optimal nominal input and the corresponding optimal state trajectory will be
found as follows

$\displaystyle\tilde{u}^{\ast}=(\tilde{u}^{\ast}(0),\ldots,\tilde{u}^{\ast}(T-1))$
$\displaystyle=[\tilde{u}^{\ast}_{0|0},\ldots,\tilde{u}^{\ast}_{T-N|T-N},\tilde{u}^{\ast}_{T-N+1|T-N},\ldots,\tilde{u}^{\ast}_{T-1|T-N}]$
$\displaystyle\tilde{x}^{\ast}=(\tilde{x}^{\ast}(0),\ldots,\tilde{x}^{\ast}(T-1),\tilde{x}^{\ast}(T))$
$\displaystyle=[\tilde{x}^{\ast}_{0|0},\ldots,\tilde{x}^{\ast}_{T-N+1|T-N},\ldots,\tilde{x}^{\ast}_{T-1|T-N},\tilde{x}^{\ast}_{T|T-N}].$
(14)

Clearly, the real trajectory of the system, $x$, will be different from the
nominal trajectory, $\tilde{x}$, due to the presence of uncertainties and
disturbances. The real control input of the system, $u$, is designed in a way
that the real trajectory of the system lies close to the nominal trajectory as
possible, using the following feedback policy

$u(k)=\tilde{u}^{\ast}(k)+K(x(k)-\tilde{x}^{\ast}(k)),$ (15)

where the feedback gain $K\in R^{2\times 4}$ can be obtained by pole
placement, LQR, and LMI technique to make the matrix $A_{K}=A+BK$ Schur stable
under the controllability of the pair (A,B). In (15), the gain $K$ compensates
deviations from nominal trajectories. By calculating the difference between
the real system and nominal system and substituting the feedback policy (15)
in (3) the error dynamic becomes

$e(k+1)=A_{K}e(k)+w(x,u,k),$ (16)

where $e:=x-\tilde{x}^{\ast}$ and $w(x,u,k)=w_{1}(k)+w_{2}(x,u,k)$ whose set
is denoted by $W$. As it was mentioned and also will be discussed in the next
subsection in more detail, this set is estimated at each iteration using the
GP estimator of $w_{1}(k)$ and the worst-case bound on $w_{2}(x,u,k)$.

The accumulative set of disturbances, $S_{K}(k)$, is defined as

$\displaystyle S_{K}(k):=\sum_{j=0}^{k-1}A_{K}^{j}W=W\oplus
AW\oplus\dots\oplus A^{k-1}W,$ (17)

As shown in [24], by setting $x(0)=\tilde{x}(0)=x_{0}$, since $A_{K}$ is
Schur, $S_{K}(\infty)$ exists and is positive invariant for the error dynamics
in (16). The feedback policy (15) assures that the state of the real system
(3) is compelled to be close to that of the nominal system (12).

As $x=\tilde{x}^{\ast}+e$, and knowing that $K$ is constant, the state of the
nominal trajectories will be the center of the tube generated by feedback
policy (15) is defined as follows

$\displaystyle x\in Z(k)\coloneqq{\tilde{x}^{\ast}(k)\oplus S_{K}(\infty)}$
(18)

Finally, to maintain the state and control input in $X$ and $U$, the state and
control input set of the nominal system need to be defined as follows

$\displaystyle\hskip 2.84526pt\bar{U}\overset{\Delta}{=}U\ominus
KS_{K}(\infty),\hskip 2.84526pt\bar{X}\overset{\Delta}{=}X\ominus
S_{K}(\infty),\hskip 2.84526pt\tilde{X}_{f}\subset X\ominus S_{K}(\infty),$
(19)

###### Proposition 1

Considering $x(k)\in Z(k)$ and
$u(k)=\tilde{u}^{\ast}(k)+K(x(k)-\tilde{x}^{\ast}(k))$, then $x(k+1)\in
Z(k+1)$ for all $w(x,u,k)\in W$ [26].

As a result, the solution of (3) using control policy (15) as its input lies
in the tube $Z(k)$ for every admissible disturbance sequence [24].

##  IV Learning Tube-Based Robust MPC for Islanded Microgrid

In this subsection, we show how the set $W$ and accordingly $S_{K}(\infty)$
can be estimated using GP and the available bound on $w_{2}$. This estimation
will then be used for tube shaping and formulating the MPC controller of the
real system.

Since the disturbance vector $w_{1}(k)$ is measured directly, it will be used
to collect the training data. We get the sample-time vector of training data
as $[k_{1},k_{2},\cdots,k_{n}]=[0,\dots,k-1]$ and the prediction sample-time
as $k_{n+1}=k$. Therefore, by estimating $\mu_{\ast}(k)$ and
$\sigma^{2}_{\ast}(k)$ using (10), the distribution of $w_{1}(k)$ is obtained
as GP $\sim\mathcal{N}(\mu_{\ast}(k),\sigma^{2}_{\ast}(k))$. By assuming
${\left\|{\mu_{\ast}(t+1)-\mu_{\ast}(t)}\right\|_{2}}\leq\Delta\mu$, the set
$\hat{W}$ can be defined as follows which represents an estimation of $W$.

$\displaystyle\hat{W}(k)=\\{$ $\displaystyle w(x,u,k)\in R^{n}:$ (20)
$\displaystyle\|\hat{w}_{1}(k)-{\mu_{\ast}(k)}\|_{{\sigma^{2}_{\ast}(k)}}\leq\chi^{2}_{n}(\vartheta)+\Delta\mu,$
$\displaystyle\|w_{2}(x,u,k)\|_{\infty}\leq L_{2}\\},$

where $\chi^{2}_{n}$ and $\vartheta$ denote the chi-distribution with $n$
degree of freedom and the confidence interval. More details are discussed in
[27]. By substituting $\hat{W}(k)$ into (17), the set $\hat{S}_{K}(\infty)$
can be defined as follows

$\displaystyle\hat{S}_{K}{(k)}:=\sum_{j=0}^{k-1}A_{K}^{j}\hat{W}=\hat{W}\oplus
A\hat{W}\oplus\dots\oplus A^{k-1}\hat{W},$ (21)

Accordingly, the new state and control input set of the nominal system will be

$\displaystyle\hskip 2.84526pt\hat{U}\overset{\Delta}{=}U\ominus
K\hat{S}_{K}(\infty),\hskip
2.84526pt\hat{X}\overset{\Delta}{=}X\ominus\hat{S}_{K}(\infty),\hskip
2.84526pt\hat{X}_{f}\subset X\ominus\hat{S}_{K}(\infty),$ (22)

In this way, we can replace constraint (13e) and (13c) in optimization problem
(13) with $(\tilde{x}_{N|t}-r(t))\in\hat{X}_{f}\subseteq\hat{X}$ and
$\tilde{x}_{k|t}\in\hat{X},\tilde{u}_{k|t}\in\hat{U}$, Then the Learning Tube-
Based RMPC optimization problem can be defined as follows:

$\displaystyle\min_{\tilde{x}_{0|0},\tilde{u}_{t}}\quad$
$\displaystyle\sum_{k=0}^{N-1}{\|\tilde{x}_{k|t}-r(t)\|^{2}_{Q}+\|\tilde{u}_{k|t}\|^{2}_{R}+\|\tilde{x}_{N|t}-r(t)\|^{2}_{P}}$
(23a) $\displaystyle\tilde{x}_{k+1|t}=A\tilde{x}_{k|t}+B\tilde{u}_{k|t},$
(23b) $\displaystyle\tilde{x}_{k|t}\in\hat{X},\tilde{u}_{k|t}\in\hat{U},$
(23c) $\displaystyle(\tilde{x}_{N|t}-r(t))\in\hat{X}_{f}\subseteq\hat{X},$
(23d) $\displaystyle(\tilde{x}_{0|t}-x_{t})\in\hat{S}_{K}(\infty),$ (23e)

where $\hat{X}_{f}$ is the invariant terminal set which is computed using the
disturbance set introduced in (20) and the method presented in [28]. The
procedure of designing Learning RMPC is illustrated in Algorithm 1.

1 Input: initial DataSet and GP, $t=t_{0}$ Measure $x_{t}$ ;

2 Update GP’s DataSet and compute ${\mu}_{*}$ and ${\sigma}_{*}$ ;

3 Update tube and constraint by using Equ. (20), (21), (22) ;

4 Set ${x}_{0|t}=x_{t}$ ;

5 Solve Problem (23) ;

6 Set $\tilde{u}_{t}=\tilde{u}_{0|t}$ ;

7 Implement $u_{t}=\tilde{u}_{t}+K(x_{t}-\tilde{x}_{0|t})$ ;

8 Set $t=t+1$ and go to step 1;

Algorithm 1 Learning Tube-Based RMPC

##  V Analytical Results

###  V-A Recursive Feasibility

An MPC is called recursively feasible if it always keeps the states in a
region from where the online optimization problem has a feasible solution
[29]. In other words, if the problem (23) has a feasible solution for the
initial condition $x_{0}$ and it remains feasible for any subsequent states
$x_{i}$ of the controlled system (3), then the MPC problem is recursively
feasible [30].

###### Lemma 1

Let’s define ${u}_{r}=B^{\dagger}(r-A{r})$. Then by having
$\tilde{u}_{N|t}=K(\tilde{x}_{N|t}-r)+u_{r}$, an invariant terminal set
$\hat{X}_{f}\ $ exists such that if $\ \tilde{x}_{N|t}-r\in\hat{X}_{f}\ $ then
$A\tilde{x}_{N|t}+B\tilde{u}_{N|t}-r\in\hat{X}_{f}$.

Proof: proof is provided in [28].

Now we can express the Recursive Feasibility theorem.

###### Theorem 1

. Assume that $\hat{X}_{f}$ is a terminal invariant set given in Lemma 1. If
the MPC problem has a feasible solution for the initial condition $x_{0}$,
then the solution of MPC is feasible for all the times.

Proof: The proof is based on mathematical induction. For this purpose, it is
assumed that the optimization problem (23) at time $t$ has a solution like
([$\tilde{x}^{\ast}_{0|t},\tilde{x}^{\ast}_{1|t},\cdots,\tilde{x}^{\ast}_{N|t}$],[$\tilde{u}^{\ast}_{0|t},\tilde{u}^{\ast}_{1|t},\cdots,\tilde{u}^{\ast}_{N-1|t}$]).

Since the tube $Z(k)$ is an invariant set for the error dynamics (16),
therefore,
$([\tilde{x}^{\ast}_{1|t},\tilde{x}^{\ast}_{2|t},\tilde{x}^{\ast}_{3|t},\cdots,\tilde{x}^{\ast}_{N|t}$]
,
[$\tilde{u}^{\ast}_{1|t},\tilde{u}^{\ast}_{2|t},\cdots,\tilde{u}^{\ast}_{N-1|t}])$
is feasible solution for MPC problem at $t+1$ from $1$ to $N-1$ prediction
horizon. Also, according to the constraint $(\tilde{x}_{N|t}-r)\in\hat{X}_{f}$
and Lemma 1, we can say that there is an input $\tilde{u}$ and final state
$\tilde{x}_{N|t+1}$ that satisfies the equation (23e). Hence, a feasible
solution for optimization (23) at time $t+1$ can be suggested as follows

$\displaystyle(\boldsymbol{x}_{t+1},\boldsymbol{u}_{t+1})=$
$\displaystyle([\tilde{x}^{\ast}_{1|t},\tilde{x}^{\ast}_{2|t},\cdots,\tilde{x}^{\ast}_{N|t},$
(24) $\displaystyle A\tilde{x}_{N|t}+B({K}(\tilde{x}_{N|t}-r)+u_{r})]$
$\displaystyle,[\tilde{u}^{\ast}_{0|t},\tilde{u}^{\ast}_{1|t},\cdots,\tilde{u}^{\ast}_{N-1|t},{K}(\tilde{x}_{N|t}-r)+u_{r}])$

Therefore, the proof is completed by mathematical induction.

###  V-B Stability

It is also important to verify the stability of the system under the proposed
Learning RMPC algorithm. The following theorem expresses the stability of the
system.

###### Theorem 2

Using control policy (15) with the nominal input derived from the optimization
problem (13) subject to the system (3), the following properties are
satisfied:

  1. 1.

$x(k)\in X$, $u(k)\in U$ for all $k\geq 0$.

  2. 2.

The closed-loop system is Input-to-State-Stable (ISS).

Proof. In the following, by using the Lyapunov method, the stability of the
system will be expressed.

$J(t)$ is defined as the cost function used in the learning tube-Based RMPC
controller as follows

$\displaystyle J(t)\
=\sum_{k=0}^{N-1}{\|\tilde{x}_{k|t}-r\|^{2}_{Q}+\|\tilde{u}_{k|t}\|^{2}_{R}+\|\tilde{x}_{N|t}-r\|^{2}_{P}}$
(25)

then, the Lyapunov function can be defined as the optimal cost function
obtained from (23e) as follows

$\displaystyle v_{t}=J_{\ast}(t)$ (26)

$v_{t}$ is clearly positive. According to the Lyapunov theorem, if the
Lyapunov function can be proved to be decreasing over the time, then the
system is stable.

$\displaystyle v_{t+1}-v_{t}=J_{\ast}(t+1)-J_{\ast}(t)$ (27)

where $J_{\ast}(t+1)$ is the optimal cost at $t+1$. Let $\hat{J}(t+1)$ be the
cost function obtained by substituting $u_{t+1}$ in (25), then:

$\displaystyle v_{t+1}-v_{t}$ $\displaystyle\leq\hat{J}(t+1)-J_{\ast}(t)$ (28)
$\displaystyle\leq\|\tilde{x}_{N|t}-r\|^{2}_{Q}+\|\tilde{u}_{N|t}\|^{2}_{R}+\|\tilde{x}_{N|t+1}-r\|^{2}_{P}$
$\displaystyle-(\|\tilde{x}_{0|t}-r\|^{2}_{Q}+\|\tilde{u}_{0|t}\|^{2}_{R}+\|\tilde{x}_{N|t}-r\|^{2}_{P})$

The $\|\tilde{x}_{N|t+1}-r\|^{2}_{P}$ can be simplified as follows

$\displaystyle\tilde{x}_{N|t+1}-r$ (29)
$\displaystyle=A\tilde{x}_{N|t}+B({K}(\tilde{x}_{N|t}-r)+u_{r})-r$
$\displaystyle=(A+B{K})(\tilde{x}_{N|t}-r_{k})+Ar+Bu_{r}-r$
$\displaystyle=(A+B{K})(\tilde{x}_{N|t}-r)$
$\displaystyle\implies\|\tilde{x}_{N|t+1}-r\|^{2}_{P}\leq\|(A+B{K})(\tilde{x}_{N|t}-r)\|^{2}_{P}$

by considering (28) and (29) :

$\displaystyle v_{t+1}-v_{t}\leq$
$\displaystyle\|\tilde{x}_{N|t}-r\|^{2}_{Q}+\|{K}(\tilde{x}_{N|t}-r)\|^{2}_{R}$
(30) $\displaystyle+$
$\displaystyle{\|u_{r}\|^{2}_{R}}+\|A_{{K}}(\tilde{x}_{N|t}-r)\|^{2}_{P}$
$\displaystyle-$
$\displaystyle(\|\tilde{x}_{0|t}-r\|^{2}_{Q}+\|\tilde{u}_{0|t}\|^{2}_{R}+\|\tilde{x}_{N|t}-r\|^{2}_{P})$
$\displaystyle\leq$
$\displaystyle\|\tilde{x}_{N|t}-r\|^{2}_{Q+{K}^{T}R{K}+(A_{{K}})^{T}P(A_{{K}})-P}$
$\displaystyle+$
$\displaystyle{\|u_{r}\|^{2}_{R}}-(\|\tilde{x}_{0|t}-r\|^{2}_{Q}+\|\tilde{u}_{0|t}\|^{2}_{R})$

${K}$ should be chosen in such a way that it applies to the discrete Lyapunov
equation $P=Q+{K}^{T}R{K}+(A_{{K}})^{T}P(A_{{K}})$, then:

$\displaystyle
v_{t+1}-v_{t}\leq{\|u_{r}\|^{2}_{R}}-(\|\tilde{x}_{0|t}-r\|^{2}_{Q}+\|\tilde{u}_{0|t}\|^{2}_{R})$
(31)

(31) is negative when $(\|\tilde{x}_{0|t}-r\|^{2}_{Q})\leq{\|u_{r}\|^{2}_{R}}$
. Therefore, the system (3) with the proposed control design within the region
of attraction $(\|\tilde{x}_{0|t}-r\|^{2}_{Q})\leq{\|u_{r}\|^{2}_{R}}$ is
stable.

##  VI Simulation Results

To evaluate the performance of the proposed Learning Tube-Based RMPC, first, a
single-DG microgrid shown in Fig. 1 is simulated to show the efficiency of
voltage control. Then, two DGs are connected in parallel and power sharing has
been done using droop control. The microgrid parameters are listed in Table I.

TABLE I: Parameters of DG Parameter | Value | Remark
---|---|---
$R_{f}$ | 1.5 m$\Omega$ |
$L_{f}$ | 100 mH |
$C_{f}$ | 100 $\mu$F |
$f_{0}$ | 60 Hz | Converter and
$S_{base}$ | 3 MVA | LC Filter
$V_{base}$ | 600 V |
$V_{dc}$ | 2000 V |
Transformer ratio ($Y$/$\Delta$) | 600/13800 V |
$T_{s}$ | 250 $\mu$s | MPC
$Z_{Line1}$, $Z_{Line2}$ | 0.35+j1.16 $\Omega$ |
$n_{1}$,$n_{2}$ | 0.5, 0.87 V/MVAr | Droop Parameters
$m_{1}$,$m_{2}$ | 0.6, 0.9 Hz/MW |

Tube-Based RMPC and Learning Tube-Based RMPC are used for the DG in the same
operational condition. A comparison between methods has been made to
illustrate the priority of the proposed method. During all the simulations,
The upper bound of parametric uncertainties is considered in parameters while
it is assumed that $\Delta R_{f}$ = 0.1$R_{f}$, $\Delta C_{f}$ = 0.1$C_{f}$
and $\Delta L_{f}$ = 0.2$L_{f}$. THD, tracking with low steady-state error and
proper dynamic response as evaluation features, has been used in this
comparison. Also, parasitic elements of switches and filter, besides the
sampling and computational delays of 202$\mu$s, have been considered. THD
permissible range is defined by the IEEE standard [31] (i.e., 5%). Simulations
have been done by the use of Matlab, Python, YALMIP, and SimPowerSystems
toolbox.

The proposed method shows how much variation exists in the worst case by
assuming a bound for disturbances and parametric uncertainties shaping the
tubes. The tubes are calculated online during the operation, which will be
advantageous for precise analysis. Shaping the control policy mentioned in
(15) will be a crucial step in Tube-Based RMPC. It can be proven that, under
the controllability of (A,B) and observability of (A,$\sqrt{Q}$), the gain $K$
obtained from above LQR ensures the optimality and stability of the system
[32]. Compared to other gains, the chosen $K$ will minimize the control effort
and give a better transient response. Using the LQR technique, matrix $K$ is
obtained as follows

$\displaystyle K$ $\displaystyle=\begin{bmatrix}-2.79\times
10^{-4}&-1.14\times 10^{-4}&0.0369&0.018\\\ -0.0028&0.0057&0.018&0.1141\\\
\end{bmatrix},$ (32)

If $K$ does not determine correctly, at least two problems will happen: First,
the controller will not be able to track the voltage precisely. Second, it is
probable that the system will become unstable. noisy observations are
considered, which can represent the measurement noise with
$\delta^{2}_{\varepsilon}$ = 0.01 variance. Moreover, a 95% confidence
interval will be shaped for the noisy data, which will be used as a bound to
shape the tubes. It is worth mentioning that the proposed method can be
implemented practically using MBE.300.E500 PMSM, commercially available by
Technosoft SA. Considering the prediction horizon of 1.25ms, and the control
horizon of 85ms, it takes 0.027455s to run Learning Robust Model Predictive
Control for the proposed model on average.

###  VI-A Harmonic Loads

Figure 2: Voltage regulation using Learning Tube-Based RMPC: voltages in d-q
frame and abc-frame while the constant impedance and harmonic loads are
connected at $t$ = 50 ms and $t$ = 70 ms.

In this section, the performance of the proposed method is evaluated. At $t$ =
50ms, a 340kVA constant impedance load with PF = 0.9 is connected to the
islanded microgrid, and a nonlinear load consisting of 5th and 7th harmonics
is connected at $t$ = 70ms. The THD of the nonlinear load’s current is
approximately 38%.

As shown in Fig. 2, voltage regulation has been done efficiently using
Learning Tube-Based RMPC. The THD of the voltage is 2.01%, which is desirable.
The tubes get tighter when the learning method is used at $t$ = 70ms.

Figure 3: Gaussian Process Regression is used to predict the harmonic current
and shape a 95% confidence interval with noisy observations.

In Fig. 3, estimating the current of the load has been shown with the 95%
confidence interval. The performance is enhanced because $\mu_{\ast}$, as an
estimation, and $\sigma_{\ast}^{2}$ is used to give the bounds besides the
parametric uncertainties upper bound, which is calculated and implemented in
tubes.

TABLE II: Comparison Among Methods for Different Harmonic Loads |  Voltage THDs $\%$
---|---
Harmonics |  Current THDs $\%$ | LRMPC | RMPC | MPC
5th, 7th, 11th | 30.83 | 1.97 | 2.02 | 3.58
41.69 | 2.03 | 2.07 | 3.67
51.39 | 2.07 | 2.15 | 3.91
5th, 11th, 13th | 31.49 | 2.02 | 2.11 | 3.62
44.62 | 2.12 | 2.24 | 3.76
52.71 | 2.21 | 2.36 | 4.02

In Table II, various harmonic loads have been connected to the islanded
microgrid, and regulating the voltage has been done. Voltage THDs have been
calculated and reported. As the table illustrates, voltage THDs are desirable
for the proposed method. It is worth mentioning that the conventional
controller, PI, have been implemented to handle the above scenario while the
harmonic load with $5^{th}$, $7^{th}$, and $11^{th}$ and the current THD of
38% is connected. The results have been shown in the following table

TABLE III: Comparison Among Methods Based on Voltage THDs - Harmonic Loads Learning RMPC | Tube-Based RMPC | MPC | Conventional PI
---|---|---|---
2.01% | 2.08% | 3.62% | 3.74

Moreover, the relationship between the learning results and the control input
have been illustrated in the following figure

Figure 4: Voltage regulation using Learning Tube-Based RMPC: learning results
and control input in d-q frame and abc-frame while the constant impedance and
harmonic loads are connected at $t$ = 50 ms and $t$ = 65 ms.

In Fig. 4, we have considered a constant impedance load besides the harmonic
loads with $5^{th}$, $11^{th}$, and $13^{th}$ and the current THD of 44.62%.

###  VI-B Constant Power Loads

In this section, a 340kVA constant impedance load with PF = 0.9 is connected
to the islanded microgrid at $t$ = 50ms, and a 500kVA constant power load with
PF = 1 parallel with a harmonic load consist of $5^{th}$ and $7^{th}$
harmonics is connected to the output terminal at $t$ = 70ms. Learning RMPC is
used to handle the voltage regulation.

Figure 5: Voltage regulation using Learning Tube-Based RMPC: voltages in d-q
frame and abc-frame while the constant impedance and the constant power load
parallel with harmonic load are connected at $t$ = 50 ms and $t$ = 70 ms.

In Fig 5, Learning RMPC will track the voltage with desirable voltage THD of
2.39% which is proper. A comparison has been made in the following table among
different methods as follows

TABLE IV: Comparison Among Methods Based on Voltage THDs - Constant Power Load Learning RMPC | Tube-Based RMPC | Conventional PI
---|---|---
2.39% | 2.48% | 3.95%
Figure 6: Gaussian Process Regression is used to predict the current of the
loads and shape a 95% confidence interval with noisy observations.

In Fig. 6, GP is used to predict the current of the loads. The variance has
been used to shape the tube with less conservativeness.

###  VI-C Droop Control for Power Sharing

Learning RMPC is proposed at the primary level of control for a single-DG
unit. The proposed method can be also used in a microgrid with multiple DG
units. two DGs are considered for adding the power sharing layer while they
are working in parallel using droop control method.

The instantaneous voltage and the current of the output terminal will be used
to calculate the active power and reactive power, respectively. Then, the
powers will be applied to a low pass filter with a cutoff frequency of 10Hz to
remove the ripples. $P_{DG}$ and $Q_{DG}$ as the average powers will be the
input of droop controller to produce the reference signals as follows

$\displaystyle f_{i}=f_{n_{i}}-m_{i}P_{DG_{i}},$ (33) $\displaystyle
v^{*}_{d_{i}}=V_{n_{i}}-n_{i}Q_{DG_{i}},$ $\displaystyle v^{*}_{q_{i}}=0,$

where $m_{P_{i}}$ and $n_{Q_{i}}$ are the frequency and voltage droop gains,
respectively. $V_{n_{i}}$ and $f_{n_{i}}$ are the desired frequency and
voltage for i-th DG unit, respectively.

Droop control is a vital aspect of parallel inverter operation, while phase-
lock loops (PLL) are commonly used for synchronization. Notably, droop control
possesses an inherent synchronization mechanism, obviating the need for a
dedicated synchronization unit. This allows droop control to effectively
achieve the same function as a PLL, as detailed in the study by Zhong and
Zhang (2020) [33].

Figure 7: Instantaneous power and reactive power of two DG units while a
constant impedance and a constant power load are connected.

Fig. 7 shows the instantaneous active and reactive powers while a 340kVA
constant impedance with PF=0.9 and a 500kVA constant load with PF=1 are
connected at $t$ = 350ms and $t$=700ms, respectively. As observed, the LPF
with $f_{L}$ = 10Hz removes the ripples, and the average of powers is shared
among DG units among them according to their droop charactristics.

Figure 8: Voltage and frequency of the load terminal.

Fig. 8 shows the voltage and frequency of the load terminal. The references
from droop controller will be used in Learning RMPC to force to track them. As
the dynamics are slow, the performance of LRMPC is not affected significantly.

##  VII Conclusions

This paper proposed a novel Learning Tube-Based RMPC to track the voltage of
an AC inverter-based islanded microgrid at the primary level of control. This
technique could handle nonlinear loads as disturbances—besides the upper bound
of parametric uncertainties. As a powerful learning method, GPR assisted in
designing the Tube-Based RMPC with less conservativeness by estimating the
load and shaping the tubes online using available data. The advantages and
disadvantages of the proposed method are mentioned as follows: (i): By using
Tube-Based RMPC, the performance of voltage regulation of an islanded
microgrid is enhanced in the presence of parametric uncertainties and load
disturbances to meet PQ standards. (ii) The performance of Tube-based RMPC is
enhanced by incorporating a Gaussian Process method to estimate the mean value
and the bound of uncertainties based on real data, which leads to less
conservativeness compared to the worst-case design. (iii): Harmonic loads as a
disturbance are handled using Learning Tube-Based MPC. (iv): The computational
cost of Learning Tube-Based MPC is increased compared to benchmarking methods
like PI, MPC, and also Tube-Based MPCs.

Moreover, we showed that our method can be used in an islanded microgrid with
multiple DG units, while droop control was used for power sharing among DGs.
Besides, the recursive feasibility and stability analysis of Learning Tube-
Based RMPC has been verified analytically. Finally, MPC, Tube-Based RMPC, and
Learning Tube-Based RMPC were compared with each other in terms of THD,
performance, and steady-state error in different conditions, which verified
the effectiveness of the proposed Learning Tube-Based RMPC.

## References

  * [1] N. M. Dehkordi, N. Sadati, and M. Hamzeh, “Fully distributed cooperative secondary frequency and voltage control of islanded microgrids,” _IEEE Transactions on Energy Conversion_ , vol. 32, no. 2, pp. 675–685, 2016.
  * [2] A. Singh and B. S. Surjan, “Microgrid: A review,” _International Journal of Research in Engineering and Technology_ , vol. 3, no. 2, pp. 185–198, 2014\.
  * [3] A. Bidram, A. Davoudi, F. L. Lewis, and S. S. Ge, “Distributed adaptive voltage control of inverter-based microgrids,” _IEEE Transactions on Energy Conversion_ , vol. 29, no. 4, pp. 862–872, 2014.
  * [4] X. Shen, J. Liu, H. Lin, Y. Yin, A. M. Alcaide, and J. I. Leon, “Cascade control of grid-connected npc converters via sliding mode technique,” _IEEE Transactions on Energy Conversion_ , 2023.
  * [5] V. Thomas, S. Kumaravel, and S. Ashok, “Fuzzy controller-based self-adaptive virtual synchronous machine for microgrid application,” _IEEE Transactions on Energy Conversion_ , vol. 36, no. 3, pp. 2427–2437, 2021.
  * [6] N. R. Merritt, C. Chakraborty, and P. Bajpai, “New voltage control strategies for vsc-based dg units in an unbalanced microgrid,” _IEEE Transactions on Sustainable Energy_ , vol. 8, no. 3, pp. 1127–1139, 2017.
  * [7] M. Huang, Y. Kang, and F. Liu, “A model predictive control strategy for voltage regulation in islanded microgrids,” _IEEE Transactions on Industrial Informatics_ , vol. 14, no. 7, pp. 3013–3024, 2018.
  * [8] M. Davari, W. Gao, Z.-P. Jiang, and F. L. Lewis, “An optimal primary frequency control based on adaptive dynamic programming for islanded modernized microgrids,” _IEEE Transactions on Automation Science and Engineering_ , vol. 18, no. 3, pp. 1109–1121, 2021.
  * [9] A. La Bella, S. Raimondi Cominesi, C. Sandroni, and R. Scattolini, “Hierarchical predictive control of microgrids in islanded operation,” _IEEE Transactions on Automation Science and Engineering_ , vol. 14, no. 2, pp. 536–546, 2017.
  * [10] A. V. Doostdar, A. Khodabakhshian, and S. H. Hosseini, “Comparing sliding mode and model predictive control for a three-phase four-leg inverter in islanded microgrid applications,” _IEEE Transactions on Industrial Electronics_ , vol. 64, no. 10, pp. 8112–8122, 2017.
  * [11] A. Ouammi, Y. Achour, D. Zejli, and H. Dagdougui, “Supervisory model predictive control for optimal energy management of networked smart greenhouses integrated microgrid,” _IEEE Transactions on Automation Science and Engineering_ , vol. 17, no. 1, pp. 117–128, 2020.
  * [12] N. Sockeel, J. Gafford, B. Papari, and M. Mazzola, “Virtual inertia emulator-based model predictive control for grid frequency regulation considering high penetration of inverter-based energy storage system,” _IEEE Transactions on Sustainable Energy_ , vol. 11, no. 4, pp. 2932–2939, 2020.
  * [13] W. Zhang, Y. Hu, H. He, Y. Wu, and M. Li, “Model predictive control for voltage regulation in islanded microgrid based on transformer-less cascaded h-bridge inverter,” _IEEE Transactions on Industrial Informatics_ , vol. 15, no. 3, pp. 1476–1486, 2019.
  * [14] A. Alessio and S. Di Cairano, “Tube MPC for linear systems with bounded disturbances,” _IEEE Transactions on Automatic Control_ , vol. 60, no. 3, pp. 805–810, 2015.
  * [15] A. Oshnoei, S. Peyghami, F. Blaabjerg _et al._ , “Robust control of voltage source converters: A tube-based model predictive approach,” _IEEE Transactions on Circuits and Systems II: Express Briefs_ , 2023.
  * [16] J. Paulson and J. Luedtke, “Tube-based robust model predictive control via sample-based gaussian process regression,” _Automatica_ , vol. 47, no. 8, pp. 1663–1671, 2011.
  * [17] P. Zhang, B. Li, X. Liu, and Q. Liu, “Data-driven model predictive control using gaussian process regression for industrial processes,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 8, pp. 4957–4966, 2020.
  * [18] P. Li, J. Ji, H. Ji, J. Jian, F. Ding, J. Wu, and C. Wang, “Mpc-based local voltage control strategy of dgs in active distribution networks,” _IEEE Transactions on Sustainable Energy_ , vol. 11, no. 4, pp. 2911–2921, 2020.
  * [19] L. Hewing, K. P. Wabersich, M. Menner, and M. N. Zeilinger, “Learning-based model predictive control: Toward safe learning in control,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 3, pp. 269–296, 2020\.
  * [20] C. E. Rasmussen and C. K. I. Williams, _Gaussian processes for machine learning._ , ser. Adaptive computation and machine learning. MIT Press, 2006.
  * [21] S. Roberts, M. Osborne, M. Ebden, S. Reece, N. Gibson, and S. Aigrain, “Gaussian processes for time-series modelling,” _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , vol. 371, no. 1984, p. 20110550, 2013.
  * [22] B. A. Swastanto, “Gaussian process regression for long-term time series forecasting,” 2016.
  * [23] D. Limón, I. Alvarado, T. Alamo, and E. F. Camacho, “Robust tube-based mpc for tracking of constrained linear systems with additive disturbances,” _Journal of Process Control_ , vol. 20, no. 3, pp. 248–260, 2010.
  * [24] J. B. Rawlings, D. Q. Mayne, and M. Diehl, _Model predictive control: theory, computation, and design_. Nob Hill Publishing Madison, WI, 2017, vol. 2.
  * [25] J. Lofberg, “Yalmip: A toolbox for modeling and optimization in matlab,” in _2004 IEEE international conference on robotics and automation (IEEE Cat. No. 04CH37508)_. IEEE, 2004, pp. 284–289.
  * [26] D. Q. Mayne, M. M. Seron, and S. Raković, “Robust model predictive control of constrained linear systems with bounded disturbances,” _Automatica_ , vol. 41, no. 2, pp. 219–224, 2005.
  * [27] A. D. Bonzanini, D. B. Graves, and A. Mesbah, “Learning-based smpc for reference tracking under state-dependent uncertainty: An application to atmospheric pressure plasma jets for plasma medicine,” _IEEE Transactions on Control Systems Technology_ , vol. 30, no. 2, pp. 611–624, 2021\.
  * [28] L. Fagiano and A. R. Teel, “Model predictive control with generalized terminal state constraint,” _IFAC Proceedings Volumes_ , vol. 45, no. 17, pp. 299–304, 2012.
  * [29] L. Yang, A. Katnik, and N. Ozay, “Quickly finding recursively feasible solutions for mpc with discrete variables,” in _2019 IEEE Conference on Control Technology and Applications (CCTA)_ , 2019, pp. 374–381.
  * [30] J. Rawlings, D. Mayne, and M. Diehl, _Model Predictive Control: Theory, Computation, and Design_ , 01 2017.
  * [31] “Ieee recommended practice and requirements for harmonic control in electric power systems,” _IEEE Std 519-2014 (Revision of IEEE Std 519-1992)_ , pp. 1–29, 2014.
  * [32] D. P. Bertsekas, _Dynamic Programming and Optimal Control_ , 3rd ed. Belmont, MA, USA: Athena Scientific, 2005, vol. I.
  * [33] Q.-C. Zhong, “Synchronization mechanism of droop control,” 2020.


# An Evolving Population Approach to Data-Stream Classification with Extreme
Verification Latency

Conor Fahy and Shengxiang Yang  Institute of Artificial Intelligence, School
of Computer Science and Informatics
De Montfort University, Leicester LE1 9BH, UK
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Recognising and reacting to change in non-stationary data-streams is a
challenging task. The majority of research in this area assumes that the true
class label of incoming points are available, either at each time step or
intermittently with some latency. In the worse case this latency approaches
infinity and we can assume that no labels are available beyond the initial
training set. When change is expected and no further training labels are
provided the challenge of maintaining a high classification accuracy is very
great. The challenge is to propagate the original training information through
several timesteps, possibly indefinitely, while adapting to underlying change
in the data-stream. In this paper we conduct an initial study into the
effectiveness of using an evolving, population-based approach as the mechanism
for adapting to change. An ensemble of one-class-classifiers is maintained for
each class. Each classifier is considered as an agent in the sub-population
and is subject to selection pressure to find interesting areas of the feature
space. This selection pressure forces the ensemble to adapt to the underlying
change in the data-stream.

###### Index Terms:

Data Stream Analysis, Scarcity of Labels , Extreme Verification Latency ,
Evolutionary Algorithms

##  I Introduction

A data-stream is a continuously arriving sequence of data. Let
$S=[(x_{i},y_{i})^{t}]_{t=0}^{\infty}$ represent a stream where $x_{i}$ is a
point which describes a class $y_{i}$, at time $t$. Often the underlying
operation which generates $S$ is non-stationary and distributions within $S$
will drift as $S$ progresses. Classification models which are static and not
adaptive to these changes are likely to degrade over time and their
predictions will become increasingly less reliable. In the worst case, change
in the stream can be so severe that our non-adaptive model becomes entirely
obsolete.

This change can happen in a number of ways. Concept drift occurs when the
characteristics of the incoming data points ($x_{i}$) change, or the
relationship between the data ($x_{i}$) and the target class ($y_{i}$)
changes. Formally, we describe this relationship as a conditional probability
$P^{t}(y_{i}|x_{i})$ where vector $x_{i}$ describes class $y_{i}$ at time $t$
. Concept drift can be described as being real or virtual. In virtual concept
drift there is a change in $P(x_{i})$ which does not result in a change in
$P(y_{i}|x_{i})$. Here, we could say that the underlying distribution, or
features of $x_{i}$ have changed but that $x_{i}$ still describes a class
$y_{i}$. When real drift occurs $P(y_{i}|x_{i})$ changes. This change might
happen gradually, maybe as virtual drift initially, or the change could be
sudden without any change in $P(x_{i})$. When real drift occurs,
$P^{t}(y_{i}|x_{i})\neq P^{t+\delta}(y_{i}|x_{i})$ after some time $\delta$.

In a fully supervised environment, the true class label $y_{i}$ of $x_{i}$, is
available and there is no _scarcity of labels_ , however classification in
non-stationary settings when there is a scarcity of labels represents a
greater challenge but is often a practical reality when processing data
streams. Here, the true class labels of stream samples are not available
immediately and there is a verification latency.

Latency $L$ is in the range $[0$, $\infty]$ and three categories of latency
are outlined in [20]: null latency when $L=0$. Intermediate latency where
$0<L<\infty$, and extreme latency where $L\to\infty$, which is the worst case
where the true label is never available.

This study focuses on the third, and worst, case of verification latency. We
assume that some label information is provided in an initialisation phase but
subsequently never again. Furthermore, some form of change is expected and the
classification model must adapt to this change without access to labels. Such
a setting is referred to as Initially Labelled Non-stationary Streaming (ILNS)
[4] or Extreme Verification Latency (EVL) [5]. For the remainder of this paper
we adopt the latter, $EVL$.

This study evaluates population-based approaches for updating and maintaining
the classification model. In the proposed method the classification model is a
set of one-class-classifier ensembles. Each ensemble is considered to be a
sub-population and each individual classifier within the ensemble is an agent
within the population. Each agent is subject to selection pressure which
forces the population to evolve and adapt to the changing data distribution.
We experiment with two bio-inspired approaches, namely: Genetic Algorithm
(GA), and Particle Swarm Optimisation (PSO).

The remainder of this paper is outlined as follows: A review of the existing
EVL literature is provided in Section II. Our proposed method is outlined in
Section III. The experimental section is outlined in Section IV. This includes
the data used, metrics, results, peer-comparisons, and a complexity analysis.
Finally, Section V provides conclusions and future work.

##  II Related Work

An approach to address the problem of EVL involves representing data as a
combination of parametric distributions. This concept was initially introduced
in the Arbitrary Sub Populations Tracker (APT) [16]. APT operates based on
three assumptions: 1) the drift is gradual and systematic, 2) all sub-
populations exist in the training data, and no new sub-populations emerge
later in the data stream, and 3) the drift remains constant. However, APT is
computationally expensive and relies on assumptions that may not always be
practical. The assumption of constant drift affecting each class identically
is often unrealistic. Nevertheless, the assumption of gradual drift is
typically necessary and prevalent in most EVL algorithms.

When assuming gradual drift, a drifting distribution will exhibit overlap
between successive windows. In such cases, the central region of the
distribution, known as the core, will have the most substantial overlap with
the distribution in subsequent time steps. Core Support Extraction (CSE) is a
popular strategy for handling EVL. CSE involves identifying instances within
the core region of each distribution at time $t$, referred to as Core samples
(CS), while discarding non-core samples.

COMpacted POlytope Sample Extraction (COMPOSE) [5] is an example of an EVL
approach employing CSE. COMPOSE employs the $\alpha$-shape [6], a
generalization of the convex hull, to identify the core support region for
each class. However, the computational complexity associated with calculating
the $\alpha$-shape makes COMPOSE impractical for high-dimensional data
streams. To address this, Learning Extreme Verification Latency with
Importance Weighting (LEVELIW) [21] proposed an importance weighting scheme as
a substitute for the $\alpha$-shape. Both methods achieve similar results, but
LEVELIW offers improved speed. Another proposal to replace the $\alpha$-shape
is discussed in [3].

FAST COMPOSE [22] is a novel method that eliminates Core Support Extraction
entirely. The idea behind FAST COMPOSE is to propagate all labeled points from
one window to the next and transfer the labels from previously labeled points
to unknown points in the latest window.

AMANDA [11] presents another CSE approach. After a training phase, incoming
unlabeled data is classified using a Semi-Supervised Learning (SSL)
classifier. Samples in each class are then evaluated using Kernel Density
Estimation (KDE) and subsequently filtered based on their KDE scores. Only the
densest samples from each class are propagated to the next time step.

SLAYER (Semi-Supervised Stream Learning with density based drift) [2] proposes
an alternative to CSE by utilizing clustering. A similar approach is outlined
in Stream Classification Algorithm Guided by Clustering [18]

A micro-cluster is an $n$-dimensional sphere characterized by its center, $c$,
and radius, $r$. It was initially proposed as a density clustering mechanism
[1] and has been widely employed in the stream clustering literature [13, 7,
12]. The supervised counterpart of a micro-cluster is called a micro-
classifier (MClassification), which also records a class label [19]. The label
of an incoming point $p$ is predicted by its nearest micro-cluster $m\phi_{i}$
and MClassification updates $m\phi_{i}$ with $p$. Bio-Inspired approaches to
data stream analysis have been explored in FlockStream [12] (inspired by the
flocking behaviour of birds) and there have been approaches based on the
observed behaviour of ants [9, 10].

A survey on EVL methods is given in [23] and an overview of classification
with scarcity of labels, which includes but not limited to the EVL setting, is
given in [8].

##  III Proposed Method

An incoming Stream $S$ is processed in tumbling windows of size $w$ (in a
tumbling window, a fixed size, non-overlapping chunk of data is considered).
Each incoming point $p_{i}$ in $w$ is passed to the set of ensembles
($e_{1},\dots,e_{n}\in E$) and a prediction $\hat{y}$ is returned. At the end
of each window, the fitness of each base-classifier is calculated and the
population of classifiers is updated to adapt to the underlying stream.

###  III-A Classification

The proposed method consists of a set of ensembles $E=\\{e_{1},\dots,e_{n}\\}$
where $n$ is the number of classes present in the training data and
$e_{i}=\\{\phi_{1},\dots,\phi_{j}\\}$ where $j$ is a user supplied parameter
which determines the ensemble size and $\phi$ is a base-classifier. In the
following, we use micro-classifiers as the base one-class-classifier (though
any distance based or density based one-class classifier could be used). A
micro-classifier is an $n$-dimensional sphere, represented as $[{c,r,y,f}]$,
where $c$ is the centre of the micro-classifier in $D$ dimensions, $r$ is the
radius of the micro-classifier, $y$ is the associated label, and $f$ is the
classifier’s current fitness. A point $p$ is recognised by a classifier if it
falls inside the classifier’s boundary; otherwise, the point is unrecognised
by this component of the ensemble. We estimate the confidence of the
prediction by using the distance from the centre $c$ to point $p$, formally:

$\begin{split}dist({\phi_{i}}_{c},p)\leq{\phi_{i}}_{r}=\\\
\begin{cases}1:prediction~{}{\phi_{i}}_{y},&confidence:dist({\phi_{i}}_{c},p)\\\
0:unrecognised,&confidence:dist({\phi_{i}}_{c},p)\end{cases}\end{split}$ (1)

Point $p$ is passed to each member of the ensemble and the most confident
classifier’s prediction is returned. If $p$ falls outside the recognised area
of the feature space, then the point is classified as ‘unrecognised’. In these
instances, $p_{i}$ might represent noise, outlier or could be an indication of
change in the stream. In some cases, it might be more appropriate to make no
prediction if a point is unrecognised, in other instances it might be better
to predict every point. In the latter, the ensemble can return the prediction
of the classifier most similar to $p_{i}$. This process is outlined in
Algorithm 1.

Algorithm 1 Classification

$\boldmath{Input:}$ Ensemble E, point $p$
$\boldmath{Output:}$ Prediction $\hat{y}$

1:classifiersWhichRecognisePoint $\leftarrow\\{\\}$

2:LabelsSet $\leftarrow\\{\\}$

3:for $e_{i}$ in $E$ do

4: for $\phi_{i}$ in $e_{i}$ do // for each classifier in ensemble

5: $d\leftarrow$ distance($\phi_{i_{centre}},p$)

6: if $(d\leq\phi_{i_{rad}})$ then // if point recognised

7: $\phi_{i_{dist}}\leftarrow d$

8: $\phi_{i_{fitness}}\leftarrow\phi_{i_{fitness}}+1$

9: classifiersWhichRecognisePoint $\leftarrow$ $\phi_{i}$

10: if $\phi_{y}$ not in LabelsSet then

11: $LabelsSet$ $\leftarrow$ $\phi_{y}$

12:$C\leftarrow\arg\min_{dist}$(classifiersWhichRecognisePoint)

13:$\hat{y}\leftarrow C_{y}$

14:if ($|LabelsSet|>1$) then // there is conflict

15: for ($\phi_{i}$ in classifiersWhichRecognisePoint) do

16: $\phi_{i_{fitness}}\leftarrow 0$

17:return $\hat{y}$

###  III-B Initialising the Ensemble

Given a training set containing $n$ classes where $y_{1},\dots,y_{n}\in Y$,
the trained classifier $E$ will be composed of $n$ ensembles, each ensemble
$e_{i}$ is a one class classifier specialised in one area of the feature
space. To initialise $E$ two user supplied variables are required; the number
of classifiers in each ensemble ($numAgents$) and the radius of each base
classifier (the sensitivity of these values will be discussed in Section
IV-E). First, the training data is split by class and an ensemble is created
for each class. Training samples in each class $y_{i}\in Y$ are added
sequentially. The first instance $x_{i}$ creates the first micro-classifier
with c = $x_{i}$, r = radius, y = $y_{i}$, and $f$ = 0. Subsequent points
$x_{i}$ are passed to existing micro-classifiers to check if they are
recognised. If so, the training point is ignored. However, if the point is
unrecognised by any of the existing classifiers, it forms a new classifier.
This process continues while the number of classifiers is less than or equal
to $numAgents$. If this limit is reached before all training samples have been
considered, it is typically a sign that the value for $\phi_{r}$ is too small.
Conversely, if training is complete and $|e_{i}|<numAgents$, we generate
additional $k$ classifiers ($k=numAgents-|e_{i}|$) based on the distribution
of existing classifiers. We sample (without replacement) an existing
classifier’s centre $\phi_{rand{{}_{c}}}$ and alter the value of a random
dimension using a Gaussian distribution, a new classifier is created from this
point and added to $e_{i}$.

###  III-C Calculating the fitness function

If we had access to the ground truth, we could compare each classifier’s
prediction with the correct label and use this measure as the classifier’s
fitness. However, in the EVL setting, we do not have access to the underlying
truth.

In the proposed method an assumption is made that change in the stream is
gradual. $P(x)$ will drift, potentially leading to a change in $P(y|x)$, but
$P(y|x)$ will not change without a change in $P(x)$. In the gradual-change
assumption, interesting areas of the feature space in one window will have the
most overlap with interesting areas in subsequent windows. Once the ensemble
has been initialised with labelled data, these interesting areas have been
identified and subsequently a $local$ exploration and adaption to the gradual
underlying change is required.

Classifiers which recognise incoming points are deemed fit and classifiers
which are no longer recognising points (or fewer points) are deemed unfit.
This will happen when the underlying distribution has drifted from this part
of the feature space. When a classifier makes a prediction (even if this
prediction is not used by the ensemble), its fitness increases and after a
window the fitness of each classifier decreases linearly
($\phi_{f}(t+1)=\phi_{f}(t)-1$)), so classifiers which we once deemed to be
fit gradually become less fit and have fewer opportunities to reproduce or
continue to subsequent populations. Alternately, classifiers which are
recognising a lot of points are interpreted as being in an interesting part of
the feature space and so have a greater chance to reproduce and explore their
locality. To prevent overlap and situations where classifiers $A$ and $B$ with
associated labels $y_{a}$ and $y_{b}$ both recognise a point and return
contradicting predictions, we penalise the classifiers which make conflicting
predictions. These classifiers have a fitness of zero which promotes
separation.

We evaluate two population-based approaches to update the ensemble: A swarm-
based approach in Particle Swarm Optimsation (PSO), and an evolutionary
approach with a Genetic Algorithm (GA). PSO is an optimisation technique which
draws inspiration from the flocking behaviour observed in nature. Each
classifier is considered a ‘particle’ and these particles move through the
feature space following the fittest members of the swarm and biasing their
movement towards historically interesting areas of the space. The position of
a particle is influenced by its own best past location and the best past
location of the swarm. In the proposed method, the number of swarms is equal
to the number of classes in the initialisation phase, so each class is
described by a swarm consisting of $n$ particles, where each particle is a
classifier.

The GA is inspired by population genetics and belongs to the family of
evolutionary computation methods. In the proposed method, each class in the
initialisation phase is considered a separate sub-population consisting of $n$
agents. The fittest members of each sub-population are selected (proportional
to their fitness) as ‘parents’. In this phase only the centre of the
classifier is considered (not the radius or label). Selected parent solutions
are recombined, mutated and the ‘children’ solutions are passed to the next
generation. This recombination and mutation constitute a local search and
solutions from sub-population $A$ can not combine with solutions from sub-
population $B$. The aim is that fitter classifiers will have more of an
opportunity to explore their

TABLE I: Training Data Used in Experiments Dataset | Classes | Features | Samples
---|---|---|---
$1CDT$ | 2 | 2 | 16,000
$4CR$ | 4 | 2 | 144,000
$1CSurr$ | 2 | 2 | 55,283
$GEARS$ | 2 | 2 | 200,000
$UG-2C-2D$ | 2 | 2 | 100,000
$UG-2C-3D$ | 2 | 3 | 200,000
$UG-2C-2D$ | 2 | 5 | 200,000
$4CRE-V1$ | 4 | 2 | 125,000
$4CRE-V2$ | 4 | 2 | 183,000
$KeyStroke$ | 4 | 10 | 1,600

local space and can adapt to gradual change in this way.

Parent solutions are selected using 2-way tournament selection. Two candidates
are randomly selected from the population and a copy of the fittest is passed
to a pool of parent solutions which are then combined with fixed-point
crossover.

##  IV Experimental Study

The following section outlines the experimental data-streams used to validate
the proposed method, along with the metrics used. A comparison of both bio-
inspired population based approaches is presented along with a comparison with
peer algorithms. A complexity analysis is offered in Section IV-E.

###  IV-A Data Streams

Ten dynamic data streams are used to validate the proposed methods. These are
taken from the non-stationary archive
page111https://sites.google.com/site/nonstationaryarchive/ and have been used
in previous studies on non-stationary data streams [7, 18] and also in the EVL
setting [11, 2].

The characteristics of each stream is presented in Table I (and a
visualisation of the streams can be seen from the non-stationary archive page,
footnote 1) . The number of classes ranges from 2 to 4, and the dimensionality
ranges from 2 to 10. Although some of the streams appear simple, they are non
trivial. For example, synthetic stream 4CR consists of four rotating classes.
Each class gradually moves into an area of the feature space previously
occupied by a different class. So, along with adaptation, some form of
‘forgetting’ is required. Most are synthetic streams and a real stream is
included in Keystroke data. This stream is based on keystroke dynamics, users
are recognised by their typing rhythm which is likely to change over time as
concept drift. The stream describes 4 different users typing a 10-key password
400 times. The 10 variables measure the time difference between a key being
released and the next one being pressed.

TABLE II: Comparison with Peer EVL methods

Dataset | COMPOSE | LEVELiw | PSO | GA
---|---|---|---|---
$1CDT$ | 0.99 | 0.99 | 0.99(.002) | 0.99(.002)(=)
$4CR$ | 0.99 | 0.99 | 0.99(.001) | 0.99(.001)(=)
$1CSurr$ | 0.90 | 0.63 | 0.95(.004) |  0.95(.024)(+)
$GEARS$ | 0.96 | 0.93 | 0.82(.001) | 0.91(.003)(-)
$UG\\_2C\\_2D$ | 0.94 | 0.73 | 0.91(.000) |  0.96(.001)(+)
$UG\\_2C\\_3D$ | 0.67 | 0.64 | 0.96(.001) |  0.97(.001)(+)
$UG\\_2C\\_5D$ | 0.90 | 0.60 | 0.86(.001) |  0.92(.03(+)
$4CRE\\_V1$ | 0.20 | 0.24 | 0.28(.031) |  0.32(.061)(+)
$4CRE\\_V2$ | 0.19 | 0.24 | 0.25(.061) |  0.33(.067)(+)
$Keystroke$ | 0.85 | 0.78 | 0.73(.023) | 0.76(.004) (-)

###  IV-B Metrics

The performance of the proposed method is measured using the F1-Score
(sometimes called the F-Score or F1-measure). It is the mean of the precision
and recall scores of the algorithm over a given window $w_{i}$. We use the
macro-F1 score as it is suitable to use for balanced and unbalanced data-sets
(the alternative micro-F1 score favours the majority class). The final result
is the average of all of the windows. An incoming point $x_{i}$ with a true
label $y_{i}$ is passed to the ensemble E and a prediction $\hat{y}_{i}$ is
returned; $\hat{y}_{i}=E(x_{i})$.

The predicted values are compared with the actual values and precision and
recall scores calculated. Combining these metrics gives the F1-Score:

$F1=2*\frac{Precision*Recall}{Precision+Recall}$ (2)

When evaluating the population-based algorithms, we run each experiment 20
times. To evaluate whether the resulting distributions are significantly
different, we use the Wilcoxen Signed-Rank test. We reject the null hypothesis
that two distributions are equal with $p<0.05$. When comparing the population
algorithm results with peer results published in the literature we use the
same statistical test and reject the null hypothesis that the distribution of
results are symmetric around the peer results with $p<0.05$

###  IV-C Evaluation

We begin by evaluating the datastreams on a $static$ ensemble. In these
experiments the ensemble is initialised with a labelled training set but is
never updated as the stream progresses. In these instances the ensemble begins
with a high accuracy but does not react to change. As expected, the ensemble
performs poorly over the course of the non-stationary stream. This is because
the constituent classifiers remain specialised in areas of the feature space
which are no longer relevant to the underlying distribution. Virtual concept
drift or real drift has occurred and the ensemble has not reacted resulting in
misclassified instances.

Figure 1: Static Ensemble vs adaptive ensemble

t = 0

t = 350

t = 700

t = 1050

t = 1400

Figure 2: Stream 4CR, 4 classes rotate. Ensemble adapts to underlying change
using Genetic Algorithm to update population of classifiers

As an illustrative example, the performance of the static ensemble is
presented alongside the performance of the adaptive ensemble in Fig. 1. Here,
a GA is used to adapt to the underlying change. The left plot shows the
progression on the 1CDT stream. Here, two classes are present. One is
stationary and the other dynamic. The static ensemble maintains $\approx 50\%$
accuracy because of the stationary class but does not track the non-stationary
class. The GA does and maintains a high accuracy through the stream. In the
right plot, stream 4CR is used. This stream consists of four classes which
rotate and move into positions previously occupied by other classes. The
static ensemble starts with high accuracy but as the stream progresses the
accuracy drops. It rises again once the classes complete their rotation and
return to the starting position. In Fig. 2, this process is visualised; the
adaptation of the ensemble with a Genetic Algorithm is presented. The drift of
the 4 classes can be observed as the algorithm adapts to this drift.

###  IV-D Performance Comparison with Peer Methods

In this subsection, we tabulate the results of the ensembles using a GA and
PSO adaptation method and compare their performance with other EVL methods. We
evaluate the proposed method with COMPOSE [5] and LEVELiw [21]. The peer-
results are taken from results previously published in the literature. The
results are presented in Table II. The proposed approach is favourable to both
peer algorithms.

###  IV-E Complexity Analysis

The memory requirements of the proposed method is $O(NC)$ where $N$ is the
number of agents per class and $C$ is the number of classes in the stream
(assuming this number does not increase during the stream). The time
complexity of the solution is similar. In the classification phase, the time
requirement is also $O(NC)$ and the population update requires a further
$O(NC)$ as this depends on the number of ensembles and agents per ensemble.
This gives a total time complexity of $O(NC^{2})$.

Of both methods evaluated, PSO is the fastest because it requires the least
amount of computation to update the population. However, it returns the lowest
accuracy. Conversely, the GA is slower but returns the highest accuracy. The
time requirement is a function of the number of agents, the more agents the
more calculations and the longer the algorithm takes to process the stream. As
an illustrative example the time requirement for two streams is presented. In
Table III, 1CDT stream is presented and the number of seconds along with the
F1 score is presented as the number of agents in each ensemble increases. PSO
is faster than the GA with the same F1 accuracy and increasing the number of
agents does not improve accuracy despite taking longer to process. In Table
IV, it can be seen that PSO is again much faster but achieves a lower F1
score. The GA takes more time but gives a better accuracy and this accuracy
increases with a larger number of agents.

TABLE III: Number of agents per class in 1CDT Stream | PSO |  | GA
---|---|---|---
Agents | F1 | Secs. |  | F1 | Secs.
10 | 0.98 | 0.3 |  | 0.99 | 0.4
20 | 0.99 | 0.57 |  | 0.99 | 0.8
50 | 0.99 | 1.12 |  | 0.99 | 2.0
100 | 0.99 | 2.4 |  | 0.99 | 4.5
TABLE IV: Number of agents per class in UG_2C_5D Stream | PSO |  | GA
---|---|---|---
Agents | F1 | Secs. |  | F1 | Secs.
10 | 0.81 | 5.5 |  | 0.81 | 7.0
20 | 0.84 | 10.3 |  | 0.85 | 12.6
50 | 0.86 | 19.6 |  | 0.91 | 31.7
100 | 0.86 | 46.1 |  | 0.92 | 73.0

##  V Conclusions and Future Work

In this paper we proposed an evolving population-based approach to the problem
of classification in data-streams with extreme verification latency (EVL). Two
population based approaches were examined; PSO and GA with the GA giving the
best results. These algorithms were evaluated on a range of data-streams and
the population based approach shows promise when compared to peer EVL
algorithms. An ensemble of one-class classifiers is maintained for each class
whereby each classifier is considered an agent in the sub-population. After
each window the population is subject to selection pressure and in this way
the ensemble reacts to underlying change in the stream and maintains high
classification rate throughout the dynamic stream without access to any
further labels beyond the initial training set.

The research serves as an initial study in the use of evolutionary algorithms
for addressing the problem of non-stationary and EVL in data-streams. It also
reveals a number of future directions. In this study, each ensemble is
composed of relatively simple base-classifiers (a micro-classifier, similar to
$k$-NN). Performance could be improved with more sophisticated one-class
classifiers such as Minimum Spanning Trees [14], One Class Support Vector
Machine [17], or others [15], or perhaps a heterogeneous combinations thereof
(for example more heavy-weight sophisticated classifiers in stable areas of
the feature space with simpler lightweight models which can be deployed
quickly and provide a better response to change). There is also potential for
improving/specialising the optimisation techniques for updating the ensemble.
For example, PSO requires fewer calculations but does not provide as accurate
performance. These are considerations which might be stream dependant and this
is an area for future research. The study makes some assumptions which are
common in the EVL literature. Firstly, it is assumed that the data-stream will
not exhibit concept evolution, so classes which are not present in the initial
training data will not appear later in the stream. Secondly, change in the
stream is assumed to be gradual. The proposed method can capture real drift
assuming that the change in $P(y|x)$ begins as virtual drift in $P(x)$.
However, sudden real drift (a change in $P(y|x)$ without change in $P(x)$) is
difficult to recognise without access to the true class label (or a human
expert) and is an important future research direction.

## References

  * [1] C.C. Aggarwal, J. Han, J. Wang, and P.S. Yu. “On demand classification of data streams” In Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining (pp. 503-508). ACM. 2004.
  * [2] M. Arostegi, J.L. Lobo, and J. Del Ser, “SLAYER: A Semi-supervised Learning Approach for Drifting Data Streams under Extreme Verification Latency”. In Workshop on Interactive Adaptive Learning (p. 50). 2021
  * [3] R. Capo, A. Sanchez, and R. Polikar, “Core support extraction for learning from initially labeled nonstationary environments using compose”. In 2014 International Joint Conference on Neural Networks (IJCNN) (pp. 602-608). IEEE. July, 2014.
  * [4] G. Ditzler, M. Roveri, C. Alippi, and R. Polikar. “Learning in nonstationary environments: A survey” IEEE Computational Intelligence Magazine, 10(4), pp.12-25. 2015.
  * [5] K.B. Dyer, R., Capo, and R. Polikar, R.“Compose: A semisupervised learning framework for initially labeled nonstationary streaming data”. IEEE transactions on neural networks and learning systems, 25(1), pp.12-26. 2013.
  * [6] H. Edelsbrunner, D., Kirkpatrick, and R. Seidel, “On the shape of a set of points in the plane”. IEEE Transactions on information theory, 29(4), pp.551-559. 1983.
  * [7] C. Fahy, S. Yang, and M. Gongora, “Classification in Dynamic Data Streams with a Scarcity of Labels”. IEEE Transactions on Knowledge and Data Engineering. 2021
  * [8] C. Fahy.,S. Yang, and M. Gongora, “Scarcity of labels in non-stationary data streams: A survey”. ACM Computing Surveys (CSUR), 55(2), pp.1-39. 2022
  * [9] C. Fahy and S. Yang, “Ant colony stream clustering: A fast density clustering algorithm for dynamic data streams”. IEEE transactions on cybernetics, 49(6), pp.2215-2228.
  * [10] C. Fahy and S. Yang, “Dynamic stream clustering using ants”. In Advances in Computational Intelligence Systems: Contributions Presented at the 16th UK Workshop on Computational Intelligence, September 7–9, 2016, Lancaster, UK (pp. 495-508). Cham: Springer International Publishing.2016
  * [11] R.S. Ferreira, G. Zimbrão, G. and L.G. Alvim, “AMANDA: Semi-supervised density-based adaptive model for non-stationary data with extreme verification latency”. Information Sciences, 488, pp.219-237. 2019.
  * [12] A. Forestiero, C. Pizzuti, and G. Spezzano, “A single pass algorithm for clustering evolving data streams based on swarm intelligence,” Data Mining and Knowledge Discovery, vol. 26, no. 1, pp. 1–26, Nov. 2011.
  * [13] R. Hyde, P. Angelov, and A. R. MacKenzie, “Fully online clustering of evolving data streams into arbitrarily shaped clusters,” Information Sciences, vol. 382-383, pp. 96–114, 2017.
  * [14] P. Juszczak, D.M. Tax, E. Pe, and R.P. Duin, “Minimum spanning tree based one-class classifier”. Neurocomputing, 72(7-9), pp.1859-1869. 2009
  * [15] S.S. Khan, and M.G. Madden, “A survey of recent trends in one class classification”. In Irish conference on artificial intelligence and cognitive science (pp. 188-197). Springer, Berlin, Heidelberg. August, 2009
  * [16] G. Krempl, “The algorithm APT to classify in concurrence of latency and drift,” in Advances in Intelligent Data Analysis (Lecture Notes in Computer Science), vol. 7014, J. Gama, E. Bradley, and J. Hollmen, Eds. Berlin Heidelberg, Germany: Springer-Verlag, pp. 222–233. 2011.
  * [17] J. Mũnoz-Marí, F. Bovolo, L. Gómez-Chova, L. Bruzzone, and G. Camp-Valls, “Semisupervised one-class support vector machines for classification of remote sensing data.” IEEE transactions on geoscience and remote sensing, 48(8), pp.3188-3197. 2010
  * [18] V.M. Souza, D.F. Silva, J. Gama and E. Gustavo. “Data stream classification guided by clustering on nonstationary environments and extreme verification latency” SIAM-ICDM, pp. 873-881. 2015.
  * [19] V.M. Souza, D.F Silva, G.E., Batista, and J. Gama, J., “Classification of evolving data streams with infinitely delayed labels”. In 2015 IEEE 14th International Conference on Machine Learning and Applications (ICMLA) (pp. 214-219). IEEE. December, 2015.
  * [20] V. Souza, T. Pinho, and G. Batista, “Evaluating stream classifiers with delayed labels information”. In 2018 7th Brazilian Conference on Intelligent Systems (BRACIS) (pp. 408-413). IEEE. 2018.
  * [21] M. Umer, R. Polikar, and C. Frederickson, “Level iw: Learning extreme verification latency with importance weighting”. In 2017 International Joint Conference on Neural Networks (IJCNN) (pp. 1740-1747). IEEE. May, 2017.
  * [22] M. Umer, C. Frederickson, C. and R. Polikar, “Learning under extreme verification latency quickly: Fast compose.” In 2016 IEEE Symposium Series on Computational Intelligence (SSCI) (pp. 1-8). IEEE. December, 2016.
  * [23] M. Umer and R. Polikar, ”Comparative analysis of extreme verification latency learning algorithms”. arXiv preprint arXiv:2011.14917. 2020