Estimating the saddle-node bifurcation point of static power systems using the holomorphic embedding method

Electrical Power and Energy Systems 84 (2017) 1–12

Contents lists available at ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

Estimating the saddle-node bifurcation point of static power systems using the holomorphic embedding method Shruti D. Rao a,⇑, Daniel J. Tylavsky a, Yang Feng b a b

School of Electrical Computer and Energy Engineering, Arizona State University, Tempe, AZ 85287, USA Siemens Industry, Inc., RC-US EM SG PTI SW ARC, 4920 Westway Park Boulevard, Suite 150, Houston, TX 77041, USA

a r t i c l e

i n f o

Article history: Received 24 November 2015 Received in revised form 29 March 2016 Accepted 13 April 2016

Keywords: Analytic continuation Saddle-node bifurcation point Voltage stability Holomorphic embedding Holomorphically embedded power flow method Voltage collapse point

a b s t r a c t Voltage stability studies have been progressively gaining importance in the power engineering community. Predicting the saddle-node bifurcation point (SNBP) of a power system has become more critical as the power-system loading has increased in many places without a concomitant increase in transmission resources. Since a Newton–Raphson power-flow method is inherently unstable near the SNBP, adaptations such as continuation methods have been used as stabilizers. A new class of nonlinear equation solvers known as the holomorphic embedding method (HEM) is theoretically guaranteed to find the high-voltage solution to the power-flow problem, if one exists, up to the SNBP, provided sufficient precision is used and the conditions of Stahl’s theorem are satisfied by the equation set. In this paper, four different HEM-based methods to estimate the saddle-node bifurcation point of a power system, are proposed and compared in terms of accuracy as well as computational efficiency. Ó 2016 Elsevier Ltd. All rights reserved.

Introduction Because of the difficulty of siting transmission lines, utilities are often forced to serve increased electric power demand, without a concomitant expansion of infrastructure. This can lead to the system being operated closer to its saddle-node bifurcation point (SNBP) and therefore closer to voltage collapse than desired. There are many examples of black-outs occurring because of a slow reduction in voltage magnitudes at buses over a time scale of a few minutes to hours followed by a sudden sharp fall in the voltage magnitudes, e.g., [1]. One goal of a voltage stability study is to determine the voltage stability margin, e.g., the amount of real and/or reactive power that can be added before the system experiences voltage collapse, with the distance to the SNBP being a quick indicator of stability margin. Significant work has been done to analyze the voltage stability of systems. It has been shown that no dynamics are required to be modeled in order to obtain the SNBP of a system and that the small signal voltage stability limit depends only on the steady-state characteristics of the system [2–4]. New techniques of analyzing the voltage stability in steady-state systems while considering system limits such as generator reactive power limits, voltage magnitude ⇑ Corresponding author. E-mail addresses: [email protected] (S.D. Rao), [email protected] (D.J. Tylavsky), [email protected] (Y. Feng). http://dx.doi.org/10.1016/j.ijepes.2016.04.045 0142-0615/Ó 2016 Elsevier Ltd. All rights reserved.

constraints, and AVR constraints have been presented in numerous papers, ([5–13] to cite only a few.) Optimization algorithms such as the genetic algorithm, particle swarm optimization and their variants have been used to detect the closest SNBP in [14–17]. The analysis of voltage stability of power systems and computation of maximum loadability while incorporating dynamic constraints along with steady-state constraints has been explored in [18–24]. Wide-area-measurement-based voltage stability analysis using modified coupled single-port models has been examined in [25]; while physical constraints such as var limits are not considered in that work, they are shown to be important. Ultimately, the power-flow (PF) equations and their solution is at the core of voltage stability analyses and the aforementioned voltage stability analysis methods use iterative solution methods which suffer from non-convergence issues. Two classes of iterative methods for solving the PF problem are well known, Gauss–Seidel and Newton Raphson (NR) (notable is the Fast Decoupled Load Flow) [26–30], both of which have many variants developed to improve the convergence properties ([31–34] to cite only a few.) The problems with these classes of solvers are also well known: initial estimate dependence [27,28]; non-convergence when a solution exists and convergence to the wrong (low-voltage/large angle) solution [35]. To date, no variants within these classes have been shown to deal with these problems in a consistent way. Often, the convergence issues are exacerbated near the bifurcation point. The most popular class of industrial methods, the NR class,

2

S.D. Rao et al. / Electrical Power and Energy Systems 84 (2017) 1–12

converges reliably for systems with nominal voltage profiles [36], but often struggles near the SNBP. The method used in most PF applications to estimate the SNBP is the continuation power flow (CPF), which is an NR-based method [37]. The computational complexity of the CPF is much higher than a simple NR PF since it requires solving a new PF problem at each step as one moves along the P–V curve toward the SNBP [37]. The problems of a singular Jacobian matrix occurring at the SNBP are eliminated in [38] by reformulating the powerflow problem and introducing a new bus type called AQ bus which has the voltage angle and reactive power consumption specified. The holomorphic embedding method (HEM) is a more recently developed class of nonlinear equation solvers that is guaranteed to converge to the operable solution for the basic PF problem, if it exists, provided the conditions of Stahl’s theorem are satisfied [39,40,42]. It is recursive, rather than iterative, and can be configured so that the initial state needed to guarantee the above properties can be easily calculated. The advantages of the HEM can be exploited to develop methods that can reliably estimate the SNBP of a system. The HEM uses the concept of holomorphic embedding (HE) to convert the non-holomorphic power-balance equations (PBE’s) of the PF problem into a set of holomorphic functions. A holomorphic function is a complex-valued analytic function, which has the property that it is infinitely complex differentiable around every point within its domain. One property of holomorphic functions important for the purpose here is that they can be represented by their Taylor series in a neighborhood of each point in their domain [47]. Using the HE formulation, the voltages at all buses and the reactive power at the PV buses are expressed as Maclaurin series of a complex embedding parameter, a. With the properties of holomorphic functions and the use of analytic continuation, Stahl’s theorem [40] guarantees that, if an operable solution exists at the given loading level, the correct voltage solution will be obtained using Padé approximants of the holomorphic series as long as the correct germ is used [39,43]. (The germ for the HEM is analogous to the initial estimate of the solution in the NR method and will be explained in later sections. However, unlike the NR initial estimate, the germ must be systematically obtained by solving a set of equations.) The methods in the literature that calculate the SNBP rely on solving successive problems, each of which has the complexity of the PF problem. This paper presents four HEM-based methods that estimate the SNBP of a system, three of which do not require multiple PF problems to be solved. All of these methods rely on an important property of a Padé approximant, that it is the maximal analytic continuation of the given function [40]. The paper is organized as follows: Discussed in Section ‘‘Two HEM formulations” is a formulation that allows load and realpower-generation extrapolation (a property critical to the three approaches proposed here) for the power-flow problem. Also discussed in this section is a formulation previously published and why it cannot be used for load extrapolation. Section ‘‘Padé approximants and branch cuts” contains the method of using the roots of Padé approximants to find the SNBP of the system and the fundamental theory behind this approach. In Section ‘‘The sigma methods”, two so-called ‘‘sigma methods” used to identify weak nodes of the power system and to estimate the SNBP of the system are introduced [44]. In Section ‘‘Numerical results for scaling all loads uniformly”, the results of numerical experiments that compare the SNBP’s predicted by different methods when all loads are scaled uniformly are presented. Section ‘‘Formulation to allow loads at different buses to be scaled by different amounts” contains the development of a formulation that can be used to scale the load vector in an arbitrary direction and, consequently, be used

to obtain the SNBP. Section ‘‘Numerical results with loads at different buses scaled by different amounts” provides the numerical results for the formulation described in Section ‘‘Formulation to allow loads at different buses to be scaled by different amounts”. Section ‘‘Incorporating var limits in the SNBP estimation” presents the results when var limits are accounted for and Section ‘‘Proposed ZIP load model for the HEM” provides the formulation for polynomial ZIP load models. Finally, conclusions are presented in Section ‘‘Conclusion”. Two HEM formulations To apply the HEM to a complex-valued problem requires that the system of equations to be solved be holomorphic. Because of the presence of the complex conjugate operator, the traditional PF equations are not holomorphic. Hence the first step in developing a proper HEM formulation is to render the PBE’s holomorphic. Consider a general ðNÞ-bus system consisting of a slack bus, called slack, a set m consisting of PQ buses, a set p consisting of PV buses which are not on var limits and a set q consisting of PV buses on maximum/minimum var limits. The PBE for a PQ bus with a constant power load is given by N X S Y ik V k ¼ i ; Vi k¼1

i2m

ð1Þ

