# Dynamic State Estimation for Radial Microgrid Protection

Arthur K. Barnes 1*  — and Adam Mate 1*  — Manuscript submitted: December 14,
2020. 1 The authors are with the Advanced Network Science Initiative at Los
Alamos National Laboratory, Los Alamos, NM 87544 USA. Email:{abarnes,
<EMAIL_ADDRESS>versions of one or more of the figures in this paper are
available online at https://ieeexplore.ieee.org.LANL ANSI LA-UR-20-30126.

###### Abstract

Microgrids are localized electrical grids with control capability that are
able to disconnect from the traditional grid to operate autonomously. They
strengthen grid resilience, help mitigate grid disturbances, and support a
flexible grid by enabling the integration of distributed energy resources.
Given the likely presence of critical loads, the proper protection of
microgrids is of vital importance; however, this is complicated in the case of
inverter-interfaced microgrids where low fault currents preclude the use of
conventional time-overcurrent protection. This paper introduces and
investigates the application of dynamic state estimation, a generalization of
differential protection, for the protection of radial portions of microgrids
(or distribution networks); both phasor-based and dynamic approaches are
investigated for protection. It is demonstrated through experiments on three
case-study systems that dynamic state estimation is capable of correctly
identifying model parameters for both normal and faulted operation.

###### Index Terms:

power system operation, dynamic state estimation, microgrid, distribution
network, protection.

##  I Introduction

Dynamic state estimation (abbr. DSE) is a generalization of differential
protection, which offers a reduced likelihood of misoperation, particularly in
the case of assets with nonlinear characteristics (e.g., transformers that are
being energized) [1]. It is also useful in cases where distance protection
performs poorly (e.g., transmission lines with series compensation [2] or
mutually coupled transmission lines [3]). DSE has previously been introduced
to microgrid branch protection [4, 5, 6].

This paper investigates the application of DSE for the protection of radial
portions of a microgrid (or a distribution network). This can be a challenge
in electrical grids with distributed generation on account of lack of fault
current from inverter-interfaced generation [7], varying fault current between
grid-connected and islanded modes [7], and the potential for normally-meshed
operation [8] and unbalanced operation due to single-phase loads [8].
Admittance relaying has been investigated as a solution for the protection of
microgrids [9], however, it has been observed to present issues with grounded-
wye connected loads [10]; consequently, additional relaying is necessary to
prevent misoperation [8].

This paper treats the radial portions of a microgrid as load busses. It is
assumed that these portions contain no loops or downstream generation; they
are modeled as constant-impedance networks with unknown impedances but known
connectivity. To ensure that the number of measured variables is greater than
approximately 1.6 times the number of free parameters (where 1.6 is a commonly
selected number for redundancy to ensure sufficient measurements for system
identification [11, 12]), most models presented here make the assumption that
the loads are balanced.

Every load and fault configuration requires a separate model. For a given load
configuration, a model for each fault configuration is fit to measured values;
the model with the lowest error, in terms of fitting the observed variables,
is assumed to be the correct one. On a grounded-wye-connected load, the
following models are necessary to distinguish between normal operation, line-
ground faults and line-line faults:

  1. 1.

Normal operation: each branch of the load has the same impedance, which is
modeled as a series resistive-inductive (abbr. RL) network.

  2. 2.

Phase A-ground fault: the faulted branch A is modeled as a resistance, while
the unfaulted branches B and C are modeled as series RL networks with equal
parameters.

  3. 3.

Phase B-ground fault: the faulted branch B is modeled as a resistance, while
the unfaulted branches C and A are modeled as series RL networks with equal
parameters.

  4. 4.

Phase C-ground fault: the faulted branch C is modeled as a resistance, while
the unfaulted branches A and B are modeled as series RL networks with equal
parameters.

  5. 5.

Phase A-B fault: the fault impedance is modeled as a resistance across the
load terminals A and B, while each branch of the load is modeled as a series
RL network.

  6. 6.

Phase B-C fault: the fault impedance is modeled as a resistance across the
load terminals B and C, while each branch of the load is modeled as a series
RL network.

  7. 7.

Phase C-A fault: the fault impedance is modeled as a resistance across the
load terminals C and A, while each branch of the load is modeled as a series
RL network.

On a delta-connected system, the following models are necessary to distinguish
between normal operation, line-ground and line-line faults:

  1. 1.

Normal operation: each branch of the load has the same impedance, which is
modeled as series RL network.

  2. 2.

Phase A-ground fault: the fault impedance is modeled as a resistance between
load terminal A and ground, while the load branches are modeled as series RL
networks.

  3. 3.

Phase B-ground fault: the fault impedance is modeled as a resistance between
load terminal B and ground, while the load branches are modeled as series RL
networks.

  4. 4.

Phase C-ground fault: the fault impedance is modeled as a resistance between
load terminal C and ground, while the load branches are modeled as series RL
networks.

  5. 5.

Phase A-B fault: the fault impedance is modeled as a resistance across the
load terminals A and B, while the branches across load terminals B-C and C-A
are modeled as series RL networks.

  6. 6.

Phase B-C fault: the fault impedance is modeled as a resistance across the
load terminals B and C, while the branches across load terminals C-A and A-B
are modeled as series RL networks.

  7. 7.

Phase C-A fault: the fault impedance is modeled as a resistance across the
load terminals C and A, while the branches across load terminals A-B and B-C
are modeled as series RL networks.

Both phasor-based and dynamic approaches are investigated for radial microgrid
protection. Section II describes the implementation of phasor-based state
estimation: it is conceptually similar to DSE but more straightforward to
derive and implement as it only requires a single time period. Section III
describes the implementation of DSE. Section IV describes how two different
transient models of loads are developed as test cases and run to test both
phasor and dynamic state estimation; next, Section V presents the performance
of state estimation on the test cases. Finally, Section VI summarizes the
conclusions of this paper.

##  II Phasor Implementation

The phasor implementation of state estimation-based protection is simpler:
only a single time period is used, which limits the number of measurements and
therefore the number of parameters that can be estimated. In this section, it
is applied to single-phase, grounded-wye and delta-connected load
configurations.

###  II-A Single-Phase Impedance

The output of the system (illustrated in Fig. 1a):

$\mathbf{y}=\begin{bmatrix}V\\\ I\end{bmatrix}$

where $V$ and $I$ are phasor quantities.

The state of the system:

$\mathbf{x}=\begin{bmatrix}Z\\\ I_{z}\end{bmatrix}$

The output-state mapping for the system is the following vector-valued
function:

$\mathbf{y}=\mathbf{h}(\mathbf{x})$

where

$\displaystyle h_{1}(\mathbf{x})$ $\displaystyle=V_{z}=ZI_{z}$ $\displaystyle
h_{2}(\mathbf{x})$ $\displaystyle=I_{z}$

The Jacobian of $\mathbf{h}(\mathbf{x})$ is determined as follows:

$\displaystyle\frac{\partial V_{z}}{\partial Z}$
$\displaystyle=\frac{\partial}{\partial Z}ZI_{z}=I_{z}$
$\displaystyle\frac{\partial V_{z}}{\partial I_{z}}$
$\displaystyle=\frac{\partial}{\partial I_{z}}ZI_{z}=Z$
$\displaystyle\frac{\partial I_{z}}{\partial Z}$
$\displaystyle=\frac{\partial}{\partial Z}I_{z}=0$
$\displaystyle\frac{\partial I_{z}}{\partial I_{z}}$
$\displaystyle=\frac{\partial}{\partial I_{z}}I_{z}=1$

The mapping between variables and the state vector:

$\displaystyle Z$ $\displaystyle=x_{1}$ $\displaystyle I_{z}$
$\displaystyle=x_{2}$

Given the variable and state mapping, the Jacobian can be built as follows:
$H(n,m)=0$, unless specified below.

$H=\begin{bmatrix}[1.5]\frac{\partial V_{z}}{\partial Z}&\frac{\partial
V_{z}}{\partial I_{z}}\\\ \frac{\partial I_{z}}{\partial Z}&\frac{\partial
I_{z}}{\partial I_{z}}\end{bmatrix}$

Given the Jacobian, the state of the system can be solved for iteratively:

$\displaystyle{\bm{\epsilon}}_{i}$
$\displaystyle=\mathbf{y}-\mathbf{h}(\mathbf{x}_{i})\hskip 36.135pt$
$\displaystyle J_{i}$ $\displaystyle=||{\bm{\epsilon}}_{i}||^{2}$
$\displaystyle\mathbf{x}_{i+1}$
$\displaystyle=\mathbf{x}_{i}+(H^{\prime}_{i}H_{i})^{-1}H^{\prime}_{i}{\bm{\epsilon}}_{i}$

###  II-B Grounded-Wye with Line-Ground Fault

The output of the system (illustrated in Fig. 1b):

$\mathbf{y}=\begin{bmatrix}I_{a}&I_{b}&I_{c}&V_{a}&V_{b}&V_{c}\end{bmatrix}^{T}$

The easiest way to model this is as an unbalanced load where the fault
impedance is not treated specially. The state of the system therefore:

$\mathbf{x}=\begin{bmatrix}Y_{a}&Y_{b}&Y_{c}&V_{za}&V_{zb}&V_{zc}\end{bmatrix}^{T}$

The output function $\mathbf{h}(\mathbf{x})$ can be written as:

$\displaystyle\mathbf{h}_{1}(\mathbf{x})$ $\displaystyle=I_{a}=y_{a}V_{a}$
$\displaystyle\mathbf{h}_{2}(\mathbf{x})$ $\displaystyle=I_{b}=y_{b}V_{b}$
$\displaystyle\mathbf{h}_{3}(\mathbf{x})$ $\displaystyle=I_{c}=y_{c}V_{c}$
$\displaystyle\mathbf{h}_{4}(\mathbf{x})$ $\displaystyle=V_{a}=V_{za}$
$\displaystyle\mathbf{h}_{5}(\mathbf{x})$ $\displaystyle=V_{b}=V_{zb}$
$\displaystyle\mathbf{h}_{6}(\mathbf{x})$ $\displaystyle=V_{c}=V_{zc}$

###  II-C Grounded-Wye with Line-Line Fault

The output of the system (illustrated in Fig. 1c):

$\mathbf{y}=\begin{bmatrix}I_{a}&I_{b}&I_{c}&V_{a}&V_{b}&V_{c}\end{bmatrix}^{T}$

The state of the system:

$\mathbf{x}=\begin{bmatrix}Y_{l}&Y_{f}&V_{za}&V_{zb}&V_{zc}\end{bmatrix}^{T}$

The output function $\mathbf{h}(\mathbf{x})$ can be written as:

$\displaystyle h_{1}(\mathbf{x})$
$\displaystyle=I_{a}=(y_{l}+y_{f})V_{za}-y_{f}V_{zb}$ $\displaystyle
h_{2}(\mathbf{x})$ $\displaystyle=I_{b}=-y_{f}V_{za}+y_{l}V_{zb}$
$\displaystyle h_{3}(\mathbf{x})$ $\displaystyle=I_{c}=y_{l}V_{zc}$
$\displaystyle h_{4}(\mathbf{x})$ $\displaystyle=V_{a}=V_{za}$ $\displaystyle
h_{5}(\mathbf{x})$ $\displaystyle=V_{b}=V_{zb}$ $\displaystyle
h_{6}(\mathbf{x})$ $\displaystyle=V_{c}=V_{zc}$

###  II-D Delta-Connected Load with Line-Line Fault

The output of the system (illustrated in Fig. 1d):

$\mathbf{y}=\begin{bmatrix}I_{a}&I_{b}&I_{c}&V_{a}&V_{b}&V_{c}\end{bmatrix}^{T}$

The state of the system:

$\mathbf{x}=\begin{bmatrix}Y_{f}&Y_{ll}&V_{za}&V_{zb}&V_{zc}\end{bmatrix}^{T}$

The output function $\mathbf{h}(\mathbf{x})$ can be written as:

$\displaystyle h_{1}(\mathbf{x})$
$\displaystyle=I_{a}=y_{aa}V_{za}-y_{ab}V_{zb}-y_{ca}V_{zc}$
$\displaystyle=(y_{ab}+y_{ca})V_{za}-y_{ab}V_{zb}-y_{ca}V_{zc}$
$\displaystyle=(y_{f}+y_{ll})V_{za}-y_{f}V_{sb}-y_{ll}V_{zc}$ $\displaystyle
h_{2}(\mathbf{x})$
$\displaystyle=I_{b}=-y_{ab}V_{za}+y_{bb}V_{zb}-y_{bc}V_{zc}$
$\displaystyle=-y_{ab}V_{za}+(y_{ab}+y_{bc})V_{zb}-y_{bc}V_{zc}$
$\displaystyle=-y_{f}V_{za}+(y_{f}+y_{ll})V_{zb}-y_{ll}V_{zc}$ $\displaystyle
h_{3}(\mathbf{x})$
$\displaystyle=I_{c}=-y_{ca}V_{za}-y_{bc}V_{zb}+y_{cc}V_{zc}$
$\displaystyle=-y_{ca}V_{za}-y_{bc}V_{zb}+(y_{ac}+y_{bc})V_{zc}$
$\displaystyle=-y_{ll}V_{za}-y_{ll}V_{zb}+2y_{ll}V_{zc}$ $\displaystyle
h_{4}(\mathbf{x})$ $\displaystyle=V_{a}=V_{za}$ $\displaystyle
h_{5}(\mathbf{x})$ $\displaystyle=V_{b}=V_{zb}$ $\displaystyle
h_{6}(\mathbf{x})$ $\displaystyle=V_{c}=V_{zc}$

###  II-E Delta-Connected Load with a Line-Ground Fault

The output of the system (illustrated in Fig. 1e):

$\mathbf{y}=\begin{bmatrix}I_{a}&I_{b}&I_{c}&V_{a}&V_{b}&V_{c}\end{bmatrix}^{T}$

The state of the system:

$\mathbf{x}=\begin{bmatrix}Y_{ll}&Y_{f}&V_{za}&V_{zb}&V_{zc}\end{bmatrix}^{T}$

The output function $\mathbf{h}(\mathbf{x})$ can be written as:

$\displaystyle h_{1}(\mathbf{x})$
$\displaystyle=I_{a}=y_{aa}V_{za}-y_{ab}V_{zb}-y_{ca}V_{zc}$
$\displaystyle=(y_{ab}+y_{ca}+y_{ag})V_{za}-y_{ab}V_{zb}-y_{ca}V_{zc}$
$\displaystyle=(y_{f}+2y_{ll})V_{za}-y_{ll}V_{sb}-y_{ll}V_{zc}$ $\displaystyle
h_{2}(\mathbf{x})$
$\displaystyle=I_{b}=-y_{ab}V_{za}+y_{bb}V_{zb}-y_{bc}V_{zc}$
$\displaystyle=-y_{ab}V_{za}+(y_{ab}+y_{bc})V_{zb}-y_{bc}V_{zc}$
$\displaystyle=-y_{ll}V_{za}+2y_{ll}V_{zb}-y_{ll}V_{zc}$ $\displaystyle
h_{3}(\mathbf{x})$
$\displaystyle=I_{c}=-y_{ca}V_{za}-y_{bc}V_{zb}+y_{cc}V_{zc}$
$\displaystyle=-y_{ca}V_{za}-y_{bc}V_{zb}+(y_{ac}+y_{bc})V_{zc}$
$\displaystyle=-y_{ll}V_{za}-y_{ll}V_{zb}+2y_{ll}V_{zc}$ $\displaystyle
h_{4}(\mathbf{x})$ $\displaystyle=V_{a}=V_{za}$ $\displaystyle
h_{5}(\mathbf{x})$ $\displaystyle=V_{b}=V_{zb}$ $\displaystyle
h_{6}(\mathbf{x})$ $\displaystyle=V_{c}=V_{zc}$

(a) Single-phase load

(b) Grounded-wye load with line-ground fault fault

(c) Grounded-wye load with line-line fault

(d) Delta-connected load with a line-line fault

(e) Delta-connected load with a line-ground fault

Figure 1: Phasor-based implementation of state estimation-based protection,
applied to single-phase, grounded-wye and delta-connected load configurations.

##  III Dynamic Implementation

While the phasor implementation uses a single time period for state
estimation, with the dynamic implementation several periods are used; in this
paper, 12 cycles are sampled at a 2 [kHz] sample rate. As in Section II, here
the DSE-based protection is applied to single-phase, grounded-wye and delta-
connected load configurations.

###  III-A Single-Phase Series RL Load

The output of the system (illustrated in Fig. 2a):

$y(t)=\begin{bmatrix}v(t)\\\ i(t)z(t)\end{bmatrix}$

For the purposes of state estimation, this is sampled at points
$n\in\\{1,...,N\\}$ giving the following vector-value equation:

$\mathbf{y}=\begin{bmatrix}\mathbf{v}\\\ \mathbf{i}\end{bmatrix}$

where

$\displaystyle\mathbf{v}$
$\displaystyle=\begin{bmatrix}v(1)&v(2)&\cdots&v(N)\end{bmatrix}^{T}$
$\displaystyle\mathbf{i}$
$\displaystyle=\begin{bmatrix}i(1)&i(2)&\cdots&i(N)\end{bmatrix}^{T}$
$\displaystyle\mathbf{z}$
$\displaystyle=\begin{bmatrix}z(1)&z(2)&\cdots&z(N)\end{bmatrix}^{T}$

The state of the system:

$\mathbf{x}=\begin{bmatrix}R&L&\mathbf{v}_{r}&\mathbf{v}_{l}\end{bmatrix}^{T}$

The output-state mapping for the system is the following vector-valued
function:

$\mathbf{y}=\mathbf{h}(\mathbf{x})$

where

$\displaystyle h_{n}(\mathbf{x})=v_{r}(n)+v_{l}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{N+n}(\mathbf{x})=Gv_{r}(n)\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{2N+n}(\mathbf{x})=Gv_{r}(n)-Gv_{r}(n-2)+\frac{2\Lambda\Delta
t}{6}(v_{l}(n)+$ $\displaystyle+4v_{l}(n-1)+v_{l}(n-2)),\quad\forall
n\in\\{3,4,\ldots,N\\}$

In the above, $v_{R}(n)=Ri_{L}(n)$ follows from discretizing
$v_{R}(t)=Ri_{L}(t)$, and

$v_{l}(n)=\frac{2\Lambda\Delta t}{6}(v_{l}(n)+4v_{l}(n-1)+v_{l}(n-2))$

follows from discretizing

$i_{l}(t)=\frac{1}{L}\int_{t-\Delta t}^{t}v_{l}(\tau)d\tau$

via Simpson’s 1/3 rule [13].

Given the variable and state vector mapping, the Jacobian can be built as
follows: $H(n,m)=0$, unless specified otherwise below.

$\displaystyle H(n,2+n)=\frac{\partial v(n)}{\partial v_{r}(n)}=1\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle H(n,2+N+n)=\frac{\partial
v(n)}{\partial v_{l}(n)}=1\quad\forall n\in\\{1,2,\ldots,N\\}$ $\displaystyle
H(N+n,1)=\frac{\partial i(n)}{\partial G}=v_{r}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle H(N+n,2+n)=\frac{\partial
i(n)}{\partial v_{r}(n)}=G\quad\forall n\in\\{1,2,\ldots,N\\}$ $\displaystyle
H(2N+n-2,1)=\frac{\partial z(n-2)}{\partial R}=v_{r}(n)-v_{r}(n-2)$
$\displaystyle\quad\forall n\in\\{3,4,\ldots,N\\}$ $\displaystyle
H(2N+n-2,n)=\frac{\partial z(n-2)}{\partial v_{r}(n)}=G\quad\forall
n\in\\{3,4,\ldots,N\\}$ $\displaystyle H(2N+n-2,2)=\frac{\partial
z(n-2)}{\partial\Lambda}=\frac{2\Delta t}{6}(v_{l}(n)+$
$\displaystyle+4v_{l}(n-1)+v_{l}(n-2))\quad\forall n\in\\{3,4,\ldots,N\\}$
$\displaystyle H(2N+n-2,2+N+n)=\frac{\partial z(n-2)}{\partial
v_{l}(n)}=\frac{\Delta t\Lambda}{3}$ $\displaystyle\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle H(2N+n-2,1+N+n)=\frac{\partial
z(n-2)}{\partial v_{l}(n-1)}=\frac{4\Delta t\Lambda}{3}$
$\displaystyle\quad\forall n\in\\{1,2,\ldots,N\\}$ $\displaystyle
H(2N+n-2,N+n)=\frac{\partial z(n-2)}{\partial v_{l}(n-2)}=\frac{\Delta
t\Lambda}{3}$ $\displaystyle\quad\forall n\in\\{1,2,\ldots,N\\}$

Given the Jacobian, the state of the system can be solved for iteratively, by
applying the same equations as in Section II-A.

###  III-B Grounded-Wye Load without Fault

The sampled output of the system (illustrated in Fig. 2b):

$\mathbf{y}=\begin{bmatrix}\mathbf{v}_{a}&\mathbf{v}_{b}&\mathbf{v}_{c}&\mathbf{i}_{a}&\mathbf{i}_{b}&\mathbf{i}_{c}&\mathbf{z}_{a}&\mathbf{z}_{b}&\mathbf{z}_{c}\end{bmatrix}^{T}$

where

$\displaystyle\mathbf{v}_{\phi}$
$\displaystyle=\begin{bmatrix}v_{\phi}(1)&v_{\phi}(2)&\cdots&v_{\phi}(N)\end{bmatrix}^{T}$
$\displaystyle\mathbf{i}_{\phi}$
$\displaystyle=\begin{bmatrix}i_{\phi}(1)&i_{\phi}(2)&\cdots&i_{\phi}(N)\end{bmatrix}^{T}$
$\displaystyle\mathbf{z}_{\phi}$
$\displaystyle=\begin{bmatrix}z_{\phi}(1)&z_{\phi}(2)&\cdots&z_{\phi}(N-2)\end{bmatrix}^{T}for\hskip
3.61371pt\phi\in\\{a,b,c\\}$

The state of the system:

$x(t)=\begin{bmatrix}G&\Lambda&\mathbf{v}_{ra}&\mathbf{v}_{rb}&\mathbf{v}_{rc}&\mathbf{v}_{la}&\mathbf{v}_{lb}&\mathbf{v}_{lc}\end{bmatrix}^{T}$

where $G=R^{-1}$ is the conductance, $\Lambda=L^{-1}$ is the reciprocal of the
inductance, $\mathbf{v}_{r\phi}$ is the voltage across the resistance on phase
$\phi$ at each time period $1,\ldots,N$ and $\mathbf{v}_{l\phi}$ is the
voltage across the inductance on phase $\phi$ at each time period
$1,\ldots,N$.

The output function $\mathbf{h}(\mathbf{x})$ can be written as:

$\displaystyle h_{n}(\mathbf{x})=v_{ra}(n)+v_{la}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+N}(\mathbf{x})=v_{rb}(n)+v_{lb}(n)\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+2N}(\mathbf{x})=v_{rb}(n)+v_{lb}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+3N}(\mathbf{x})=Gv_{ra}(n)\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+4N}(\mathbf{x})=Gv_{rb}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+5N}(\mathbf{x})=Gv_{rc}(n))\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+6N}(\mathbf{x})=G(v_{ra}(n)-v_{ra}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{la}(n)+$ $\displaystyle+4v_{la}(n-1)+v_{la}(n-2))\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+7N}(\mathbf{x})=G(v_{rb}(n)-v_{rb}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lb}(n)+$ $\displaystyle+4v_{lc}(n-1)+v_{lb}(n-2))\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+8N}(\mathbf{x})=G(v_{rc}(n)-v_{rc}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lc}(n)+$ $\displaystyle+4v_{lc}(n-1)+v_{lc}(n-2))\quad\forall
n\in\\{1,2,\ldots,N\\}.$

