PMU-based Distributed Non-iterative Algorithm for Real-time Voltage Stability Monitoring

TL;DR

This paper introduces a novel distributed, non-iterative voltage stability index using PMU data, enabling real-time, secure, and low-cost monitoring of power system stability with minimal communication.

Contribution

It proposes a new distributed non-iterative VSI based on recasting power flow equations as circles, improving speed, interpretability, and robustness over existing methods.

Findings

Effective in IEEE test systems

Faster and more interpretable than existing methods

Robust to measurement noise

Abstract

The Phasor measurement unit (PMU) measurements are mandatory to monitor the power system's voltage stability margin in an online manner. Monitoring is key to the secure operation of the grid. Traditionally, online monitoring of voltage stability using synchrophasors required a centralized communication architecture, which leads to high investment cost and cyber-security concerns. The increasing importance of cyber-security and low investment cost have recently led to the development of distributed algorithms for online monitoring of the grid that are inherently less prone to malicious attacks. In this work, a novel distributed non-iterative voltage stability index (VSI) is proposed by recasting the power flow equations as circles. The online computations of VSI are simultaneously performed by the processors embedded at each bus in the smart grid with the help of PMUs and communication of voltage phasors between neighboring buses. The distributed nature of the index enables the real-time identification of the critical bus of the system with minimal communication infrastructure. The effectiveness of the proposed distributed index is demonstrated on IEEE test systems and contrasted with existing methods to show the benefits of the proposed method in speed, interpretability, identification of outage location, and low sensitivity to noisy measurements.

Tables

Table 1. TABLE II : Comparison between centralized, proposed distributed VSI and decentralized VCPIs to monitor buses 14 , 29 14 29 14,29 and 30 30 30 in IEEE- 30 30 30 bus system [ 36 ] .

VCPI calculation method	No. of PMUs required	Sensitivity to noise	Full admittance matrix	Single point failure
Centralized	$30$	Low	Required	Yes
Distributed : iterative index [26]	$30$	Low	Not required	No
Distributed : proposed index	$5$	Low	Not required	No
Decentralized	$3$	High	Not required	No

Table 2. TABLE III : The standard deviation ( σ 𝜎 \sigma ) for the proposed VSI, LTI and distributed sensitivity index [ 26 ] at bus 30 30 30 in IEEE- 30 30 30 bus system when Gaussian noise is introduced in the voltage phasor measurements.

Noise in voltage angle ($\sigma V_a$) and magnitude ($\sigma V_m$ )	Proposed index (VSI)	Local Thevenin index (LTI)	Distributed iterative sensitivity index [26]
$\sigma V_m=0.001p.u.$ $\sigma V_a={0.01}^{\circ }$, [41]	$\sigma =0.0043$	$\sigma =0.0184$	$\sigma =0.0057$ when convergence tolerance = $0.001$
$\sigma V_m=0.001p.u.$ $\sigma V_a={0.5}^{\circ }$, [39]	$\sigma =0.0058$	$\sigma =0.1942$	$\sigma =0.0084$ when convergence tolerance = $0.1$
$\sigma V_m=0.001p.u.$ $\sigma V_a={0.5}^{\circ }$, [39]	$\sigma =0.0058$	$\sigma =0.1942$	algorithm diverges when convergence tolerance = $0.001$

Equations

$$p_d=t_{d,1}\cdot v_{d,r}^2+t_{d,2}\cdot v_{d,r}+t_{d,1}\cdot v_{d,i}^2+t_{d,3}\cdot v_{d,i},$$

$$q_d=t_{d,4}\cdot v_{d,r}^2-t_{d,3}\cdot v_{d,r}+t_{d,4}\cdot v_{d,i}^2+t_{d,2}\cdot v_{d,i}.$$

$$t_{d,1}=-\sum _{k\in \mathcal{N}(d)}{g}_{k,d},t_{d,2}=\sum _{k\in \mathcal{N}(d)}(v_{k,r}{g}_{k,d}-v_{k,i}{b}_{k,d}),$$

$$t_{d,3}=\sum _{k\in \mathcal{N}(d)}(v_{k,r}{b}_{k,d}+v_{k,i}{g}_{k,d}),t_{d,4}=\sum _{k\in \mathcal{N}(d)}{b}_{k,d}.$$

$$o_p=\left(\frac{-t_{d,2}}{2t_{d,1}},\frac{-t_{d,3}}{2t_{d,1}}\right),o_q=\left(\frac{t_{d,3}}{2t_{d,4}},\frac{-t_{d,2}}{2t_{d,4}}\right),$$

$$r_p=\sqrt{\frac{p_d}{t_{d,1}}+\frac{{\left(t_{d,2}\right)}^2+{\left(t_{d,3}\right)}^2}{4t_{d,1}^2}},$$

$$r_q=\sqrt{\frac{q_d}{t_{d,4}}+\frac{{\left(t_{d,3}\right)}^2+{\left(t_{d,2}\right)}^2}{4t_{d,4}^2}}.$$

$$C_k\left(z,z^{\ast }\right)$$

$$\triangleq \left[\begin{array}{cc}A& B\\ C& D\end{array}\right]\triangleq \left[\begin{array}{cc}1& -{\gamma }^{\ast }\\ -\gamma & (\gamma \cdot {\gamma }^{\ast }-{\rho }^2)\end{array}\right],$$

$${\mathcal{C}}_p\triangleq \left[\begin{array}{cc}{A}_p& B_p\\ C_p& D_p\end{array}\right],$$

$$A_p=1,B_p=-{\left(o_p\right)}^{\ast }=-{\left(\frac{-{\mathbf{b}}_p}{2}\right)}^{\ast }={\left(\frac{{\mathbf{b}}_p}{2}\right)}^{\ast },$$

$$C_p=-\left(o_p\right)=-\left(\frac{-{\mathbf{b}}_p}{2}\right)=\left(\frac{{\mathbf{b}}_p}{2}\right),$$

$$D_p=\left(o_p\right)\cdot {\left(o_p\right)}^{\ast }-r_p^2=c_p,$$

$${\mathcal{C}}_p$$

$${\Delta }^{\ast }$$

$${\Delta }_p$$

$${\Delta }_{pq}$$

$${\frac{1}{2}\cdot \left(\frac{{\Vert {\mathbf{b}}_q\Vert }^2}{4}-c_q\right))}^2.$$

