A Passivity Interpretation of Energy-Based Forced Oscillation Source Location Methods

TL;DR

This paper introduces a passivity-based framework to analyze and improve energy-based methods for locating forced oscillation sources in power systems, highlighting limitations and proposing a simulation-free performance prediction algorithm.

Contribution

It provides a systematic passivity interpretation of the DEF method, identifies its deficiencies, and develops a new algorithm for predicting its performance without simulations.

Findings

The DEF method's deficiencies are mathematically characterized.

No universal passivity transformation exists for all power system components.

The proposed algorithm accurately predicts DEF performance on standard test systems.

Abstract

This paper develops a systematic framework for analyzing how low frequency forced oscillations propagate in electric power systems. Using this framework, the paper shows how to mathematically justify the so-called Dissipating Energy Flow (DEF) forced oscillation source location technique. The DEF's specific deficiencies are pinpointed, and its underlying energy function is analyzed via incremental passivity theory. This analysis is then used to prove that there exists no passivity transformation (i.e. quadratic energy function) which can simultaneously render all components of a lossy classical power system passive. The paper goes on to develop a simulation-free algorithm for predicting the performance of the DEF method in a generalized power system, and it analyzes the passivity of three non-classical load and generation components. The proposed propagation framework and performance algorithm are both tested and illustrated on the IEEE 39-bus New England system and the WECC 179-bus system.

Figures

Tables

Table 1. TABLE I: Eigenvalues of System with Droop-Controlled Inverters

	$\alpha =0$	$\alpha ={10}^0$	$\alpha ={10}^1$	$\alpha ={10}^2$	$\alpha ={10}^3$	$\alpha ={10}^4$
${\lambda }_1$	0.007	0.007	0.008	0.01	0.05	0.4
${\lambda }_2$	2.11	2.11	2.12	2.17	2.80	28.04

Equations

$$\dot{\mathbf{x}}$$

$$\mathbf{y}$$

$$\Delta \dot{\mathbf{x}}$$

$$\Delta \mathbf{y}$$

$$\tilde{\mathbf{y}}=\underset{\mathcal{Y}(j\Omega )}{\underbrace{\left(C(j\Omega 1-A{)}^{-1}B+D\right)}}{\tilde{\mathbf{u}}}_v,$$

$$\displaystyle E=\left[\!\begin{array}[]{ccc}1&-1&\!\!\!\cdots\\ 0&1\vspace{-0.15cm}\\ \vdots&&\!\!\!\ddots\end{array}\!\right]\,\Leftrightarrow\,\,\,E_{a}=\left[\!\begin{array}[]{ccc}\mathbb{1}_{2}&-\mathbb{1}_{2}&\!\!\!\cdots\\ {\bf 0}&\mathbb{1}_{2}\vspace{-0.15cm}\\ \vdots&&\!\!\!\ddots\end{array}\!\right].$$

$${\bf V}_{b}=\left[\!\begin{array}[]{c}\tilde{V}_{r,1}\\ \tilde{V}_{i,1}\\ \vdots\\ \tilde{V}_{r,n}\\ \tilde{V}_{i,n}\end{array}\!\right],$$

$${\bf I}_{l}=\left[\!\begin{array}[]{c}\tilde{I}_{r,1}\\ \tilde{I}_{i,1}\\ \vdots\\ \tilde{I}_{r,m}\\ \tilde{I}_{i,m}\end{array}\!\right].$$

$${\mathcal{Y}}_L$$

$${\mathcal{Y}}_L{E}_a{\mathbf{V}}_b$$

$$\displaystyle\mathcal{Y}_{S}=\left[\begin{array}[]{ccc}\mathcal{Y}_{s,1}\!\!\!&&\!\!\!{\bf 0\vspace{-0.15cm}}\\ &\ddots\vspace{-0.15cm}\\ {\bf 0}\!\!\!&&\!\!\!\mathcal{Y}_{s,n}\end{array}\right].$$

$${\mathbf{I}}_s={\mathcal{Y}}_S{\mathbf{V}}_b.$$

$$\mathbf{J}=[0,\dots 0,\mathcal{I},0,\dots 0{]}^{\top }.$$

$$\mathbf{J}+{\mathbf{I}}_s+E_a^{†}{\mathbf{I}}_l=0.$$

$$-\mathbf{J}$$

$$=E_a^{†}{\mathbf{I}}_l+{\mathcal{Y}}_S{\mathbf{V}}_b$$

$$=\left(E_a^{†}{\mathcal{Y}}_L{E}_a+{\mathcal{Y}}_S\right){\mathbf{V}}_b.$$

$${\mathcal{Y}}_B=E_a^{†}{\mathcal{Y}}_L{E}_a+{\mathcal{Y}}_S.$$

$$-\mathbf{J}$$

$$\displaystyle\left[\begin{array}[]{c}-{\bm{\mathcal{I}}}\\ \hline\cr{\bf 0}\end{array}\right]$$

$${\mathcal{I}}^{\prime }=\mathcal{I}+{\mathcal{Y}}_{s,1}{\mathbf{V}}_s.$$

$$-{\mathbf{J}}^{\prime }$$

$$\displaystyle\left[\begin{array}[]{c}-{\bm{\mathcal{I}}}^{\prime}\\ \hline\cr{\bf 0}\end{array}\right]$$

$$-{\mathcal{I}}^{\prime }$$

$${\mathbf{V}}_s=-{\mathcal{Z}}_N{\mathcal{I}}^{\prime },$$

$$0$$

$$={\mathbf{V}}_b^{†}\left(E_a^{†}{\mathbf{I}}_l+{\mathbf{I}}_s+\mathbf{J}\right)$$

$$={\mathbf{V}}_b^{†}\left({\mathcal{Y}}_B-{\mathcal{Y}}_B\right){\mathbf{V}}_b,$$

$$\mathcal{Q}$$

$$-{\mathbf{V}}_b^{†}\mathcal{Q}{\mathbf{J}}^{\prime }$$

$$-{\mathbf{V}}_s^{†}{\mathcal{Q}}_b{\mathcal{I}}^{\prime }$$

$$\displaystyle{\mathcal{Q}}_{b}=\left[\!\begin{array}[]{cc}0\!&-\frac{1}{j\Omega}\\ \frac{1}{j\Omega}\!&0\end{array}\!\right],$$

$$\underset{\text{Source Energy}}{\underbrace{\mathrm{Re}\{{\mathbf{V}}_s^{†}{\mathcal{Q}}_b{\mathcal{I}}^{\prime }\}}}+\sum _{i=2}^n\underset{\text{Generator Damping Contributions}}{\underbrace{\mathrm{Re}\{{\mathbf{V}}_{ns,i}^{†}\left({\mathcal{Q}}_b{\mathcal{Y}}_{s,i}\right){\mathbf{V}}_{ns,i}\}}}=0,$$

$$\underset{\text{Source Energy}}{\underbrace{\mathrm{Re}\{{\mathbf{U}}_s^{†}{\mathcal{Q}}_b{\mathcal{I}}^{\prime }\}}}+\sum _{i=2}^n\underset{\text{Generator Damping Contributions}}{\underbrace{\text{Re{}{\mathbf{U}}_{ns,i}^{†}\left({\mathcal{Q}}_b{\mathcal{Y}}_{s,i}{\mathcal{P}}_b\right){\mathbf{U}}_{ns,i}\}}}=0.$$

$${\overline{T}}_1$$

Full text

I Introduction

Forced oscillations (FOs) are still a problematic reality in modern electric power systems. Caused by extraneous periodic perturbations, FOs can compromise system security, degrade system performance, and resonantly excite poorly damped interarea modes [1, 2, 3]. Despite widescale deployment of Phasor Measurement Units (PMUs) across the US high voltage transmission network, locating the sources of these FOs remains a challenging task due to their sporadic nature, speed of propagation, and inability to be predicted by the system operators’ dynamical models. Since the most effective way for dealing with a FO is to locate the source and disconnect it from service [4], an effective source location technique is an indispensable smart grid application.

Of the many source location techniques currently available in the academic marketplace [5], the so-called Dissipating Energy Flow (DEF) method has enjoyed some of the most successful testing results, both in simulation environments [4] and in real-time applications [6]. The method was originally developed by Chen et al. [7], but its underlying mathematics leverage the Lyapunov functions from [8]. Despite its success, when its underlying modeling assumptions are violated, the method may perform poorly [9, 10, 11]. Additionally, the method has lacked a generalized framework which is capable of providing a system-wide justification for its usage.

The primary challenge to reliable DEF performance is the contribution of dissipating energy from non-FO sources, such as lossy transmission lines, lossy or negatively damped loads, and generator dynamics dominated by non-passive controllers. An actual example of this contribution may be found in Fig. 1 of [11]. If these contributions are large enough, the FO source can appear to be a dissipating energy sink and the DEF method can fail. Despite its inadequacies, the DEF’s excellent performance in real-time application at ISO New England strongly implies that further research should be performed in order to more systematically characterize the method. Shortcomings of the DEF method have been analyzed in [9], and [11] has recommended using passivity theory to interpret the method from a new mathematical perspective, but no theoretical methods have been devised for testing how the DEF method will perform in an arbitrary network. Such testing is essential in order to ensure that the DEF can perform adequately in new environments, such as in microgrids where “R/X” line ratios are high and voltage control is fast, or in networks that have particularly resistive load pockets. To make such predictions, a systematic framework is needed in order to thoroughly study the DEF.

Accordingly, this paper develops a framework which is suitable for analyzing the DEF method via passivity theory in the context of a full-scale system. Initially, this framework is used to analyze a lossy classical power system, but the methods are then generalized to include arbitrary system models; these generalized methods are thus capable of investigating all of the problematic effects presented by lines, loads and generator controllers. The primary contributions of this paper are as follows:

Using the Frequency Response Function (FRF) analysis proposed in [10], we leverage a variety of tools from AC circuit theory in order to develop a linearized framework for analyzing oscillation propagation at the system level. 2. 2.

