# Robust Dynamic State Estimation of Multi-Machine Power Networks with Solar
Farms and Dynamics Loads

Muhammad Nadeem and Ahmad F. Taha This work is supported by National Science
Foundation under Grants 2152450 and 2151571. The authors are with the Civil
and Environmental Engineering Department, Vanderbilt University, 2201 West End
Ave, Nashville, Tennessee 37235<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Conventional state estimation routines of electrical grids are mainly reliant
on dynamic models of fossil fuel-based resources. These models commonly
contain differential equations describing synchronous generator models and
algebraic equations modeling power flow/balance equations. Fuel-free power
systems that are driven by inertia-less renewable energy resources will hence
require new models and upgraded estimation routines. To that end, in this
paper we propose a robust estimator for an interconnected model of power
networks comprised of a comprehensive ninth order synchronous generator model,
advanced power electronics-based models for photovoltaic (PV) power plants,
constant power loads, constant impedance loads, and motor loads. The presented
state estimator design is based on Lyapunov stability criteria for nonlinear
differential algebraic equation (DAE) models and is posed as a convex semi-
definite optimization problem. Thorough simulations studies have been carried
out on IEEE-39 bus test system to showcase the robustness of the proposed
estimator against unknown uncertainty from load demand and solar irradiance.

##  I Introduction and Paper Contributions

Power systems are moving rapidly toward a greener and more sustainable future
through the widespread installation of renewable energy resources (RERs). The
intermittent and volatile nature of RERs has created challenges in modern
power systems and made the secure and reliable operation of the grid more
difficult. The solution to these problems can be found in dynamic state
estimation (DSE) for the next generation of smart power grids. In particular
DSE can provide accurate realtime estimates of the physical states of
electrical grid which can be highly beneficial for the realtime feedback
control and health monitoring of electric grid [1].

In this regard, thorough research has been carried out in the past two decades
in power system DSE mainly focusing on estimating generator internal states
(rotor angle and frequency) and algebraic states (voltage and current
phasors). All of the current DSE algorithms can be divided into two broad
categories: (1) deterministic observers and (2) Kalman filters (also known as
stochastic estimators). In deterministic observers the basic idea is as
follows: First, a dynamical model for the error dynamics is computed (this
model explains the evolution of error between original and estimated state
variables), then observer gain is designed based on Lyapunov stability notion
such that it drives the error model near zero [2, 3, 4].

On the other hand Kalman filters (KF) are iterative and they primarily exploit
the statistical properties (such as normal distribution) of the
noise/disturbances to minimize their impact on the performance of DSE [5, 6,
7].

A summary of all the techniques used for power system DSE can be found in
recent survey paper [1], while a thorough comparative study between different
derivatives of Kalman filter can be found in [8]. Similarly, a comparative
analysis between stochastic estimators and deterministic observers can also be
seen in [9]. An extensive summary of these methods is beyond the scope of this
paper.

In the current literature, the DSE algorithms are primarily based on a simple
power system model capturing only dynamical model of generators and few
algebraic constraints showcasing power balance equations. Thus the monitoring
scope of DSE is mainly limited to the estimation of algebraic states and
generator internal states. However, with the rapid deployment of solar PV
power plants and dynamic loads in the electrical grid these traditional power
systems models are gradually becoming obsolete and out of context [1]. Hence
to achieve grid-wise situational awareness in the future power systems it is
important to model solar PV power plants and dynamic loads in the nonlinear
descriptor (NL-D) model of power system while performing DSE.

Unfortunately, the literature on state estimation of multi-machine power
network models with synchronous machines, solar PV power plants, and motor
loads while considering PMU measurements models and algebraic constraints is
virtually non-existent. In the recent study [10] DSE has been extended and
state estimation for a detailed PV power plant model has been reported,
however, the study is limited to distribution network only and does not
estimate algebraic variables and generators’ internal states.

Based on the above discussion the contributions of this paper are as follows:

  * •

We propose DSE technique for a holistic model of a power system having ninth
order synchronous generator model, advanced power electronics-based models of
solar farms, motor loads, constant power loads, and constant impedance loads.
The proposed DSE algorithm can simultaneously estimate all the states of the
power system that include dynamic states of solar farm, loads and algebraic
states.

  * •

To design a robust estimator that can handle load fluctuations and uncertainty
from RERs, we propose an $H_{\infty}$ based state estimation technique. The
main advantage of $H_{\infty}$ based state estimator over Kalman filters is
that no prior knowledge about the statistic of uncertainty is required.

  * •

We showcase the performance of the proposed estimator on IEEE case39 test
power system which is widely used for DSE studies in power systems.

The rest of the paper is organized as follows: Section II summarizes the NL-D
model; Section III discuses the proposed estimator design; Section IV
showcases simulation studies and the paper is concluded in Section V.

##  II Power Network Model

We consider a graphical model of the electrical grid with $N$ number of buses,
$\mathcal{N}=\left\\{1,...,N\right\\}$ set of nodes and
$\mathcal{E}\subseteq\mathcal{N}\times\mathcal{N}$ set of transmission lines.
The overall power system is assumed to have: $G$ number of synchronous
machines, $J$ number of PV power plants, and $L_{p}$, $L_{z}$, and $K$ number
of constant power, constant impedance, and motor loads respectively. Notice
that $\mathcal{N}=\mathcal{G}\cup\mathcal{R}\cup\mathcal{L}$, where
$\mathcal{G},\mathcal{R}$, and $\mathcal{L}$ denote a set of buses containing
synchronous machines, renewables, and loads respectively. To that end, we
model the interconnected model of power systems using a set of differential-
algebraic equations as follows:

$\displaystyle\vspace{-0.5cm}\textit{system
dynamics:}\;\;\;\dot{\boldsymbol{x}}(t)$
$\displaystyle=\boldsymbol{f}(\boldsymbol{x}_{d},\boldsymbol{x}_{a},\boldsymbol{u},\boldsymbol{w})~{}$
(1a) $\displaystyle\textit{algebraic constraints:}\;\;\;\boldsymbol{0}$
$\displaystyle=\boldsymbol{h}(\boldsymbol{x}_{d},\boldsymbol{x}_{a},\boldsymbol{u},\boldsymbol{w})~{}$
(1b)