$$\text{Proposed VSI}\equiv {\Delta }_{\text{norm}}^{\ast }=\frac{{\Delta }^{\ast }}{{\Delta }_{\text{no-load}}^{\ast }},$$

$${\left(x-\alpha \right)}^2+{\left(y-\beta \right)}^2={\rho }^2,$$

$${\left|z-\gamma \right|}^2={\rho }^2,$$

$$\left(z-\gamma \right)\cdot \left(z^{\ast }-{\gamma }^{\ast }\right)={\rho }^2,$$

$$z\cdot z^{\ast }-{\gamma }^{\ast }\cdot z-\gamma \cdot z^{\ast }+\gamma \cdot {\gamma }^{\ast }-{\rho }^2=0,$$

$$\triangleq \left[\begin{array}{cc}A& B\\ C& D\end{array}\right]\triangleq \left[\begin{array}{cc}1& -{\gamma }^{\ast }\\ -\gamma & (\gamma \cdot {\gamma }^{\ast }-{\rho }^2)\end{array}\right],$$

$${\Delta }_{\text{norm}}^{\ast }=\frac{{\Delta }^{\ast }}{{\Delta }_{\text{no-load}}^{\ast }},$$

$${\Delta }^{\ast }$$

$${-\frac{1}{2}\cdot \left(\frac{{\Vert {\mathbf{b}}_p\Vert }^2}{4}-c_p\right)-\frac{1}{2}\cdot \left(\frac{{\Vert {\mathbf{b}}_q\Vert }^2}{4}-c_q\right))}^2,$$

$${\mathbf{b}}_p$$

$$c_p$$

$${\mathbf{b}}_p$$

$$c_p$$

$${\Delta }_{norm}^{\ast }=\frac{b^4+2b^2{g}^2-4b^{2}gP-4b^2{P}^2+g^4-4g^{3}P}{{\left(b^2+g^2\right)}^2},$$

$$P_{max}=\frac{-(b^2+g^2)g+(b^2+g^2{)}^{\frac{3}{2}}}{(2b^2)},$$

$$\text{No load condition:}{\Delta }_{norm}^{\ast }=\frac{b^4+2b^2{g}^2+g^4}{{\left(b^2+g^2\right)}^2}=1,$$

$$\text{Full load condition (P=Pmax):}{\Delta }_{norm}^{\ast }=0.$$

$$\displaystyle\eqref{eq:simplified_eq_0_rev_4}\ \text{(when $x/r=5$)}\ =-\frac{P^{2}}{169}-\frac{10\,P}{13}+1\approx-\frac{10\,P}{13}+1.$$

$$\displaystyle\eqref{eq:simplified_eq_0_rev_4}\ \text{(when $x/r=10$)}\ =-\frac{4\,P^{2}}{10201}-\frac{40\,P}{101}+1\approx-\frac{40\,P}{101}+1.$$

$$\frac{d^2({\Delta }_{norm}^{\ast })}{d{P}^2}=\frac{d^2\left(-\frac{4P^2}{10201}-\frac{40P}{101}+1\right)}{d{P}^2}=-7.8424e^{-4}.$$

Full text

I Introduction

The advent of phasor measurement units (PMUs) and their wide-spread acceptance has made it possible to obtain real-time and time-synchronized information of voltage and current phasors for use in diversified applications [1, 2]. Consequently, the proliferation of PMUs transformed the power grid into a tightly integrated cyber-physical system, where local computation is viable for distributed monitoring and coordinated control actions in real-time for wide-area monitoring, protection, and control (WAMPAC).

Among WAMPAC applications, online monitoring of long-term voltage instability (LTVI) is an area of interest for industry [3]. Specifically, LTVI is a quasi-static bifurcation (nose point of the PV curve) [4], caused by the inability of the generation and transmission system to provide sufficient power to loads, e.g., due to increased demand, generation outage, or generators on VAR limits [5, 6]. If left unattended, LTVI results in a system-wide voltage collapse and blackout [7, 8]. Using PMU measurements, we can compute the margin (a measure of distance) from the current operating condition to LTVI. This measure of distance is quantified using metrics known as voltage collapse proximity indices (VCPIs). The accuracy of VCPIs depends on the theoretical basis used to derive them, e.g., Thevenin circuit, sensitivity methods [9, 10]. These methods compute the VCPIs by either decentralized(local) or centralized communication architectures.

Conventional centralized VCPIs typically use Jacobian-based sensitivity indices. These methods use Jacobian to identify LTVI since the Jacobian matrix becomes singular at LTVI [11, 12]. However, the calculation of the Jacobian requires voltage phasor information at every bus in the power system, which is expensive [13, 14, 15, 16]. In addition to this, the interpretability of sensitivity based VCPIs is harder than decentralized(local) VCPIs. It is hard because, when the operating point is close to the limit, the corresponding centralized VCPI value suddenly increases, becoming unbounded. The unbounded nature and non-linear behavior make the centralized VCPI hard to interpret as a fair distance to the critical loading. References [17, 13, 14] address the non-interpretability partially in a centralized manner using measurements from all the nodes in the system. Specifically, the original $n$-bus system is simplified into a $2$-bus multi-port coupled network equivalents. However, there is no theoretical proof that the LTVI of the full-system can always be identified using these simplifications.

To overcome the drawbacks of centralized methods, methods (indices) that utilize only local PMU measurements at a substation are developed to monitor the LTVI locally, and these indices are referred to as decentralized VCPIs [18]. Such VCPIs employ Thevenin-based methods by utilizing the quasi-steady-state behavior of the power system to calculate a $2$-bus Thevenin equivalent at a bus [18]. In such an equivalent model, one node represents the bus of interest while the other accounts for the rest of the system. The equivalent circuit estimation only requires local measurements (i.e., requires only a single PMU) at a bus of interest over a time period, making it a decentralized architecture. Once the equivalent circuit is estimated, an index value is computed at the bus of interest using maximum power transfer theorem [18]. Therefore, these Thevenin-based methods have a nice property of decentralized computation of the index values [19, 18, 20, 21]. It is also shown that the local index can theoretically (in the presence of no measurement noise) always identify the LTVI of the full system by correlating it with the Jacobian [11]. Interpretability is an additional advantage of decentralized VCPIs, as the values are usually between [math] and $1$. Unfortunately, the decentralized method suffers from errors in practice, especially due to the presence of measurement noise, e.g., the local VCPIs can be inaccurate with large errors [22].