###  III-C Grounded-Wye Load with Line-Ground Fault

The sampled output of the system (illustrated in Fig. 2c):

$\mathbf{y}=\begin{bmatrix}\mathbf{v}_{a}&\mathbf{v}_{b}&\mathbf{v}_{c}&\mathbf{i}_{a}&\mathbf{i}_{b}&\mathbf{i}_{c}&\mathbf{z}_{b}&\mathbf{z}_{c}\end{bmatrix}^{T}$

Note that there are no $z_{a}(n)$ output variables as the reactive impedance
on phase A is large compared to the parallel fault conductance $G_{f}$.

The state of the system:

$x(t)=\begin{bmatrix}G&\Lambda&G_{f}&\mathbf{v}_{ra}&\mathbf{v}_{rb}&\mathbf{v}_{rc}&\mathbf{v}_{lb}&\mathbf{v}_{lc}\end{bmatrix}^{T}$

where $G_{f}=R^{-1}$ is the conductance and the remaining states are the same
as those in that of the grounded-wye no-fault state.

The output function $\mathbf{h}(\mathbf{x})$ can be written as:

$\displaystyle h_{n}(\mathbf{x})=v_{ra}(n)\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+N}(\mathbf{x})=v_{rb}(n)+v_{lb}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+2N}(\mathbf{x})=v_{rb}(n)+v_{lb}(n)\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+3N}(\mathbf{x})=Gv_{ra}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+4N}(\mathbf{x})=Gv_{rb}(n)\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+5N}(\mathbf{x})=Gv_{rc}(n))\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+6N}(\mathbf{x})=G(v_{rb}(n)-v_{rb}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lb}(n)+$ $\displaystyle+4v_{lc}(n-1)+v_{lb}(n-2))\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+7N}(\mathbf{x})=G(v_{rc}(n)-v_{rc}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lc}(n)+$ $\displaystyle+4v_{lc}(n-1)+v_{lc}(n-2))\quad\forall
n\in\\{1,2,\ldots,N\\}$