where $\boldsymbol{x}_{d}\in\mathbb{R}^{n_{d}}$ denotes differential
variables, $\boldsymbol{x}_{a}\in\mathbb{R}^{n_{a}}$ denotes algebraic
variables, $\boldsymbol{w}\in\mathbb{R}^{n_{w}}$ models the exogenous
disturbances to the system, and $\boldsymbol{u}\in\mathbb{R}^{n_{u}}$ models a
set of commands that are used to drive the grid to a desirable state.

In (1) we model $\boldsymbol{x}_{a}$ as:

$\displaystyle\boldsymbol{x}_{a}(t):=\boldsymbol{x}_{a}=\begin{bmatrix}\boldsymbol{V}_{Re}^{\top}&\boldsymbol{V}_{Im}^{\top}&\boldsymbol{I}_{Re}^{\top}&\boldsymbol{I}_{Im}^{\top}\end{bmatrix}^{\top}\vspace{-0.3cm}$
(2)

where
$\boldsymbol{V}_{Re}=\\{V_{Re_{i}}\\}_{i\in\mathcal{N}},\boldsymbol{V}_{Im}=\\{V_{Im_{i}}\\}_{i\in\mathcal{N}},\boldsymbol{I}_{Re}=\\{I_{Re_{i}}\\}_{i\in\mathcal{N}}$
and $\boldsymbol{I}_{Im}=\\{I_{Im_{i}}\\}_{i\in\mathcal{N}}$ denotes the real
and imaginary part of voltage and current phasors respectively.

The vector $\boldsymbol{u}$ contains reference inputs to the generators and PV
plant and is modeled as:

$\displaystyle\boldsymbol{u}=\begin{bmatrix}\boldsymbol{u}_{G}^{\top}&\boldsymbol{u}_{R}^{\top}\end{bmatrix}^{\top}$
(3)

where $\boldsymbol{u}_{G}$ contains voltage reference set points
$\boldsymbol{V}_{ref}$ and steam/hydro valve reference positions
$\boldsymbol{P}_{v_{ref}}$ for the synchronous machines, such that
$\boldsymbol{u}_{G}=\begin{bmatrix}\boldsymbol{V}_{ref}^{\top}&\boldsymbol{P}_{v_{ref}}^{\top}\end{bmatrix}$.
Similarly
$\boldsymbol{u}_{R}=\begin{bmatrix}\boldsymbol{V}_{ref}^{\top}&\boldsymbol{P}_{ref}^{\top}\end{bmatrix}$,
where $\boldsymbol{V}_{ref}$ is the voltage reference and
$\boldsymbol{P}_{ref}$ is the power reference set points for PV power plants.

Similarly, vector $\boldsymbol{w}$ in Eq. (1) is expressed as follows:

$\displaystyle\boldsymbol{w}=\begin{bmatrix}\boldsymbol{P}_{d}^{\top}&\boldsymbol{I}_{r}^{\top}\end{bmatrix}^{\top}$
(4)

where $\boldsymbol{P}_{d}$ is the disturbance in load demand and
$\boldsymbol{I}_{r}$ is the disturbance in the irradiance from the sun.
Furthermore in (1) we model $\boldsymbol{x}_{d}$ as:

$\displaystyle\boldsymbol{x}_{d}(t):=\boldsymbol{x}_{d}=\begin{bmatrix}\boldsymbol{x}_{G}^{\top}&\boldsymbol{x}_{R}^{\top}&\boldsymbol{x}_{m}^{\top}\end{bmatrix}^{\top}\vspace{-0.1cm}$
(5)

where $\boldsymbol{x}_{G}$ denote states of generators, $\boldsymbol{x}_{R}$
denote dynamic states of PV plant, and $\boldsymbol{x}_{m}$ represents dynamic
states of motor load. To that end, we model synchronous machines using a
comprehensive ninth order generator model, thus $\boldsymbol{x}_{G}$ can be
given as follows [11, 12]:

$\displaystyle\boldsymbol{x}_{G}=\begin{bmatrix}\boldsymbol{e}_{dq}^{\top}\;\;\boldsymbol{\omega}_{sg}^{\top}\;\;\boldsymbol{\delta}_{sg}^{\top}\;\;\boldsymbol{T}_{m}^{\top}\;\;\boldsymbol{P}_{v}^{\top}\;\;\boldsymbol{E}_{fd}^{\top}\;\;\boldsymbol{v}_{a}^{\top}\;\;\boldsymbol{r}_{f}^{\top}\end{bmatrix}^{\top}\in\mathbb{R}^{9G}$
(6)

where
$\boldsymbol{e}_{dq}=[\boldsymbol{e}_{d}^{\top}\,\,\boldsymbol{e}_{q}^{\top}]$
denotes transient voltages along q- and d-axis, $\boldsymbol{\omega}_{sg}$ is
the rotor speed, $\boldsymbol{\delta}_{G}$ is the rotor angle,
$\boldsymbol{T}_{m}$ represent prime mover torque, $\boldsymbol{P}_{v}$ is the
valve position, $\boldsymbol{E}_{fd}$ denotes field voltage,
$\boldsymbol{v}_{a}$ represents amplifier voltage, and $\boldsymbol{r}_{f}$
denotes stabilizer output.

The models for the PV power plant has been obtained from [12] and thus the
state vector $\boldsymbol{x}_{R}$ can be expressed as follows:

$\displaystyle\small\boldsymbol{x}_{R}=\begin{bmatrix}\boldsymbol{i}_{dq_{f}}^{\top}\;\;\boldsymbol{v}_{dq_{c}}^{\top}\;\;\boldsymbol{E}_{dc}^{\top}\;\;\boldsymbol{P}_{e}^{\top}\;\;\boldsymbol{Q}_{e}^{\top}\;\;\boldsymbol{\delta}_{inv}^{\top}\;\;\boldsymbol{z}_{dq_{0}}^{\top}\;\;\boldsymbol{z}_{dq_{f}}^{\top}\end{bmatrix}^{\top}\in\mathbb{R}^{12J}$
(7)