As the centralized and decentralized methods have trade-offs between accuracy, cost, and interpretability, we are interested in an intermediate scheme to include these advantages without minimal trade-offs. Specifically, the centralized method requires measurements from the entire system, while the decentralized method only requires the PMU measurement on the bus of interest. In this paper, we propose to use an intermediate approach that requires PMU measurements from the neighboring buses of the interested bus to calculate the proposed VCPI. This set up is known as a distributed monitoring scheme, as shown in Fig. 1. In this scheme, PMUs communicate their measurements at a bus in the physical layer with its adjacent buses via the cyber layer. The processor embedded at each bus only uses its adjacent buses’ voltage phasor measurements to compute the proposed VCPI referred to as voltage stability index (VSI).

The structure of the power system in Fig. 1 is for a future grid communication architecture. Some of these communication advantages are robustness to cyber-attacks on the communication nodes, and improved reliability of communications due to shorter communication links [23, 24]. Hence it is highly desired to have a distributed communication architecture due to its various advantages in contrast to star/centralized architecture. Hence it has been of interest in the community recently for distributed implementation of various applications. For example, WAMPAC [25, 26, 27], load shedding [28, 29], optimal power flow [30], economic dispatch [31], transfer capability assessment [32], and control techniques in microgrids [33, 34]. Even for VCPI, [26] uses such a distributed communication setup, where the Jacobian is used to calculate a sensitivity index iteratively. Unfortunately, the advantage of the trade-off between centralized method and the decentralized method is not utilized fully, making [26] suffer from $1$) necessity for PMUs to cover all buses, $2$) increasing computational time as the power system operating point gets closer to critical loading, $3$) non-interpretable nature of the index, $4$) theoretical flaw of assuming eigenvalues of the Jacobian to be non-negative (invalid in the IEEE $300$-bus network).

Contributions: In this paper, we propose a novel distributed index that overcomes the drawbacks of [26] mentioned above. Instead of exploiting the properties of the Jacobian, we analyzed the power flow equations themselves with a new perspective as circles in a distributed framework. This leads to a new voltage stability index (VSI) that $1$) does not need PMUs at all buses $2$) is non-iterative leading to fast computation irrespective of system loading, $3$) is interpretable (normalized between $1$ and [math]), where $1$ indicates no-load and [math] indicates maximum load, $4$) is exact (assumption free) unlike other indices in literature [18, 17, 13, 14]. This is the first time that a non-iterative, approximation free as well as distributed communication based method has been proposed for monitoring the LTVI.

The paper is organized as follows: Section II explains the power flow equations as circles in a distributed framework. Section III quantifies the distance between the power flow circles as a measure of margin to LTVI. Section IV shows the efficacy of the proposed scheme and VSI via simulations on different IEEE systems. Section V concludes the paper.

II Power Flow Equations as Intersection of Circles

Let $p_{d}$ and $q_{d}$ be the active and reactive power injections at bus $d$. Let $\mathit{g}_{k,d}+j\cdot\mathit{b}_{k,d}$ be the $(k,d)^{th}$ element of the admittance matrix $Y$. $\mathit{v}_{d,r}$ and $\mathit{v}_{d,i}$ be the real and imaginary part of the voltage phasor at bus $d$, respectively. We represent bus $d$’s neighboring bus set as $\mathcal{N}(\mathit{d})$. The power flow equations used in this paper are in rectangular coordinate form. They are given by

[TABLE]

The parameters $t_{d,1},t_{d,2},t_{d,3}$ and $t_{d,4}$ are given by

[TABLE]

II-A Distributed Nature of Power Flow Equations as Circles

The advantage of using power flow equations (1) comes from their visualization as circles. Specifically, for fixed constants $t_{d,1},t_{d,2},t_{d,3},t_{d,4}$, equations (1a) and (1b) represents two circles at bus $d$ in $v_{d,r}$ and $v_{d,i}$ space [12]. The real power equation (1a) at bus $d$ can be represented as a circle with center ${o}_{p}$ and radius $\mathit{r}_{p}$. The reactive power equation (1b) at bus $d$ can be represented as a circle with center ${o}_{q}$ and radius $\mathit{r}_{q}$. These centers and radii are given by

[TABLE]

Using the centers and radii (3); the intersection points of real and reactive power circles at bus $d$ provides the voltage solution $v_{d}$ at bus $d$ where $v_{d}=v_{d,r}+j\cdot v_{d,i}$. It is important to note that the expressions for the centers and radii of the power flow circles at bus $d$ contain $t_{d,2}$ and $t_{d,3}$ which contain the neighboring bus voltages $v_{k,r}+j\cdot v_{k,i}$ where $\mathit{k}\in\mathcal{N}(\mathit{d})$. To plot the power flow circles at bus $d$, the centers and radii in (3) are calculated locally by the embedded processor at bus $d$ using the voltage phasor measurements from the PMUs located at the immediate neighboring buses to bus $d$. For example, to plot the power flow circles at bus $3$ in Fig. 2, we only need the voltage phasors from the PMUs at buses $2$ and $4$; branch admittance of lines joining bus $3$ to $2$ and $4$. Fortunately, in the proposed distributed scheme, we do not need the admittance matrix of the entire system.

The discussion above is for PQ buses with constant real and reactive power constraints. In the case of PV buses, we get constant real power and voltage constraints, which also lead to two circles, namely the real power circle (1a) and voltage circle centered at the origin with radius $\left|V_{\textrm{specified}}\right|$. We will next discuss how the power flow circles are impacted due to load increase.

II-B Impact of Load Increase on Power Flow Circles

The long-term voltage instability is caused when the generation and transmission system cannot supply power demanded by the load. This situation corresponds to the lack of a solution for the power flow equations. This can be identified by looking at the intersection of power flow circles constructed using PMU measurements at bus $d$. Geometrically, two circles can have two, one or no points of intersections. Hence, the voltage solution at bus $d$ represented as an intersection of power flow circles can have (a) two common points to indicate multiple feasible voltage solutions, or (b) one common point to indicate a single feasible voltage solution corresponding to the nose point (LTVI) of PV curve, or (c) no common points to represent the in-feasibility of operating conditions.