###  III-D Grounded-Wye Load with Line-Line Fault

The sampled output of the system (illustrated in Fig. 2d):

$\mathbf{y}=\begin{bmatrix}\mathbf{v}_{a}&\mathbf{v}_{b}&\mathbf{v}_{c}&\mathbf{i}_{a}&\mathbf{i}_{b}&\mathbf{i}_{c}&\mathbf{z}_{a}&\mathbf{z}_{b}&\mathbf{z}_{c}\end{bmatrix}^{T}$

The state of the system:

$x(t)=\begin{bmatrix}G&\Lambda&\mathbf{v}_{ra}&\mathbf{v}_{rb}&\mathbf{v}_{rc}&\mathbf{v}_{la}&\mathbf{v}_{lb}&\mathbf{v}_{lc}\end{bmatrix}^{T}$

The output function $\mathbf{h}(\mathbf{x})$ can be written as:

$\displaystyle h_{n}(\mathbf{x})=v_{ra}(n)+v_{la}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+N}(\mathbf{x})=v_{rb}(n)+v_{lb}(n)\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+2N}(\mathbf{x})=v_{rb}(n)+v_{lb}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+3N}(\mathbf{x})=Gv_{ra}(n)+G_{f}(v_{ra}(n)+v_{la}(n)-v_{rb}(n)-$
$\displaystyle-v_{lb}(n))\quad\forall n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+4N}(\mathbf{x})=Gv_{rb}(n)-G_{f}(v_{ra}(n)+v_{la}(n)-v_{rb}(n)-$
$\displaystyle-v_{lb}(n))\quad\forall n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+5N}(\mathbf{x})=Gv_{rc}(n))\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+6N}(\mathbf{x})=G(v_{ra}(n)-v_{ra}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{la}(n)+$ $\displaystyle+4v_{la}(n-1)+v_{la}(n-2))\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+7N}(\mathbf{x})=G(v_{rb}(n)-v_{rb}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lb}(n)+$ $\displaystyle+4v_{lc}(n-1)+v_{lb}(n-2))\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+8N}(\mathbf{x})=G(v_{rc}(n)-v_{rc}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lc}(n)+$ $\displaystyle+4v_{lc}(n-1)+v_{lc}(n-2))\quad\forall
n\in\\{1,2,\ldots,N\\}$