where $\boldsymbol{i}_{dq,f}$ denotes the _dq_ -axis current at the terminals
of the inverter of the PV plant, $\boldsymbol{v}_{dq,c}$ is the _dq_ -axis
voltage across the AC capacitor, $\boldsymbol{E}_{dc}$ the energy stored in
the dc-link capacitor of the PV plant, $\boldsymbol{{P}_{e}}$ and
$\boldsymbol{{Q}_{e}}$ are the filtered real and reactive power at the
terminals of the PV power plant, $\boldsymbol{\delta_{inv}}$ is the relative
angle of the PV power plant, $\boldsymbol{z}_{dq,o}$ and
$\boldsymbol{z}_{dq,f}$ are the states of the voltage and current regulators
used in the grid-forming controller of the PV plant. Readers are referred to
[12] for further details and complete in-depth description about the PV plant
model.

The dynamics of the motor loads are detailed as [13]:

$\dot{\omega}_{mot_{k}}=\frac{1}{2H_{m_{k}}}(T_{e_{k}}-T_{m_{k}})$ (8)

where $H_{m_{k}}$ is the motor inertia constant and $T_{e_{k}}$, $T_{m_{k}}$
denote electromagnetic and mechanical torques in the $k$-th motor [13, p.
244]. Thus $\boldsymbol{x}_{m}=[\boldsymbol{w}_{mot_{k}}]$.

To take into account the topological effect of power systems, the power
flow/balance or the current balance equations need to considered. These can be
expressed as follows [11]:

$\displaystyle\underbrace{\begin{bmatrix}\boldsymbol{\widetilde{I}}_{G}\\\
\boldsymbol{\widetilde{I}}_{R}\\\
\boldsymbol{\widetilde{I}}_{L}\end{bmatrix}}_{\boldsymbol{I}}-\underbrace{\begin{bmatrix}\boldsymbol{Y}_{GG}&\boldsymbol{Y}_{GR}&\boldsymbol{Y}_{GL}\\\
\boldsymbol{Y}_{RG}&\boldsymbol{Y}_{RR}&\boldsymbol{Y}_{RL}\\\
\boldsymbol{Y}_{LG}&\boldsymbol{Y}_{LR}&\boldsymbol{Y}_{LL}\end{bmatrix}}_{\boldsymbol{Y}}\underbrace{\begin{bmatrix}\boldsymbol{\widetilde{V}}_{G}\\\
\boldsymbol{\widetilde{V}}_{R}\\\ \boldsymbol{\widetilde{V}}_{L}\\\
\end{bmatrix}}_{\boldsymbol{V}}=\boldsymbol{0}$ (9)

where $\boldsymbol{I}$ is the net injected current vector, $\boldsymbol{Y}$ is
the admittance matrix, and $\boldsymbol{V}$ is the bus voltage vector. In (9),
$\boldsymbol{\widetilde{I}}_{G}=\\{I_{Re_{i}}\\}_{i\in\mathcal{G}}+j\\{I_{Im_{i}}\\}_{i\in\mathcal{G}}$
denotes phasor currents injected by synchronous generators and
$\boldsymbol{\widetilde{V}}_{G}=\\{V_{Re_{i}}\\}_{i\in\mathcal{G}}+j\\{V_{Im_{i}}\\}_{i\in\mathcal{G}}$
represents voltage phasors at the terminal of generator buses. Similarly
$\boldsymbol{\widetilde{V}}_{R}$ , $\boldsymbol{\widetilde{V}}_{L}$,
$\boldsymbol{\widetilde{I}}_{R}$, and $\boldsymbol{\widetilde{I}}_{L}$ denotes
voltage and current phasors of all loads and PV power plants.

To that end, by considering (2)–(9) and incorporating the associated dynamics
then the overall interconnected model of power systems can be represented in a
following compact form:

$\displaystyle\textbf{NL-D:}\;\;\;\boldsymbol{E}\dot{{\boldsymbol{x}}}$
$\displaystyle={\boldsymbol{A}}{\boldsymbol{x}}+{\boldsymbol{f}}\left({\boldsymbol{x}},{\boldsymbol{u}},{\boldsymbol{w}}\right)+{\boldsymbol{B}}{{\boldsymbol{u}}}+{\boldsymbol{B}}_{w}\boldsymbol{w}$
(10a) $\displaystyle\boldsymbol{y}$
$\displaystyle=\boldsymbol{C}\boldsymbol{x}+\boldsymbol{v}$ (10b)

where
$\boldsymbol{x}=\begin{bmatrix}\boldsymbol{x}_{d}^{\top}&\boldsymbol{x}_{a}^{\top}\end{bmatrix}^{\top}\in\mathbb{R}^{n}$
denotes the overall state vector; $\boldsymbol{E}$ encodes algebraic
constraints with rows of zeros; $\boldsymbol{C}$ maps state vector
$\boldsymbol{x}$ to what PMUs usually measures (i.e., current and voltage
phasors) and $\boldsymbol{v}\in\mathbb{R}^{p}$ denotes measurement noise on
PMU measurements $\boldsymbol{y}\in\mathbb{R}^{p}$. Since vector
$\boldsymbol{x}_{a}$ contains voltage and current phasors, we define
$\boldsymbol{C}$ as:
$\boldsymbol{C}=\begin{bmatrix}\boldsymbol{O}&\tilde{\boldsymbol{C}}\end{bmatrix}$,
where $\tilde{\boldsymbol{C}}$ is a diagonal binary matrix with ones only at
those locations where PMUs are connected and voltage and current phasors are
measured. The other state-space matrices $\boldsymbol{A},\boldsymbol{B}$, and
$\boldsymbol{B}_{w}$ are obtained by capturing the linear components of the
SLS-DAE model. The vector-valued function $\boldsymbol{f}(\cdot)$ captures the
encompassed nonlinearities in the dynamics.

##  III Estimator for NL-D Model of Power Systems

In this section, we propose a Luenberger type observer design for the NL-D
model depicted in (11). First, we focus on modeling the uncertainties in loads
and renewables and present a technique to parameterize the nonlinear function
$\boldsymbol{f}(.)$ in the power system dynamics.

###  III-A Modeling uncertainty from loads and renewables