Before we provide an example, here we present a high-level data-flow of any voltage stability monitoring method and explain the overall procedure for the calculation of a VCPI. The sequential steps for any online stability monitoring method are: (a) Collection of real-time PMU measurements from the grid, (b) Estimation of power system states from the PMU measurements, and (c) Utilizing the estimated states to calculate an index that provides a metric of the system stability (VCPI). At a given time step or “snapshot”, we are monitoring the voltages in the power grid directly from the PMUs. When there is a load change in the power grid changes in the next snapshot, the PMU voltage measurements in the next snapshot are automatically updated since PMU device voltages are functions of loading on the power grid. Using these real-time PMU measurements directly the proposed index is calculated. Alternatively, the input measurements to calculate the proposed index can also be obtained from the output of a preprocessor such as a state estimator.

For illustration purpose, we use a completely connected $3$-bus system with the same branch admittance value of $1-j\cdot 0.5$ p.u. To showcase the behavior of power flow circles due to the impact of load change, we will increase the loading on the power grid until there is a blackout and observe the power flow circles’ behavior. We also assume the same apparent load at bus $2$ and $3$, e.g., $S_{2}=S_{3}$ for simplicity. The loads at bus $2$ and $3$ are increased in $5$ different time instants with the following values (in p.u.) $[-0.01+j\cdot 0.33,-0.04+j\cdot 0.40,-0.13+j\cdot 0.44,-0.28+j\cdot 0.45,-0.49+j\cdot 0.43]$ until the power flow circles at bus $3$ had only one common point. To calculate the power flow circles’ radius and center corresponding to the $5$ different time instants in real-time, we use the distributed communication architecture described in Fig. 2 to directly collect the real-time PMU voltage measurements for the $5$ different time instants. Using these neighboring bus PMU voltages and (3), the circles can be drawn in real-time. The corresponding power flow circles’ intersection is plotted as shown in Fig. 3, where the solid and dashed circles indicate the reactive and real power circles respectively.

In Fig. 3, the common points “A”, “B”, “C”, “D”, and “E” indicate the higher magnitude voltage solutions corresponding to an apparent power load of $[-0.01+j\cdot 0.33,-0.04+j\cdot 0.40,-0.13+j\cdot 0.44,-0.28+j\cdot 0.45,-0.49+j\cdot 0.43]$ p.u. respectively at the bus $3$. It can be observed that the reactive power circle steadily becomes smaller, and the real power circle steadily becomes larger as the system load increases. This behavior of power flow circles shows that the circles move further away from each other until they have only one common point (nose point) as the load increases to the critical value. This occurs in Fig. 3 at the intersection of black power flow circles. It corresponds to the critical loading in the system and has only one intersection point “E”. This point represents the nose point of the PV curve. It is to be noted that all the power flow circles at every bus will not touch each other externally at a given snapshot, but only the power flow circles at the bus that is critically loaded will touch externally indicating the LTVI.

Thus, the distance between the power flow circles at the critical bus decreases to zero as the load increases to the critical value. Hence, this distance can be used as an indicator of LTVI. Quantifying the distance between the circles for deriving a VSI is described in the next section.

III Characterizing the Distance Between Power Flow Circles for VSI

The distance between the power flow circles constructed using the PMU measurements indicates the margin to reach LTVI. There are several ways to quantify the distance between two circles. One such indicator may consider the distance between the circles using their centers and radii. This distance is positive, zero and negative when the circles intersect, touch externally and do not intersect with each other, respectively. Another indicator is the Euclidean distance between the two common points when the circles intersect each other. When the circles intersect, touch externally, and do not intersect with each other, the Euclidean distance between the two voltage solutions is positive, zero, and does not exist, respectively. We observed that the proposed index (VSI) derived below in contrast to the above mentioned possible formulations of distributed indices has a relatively better behavior with changing loading.

Different than the indicators mentioned above, the proposed index is derived by using the concept of the family of circles instead of directly using (1a) and (1b). In this section first, we introduce the concept of determinant to infer the intersection of any two circles. Second, using the determinant, we derive the expression for the proposed distributed voltage stability index.

III-A Identifying the Intersection of Power Flow Circles Using Determinant

Given real and reactive power circles, two important aspects have to be inferred to derive the proposed indicator. They are 1) the existence of a feasible solution, i.e., does the power flow circles intersect? and 2) if there is a feasible solution, then what is the distance to voltage collapse point?. These two aspects are analyzed by representing the power flow circles as a family of circles, and this analysis is facilitated by using the determinant concept [35].

To study (1) and (2) aspects mentioned in the above paragraph, we need to analyze how the real and reactive power circles interact together i.e., family of circles (parameter of circles). Given two circles $C_{1}$ and $C_{2}$, their family of circles $C_{f}$ is represented by the equation $C_{f}=\lambda_{1}\cdot C_{1}+\lambda_{2}\cdot C_{2}$ $\forall\ \lambda_{1},\lambda_{2}\ \in\mathbb{R}$. For a given $(\lambda_{1},\lambda_{2})$ pair, $C_{f}$ represents a single circle that passes through the common points of circles $C_{1}$ and $C_{2}$. Fig. 4 shows a typical illustration of family of circles $C_{f}$ with radical axis. In Fig. 4, the radical axis is the line that passes through the common points of circles $C_{1}$ and $C_{2}$. Table I presents some important geometrical properties of $C_{f}$.

From Tab. I, [35] shows that the determinant can be used to study the property of symmetry of $C_{f}$ around the radical axis. We can see from Tab. I that the symmetry around the radical axis can directly indicate the common points between the two given circles. The following provides the inferences for the approach selected in this paper.

For example, when the power flow circles intersect, then we have common points between them and vice versa. From Table I, the mirror image (symmetry) around the radical axis is a good indicator to identify if the power flow circles intersect. Mathematically, property of symmetry around the radical axis in Fig. 4 is indicated by the determinant value of $C_{f}$ [35]. For example, a non-negative determinant value of $C_{f}$ indicates the symmetry around radical axis, i.e., power flow circles intersect, and thereby a feasible solution exists. Similarly, a negative determinant value of $C_{f}$ indicates that the power flow circles do not intersect. In Section III-B, we show that family of circles $C_{f}$ corresponding to real and reactive power circles can be represented in a square matrix form. This square matrix form enables us to calculate its determinant.

III-B Proposed Voltage Stability Indicator