###  III-E Delta Load without Fault

The sampled output of the system (illustrated in Fig. 2f):

$\mathbf{y}=\begin{bmatrix}\mathbf{v}_{ab}&\mathbf{v}_{bc}&\mathbf{v}_{ca}&\mathbf{i}_{a}&\mathbf{i}_{b}&\mathbf{i}_{c},&\mathbf{z}_{ab}&\mathbf{z}_{bc}&\mathbf{z}_{ca}\end{bmatrix}^{T}$

The state of the system:

$x(t)=\begin{bmatrix}G&\Lambda\mathbf{v}_{rab}&\mathbf{v}_{rbc}&\mathbf{v}_{rca}&\mathbf{v}_{lab}&\mathbf{v}_{lbc}&\mathbf{v}_{lca}\end{bmatrix}^{T}$

The output function $\mathbf{h}(\mathbf{x})$ can be written as:

$\displaystyle h_{n}(\mathbf{x})=v_{rab}(n)+v_{lab}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+N}(\mathbf{x})=v_{rbc}(n)+v_{lbc}(n)\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+2N}(\mathbf{x})=v_{rca}(n)+v_{lca}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+3N}(\mathbf{x})=G(v_{rab}(n)-v_{rca}(n))\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+4N}(\mathbf{x})=G(v_{rbc}(n)-v_{rab}(n))\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+5N}(\mathbf{x})=G(v_{rca}(n)-v_{rbc}(n))\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+6N}(\mathbf{x})=G(v_{rab}(n)-v_{rab}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lab}(n)+$
$\displaystyle+4v_{lab}(n-1)+v_{lab}(n-2))\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+7N}(\mathbf{x})=G(v_{rbc}(n)-v_{rbc}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lbc}(n)+$
$\displaystyle+4v_{lbc}(n-1)+v_{lbc}(n-2))\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+8N}(\mathbf{x})=G(v_{rca}(n)-v_{rca}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lca}(n)+$
$\displaystyle+4v_{lca}(n-1)+v_{lca}(n-2))\quad\forall n\in\\{1,2,\ldots,N\\}$