In Eq. (11) the vector $\boldsymbol{w}$ encapsulates load demand and
irradiance, note that both of these quantities are time-varying and
fluctuating. To that end, herein we assume that hour- or minute-ahead
predictions of loads and irradiance are available while the
fluctuation/disturbances in these quantities are unknown. Notice that this is
realistic as power systems operators record and publish these quantities on a
daily basis (see daily hour- and minute-ahead predictions of load and
renewables published by California independent system operator [14]). However,
the prediction may be inaccurate. In particular, high fidelity estimate of
RERs are difficult to obtain. Accordingly, we can write
$\boldsymbol{w}=\bar{\boldsymbol{w}}+\Delta\boldsymbol{w}$, where
$\bar{\boldsymbol{w}}$ is the predicted/known values and
$\Delta\boldsymbol{w}$ defines all the disturbances/uncertainties in these
quantities. The goal of the estimator is to provide accurate state estimation
results under unknown uncertainty $\Delta\boldsymbol{w}$. To that end, the
NL-D model (11) can be rewritten as follows:

$\displaystyle\begin{split}\boldsymbol{E}\dot{{\boldsymbol{x}}}&={\boldsymbol{A}}{\boldsymbol{x}}+{\boldsymbol{f}}\left({\boldsymbol{x}},{\boldsymbol{u}},\boldsymbol{w}\right)+{\boldsymbol{B}}{{\boldsymbol{u}}}+{\boldsymbol{B}}_{w}\bar{\boldsymbol{w}}+{\boldsymbol{B}}_{w}\Delta\boldsymbol{w}\end{split}$
(11a)
$\displaystyle\begin{split}\boldsymbol{y}&=\boldsymbol{C}\boldsymbol{x}+\boldsymbol{v}.\end{split}$
(11b)

The main reason for splitting the vector $\boldsymbol{w}$ into known/predicted
and unknown parts is because later on in the estimator dynamics to achieve
robust performance the estimator only has access to the predicted value of
$\boldsymbol{w}$ and disturbances $\Delta\boldsymbol{w}$ is kept unknown to
the estimator.

In the following section we discuss Lipschitz continuity condition to classify
the nonlinear vector valued function $\boldsymbol{f}(.)$ in the system
dynamics.

###  III-B Parameterizing nonlinearities in the NL-D model

To synthesize robust observers it is crucial to identify and parameterize the
nonlinearities in the dynamical system. Hence, to express the nonlinear
function $\boldsymbol{f}(.)$ in a better way in the NL-D model we assume that
$\boldsymbol{f}(.)$ is Lipschitz continuous, which means $\boldsymbol{f}(.)$
is differentiable everywhere and the rate of change is bounded above by a real
number. Smallest such real number is called Lipschitz constant. With that in
mind we assume $\boldsymbol{x}\in\mathcal{X}$, where $\mathcal{X}$ define the
operating region of state vector $\boldsymbol{x}$, then the Lipschitz bounding
condition for $\boldsymbol{f}(.)$ can be written as follows:

$\displaystyle\lVert\boldsymbol{f}(\boldsymbol{x},\boldsymbol{u},\boldsymbol{w})-\boldsymbol{f}(\hat{\boldsymbol{x}},\boldsymbol{u},\boldsymbol{w})\rVert_{2}\leq\alpha\lVert(\boldsymbol{x}-\hat{\boldsymbol{x}})\rVert_{2}$
(12)

where $\alpha$ is the Lipschitz constant. In the case of power systems the
Lipschitz constant can be assumed based on operator knowledge and by looking
at the overall size/topology of the grid or it can also be computed using much
more sophisticated ways as shown in [15, 16]. Also, since the Lipschitz
constant provides an operational bound on grid state variables, its value can
be varied to consider anomalies and deviations from steady state values and
vice versa [2]. Showcasing methods to compute Lipschitz constant is beyond the
scope of this paper.

###  III-C $H_{\infty}$ based state estimator design

In this study we are designing a Luenberger type state estimator for NL-D
model of power systems depicted in (11). The overall estimator design is
primarily based on Lyapunov stability criterion and $H_{\infty}$ estimation
concept is used to achieve robust performance under unknown
fluctuations/disturbances from loads and RERs.

In state estimation literature, the $H_{\infty}$ concept was initially
introduced in [17] to synthesize a robust estimator for linear systems subject
to unknown uncertainty. The main advantage of $H_{\infty}$ based state
estimation is that no prior knowledge or assumptions are required about the
statistical properties of uncertainty. In $H_{\infty}$ based state estimation
the uncertainty is considered as a random bounded signal and then the observer
is designed such that it ensures a particular $H_{\infty}$ performance for the
error dynamics for all such bounded uncertainty, such that:
$\lVert{\boldsymbol{e}}\rVert^{2}_{L_{2}}\leq\gamma^{2}\lVert{\Delta\boldsymbol{w}}\rVert^{2}_{L_{2}}$
with performance level $\gamma$.

That being said, we now focus on presenting the structure of our observer
design for the NL-D model of power system depicted in (11). To begin, let
$\hat{\boldsymbol{y}}$ be the estimated outputs and $\hat{\boldsymbol{x}}$ be
the estimated states variables. Then Luenberger type state estimator for Eq.
(11) can be expressed as follows:

$\displaystyle\begin{split}\boldsymbol{E}\dot{\hat{\boldsymbol{x}}}&={\boldsymbol{A}}{\hat{\boldsymbol{x}}}+{\boldsymbol{f}}\left({\hat{\boldsymbol{x}}},{\boldsymbol{u}},{\bar{\boldsymbol{w}}}\right)+\boldsymbol{L}(\boldsymbol{y}-{\hat{\boldsymbol{y}})}+{\boldsymbol{B}}{{\boldsymbol{u}}}+{\boldsymbol{B}_{w}}{{\bar{\boldsymbol{w}}}}\end{split}$
(13a)
$\displaystyle\begin{split}\hat{\boldsymbol{y}}&=\boldsymbol{C}\hat{\boldsymbol{x}}\end{split}$
(13b)

where $\boldsymbol{L}\in\mathbb{R}^{n\times p}$ is the estimator gain matrix.
Even if the estimator starts from different initial conditions, using
measurement provided by PMUs $\boldsymbol{y}$, the gain matrix
$\boldsymbol{L}$ ensures the convergence of estimated states
$\hat{\boldsymbol{x}}$ to their true values $\hat{\boldsymbol{x}}$. The main
objective of this paper is to design an appropriate gain matrix
$\boldsymbol{L}$ such that robust state estimation can be achieved from the
estimator dynamics (13).