We subsequently use the proposed framework, along with the passivity observations from [11], to theoretically justify the DEF method and show that there exists no other quadratic passivity transformation which will render all components of the entire network passive in a classical power system. 3. 3.

A simulation-free algorithm is developed which predicts the performance of the DEF method in a generalized power system. 4. 4.

The proposed passivity transformation is used to analytically investigate the dissipating energy flow properties of three common, yet non-classical, grid components.

The remainder of this paper is structured as follows. In Section II, we derive the linearized propagation framework and show how quadratic energy is conserved in the system. In Section III, we extend the energy function analysis and prove that there is no quadratic energy function which will render all components of a classical power system passive. We then leverage the proposed framework in Section IV in order to develop an algorithm which analytically predicts the performance of the DEF method for a generalized power system. We present illustrative and supportive test results in Section V and offer concluding remarks in Section VI.

II Perturbative Network Model Derivation

In this section, we introduce a linearized network model which is particularly useful for analyzing FO propagation in power systems. We refer to it as a “perturbative” model since it captures the network’s response to small111The framework presented in this paper is valid when FOs are sufficiently “small”, such that quadratic nonlinearities may be neglected. perturbations.

II-A Complex Admittance Matrices

We consider a power system component, shown in Fig. 1, whose dynamics are governed according to the DAE set

[TABLE]

where inputs ${\bf u}_{v}=\left[V_{r},\,V_{i}\right]$ and outputs ${\bf y}=\left[I_{r},\,I_{i}\right]$ are vectors of real and imaginary voltages and currents, respectively.

We linearize this model around a steady state operating point to linearly relate the voltage and current perturbations:

[TABLE]

Assuming (2) is BIBO stable, its Fourier transform admits the Frequency Response Function (FRF) of the component [10]:

[TABLE]

where $\Omega$ is the angular frequency of the input and output signals, and $\mathcal{Y}\equiv\mathcal{Y}(j\Omega)\!\in\!{\mathbb{C}}^{2\times 2}$ is referred to as the admittance matrix relating voltage $\tilde{{\bf u}}_{v}\!\!\in\!\!{\mathbb{C}}^{2\times 1}$ and current $\tilde{{\bf y}}\!\!\in\!\!{\mathbb{C}}^{2\times 1}$ perturbations. These perturbations can be given in rectangular (${\tilde{V}}_{r}$, ${\tilde{V}}_{i}$) or polar (${\tilde{\rm V}}$, ${\tilde{\theta}}$) coordinates depending on convenience. In this paper, we will primarily consider the effects of the FRF matrix at the relevant forcing frequency $\Omega_{d}$ of the FO. When referring to $\mathcal{Y}$, we will use the terms FRF and admittance interchangeably, depending on the context.

II-B Network Modeling

Consider a power system network whose graph ${\mathcal{G}}({\mathcal{V}},{\mathcal{E}})$ has edge set $\mathcal{E}$, $|\mathcal{E}|=m$, vertex set $\mathcal{V}$, $|\mathcal{V}|=n$, and directed nodal incidence matrix $E\in{\mathbb{R}}^{m\times n}$ [12]. In power system modeling, it is standard practice to use the incidence matrix to build the nodal admittance matrix $Y_{\rm bus}$ from the primitive admittance matrix $Y_{p}$ via $Y_{\rm bus}=E^{\dagger}Y_{p}E$ [13]. This procedure is possible because the admittances in $Y_{p}$ are complex values $\in{\mathbb{C}}^{1}$. In the following analysis, all admittances in the form of (3) must be expressed in $2\times 2$ matrix form (i.e. $\mathcal{Y}\in{\mathbb{C}}^{2\times 2}$). Accordingly, we need an incidence matrix which can relate $2\times 2$, rather than $1\times 1$, admittance values. We thus build a so-called “augmented” incidence matrix $E_{a}\in{\mathbb{R}}^{2m\times 2n}$. This matrix is constructed by taking $E$ and replacing all values of 1 with the $2\times 2$ identity matrix $\mathbb{1}_{2}$ and all values of 0 with a $2\times 2$ zero matrix ${\bf 0}$:

[TABLE]

Considering voltage and current perturbation phasors such as $\tilde{{\bf u}}_{v}$ and $\tilde{{\bf y}}$ from (3), we define the vector ${\bf V}_{b}\in{\mathbb{C}}^{2n\times 1}$ as the vector of rectangular bus voltage perturbation phasors, and we define the vector ${\bf I}_{l}\in{\mathbb{C}}^{2m\times 1}$ as the vector of rectangular line current perturbation phasors, where the convention of the positive line current flows agrees with the direction of the augmented incidence matrix:

{\bf V}_{b}=\left[\!\begin{array}[]{c}\tilde{V}_{r,1}\\ \tilde{V}_{i,1}\\ \vdots\\ \tilde{V}_{r,n}\\ \tilde{V}_{i,n}\end{array}\!\right],

(11)

{\bf I}_{l}=\left[\!\begin{array}[]{c}\tilde{I}_{r,1}\\ \tilde{I}_{i,1}\\ \vdots\\ \tilde{I}_{r,m}\\ \tilde{I}_{i,m}\end{array}\!\right].

(12)

Admittance matrices $\mathcal{Y}_{l,1},...,\mathcal{Y}_{l,m}$, $\mathcal{Y}_{l,i}\in\mathbb{R}^{2\times 2}$, associated with the $m$ transmission lines, are placed diagonally in the matrix ${\mathcal{Y}}_{L}\in{\mathbb{R}}^{2m\times 2m}$ such that

[TABLE]

Transformers with off-nominal tap ratios, such as tap changing transformers, are discussed in Appendix A. Line current and bus voltage perturbations obey Ohm’s law:

[TABLE]

Admittance matrices $\mathcal{Y}_{s,1},...,\mathcal{Y}_{s,n}$, $\mathcal{Y}_{s,i}\in\mathbb{C}^{2\times 2}$, associated with shunt elements at each of the $n$ buses, are placed diagonally in the matrix ${\mathcal{Y}}_{S}\in{\mathbb{C}}^{2n\times 2n}$, such that

[TABLE]

These shunt admittances are not simply capacitors or inductors; they can represent generators, loads, or any other terminal element in the system and can be constructed via (3). If multiple elements are connected in parallel, such as a generator and its station load, their admittances can be modeled independently and summed to compute the aggregate shunt admittance. The shunt matrix ${\mathcal{Y}}_{S}$ may be used to compute the shunt current injections222Shunt currents flowing out of the circuit to ground are defined as positive. ${\bf I}_{s}\!\in\!{\mathbb{C}}^{2n\times 1}$ via

[TABLE]

As outlined in [10], when analyzing a network with this representation, FOs show up like current sources at their respective source buses. For a system experiencing a single FO, there will be a single current source ${\bm{\mathcal{I}}}\in{\mathbb{C}}^{2\times 1}$ driving the entire network. We define a sparse FO vector of current injections ${\bf J}\in\mathbb{C}^{2n\times 1}$ whose structure will take the form

[TABLE]

If bus $k$ is the source of the FO, then ${\bm{\mathcal{I}}}$ will be located in elements ${\bf J}_{2k-1}$ and ${\bf J}_{2k}$. Source injections in ${\bf J}$ obey the same current convention as ${\bf I}_{s}$. The network obeys KCL, i.e., all nodal currents sum to 0:

[TABLE]

We define ${\bf I}_{I}\equiv E_{a}^{\dagger}{\bf I}_{l}$ to be the aggregate current injection at each node: it represents the sum of the source current injection at each bus plus the shunt current flowing to ground. Via conservation of current at each bus, we have

[TABLE]

The block-Hermitian dynamic nodal admittance (or augmented dynamic Y-bus) matrix $\mathcal{Y}_{B}\in{\mathbb{C}}^{2n\times 2n}$ is this defined to be

[TABLE]

Assuming there is a single FO in the system, it is instructive to rewrite (25c) with partitioned matrices and vectors, where the system has been renumbered such that the source bus is bus 1, and the current injection has a value of ${\bm{\mathcal{I}}}\in{\mathbb{C}}^{2\times 1}$:

[TABLE]

where ${\bf V}_{s}\!\in\!{\mathbb{C}}^{2\times 1}$ represents voltage perturbations at the source bus, ${\bf V}_{ns}\!\in\!{\mathbb{C}}^{(2n-2)\times 1}$ represents voltage perturbations at all other buses, and ${\bf V}_{b}={\bf V}_{s}^{\frown}{\bf V}_{ns}$. While ${\bm{\mathcal{I}}}$ represents the true source current injection, we may also define ${\bm{\mathcal{I}}}^{\prime}$ as the sum of the source current at the source bus plus its shunt injection:

[TABLE]

Correspondingly, we say that ${\mathcal{Y}}_{B1}^{\prime}$ contains no shunt element, and ${\bm{\mathcal{I}}}^{\prime}$ is the current directly measured at the source bus flowing into the network; we note that it is equal to the first two elements of ${\bf I}_{I}$. We now restate (27a)-(27h) with this update:

[TABLE]

A simple Kron reduction can be performed in order to determine the effective admittance “seen” by the current source:

[TABLE]

where $\mathcal{Y}_{N}\!\in\!{\mathbb{C}}^{2\times 2}$ is a $2\times 2$ complex aggregate network admittance matrix and ${\bf V}_{s}$ is the resulting voltage caused by the current injection ${\bm{\mathcal{I}}}^{\prime}$ interacting with the aggregate network dynamics codified in $\mathcal{Y}_{N}$. In this model, since voltage perturbations are considered a response to rouge current injections, it is helpful to rewrite (30) with voltage as a function of current:

[TABLE]

where $\mathcal{Z}_{N}\!=\!\mathcal{Y}_{N}^{-1}$ is the aggregate network impedance. In other words, the current injection ${\bm{\mathcal{I}}}^{\prime}$ gives rise to the network voltages, and the vector ${\bf V}_{b}$ in (29a) is not arbitrary: the Kron reduction of (30) is only meaningful when ${\bf V}_{b}$ acts as a solution to the linear system $-{\bf J}^{\prime}=\mathcal{Y}_{B}^{\prime}{\bf V}_{b}$, i.e. ${\bf V}_{b}=-(\mathcal{Y}_{B}^{\prime})^{-1}{\bf J}^{\prime}$.