###  III-F Delta Load with Line-Line Fault

The sampled output of the system (illustrated in Fig. 2d):

$\mathbf{y}=\begin{bmatrix}\mathbf{v}_{ab}&\mathbf{v}_{bc}&\mathbf{v}_{ca}&\mathbf{i}_{a}&\mathbf{i}_{b}&\mathbf{i}_{c},&\mathbf{z}_{bc}&\mathbf{z}_{ca}\end{bmatrix}^{T}$

Note that there are no $z_{ab}(n)$ output variables as the reactive impedance
on phase A is large compared to the parallel fault conductance $G_{f}$.

The state of the system:

$x(t)=\begin{bmatrix}G&\Lambda&G_{f}&\mathbf{v}_{rab}&\mathbf{v}_{rbc}&\mathbf{v}_{rca}&\mathbf{v}_{lbc}&\mathbf{v}_{lca}\end{bmatrix}^{T}$

where $G_{f}=R^{-1}$ is the conductance and the remaining states are the same
as those in that of the delta no-fault state.

The output function $\mathbf{h}(\mathbf{x})$ can be written as:

$\displaystyle h_{n}(\mathbf{x})=v_{rab}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+N}(\mathbf{x})=v_{rbc}(n)+v_{lbc}(n)\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+2N}(\mathbf{x})=v_{rca}(n)+v_{lca}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+3N}(\mathbf{x})=Gv_{ra}(n)\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+4N}(\mathbf{x})=Gv_{rb}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+5N}(\mathbf{x})=Gv_{rc}(n))\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+6N}(\mathbf{x})=G(v_{rbc}(n)-v_{rbc}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lbc}(n)+$
$\displaystyle+4v_{lbc}(n-1)+v_{lbc}(n-2))\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+7N}(\mathbf{x})=G(v_{rca}(n)-v_{rca}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lca}(n)+$
$\displaystyle+4v_{lca}(n-1)+v_{lca}(n-2))\quad\forall n\in\\{1,2,\ldots,N\\}$