To that end, we define the error between actual and true values of state
variables as $\boldsymbol{e}=\boldsymbol{x}-\hat{\boldsymbol{x}}$. Then using
(11) and (13) the model for the error dynamics can be computed as:

$\displaystyle\boldsymbol{E}\dot{{\boldsymbol{e}}}$
$\displaystyle=\boldsymbol{A}_{c}\boldsymbol{e}+\Delta\boldsymbol{f}+\boldsymbol{B}_{w}\Delta\boldsymbol{w}\vspace{-0.2cm}$
(14)

where $\boldsymbol{A}_{c}=(\boldsymbol{A}-\boldsymbol{L}\boldsymbol{C})$ and
$\Delta\boldsymbol{f}=\boldsymbol{f}\left({\boldsymbol{x}},{\boldsymbol{u}},{\boldsymbol{w}}\right)-{\boldsymbol{f}}\left({\hat{\boldsymbol{x}}},{\boldsymbol{u}},{\boldsymbol{w}}\right)$.
The primary objective of the estimator is to drive the differential equation
of the error (14) near zero under unknown uncertainty $\Delta\boldsymbol{w}$.
Before we proceed to the theory of designing the gain matrix $\boldsymbol{L}$,
we present key assumption that is critical in developing the proposed
methodology.

###### Assumption 1

The effect of uncertainty $\Delta\boldsymbol{w}$ on the norm of
$\boldsymbol{f}(.)$ is negligible, meaning:

$\displaystyle\lVert\boldsymbol{f}(\hat{\boldsymbol{x}},\boldsymbol{u},\bar{\boldsymbol{w}})\rVert_{2}\approx\lVert\boldsymbol{f}(\hat{\boldsymbol{x}},\boldsymbol{u},{\boldsymbol{w}})\rVert_{2}.$
(15)

Assumption 1 is crucial in developing the theory of the observer design as it
is required while applying the S-procedure lemma [18] in the derivation of
estimator gain in a tractable way. Notice that Assumption 1 is mild and holds
in the case of power networks as the nonlinear function in the power systems
models are commonly independent of uncertainties and are mostly dependent only
on state vector, see [19] and references therein. To validate this on the NL-D
model (13) we used different IEEE test power systems such as WECC-9 bus system
and IEEE-39 bus system and we note that by varying $\boldsymbol{P}_{d}$ in
$\boldsymbol{w}$ the norm of $\boldsymbol{f}(.)$ does not change while by
adding uncertainty in $\boldsymbol{I}_{r}$ the norm of $\boldsymbol{f}(.)$
changes slightly, such as by decreasing the irradiance by $40\%$ the
percentage change in $\lVert\boldsymbol{f}(.)\rVert_{2}$ is $0.09\times
10^{-3}\%$. Hence Assumption 1 is empirically satisfied. Proving Assumption 1
theoretically is an interesting research question beyond the scope of this
paper.

With that in mind, we can find necessary and sufficient conditions for the
existence of observer gain matrix $\boldsymbol{L}$. To that end, let us assume
a Lyapunov candidate function
$V(\boldsymbol{e})=\boldsymbol{e}^{\top}\boldsymbol{E}^{\top}\boldsymbol{P}\boldsymbol{e}$,
where $V:\mathbb{R}^{n}\rightarrow\mathbb{R}_{+}$,
$\boldsymbol{P}\in\mathbb{R}^{n\times n}$, and
$\boldsymbol{E}^{\top}\boldsymbol{P}=\boldsymbol{P}^{\top}\boldsymbol{E}\succeq
0$, then its derivative along the trajectories of Eq. (14) yields:

$\displaystyle{\dot{V}}(\boldsymbol{e})$
$\displaystyle=(\boldsymbol{E}\dot{\boldsymbol{e}})^{\top}\boldsymbol{P}\boldsymbol{e}+(\boldsymbol{P}\boldsymbol{e})^{\top}(\boldsymbol{E}\dot{\boldsymbol{e}}).$

For any bounded disturbances $\Delta\boldsymbol{w}$ the $H_{\infty}$ stability
condition can be satisfied as:
${\dot{V}}(\boldsymbol{e})+\boldsymbol{e}^{\top}\boldsymbol{e}-\gamma^{2}\Delta\boldsymbol{w}^{\top}\Delta\boldsymbol{w}<0$,
which can be written as:
$\boldsymbol{\Upsilon}^{\top}\boldsymbol{\Theta}_{1}\boldsymbol{\Upsilon}<0$,
where

$\displaystyle\boldsymbol{\Upsilon}=\begin{bmatrix}\boldsymbol{e}\\\
\Delta\boldsymbol{f}\\\
\Delta\boldsymbol{w}\end{bmatrix},\boldsymbol{\Theta}_{1}=\begin{bmatrix}\boldsymbol{A}_{c}^{\top}\boldsymbol{P}+\boldsymbol{P}^{\top}\boldsymbol{A}_{c}+\boldsymbol{I}&\boldsymbol{P}^{\top}&\boldsymbol{P}^{\top}\boldsymbol{B}_{w}\\\
\boldsymbol{P}&\boldsymbol{O}&\boldsymbol{O}\\\
\boldsymbol{B}_{w}^{\top}\boldsymbol{P}&\boldsymbol{O}&-{\gamma}^{2}\boldsymbol{I}\end{bmatrix}.$

Now from Lipschitz continuity assumption (12) we obtain:
$\Delta\boldsymbol{f}^{\top}\Delta\boldsymbol{f}-\alpha^{2}\boldsymbol{e}^{\top}\boldsymbol{e}\leq
0$, which can be written as
$\boldsymbol{\Upsilon}^{\top}\boldsymbol{\Theta}_{2}\boldsymbol{\Upsilon}\leq
0$, with
$\boldsymbol{\Theta}_{2}=\mathrm{diag}\left(\begin{bmatrix}-\alpha^{2}\boldsymbol{I}&\boldsymbol{I}&\boldsymbol{O}\end{bmatrix}\right)$,
where $\mathrm{diag}$ denotes a diagonal matrix. Now applying the S-procedure
lemma [18], then
$\boldsymbol{\Theta}_{1}-(\epsilon)\boldsymbol{\Theta}_{2}\prec 0$ for some
scalar $\epsilon>0$, which is equivalent to:

$\displaystyle\begin{bmatrix}\boldsymbol{A}_{c}^{\top}\boldsymbol{P}+\boldsymbol{P}^{\top}\boldsymbol{A}_{c}+\boldsymbol{I}+\epsilon\alpha^{2}\boldsymbol{I}&\boldsymbol{P}^{\top}&\boldsymbol{P}^{\top}\boldsymbol{B}_{w}\\\
\boldsymbol{P}&-\epsilon\boldsymbol{I}&\boldsymbol{O}\\\
\boldsymbol{B}_{w}^{\top}\boldsymbol{P}&\boldsymbol{O}&-{\gamma}^{2}\boldsymbol{I}\end{bmatrix}\prec
0.$ (16)

Finally, by rewriting $\boldsymbol{P}$ as
$\boldsymbol{P}=\boldsymbol{XE}+{\boldsymbol{E}^{\perp\top}}{\boldsymbol{Y}}$
for some $\boldsymbol{Y}\in\mathbb{R}^{n_{a}\times n}$ and
$\boldsymbol{X}\in\mathbb{S}^{n\times n}_{++}$, where
$\boldsymbol{E}^{\perp\top}$ is the orthogonal complement of $\boldsymbol{E}$
[20, ch. 2] and defining
$\boldsymbol{N}=\boldsymbol{L}^{\top}\boldsymbol{P}\in\mathbb{R}^{p\times n}$
then the necessary and sufficient condition for the existence of observer gain
$\boldsymbol{L}$ can be written in term of strict linear matrix inequality
(LMI) as follows:

$\displaystyle\begin{bmatrix}\boldsymbol{\Omega}_{11}&*&*\\\
\boldsymbol{X}\boldsymbol{E}+\boldsymbol{E}^{\perp\top}\boldsymbol{Y}&-{\epsilon}\boldsymbol{I}&*\\\
\boldsymbol{B_{w}}^{\top}\boldsymbol{X}\boldsymbol{E}+\boldsymbol{B_{w}}^{\top}\boldsymbol{E}^{\perp\top}\boldsymbol{Y}&\boldsymbol{O}&-{\gamma}^{2}\boldsymbol{I}\end{bmatrix}\prec
0$ (17)

where $\boldsymbol{\Omega}_{11}$ is given as:

$\displaystyle\begin{split}\boldsymbol{\Omega}_{11}=\boldsymbol{A}^{\top}\boldsymbol{X}\boldsymbol{E}+\boldsymbol{A}^{\top}\boldsymbol{E}^{\perp\top}\boldsymbol{Y}+\boldsymbol{E}^{\top}\boldsymbol{X}^{\top}\boldsymbol{A}+\boldsymbol{Y}^{\top}\boldsymbol{E}^{\perp}\boldsymbol{A}-\\\
\boldsymbol{C}^{\top}\boldsymbol{N}-\boldsymbol{N}^{\top}\boldsymbol{C}+\epsilon\alpha^{2}\boldsymbol{I}.\end{split}$

If LMI (17) is solved and there exist matrices $\boldsymbol{N}$,
$\boldsymbol{Y}$, $\boldsymbol{X}$, and $\gamma$ and $\epsilon>0$, then the
estimator gain can be retrieved as
$\boldsymbol{L}=\left({\boldsymbol{N}\boldsymbol{P}^{-1}}\right)^{\top}$.

Notice that, since $\gamma$ denotes performance level, then one can minimize
$\gamma$ to achieve robust performance. Also, minimizing
$\lVert\boldsymbol{N}\rVert_{2}$ can provide estimator gain of reasonable
magnitude. Furthermore, as our Lyapunov candidate function can be written as
$V(e)=\boldsymbol{e}^{\top}\boldsymbol{E}^{\top}\boldsymbol{X}\boldsymbol{E}\boldsymbol{e}$
then to achieve quick convergence of $\hat{\boldsymbol{x}}$ to
$\boldsymbol{x}$ one can minimize the maximum eigenvalue of
$\boldsymbol{E}^{\top}\boldsymbol{X}\boldsymbol{E}$ which can be written as a
convex SDP optimization problem:

$\displaystyle\mathbf{\left(P_{1}\right)}\;\operatorname*{minimize}_{\kappa,\boldsymbol{X}}\;\;\;$
$\displaystyle
a_{3}\kappa\;\;\;\;\mathrm{subject\;to}\,\,\;\;\;\kappa\boldsymbol{I}-\boldsymbol{E}^{\top}\boldsymbol{X}\boldsymbol{E}\succ
0,\kappa>0.$

Hence, the necessary and sufficient condition of computing $\boldsymbol{L}$
given in Eq. (17) can be converted to a convex semidefinite optimization
problem expressed as follows:

$\displaystyle\mathbf{\left(P_{2}\right)}\;\operatorname*{minimize}_{\kappa,{\epsilon},\gamma,\boldsymbol{N},\boldsymbol{X},\boldsymbol{Y}}\;\;\;$
$\displaystyle a_{1}\lVert\boldsymbol{N}\rVert_{2}+a_{2}\gamma+a_{3}\kappa$
$\displaystyle\operatorname*{subject\,to}\;\;\;$
$\displaystyle\kappa\boldsymbol{I}-\boldsymbol{E}^{\top}\boldsymbol{X}\boldsymbol{E}\succ
0,\;\mathrm{LMI}\;\eqref{eq:LMI_Hinf},\;$ $\displaystyle\boldsymbol{X}\succ
0,\;\kappa>0,\;{\epsilon}>0,\;\gamma>0$