In this subsection to derive the proposed VSI using determinant, first at given a loading condition, the power flow circles are computed locally at a bus using PMU voltage phasor measurements from its adjacent buses and these circles are then represented as square matrices to enable the determinant calculation. Second, we derive the proposed VSI by utilizing the properties of determinant inferred in Section III-A.

Equation (4) represents the standard homogeneous form of a circle ($C_{k}$) with locus $z=x+j\cdot y$ and $\mathit{z^{*}}$ is complex conjugate of $z$. The center and radius of the circle $C_{k}$ are denoted by $\gamma=\alpha+j\cdot\beta$ and $\rho$ respectively. Reference [35] shows the idea of representing circle as a hermitian matrix as shown in (5). The derivation from (4) to (5) is provided in Appendix -A.

[TABLE]

where $\mathit{A}\neq 0$. $\mathit{A}$ and $\mathit{D}$ are always real, $\mathit{B}$ and $\mathit{C}$ are always complex conjugates for a circle. Since the power flow equations are also circles, they can be represented as hermitian matrices and such matrix representation of power flow equations enable the determinant calculation which helps to derive the proposed index.

Matrix Representation of Power Flow Equations and Proposed VSI: Once we obtain the matrix form of power flow equations, its determinant is useful to 1) identify the existence of power flow solution and 2) distance to voltage collapse point. Similar to the above calculations, the elements of real power circle-matrix $\mathcal{C}_{p}$ are calculated using its center (3a) and radius (3b) as shown below.

[TABLE]

where $\mathbf{b}_{p}=\begin{bmatrix}\dfrac{t_{d,2}}{t_{d,1}}&\dfrac{t_{d,3}}{t_{d,1}}\end{bmatrix}^{T}$, $\mathbf{b}_{q}=\begin{bmatrix}\dfrac{-t_{d,3}}{t_{d,4}}&\dfrac{t_{d,2}}{t_{d,4}}\end{bmatrix}^{T}$, $\mathit{c}_{p}=\dfrac{-\mathit{p}_{d}}{t_{d,1}}$ and $\mathit{c}_{q}=\dfrac{-\mathit{q}_{d}}{t_{d,4}}$. Similarly, the reactive power circle-matrix $\mathcal{C}_{q}$ is also shown below.

[TABLE]

Finally, the power flow equations are represented as matrices $\mathcal{C}_{p}$ and $\mathcal{C}_{q}$ respectively. As discussed in Section III-A, to study the voltage collapse phenomenon one must analyze the real and reactive power circles together i.e., family of circles ($\lambda_{1}\cdot\mathcal{C}_{p}+\lambda_{2}\cdot\mathcal{C}_{q}$). Using the inferences and property of determinant described in Section. III-A, at a given loading condition $C_{f}$ is given by $\mathcal{C}_{p}+\mathcal{C}_{q}$ and the expression for its determinant $\left|C_{f}\right|$ (say $\Delta^{*}$) is given by (7).

[TABLE]

where

[TABLE]

We determine $\Delta^{*}$ by calculating $\Delta_{p}$, $\Delta_{q}$ and $\Delta_{pq}$ at a given loading condition using the centers and radii. The power flow circles touch each other at the voltage collapse point (nose point of the PV curve). In such a scenario, Tab. I shows that the determinant $\Delta^{*}$ should become zero. Hence the lower bound of $\Delta^{*}$ is zero, which occurs at the nose point of the PV curve (externally touching circles). However, to interpret the results of the proposed index, it should range between $1$ and [math]. In order to bound the upper limit of $\Delta^{*}$, we normalize it using the no-load value of the index. The normalized form of $\Delta^{*}$ is the proposed VSI, as shown below.

[TABLE]

where $\Delta^{*}_{\textrm{no-load}}$ is the value of (7) with zero load. Appendix -B provides an explanation of the near-linear behaviour of the proposed index for a $2$-bus system.

Remark 1**.**

(VSI for PV Buses): (8) is the distributed VSI for a PQ bus, and a similar VSI can be derived for a PV bus. For PV buses, the voltage solution is determined by the intersection of voltage and real power circles, and the distance between them indicates the VAR limit violation. These circles are again computed using the PMU voltage phasor measurements of adjacent buses.

IV Simulations and Discussion

The proposed distributed non-iterative VSI is tested on various test cases such as IEEE $30$, $300$ and $2383$-bus systems [36]. To obtain the voltage phasor measurements and true voltage stability margin, we use an adaptive Newton-Raphson step size based power flow solver known as continuation power flow (CPF) [37] from MATPOWER [38]. CPF calculates the voltage solution on the PV curve by increasing the load and generation along a load/generator increase direction parametrized by continuation parameter $\lambda$ [37].

In this section, we show that the proposed distributed VSI is advantageous compared to both centralized and decentralized VCPIs. We also compare the proposed distributed index with the only other distributed method [26] in the literature to show its merits. Finally, we show the behavior of the proposed index for line & generator outages.

IV-A Non-iterative Distributed VSI Versus Centralized VCPIs

The power system operator generally monitors the power system from a control room. In practice, there are specific regions of the power system that are critical for LTVI and can be identified by offline means. Thus, the operator is interested in online LTVI monitoring in a few strategic locations/regions. For instance, in the IEEE-$30$ bus system, the operator would like to monitor the buses $14,29$ and $30$. To calculate centralized VCPIs at buses $14,29$ and $30$, the voltage measurements at all the buses in the power system are required, i.e., full observability. To calculate the proposed distributed index at buses $14,29$ and $30$, the voltage measurements at only their immediate neighboring buses are required. The adjacent nodes connected to buses $14$, $29$ and $30$ are buses $(12,15)$, $(27,30)$ and $(27,29)$ respectively. Thus we only require $5$ PMUs at buses $12,15,27,29$ and $30$ to calculate the proposed distributed VSI.

Additionally, the proposed VSI does not require the complete knowledge of the system admittance matrix. Instead, it takes advantage of the sparse nature of power system graph by using only branch admittance values of branches connecting to buses $14,29$ and $30$ to calculate the proposed distributed index. Table II presents the requirements of centralized, decentralized, and distributed methods such as proposed index, existing iterative index [26] to monitor buses $14,29$ and $30$. It is important to note that even though the iterative distributed index proposed in [26] uses communication links between the neighboring buses, [26] mandates requirement of PMUs at all buses in the grid for convergence of its iterative algorithm. Whereas the proposed method does not require PMUs to cover all the buses.