###  III-G Delta Load with Line-Ground Fault

The sampled output of the system (illustrated in Fig. 2g):

$\mathbf{y}=\begin{bmatrix}\mathbf{v}_{ab}&\mathbf{v}_{bc}&\mathbf{v}_{ca}&\mathbf{v}_{a}&\mathbf{i}_{a}&\mathbf{i}_{b}&\mathbf{i}_{c},&\mathbf{z}_{ab}&\mathbf{z}_{bc}&\mathbf{z}_{ca}\end{bmatrix}^{T}$

The state of the system:

$\displaystyle
x(t)=\begin{bmatrix}G&\Lambda&G_{f}&\mathbf{v}_{rab}&\mathbf{v}_{rbc}&\mathbf{v}_{rca}&...\end{bmatrix}^{T}$
$\displaystyle\begin{bmatrix}...&\mathbf{v}_{lab}&\mathbf{v}_{lbc}&\mathbf{v}_{lca}&\mathbf{v}_{f}\end{bmatrix}^{T}$

where $G_{f}=R^{-1}$ is the fault conductance, $\mathbf{v}_{f}$ is the voltage
across the fault and the remaining states are the same as those in that of the
delta no-fault state.

The output function $\mathbf{h}(\mathbf{x})$ can be written as:

$\displaystyle h_{n}(\mathbf{x})=v_{rab}(n)+v_{lab}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+N}(\mathbf{x})=v_{rbc}(n)+v_{lbc}(n)\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+2N}(\mathbf{x})=v_{rca}(n)+v_{lca}(n)\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+3N}(\mathbf{x})=v_{f}(n)\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+4N}(\mathbf{x})=G(v_{rab}(n)-v_{rca}(n))+G_{f}v_{f}(n)$
$\displaystyle\quad\forall n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+5N}(\mathbf{x})=G(v_{rbc}(n)-v_{rab}(n))\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+6N}(\mathbf{x})=G(v_{rca}(n)-v_{rbc}(n))\quad\forall
n\in\\{1,2,\ldots,N\\}$ $\displaystyle
h_{n+7N}(\mathbf{x})=G(v_{rab}(n)-v_{rab}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lab}(n)+$
$\displaystyle+4v_{lab}(n-1)+v_{lab}(n-2))\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+8N}(\mathbf{x})=G(v_{rbc}(n)-v_{rbc}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lbc}(n)+$
$\displaystyle+4v_{lbc}(n-1)+v_{lbc}(n-2))\quad\forall n\in\\{1,2,\ldots,N\\}$
$\displaystyle h_{n+9N}(\mathbf{x})=G(v_{rca}(n)-v_{rca}(n-2))-\frac{2\Delta
t\Lambda}{6}(v_{lca}(n)+$
$\displaystyle+4v_{lca}(n-1)+v_{lca}(n-2))\quad\forall n\in\\{1,2,\ldots,N\\}$

(a) Single-phase RL series load

(b) Grounded-wye-connected RL load

(c) Grounded-wye-connected RL load with a line-ground fault

(d) Grounded-wye-connected RL load with a line-line fault

(e) Delta-connected RL load

(f) Delta-connected RL load with a line-line fault

(g) Delta-connected RL load with a line-ground fault

Figure 2: Dynamic implementation of state estimation-based protection, applied
to single-phase, grounded-wye and delta-connected load configurations.

##  IV Experiments