Definition 1**.**

We refer to the admittance matrix ${\mathcal{Y}}_{N}$ of (30) as the system’s dynamic Ward equivalent (DWE) admittance.

II-C Quadratic Energy Considerations

As with any network which obeys Kirchhoff’s laws, Tellegen’s theorem is also obeyed: the sum of the products of branch (including shunt branches) potential differences and branch flows is equal to 0. Accordingly,

[TABLE]

where (32a) is the statement of Tellegen’s theorem, (32b) is the conservation of current, and (32c) is the resulting proof. As a consequence of this theorem, there exist a family of quadratic functionals for which conservation laws can be formulated. An obvious one is the “real power”, $\rm{Re}\{{\bf V}{\bf I}^{\dagger}\}$, that is consumed only on the elements with positive resistance. In the context of FOs, an alternative interpretation of the conservation of power can be acquired by manipulating (29h) in order to define another (arbitrary) type of quadratic power333The term “quadratic power” is used since we are multiplying voltages and currents, but the quadratic quantity doesn’t necessarily have the interpretation of physical power. It can also be interpreted as an “energy function”. Quadratic energy and quadratic power are therefore used interchangeably.. The key observation is that this new quadratic power will be conserved throughout the network. To show why, we consider matrix ${\mathcal{Q}}\in{\mathbb{C}}^{2n\times 2n}$ with block diagonal sub-matrices ${\mathcal{Q}}_{b}\in{\mathbb{C}}^{2\times 2}$:

[TABLE]

We now left multiply (29a) by ${\bf V}_{b}^{\dagger}{{\mathcal{Q}}}$, which represents the application of a quadratic energy function:

[TABLE]

Different choices for matrix ${{\mathcal{Q}}}_{b}$ correspond to different energy function applications, but in each case, the quadratic quantity is conserved. For example, if ${\mathcal{Q}}_{b}$ is chosen such that

[TABLE]

then the associated energy function corresponds to the DEF method [11]. Under DEF assumptions, lines and loads are rendered lossless, i.e $(\mathcal{Q}_{b}\mathcal{Y})+(\mathcal{Q}_{b}\mathcal{Y})^{\dagger}={\bf 0}$. Thus, in taking the real part, (37b) simplifies to

[TABLE]

where ${\bf V}_{s}$ and ${\bf V}_{ns,i}$ are the source and $i^{\rm th}$ non-source bus voltage perturbation vectors, respectively. The formulation of (41) further clarifies the DEF’s functionality: since all the damping energy consumed by generators is positive [11], the source energy is necessarily negative and can be traced back to the single, negative source. The DEF technique, therefore, is based on tracking a particular type of quadratic power in the network. When constructing the system’s quadratic energy function, we are not restricted to choosing just a ${\mathcal{Q}}$ matrix. We may also introduce matrix ${\mathcal{P}}$ whose structure matches ${\mathcal{Q}}$. For example, we may set ${\bf U}_{b}={\mathcal{P}}^{-1}{\bf V}_{b}$, left multiply (29a) by ${\bf U}_{b}^{\dagger}{{\mathcal{Q}}}$, and insert a ${\mathcal{P}}{\mathcal{P}}^{-1}$ term. Updating (41) yields:

[TABLE]

III Energy Function Analysis

The DEF method can be interpreted as the application of a specific quadratic energy function to all elements of a network. Reference [11], which reviews relevant passivity concepts, explains how this quadratic energy function can be interpreted as a matrix transformation (called a “passivity transformation”) which attempts to render system elements incrementally passive. Generally, a passive component is one which can only dissipate and store, but not produce, physical power. To avoid confusion, we point out that this paper discusses the so-called incremental passivity, which refers to the passive nature of a system’s incremental change from its equilibrium, as defined by [14, eq. (4.159)]. This is also known as “shifted passivity”. We stress that physical passivity of an element does not imply incremental passivity when some transformation, such as (36), has been applied. In the remainder of the paper, the term passivity always refers to incremental passivity.

Reference [11] further shows that the DEF energy function is inadequate for lossy network elements. While it may be tempting to develop a new passivity transformation which is suitable for lossy networks, in this section, we use passivity theory to prove that no constant quadratic energy function exists for the classical model of a multimachine power system. For coherency, we first offer the following clarification of how we use the terms “lossy” and “lossless” in this paper.

Definition 2**.**

In keeping with the conventional nomenclature for both power systems and passivity literature, we offer the following definitions for lossy and lossless:

• •

Lossy*: A lossy element is an element which contains resistance; a lossy system is a system which contains resistive transmission lines.*

• •

Lossless444While “lossless” certainly can also refer to a transmission network with purely reactive lines, we do not invoke that definition anywhere in this paper.: The term lossless describes a FRF whose Hermitian part is $\bf 0$. If ${\mathcal{Y}}+{\mathcal{Y}}^{\dagger}={\bf 0}$, ${\mathcal{Y}}$ is said to be lossless; if ${\bf M}{\mathcal{Y}}+({\bf M}{\mathcal{Y}})^{\dagger}={\bf 0}$, ${\bf M}{\mathcal{Y}}$ is said to be lossless.

III-A Basis Matrices

To aide in the energy function analysis and inference techniques, we define a useful set of basis matrices.

Definition 3**.**

We define orthogonal ($A^{-1}\!=\!A^{\dagger}$) basis matrices

[TABLE]

and the set $\mathbb{T}=\{\overline{T}_{1},\overline{T}_{2},\overline{T}_{3},\overline{T}_{4}\}$. Set $\mathbb{T}$ spans region $\mathbb{R}^{2\times 2}$:

[TABLE]

Lemma 1**.**

There exists no non-singular matrix $\Gamma\in\mathbb{C}^{2\times 2}$ for which, simultaneously,

[TABLE]

Proof.

We write $\Gamma$ as the sum of its diagonal ($\Gamma_{d}$) and off-diagonal ($\Gamma_{o}$) component matrices: $\Gamma=\Gamma_{d}+\Gamma_{o}$. Since ${\rm Re}\{{{\bf v}}^{\dagger}(\overline{T}_{i}\Gamma){{\bf v}}\}=0,\;\forall{{\bf v}}\in\mathbb{C}^{2\times 2}$ is equivalent to stating that $\overline{T}_{i}\Gamma+(\overline{T}_{i}\Gamma)^{\dagger}={\bf 0}$, the constraints on $\Gamma$ caused by ${\overline{T}}_{2}$, ${\overline{T}}_{3}$, and ${\overline{T}}_{4}$ from (45) may be stated as

[TABLE]

where $K\equiv{\overline{T}}_{3}$ is defined to be the reversal matrix in [15]. Accordingly, $\Gamma_{d}$ must be simultaneously skew-Hermitian, skew-perhermitian and perhermitian, respectively [15]. Necessarily, $\Gamma_{d}={\bf 0}$. The matrix $\Gamma_{o}$ must be simultaneously Hermitian and skew-perhermitian. Necessarily, $j\beta{\overline{T}}_{4},\,\beta\!\in\mathbb{R}^{1}$, is the only matrix which fits this description. We define

[TABLE]

as the only matrix which uniformly satisfies (45). We apply $\Gamma=\Gamma^{\star}$ to (44) and consider the eigenvalues of the matrix ${\overline{T}}_{1}\Gamma^{\star}+\left({\overline{T}}_{1}\Gamma^{\star}\right)^{\dagger}=2\beta j{\overline{T}}_{4}$:

[TABLE]

which violates (44). Since ${\overline{T}}_{1}\Gamma^{\star}\!+\!({\overline{T}}_{1}\Gamma^{\star})^{\dagger}$ is an indefinite matrix but $\Gamma^{\star}$ is the only matrix which satisfies (45), the theorem has been proved. ∎

Corollary 1**.**

$\Gamma^{\star}$* is a solution to $P\Gamma^{\star}+(P\Gamma^{\star})^{\dagger}={\bf 0}$ for any matrix $P$ which may be written as $P=\sum_{i=2}^{4}\alpha_{i}{\overline{T}}_{i},\,\alpha_{i}\!\in\!{\mathbb{R}}^{1}$.*

Corollary 2**.**

The results of Lemma 1 stand if ${\overline{T}}_{i}$ is right multiplied by transformation matrix ${\bf M}$ instead of left multiplied by transformation matrix $\Gamma$. There exists no non-singular matrix ${\bf M}\in\mathbb{C}^{2\times 2}$ for which the following simultaneously hold:

[TABLE]

Corollary 3**.**

By employing both non-singular matrices $\Gamma\in\mathbb{C}^{2\times 2}$ and ${\bf M}\in\mathbb{C}^{2\times 2}$, the solution to

[TABLE]

must take the form $\Gamma=(j\beta\overline{T}_{4}){\bf M}^{\dagger}$ for any ${\bf M}\in{\mathbb{C}}^{2\times 2}$. This may be seen via the following manipulation:

[TABLE]

Since (55b) may only be solved by $\Gamma{\bf M}^{\dagger^{-1}}=j\beta\overline{T}_{4}$, per Lemma 1, we have that $\Gamma=(j\beta\overline{T}_{4}){\bf M}^{\dagger}$ must be satisfied.

III-B Quadratic Energy Functions in a Classical Power System

We now assume the classical model of a lossy multimachine power system model [16] and allow for constant power loads to be present. The forms of the FRFs associated with constant power loads (${\mathcal{Y}}_{p}$), constant impedance lines/shunts (${\mathcal{Y}}_{z}$), and classical generators (${\mathcal{Y}}_{g}$) are given in [11]. The set of plausible FRFs associated with these three elements may be constructed according to the following basis matrix combinations:

[TABLE]

The FRF of a classical generator is derived and fully explained in [10], and, for convenience, is explicitly re-stated here:

[TABLE]

where $\delta$ is the generator’s absolute rotor angle and $\gamma\equiv\gamma(\Omega)\!\in\!{\mathbb{C}}^{1}$ is a complex frequency dependent parameter. Matrix (61) has many useful properties, one of which will be addressed in subsection III-C.

Given that any admittance matrix $\mathcal{Y}$ may be written as the complex sum of the four weighted basis matrices, we introduce the following useful definition:

Definition 4**.**

We assume some admittance $\mathcal{Y}$ may be written as ${\mathcal{Y}}=\sum_{i=1}^{4}\left(a_{i}+jb_{i}\right)\overline{T}_{i}$. Using matrices $\bf M$ and $\Gamma$ from Corollary 3, where $\Gamma=(j\beta\overline{T}_{4}){\bf M}^{\dagger}$, we define $P^{\star}={\rm Re}\{{{\bf u}}^{\dagger}\left({\bf M}\mathcal{Y}\Gamma\right){{\bf u}}\}$ as the dissipating power for input vector ${\bf u}$. We further define two other types of quadratic power: resistive power $P_{r}^{\star}$ and damping power $P_{d}^{\star}$, where $P^{\star}=P_{r}^{\star}+P_{d}^{\star}$, and

• •

$P_{r}^{\star}={\rm Re}\left\{{{\bf u}}^{\dagger}\left({\bf M}\left(a_{1}\overline{T}_{1}\right)\Gamma\right){{\bf u}}\right\}$**

• •

$P_{d}^{\star}={\rm Re}\left\{{{\bf u}}^{\dagger}\left({\bf M}\left(jb_{2}\overline{T}_{2}+jb_{3}\overline{T}_{3}+jb_{4}\overline{T}_{4}\right)\Gamma\right){{\bf u}}\right\}$.

We now prove that a perturbative system model containing elements (56a)-(56c) cannot be rendered passive under any quadratic passivity transformation. In other words, there is no quadratic quantity that is dissipated by all elements present in common networks.

Theorem 1**.**

There exist no non-singular matrices ${\bf M}\in\mathbb{C}^{2\times 2}$ and $\Gamma\in\mathbb{C}^{2\times 2}$ for which

[TABLE]

Proof.

The FRF of a strictly reactive element, such as matrix $T_{X}$ in (61), is $\propto\overline{T}_{4}$ while the FRF of a strictly capacitive element is $\propto-\overline{T}_{4}$. The only way for ${\bf M}\overline{T}_{4}\Gamma+\left({\bf M}\overline{T}_{4}\Gamma\right)^{\dagger}\succeq 0$ and $-{\bf M}\overline{T}_{4}\Gamma-\left({\bf M}\overline{T}_{4}\Gamma\right)^{\dagger}\succeq 0$ to be simultaneously true is for ${\bf M}\overline{T}_{4}\Gamma+\left({\bf M}\overline{T}_{4}\Gamma\right)^{\dagger}\equiv{\bf 0}$. Since both reactive and capacitive elements appear in classical power systems, $\overline{T}_{4}$ must be lossless under the desired passivity transformation.

We now consider some classical generator whose damping characteristics are sufficiently small ($D\approx 0$), such that $\gamma$ is a real parameter. In this case, the matrix ${\bf M}\mathcal{Y}_{g}\Gamma+\left({\bf M}\mathcal{Y}_{g}\Gamma\right)^{\dagger}$ reduces to ${\bf M}\gamma T_{\delta}\Gamma+\left({\bf M}\gamma T_{\delta}\Gamma\right)^{\dagger}$ since ${T_{X}}$ must be a lossless element according to the previous conclusion about $\overline{T}_{4}$.

We define the squared electromechanical resonant frequency associated with the classical generator as $\Omega^{2}_{r}=\frac{{\rm V}{\rm E}^{\prime}}{MX^{\prime}_{d}}\cos(\varphi)$. For some $\epsilon$, $\gamma(\Omega_{r}-\epsilon)=-\gamma(\Omega_{r}+\epsilon)$. We must therefore ensure that ${\bf M}\gamma T_{\delta}\Gamma+\left({\bf M}\gamma T_{\delta}\Gamma\right)^{\dagger}\succeq 0$ and $-{\bf M}\gamma T_{\delta}\Gamma-\left({\bf M}\gamma T_{\delta}\Gamma\right)^{\dagger}\succeq 0$, respectively, when ensuring the generator’s passivity on either side of the resonant peak. The only way for these statements to be simultaneously true is for ${\bf M}T_{\delta}\Gamma+\left({\bf M}T_{\delta}\Gamma\right)^{\dagger}\equiv{\bf 0}$. To accomplish this, we consider the numerical structure of $T_{\delta}$ for two plausible rotor angle values: $\delta_{1}=0$ and $\delta_{2}=\frac{\pi}{4}$:

[TABLE]

Since ${\bf M}\overline{T}_{4}\Gamma+\left({\bf M}\overline{T}_{4}\Gamma\right)^{\dagger}\equiv{\bf 0}$, then we must also require ${\bf M}\overline{T}_{3}\Gamma+\left({\bf M}\overline{T}_{3}\Gamma\right)^{\dagger}\equiv{\bf 0}$ from (66) and ${\bf M}\overline{T}_{2}\Gamma+\left({\bf M}\overline{T}_{2}\Gamma\right)^{\dagger}\equiv{\bf 0}$ from (69) in order to ensure that ${\bf M}T_{\delta}\Gamma+\left({\bf M}T_{\delta}\Gamma\right)^{\dagger}\equiv{\bf 0}$. We are thus requiring that ${\rm Re}\left\{\tilde{{\bf v}}^{\dagger}\left({\bf M}\overline{T}_{i}\Gamma\right)\tilde{{\bf v}}\right\}=0,\;\forall\tilde{{\bf v}}\in\mathbb{C}^{2\times 1},\,\forall i\in\{2,3,4\}$. As stated in Corollary 3, the only way to achieve losslessness for basis matrices 2, 3 and 4 is for $\Gamma=(j\beta\overline{T}_{4}){\bf M}^{\dagger}$. By employing this transformation, Lemma 1 proves that the quadratic energy associated with any element containing $\overline{T}_{1}$ will be rendered indefinite in sign. Since (56b) contains $\overline{T}_{1}$ when resistance is present in the network, then there exists no nonsingular matrices ${\bf M}$ and $\Gamma$ for which ${\bf M}\mathcal{Y}\Gamma+\left({\bf M}\mathcal{Y}\Gamma\right)^{\dagger}\succeq 0,\quad\forall\mathcal{Y}\in\left\{\mathcal{Y}_{p},\;\mathcal{Y}_{z},\;\mathcal{Y}_{g}\right\}$. ∎

Corollary 4**.**

By choosing ${\bf M}=\overline{T}_{1}$, $\beta=\frac{-1}{\Omega}$, and $\Gamma=(j\beta\overline{T}_{4}){\bf M}^{\dagger}$, we arrive at the passivity transformation implicitly employed by the DEF, as given by[11, eq. (38a,b)].

Corollary 5**.**

The passivity transformation $\Gamma=(j\beta\overline{T}_{4}){\bf M}^{\dagger}$ renders constant power and reactive (capacitive or inductive) elements of the system lossless, classical generators passive (see[11, eq. (49)]), and resistive elements non-passive. Without resistive elements, the DEF method is a fully reliable source location technique in classical power systems.

Corollary 6**.**

Since the linearized admittance matrix ${\mathcal{Y}}$ associated with any ZIP load is purely real, i.e. ${\mathcal{Y}}={\rm Re}\{{\mathcal{Y}}\}$, the eigenvalues of its transformed Hermitian part will be equal and opposite in value: $\lambda({\bf M}\mathcal{Y}\Gamma+\left({\bf M}\mathcal{Y}\Gamma\right)^{\dagger})=\pm\alpha$, where $\alpha=0$ for constant power or purely reactive loads.

If a network has no resistive elements and is truly passive, no quadratic energy production can occur on regular network elements, so injections of energy have to be related to external sources, like FOs. This is why finding a passive energy-like form is important, and why the almost-passive form used by DEF has had so much success. To further show why the results of Theorem 1 are problematic for the DEF method, we consider the structure of (41): since the generator damping contributions are positive definite, the source energy is necessarily negative in a power system with no resistance. This is not true in a lossy power system. To show why not, we define block matrices $\underline{\bf M}\in{\mathbb{C}}^{2n\times 2n}$ and $\underline{\Gamma}\in{\mathbb{C}}^{2n\times 2n}$, where

[TABLE]

and whose block diagonal matrices are given by ${\bf M}$ and $\Gamma=(j\beta\overline{T}_{4}){\bf M}^{\dagger}$, respectively. We left multiply (29a) by $\underline{{\bf M}}$ and insert $\underline{\Gamma}\underline{\Gamma}^{-1}$ on the RHS:

[TABLE]

Defining the transformed voltage vectors ${\bf U}_{b}=\underline{\Gamma}^{-1}{\bf V}_{b}$ and ${\bf U}_{s}={\Gamma}^{-1}{\bf V}_{s}$, we left multiply (70) by ${\bf U}_{b}^{\dagger}$ and simplify. We may group the dissipating power injections into their respective contributing groups (assuming lossless loads):

[TABLE]

The FO source term can produce only negative damping energy, i.e. ${\rm Re}\{{\bf U}_{s}^{\dagger}{\bf M}{\bm{\mathcal{I}}}^{\prime}\}{<}0$, if the condition

[TABLE]

is met. If it is not met, indefinitely signed resistive energy can dominate damper winding energy absorption and the source term can, in fact, appear as a positively damped element. In this plausible situation, the DEF method will fail. We note that violation of (72) is a necessary but not sufficient condition for DEF failure. Next, we show that the sign of the injected resistive energy associated with a transmission network can be negated if all system voltages are complex conjugated.

Theorem 2**.**