where $a_{1}$, $a_{2}$, $a_{3}$ $\in\mathbb{R_{++}}$ are predefined weighting
constants. Calculating estimator gain matrix $\boldsymbol{L}$ by solving
$\mathbf{P_{2}}$ always ensures that the performance of state error dynamics
(14) is bounded in $H_{\infty}$ sense, such that
$\lVert{\boldsymbol{e}}\rVert^{2}_{L_{2}}\leq\gamma^{2}\lVert{\Delta\boldsymbol{w}}\rVert^{2}_{L_{2}}$.
In other words, solving $\mathbf{P_{2}}$ guarantees that
$\lVert{\boldsymbol{e}}\rVert^{2}_{L_{2}}$ always lies in a tube of circle
with center at the origin and having a radius of
$\gamma^{2}\lVert{\boldsymbol{w}}\rVert^{2}_{L_{2}}$.

In the following section we present thorough numerical studies to assess the
performance of the estimator under various uncertainties from loads and
renewables.

##  IV System Setup and Simulation Studies

The proposed estimator has been tested on modified IEEE-39 bus system [21].
This system consists of two PV power plants at Buses 34 and 36, a motor load
at Bus 14, eight synchronous generators, constant power loads and constant
impedance loads. The one-line diagram of this system is depicted in Fig. 1.
All the parameters for the synchronous generator and its excitation systems
are detailed in [11], while the parameters of PV power plants and other
complete in-depth description about the system model can be found in [12].

Figure 1: IEEE 39 bus system with two PV power plants at Buses 34, 36 and a
motor load at Bus 14.

We deploy 13 PMUs at Buses $[2,6,9,10,13,14,17,19,20,22,23,25,29]$ as given in
[22] to ensure that the power system is observable. Each PMU at a bus is
measuring the total current demanded/injected and voltage phasors of that bus.
The dynamic states that we are estimating are $\boldsymbol{x}_{a}$ and
$\boldsymbol{x}_{d}$ as given in (2) and (5). The uncertain quantities for the
estimator are assumed to be the overall load demand and irradiance from the
sun. Moreover, for all the case studies Gaussian noise with diagonal
covariance matrix and variance of $0.001^{2}$ has also been added to the PMU
measurements.

Both the estimator and power system NL-D models are simulated using MATLAB
differential algebraic equations solver ode15s, whereas all the SDPs are
solved in YALMIP [23] with MOSEK [24] as optimization solver. The initial
conditions for the power system are computed using power flow solutions, while
the estimator is initialized from random initial conditions having $10\%$
maximum variation from the steady state values of power system. Through out
the simulations we consider weighting constant in $\mathbf{P_{2}}$ as
$a_{1}=a_{2}=a_{3}=1$, power system base as $S_{b}=100\mathrm{MVA}$, and
$\omega_{0}=2\pi 60$ $\mathrm{{rad/sec}}$.

###  IV-A Case 1: Estimation under demand disturbance

In power systems the overall load demand is usually fluctuating and the exact
estimate of the load demand of any power network is difficult to predict. To
that end, herein we demonstrate the performance of the proposed estimator in
performing DSE under unknown disturbances in power demand. To showcase the
robustness of the estimator we provide only steady state values of loads to
the estimator and the actual transient/fluctuating values are kept unknown.

With that in mind, to generate random fluctuations in the constant power loads
we provide step disturbances at $t=3$s to the loads connected at Buses
$[7,8,12,14,15,16,18,20]$, while the loads connected at Buses $3$ and $4$ are
assumed to be changing sinusoidally with some Gaussian noise as shown in Fig.
2.

(a)  (b)

Figure 2: (a) Estimation error norm for Case 1, Case 2 and (b) Variation in
real power demand at Bus 3 and 4.

(a)  (b)  (c)  (d)

Figure 3: State estimation results for PV power plant 1 under disturbances in
real power demand.

(a)  (b)  (c)  (d)

Figure 4: Estimation results for PV power plant 2 and algebraic variables of
Bus 25 under irradiance and reactive power disturbances.

The estimation results for the PV power plant 1 are shown in Fig. 3. We notice
that although the estimator starts from initial conditions different than
system steady state values and after $t=3$s the fluctuations in load demand
are not known to the estimator, it can still accurately track the different
states of the PV power plant. This can also be corroborated from the
estimation error norm given in Fig. 2 from which we can validate that the
estimator is successfully driving the error between true and estimated values
near zero for all of the system dynamic and algebraic states, under
uncertainty. Similar performance from the estimator has been achieved in
predicting the internal states of synchronous machines and algebraic variables
of all the buses as presented in Fig. 5.

###  IV-B Case 2: Estimation under irradiance & power uncertainty

To further assess the robustness of the estimator toward different sources of
uncertainty, in this section we demonstrate the performance of the estimator
under disturbances from renewables. To that end, right after $t>0$ we
decreased the irradiance from the sun on both of the PV power plants by
$20\%$. Also, to create further unknown transient conditions we added step
disturbances in the overall reactive power demand of the system. Note that,
the estimator is completely unaware of these disturbances and only knows the
steady state value of these quantities. This can be validated from the
structure of the proposed estimator design (13). We can see that the estimator
only has access to ${\bar{\boldsymbol{w}}}$ which contains the steady
state/predicted values of load demand and irradiance.

The state estimation results are depicted in Fig. 4 and 6. For brevity
estimation results for only a few of the states are shown. We notice that the
estimator is still able to track the true values of state variables with good
accuracy and is also driving the error norm near zero as shown in Fig. 2.

(a)  (b)  (c)  (d)

Figure 5: Estimation results for Generator-6, real current at Bus 21, and
motor speed under disturbance in real power demand.

(a)  (b)

(c)  (d)

Figure 6: Estimation results for Case 2: States of Generator-1 and motor
speed.

##  V Conclusions and Future Work

In this paper, we proposed a robust estimator for an interconnected model of
power systems. The observer design is posed as a convex optimization and work
as a one step predictor. Using a few measurements provided by PMUs the
proposed estimator can provide an accurate estimate of all the states of a
power system including dynamic states of PV power plants and motor loads. The
presented estimator does not require any prior knowledge about the statistical
properties of the uncertainty and can provide accurate estimation results as
long as the uncertainty/disturbance is bounded.

The limitations of this study are twofold: First, the proposed estimator
performs DSE in a centralized fashion, thus for a very big power network it
needs to be extended to a decentralized framework. Second, the theory of the
proposed estimator is based on continuous models, however PMUs measurements
are commonly transmitted via a digitized network, thus a discrete time version
of this estimator will be more appropriate. To that end, as a future work the
proposed estimator will be extended to a discretized and decentralized
framework.