Three different case-study systems are considered: 1) a single-phase load, 2)
a grounded-wye constant-impedance load, and 3) a delta-connected constant
impedance load. These load configurations are studied for both phasor-based
state estimation and DSE. To verify the noise immunity of the methods, random
noise with an amplitude of approximately 10% of the signal peak is added to
the measurements in both cases. Experiments are performed in Julia v1.5 [14]
and 64-bit MATLAB R2019b [15].

###  IV-A Phasor Implementation

For the single-phase phasor model, it is assumed that the source voltage is
240 [V] and the load impedance $R+jX$ is such that it draws a current of
$10-j5$ [A]. For the three-phase phasor models, both line-ground and line-line
fault configurations are considered. These assume that the voltage source is
480 [V] rms line-line and the load impedance $R+jX$ is such that it draws
$30-i15$ [A] per phase. The fault resistance $R_{f}$ is selected such that
$R_{f}=R/10$.

Measured data is obtained by assuming a balanced input voltage and calculating
the current by multiplying the input voltage phasor vector with the admittance
matrix of the load-fault network. This is also the case for the single-phase
dynamic load, however, in that case the measured phasor voltage is converted
to instantaneous voltage to obtain the input for DSE.

###  IV-B Dynamic Implementation

The first case-study system is solved ad-hoc, assuming an ideal source with
the parameters listed in Table I.

Variable | Symbol | Value | Units
---|---|---|---
Total load real power | $P$ | 10 | kW
Total load reactive power | $Q$ | 5 | kVAR
Line-line RMS source voltage | $V_{ll}$ | 480 | V
Simulation time | T | 10 | ms
Sample rate | $T_{s}$ | 100 |  $\mu$s
TABLE I: Parameters for Single-Phase Dynamic Load

The latter two case-study systems are modeled in the MATLAB/Simulink® SimScape
multi-physics simulation environment, using the Specialized Power Systems
library with the parameters listed in Table II and Table III.

Variable | Symbol | Value | Units
---|---|---|---
Total load real power | $P$ | 10 | kW
Total load reactive power | $Q$ | 5 | kVAR
Line-line RMS source voltage | $V_{ll}$ | 240 | V
Source resistance | $R_{s}$ | 19.2 | $\Omega$
Source inductance | $L_{s}$ | 25.465 | mH
Fault resistance | $R_{f}$ | 1 | m$\Omega$
Ground resistance | $R_{g}$ | 10 | m$\Omega$
Cable positive-sequence resistance | $R_{c}$ | 183.7 | m$\Omega$
Cable positive-sequence reactance | $L_{c}$ | 26.6 | m$\Omega$
Simulation time | T | 200 | ms
Fault start time | $T_{f}$ | 50 | ms
TABLE II: Common Parameters for Three-Phase Dynamic Models Variable | Grounded-Wye | Delta
---|---|---
Load resistance R ($\Omega$) | 18.432 | 55.296
Load inductance L (mH) | 24.457 | 73.3
TABLE III: Varying Parameters for Three-Phase Dynamic Models

In the systems, the load is connected to a 480 [V] rms line-line source
through 1000 [ft] of 1/0 AWG quadruplex overhead service drop cable. Three
different cases are considered: 1) no-fault, 2) line-ground fault, and 3)
line-line fault.

###  IV-C Grounded-Wye Load

For the grounded-wye case, the system used is depicted in Fig. 3. The load
consists of three balanced series RL branches wired in a grounded-wye
configuration. This system has the parameters listed in Tables II and III.
Note that the total fault resistance for the line-ground fault is
$R_{f}+R_{g}=110$ [m$\Omega$], while the total fault resistance for the line-
line fault is $2R_{f}=200$ [m$\Omega$].

Figure 3: MATLAB Simulink model for a grounded-wye load with faults.

###  IV-D Delta Load

For the delta-connected case, the system used is depicted in Fig. 4. The load
consists of three balanced series RL branches wired in a delta configuration.
This system has the parameters listed in Tables II and III. Note that the
total fault resistance for the line-ground fault is $R_{f}+R_{g}=110$
[m$\Omega$], while the total fault resistance for the line-line fault is
$2R_{f}=200$ [m$\Omega$].

(a) Main model (b) Delta load

Figure 4: MATLAB Simulink model for a delta load with faults.

##  V Results

Results for the phasor models are presented in Table IV, while results for the
dynamic models are presented in Table V. The phasor models provide an
excellent estimate of the system parameters. DSE has difficulty estimating the
fault resistance for the dynamic model of grounded-wye network with a line-
line fault; a potential solution is to reduce the model order by neglecting
the load impedance on the faulted phases. Some moderate error is observed for
the case of the delta-connected load with a line-ground fault; again, it may
be possible to improve performance by neglecting load impedance on the faulted
phases.

Case | $R$ | $\hat{R}$ | $X$ | $\hat{X}$ | $R_{f}$ | $\hat{R}_{f}$
---|---|---|---|---|---|---
Single-Phase RL Load | 19.200 | 19.200 | 9.600 | 9.600 | – | –
Grounded-Wye Line-Ground Fault | 7.387 | 7.387 | 3.693 | 3.693 | 0.923 | 0.923
Grounded-Wye Line-Line Fault | 7.387 | 5.184 | 3.693 | 4.787 | 0.923 | 0.935
Delta Line-Line Fault | 22.160 | 22.160 | 11.080 | 11.080 | 2.770 | 2.770
Delta Line-Ground Fault | 22.160 | 22.160 | 11.080 | 11.080 | 2.770 | 2.770
TABLE IV: Results for Phasor State Estimation Case |  $R$ ($\Omega$) |  $\hat{R}$ ($\Omega$) |  $L$ (mH) |  $\hat{L}$ (mH) |  $R_{f}$ ($m\Omega$) |  $\hat{R}_{f}$ ($m\Omega$)
---|---|---|---|---|---|---
Single-Phase RL Load | 19.200 | 19.265 | 25.465 | 25.988 | – | –
Grounded-Wye No Fault | 18.432 | 18.404 | 24.446 | 24.485 | – | –
Grounded-Wye Line-Ground Fault | 18.432 | 18.432 | 24.446 | 24.446 | 11.000 | 10.997
Grounded-Wye Line-Line Fault | 18.432 | 18.432 | 24.446 | 24.446 | 11.000 | 3.165
Delta No Fault | 55.296 | 55.412 | 73.339 | 73.495 | – | –
Delta Line-Line Fault | 55.296 | 55.895 | 73.339 | 73.666 | 2.000 | 2.001
Delta Line-Ground Fault | 55.296 | 55.479 | 73.339 | 73.495 | 11.000 | 11.405
TABLE V: Results for Dynamic State Estimation