Consider a lossy transmission network (just the network). For any transformed voltage vector ${\bf U}_{b}$ which yields quadratic energy ${\rm Re}\{{\bf U}_{b}^{\dagger}(\underline{{\bf M}}(E_{a}^{\dagger}\mathcal{Y}_{L}E_{a})\underline{\Gamma}){\bf U}_{b}\}\!=\!P_{r}^{\star}$, there exists conjugated vector ${\bf U}^{*}_{b}$ which yields an equal and opposite quadratic energy ${\rm Re}\{{\bf U}_{b}^{*\dagger}(\underline{{\bf M}}(E_{a}^{\dagger}\mathcal{Y}_{L}E_{a})\underline{\Gamma}){\bf U}_{b}^{*}\}=-P_{r}^{\star}$, where $\underline{{\bf M}}$ and $\underline{\Gamma}$, with submatrices $\bf M$ and $\Gamma$, yield from Corollary 3.

Proof.

We split the transmission line matrix into its conductive and susceptive parts: $\mathcal{Y}_{L}=\mathcal{Y}_{\bar{G}}+\mathcal{Y}_{\bar{B}}$, where $\mathcal{Y}_{{\bar{G}},i}=G_{i}\overline{T}_{1}$ and $\mathcal{Y}_{{\bar{B}},i}=B_{i}\overline{T}_{4}$. We also define block matrices ${\bf M}=j\overline{T}_{4}\Gamma$ and $\Gamma=\overline{T}_{1}$ from $\underline{\bf M}$ and $\underline{\Gamma}$. Therefore, $\underline{{\bf M}}\left(E_{a}^{\dagger}\mathcal{Y}_{L}E_{a}\right)\underline{\Gamma}=\underline{{\bf M}}\left(E_{a}^{\dagger}\left(\mathcal{Y}_{\bar{G}}+\mathcal{Y}_{\bar{B}}\right)E_{a}\right)$. The Hermitian part (termed ${\bf H}$) is

[TABLE]

where the matrices of (73b) commute since the product of Hermitian matrices is also Hermitian. We define ${\overline{G}}_{ij}\equiv jG_{ij}{\overline{T}}_{4}$, where $G_{ij}$ is the scalar line conductance connecting buses $i$ and $j$, and ${\bf u}_{i}\!\subset\!{\bf U}_{b}$ is the voltage element of ${\bf U}_{b}$ associated with bus $i$. From (73c), the quadratic power is

[TABLE]

The quadratic quantity ${\bf x}^{\dagger}{\overline{G}}_{ij}{\bf x}=\epsilon$ may be negated by conjugating the input (proof trivial): ${\bf x}^{*\dagger}{\overline{G}}_{ij}{\bf x}^{*}\!=\!-\epsilon$. By taking ${\bf U}_{b}^{*}$ as an input to ${\rm Re}\{\underline{\bf M}{\bf U}_{b}^{*\dagger}(E_{a}^{\dagger}\mathcal{Y}_{G}E_{a}){\bf U}_{b}^{*}\}$, by (74c) we thus have $\sum_{i,j\in\mathcal{E}}{\bf u}_{ij}^{*\dagger}\overline{G}_{ij}{\bf u}_{ij}^{*}=-P_{r}^{\star}$. ∎

We note that ${\bf U}_{b}$ cannot be chosen arbitrarily; it must represent a valid solution to the linear system of (29a). Since ${\bf U}_{b}$ is not itself a degree of freedom but rather a response to some current injection, statements about the mathematical characteristics of ${\bf M}\mathcal{Y}_{N}\Gamma+({\bf M}\mathcal{Y}_{N}\Gamma)^{\dagger}$ are difficult to prove using energy-based arguments. Despite this fact, the following two theorems offer useful insights into these characteristics. First, we offer a definition of an associated system.

Definition 5**.**

Consider classical power system $\Sigma_{c}$ with ZIP loads, passive generators, and a lossy transmission network. This system’s voltage vector ${\bf V}_{b}$ is ${\bf V}_{b}={\bf V}_{s}^{\frown}{\bf V}_{ns}$, where ${\bf V}_{ns}=-{\mathcal{Y}}_{4}^{-1}{\mathcal{Y}}_{3}{\bf V}_{s}$ from (29h). All voltage vectors are transformed such that ${\bf U}=\Gamma^{-1}{\bf V}$. The Kron reduced admittance (DWE) seen by a FO source is ${\mathcal{Y}}_{N}$ as in (30). The matrix ${\bf N}_{c}=\frac{1}{2}({\bf M}\mathcal{Y}_{N}\Gamma)+\frac{1}{2}({\bf M}\mathcal{Y}_{N}\Gamma)^{\dagger}$, with matrices $\bf M$ and $\Gamma$ from Corollary 3, will have two eigenvalues: $\lambda_{1}$ and $\lambda_{2}$.

Theorem 3**.**

Consider ${\bf N}_{c}$ from $\Sigma_{c}$:

• $(a)$

If condition (72) is met, $\lambda_{1}$, $\lambda_{2}>0$.

• $(b)$

If ${\rm Re}\{{\bf U}_{s}^{\dagger}{\bf M}{\bm{\mathcal{I}}}^{\prime}\}>0$ is measured, i.e. the DEF method has failed, condition (72) is violated and $\lambda_{1}\lor\lambda_{2}<0$.

Proof.

By modifying (71), we have $-{\bf U}_{s}^{\dagger}{\bf M}{\bm{\mathcal{I}}}^{\prime}\!=\!{\bf U}_{b}^{\dagger}\underline{{\bf M}}(\mathcal{Y}_{S}^{\prime}\!+\!E_{a}^{\dagger}\mathcal{Y}_{L}E_{a})\underline{\Gamma}{\bf U}_{b}$, and with condition (72) met, $-{\bf U}_{s}^{\dagger}{\bf M}{\bm{\mathcal{I}}}^{\prime}>0$. Since $-{\bm{\mathcal{I}}}^{\prime}={\mathcal{Y}}_{N}{\bf V}_{s}$, we have that ${\bf U}_{s}^{\dagger}{\bf M}\mathcal{Y}_{N}\Gamma{\bf U}_{s}>0$ by substitution, implying that $\lambda\{({\bf M}\mathcal{Y}_{N}\Gamma)+({\bf M}\mathcal{Y}_{N}\Gamma)^{\dagger}\}>0,\;\forall\lambda\in\{\lambda_{1},\lambda_{2}\}$, proving proposition $(a)$.

If ${\rm Re}\{{\bf U}_{s}^{\dagger}{\bf M}{\bm{\mathcal{I}}}^{\prime}\}>0$, then ${\bf U}_{s}^{\dagger}{{\bf N}_{c}}{\bf U}_{s}{<}0$ is directly implied. Accordingly, $\lambda_{1}\lor\lambda_{2}<0$. By the conservation argument presented in the proof of proposition $(a)$, if ${\bf U}_{s}^{\dagger}{{\bf N}_{c}}{\bf U}_{s}{<}0$, then ${\rm Re}\{{\bf U}_{b}^{\dagger}\underline{{\bf M}}(\mathcal{Y}_{S}^{\prime}\!+\!E_{a}^{\dagger}\mathcal{Y}_{L}E_{a})\underline{\Gamma}{\bf U}_{b}\}<0$ implying the violation of condition (72), proving proposition $(b)$. ∎

Our final theorem characterizes how the eigenvalues change as resistance is added to $\Sigma_{c}$.

Theorem 4**.**

Consider ${\bf N}_{c}$ from $\Sigma_{c}$. No amount of additional resistance to lines or loads can cause the eigenvalues of ${\bf N}_{c}$ to both become negative: if ${\rm det}({\bf N}_{c})\geq 0$ then ${\rm trace}({\bf N}_{c})\geq 0$.

Proof.

Consider some altered version of $\Sigma_{c}$ where all resistance has been removed from the system. In this situation, $\lambda_{1},\lambda_{2}\geq 0$ according to Theorem 3.

Consider some secondary altered version of $\Sigma_{c}$ where generators have no damping (making them “lossless”). In this situation, matrix ${\mathcal{Y}}_{N}$ will be purely real (${\mathcal{Y}}_{N}\!\in\!{\mathbb{R}}^{2\times 2}$) via (30). Therefore, assuming lossy lines, $\lambda_{1}=-\lambda_{2}\neq 0$.

All real systems exist between these two alternatives. We now consider a system with nominal resistance and some value $\alpha$ which parameterizes the amount of damping in the system generators ($\alpha=0$ corresponds to no damping). When $\alpha=0$, $\lambda_{1}=-\lambda_{2}\neq 0$. As $\alpha\rightarrow\infty$, then $\lambda_{1},\,\lambda_{2}\rightarrow{\mathbb{R}}^{+}$. In between these extremes, if for some value of $\alpha>0$, there is $\lambda_{1},\,\lambda_{2}\in{\mathbb{R}}^{-}$, then it would imply that the addition of positive damping causes the system to lose the ability to dissipate the quadratic energy which the generator damping consumes. This is a contradiction, so at least one eigenvalue must always remain positive for all levels of damping and resistance. ∎

Remark 1**.**

Since the DWE ${\mathcal{Y}}_{N}$ is the admittance of the network “seen” by the source, the negative DWE $-{\mathcal{Y}}_{N}$ is the admittance “seen” by the network behind the source bus. By Theorem 4, the eigenvalues of ${\bf N}_{c}=\frac{1}{2}({\bf M}\mathcal{Y}_{N}\Gamma)+\frac{1}{2}({\bf M}\mathcal{Y}_{N}\Gamma)^{\dagger}$ can never be simultaneously negative. Stated differently, the eigenvalues of $-{\bf N}_{c}$ can never be simultaneously positive. For this reason, the admittance seen by the network behind the FO source can never be truly passive.

III-C Inference Analysis Observations