where, Y ik is the (i; kÞ element of the bus admittance matrix, and Si , and V i are the complex power injection and voltage at bus i, respectively. (The HEM model for polynomial ZIP loads will be discussed in Section ‘‘Proposed ZIP load model for the HEM”.) The traditional defining equations for a PV bus are given by (2) and (3).

! N X   Pi ¼ Re V i Y ik V k ;

8i 2 p

ð2Þ

k¼1

jV i j ¼ V sp i ;

8i 2 p

ð3Þ

where Pi denotes the real power injection and V sp i is the specified voltage magnitude at bus i. PV buses on var limits are treated similar to PQ buses with their reactive power generation fixed at the appropriate limit and the real power generation given by (2). Formulation that allows extrapolation of the load The above non-holomorphic equations can be holomorphically embedded in an infinite number of ways. It is possible to embed (1)–(3) in such a way that the solution obtained at different values of real a, represents the solution (if it exists) when the complex power injections at the load buses and real power at generation buses are scaled by a factor of a. It is necessary to have such a formulation in order to be able to estimate the SNBP of the system without having to solve a new PF problem at different loading conditions. The HEM formulations published in the past, do not allow one to scale the load by a factor of a, since they solve the given powerflow problem at a = 1.0, because they are consistent with the power system equations only at a = 1.0. This will be explained in more detail in Section ‘‘Formulation with a simple germ”. Consider the following set of holomorphically embedded equations, where (4) represents the PBE for PQ buses, (5) represents the voltage magnitude constraint for the slack bus, (6) represents the PBE for the PV buses, (7) represents the voltage magnitude constraint for the PV buses and (8) represents the PBE for PV buses on var limits. N X aS Y ik V k ðaÞ ¼  i  ; V i ða Þ k¼1

i2m

ð4Þ

3

S.D. Rao et al. / Electrical Power and Energy Systems 84 (2017) 1–12

V i ðaÞ ¼ V sp i ;

i 2 slack

ð5Þ

N X aðPgi  Pli Þ  jðQ gi ðaÞ  aQ li Þ Y ik V k ðaÞ ¼ ; V i ða Þ k¼1

 2  ; V i ðaÞ  V i ða Þ ¼ V sp i

k¼1

i2p

i2p

N X aðPgi  Pli Þ  jðQ gi Y ik V k ðaÞ ¼ V i ða Þ k¼1

 aQ li Þ

;

i2q

ð8Þ

where P gi denotes the real-power generation, Pli denotes the realpower load, Q gi (a) denotes the reactive-power generation, Q li denotes the reactive-power load at bus i and Q gi lt denotes the respective maximum/minimum var limit for a generator bus on maximum/minimum var limits. The real- and reactive-power injections at bus i are given by:

Pi ¼ Pgi  Pli ;

8i 2 p

Q i ¼ Q gi ðaÞ  Q li ;

ð9Þ

8i 2 p

ð10Þ

Using (9), (6) and (8) can be written as: N X

Y ik V k ðaÞ ¼

k¼1

aPi þ jaQ li  jQ gi ðaÞ ; i2p V i ða Þ

N X aPi þ jaQ li  jQ gi lt Y ik V k ðaÞ ¼ ; V i ða Þ k¼1

ð11Þ





 V i ½0 þ V i ½1a þ V i ½2a2 þ       2  ;  V i ½0 þ V i ½1a þ V i ½2a2 þ    ¼ V sp i

i2p

ð17Þ

ð18Þ

N X   Y ik V k ½0 þ V k ½1a þ V k ½2a2 þ    k¼1

¼

aP i þ jaQ li  jQ gi lt ; V i ½0 þ V i ½1a þ V i ½2a2 þ   

i2q

ð19Þ

Since it is required that (15)–(19) hold true for any value of a, the coefficients of the respective powers of a on both sides of the equations must be equal. In order to find the series coefficients that satisfy these equations, the inverse of the voltage power series on the RHS has to be represented as a power series. To achieve this, let the inverse of the voltage function V(a), be represented by another power series, W(a), defined by W(a) = W[0] + W[1]a + W[2]a2 +    where the relationship between W(a) and V(a) is given by (20).

WðaÞ ¼

1 VðaÞ

ð20Þ

Thus Eqs. (15)–(19) are now modified to:

i2q

ð12Þ

Since V(a) and Q gi (a) are holomorphic functions of the parameter a, they can be expressed as Maclaurin series given by:

VðaÞ ¼ V½0 þ V½1a þ    þ V½nðaÞ

n

Q gi ðaÞ ¼ Q gi ½0 þ Q gi ½1a þ    þ Q gi ½nðaÞn

ð13Þ

where Q gi (a) is a real-valued series. Note that the variable V  in (4), (7), (11) and (12) is embedded with a instead of a in order to ensure that the Cauchy-Riemann conditions are satisfied, thereby retaining the holomorphicity of the function. For this reason, it is important to keep in mind, that (4), (7), (11) and (12) represent the given power system only for real values of a [45]. The Maclaurin series for V  ða Þ is given by:

V  ða Þ ¼ V  ½0 þ V  ½1a þ    þ V  ½nðaÞn

ð14Þ

Eqs. (4), (5), (7), (11) and (12), represent a new formulation that allows the loads at all buses, and the real-power generation at PV buses, to be scaled by a factor of a. Note that in (11), only the reactive-power generation Q gi is written as a function of a instead of the net reactive-power injection Q i . This allows the reactive-power load to be scaled by a factor of a, while the variable value of reactive-power generation needed to maintain bus voltage control, is calculated from the power series. Also in (12), Q gi lt is not multiplied by a, since this is fixed for a bus on var limits and cannot be scaled with a. This is the only difference between the equation for PQ buses and that for PV buses on var limits. By substituting the series V(a), V  (a Þ and Q(a) given by (13) and (14) into (4), (5), (7) (11) and (12), we obtain:

k¼1

  ¼ aSi W i ½0 þ W i ½1a þ W i ½2a2 þ    ;





V i ½0 þ V i ½1a þ V i ½2a2 þ    ¼

aSi ; i 2 m ¼   V i ½0 þ V i ½1a þ V i ½2a2 þ    i 2 slack

ð15Þ ð16Þ

V sp i ;

i2m

i 2 slack

ð21Þ ð22Þ

N X   Y ik V k ½0 þ V k ½1a þ V k ½2a2 þ    k¼1

  ¼ ðaPi þ jaQ li  j Q gi ½0 þ Q gi ½1a þ Q gi ½2a2 þ    Þ     W i ½0 þ W i ½1a þ W i ½2a2 þ    ; i 2 p



ð23Þ



V i ½0 þ V i ½1a þ V i ½2a2 þ       2  ;  V i ½0 þ V i ½1a þ V i ½2a2 þ    ¼ V sp i

i2p

ð24Þ

N X   Y ik V k ½0 þ V k ½1a þ V k ½2a2 þ    ¼ ðaPi þ jaQ li  jQ gi kt Þ k¼1

   W i ½0 þ W i ½1a þ W i ½2a2 þ    ;

i2q

ð25Þ

In order to guarantee that the HEM can find an operable solution (if one exists), it is critical to have the correct germ [39,40,43]. The germ is the solution obtained by solving the holomorphically embedded equations representing the PF problem at a = 0. This extrapolation-capable formulation presented above however, poses some challenges when calculating the germ. The system of equations to be solved for the germ is given by (26a)–(26d). N X Y ik V k ½0 ¼ 0;

i2m

ðaÞ

k¼1

k¼1

k¼1

 V i ½0 þ V i ½1a þ V i ½2a2 þ    ¼ V sp i ;

N X   Y ik V k ½0 þ V k ½1a þ V k ½2a2 þ . . .

N X jQ ½0 Y ik V k ½0 ¼ V  gi½0 ;

N X   Y ik V k ½0 þ V k ½1a þ V k ½2a2 þ   





aPi þ jaQ li  j Q gi ½0 þ Q gi ½1a þ Q gi ½2a2 þ       ; i2p ¼ V i ½0 þ V i ½1a þ V i ½2a2 þ   

ð6Þ ð7Þ

lt

N X   Y ik V k ½0 þ V k ½1a þ V k ½2a2 þ   

i

 2  ; V i ½0  V i ½0 ¼ V sp i

i 2 p ðbÞ i2p

ðcÞ

V i ½0 ¼ V sp i 2 slack i ; N X jQ gi lt Y ik V k ½0 ¼ V  ½0 ; i2q

ðdÞ

k¼1

i

ðeÞ

ð26Þ

4

S.D. Rao et al. / Electrical Power and Energy Systems 84 (2017) 1–12

Note that (26b) and (26c), representing the PBE and the voltage magnitude constraints for PV buses and (26e) representing the PBE for PV buses on var limits, cause the system of equations that need to be solved to obtain the germ to be nonlinear. Since the NR is known to have good convergence properties for lightly loaded systems and the germ represents the solution for a no-load, no-realpower-generation scenario, this germ can be obtained by using the NR for this special case. However, it is possible, though unlikely, that NR might not give the right germ; it could converge to the low-voltage solution for this lightly-loaded problem (the only loads are the shunts). However, if one wants an algorithm that is guaranteed to find the correct germ, the HEM can be applied to solve (26) as explained in Appendix A. Once the germ is obtained, the recursive relations to find the remaining terms of the power series, are obtained by equating the coefficients of same powers of a on both sides of Eqs. (4), (5), (7), (11) and (12). Eqs. (27) and (28) provide the recurrence relations for the PQ buses and the slack bus respectively. N X Y ik V k ½n ¼ Si W i ½n  1;

i2m

ð27Þ

k¼1

V i ½n ¼ 0; n > 0 i 2 slack

ð28Þ

The recurrence relation for the PBE of PV buses is obtained as: N n X X Y ik V k ½n ¼ ðP i þ jQ li ÞW i ½n  1  j Q gi ½kW i ½n  k k¼1

!

ð29Þ

k¼0

Eq. (29) can be further simplified to: N X Y ik V k ½n ¼ ðP i þ jQ li ÞW i ½n  1 k¼1

j

n1 X

! Q gi ½kW i ½n  k þ Q gi ½nW i ½0 þ Q gi ½0W i ½n

k¼1

ð30Þ The unknowns Q gi [n] and W i [n] can be moved to the LHS, leaving all known quantities on the RHS to obtain: N X Y ik V k ½n þ jQ gi ½nW i ½0 þ jQ gi ½0W i ½n k¼1

¼ ðPi þ jQ li ÞW i ½n  1  j

n1 X

! Q gi ½kW i ½n  k

ð31Þ

k¼1

Notice that in the above equation, a new unknown, W[n], is introduced. In order to calculate the coefficients of the inverse power series, W[n], both sides of (20) are multiplied by V(a).

WðaÞVðaÞ ¼ 1    W½0 þ W½1a þ W½2a2 þ    V½0 þ V½1a þ V½2a2 þ    ¼ 1 ð32Þ By equating the coefficients of the same powers of a on both sides of (32) a relation between W[n] and V[n] is obtained as given in (33). 1 W½0 ¼ V½0

W½1V½0 þ W½0V½1 ¼ 0 .. .

ð33Þ n1 X

W½0V½n þ V½0W½n ¼ 

W½kV½n  k;

nP1

k¼1

Equating the coefficients of same powers of a, on both sides of the embedded voltage magnitude constraint equation given by (18), we obtain

V i ½0V i ½n þ V i ½nV i ½0 ¼ ðV i ½1V i ½n  1 þ    þ V i ½n  1V i ½1Þ n1 X ¼  V i ½kV i ½n  k;

nP1

ð34Þ

k¼1

Note that all the quantities on the RHS of (34) are known. The recurrence relation for PV buses on var limits is obtained as: N X Y ik V k ½n ¼ ðPi þ jQ li ÞW i ½n  1  jQ gi

lt

 W i ½n

k¼1

)

N X Y ik V k ½n þ jQ gi

lt

 W i ½n ¼ ðPi þ jQ li ÞW i ½n  1

ð35Þ

k¼1

Note that in (35) also, there is an extra unknown W[n], and thus an additional equation given by (33) needs to be added to the system of equations for PV buses on var limits in order to ensure that the problem is not under-determined. The system of linear Eqs. (27), (28), (31), and (33)–(35) represent the recurrence relations for the formulation given by (4), (5), (7), (11) and (12). It is important to keep in mind that Q gi (a) is a purely real series. In order to obtain the power series coefficients using a single linear matrix equation, the unknowns need to be split into real and imaginary parts. Since the power series resulting from solving these recurrence relationships will only converge within its radius of convergence, Padé approximants are used to evaluate these series everywhere in the function’s domain, including where the series is nonconvergent. A Padé approximant is a rational-function (ratio of two polynomials) approximation to a power series that can yield a more accurate approximation to the defining function than the truncated series of the same length. A common Padé approximant notation to indicate that the degree of the numerator polynomial is L and the degree of the denominator polynomial is M, is ½L=Mf ðaÞ , given in (36).

f ðaÞ  ½L=Mf ðaÞ ¼

a0 þ a1 a þ    aL aL b0 þ b1 a þ    bM aM

ð36Þ

The direct or matrix method used in this work for finding the rational approximant [49] is straight forward and involves solving a dense set of linear equations of dimension M and then a forward substitution through a dense lower triangular matrix of dimension L + 1. Since the Padé approximant is the maximal analytic continuation of the given power series [40], its value for any a within the function’s domain is guaranteed to be the HV solution for the above embedding, provided the conditions of Stahl’s theorem are satisfied [40,41]. Stahl’s theorem states that if the function f which is analytic at infinity and has all its singularities in a compact set of logarithmic capacity 0; then the convergence domain where the sequence {[mj/nj]} (j2N) of Padé approximants converges in capacity to f for any sequence {[mj/nj]} satisfying the condition given by (37), is identical with the extremal domain of the function f up to a set of capacity zero.

mj þ nj ! 1;

mj ! 1 as j ! 1 nj

ð37Þ

Note that once the Padé approximants for the voltages are obtained, the SNBP may be found by varying a in the Padé approximants (not re-solving the PF problem), using a binary search approach, until the boundary is reached at which the Padé approximant no longer obeys the PBE’s. This point is the boundary of the function’s domain. Thus this method may be viewed as a curvetracing of the P–V curve for a system. It is important to point out that it is not necessary to evaluate the Padé approximants for all buses in order to check for convergence, but need be done only for a few arbitrarily chosen buses.

S.D. Rao et al. / Electrical Power and Energy Systems 84 (2017) 1–12

The Padé approximants will have converged when the solution is reached; however, if there is no solution to the given problem, the Padé approximants will oscillate. If the given loading on the system is too low, the solution at a = 1.0 may converge with very few terms in the power series. However, care must be taken to ensure that a sufficient number of terms is included so that an accurate solution is obtained at higher values of a, especially near the SNBP. Formulation with a simple germ While the germ calculation for the formulation presented in S ection ‘‘Formulation that allows extrapolation of the load” requires solving a nonlinear problem, this section will present a nonextrapolating formulation whose germ may be found by inspection. One of the ways to holomorphically embed (1)–(3) is given by (38)–(41) where Y ik trans corresponds to the ‘‘series branch” part of the admittance matrix and Y i shunt corresponds to the shunt part of the admittance matrix [42,45]. Since this formulation is nonextrapolating and is valid at only a = 1.0, (38) can be used to represent PV buses on var limits as well, by setting the reactive power generation to be at its appropriate limit. N X a S Y ik trans V k ðaÞ ¼  i   aY i V i ða Þ k¼1

V i ðaÞ ¼ 1 þ ðV sp i  1Þa;

shunt V i ð

aÞ; i 2 m

i 2 slack

N X aPi  jQ i ðaÞ Y iktrans V k ðaÞ ¼  aY i V i ða Þ k¼1

   2  1 ; V i ðaÞ  V i ða Þ ¼ 1 þ a V sp i

ð38Þ ð39Þ

shunt V i ð

i2p

aÞ; i 2 p

ð40Þ

ð41Þ

The germ for this method can be obtained easily and uses the same process as given in Appendix A. Once the germ is obtained, the approach to obtaining the solution for this formulation is very similar to that presented in Section ‘‘Formulation that allows extrapolation of the load”. It should be mentioned that this formulation not only provides an obvious germ, but also is in general a bit more advantageous in terms of numerical stability and precision. This is because by doing so, the magnitude of the voltages does not vary as much when going from a = 0 to a = 1.0, and thus the Pade approximants behave slightly better numerically in general [45]. A drawback of the above formulation is that it represents the given power system only at a = 1.0 and has no meaningful interpretation at any other value of a. Since Y i shunt is also multiplied by a, at any value of a – 1, the solution corresponds to a network whose shunt-branch parameters differ from those of the original network. Also, the voltage magnitude constraints, embedded as given in (39) and (41) will hold true only at a = 1.0. For any other value of a, the voltage magnitude constraints for the PV buses and the slack bus will not be obeyed. However, because the HEM, with this formulation, is theoretically guaranteed to converge to the HV solution (precision issues notwithstanding) this formulation can be used to calculate the SNBP using a binary search approach. This involves solving multiple PF problems and is of the order of complexity of the CPF, where complexity here is taken as the number of PF’s that must be solved. The computational complexity of the method proposed in Section ‘‘Formulation that allows extrapolation of the load” is much less than that of the repeated solution approach proposed here, since for that method, only the Padé approximants for the voltages must be evaluated until the function’s domain boundary is reached, which is indicated by high bus-power mismatches.

5

Padé approximants and branch cuts While the power series obtained may or may not converge depending on the value of a, i.e., for various loading conditions, there are cases where the power series diverges but the solution for the PF problem exists. The power series represents one of the branches of the solution to the embedded problem. The radius of convergence of the power series is determined by the closest singularity of the function to the origin [47]. Since the full solution V(a) is a multi-valued function, more precisely an algebraic curve, it has no poles; therefore the only singularities are its branch points on the s-plane (i.e. points where one or more branches coincide) [48]. The first branch point encountered along the real axis is precisely the bifurcation point that we seek [45]. However, in general, other branch points on the complex a-plane may be closer to a = 0 [50]. In that case, the power series has a reduced convergence radius, and that is why Padé approximants are needed. If the point of interest (given loading condition) is outside the radius of convergence, then the voltage series will diverge. The series therefore needs to be ‘redefined’ such that a converged solution is obtained for all loading levels where a solution exists. The sequence of near-diagonal Padé approximants, proven to be the maximal analytic continuation of the power series by Stahl’s theorem [40], can be used to extend the convergence region of the power series over the domain of the function. There are several ways to calculate the Padé approximant and one commonly used is the direct/matrix method, which is well-described in [49]. A common Padé approximant notation is to use L as the degree of the numerator polynomial and M as the degree of the denominator polynomial. Though the matrix method described in [49] allows the calculation of a rational approximant of an arbitrary degree, Stahl’s theory [51] indicates that the close-to-diagonal Padé approximants (limM!1 L=M ¼ 1Þ yield the best accuracy when evaluating the power series outside its radius of convergence. Stahl’s work [40,41], originated from Nuttall’s work [48], which proved that the poles of the close-to-diagonal Padé approximant accumulate on Stahl’s compact set. When a scaling formulation is used, the distance, from the origin to the closest point on the compact set obtained from the voltage power series gives the load-scaling value at the SNBP. The zeros of the Padé approximant simulate Stahl’s compact set in a manner similar to that of the poles; therefore the smallest real zero or the smallest real pole of the Padé approximant can be used to determine the proximity of the system to its SNBP along with the poles. Some interesting results on the complexity of Stahl’s compact set have been shown in [50]. While the SNBP can only be precisely obtained if computing the power series with an infinite number of terms and infinite precision, promising numerical results, given in Sections ‘‘Numerical results for scaling all loads uniformly”, ‘‘Numerical results with loads at different buses scaled by different amounts” and ‘‘Incorporating var limits in the SNBP estimation”, show that the method works well for predicting the SNBP with a reasonable number of terms and finite precision.

The sigma methods The method described in this section to calculate the SNBP of a system is based on calculating a set of complex-valued indices [44]. The idea behind the method is to develop a two-bus equivalent system for each node of the power system consisting of that node along with the slack bus, an equivalent that preserves only the slack bus voltage and the voltage at the retained bus. It is a local nonlinear Thevenin-type equivalent at each node. This section will describe the derivation of the sigma indices, starting with a twobus example.

6

S.D. Rao et al. / Electrical Power and Energy Systems 84 (2017) 1–12

Consider a two-bus system with a slack bus and a PQ bus as shown in Fig. 1, where Z is the line impedance, S is the complexpower injection at the PQ bus (and hence P will be negative for real-power loads), V 0 is the slack bus voltage and V is the PQ bus voltage. It is easy to show that

U ¼1þ

r

ð42Þ

U

where U ¼ V/V 0 and r is defined as

r¼

ZS

ð43Þ

jV 0 j2

The roots of (42), which is a quadratic equation, are:

U¼

1  2

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ rR  r2I þ jrI 4

ð44Þ

where rI and rR are the imaginary and real parts of r, respectively. If the slack bus voltage is assumed to be 1.0 p.u., the two roots represent the high- and low-voltage solutions for the given twobus network. The two solutions meet at the point at which the radicand becomes zero, the SNBP. For the high-voltage solution to exist, it is necessary that the radicand in (44) be positive. Thus the condition to ensure that the system is short of or at its static voltage collapse point, called the ‘r condition’ is given by:

1 þ rR  r2I P 0 4

ð45Þ

Thus the index r can be used to estimate the SNBP of the system. Sigma method—for weak bus determination The method, given in [44], to estimate the voltage stability margin for a system, is an extension of the two-bus derivation to a general (N þ 1Þ-bus system. Eq. (42) can be extended to represent each node of a general (N þ 1Þ-bus system consisting of PV and PQ buses as follows,

ari ðaÞ V i ðaÞ ¼ 1 þ   V i ða Þ

ri ½n  1 ¼

V i ðaÞV i ða Þ ¼ V i ða Þ þ ari ðaÞ Solving this for the power series

ð47Þ

ri (a) yields:

ari ðaÞ ¼ V i ðaÞV i ða Þ  V i ða Þ

ð48Þ

where ri (a)=ri [0]+ri [1]a þ ri [2]a2 . . . . By equating the constant terms on both sides of (48), we obtain,

V i ½0V i ½0 ¼ V i ½0

V i ½0ðV i ½0  1Þ ¼ 0

ð49Þ

Thus the germ corresponding to the high-voltage solution, must be V i [0] = 1.0. Equating higher powers of a on both sides of (48), the r terms can be obtained as:

Fig. 1. Two-bus system diagram.

n>1

ð50Þ

k¼0

Note that (49) implies that (46), in its given form, is valid only if the voltage germ for the PF problem is 1.0 p.u. for the high voltage solution, which is not generally true for the extrapolation formulation given in Section ‘‘Formulation that allows extrapolation of the load”. Such a germ is obtained for the formulation given by (38)–(41); however this formulation cannot be used for extrapolation, which is problematic for the purposes here since the objective when using (46) is to increase a (extrapolate) until the r condition of (45) is met. A modification to (46), which would allow extrapolation of a, is given below (Section ‘‘Extrapolating sigma method–for SNBP estimation”) Regardless of its extrapolation potential, it has been observed that the margin by which r satisfies the r-condition is a reliable indicator of which buses are the weak nodes of the system, with weak nodes closer to violating the condition [44]. Extrapolating sigma method–for SNBP estimation If (46) is viewed as the mathematical analog of a two-bus equivalent of the original system reduced to a slack bus and bus of interest, then such an approximate equivalent might be used to predict the SNBP by scaling a until the sigma condition is violated [44]; however (46), in its present embedding, is valid for extrapolation to the SNBP only for voltage series obtained from HE formulations that allow a extrapolation and where the germ for the voltages is 1.0. For the formulation presented in Section ‘‘Formulation that allows extrapolation of the load”, which allows extrapolation of a, the voltage germs are not necessarily 1.0 p.u. for practical systems. Hence (46) is only strictly valid for a = 1.0 and extrapolation is approximate. We can modify (46) for use with the proposed extrapolation formulation, where the voltage germ is not 1.0. This can be done by eliminating a from the second term on the RHS of (46) as follows:

V i ðaÞ ¼ 1 þ

ð46Þ

where V i (a) represent the bus i voltage series and ri (a) is the twobus-model equivalent parameter to be determined. Multiplying both sides of (46) by V i (a Þ yields:

n X V i ½kV i ½n  k  V i ½n;

ri ðaÞ V i ða Þ

ð51Þ

Thus the equation for

ÞV i ð  Þ

V i ða )

ri (a) can be obtained as:

V i ð  Þ þ i ð Þ ÞV i ð  Þ  V i ð  Þ

a ¼ ri ðaÞ ¼ V i ða

a a

r a

a

ð52Þ

The value of r[0] for this formulation is obtained as:

ri ½0 ¼ V i ½0V i ½0  V i ½0

ð53Þ

The r[n] terms in the modified formulation are given by:

ri ½n ¼

n X V i ½kV i ½n  k  V i ½n

ð54Þ

k¼0

Thus ri (a) can now be used as the parameter of an approximate two-bus equivalent network at different values of a, up to the SNBP. The SNBP of the system can be estimated by evaluating the ri (a) at escalating values of a, using its Padé approximant, until the r condition is no longer satisfied. The value of a thus obtained is the scaling factor of the loads and real-power generation corresponding to the SNBP. The extrapolating sigma method can also be used to determine the weakness of the buses in the system similar to the method explained in Section ‘‘Sigma method—for weak bus determination”. The Extrapolating Sigma Method has complexity advantages over the Non-Extrapolating Sigma Method of Section ‘‘Sigma method—for weak bus determination”. The non-extrapolating form needs a binary search, requiring multiple PF solutions to be

S.D. Rao et al. / Electrical Power and Energy Systems 84 (2017) 1–12

performed and a sigma series generated for each PF solution in order to home in on the SNBP. However with the extrapolating form, one does not need to solve multiple power-flow problems but only needs the sigma series at the base-case PF solution to be evaluated at different values of a and thus needs many fewer computations. Both sigma methods can be used to detect the weak buses of the system in a relatively quick manner and thereby identify the buses that could benefit from more refined load models. Such identification of weak buses might be used to inform the model-building step when dynamic voltage stability simulations are to be performed.

Numerical results for scaling all loads uniformly Given in this section are the results of some numerical experiments to obtain the SNBPs of different systems, estimated by different HEM-based methods and compared with the results obtained from different NR-based applications. It is difficult to assess SNBP accuracy as well-accepted industry applications often do not agree among themselves. The SNBP using the NR-based method is obtained from VSAT [53], MATPOWER [54] and PSAT [55]. While using these methods, the loads at all the buses and the real-power-generation at the PV buses are scaled uniformly until the SNBP is reached except VSAT which does not allow the generation to be scaled uniformly with the load, but instead adds an incremental amount to the generation. In order to assess the accuracy of the HEM methods in predicting the SNBP, without having the added complexity of bus-type switching, generator var limits were ignored but will be considered in a later section. The ratios of the loads at the SNBP to the base-case loads is reported as the load-scaling factor, a, corresponding to the SNBP. No execution time comparison with the other SNBP-identifying applications cited in this paper is provided for the following reasons: The HEM-based methods, PSAT and MATPOWER are implemented using MATLAB, an interpretive language which can be more than an order of magnitude slower than compiled languages, such as that used to implement VSAT. A timing comparison even among the MATLAB based application is unwise as execution times will be significantly affected by if and how sparsity is exploited and whether pivoting (diagonal, partial, complete or none) is used. While execution time is an important metric, it is not important in all situations. The primary advantage of the HEM-based methods is that they can eliminate the non-convergence issues of the traditional iterative methods, upon which the CPF and other methods of obtaining the SNBP are based. For instance, the authors have observed that for a 6057 bus model of the ERCOT system, the NR fails to converge to an operable solution if a flat-start is used, whereas it was possible to find the solution using the HEM using the algorithm described here [42]. The guarantee that HEM finds the high-voltage solution, if one exists, is usually more important than execution time when ill-conditioned cases arise. Further, HEM-based methods provide additional information about the weak buses in the system and an analytical approximation of solution over a range of loads, which other methods do not. To compare the performance of HEM and NR-based methods which are used to calculate the SNBP, the load scaling factor at the SNBP was estimated using four different HEM-based methods: i. Power-Flow Search Method (PFSM)—In this approach, a binary search is performed with the non-extrapolating formulation given by (38)–(41) until the SNBP is reached. This method is computationally the most expensive method and requires a power-flow problem to be solved for each scaled-load value, a process which is similar to the CPF.

7

ii. Padé Approximant Search (PAS)—In this approach, using the extrapolation formulation given by (4), (5), (7) and (11) the Padé approximants are calculated for each bus voltage once and then a binary search is conducted using a as the search parameter, until the SNBP is reached. This involves solving the PF problem only once and evaluating the Padé approximants for the voltage series at all buses for various a values until the binary search criteria for accuracy is satisfied. Note that it is necessary to calculate the bus power-mismatches to determine the point at which the Padé approximant breaks down. Additionally, if a sufficient number of terms to make accurate predictions near the SNBP are not included in the series, this method may under-predict the actual SNBP. iii. Roots Method—In the Roots Method, the scaling formulation given by (4), (5), (7), (11) and (12), is used to obtain the SNBP using a. The zeros of the Padé approximants. b. The poles of the Padé approximants. The zeros/poles for the voltage series of an arbitrary bus are calculated and the smallest real zero/pole is taken as the load-scaling factor at the SNBP. Note that this requires a single PF solution and the calculation of the roots of the Padé approximant for only one bus. Unlike (i) and (ii), it is unnecessary to check if the bus mismatches are being satisfied. Like (i) and (ii), the requirement exists that a sufficient number of terms in the series is needed to guarantee accuracy. Computationally, this is the most efficient of all of the HEM methods. iv. Extrapolating Sigma Method (ESM)—In the ESM, a binary search, with a as the search parameter, is conducted until the SNBP is reached using the Extrapolating Sigma Method. The value of a reported corresponds to the load-scaling factor at the SNBP. Here, the PF problem is solved only once and the Padé approximants for all the ri (a) are evaluated for multiple values of a until the r condition is violated. Note that checking for the r condition takes fewer computations than checking for the mismatch, as required in (i) and (ii). Fig. 2 provides the results for the predicted SNBP for the IEEE118 and the IEEE-300 bus systems, when the loads at all buses and the real-power generation are scaled uniformly and the var limits are ignored. Note that the results from VSAT differ slightly from the other NR-based methods because VSAT does not allow the generation to be scaled uniformly with the load, but instead adds an incremental amount to the generation, which leads to the difference being observed. A total of 41 terms were used in the power series for the methods based on the HEM. It can be seen from the results that the four methods of predicting the SNBP using HEM, compare fairly well with the NR-based methods. One of the biggest advantages of the HEM methods based on extrapolation is that the PF problem needs to be solved only once, and then only the Padé approximants have to be evaluated at different values of a. In order to determine the weakest bus of the system for the IEEE 118-bus system, the sigma series for all the buses were evaluated at increasing values of a until the sigma condition (45) was violated. It was observed that bus number 8 was the first bus to violate the sigma condition described in Section ‘‘The sigma methods”. The absolute value of the radicand for sigma given by (44) was plotted on a log scale for bus number 8 for all the values of a at which bus 8 violated the r condition (the radicand became negative), and is shown in Fig. 3. It can be seen that, although the first instance of violation occurs at a = 3.12, the value of the radicand’s small excursion into negative territory may be attributed to numerical roundoff error. The value of the radicand suddenly

8

S.D. Rao et al. / Electrical Power and Energy Systems 84 (2017) 1–12

Fig. 2. BP obtained using different methods.

increased in magnitude at a = 3.2 and, beyond that the magnitude remained large. Thus, in order to have more accurate predictions of the SNBP, the value of a at which the value of the radicand decreases below 0.01 is taken as the SNBP. Based on this criterion, a load-scaling factor of 3.2 is estimated as the SNBP for the 118-bus system. Using this same approach for the IEEE 300-bus system, 1.43 is calculated as the a corresponding to the SNBP. Of all the HE methods discussed in this paper, the roots method is computationally the least expensive and hence only the roots method will be used in subsequent sections of this paper. Fig. 4 shows the behavior of the SNBP predicted by the zeros of the Padé approximants as the number of terms in the series is increased for the IEEE 118-bus system. It can be seen that after nearly 33 terms, not much is gained from adding more terms to the series. Based on the computational complexity and the limited numerical results presented here, the Roots Method for determining the SNBP is the most reliable and efficient.

a. However, scaling all loads uniformly is an unacceptable limitation in some cases. In voltage stability studies, it is often desired to analyze cases when the loads at only a few buses are increased and by different amounts, for instance, increases in one load pocket. Methods to estimate the direction for which the scaling yields the worst-case scenario are discussed in [56,57]. In such situations, a formulation based on that described in Section ‘‘Formulation that allows extrapolation of the load” can be derived. Consider the formulation given by: N X S þ aDS Y ik V k ðaÞ ¼ i   i ; V i ða Þ k¼1

V i ½0 ¼ V sp i ;

The formulation described in Section ‘‘Formulation that allows extrapolation of the load”, allows the loads at all buses and the real-power-generation at the PV buses to be scaled by a factor of

Fig. 3. Radicand of (44) for bus 8 vs. number of terms in the series.

ð55Þ

i 2 slack

ð56Þ

N X ðPi þ jQ li Þ þ aðDP i þ jDQ li Þ  jQ gi ðaÞ Y ik V k ðaÞ ¼ ; V i ða Þ k¼1

 2  ; V i ðaÞ  V i ða Þ ¼ V sp i Formulation to allow loads at different buses to be scaled by different amounts

i2m

i2p

i2p

N X ðPi þ jQ li Þ þ aðDP i þ jDQ li Þ  jQ gi Y ik V k ðaÞ ¼ V i ða Þ k¼1

ð57Þ ð58Þ

lt

;

i2q

ð59Þ

Fig. 4. Predicted BP using roots of the numerator vs. number of terms in the series.

9

S.D. Rao et al. / Electrical Power and Energy Systems 84 (2017) 1–12

where DSi denotes the incremental complex-power injection at the PQ buses, DPi denotes the incremental real-power injection at the PV buses and DQ li represents the incremental reactive-power load at the PV buses. At a = 0, (55)–(59) reduce to: N X S Y ik V k ½0 ¼  i ; V i ½0 k¼1

V i ½0 ¼ V sp i ;

i2m

ð60Þ

i 2 slack

ð61Þ

N X ðP i þ jQ li Þ  jQ gi ½0 Y ik V k ½0 ¼ ; V i ½0 k¼1

 2  ; V i ½0  V i ½0 ¼ V sp i

i2p

ð62Þ

i2p

ð63Þ

N X ðP i þ jQ li Þ  jQ gi Y ik V k ½0 ¼ V i ½0 k¼1

lt

;

i2q

ð64Þ

Notice that solution to (60)–(64) represents the voltages and reactive-power generation at the buses for the power system with its given loading level. It is clear that (60)–(64), is a system of nonlinear equations. This nonlinear problem can be handled in the same way suggested in Appendix A, i.e., holomorphically embed the system of equations representing the germ solution, with the unknown germ voltages and reactive powers expressed as Maclaurin series and solve it using the formulation given by (38)–(41). The converged solution to this new HEM problem represents the germ. Once the germ is obtained, the remaining terms of the series can be obtained by equating the coefficients of a on both sides of (55)– (59). Since the procedure is the same as that derived in Section ‘‘Two HEM formulations”, the details of the derivation will not be presented. The recurrence relations for this formulation are given by: N X Y ik V k ½n  Si W i ½n ¼ DSi W i ½n  1;

i 2 m; n P 1

ð65Þ

k¼1

V i ½n ¼ 0;

n > 0 i 2 slack

ð66Þ

N X Y ik V k ½n  ðPi þ jQ li ÞW i ½n þ jQ gi ½nW i ½0 þ jQ gi ½0W i ½n k¼1

¼ ð DP i þ

jDQ li ÞW i ½n

 1  j

n1 X

! Q gi ½kW i ½n

 k ;

i2p

ð67Þ

k¼1

n1 X W½0V½n þ V½0W½n ¼  W½kV½n  k;

nP1

ð68Þ

k¼1

V i ½0V i ½n þ V i ½nV i ½0 ¼ 

n1 X V i ½kV i ½n  k;

n P 1; i 2 p

ð69Þ

k¼1

N X Y ik V k ½n  ðPi þ jQ li ÞW i ½n þ jQ gi lt W i ½n k¼1

¼ ðDPi þ jDQ li ÞW i ½n  1;

i2q

ð70Þ

These equations will have to be split into real and imaginary parts as in Section ‘‘Two HEM formulations”. Once the required number of terms of the power series is obtained, Padé approximants can be used to obtain the converged voltages.

Numerical results with loads at different buses scaled by different amounts This section presents results of experiments designed to study the accuracy of the SNBP, obtained from the formulation described in Section ‘‘Formulation to allow loads at different buses to be scaled by different amounts”. Since this method is intended to obtain the SNBP when the loads on only certain buses are increased, for instance in studying the behavior of unusually heavy loads in a load pocket, this method was tested by increasing loads on five weakest buses of the IEEE 118-bus and 300-bus systems, three of which were electrically close to each other. For this experiment, the PF problem was solved for the base-case scenario. Next five buses were identified by using the Extrapolating Sigma Method as weak buses by increasing the value of alpha, until at least five buses violated the sigma condition, (45), of which at least three buses were electrically close. Incremental loads were then added to these buses with the increment being proportional to the existing load at the respective buses. Generators in the vicinity of these buses (largely) supplied the incremental load by assigning an incremental generation variable to these generators. To check the accuracy of this approach, the same increments were used in VSAT and these incremental powers were scaled until the SNBP was reached. Of the discussed NR-based SNBP estimation methods, only VSAT allows scaling the incremental load and generation on a subset of buses and hence only VSAT was used to compare the performance. The SNBPs obtained from VSAT and the roots of Padé approximants described in Section ‘‘Padé approximants and branch cuts” were then compared. Table 1 provides the results for the comparison of the loadscaling factors corresponding to the SNBP obtained from VSAT and from the HEM-based methods, for the described experiment. It can be seen that the prediction of the SNBP by all the three methods is very consistent for the 300- bus system whereas for the 118bus system, VSAT predicts a SNBP at a somewhat higher loading level as compared to the HEM methods. Note that VSAT solved a number of PF problems (based on the step-size) of incremental power to reach its predicted SNBP, whereas the HEM-based methods were able to estimate the SNBP fairly accurately, by solving only one PF with 41 terms present in the power series. From a computational complexity standpoint, the extra computational expense of the HEM is offset by the repetitive PF solutions required by VSAT.

Incorporating var limits in the SNBP estimation Since lack of reactive power support can lead to a voltage collapse, it is important to account for generator var limits while estimating the SNBP. Note that among the HEM methods discussed earlier to estimate the SNBP, var limits can be incorporated in PFSM and PAS methods by performing bus-type switching while solving each power-flow problem as is done in the CPF. Currently, the reactive power limits on generators are handled in the HEM methods as follows. Let Q Gi MIN ; Q Gi MAX represent the minimum and maximum reactive power limits on a generator at bus i, respectively. The

Table 1 Comparison of VSAT with the roots method. System name

VSAT

Zeros of Padé approximants

Poles of Padé approximants

IEEE-118 IEEE-300

27.74 10.24

26.06 10.26

26.06 10.28

10

S.D. Rao et al. / Electrical Power and Energy Systems 84 (2017) 1–12

bus voltages are obtained using the algorithm presented with 41 terms in the series. Any reactive power load at bus i is added to the calculated net reactive power injection at that bus to obtain the net reactive power generated, Q Gi . If the reactive power limits are violated (Q Gi > Q Gi MAX or Q Gi < Q Gi MIN Þ, the bus type is changed from a PV bus to PQ bus with appropriate limits. For PV buses on var limits, if by reacquiring voltage control the net reactive power generated is brought within reactive power limits, then the generator bus model is switched back to a PV bus model. The authors are currently working on more elegant and efficient methods for bus type switching and tap changing that take advantage of the HEM. For the roots method however, an iterative approach has to be adopted to account for var limits which is as follows: The basecase problem (a = 1.0) is solved using HEM while ensuring that var limits are obeyed and then the roots of the Padé approximants are used to estimate the SNBP, say SNBP1. At this loading condition given by the load scaling parameter a = SNBP1, the generator var injections are calculated and a bus-type switching iteration is performed. Using these new bus-type assignments, a new HEM-based series is obtained and the roots method is used to estimate the SNBP. This process is continued until no additional bus type switching is required. The results of this approach for the scalable form given in Section ‘‘Formulation that allows extrapolation of the load” were compared with the results from PSAT when all loads and realpower generations were uniformly scaled. It was observed, for these test systems, that the HEM bus-types at the estimated SNBP differed from the PSAT bus-types at the SNBP for one or two buses, which is expected to occur since bus type switching is an ad hoc process and it is not uncommon for well accepted algorithms to lead to different sets of generators on var limits. Hence, in order to make the SNBP prediction from all algorithms comparable, the same set of buses on var limits at the SNBP—taken from the PSAT solution—was used for all methods tested in our research. It can be seen in Table 2 that the results of the roots method are consistent with the PSAT results for the IEEE 118-bus and 300-bus systems.

Table 2 Comparison of PSAT with the roots method with active var limits. System name

PSAT

Zeros of Padé approximants

Poles of Padé approximants

IEEE-118 IEEE-300

2.08 1.058

2.1 1.03

2.1 1.03

where jV i j (a) denotes the power series for the voltage magnitude at the bus given by (74).

jV i jðaÞ  jV i jðaÞ ¼ V i ðaÞ  V i ða Þ;

i2m

ð74Þ

Note that the constant impedance part of the load is not scaled by a and can be moved to the LHS by modifying the corresponding diagonal element of the admittance matrix leading to a modified ^ ik . Eq. (32) can then be substituted to admittance matrix called Y obtain N X   Y^ ik V k ðaÞ ¼ P li aðp2 jV i jðaÞ þ p3 ÞW i ða Þ k¼0

  þ jQ li aðq2 jV i jðaÞ þ q3 ÞW i ða Þ ;

i2m

ð75Þ

The recurrence relations can be obtained by equating the powers of

a on both sides of (74) and (75), given by (76) and (77). Note that V [n] depends on lower-indexed terms of jVj(a) and W(a). Once, V[n]

is obtained, jVj[n] can be obtained using (77).

! N n1 X X Y^ ik V k ½n ¼  Pli  p2 jV i j½kW i ½n  1  k þ Pli  p3 W i ½n  1 k¼0

k¼0

þ j Q li  q2

n1 X

! jV i j½kW i ½n  1  k þ Q li  q3 W i ½n  1 ;

i2m

ð76Þ

k¼0

n1 n X X 2  jV i j½0jV i j½n ¼  jV i j½kjV i j ½n  k þ V i ½kV i ½n  k k¼1

ð77Þ

k¼0

Proposed ZIP load model for the HEM

Conclusion

It is often desirable to use ZIP load models in power system analysis, instead of constant P/Q loads, in order to better represent the load behavior. This section provides a HEM formulation that incorporates polynomial ZIP load models. Consider a load bus with a ZIP load given by:

In this paper, two HEM-based full PF formulations have been proposed for identifying the SNBP and a model has been proposed for using ZIP load models. Using these formulations, four short-cut methods are proposed to estimate the SNBP. Another formulation that allows the loads at different buses to be scaled along different directions has also been proposed. Numerical results are shown that demonstrate the accuracy with which these methods predict the SNBP with and without var limits. The Roots Method is considered to be the most effective HEM-based method to estimate the SNBP. The biggest advantage of the HEM methods is that convergence is guaranteed, even at the SNBP, precision issues notwithstanding, and provided the conditions of Stahl’s theorem are satisfied [40,41]. While the traditional NR-based methods to estimate the SNBP have worked well for the cases discussed in this paper, they do not offer the guarantee of converging to a solution if one exists and thus are not guaranteed to converge at the SNBP. It is important to keep in mind that since no dynamics have been modeled, the maximum allowable loading predicted by this method does not consider Hopf bifurcations. Further research should be performed to include the equations associated with HVDC links and implement a more sophisticated algorithm to perform discrete changes such as tap-changing transformers, bus-type switching, phase-shifting transformers, and switched capacitances which is not a trivial task; however, the work presented in this paper we believe represents an important step towards using the HEM to estimate the saddle-node bifurcation point of a power system.

Pli ¼ Pli ðp1 jV i j2 þ p2 jV i j þ p3 Þ; Q li ¼ Q li ðq1 jV i j2 þ q2 jV i j þ q3 Þ;

i2m

ð71Þ

i2m

where p1 ; p2 and p3 are the constant impedance, constant current and constant power components of the real load Pli and q1 ; q2 and q3 are the constant impedance, constant current and constant power components of the reactive load Q li , respectively, at the bus i. The PBE for this load model is given by: N X ðP li ðp1 jV i j2 þp2 jV i jþp3 ÞÞ jðQ ðq jV j2 þq jV jþq ÞÞ Y ik V k ¼ þ li 1 i V  2 i 3 ; V k¼0

)

N X k¼0

i

i2m

i

    Y ik V k ¼ Pli p1 V i þ p2 jVVi jþp3 þ jQ li q1 V i þ q2 jVVi jþq3 i

ð72Þ

i

Eq. (72) can be holomorphically embedded as follows:

 N X p jV i jðaÞ þ p3 Y ik V k ðaÞ ¼ Pli p1 V i ðaÞ þ a 2   V i ða Þ k¼0  q2 jV i jðaÞ þ q3 þ jQ li q1 V i ðaÞ þ a ;  V i ða Þ

i2m

ð73Þ

11

S.D. Rao et al. / Electrical Power and Energy Systems 84 (2017) 1–12

Appendix A The nonlinear problem given by (26) can be solved using the HEM, wherein the voltages and reactive powers at a = 0 can be represented as different Maclaurin series of the complex embedding parameter b, given by (78)–(81), N X Y ik trans V k 0 ðbÞ ¼ bY i

shunt V i 0 ðbÞ;

i2m

ð78Þ

k¼1

V i 0 ðbÞ ¼ 1 þ ðV sp i  1Þb; N X Y ik

trans V k 0 ðbÞ

¼

k¼1

i 2 slack

jQ gi 0 ðbÞ  bY i V i 0 ðb Þ

ð79Þ shunt V i 0 ðbÞ;

   2  1 ; V i 0 ðbÞ  V i 0 ðb Þ ¼ 1 þ b V sp i N X Y ik k¼1

jbQ gi lt  bY i trans V k 0 ðbÞ ¼ V i 0 ðb Þ

i2p

ð80Þ

i2p

ð81Þ

V i ½0V i ½0 ¼ 1 ) V i ½0 ¼ 1 ðgermÞ  2  1 V i ½0V i ½1 þ V i ½1V i ½0 ¼ V sp i  2   1 ) V i ½1 þ V i ½1 ¼ 2V i re ½1 ¼ V sp i .. . V i ½0V i ½n þ V i ½1V i ½n  1 þ    þ V i ½n  1V i ½1 þ V i ½nV i ½0 ¼ 0 n1 X ) V i ½n þ V i ½n ¼ 2V i re ½n ¼  V½kV  ½n  k k¼1

ð89Þ shunt V i 0 ðbÞ;

i2q

ð82Þ

where V i 0 ðb ¼ 1Þ and Q gi 0 ðb ¼ 1Þ, represent the voltage and PV-bus-reactive-power injections under no-load conditions, respectively. The set of equations to obtain the germ for the above formulation is given by (83). N X Y ik trans V k 0 ½0 ¼ 0;

Equating the coefficients of the powers of a on both sides of (81) yields (90), where V i re [n] represents nth coefficient of the real part of the voltage power series from the voltage-magnitude constraint. The notation dni , as defined in (86), may be used to write a generalized expression to evaluate V i re [n]. The derivation of V i re [n], i 2 p, for an arbitrary value of n is given in (89).

i 2 m; q

ðaÞ

2 n1 V sp  1 1 X V i ½kV i ½n  k  V i re ½n ¼ dn0 þ dn1 i 2 k¼1 2

! ð90Þ

Eqs. (84), (85), (87), (88) and (90) are solved to obtain the remaining terms of the power series. Once a sufficient number of terms are obtained to meet the desired convergence tolerance, Padé approximants can be used to obtain the converged voltages and PV-bus-generated reactive powers for the germ.

k¼1 N X Y ik

trans

i V k 0 ½0 ¼ jQ V

i 0

k¼1

V i 0 ½0  V i 0 ½0 ¼ 1; V i 0 ½0 ¼ 1;

0 ½0 ; ½0

ð83Þ

i 2 p ðbÞ

i2p

ðcÞ

i 2 slack

ðdÞ

Since only the series branches are included in Y ik trans , the sum of all off-diagonal elements in any row is equal to the negative of the diagonal element. Thus a voltage of 1.0 p.u. for all of the PQ buses is consistent with (83a). Eqs. (83c) and (83d), are likewise satisfied when voltage germs for the PV buses and the slack bus are assumed to be 1.0 p.u. Finally, (83b) must be solved for the generator reactive power injection at a = 0. From inspection, we see that the reactive-power germ is 0.0. Thus, the germ for this formulation can be obtained through observation rather than calculation. Using this germ, a recursive relation can be developed to find the remaining terms of the power series by equating the coefficients of the same powers of a on both sides of (78)–(82). Equating the coefficients of the same powers of a on both sides of (78)–(80) and (82), we obtain (84), (85) and (87), (88). N X Y ik

trans V k

½n ¼ Y i

shunt V i 0 ½n

 1;

i2m

ð84Þ

k¼1

V i 0 ½n ¼ dn0 þ dn1 ðV sp i  1Þ;

i 2 slack

ð85Þ

where dni is the Kronecker delta, defined by (86).



dni ¼ N X Y ik

1; if n ¼ i 0;

ð86Þ

otherwise

n1 X Q gi 0 ½kW i 0 ½n  k trans V k 0 ½n þ jQ gi 0 ½n ¼ j

k¼1

k¼1

 Y ishunt V i 0 ½n  1; N X Y ik

!

trans V k 0 ½n

i2p

¼ jQ gi ly W i 0 ½n  1  Y ishunt V i 0 ½n  1;

ð87Þ i2q

k¼1

ð88Þ

References [1] Barbier C, Barret JP. An analysis of phenomena of voltage collapse on a transmission system. Rev Gen Elec CIGRE Special Issue; July 1980. p. 3–21. [2] Dobson Ian. The irrelevance of electric power system dynamics for the loading margin to voltage collapse and its sensitivities. Nonlinear Theory Appl, IEICE 2011;2(3):263–80. [3] Bompard E, Chicco G, Napoli R. A dynamic interpretation of the load-flow Jacobian singularity for voltage stability analysis. Int J Electr Power Energy Syst 1996;18(6):385–95. [4] Cañizares CA. Conditions for saddle-node bifurcations in AC/DC power systems. Int J Electr Power Energy Syst 1995;17(1):61–8. [5] Salgado R, Takashiba J. A framework to study QV-constraint exchange points in the maximum loadability analysis. Int J Electr Power Energy Syst 2015;64:347–55. [6] Razmi H, Shayanfar HA, Teshnehlab M. Steady state voltage stability with AVR voltage constraints. Int J Electr Power Energy Syst 2012;43(1):650–9. [7] Yang Xiaoyu, Zhou Xiaoxin, Mac Yichen, Du Zhengchun. Asymptotic numerical method for continuation power flow. Int J Electr Power Energy Syst 2012;43 (1):670–9. [8] Yang Xiaoyu, Zhou Xiaoxin, Du Zhengchun. Efficient solution algorithms for computing fold points of power flow equations. Int J Electr Power Energy Syst 2011;33(2):229–35. [9] Echavarren FM, Lobato E, Rouco L, Gómez T. Formulation, computation and improvement of steady state security margins in power systems. Part I: Theoretical framework. Int J Electr Power Energy Syst 2011;33(2):340–6. [10] Lee Ching-Yin, Tsai Shao-Hong, Wu Yuan-Kang. A new approach to the assessment of steady-state voltage stability margins using the P–Q–V curve. Int J Electr Power Energy Syst 2010;32(10):1091–8. [11] Costa Vander Menengoy da, Guedes Magda Rocha, Rosa Arlei Lucas de Sousa, Cantarino Marcelo. A new modeling of loading margin and its sensitivities using rectangular voltage coordinates in voltage stability analysis. Int J Electr Power Energy Syst 2010;32(4):290–8. [12] Echavarren FM, Lobato E, Rouco L. Steady-state analysis of the effect of reactive generation limits in voltage stability. Int J Electr Power Energy Syst 2009;79 (9):1292–9. [13] Hongjie Jia, Xiaodan Yu, Yixin Yu. An improved voltage stability index and its application. Int J Electr Power Energy Syst 2005;27(8):567–74. [14] Acharjee P. Identification of maximum loadability limit and weak buses using security constraint genetic algorithm. Int J Electr Power Energy Syst 2012;36 (1):40–50. [15] Phadke AR, Fozdar Manoj, Niazi KR. A new technique for computation of closest saddle-node bifurcation point of power system using real coded genetic algorithm. Int J Electr Power Energy Syst 2011;33(5):1203–10. [16] Jabr RA, Pal BC. Computing closest saddle node bifurcations in a radial system via conic programming. Int J Electr Power Energy Syst 2009;31(6):243–8. [17] Arya LD, Choube SC, Shrivastava M, Kothari DP. Particle swarm optimization for determining shortest distance to voltage collapse. Int J Electr Power Energy Syst 2007;29(10):796–802.

12

S.D. Rao et al. / Electrical Power and Energy Systems 84 (2017) 1–12

[18] Zhang JF, Tse CT, Wang KW, Bian XY, Wonge KP, Ho SL, Lock FS. Determination and enhancement of probabilistic stability margin with load uncertainties. Int J Electr Power Energy Syst 2013;45(1):19–27. [19] Calle IA, Castronuovo ED, Ledesma. Maximum loadability of an isolated system considering steady-state and dynamic constraints. Int J Electr Power Energy Syst 2013;53:774–81. [20] Pereira LES, da Costa VM. Interval analysis applied to the maximum loading point of electric power systems considering load data uncertainties. Int J Electr Power Energy Syst 2014;54:334–40. [21] Gu Wei, Wan Qiulan. Linearized voltage stability index for wide-area voltage monitoring and control. Int J Electr Power Energy Syst 2010;32(4):333–6. [22] Arya LD, Titare LS, Kothari DP. Probabilistic assessment and preventive control of voltage security margins using artificial neural network. Int J Electr Power Energy Syst 2007;29(2):99–105. [23] Zhao Jinquan, Chiang Hsiao-Dong, Li Hua. Enhanced look-ahead load margin estimation for voltage security assessment. Int J Electr Power Energy Syst 2004;26(6):431–8. [24] Padiyar KR, Rao SS. Dynamic analysis of voltage instability in AC–DC systems. Int J Electr Power Energy Syst 1996;18(1):11–8. [25] Liu JH, Chu CC. Wide-area measurement-based voltage stability indicators by modified coupled single-port models. IEEE Trans Power Syst 2014;29 (2):756–64. [26] Ward JB, Hale HW. Digital computer solution of power-flow problems. Power Appar Syst Part III, Trans Am Inst Electr Eng 1956;75(3):398–404. [27] Tinney W, Hart C. Power flow solution by Newton’s method. IEEE Trans Power Appar Syst 1967;PAS-86(11):1449–60. [28] Tinney W, Walker J. Direct solution of sparse network equations by optimally ordered triangular factorization. Proc IEEE 1967;55(11):1801–9. [29] Stott B, Alsac O. Fast decoupled load flow. IEEE Trans Power Appar Syst 1974; PAS-93(3):859–69. [30] Stott B. Decoupled Newton load flow. IEEE Trans Power Appar Syst 1972;PAS91(5):1955–9. [31] Stott B. Effective starting process for Newton–Raphson load flows. In: IEEE trans power app proc inst electr eng; November 1971. p. 983–7. [32] Iwamoto S, Tamura T. A load flow method for ill-conditioned power systems. IEEE Trans Power Appar Syst 1981;PAS-100(April):1736–43. [33] Schaffer MD, Tylavsky DJ. A nondiverging polar-form Newton-based power flow. IEEE Trans Indus Appl 1988;24(5):870–7. [34] Tylavsky DJ, Schaffer MD. Non-diverging power flow using a least-power type theorem. In: IEEE trans on industry applications, high voltage. October 1987. p. 944–51. [35] Thorp J, Naqavi S. Load-flow fractals draw clues to erratic behavior. IEEE Comput Appl Power 1997;10(1):59–62. [36] Tamura Y, Iba K, Iwamoto S. A method for finding multiple load-flow solutions for general power systems. In: IEEE paper no. A80 043-0 presented at the IEEE PES winter meeting. New York, NY, February 3–8; 1980.

[37] Ajjarapu V, Christy C. The continuation power flow: a tool for steady state voltage stability analysis. IEEE Trans Power Syst 1992;7(1):416–23. [38] Ghiocel SG, Chow JH. A power flow method using a new bus type for computing steady-state voltage stability margins. IEEE Trans Power Syst 2014;29(2):958–65. [39] Trias A. The holomorphic embedding load flow method. In: Power and energy society general meeting; July 2012. p. 1–8 [40] Stahl H. On the convergence of generalized Padé approximants. Construct Approx 1989;5:221–40. [41] Stahl H. The convergence of Padé approximants to functions with branch points. J Approx Theory 1997;91(2):139–204. [42] Rao S, Feng Y, Tylavsky DJ, Subramanian MK. The holomorphic embedding method applied to the power-flow problem. IEEE Trans Power Syst 2015;PP (99):1–13. [43] Trias A. System and method for monitoring and managing electrical power transmission and distribution networks. US Patents 7,519,506; 2009. 7,979,239 (2011). [44] Trias A. Sigma algebraic approximants as a diagnostic tool in power networks. US Patent 2014/0156094; 2014. [45] Trias A. Fundamentals of the holomorphic embedding load-flow method. Available from: arXiv:1509.02421v1 [cs.SY]; 2015. [47] Gamelin T. Complex analysis. Springer; 2001. [48] Nuttall J. On convergence of Padé approximants to functions with branch points. In: Saff EB, Varga RS, editors. Padé and rational approximation. New York: Academic Press; 1977. p. 101–9. [49] Baker G, Graves-Morris P. Padé approximants, ser. Encyclopedia of mathematics and its applications. Cambridge University Press; 1996. [50] Baghsorkhi SS, Suetin SP. Embedding AC power flow with voltage control in the complex plane: the case of analytic continuation via Padé approximants. Available from: arXiv:1504.03249v1 [cs.SY]; 2015. [51] Stahl H. Extremal domains associated with an analytic function I, II. Complex Variables Theory Appl: Int J 1985;4(4):311–24. 325–38. [53] Voltage security assessment tool, DSATools. [accessed in April 2015]. [54] Zimmerman RD, Murillo-Sánchez CE, Thomas RJ. MATPOWER: steady-state operations, planning and analysis tools for power systems research and education. IEEE Trans Power Syst 2011;26(1):12–9. [55] Milano F. Power systems analysis toolbox. [accessed, April 2015]. [56] Dobson I, Lu L. New methods for computing a closest saddle node bifurcation and worst case load power margin for voltage collapse. IEEE Trans Power Syst 1993;8(3):905–13. [57] Alvarado F, Dobson I, Hu Yi. Computation of closest bifurcations in power systems. IEEE Trans Power Syst 1994;9(2):918–28.