##  VI Conclusions

The results of performed experiments prove that DSE is capable of correctly
identifying model parameters of three different load configurations for both
normal and faulted operations. These load configurations model a lumped
equivalent of a radial electrical grid supplying multiple loads.

Several models showed sensitivity to inital conditions, particularly the
delta-connected load, so it is important that consideration be given to
providing the presented methods with good initial conditions. One issue is
making sure that there is a sufficient number of measurements to estimate
model states; for example, it is not possible to infer impedances for an
unbalanced delta-connected load given a single time snapshot. To reduce the
number of states, the models presented here assume that loads are balanced;
this assumption can be an issue for systems with a high degree of load
imbalance.

Existing work has demonstrated that DSE can operate with nonlinear elements
[16]. One option for future work is to expand the here presented methods with
other load models. These could include nonlinear voltage-dependent models
where power is a polynomial function of voltage (ZIP loads) or where power is
a polynomial function of both voltage and frequency (e.g. the WSCC load model
[17]). Alternately, these could include dynamic load models (e.g. an induction
motor model (MOTORW) or a composite load model (CMPLDW) [18]). Last, there is
the possibility of protecting line sections that include loads with
coordinated breakers at both ends; this could correspond to a distributed
parameter line or a Pi/Tee lumped equivalent model [19].

## References

  * [1] A. P. S. Meliopoulos, G. J. Cokkinides, P. Myrda, Y. Liu, R. Fan, L. Sun, R. Huang, and Z. Tan.  Dynamic State Estimation-Based Protection: Status and Promise.  IEEE Transactions on Power Delivery, 32(1):320–330, Feb. 2017\.
  * [2] Y. Liu, A. P. S. Meliopoulos, R. Fan, L. Sun, and Z. Tan.  Dynamic State Estimation Based Protection on Series Compensated Transmission Lines.  IEEE Transactions on Power Delivery, 32(5):2199–2209, 2017.
  * [3] Y. Liu, A. P. S. Meliopoulos, L. Sun, and R. Fan.  Dynamic State Estimation Based Protection of Mutually Coupled Transmission Lines.  CSEE Journal of Power and Energy Systems, 2(4):6–14, Dec. 2016\.
  * [4] Y. Liu, A. P. S. Meliopoulos, R. Fan, and L. Sun.  Dynamic State Estimation Based Protection of Microgrid Circuits.  In 2015 IEEE Power Energy Society General Meeting, pages 1–5, Jul. 2015.
  * [5] O. Vasios, S. Kampezidou, and A. P. S. Meliopoulos.  A Dynamic State Estimation Based Centralized Scheme for Microgrid Protection.  In Proceedings of the 2018 North American Power Symposium, pages 1–6, Sep. 2018.
  * [6] S. Choi and A. P. S. Meliopoulos.  Effective Real-Time Operation and Protection Scheme of Microgrids Using Distributed Dynamic State Estimation.  IEEE Transactions on Power Delivery, 32(1):504–514, Feb. 2017\.
  * [7] R. M. Tumilty, M. Brucoli, G. M. Burt, and T. C. Green.  Approaches to Network Protection for Inverter Dominated Electrical Distribution Systems.  In Proceedings of the 3rd IET International Conference on Power Electronics, Machines and Drives, 2006, pages 622–626, Apr. 2006.
  * [8] M. Dewadasa, A. Ghosh, and G. Ledwich.  Line Protection in Inverter Supplied Networks.  In Proceedings of the 2008 Australasian Universities Power Engineering Conference, pages 1–6, Dec. 2008.
  * [9] A. K. Barnes and A. Mate.  Implementing Admittance Relaying for Microgrid Protection.  In Proceedings of the 2021 IEEE/IAS 57th Industrial and Commercial Power Systems Technical Conference, pages 1–9, Apr. 2021.
  * [10] J. M. Dewadasa, A. Ghosh, and G. Ledwich.  Distance Protection Solution for a Converter Controlled Microgrid.  In Proceedings of the 15th National Power Systems Conference, 2008\.
  * [11] D. P. Kothari and I. J. Nagrath.  Modern Power System Analysis.  Tata McGraw-Hill Education, 1989.
  * [12] A. Monticelli.  Electric Power System State Estimation.  Proceedings of the IEEE, 88(2):262–282, Feb. 2000.
  * [13] S. Chapra and R. Canale.  Numerical Methods for Engineers.  McGraw-Hill Education, 2009.
  * [14] J. Bezanson, S. Karpinski, V. B. Shah, and A. Edelman.  Julia: A Fast Dynamic Language for Technical Computing.  arXiv:1209.5145, 2012.
  * [15] MathWorks.  MATLAB R2019b Documentation.  MATLAB – https://www.mathworks.com, 2020.
  * [16] H. F. Albinali.  State Estimation-Based Centralized Substation Protection Scheme.  PhD thesis, Georgia Institute of Technology, Aug. 2017.
  * [17] L. Pereira, D. Kosterev, P. Mackin, D. Davies, J. Undrill, and Z. Wenchun.  An Interim Dynamic Induction Motor Model for Stability Studies in the WSCC.  IEEE Transactions on Power Systems, 17(4):1108–1115, Nov. 2002\.
  * [18] North American Electric Reliability Corporation.  Technical Reference Document: Dynamic Load Modeling.  Technical report, NERC, 2016.
  * [19] D. T. Rizy and R. H. Staunton.  Evaluation of Distribution Analysis Software for DER Applications.  Technical report, Oak Ridge National Laboratory, 2002.