If a system-wide admittance matrix model, i.e. (26), is fully known at the time of a FO event, locating the source of the FO is trivial once PMU measurements are obtained. Such a model is often unknown or poorly known in real time, so model based approaches aren’t typically applicable. Using PMU inference methods to build the admittance model (in real time) are thus a tempting solution. Unfortunately, the system is inherently underdetermined. In a standard linear circuit, admittance may be computed as the simple ratio of complex inputs and complex outputs. This cannot be done in a MIMO system. Assuming complex input (${\tilde{\bf u}}_{v}$) and output (${\tilde{\bf y}}$) vectors for the system in (3) are given, a generalized admittance ${\mathcal{Y}}\in{\mathbb{C}}^{2\times 2}$ cannot be directly inferred since (3) represents four real linear equations while ${\mathcal{Y}}$ is specified with 8 coefficients ($a_{i},\,b_{i},\,i\in\{1,2,3,4\}$):

[TABLE]

The system is thus severely underdetermined. For a classical generator, though, this is not the case. This is an important observation, because the dynamics of a system experiencing a FO are dominated by the electromechanical response of synchronous generators.

Proposition 1**.**

The admittance ${\mathcal{Y}}_{g}$ of a classical generator can be uniquely specified as the weighted sum of six basis matrices with only four real coefficients.

Proof.

It may be observed that ${\mathcal{Y}}_{g}$ in (61) contains no complex combination of basis matrix $\overline{T}_{1}$. Therefore, ${\mathcal{Y}}_{g}$ may be written using six basis matrix coefficients $a_{i}$, $b_{i}$, $i\in\{2,3,4\}$ via $\mathcal{Y}_{g}=\sum_{i=2...4}\left(a_{i}+jb_{i}\right)\overline{T}_{i}.$ We further observe that for $(a_{2}+jb_{2})\overline{T}_{2}$ and $(a_{3}+jb_{3})\overline{T}_{3}$, the ratios of $a_{2}/b_{2}$ and $a_{3}/b_{3}$ must be equal, so $a_{2}=b_{2}a_{3}/b_{3}$. This eliminates one coefficient. Finally, we write $\gamma=\gamma_{r}+j\gamma_{i}$ and observe that $b_{2}=\gamma_{i}\sin(\delta)\cos(\delta)$ while $b_{3}+b_{4}=\gamma_{i}\sin^{2}(\delta)$ and $b_{3}-b_{4}=-\gamma_{i}\cos^{2}(\delta)$. Solving this system yields $b_{4}=+\sqrt{b_{2}^{2}+b_{3}^{2}}$ (assuming damping $D$ is positive). Therefore, ${\mathcal{Y}}_{g}$ may be specified with just four coefficients: $a_{3}$, $a_{4}$, $b_{2}$, and $b_{3}$. ∎

Corollary 7**.**

From Proposition 1, classical generator admittance ${\mathcal{Y}}_{g}$ can be specified with only four parameters. Thus, the inference problem $\min_{{\mathcal{Y}}_{g}}\{{\tilde{\bf y}}-{\mathcal{Y}}_{g}{\tilde{\bf u}}_{v}\}$ represents 4 equations and 4 unknowns and is thus a “determined” problem.

In other words, the admittance of a classical generator can be fully inferred from PMU data.

IV Practical Applications of Passivity Analysis

In this section, we expound three practical applications of the analysis presented in this paper. First, we explicitly state how to implement the DEF source location method via the proposed frequency domain methods. Second, we show how system operators can use the proposed framework to predict how the DEF method will perform in their respective networks without performing simulations. And third, we explain how the tools in this paper can be used to understand how the addition of new grid components will affect performance of the DEF method. In all further analysis, we assume

[TABLE]

IV-A Energy Flow Calculations in the Frequency Domain

In order to compute the dissipating power signal $P^{\star}$ in a system which is experiencing a FO event, we convert voltage and current PMU data into rectangular coordinates. Next, we take the Fast Fourier Transform (FFT) of these signals, thus constructing ${\tilde{\bf V}}(\Omega)=[{\tilde{V}}_{r}(\Omega),\,{\tilde{V}}_{i}(\Omega)]^{\top}$ and ${\tilde{\bf I}}(\Omega)=[{\tilde{I}}_{r}(\Omega),\,{\tilde{I}}_{i}(\Omega)]^{\top}$ at all relevant buses and lines. Finally, we evaluate these signals at the forcing frequency $\Omega_{d}$ in order to compute the dissipating power injection at each bus or flow along each line:

[TABLE]

Positive $P^{\star}$ flowing out of an element indicates the presence of a non-passive component (typically a FO source). We refer readers to Appendix B for a discussion on the difficulties of dealing with rectangular coordinates in a realistic power system. Using this procedure, a map, such as the one presented in Fig. 3, can be constructed.

IV-B Predicting Performance of DEF Source Location

The implicit goal of the DEF method is to locate the source of negative damping in the system. When resistive elements are introduced to the network, the goal becomes obfuscated because the source may appear passive, and thus, positively damped. By Remark 1 though, the negative DWE of the source cannot be physically passive. A key contribution of this paper recognizes that while generators continue to be physically passive, the source only appears555The transformed FRF will have one positive and one negative eigenvalue. passive when the DEF method fails in a classical power system.

Before system operators can consistently rely on using the DEF method for FO source location, it is important that they first test if the method will perform adequately in their respective networks. Algorithm 1 clarifies the exact analytical procedure which can be used to test if the DEF method will perform successfully for any given bus in the network at any given FO frequency. We note that this off-line analysis in not restricted to classical power systems: any system component can be incorporated so long as its dynamics can be analytically approximated and linearized. Additionally, the type of oscillation source is entirely arbitrary: the methods in Algorithm 1 predict the sign of the energy flowing from the source bus terminal as a function of how the system dynamically responds to the abstracted oscillation signal. It is thus similar to testing if condition (72) is violated.

IV-C Analysis of Additional Grid Components

Until this point in the paper, only the elements of a lossy classical power system have been considered, but it is important for system operators to consider how the addition of new elements will effect energy-based source location methods. If newly added elements are non-passive, it is also important for system operators to test the degree of allowed penetration from these elements before energy-based source location methods will fail. In this subsection, we present examples of passivity analysis applied to three non-classical elements: frequency dependent loads, droop-controlled inverter systems, and third order synchronous generators with first order Automatic Voltage Regulators (AVRs).

IV-C1 Frequency Dependent Load

Assuming power is an instantaneous function of voltage magnitude and frequency[17], we may write $P(t)$ and $Q(t)$ via

[TABLE]

where $\omega=\omega_{0}+\dot{\theta}$. Since we are interested in the linearized responses of (82a) and (82b), we evaluate their partial derivatives at equilibrium (${\rm V}_{0}$, $\omega_{0}{=1}$). Assuming sinusoidal perturbations, phasor notation yields

[TABLE]

where $\tilde{\omega}=j\tilde{\theta}$. Employing matrix ${\rm T}_{1}$ from Appendix B, we transform from polar to rectangular via $\tilde{{\bf V}}={\rm T}_{1}\tilde{{\bf V}}_{p}$ in (91):

[TABLE]

We next consider perturbations of $P=V_{r}I_{r}+I_{i}V_{i}$ and $Q=V_{i}I_{r}-V_{r}I_{i}$. Treating $\bar{V}_{i/r}$ and $\bar{I}_{i/r}$ as steady state values, we linearize and convert to phasor notation:

[TABLE]

By equating (92) and (103), ${\tilde{\bf I}}$ and ${\tilde{\bf V}}$ are directly related by

[TABLE]

Finally, we set $\alpha_{p}=\alpha_{q}=0$ to isolate the effects of frequency, and we compute eigenvalues $\lambda_{1,2}=\lambda\{{\bf M}\mathcal{Y}_{b}+({\bf M}\mathcal{Y}_{b})^{\dagger}\}$:

[TABLE]

where $\varphi=\theta-\phi$. In setting $\beta\equiv\beta_{p}\approx\beta_{{q}}$, (105) simplifies to

[TABLE]

Since power factor is usually close to unity, this case of frequency dependent load is primarily passive for $\beta>0$ because the positive eigenvalue will be far larger than the negative eigenvalue. Otherwise, if $\beta_{p}=0$, then the load behaves like a resistor (equal and opposite eigenvalues), but if instead $\beta_{q}=0$, then the load is truly passive for reasonable values of power factor.

IV-C2 Droop-controlled Inverter

We consider the dynamics of a droop-controlled inveter circuit; such circuits have the potential to dominate distributed energy resource interconnections. As given in [18], the circuit’s impedance may be stated as

[TABLE]

Associated parameters are explained in [18], but $k_{p}$ and $k_{q}$ are the active and reactive power droop coefficients, and $R_{c}$ is the virtual + coupling resistance. We seek to analyze the passivity of this system. Rather then invert (109), we consider the passivity of the impedance $\mathcal{Z}$ directly. Appendix C establishes the equivalence of passivity classifications between an admittance $\mathcal{Y}$ and its associated impedance ${\mathcal{Z}}={\mathcal{Y}}^{-1}$ for transformation matrix $\bf M$. We thus compute the Hermitian part of ${\bf M}\mathcal{Z}$:

[TABLE]

Since $\tau$, $\Omega$, and the droop coefficients are all positive, if $R_{c}$ is small, then the eigenvalues of (112) are simply

[TABLE]

and the system is effectively passive. If $R_{c}$ is not neglected, though, the eigenvalue expressions become more cumbersome. Using the parameter values suggested in [18], we plot $\lambda_{1}$ and $\lambda_{2}$ as $R_{c}$ is scaled via $R_{c}=\alpha X_{c}$, $\alpha\in\{0,0.05,0.1,0.15,0.2\}$, for the matrix ${\bf M}\mathcal{Y}+({\bf M}\mathcal{Y})^{\dagger}$ (we plot admittance eigenvalues rather than impedance eigenvalues solely for graphical clarify). Fig. 2 shows the eigenvalues. Clearly, as coupling resistance decreases in value, the inverter behaves more passively. The DEF method will thus perform more successfully when virtual + coupling resistance is minimized.

IV-C3 Third Order Generator with First Order AVR

We consider the dynamics of a third order synchronous machine with a first order ARV. The associated model [17] is stated by

[TABLE]