IV-B Non-iterative Distributed VSI Versus Decentralized VCPIs

In this subsection, we compared the proposed VSI with other methods (decentralized and distributed index [26]) when there is noise in the PMU measurements, and we show that the proposed VSI is less prone to noisy measurements. To understand the impact of noise on the proposed methodology, an additive Gaussian noise with zero mean and standard deviation of $0.001$ p.u. on voltage magnitude ($\sigma V_{m}=0.001p.u.$), $0.5^{\circ}$ on phase angles ($\sigma V_{a}=0.5^{\circ}$) are introduced in the measurements according to the analysis of field-tested PMUs by New England ISO [39, 40] and IEEE standard for acceptable PMU errors [41]. To demonstrate the robustness of the proposed VSI with regard to noise, it is compared to the local Thevenin index (LTI) i.e., a decentralized VCPI [18, 21, 42] and distributed sensitivity index [26]. The results are shown in Tab. III and Fig. 5 under different noise levels.

From Tab. III, it can be observed that the proposed VSI has a minimal standard deviation when compared to that of both LTI and distributed sensitivity index [26] for different noise levels. Additionally, the standard deviation of LTI increases considerably when $\sigma V_{a}=0.5^{\circ}$. However, the proposed index and distributed iterative sensitivity index [26] are not so impacted by the increase in the noise of voltage angle measurements. Despite the lower standard deviation of distributed iterative sensitivity index in Tab. III, it is important to note that the method [26] suffers from non-convergence as the noise increases. For example, when $\sigma V_{m}=0.001p.u.$ and $\sigma V_{a}=0.5^{\circ}$, the distributed sensitivity index [26] does not converge when its convergence tolerance is $0.001$ but converges for a lower level noise of $\sigma V_{a}=0.01^{\circ}$ with same tolerance setting. This creates a challenge in setting the threshold for convergence [26] when deployed in the field since the PMU measurements have variable noise levels. In contrast, the proposed index has no such drawback as it is a non-iterative index and it is also least sensitive to the noise in PMU measurements as shown in Tab. III.

From Fig. 5, LTI is more sensitive to noise than that of the proposed distributed index due to its approximation errors [18]. Similarly, from Tab. III, the distributed iterative index is more sensitive to noise than the proposed distributed index due to its iterative nature. While the proposed distributed index is the least sensitive one to noise.

IV-C Effect of PMU Voltage Phasor Data: Measurement Noise Variability, Large Measurement Outliers, and Missing Data

We further investigate the effect of bad quality phasor measurements from PMUs on the proposed index. To handle these bad quality phasor measurements, one must use a preprocessor that takes bad quality measurements like input and outputs the filtered/better measurements. These filtered measurements are used to calculate the proposed distributed index. Depending on the quality of data, this preprocessor can be a state estimator, bad data detector, low-rank matrix methods, etc. and this preprocessor is independent of the proposed methodology. Various types of errors such as clock drift errors, measurement spikes, high amount of noise, and missing data can lower the quality of the PMU data [43]. We generate PMU voltage phasor data that is very close to the real-world using the error statistics presented in [43, 39] as shown in Fig. 6. In this subsection, we use noisy PMU measurements to calculate the proposed index at buses $28,29,$ and $30$ of IEEE-$30$ bus system. Since these noisy measurements with different types of errors can directly impact the quality of the VCPIs, the proposed index is calculated by using the output measurements of a preprocessor that uses simple robust statistic such as median of the PMU measurements in the given time window (e.g., 1 sec). Fig. 6 shows that the proposed index is not affected by any of the data dropouts, measurement spikes, and the device clock errors. In contrast to [26], due to the iterations involved, it is unclear how the distributed iterative index [26] handles the different kinds of noise like missing data, time skew, measurement spikes and large white noise. Different methodology based preprocessors provide different quality of measurements. Irrespective of the preprocessor, for a given same amount of input noise signal, we show that the proposed index is robust (Tab. III, Fig. 5) and fast (Tab. IV, Tab. V) when compared to other methods [18, 21, 42, 26].

Thus, the proposed methodology is robust to various kinds of noise types and levels when compared to [18, 26, 21, 42], ensuring that the alarms triggered using this VSI will have fewer false-positive rates, providing a reliable grid monitoring scheme. For an illustration of the impact of a preprocessor on the proposed index, we present two cases i.e., 1) distributed state estimator that uses both current and voltage measurements, and 2) low-rank matrix method. 1) Fig. 7 shows a preprocessor that uses the same distributed communication architecture requirement as that of the proposed VSI can use redundant local branch current measurements to improve the noisy data and thereby increase the accuracy of the proposed method. Generally, the proposed index is robust to noise. However, in cases when there is an impact of noise, it is recommended to use a distributed state estimation based preprocessor that uses redundant branch current measurements to further improve the performance of the proposed methodology. 2) Fig. 8 shows a preprocessor that can recover large chunks of missing data where a simple median statistic will fail. Specifically, we used the OLAP-t method [44] to recover the missing data in a $2$ second window frame from a PMU located at bus $8$. In the IEEE-30 bus system, buses $6,8,27$ are the neighboring buses to bus $28$. It is important to use data only from buses $6,8,27$ to recover the missing data at bus $8$ due to the distributed communication architecture (framework) used by the proposed methodology. We used the recovered data from [44] and a simple heuristic that uses last known non-zero values to calculate the proposed index at bus $28$. Fig. 8 shows that [44] is more accurate than the simple heuristic method especially when there is more missing data.

IV-D Non-iterative Distributed VSI Versus Iterative Distributed Sensitivity Index

In this subsection, we compare the interpretability of the proposed index and iterative distributed sensitivity index [26] and show that the proposed index can provide a measure of distance to system voltage collapse. We also discuss the number of iterations taken by [26] in contrast with the non-iterative nature of the proposed index.

IV-D1 Proportional Load & Generation Increase