##  VI Acknowledgments

The power system model used in this study was developed by Dr. Hugo Villegas
Pico<EMAIL_ADDRESS>and Soummyar Roy<EMAIL_ADDRESS>at Iowa
State University. This work is supported by National Science Foundation under
Grants 2152450 and 2151571.

## References

  * [1] Y. Liu, A. K. Singh, J. Zhao, A. P. S. Meliopoulos, B. C. Pal, M. A. Bin Mohd Ariff, T. Van Cutsem, M. Glavic, Z. Huang, I. Kamwa, L. Mili, A. S. Mir, A. F. Taha, V. Terziya, and S. Yu, “Dynamic state estimation for power system control and protection,” _IEEE Transactions on Power Systems_ , pp. 1–1, 2021.
  * [2] S. A. Nugroho, A. F. Taha, and J. Qi, “Robust dynamic state estimation of synchronous machines with asymptotic state estimation error performance guarantees,” _IEEE Transactions on Power Systems_ , vol. 35, no. 3, pp. 1923–1935, 2020.
  * [3] S. A. Nugroho, A. Taha, N. Gatsis, and J. Zhao, “Observers for differential algebraic equation models of power networks: Jointly estimating dynamic and algebraic states,” _IEEE Transactions on Control of Network Systems_ , pp. 1–1, 2022.
  * [4] M. Jin, H. Feng, and J. Lavaei, “Multiplier-based observer design for large-scale lipschitz systems,” 2018.
  * [5] E. Ghahremani and I. Kamwa, “Dynamic state estimation in power system by applying the extended kalman filter with unknown inputs to phasor measurements,” _IEEE Transactions on Power Systems_ , vol. 26, no. 4, pp. 2556–2566, 2011.
  * [6] J. Zhao, “Dynamic state estimation with model uncertainties using $\mathcal{H}_{\infty}$ extended kalman filter,” _IEEE Transactions on Power Systems_ , vol. 33, no. 1, pp. 1099–1100, 2018.
  * [7] N. Zhou, D. Meng, and S. Lu, “Estimation of the dynamic states of synchronous machines using an extended particle filter,” _IEEE Transactions on Power Systems_ , vol. 28, no. 4, pp. 4152–4161, 2013.
  * [8] N. Zhou, D. Meng, Z. Huang, and G. Welch, “Dynamic state estimation of a synchronous machine using pmu data: A comparative study,” _IEEE Transactions on Smart Grid_ , vol. 6, no. 1, pp. 450–460, 2015.
  * [9] J. Qi, A. F. Taha, and J. Wang, “Comparing kalman filters and observers for power system dynamic state estimation with model uncertainty and malicious cyber attacks,” _IEEE Access_ , vol. 6, pp. 77 155–77 168, 2018.
  * [10] Z. Fang, Y. Lin, S. Song, C. Li, X. Lin, and Y. Chen, “State estimation for situational awareness of active distribution system with photovoltaic power plants,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 1, pp. 239–250, 2021.
  * [11] P. Sauer, M. Pai, and J. Chow, _Power System Dynamics and Stability: With Synchrophasor Measurement and Power System Toolbox_ , ser. Wiley - IEEE. Wiley, 2017.
  * [12] H. N. V. Pico and B. B. Johnson, “Transient stability assessment of multi-machine multi-converter power systems,” _IEEE Transactions on Power Systems_ , vol. 34, no. 5, pp. 3504–3514, 2019.
  * [13] P. C. Krause, O. Wasynczuk, S. D. Sudhoff, and S. D. Pekarek, _Analysis of Electric Machinery and Drive Systems_ , 3rd ed. John Wiley & Sons, 2013.
  * [14] CAISO, “Real time demand curves,” Nov 2021. [Online]. Available: http://www.caiso.com/outlook/SystemStatus.html
  * [15] S. A. Nugroho, A. F. Taha, and J. Qi, “Characterizing the nonlinearity of power system generator models,” in _2019 American Control Conference (ACC)_ , 2019, pp. 1936–1941.
  * [16] S. A. Nugroho, A. Taha, and V. Hoang, “Nonlinear dynamic systems parameterization using interval-based global optimization: Computing lipschitz constants and beyond,” _IEEE Transactions on Automatic Control_ , pp. 1–1, 2021.
  * [17] U. Shaked, “$\mathcal{H}_{\infty}$ minimum error state estimation of linear stationary processes,” _IEEE Transactions on Automatic Control_ , vol. 35, no. 5, pp. 554–558, 1990.
  * [18] K. Derinkuyu and M. Ç. Pınar, “On the s-procedure and some variants,” _Mathematical Methods of Operations Research_ , vol. 64, no. 1, pp. 55–77, 2006.
  * [19] J. Zhao, Z. Zheng, S. Wang, R. Huang, T. Bi, L. Mili, and Z. Huang, “Correlation-aided robust decentralized dynamic state estimation of power systems with unknown control inputs,” _IEEE Transactions on Power Systems_ , vol. 35, no. 3, pp. 2443–2451, 2020.
  * [20] L. Shengyuan, Xu andJames, _Robust Control and Filtering of Singular Systems, ser. Lecture Notes in Control and Information Sciences._ Springer-Verlag Berlin/Heidelberg., 2004.
  * [21] I. Hiskens, _IEEE PES task force on benchmark systems for stability controls_. IEEE PES, Piscataway, NJ, USA, Tech. Rep., IEEE PES, Piscataway, NJ, USA, Tech. Rep, Nov. 2013.
  * [22] S. Chakrabarti and E. Kyriakides, “Optimal placement of phasor measurement units for power system observability,” _IEEE Transactions on Power Systems_ , vol. 23, no. 3, pp. 1433–1440, 2008.
  * [23] J. Lofberg, “Yalmip : a toolbox for modeling and optimization in matlab,” in _2004 IEEE International Conference on Robotics and Automation (IEEE Cat. No.04CH37508)_ , 2004, pp. 284–289.
  * [24] E. D. Andersen and K. D. Andersen, _The Mosek Interior Point Optimizer for Linear Programming: An Implementation of the Homogeneous Algorithm_. Boston, MA: Springer US, 2000, pp. 197–232.