where ${\rm V}_{r}$ is the voltage reference set point, ${\rm V}e^{j\theta}$ is the terminal voltage phasor, $i_{d}=(e_{q}^{\prime}-e_{q})/X_{d}^{\prime}$, $i_{q}=e_{d}/X_{q}$, $e_{d}={\rm V}\sin(\delta-\theta)$ and $e_{q}={\rm V}\cos(\delta-\theta)$. Finally, $P_{e}=e_{d}i_{d}+i_{q}e_{q}$. We define output power variables $P$ and $Q$ as $P=-P_{e}$ and $Q=-(i_{d}e_{q}-i_{q}e_{d})$666In keeping with the conventions of Section II, positive currents flow from the network and into the machine..

In order to make analytical claims about the passivity of (114), we leverage the alternative DEF interpretation proposed in Appendix D. To make the analysis tractable, we choose to linearize the model about the generator’s unloaded equilibrium point (i.e. $\theta\!=\!\delta$, ${\rm V}\!=\!e_{q}$, $e_{d}\!=\!0$, etc.). Although this configuration may be uncommon, the resulting linearized dynamics are sufficiently simplified for approximating generator passivity. Additionally, the resulting admittance matrix $\mathcal{H}$ is diagonal, so the eigenvalues of ${\bf K}^{\dagger}\mathcal{H}+({\bf K}^{\dagger}\mathcal{H})^{\dagger}$ are trivial. We compute $\mathcal{H}$ from Appendix D:

[TABLE]

Taking the eigenvalues of ${\bf K}^{\dagger}\mathcal{H}+({\bf K}^{\dagger}\mathcal{H})^{\dagger}$ is equivalent to taking the imaginary parts of the diagonals of matrix $\mathcal{H}$. The relationship between $\tilde{P}$ and $\tilde{\theta}$ (i.e. $\lambda_{1}$) is thus passive if

[TABLE]

This will be true if damping $D$ is positive. Interestingly, when the generator is unloaded, the damping provided by the field winding is of second order and thus has no mathematical impact on passivity between $\tilde{P}$ and $\tilde{\theta}$. Via (117), the relationship between $\tilde{Q}$ and $\tilde{\rm V}$ (i.e. $\lambda_{2}$) is thus passive if

[TABLE]

This result has many interesting interpretations. For $K_{a}$ small or $T_{a}$ sufficiently large, the expression will always be passive. Strong and fast AVR response thus causes the relationship between $\tilde{Q}$ and $\tilde{\rm V}$ to lose passivity. For passivity to be guaranteed $\forall\Omega$, the AVR gain should have the upper bound

[TABLE]

Of course, the claim that the generator will be passive if $D\geq 0$ and (120) are both satisfied is only true when the generator is unloaded. Loading effects are likely to be small, though, so these results are useful for guiding intuition related to when generators will lose passivity.

V Test Results

In this section, we present test results which illustratively validate the framework presented in Sections II and III and the practical applications presented in Section IV. We perform these tests in two systems which are altered to engender poor DEF performance. In the 39-bus New England system, the average R/X transmission line ratio was increased from $6\%$ to $15\%$, and in the 179-bus WECC system, all constant power loads were converted into constant impedance loads. All simulation code is posted online777https://github.com/SamChevalier/Passivity-Enforcement-FOs for open source access.

V-A New England 39-bus Test System

Fig. 3 shows a diagram of the IEEE 39-bus New England system. In this system, generators are modeled as third order synchronous machines [19] with first order AVRs888The gains and time constants of these AVRs were chosen to approximate the full regulator + exciter system model from [19].. Initially, loads were modeled as constant power (PQ) and system lines were modeled with no-loss (purely reactive). Following the steps in Algorithm 1, we computed the eigenvalues associated with bus 31’s DWE ${\mathcal{Y}}_{N}$ for a frequency of $\Omega_{d}=2\pi\times 2$:

[TABLE]

Due to the two positive eigenvalues, any FO source originating at bus 31 with a frequency of 2 Hz will appear non- passive: the DEF method should perform accurately in this system since condition (72) is satisfied. We then simulated the response of the system to a 2Hz perturbation originating at bus 31. System excitation and simulation were performed in the frequency domain. From the simulated data, dissipating power $P^{\star}$ was computed via (81) on all lines. The resulting flows are plotted in Fig. 3. Clearly, the source is readily identifiable.

Next, the model was altered such that the average999One outlier line, which has $R\approx X$, was excluded from this average. R/X ratio was increased from $6\%$ to $15\%$ for all lines. Additionally, all loads were converted to the model described by (82) with $\alpha$ and $\beta$ parameters chosen randomly from ${\mathcal{U}}(0,2)$. We again considered the eigenvalues associated with bus 31’s DWE ${\mathcal{Y}}_{N}$ for a frequency of $\Omega_{d}=2\pi\cdot 2$:

[TABLE]

The oppositely signed eigenvalues show that the network can aggregately generate or consume dissipating power $P^{\star}$, indicating unreliable DEF results. After computing the system’s response to a 2Hz sinusoidal FO applied at bus 31’s generator, the dissipating power $P^{\star}$ flowing in the network was plotted in Fig. 4. Based on arrow directionality, all generators are shown to be $P^{\star}$ sinks: there is no apparent generator source in the system. Loads present only a small contribution of dissipating power. Although it may appear as though bus 6, for example, is a source, the dissipating power flows are exactly conserved at this bus: all dissipating power flowing in from lines 7-6, 5-6, and 11-6 exactly flows out on line 6-31.

The reason there are no apparent sources is due to the highly active nature of the lines. For example, in the flow from bus 11 to bus 6, the sending end $P^{\star}$ is larger than the receiving end $P^{\star}$ (based on arrow sizes), and given the flow direction, the line is clearly a source of dissipating power. Similarly, on line 21-22, both arrows point away from the line, indicating positive dissipating power flows out of both ends of the line. Becasue condition (72) is not met in this network, the FO source is able to act as a dissipating power sink. In other words, since the network is producing more dissipating power than it can consume, the source provides additional “sink” slack. Therefore, the source cannot be readily identified by a system operator.

As a second test, we reconsidered the no-loss system whose eigenvalues are characterized by (121). At load buses 4, 8 and 16, we added identical droop-controlled inverter circuits whose parameter values are specified in [18]. However, we intentionally neglected resistive losses (i.e. $R_{c}=0$). Next, we scaled the nominal droop gain values of $k_{p}=1.3\cdot 10^{-3}$ and $k_{q}=7.5\cdot 10^{-3}$ by scalar $\alpha\in[0\,\;10^{4}]$. For each system configuration, the eigenvalues of the DWE seen at the source bus were computed per Algorithm 1. The results, shown in Table I, confirm the prediction of (113). Across the full range of droop gain values, the eigenvalues of the source bus are both positive. Therefore, these components do not interfere with DEF performance. When losses are added back in, interference may occur.

As a final test on the 39-bus system, we again reconsidered the no-loss system whose eigenvalues are characterized by (121). Instead of adding loss, we tripled all AVR gain values and decreased all AVR time constants by one third. Accordingly,

[TABLE]

As evidenced by the slightly negative eigenvalue, these changes caused generator passivity to be lost in at least some cases. This effect was predicted by (120).

V-B WECC 179-bus Test System

In the second test, we employed the 179-bus WECC system which was prepared by the IEEE Task Force on FOs [20]. Specifically, we employed test case number F1, where a 0.86 Hz FO is applied to the reference signal in a generator’s AVR. We altered the system by converting all loads to constant impedance with some frequency dependence. We also applied Ornstein-Uhlenbeck noise on the load power terms $P_{0}$ and $Q_{0}$ from (82a)-(82b). Due to the load model modification, the natural frequencies of the system changed slightly. Thus, in order to engender better resonance, we increased the frequency of the FO by 0.1Hz, and we increased the damping at the source generator to ensure it acted as a FO sink in the DEF analysis. We simulated the system for 120 seconds with Power System Analysis Toolbox (PSAT) [17]. To show the strength of the resonance condition, we initially simulated the system with no load noise and plotted the time domain frequency response $\omega$ of the system generators. Panel $({\bf a})$ of Fig. 5 shows that there are multiple generators with larger frequency oscillations than the source. Panel $({\bf b})$ shows the frequency dynamics of the system generators once load stochasticisticy is added.

Before investigating the dissipating energy flows in this large system, we first sought to further experimentally validate the propagation framework proposed in Section II. To do so, we analytically constructed the full system model of (26). As predicted by this model, the complex current injection vector ${\bf J}\!\in\!{\mathbb{C}}^{2n\times 1}$ from (27) should be sparse: non-zero elements should only come from unmodeled, extraneous perturbations. We thus simulated the WECC system with noisy loads and a strong FO source at bus 4. We then collected all voltage PMU data, generated the FFT evaluated at $\Omega_{d}=2\pi\times 0.96$, constructed ${\bf V}_{b}$ from (27) and the methods outlined in Appendix B, and predicted the current injection vector $\bf J$ via ${\bf J}=-\mathcal{Y}_{B}{\bf V}_{b}$. Since there are two complex current injection components associated with each bus, we defined injection magnitude vector $\overline{\bf J}\in{\mathbb{R}}^{n\times 1}$. Its $i^{\rm th}$ element is given as the sum of the magnitudes of the current injections at the $i^{\rm th}$ bus:

[TABLE]

We plot the injection magnitudes in Fig. 6. As expected, the primary current injection is located at the source bus, but there is a plethora of small, non-zero injections as well. These come from the random load perturbations which are applied at each time step. Fig. 6 thus serves to validate (26) as a useful model for analyzing FO propagation.

As stated, the FO source is located in the generator at bus 4. By converting all loads to constant impedance (with some frequency dependence), we were able to engender a system whose DWE ${\mathcal{Y}}_{N}$ had oppositely signed eigenvalues. Via Algorithm 1,

[TABLE]