In this case, first, we verify the proposed VSI correctly identifies the critical bus in the system that causes the LTVI by comparing it with the only other distributed method [26] in the literature. Second, we show that the proposed VSI is easier to interpret compared to [26]. $\lambda$ (continuation parameter from CPF) is defined as the load scaling factor that scales the loads and generations in IEEE-$30$ bus system and continuation power flow (CPF) solver [37] from MATPOWER [38] is used to generate voltage phasor measurements for each $\lambda$. The voltage measurements for each $\lambda$ are used to calculate the proposed non-iterative distributed VSI. Fig. 9(a) shows the proposed distributed index values versus load scaling factor $\lambda$ for IEEE-$30$ bus system. The first observation is that at the critical loading (corresponding to $\lambda=2.8$), bus $30$ has the lowest voltage magnitude of $0.51$ p.u. in the entire network. The reason for bus $30$ to be critical is due to its location being electrically farthest away from the generators and synchronous condensers. When the load scaling factor $\lambda$ is greater than $2.8$, there is no feasible operating solution identified by CPF due to the LTVI phenomenon. Consequently, when $\lambda=2.8$, the proposed distributed index value calculated at bus $30$ is the smallest among all buses, and it is very close to 0, implying that this is the critical bus in the system that causes the LTVI.

This fact is also verified by comparing it with a distributed iterative sensitivity method [26]. Fig. 9(b) presents the index values of [26] versus load scaling factor ($\lambda$) and it also shows that the critical bus is bus $30$ i.e., the bus with maximum sensitivity index value at $\lambda=2.8$. Thus the VSI correctly identifies the critical bus in the system. However, It can be seen that the sensitivity index varies in an extremely non-linear manner with the load scaling (continuation) parameter ($\lambda$). Distributed iterative sensitivity index in Fig. 9(b) is unbounded at $\lambda=2.8$, making it hard to interpret the distance to voltage collapse point. Whereas, the distributed non-iterative VSI in Fig. 9(a) is bounded between $1$ and [math] corresponding to no-load and loadability limit, respectively. This makes the proposed index a good indicator to interpret the distance to voltage collapse point. Hence the non-interpretability and non-linear nature of the sensitivity methods makes it hard to set monitoring thresholds to reliably trigger controls while the proposed index solves these problems by effectively using the PMU measurements. A similar study for larger test case systems such as IEEE-$300$ and IEEE-$2383$ along with their critical buses are also presented in Fig. 10(a) and Fig. 10(b) respectively.

Even though the proposed index value at the critical bus (e.g., bus $30$) is expected to be [math] at critical loading ($\lambda=2.8$), we observe that the index value is a very small non-zero value (VSI = $0.03$) from Fig. 9(a). This behavior is due to the numerical instability (ill-conditioned Jacobian) of the power flow solver (CPF) from MATPOWER [38], which is used to generate the voltage phasor measurements. For example, when the loading condition $\lambda$ is very close to the nose point of the PV curve, it is well-known that the power flow Jacobian becomes ill-conditioned and power flow solvers diverge [45]. Due to this divergence behavior, it is not possible to accurately generate the voltage phasor measurements corresponding to the nose point of the PV curve, where the proposed index becomes zero. Hence we see very small non-zero values in Fig. 9(a). However, this limitation is not present when the proposed index is deployed on the field because the voltage phasor measurements are taken from the PMUs directly.

IV-D2 Computation Performance of VSI

In this case, we first compare the performance of the proposed distributed non-iterative index with that of the distributed iterative index [26] under two scenarios. They are a) increase in system loading and b) the presence of noise in voltage phasor measurements. a) In case of an increase in system loading $\lambda$, from Tab. IV, in the presence of both noise and no noise, the time is taken by the distributed iterative method [26] increases with increase in system loading $\lambda$. This behavior is not desirable since the updating frequency (calculation time) of the distributed iterative index is not constant, and it is dependent on the loading condition ($\lambda$). Fortunately, the time taken by the proposed index is constant and almost negligible for any loading $\lambda$ due to its time complexity being $O(1)$ (non-iterative index). b) In presence of noise in voltage phasor measurements, it can be seen from Tab. IV that for any system loading $\lambda$, the time taken by the distributed iterative index [26] is higher than that of the no noise scenario. To speed-up the distributed iterative index, its convergence tolerance has to be varied based on the noise level in measurements. However, it is not practical to adjust the convergence tolerance depending on varying noise levels in an online monitoring application. Fortunately, again the time taken by the proposed method is not affected by the noise level in measurements due to its approximation free formulation and non-iterative nature.

For a better comparison of the proposed and existing iterative method [26], Tab. V presents the totals floating-point operations per second (Flops) taken by both the methods. To calculate the flops, we use the information from [46] to determine the count of flop needed for performing mathematical operations such as addition, subtraction, multiplication, and division. The number of flops required to calculate the proposed index at bus $n$, which has “$M$” neighboring buses = $8\cdot M+37$. While, the total count of flop required to calculate the existing iterative method = $16\cdot M+(16\cdot M+360)\cdot I$, where “$I$” is the total iterations taken by the iterative method to converge. In case of the proposed method, since it is a non-iterative approach the required number of flops expression does not have “$I$” (number of iteration) in it.

IV-E Behavior of Proposed VSI after Line & Generator Outage

Network reconfiguration due to line outage (either due to faults or for maintenance) is a frequent occurrence in the power system. One key input to the proposed method is the admittance of lines connecting neighboring buses. In case of a line outage, this information is outdated and so it is important to investigate the behavior of the index in this scenario. We identified that the proposed index accurately tracks the system stress in real-time even when using the outdated line admittances. In order to study the effect of topology change on the VSI, the line between buses $15$ and $23$ in the IEEE $30$-bus system is taken out of service at time $t=138$ secs as the load is increasing. Fig. 11 plots the VSI at buses $14,15,18$ & $19$ and it is observed that VSI drops the moment the line outage occurs at time $t=138$ secs, indicating that the system is stressed and margin to voltage collapse point has reduced. It is also observed from Fig. 11 that the values of the VSI at various buses are very similar when using un-updated and updated connectivity information. Thus the distributed VSI calculation methodology provides a reasonable estimate of the VSI for the operator/relay to trigger controls/alarms in the time it takes to correct the admittance information, thus improving the situational awareness of the grid. In addition, it can also be observed that the index at bus $15$ that is closest to the line outage reduces the most and this fact can be used to identify the fault location in an unsupervised manner.