Again, these are the eigenvalues seen by the source bus when it looked into the system. Because one of the eigenvalues is negative, the source generator is capable of acting as a dissipating power sink. When the AVR oscillation was initially applied, simulation results showed that the generator acted like a source. When we simply increased the damping of the source generator though, we began to excite alternative eigenvalues of ${\mathcal{Y}}_{N}$, and the source generator became a sink. For a sufficiently high level of source damping, the dissipating power injections of Fig. 7 were observed at the generator buses. In this scenario, the DEF method fails due to the high level of resistance in the system. This unreliable DEF performance is predicted by the oppositely signed eigenvalues of (125).

It is instructive to note that the eigenvalues of (125) do not change as the damping at the source bus is altered. These eigenvalues are a product of the network, not the source. This is further confirmation that Algorithm 1 is agnostic to the type or cause of the FO; rather, Algorithm 1 predicts the network response to any perturbation originating from the selected source bus. Accordingly, frequency and topological location are the only important characteristics of the FO considered in the algorithm.

VI Conclusion and Future Work

In this paper, we have proposed a novel framework for understanding how FOs linearly propagate in an electrical power network. Using Tellegan’s theorem and passivity concepts, we leveraged this perturbative propagation framework to further motivate and analyze the successful DEF method. We proved that the DEF’s shortcomings in a classical power system cannot be avoided by simply selecting a new passivity transformation (quadratic energy function), and in Theorem 4, we showed that the effective source admittance associated with a FO cannot be passive in a classical power system. We further developed necessary (and by extension, sufficient) conditions for the failure of the DEF method. When power system operators have access to an analytical system model, they may use our proposed framework to analytically predict apriori if the DEF method might fail ($\lambda_{1}\geq 0$, $\lambda_{2}\leq 0$), will fail ($\lambda_{1},\lambda_{2}\leq 0$), or will succeed ($\lambda_{1},\lambda_{2}\geq 0$). This may be especially useful in small, microgrid networks where system architecture is known more fully and line loss ratios are much higher.

One obvious drawback to the proposed approach for the apriori testing of DEF performance is the need for an analytical system model. Load models, especially, may be highly uncertain from the perspective of a system operator. And as observed in subsection III-C, inference of the underlying (load or generator) admittances cannot be performed from a single PMU observation due to the underdetermined nature of the proplem. In future work, we hope to propose model-free and data-driven inference techniques which overcome this issue by leveraging surrounding spectral data to regularize the inference problem. Although predicting DEF performance is an important contribution in itself, future research will also focus on leveraging this framework to devise methods which enhance the source location performance of energy-based methods.

The main theorems in this paper assumed a classical lossy power system, but the underlying framework presented in Sections II and III is actually quite general and can be used to interpret a variety of linear propagation phenomena. As a further extension of this work, we hope to use our framework to investigate the passivity properties of other renewable energy sources, such as doubly-fed induction generators (DFIGs) and power electronic interfaces of HVDC lines. We also plan to use these approaches to analyze the performance of the DEF method in the case of interarea oscillations, especially those caused by the interaction of many positively damped AVRs. The performance of energy-based oscillation detection methods in this situation is still an open research question.

Appendix A

We consider the complications associated with a two port element whose current flow is not a linear function of terminal voltage differentials. To do so, we return to a standard steady state (rather then perturbative) power system model, and we consider a transmission line with complex admittance ${\tilde{y}}_{ij}\in{\mathbb{C}}^{1}$. If ${\tilde{V}}_{i}$, ${\tilde{V}}_{j}\in{\mathbb{C}}^{1}$ are the phasor voltages associated with buses $i$ and $j$, then the standard two port model associated with a tap changing transformer may be stated as

[TABLE]

where $c$ is the tap ratio and ${\tilde{I}}_{ij}$, ${\tilde{I}}_{ji}\in{\mathbb{C}}^{1}$ are the current flows. Considering the top equation, the flow $\tilde{I}_{ij}$ may be written as

[TABLE]

From (133b), we may thus interpret ${\tilde{y}_{ij}}/c$ as the series impedance between the lines, and we may interpret $\tilde{V}_{i}\tilde{y}_{ij}({1-c})/{c^{2}}$ as a shunt current injection at bus $i$ caused by the tap changer. A similar shunt injection will exist at bus $j$. In this way, any two port element whose current flow isn’t directly proportional to $\tilde{V}_{i}-\tilde{V}_{j}$ (phase shifting transformer, HVDC line, etc.) can be expressed as the sum of a linear flow term and a shunt injection term. Once these relations are appropriately converted into the perturbative system model, matrices (16) and (21) can be updated accordingly.

Appendix B

An inherent problem exists when attempting to compute (81) in a realistic power system. Because the system has natural integration effects, the phase angles of the system constantly drift; this can be seen in the (simulated or real) time series data of any system with stochastic loads. Because generator admittance can be a function of steady state rotor angle $\delta$, such as in (61), the underlying admittance will be constantly changing as system phase angle drifts. Due to the infusion of such nonlinear effects, FFT analysis of rectangular (or Cartesian) coordinate variables will yield outputs which are incompatible with the perturbative model proposed in Section II: the underlying admittances will be constantly changing as system phase angle drifts. Polar coordinate formulations, however, do not have this same drawback. As a workaround, we may introduce the concept of effective rectangular perturbation. To do so, we first introduce matrices ${\rm T}_{1}$ and ${\rm T}_{1,i}$ which linearize standard phasors $V_{r}+jV_{i}={\rm V}e^{j\theta}$ and $I_{r}+jI_{i}={\rm I}e^{j\phi}$, respectively [10, eq. (10)]:

[TABLE]

We now assume we have some admittance ${\mathcal{Y}}_{g}$ which relates polar voltage and current perturbations via $\tilde{{\bf I}}_{p}=\mathcal{Y}_{g}\tilde{{\bf V}}_{p}$. This relationship will not be affected by system phase angle drift. We thus transform the expression via

[TABLE]

where $\theta$ and $\phi$ from (140)-(147) are chosen as their respective instantaneous values at $t\!=\!0$ in the PMU data (any reference angle may be used so long as it is consistent across all transformations). Vectors $\tilde{{\bf I}}_{e}$ and $\tilde{{\bf V}}_{e}$ represent the effective rectangular perturbation vectors. Thus, (148b) is a fully consistent equation and may be used to update the $P^{\star}$ flow equation in (81).

Appendix C

Rather than dealing with the admittance $\mathcal{Y}$ from ${\tilde{\bf I}}={\mathcal{Y}}{\tilde{\bf V}}$, it may be convenient to deal with impedance ${\mathcal{Z}}={\mathcal{Y}}^{-1}$ instead. We define ${\bf N}_{\mathcal{Y}}={\bf M}\mathcal{Y}+({\bf M}\mathcal{Y})^{\dagger}$ and ${\bf N}_{\mathcal{Z}}={\bf M}\mathcal{Z}+({\bf M}\mathcal{Z})^{\dagger}$.

Theorem 5**.**

Iff ${\bf N}_{\mathcal{Y}}$ is (non)passive, ${\bf N}_{\mathcal{Z}}$ is (non)passive.

Proof.

Starting with $\tilde{{\bf I}}=\mathcal{Y}\tilde{{\bf V}}$, we multiply through by $\tilde{{\bf V}}^{\dagger}{\bf M}$:

[TABLE]

Similarly, starting with $\mathcal{Z}\tilde{{\bf I}}=\tilde{{\bf V}}$, we multiply through by $\tilde{{\bf I}}^{\dagger}{\bf M}$:

[TABLE]

Since ${\rm Re}\{\tilde{{\bf V}}^{\dagger}{\bf M}\tilde{{\bf I}}\}=\alpha\in\mathbb{R}$, then

[TABLE]

since ${\bf M}={\bf M}^{\dagger}$. Because ${\rm Re}\{\tilde{{\bf I}}^{\dagger}{\bf M}\mathcal{Z}\tilde{{\bf I}}\}={\rm Re}\{\tilde{{\bf V}}^{\dagger}{\bf M}\mathcal{Y}\tilde{{\bf V}}\}$,

[TABLE]

Thus, ${\bf N}_{\mathcal{Y}}$ and ${\bf N}_{\mathcal{Z}}$ share the same passivity classification. ∎

Appendix D

The original DEF integral, as given in [7], may be stated as

[TABLE]

where ${\rm d}{\rm V}e^{j\theta}=(\dot{{\rm V}}+j{\dot{\theta}}{\rm V})e^{j\theta}{\rm d}t$. We thus restate $W_{{\rm DE}}$ as

[TABLE]

Other publications commonly state the term $(Q{\dot{{\rm V}}}/{\rm V}){\rm d}t$ as $Q{\rm d}\ln({\rm V})$. Assuming small perturbations of (153), we may write ${\rm d}\ln({\rm V})\approx{\rm d}{\rm V}/{\rm V}_{0}=({\dot{\rm V}}/{\rm V}_{0}){\rm d}t$. We have

[TABLE]

We now suppose there is some FRF $\mathcal{H}$ which relates output oscillations $\tilde{P}$ and $\tilde{Q}$ with input oscillations $\tilde{\rm V}$ and $\tilde{\theta}$ via

[TABLE]

We define a new passivity transformation matrix

[TABLE]

We use $\bf K$ to write $\tilde{{\bf S}}=\mathcal{H}({\bf K}^{-1}{\bf K})\tilde{{\bf V}}_{p}$. The passivity of the underlying system may be tested via the following relation, where ${\tilde{\bf S}}$ and ${\tilde{\bf V}}_{p}$ are from (160) as previously defined in (91):

[TABLE]

By the proof logic presented in [11], the following conditions are equivalent after a sufficient number of perturbation cycles:

[TABLE]

If instead of $\mathcal{H}$, the dynamics of the system were described by the admittance $\mathcal{Y}$, (166) follows directly from (165b):

[TABLE]

Acknowledgment

The authors gratefully acknowledge the helpful correspondence and general guidance provided by Slava Maslennikov of ISO New England. ISONE’s efforts, under Slava’s direction, to pioneer the testing and implementation of the DEF method are to be sincerely applauded.