Another disturbance that is common is the outage of a generator. We investigated this scenario and observed that the proposed VSI accurately tracks the system stress after a generator outage and can identify the critical buses. Consider the IEEE-$30$ bus system, generator two located at bus 2 is outaged at 50 secs and generator five located at bus 23 is outaged at 93 secs. We can observe from Fig. 12 that at time t= 50 secs and t=93 secs, VSI at all the buses drop, indicating rising overall system stress. Further, the drop is more at locations closer to the generator with larger power output (gen-2), quantifying the severity of the outage.

Thus, these results validate the proposed VSI’s ability to identify LTVI and demonstrate it’s utility to monitor voltage stability using measurements in a distributed manner under various disturbances such as network reconfiguration and generator outage.

V Conclusion

This paper proposes a PMU measurement-based voltage stability index that can accurately identify long term voltage instability of the system. The key novelty of this work is the mathematical derivation of the index to reflect voltage security and its distributed nature. The index is derived by analyzing the power flow equations in the rectangular form as circles, and at the nose point, the value of the index at the critical bus is zero. In addition, the distributed communication framework between neighboring buses makes the calculation of index scalable and secure for the grid. In the various test scenarios such as noisy measurements, different load increase directions, line outage and generator outage, the index detected the critical bus and quantified the stress in the system. The proposed index behavior is also compared with distributed and decentralized methods and is shown to have a superior & reliable performance, particularly for noisy measurements and large grids. The wide area nature of the proposed index makes it robust to measurement and system noise that adversely affects similar techniques such as local Thevenin methods, etc. Furthermore, the proposed index can localize the line/gen outage locations in an unsupervised manner purely using PMU measurements, making it a promising method for event detection and localization of other events. Finally, the distributed nature of the index makes it possible to utilize cloud computing infrastructure and other recent trends in the big data analytics field for efficient computation and storage.

-A Matrix Representation of Standard Circle:

First, the matrix form of a standard circle is derived. Let $\mathit{z}=\mathit{x}+j\cdot\mathit{y}$ in complex plane be set of all the points in a circle $\mathit{{C_{k}}\left(z,z^{*}\right)}$ with radius $\rho$ and center $\gamma=\alpha+j\cdot\beta$. Then, the equation of this circle is given by (9)

[TABLE]

where $\mathit{A}=1$, $\mathit{B}=-\gamma^{*}$, $\mathit{C}=-\gamma$ and $\mathit{D}=\gamma\cdot\gamma^{*}-\rho^{2}$.

-B Proof for linear behavior of proposed index in case of transmission networks:

Fig. 13 shows a simple $2$-bus system which consists of one load fed by an infinite bus through a transmission line. The infinite bus is represented by an ideal voltage source with constant voltage given by $\overline{E}=E_{r}+j\cdot E_{i}=1+j\cdot 0$. The transmission line is represented by its admittance $\overline{Y}=g+j\cdot b$. The load consumes an apparent power of $\overline{S}=P+j\cdot Q$ with a constant power factor $cos(\phi)$ and it is represented by a constant equivalent impedance $\overline{Z}_{L}=R_{L}+j\cdot X_{L}$.

The proposed distributed non-iterative index is given by

[TABLE]

where $\Delta^{*}$ is given by (12), $\Delta^{*}_{\textrm{no-load}}$ is the value of $\Delta^{*}$ with zero load and voltage of 1 p.u. at all buses of the network.

[TABLE]

where

[TABLE]

Let load power factor be unity ($cos(\phi)=1$). Now, calculate the terms $\mathbf{b}_{p}$, $\mathbf{b}_{q}$, $c_{p}$ and $c_{q}$ as shown below for the $2$-bus example that is described above.

[TABLE]

Substitute above terms ($\mathbf{b}_{p},\mathbf{b}_{q}$,$c_{p},$ and $c_{q}$) in (11) and upon further simplification, (11) is given by

[TABLE]

Bounds of proposed index: First, we will show that the proposed index is bounded between $1$ and [math] at no load ($P=0$) and nose point of the PV curve ($P=P_{max}$) respectively. The maximum transferable power ($P=P_{max}$) for the $2$-bus system from Fig. 13 when the load power factor is unity is given by

[TABLE]

[TABLE]

Linear behavior of proposed index for transmission networks: Second, we show that even though the proposed index (13) is a quadratic equation in real power ($P$), the equation behaves linearly for various practical power system operation ranges since the coefficient of ($P^{2}$) term is negligible. Specifically, we showcase that the proposed index behaves linearly under various $r/x$ ratios of the transmission line and various load power factors.

Effect of $x/r$ ratio on the linear behavior of the index: The $x/r$ ratio is a critical power system parameter that effects the maximum power transfer capability of the grid. For example, a low $x/r$ ratio enables more maximum power transfer when compared to that of a high $r/x$ ratio power line. Fig. 14 shows that the proposed index (13) behaves linearly for various values of $x/r$ ratios corresponding to different levels of voltage ($22$ kV to $765$ kV) [47, 48]. Fig. 14 shows that the proposed index behaves linearly for high $x/r$ values (transmission networks) and behaves in a slightly curved fashion for low $x/r$ values (distribution networks).

For example, when $x/r=5$ and $x/r=10$ then (13) can be further simplified to (14) and (15) respectively. From (14) and (15), we can see that the high values of $x/r$ makes the coefficient of “$P^{2}$” term negligible and hence the linear part of the equation i.e., “$P$” term dominates resulting in the linear behavior of the proposed index.

[TABLE]

[TABLE]

Additionally, the proposed index does not deviate a lot from its original slope and hence it tends to stay fairly linear. For example, we know that the first order differentiation of (13) gives the slope of the proposed index. The second order differentiation of the (13) gives the information about the change of the slope of the proposed index. From (16), when $x/r=10$, this rate of the change of slope of the proposed index is bounded by the coefficient of “$P^{2}$” and it is negligible ($-7.8424e^{-4}\approx 0$) for transmission network lines due to high $x/r$ ratio as discussed. In other words, as the $x/r$ ratio increases, the coefficient of $P^{2}$ becomes negligibly small. Hence the proposed index tends to behave linearly for all the transmission network case studies presented in this paper.

[TABLE] 2. 2.

Effect of load power factor on the proposed index: For a given $x/r$ ratio of $10$, Fig. 15 shows the proposed index (13) for various values of load power factor. The proposed index’s linear behavior is not effected unless the power factor decreases to as low as $0.6$ which is not realistic for transmission networks.