Asymptotic numerical method for continuation power flow

Electrical Power and Energy Systems 43 (2012) 670–679

Contents lists available at SciVerse ScienceDirect

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

Asymptotic numerical method for continuation power flow Xiaoyu Yang a,b,⇑, Xiaoxin Zhou b, Yichen Ma c, Zhengchun Du a a

School of Electrical Engineering, Xi’an Jiaotong University, Xi’an 710049, PR China China Electric Power Research Institute, Beijing 100192, PR China c School of Science, Xi’an Jiaotong University, Xi’an 710049, PR China b

a r t i c l e

i n f o

Article history: Received 20 December 2011 Received in revised form 7 June 2012 Accepted 9 June 2012 Available online 7 July 2012 Keywords: Asymptotic numerical method Continuation power flow Higher-order predictor Adaptive step-length control

a b s t r a c t In this paper, we use the asymptotic numerical method (ANM) to solve continuation power flow (CPF) problems. ANM can be considered as a higher-order predictor without any corrections. The method has been applied with great success to the areas of fluids, elasticity and structural mechanics. Compared to the general predictor–corrector continuation methods used in power systems, ANM has the following advantages. Firstly, the computation time is smaller. With ANM, the nonlinear problems to be solved are transformed into a recursive sequence of linear systems with the same coefficient matrix, and only one sparse Jacobian matrix factorization is required at each continuation step. Secondly, the computational procedure is automatic. A simple criterion proposed by B. Cochelin et al. can be used to determine the step-length, which makes the continuation easy, and no special step-length control strategy is required. Thirdly, as the solution branch has been expressed into a closed analytical form, the Q-limit points on P–V curves due to reactive power limits violations and other breaking points with the control devices actions can be precisely located with ANM easily. Numerical examples in several power systems were presented to validate the method. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Continuation methods based on predictor–corrector schemes have long been used to study the solution behavior of nonlinear parameter-dependent systems [1]. Since the 1990s, it has also been used in power systems for P–V curves tracing, contingency analysis, and available transfer capability calculations, etc. [2–7]. Although continuation power flow (CPF) is a powerful and robust tool in voltage stability analysis, it’s time-consuming for large scale power systems and limits its use for online applications. Extensive studies have been carried out in the CPF computation. The earlier works [2,3] used linear tangent or secant predictor in the CPF programs. In recent years, some authors [8,9] have proposed to use nonlinear predictor based on polynomial interpolations. The nonlinear predictor reduces the number of iterations in the corrector step, and enables the CPF method to take longer prediction step-length than the linear case, thus it can increase the computational speed of the CPF method. To speed up the CPF computation, other authors [10] have suggested fast decoupled CPF, which combines the properties of CPF with the virtue of fast decoupled power flow.

⇑ Corresponding author at: School of Electrical Engineering, Xi’an Jiaotong University, Xi’an 710049, PR China. Tel.: +86 10 82812645; fax: +86 10 62913126. E-mail addresses: [email protected] (X. Yang), [email protected] (X. Zhou), [email protected] (Y. Ma), [email protected] (Z. Du). 0142-0615/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijepes.2012.06.050

In terms of computation, the general predictor–corrector CPF has to use corrector process (usually by Newton-type algorithms), which brings the predicted point back to the original solution curve at each continuation step and may slow down the computation speed owing to a succession of linearizations and iterations handled in the correction procedure. Moreover, during the predictor–corrector process, it is required that the correction converge. If this is not the case, the prediction should be rejected and a smaller step-length must be reconsidered, which may also lower the computation speed. Asymptotic numerical method (ANM) appears as a continuation method without any corrector iterations [11–13]. The method is very efficient and less computationally demanding than the conventional predictor–corrector schemes. ANM is easy to be implemented when the equations of the problems are formulated with a quadratic framework. The principle of the ANM is to use a higher-order predictor [14,15] and build up the solution branch in the form of power series expansions, which are truncated at large orders. With such a representation, one gets approximations of the solution path that are very accurate inside the region of convergence. This feature has permitted to establish efficient continuation procedures without any corrections. Nowadays, ANM has been applied with great success to a wide range of problems mainly in the areas of nonlinear fluids, elasticity and structural mechanics [11–13] (in these areas, ANM is usually combined with the finite element method (FEM), which is generally used to

671

X. Yang et al. / Electrical Power and Energy Systems 43 (2012) 670–679

discretize the nonlinear systems), but few attempts have been made in power systems. In this paper, we will apply the ANM concept to solve the CPF problems. When ANM is applied to real power systems, the special aspects associated to the limits of reactive power and control devices actions (e.g., load tap changers (LTCs), switchable shunts, etc.) must be taken into account as they have a great impact on voltage stability analysis. Those aspects may generate the breaking points on the solution curves and produce the discontinuities of the curves. For ANM applications in fluids and structural mechanics areas, however, such problems are usually not encountered. The paper is organized as follows. Section 2 presents the quadratic model for parameter-dependent power flow equations. Section 3 introduces the basic principles of the ANM and applies them to the CPF problems. A simple method using ANM to locate the breaking points on the solution curves is presented in this part. Numerical examples are given in Section 4 to test the ANM continuation, and comparison results between ANM continuation and the general predictor–corrector CPF are also presented. Some conclusions are drawn in Section 5. 2. The quadratic model for parameter-dependent power flow equations Expressing the bus voltage in rectangular coordinates as Vi\hi = ei + jfi, the power flow problem can be formulated as the solution of the following quadratic equations

PGi  PLi  ei

N N X X ðGij ej  Bij fj Þ  fi ðGij fj þ Bij ej Þ ¼ 0 j¼1

Q Gi  Q Li  fi

N N X X ðGij ej  Bij fj Þ þ ei ðGij fj þ Bij ej Þ ¼ 0 j¼1

ð1Þ

j¼1

ð2Þ

f ðu; kÞ :¼ LðuÞ þ qðu; uÞ þ kBðuÞ ¼ 0

V 2is  e2i  fi2 ¼ 0

ð3Þ

where Vis is specified known voltage at bus i. A nonlinear polynomial load model is used for all the load buses to represent the voltage sensitivity of various load components [16], and we assume that the load changes according to

PLi ¼ PLi;0 ½Api ðV i =V i;0 Þ2 þ Bpi ðV i =V i;0 Þ þ C pi ð1 þ K lpi kÞ

ð4Þ

Q Li ¼ Q Li;0 ½Aqi ðV i =V i;0 Þ2 þ Bqi ðV i =V i;0 Þ þ C qi ð1 þ K lqi kÞ

ð5Þ

where Pli,0, Qli,0 and Vi,0 are the base case values of the load and voltage magnitude at bus i; Api, Bpi, Cpi, Aqi, Bqi, Cqi are the polynomial coefficients of real power and reactive power respectively, usually the coefficients satisfy Api + Bpi + Cpi = 1, Aqi + Bqi + Cqi = 1; k is a slowly varying parameter (often referred to as the loading factor) that reflects the load increase; Klpi, Klqi designates the rate of load change. Correspondingly, a specified generation is given as follows

Remark 1. We call the continuation power flow mentioned in Section 2 the general power injection CPF. The model of above system (8) can also fit for the outage continuation power flow problem (e.g., for branch’s outage) [6]. Fig. 1 shows how the admittances of an outage branch are parameterized. The diagram can be seen as a two-port network, where Yseries and Yc indicate the series and shunt admittances respectively. To avoid confusion, we use the letter l instead of k in the outage continuation case. In this model, l = 1 means the branch is in service and l = 0 out of service. The power flow equations at bus k will be redefined as (similar modifications are applied to bus m)

PGk  PLk  G0kk e2k  G0kk fk2  G0km ek em  G0km fk fm þ B0km ek fm  B0km fk em N N X X ðGkj ej  Bkj fj Þ  fk ðGkj fj þ Bkj ej Þ  Pkm ¼ 0 ek j¼1 j¼1 j–k; m j–k; m ð9Þ Q Gk  Q Lk þ B0kk e2k þ B0kk fk2  G0km fk em þ G0km ek fm þ B0km fk fm þ B0km ek em N N X X fk ðGkj ej  Bkj fj Þ þ ek ðGkj fj þ Bkj ej Þ  Q km ¼ 0 j¼1 j¼1 j–k; m j–k; m ð10Þ G0kk

0 jBkk

where þ represents the new kth diagonal element of the bus 0 admittance matrix and G0km þ jBkm represents the new (k, m)-entry of the bus admittance matrix after the branch has been removed, and Pkm, Qkm in (9) and (10) are determined by   br 2 br br br br 2 P km ¼ l Gbr ð11Þ kk ek þ Gkk fk þ Gkm ek em þ Gkm fk fm  Bkm ek fm þ Bkm fk em   br br br 2 br br 2 ð12Þ Q km ¼ l Bbr kk ek  Bkk fk þ Gkm fk em  Gkm ek fm  Bkm fk fm  Bkm ek em br

br

br where Gbr kk þ jBkk and Gkm þ jBkm are self- and mutual-admittances of the two-port network. From Eqs. (9) to (12) we know B(u) in (8) will be a quadratic function in this case.

ð6Þ

where PGi,0 is active generation at base case and Kgi specifies the rate of generation change as k varies. For load types with non-constant power loads, the bus voltage magnitude Vi is introduced as a new state variable and must satisfy

V 2i  e2i  fi2 ¼ 0

ð7Þ

We have considered the nonlinear load model and transformed the power flow equations into the quadratic forms. For highvoltage direct current (HVDC) and flexible ac transmission systems

ð8Þ

where L() is a linear operator, q(, ) is a bilinear operator, B() is a quadratic operator. The mixed unknown vector u may include the bus voltage variables ei, fi, Vi (for non-constant power loads), etc. In (8), B(u) may degenerate into a constant or a linear function if the load model is taken as constant power or constant current type.

j¼1

where PGi, QGi, PLi, QLi are the active and reactive power of generation and load of N-bus system; Yij = Gij + jBij represents the (i, j)entry of the bus admittance matrix. As for a PV bus, Eq. (2) is replaced by

PGi ¼ PGi;0 ð1 þ K gi kÞ

(FACTS) applications in power systems, those can be resolved in general by introducing auxiliary variables and different relations between the variables as we did in this paper. Similar works have been done in [17]. Readers may refer to this reference for details. The above parameter-dependent power flow equations can be rewritten in the general simple form

Fig. 1. Parameterization of the branch.

672

X. Yang et al. / Electrical Power and Energy Systems 43 (2012) 670–679

Thus, the higher-order Taylor series based predictor for the approximation of the new step can be calculated as (cf. the first order tangent prediction)

3. ANM continuation 3.1. Higher-order Taylor series based predictor

ð2Þ

The continuation methods use predictor–corrector schemes to trace the solution curves of nonlinear systems as the parameter varies. In order to improve the efficiency of the general tangent predictor, some authors have proposed the use of higher-order Taylor series expansions (e.g., second or third order) to increase the accuracy of the prediction approximation [14,15]. The methods have close relationship with the ANM continuation we focus on in this paper. The principle of pseudo-arclength continuation method proposed by H.B.Keller is described in Fig. 2 [1]. Supposing that (uj, kj) is the jth (j = 0, 1, 2, . . .) computed point on the solution curve, the next point is determined by the following augmented system:

"

GðuðsÞ; kðsÞ; sÞ ¼

f ðu; kÞ _uTj ðu  uj Þ þ k_ j ðk  kj Þ  Ds

#

¼0

ð13Þ

where s denotes the pseudo-arclength, Ds = s  sj is the step-length along the solution branch, and ðu_ j ; k_ j Þ ¼ ðduj =ds; dkj =dsÞ is the tangent vector at step j. The successive differentiation of (13) with respect to the arclength parameter s yields the following systems to determine the derivatives of (uj, kj) at different orders [14,15]: At order p = 1, the tangent vector ðu_ j ; k_ j Þ of the solution branch satisfies the following arclength conditions [1]:

(

f ju u_ j þ f jk k_ j ¼ 0 ku_ j k2 þ jk_ j j2 ¼ 1

ð14Þ

2

  ðpÞ ðpÞ are deterAt order p P 2, the higher-order derivatives uj ; kj

mined by the linear systems with the same coefficient matrix [14,15]:

"

" # #2 ðpÞ 3 nl u 4 j 5 ¼  Fp ðpÞ k_ j 0 kj

f ju u_ Tj

f

ð15Þ

_ j u_ j uu u

F nl 3 ¼ 3f

þ 2f

_ j ujð2Þ uu u

þ 3f

þ 3f

_ j k_ 2j ukk u

þf

2

ð3Þ

ðDsÞ2 þ

kj

2

uj

ðDsÞ3 þ   

6

ð17Þ

ð3Þ

2

ðDsÞ þ

kj

6

3

ðDsÞ þ   

The computation of the higher-order predictor (17) requires accurate information of the RHS vectors F nl p of (15), which appears to be rather cumbersome and expensive for general nonlinear systems from a computational viewpoint. In fact, only few leading order terms of the series can be really determined because of the complexity of the expansion procedure [14,15]. But for quadratic nonlinearity, things may be different since the computation of F nl p is cheap and easy. It is possible to use a higher-order predictor truncated at larger orders to approximate the real solution path very precisely. In this way, there is definitely no need of any corrections in the continuation process. Hence, no problems of iterations and convergence arise, and the computation speed is highly enhanced. This may lead to the following ANM continuation. 3.2. Asymptotic numerical method Theoretically, there is no restriction on the order of B(u) in (8). In general ANM continuation in the fluids and mechanics areas [11–13], B(u) is a constant vector. Here we choose B(u) to be quadratic, since it is easier to deal with F nl p uniformly when ANM is applied to power systems. We consider the power injection continuation power flow first. Define

uj

ðpÞ

ðpÞ

¼ uj =p!;

kj

ðpÞ

¼ kj =p!

ð18Þ

The Taylor series (17) truncated at order K is

_ þf

_2 kk kj ;

_ j kjð2Þ uk u

þ 3f

_ j kj uk u

uj

ð2Þ

~kjþ1 ¼ kj þ k_ j Ds þ

ðpÞ

j k

where F nl p ; p ¼ 2; 3. . . are the right hand sides (RHSs) vectors of (15), and determined by the higher-order derivatives of function f, e.g. (for the sake of clarity, the path index j of the higher derivatives of f is omitted)

F nl 2 ¼ f

~ jþ1 ¼ uj þ u_ j Ds þ u

_ j u_ j u_ j uuu u

þf

ð2Þ _ uk uj kj

_3 kkk kj

þ 3f

þ 3f

_

_ j u_ j kj uuk u

_ ð2Þ kk kj kj

.. .

ð16Þ

~ jþ1 ¼ uj þ u

K K X X ðpÞ ðpÞ uj ðDsÞp ; ~kjþ1 ¼ kj þ kj ðDsÞp p¼1

ð19Þ

p¼1

The accuracy of this approximation is of course influenced by the order of truncation K and length interval Ds, where this polynomial representation is used. Using the results in Appendixes A and B, it can be easily verified by mathematical induction that F nl p of (15) may be simplified to nl F nl p ¼ p!F p

ð20Þ

where

F nl p ¼

p1 h X

ðrÞ

ðprÞ

qðuj ; uj

ðrÞ ðprÞ

Þ þ Bju uj kj

r¼1

h   i ðrÞ ðprlÞ ðlÞ qB uj ; uj kj ðp P 2Þ

i

þ

p2 pl1 X X l¼0 r¼1

ð21Þ

ð0Þ

(Define kj ¼ kj ). In (21), Bju is the Jacobian matrix of B(u) at step j, qB (, ) is a bilinear map corresponding to the second order term of B(u). We know from (21) that F nl p are vectors depending only on the previously computed terms of the series. With (18) and (20), the sequences of linear Eq. (15) are transformed to

"

Fig. 2. The pseudo-arclength continuation.

f ju u_ Tj

" # #2 ðpÞ 3 nl uj F p 4 5¼ ðpÞ k_ j 0 kj

f

j k

ð22Þ

where p P   2. So far, we can compute all the Taylor coefficients ðpÞ ðpÞ uj ; kj in (19) easily.

673

X. Yang et al. / Electrical Power and Energy Systems 43 (2012) 670–679

At order p = 1, an efficient scheme by Keller [1] to get the tangent vector ðu_ j ; k_ j Þ from nonlinear Eq. (14) is as follows: Solve

f

j u

v ¼ f

j k

ð23Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðu_ j ; k_ j Þ ¼ rðv ; 1Þ= ð1 þ kv k22 Þ

ð24Þ

where r = ± 1 is chosen so that the orientation of the path is preserved, i.e.

u_ Tj1 u_ j þ k_ j1 k_ j > 0

v p ¼ F nl p

ð26Þ

then ðpÞ

kj

ðpÞ uj

¼

u_ Tj v p k_ j þ

u_ Tj

¼ vp þ

v

¼ k_ j u_ Tj v p

ðpÞ kj

ð27Þ

v

Once ð1 6 p 6 KÞ at each order are computed, we ~ jþ1 ; ~ obtain ðu kjþ1 Þ from (19). With ANM, the nonlinear parameterdependent power flow equations to be solved are transformed into a recursive sequence of linear systems with the same coefficient matrix (Eq. (22)), and only the RHS vectors F nl p change with the order. The computation of the Taylor series requires only one global LUfactorization of sparse Jacobian matrix f ju at each continuation step, several assemblies of the RHS vectors F nl p , and several forward and backward substitutions for no significant extra computational costs. Finally, we discuss the specific forms of the bilinear forms in (21). In practical computations, we may define the bilinear form according to the specific forms of the equations. Take the bilinear form q(, ) for example, for the ith nodal active power of Eq. (1), the bilinear form can be expressed as

"

N N X X ðGij e00j  Bij fj00 Þ þ fi0 ðGij fj00 þ Bij e00j Þ j¼1

#

j¼1

 PLi;0 Api V 0i V 00i =V 2i;0

 For branch outage continuation power flow problem, as B(u) in (8) is quadratic, then Eq. (21) may be used to calculate F nl p . The bilinear form qB( , ) is determined by Eqs. (11) and (12). 3.3. Step-length determination 3.3.1. The criterion by B. Cochelin et al. The ANM continuation allows computing the Taylor series at rather large orders. But generally the series of (19) has a finite region of convergence, and the solution path has to be defined in a step-by-step manner. The ANM continuation is therefore able to find the entire branch as a succession of different sections of branch, i.e., the evaluation of the solution at the end step becomes the starting point of the next step (Fig. 3). And it is worthwhile to establish a criterion which fixes each step-length in an optimum way. Cochelin et al. [11–13] proposed a simple criterion for determining the step-length. The criterion requires that the difference between two consecutive order solutions remains small

   ðKÞ K uj ðDsÞ  ~ jþ1;orderK  u ~ jþ1;orderK1 k1 ku 1 
ð31Þ

ð28Þ

where e is a small accuracy control number to be chosen. By approximating the denominator as ku_ j Dsk, the maximum steplength is

  1=ðK1Þ    Dsmax ¼ eku_ j k1 =uðKÞ  j 1

Once each order ð1 6 p 6 KÞ has been found, the range of validity of the series (19) is defined by the interval Ds 2 [0, Dsmax], and residual error kf(u, k)k1 will be kept at sufficient small error in the whole interval. 3.3.2. Q-limit points consideration The ANM continuation in Section 3.2 does not take into consideration the reactive power limits violations of bus generation yet, which is an important factor in voltage stability analysis [18]. It may lead to a breaking point before the SNBP, usually called the Q-limit point. At the Q-limit point, a PV bus will switch to a PQ bus, thus may cause the discontinuities of P–V curves. Several

    ul ¼ e0i ; fi0 ji R SlackBus [ V 0i ; ur     ¼ e00i ; fi00 ji R SlackBus [ V 00i :

Remark 2.  For general power injection continuation, if the loads are of constant power model, the terms containing fuk, fuuk in (A.1) will be zero, then the RHS vectors F nl p become p1   X ðrÞ ðprÞ ðp P 2Þ q uj ; uj

ð32Þ

ðpÞ uj

If we define ul, ur as the first and second dummy vectors of the bilinear form respectively, then

F nl p ¼

ð30Þ

1

ðpÞ ðpÞ ðuj ; kj Þ

 e0i

p1 h   i X ðrÞ ðprÞ ðrÞ ðprÞ q uj ; u j þ Bju uj kj ðp P 2Þ

ð25Þ

Here, ðu_ j1 ; k_ j1 Þ is the preceding tangent vector. At the initial point (u0, k0), notice that the loads are increasing at the beginning procedure, so r = + 1 is chosen for power injection continuation (Note: for branch outage continuation, r = 1 is chosen to let the admittances of branch reduce at the start). At order p P 2, we solve the sequences of linear systems of (22) with Keller’s bordering algorithm [1] in the following manner: Solve j u

F nl p ¼

r¼1

then

f

Note that the above RHS vectors are the same as the ones in general ANM continuation [11–13]. If the loads are of constant current model, the terms like fuuk will be zero, then

ð29Þ

r¼1

Fig. 3. Principles of the continuation with ANM.

674

X. Yang et al. / Electrical Power and Energy Systems 43 (2012) 670–679

methods exist and have been used to compute the Q-limit points along the solution path, e.g., sensitivity method with Newton corrections [19], Q-limit performance index guided continuation [20], etc. Next, we present a simple method to locate the Q-limit points with ANM. The reactive power generation at bus i is determined by (2) and can be expressed as

hðu; kÞ ¼ Q Li þ fi

N N X X ðGij ej  Bij fj Þ  ei ðGij fj þ Bij ej Þ j¼1

ð33Þ

j¼1

where QLi is determined by (5). Suppose that bus i is a PV bus and reactive power generation at bus i reach its limit (upper or lower) at the point (uq, kq) between the current step j and (j + 1) on the solution curve (Fig. 3). At the jth step and (j + 1)th step, it satisfies the following inequalities

hðuj ; kj Þ < Q Gi;max ;

~ jþ1 ; ~kjþ1 Þ > Q Gi;max hðu

ð34aÞ

~ jþ1 ; ~kjþ1 Þ < Q Gi;min hðu

ð34bÞ

or

hðuj ; kj Þ > Q Gi;min ; At the point (uq, kq),

hðuq ; kq Þ ¼ Q Gi;max or hðuq ; kq Þ ¼ Q Gi;min

K X ðpÞ uj ð D s q Þ p ;

kq ¼ kj þ

p¼1

K X ðpÞ kj ðDsq Þp

ð36Þ

p¼1

where Dsq is the step-length for precisely locating (uq, kq) at the jth step. Introducing the series (36) into (35), the step-length Dsq for Qlimit point (uq, kq) can be thus obtained by solving (35) directly. Naturally, a bisection method may be used here to solve a simple polynomial equation. In case multiple generators exceed their reactive limits between the jth step and (j + 1)th step, the minimum solution Dsq of (35) should be used to determine the closest Q-limit point near to the jth step. We remark here that the closest Q-limit point may probably correspond to more than one generator which violates the reactive limit simultaneously. Remark 3. In fact, there are many voltage control devices (e.g., LTCs, switchable shunts, FACTS, etc.) which may have significant impacts on voltage stability. These control devices may also induce the breaking points on P–V curves. As the solution branch has been expressed into a closed analytical form as the series (19), the breaking points produced by the control devices actions can be precisely located with ANM similar as the method for locating Qlimit points.

Remark 4. ANM can also be used to locate the accurate SNBP on the solution curves. Usually the exact SNBP of the power flow equations is computed by solving the augmented systems (e.g., the MooreSpence system or the minimally extended system [21–23]). As the first-order derivative of k with respect to the arclength changes sign at the SNBP, the step-length for the exact SNBP can be determined very easily by solving the following polynomial equation

k_ SNBP ¼

K X ðpÞ pkj ðDsSNBP Þp1 ¼ 0

The final algorithm using ANM to solve the continuation power flow problems can be described as follows: Step 0: Input power systems parameters; Choose the ANM parameters K, e. Step 1: Set step index j = 0; Solve the conventional power flow to get the initial point (u0, k0).   ðpÞ ðpÞ of the series (19) Step 2: Compute Taylor coefficients uj ; kj up to order K with (14) and (22). Step 3: Compute the range of validity of the series (19) ~ jþ1 ; ~ Ds 2 [0, Dsmax] with (32); Compute ðu kjþ1 Þ with (19). Step 4: Check if there exist any Q-limit points in the interval Ds 2 [0, Dsmax] with (34a) and (34b). If exist, then solve (35) to get the closest one near to the step j with the bisection method, change the corresponding PV bus to PQ bus, and set (uj+1, kj+1): = (uq, kq); else, set ðujþ1 ; kjþ1 Þ ~ jþ1 ; ~ :¼ ðu kjþ1 Þ. Step 5: Judge whether the tracing process will stop (e.g., the whole path has been traced or the SNBP has been passed). If not, set j: = j + 1, continue and return to step 2.

ð35Þ

where QGi,max, QGi,min are the upper and lower reactive power limits at bus i. As any points between the jth step and (j + 1)th step can be expressed as a closed-form expression like the series (19), we have

u q ¼ uj þ

3.4. Implementation of the algorithm

ð37Þ

p¼1

where DsSNBP is the step-length for precisely locating the SNBP at step j. Also, the bisection method can be applied here.

4. Numerical examples In this section, we present some examples to show the performance of ANM continuation in several power systems. The test systems include IEEE-118, IEEE-300 systems and a large scale power system in China. The IEEE systems data were taken from [24]. The large scale power system is the 1620-bus North China power system. The 1620-bus system contains 218 generators, 1594 transmission lines, and 1136 transformers. The diagram of this system is shown in [21]. The computational program was implemented in C++ codes, using doubling precision arithmetic. All tests were performed on a Pentium (R) IV microcomputer with 2.2 GHz microprocessor and 512 M RAM. We remark that in our developed ANM program, the accuracy control parameter e in (32) is actually outside of the parentheses. In fact, when we use B. Cochelin’s original step-length formula, the CPF results are similar as those in [11–13]. 4.1. Determination of K, e The truncation order K and accuracy control parameter e are important parameters for ANM continuation. Choices of K, e may be problem-dependent. It has been extensively studied for proper choices of these parameters in the areas of fluids and structural mechanics [11–13]. Next, we discuss the choices of these parameters in power systems. For the sake of convenience, the Q-limit problem is not considered temporarily (i.e., there are no PV-PQ buses transitions and the function types and dimension of the parameter-dependent power flow Eq. (8) remain the same throughout the computation). The results of ANM continuation with different e and K are presented in Tables 1 and 2. We use the general power injection continuation in the tests. In Tables 1 and 2, the step means the total number continuation steps in the whole P–V curve tracing; the residue (Res) is the maximum residual error along the solution path; and the time represents the total computational costs. A common strategy for ANM continuation is to choose a small accuracy control parameter e to limit the increase of residual error along the solution path. In general references in the areas of fluids and structural mechanics [11–13], e is chosen to be a rather small value (e.g., e = 104 or even much smaller ones). In our numerical experiments in power systems, we notice that if very small values

675

X. Yang et al. / Electrical Power and Energy Systems 43 (2012) 670–679 Table 1 Influence of the accuracy control parameter e on ANM continuation (K = 6). Systems

IEEE-118 IEEE-300 1620-bus

e = 0.01

e = 0.02

e = 0.1

e = 0.2

Step

Res (p.u.)

Time (s)

Step

Res (p.u.)

Time (s)

Step

Res (p.u.)

Time (s)

Step

Res (p.u.)

Time (s)

396 525 367

2.20  109 1.43  106 1.65  107

7.54 13.23 39.04

198 263 184

2.20  109 1.43  106 1.65  107

3.85 6.52 20.36

40 51 37

2.62  105 1.08  104 8.41  107

0.75 1.32 4.92

19 25 19

1.83  103 9.01  103 6.23  105

0.35 0.68 2.99

Table 2 Influence of the truncation order K on ANM continuation (e = 0.1). Systems

IEEE-118 IEEE-300 1620-bus

K=3

K=5

K = 10

Step

Res (p.u.)

Time (s)

Step

Res (p.u.)

Time (s)

Step

Res (p.u.)

Time (s)

Step

Res (p.u.)

Time (s)

26 33 21

2.96  103 0.16 1.97  103

0.44 0.81 2.56

36 46 33

2.35  104 1.46  103 7.28  106

0.63 1.15 4.15

49 66 46

2.20  109 1.43  106 1.65  107

1.00 2.04 8.36

56 75 52

2.20  109 1.43  106 1.65  107

1.42 3.29 13.94

are to be used, it results in a large number of continuation steps and huge computation time, and may make the ANM continuation much more inefficient. Also, we find in our numerical tests that the residual error remains the same throughout the whole continuation procedure and the accuracy of the solution may not be improved any more if e reduces to a certain level (e.g., e = 0.01,0.02). On the other hand, with larger values of e, the continuation steps decrease, but less accuracy occurs. So, we choose e = 0.1 as a proper value for power systems. When turning to the truncation order K, it is rather large, generally between 10 and 30 or even higher in ANM continuation [11–13]. For power systems, we need not choose that larger number owing to the following facts. Firstly, we find in our numerical tests that the number of continuation steps increases when the truncation order K is raised. Secondly, the accuracy of the solution may not be improved much when larger K is used (similar as the facts when e decreases). Thirdly, as in general ANM continuation, the large orders also lead to large computation of RHS vectors and more forward and backward substitutions, and naturally lead to larger computation time. Our experiences show that the truncation order K chosen to be 6 to 10 is a good choice for power systems. We have observed that the residual error remains unchanged when K is larger and e is smaller to a certain level (e.g., in Tables 1 and 2, for the 1620-bus system, the residual error is kept at 1.65  107 p.u. when e = 0.01 or K = 15). We make a brief explanation here. Indeed, by introducing the series (17) into (8), expanding it into Taylor series and ordering like powers of Ds, we get

f ðuj þ

K K X X ðpÞ ðpÞ uj ðDsÞp =p!; kj þ kj ðDsÞp =p!Þ p¼1

K = 15

p¼1

  h ð2Þ ð2Þ ¼ f ðuj ; kj Þ þ f ju u_ j þ f jk k_ j Ds þ 12 f ju uj þ f jk kj þ f juu u_ j u_ j þ 2f h i Kþ1 1 þ    þ ðKþ1Þ! f ju uðKþ1Þ þ f jk kðKþ1Þ þ F nl þ  Kþ1 ðDsÞ j j h i Kþ1 j ðKþ1Þ j ðKþ1Þ nl þ f k kj þ F Kþ1 ðDsÞ þ  ¼ f ðuj ; kj Þ þ f u uj

j _ _ uk uj kj

i

ðDsÞ2

ð38Þ

The middle terms of 1 to K orders in (38) are zero (considering (14) and (22)). The residual error of the (j + 1)th step equals jth step if higher-order terms of (38) reduce quickly. In this case, the residual error will be kept constant at the initial value of kf(u0, k0)k1 along the continuation. We did the tests with the outage continuation power flow. The similar results were obtained. Details of the information are omitted here.

We emphasize that for some values of the chosen parameters, the accuracy of the solution obtained with ANM is not satisfactory (e.g., K = 3 in Table 2 for IEEE-300 system). In order to increase the reliability of the ANM continuation, a corrector may be included in the ANM continuation [12], although the corrections are seldom used in the continuation procedure. In our ANM program, we use the general Newton correction.

4.2. The selected examples In the framework of classical predictor–corrector continuation method, it is not easy to choose an optimal step-length with which the solution converges in an optimal time-consuming way, though it is possible to use a variable step-length in the continuation procedure [14,15]. It is not the case for ANM in which the step-length is naturally adaptive and is estimated a posteriori by using the computed terms of the series. The step-length of ANM continuation is naturally larger for poorly nonlinear solutions and becomes shorter when strong nonlinearities occur. Cases A and B of IEEE118 and IEEE-300 systems, respectively, can be used to verify this point. Case C is about the power injection continuation in the 1620-bus North China power system when generators reactive power limits and load tap changers actions are considered. In these examples, the ANM parameters are set as: the truncation order K = 6 and accuracy control parameter e = 0.1.

4.2.1. Case A Figs. 4a–c show the computational results of IEEE-118 system for power injection continuation. In IEEE-118 system, the transmission lines of 20–21 and 89–90 are replaced by two HVDC lines. The load types in this example are taken as constant power load model. Assuming that the load increase with constant power factors takes place at the buses number 1–10. The Q-limit problem is not considered in this example for the moment. The full V-k (P–V) curves are obtained in 40 continuation steps, which means with 40 Jacobian matrix LU-factorizations. Fig. 4b shows that the step-length first decreases as the V-k curves advance from the start. At the 14th step, a SNBP is passed (i.e., the first-order derivative of k changes sign, the accurate SNBP occurs at k = 1.436577, by solving (37)). The step-length keeps at low values after the SNBP until the 22th step, and then it begins to increase. The smaller step-lengths correspond to the points at the right bottom of the V-k curves just after the SNBP. The step-length is automatically determined by (32) in the whole curve tracing process. The residual error is kept well under control (105 order).

676

X. Yang et al. / Electrical Power and Energy Systems 43 (2012) 670–679

1.0

1.0

0.9 0.9

V(p.u.)

V(p.u.)

0.8 0.8 V2 V5

0.7

0.7 0.6

V9002 V9012

0.5

0.6

0.4 0.5 0.3 -0.2

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

0.0

1.6

0.2

Loading factor λ

0.6

0.8

1.0

Fig. 5a. Test results of the IEEE-300 system for branch outage continuation: V versus l.

0.6

0.30

0.5

0.25

0.4

0.20

Step-length

Step-length

Fig. 4a. Test results of the IEEE-118 system for power injection continuation: V versus k.

0.4

Branch parameter μ

0.3 0.2

0.15 0.10

0.1 0.05 0.0 0.00 0

10

20

30

40

Step number

0

Fig. 4b. Test results of the IEEE-118 system for power injection continuation: steplength versus step number.

10

20

30

40

50

60

70

Step number Fig. 5b. Test results of the IEEE-300 system for branch outage continuation: steplength versus step number.

-5

3.0

×10

-5

1.8

2.5

×10

1.6 1.4 1.5

1.2

||f|| (p.u.)

||f|| (p.u.)

2.0

1.0 0.5

1.0 0.8 0.6 0.4

0.0 0

10

20

30

40

Step number Fig. 4c. Test results of the IEEE-118 system for power injection continuation: residue versus step number.

0.2 0.0 0

10

20

30

40

50

60

70

Step number Fig. 5c. Test results of the IEEE-300 system for branch outage continuation: residue versus step number.

4.2.2. Case B This example is about the branch outage continuation of IEEE300 system. Figs. 5a–c show the results when a double-line (9002–9012) contingency occurs. The power flow solution does not exist at this time since l is larger than zero during the whole solution curve tracing. The step-length reaches its high peak at the 18th step, just near the SNBP (the exact SNBP occurs at l = 0.074212, between the 19th and 20th step). As we can see from

Fig. 5a, the problem in this case takes on a strong nonlinear peculiarity with sharp turn’s occurrence at the left bottom of the V–l curves at the 26th–48th step and smaller step-lengths (<0.05) are used there. Compared with Fig. 4b, the step-lengths in Fig. 5b change dramatically in the outage continuation. The ANM continuation still works well in this example. The residual error is of 105 order in this case (Fig. 5c).

677

X. Yang et al. / Electrical Power and Energy Systems 43 (2012) 670–679 -7

0.94

×10 5.0

0.92 4.5

A(2) B(3)

0.88

||f|| (p.u.)

V(p.u.)

0.90 HuaiRou110 QingHe110

0.86

4.0 3.5 3.0 2.5

0.84 2.0 0.82

1.5 0.0

0.5

1.0

1.5

2.0

2.5

3.0

-2

0

2

4

Loading factor λ

6

8

10

12

14

16

18

Step number

Fig. 6a. Test results of the 1620-bus system with Q-limits: V versus k (‘‘o’’-Q limit points).

Fig. 6c. Test results of the 1620-bus system with Q-limits: residue versus step number.

4.2.3. Case C This example is about the power injection continuation of the 1620-bus North China power system. Figs. 6a–f show the computational results of the 1620-bus system. The 1620-bus system covers Jing-Jin-Tang (Beijing–Tianjin–Tangshan), Shanxi, Shandong, West inner-Mongolia and South Hebei subsystems, where JingJin-Tang subsystem is the main load center. In this example, we assume that the load increase is focused on several buses in JingJin-Tang subsystem and the load types of these buses are taken as constant power, constant current and constant impedance respectively. Figs. 6a–c show the results when generator reactive power limits are enforced in the V–k curves tracing. The reactive power limits are only considered on the top half of the V–k curves. In this example, it takes 17 steps to pass the SNBP (the accurate SNBP occurs at k = 2.728941). Fig. 6b shows the step-length originally calculated with (32) and the final computed step-length in the procedure of upper V–k curves tracing. In Fig. 6b, the step-lengths which are lower than the corresponding Dsmax are those for determining the Q-limit points. The last two Q-limit points (A, B) on the upper V–k curves in Fig. 6a correspond to the multiple generators which reach the reactive limits at the same time (generator numbers are in parentheses). The residual error is of 107 order in this case (Fig. 6c). Figs. 6d–f show the results of the 1620-bus system when two parallel fixed tap transformers in QingHe 220/110 kV substation in Jing-Jin-Tang subsystem are replaced by LTCs. The LTC model used in our simulation is a discrete model with no time delay

[16]. When the controlled bus voltage drops below 0.9 p.u., the LTCs begin to alter taps. In the procedure of V–k curves tracing, the LTCs act twice (corresponding two jumps on V–k curves in Fig. 6d). As the LTCs have reached the lowest taps position, they cannot continue to act to restore the bus voltage in the end. In this example, it takes 17 steps to pass the SNBP (the accurate SNBP occurs at k = 2.535994). Fig. 6e shows the step-length originally calculated with (32) and the final computed step-length in the procedure of V–k curves tracing. Similar to Fig. 6b, the step-lengths which are lower than the corresponding Dsmax in Fig. 6e are those for determining the breaking points on V–k curves (the 0th, 1th, 2th, 4th, 6th, 10th step for the Q-limit points and the 3th, 7th step for LTCs actions). The last two Q-limit points (A, B) on the upper V– k curves in Fig. 6d are multiple Q-limit points. The residual error is of 106 order in this case (Fig. 6f).

4.3. Comparison with the predictor–corrector methods Finally, we compare the performance of ANM continuation with the general predictor–corrector CPF. The predictor–corrector we used here is the pseudo-arclength continuation by Keller [1]. The reactive power limits constraints are considered in the computation. The predictor–corrector CPF uses prescribed step-lengths, and the step-lengths are adjusted properly so that the final errors for the SNBP’s approximation can be compatible with ANM continuation. The computation procedure stops when the SNBP is justly passed.

0.92

1.6

0.91

1.4 Calculated with (32) Actual step-length

B(3)

0.89

1.0

V(p.u.)

Step-length

1.2

A(2)

0.90

0.8 0.6

LTCs+Q_limts Q_limits

0.88 0.87 0.86 0.85

0.4

0.84

0.2

0.83 0.82

0.0

0.81 -2

0

2

4

6

8

10

12

14

16

18

Step number Fig. 6b. Test results of the 1620-bus system with Q-limits: step-length versus step number.

0.0

0.5

1.0

1.5

2.0

2.5

3.0

Loading factor λ Fig. 6d. Test results of the 1620-bus system with Q-limits and LTCs actions: V versus k (‘‘o’’-Q limit points).

678

X. Yang et al. / Electrical Power and Energy Systems 43 (2012) 670–679 Table 4 Comparison between ANM Continuation and predictor–corrector methods in the 1620-bus system (11 Q-limit points).

1.6 1.4

Calculated with (32) Actual step-length

Step-length

1.2

Items

ANM (K = 6, e = 0.1)

1.0 Steps number n.b.f. Q-limit points kSNBP Residue (p.u.) Time (s)

0.8 0.6 0.4 0.2

14 14 11 1.8644 1.62  106 2.50

Predictor–corrector

Ds = 0.05

Ds = 0.1

Ds = 0.2

60 75 11 1.8659 9.51  105 6.41

31 51 11 1.8631 8.18  105 4.37

18 40 12 1.9324 3.11  105 3.34

0.0 -2

0

2

4

6

8

10

12

14

16

18

points on the solution path in Table 4. We can observe from Tables 3 and 4 that:

Step number Fig. 6e. Test results of the 1620-bus system with Q-limits and LTCs actions: steplength versus step number.

-6

4.0

×10

3.5

||f|| (p.u.)

3.0 2.5 2.0 1.5 1.0 0.5 0.0 -2

0

2

4

6

8

10

12

14

16

18

Step number Fig. 6f. Test results of the 1620-bus system with Q-limits and LTCs actions: residue versus step number.

For both ANM continuation and the predictor–corrector CPF, if Keller’s bordering algorithm [1] is applied in the computation, the computational effort is basically due to the number of Jacobian matrix’s factorization (n.b.f.). Therefore, the number of Jacobian matrix’s factorization can be used as a benchmark for the comparison (but not the continuation steps). Tables 3 and 4 show the comparison results for ANM continuation and the predictor–corrector CPF. They correspond to different quantities of Q-limit points on upper V–k curves. In the Tables, the continuation steps, number of Jacobian matrix’s factorization, number of computed Q-limit points on upper V–k curves, approximate kSNBP, the residual error, and the computation times are presented. Table 3 corresponds to Case C in the previous subsection (only taking into consideration Q-limit points constraints), while we make minor modifications to the reactive limits of the generators to produce more Q-limit

Table 3 Comparison between ANM continuation and predictor–corrector methods in the 1620-bus system (6 Q-limit points). Items

ANM (K = 6, e = 0.1)

Steps number n.b.f. Q-limit points kSNBP Residue (p.u.) Time (s)

17 17 6 2.7278 4.90  107 2.68

Predictor–corrector

Ds = 0.2

Ds = 0.6

Ds = 1.0

24 48 6 2.7252 3.10  105 3.87

10 34 6 2.7226 1.55  106 3.02

11 69 6 2.7275 8.73  105 4.51

 ANM continuation is faster. Because only one sparse Jacobian matrix has to be factorized at each continuation step, ANM continuation requires much less matrix factorizations than the predictor–corrector CPF. In the examples, the CPU times of ANM continuation are smaller than those of the predictor–corrector CPF.  Step-length determination for ANM continuation is easier. The efficiency of CPF is closely related to the step-length control. In fact, ANM continuation does not require any step-length control strategy. Also, we know from above that ANM continuation may locate the Q-limit points precisely as the solution path has been expressed into a closed analytical form; while for the predictor–corrector CPF, it is more difficult to choose an optimal step-length to ensure the proper sequences of reactive limits enforcement because the shape of the solution curves and limits enforcement sequences are not known in advance. In Table 4, for example, we have to choose smaller step-lengths for the predictor–corrector CPF when more Q-limit points are on the path. When a larger step-length, e.g., Ds = 0.2 is used, this may generate a large computation error (e.g., number of Q-limit points and approximate kSNBP in the last column of Table 4).

5. Conclusions In this paper, we have applied asymptotic numerical method to the continuation power flow problems. We consider the power injection continuation and the branch outage continuation used in power systems. The load types can be nonlinear load model. In comparison with the classical predictor–corrector schemes, ANM has the following important features:  The computation time is smaller. ANM can be considered as a higher-order predictor without corrections. With ANM, the nonlinear problems to be solved are transformed into a sequence of linear systems with the same coefficient matrix, and only one sparse Jacobian matrix factorization is required at each continuation step. The RHS vectors change with the order and can be calculated with the bilinear maps with no difficulties.  Unlike the classical predictor–corrector continuation, ANM continuation does not require any special step-length control strategy. The step-lengths of ANM are naturally adaptive.  Very few parameters have to be defined in ANM continuation. Indeed, only two parameters are needed to be set, i.e., the truncation order K and the accuracy control parameter e, which may be easily determined in power systems.  With ANM, as the solution branch has been expressed into a closed analytical form, the special points like the breaking points and the SNBP on the solution curves can be easily located.

679

X. Yang et al. / Electrical Power and Energy Systems 43 (2012) 670–679

Acknowledgment

Q u ðu þ taÞ  Q u ðuÞ t qðu þ ta; Þ þ qð; u þ taÞ  qðu; Þ  qð; uÞ ¼ lim t!0 t ¼ qða; Þ þ qð; aÞ

Q uu a ¼ lim t!0

The authors thank Professor B. Cochelin and I. Charpentier from Laboratory of Mechanics and Acoustics (LMA-CNRS) in France for helpful discussions about ANM.

Hence, we obtain

Appendix A For the parameter-dependent power flow Eq. (8), the RHS vectors of (15) are given by the following equations. Note that the terms like fuuu, fuuu, fkk, fkk will be zero. We detail the computations for p = 5. F nl 2 ¼f F nl 3 F nl 4

_ j u_ j uu u

þ 2f

¼ 3f

_ j uð2Þ uu u j

¼ 4f

_ j uð3Þ uu u j

þ12f F nl 5 ¼ 5f

þ5f

þ 3f

_ j kð2Þ uk u j

þ 3f

ð2Þ ð2Þ uu uj uj

þ 10f

ð4Þ _ uk uj kj

þ10f

_

_ j kj ; uk u

_ _ j uð2Þ uuk u j kj

_ j uð4Þ uu u j

þ 6f

þ 3f

þ 4f

þ 3f uuk u_ j u_ j k_ j ;

_ j kð3Þ uk u j

þ 6f

ð2Þ ð2Þ uk uj kj

þ 4f

ð3Þ _ uk uj kj

_ j u_ j kð2Þ uuk u j ;

ð2Þ ð3Þ uu uj uj

þ 20f

ð2Þ _ uk uj kj

þ 5f

_ j kð4Þ uk u j

þ 10f

ð2Þ ð3Þ uk uj kj

þ 10f

ð3Þ ð2Þ uk uj kj

ð2Þ ð2Þ _ _ _ j uð3Þ _ ð2Þ ð2Þ uuk u j kj þ 15f uuk uj uj kj þ 30f uuk uj uj kj

_ j u_ j kð3Þ uuk u j ðA:1Þ

Appendix B The following proposition is about the second derivative of the quadratic functions and is fundamental in the procedure to simplify the expression F nl p of (15). Proposition 1. Assuming that Q(u) = q(u, u) is a set of quadratic functions, where q:Rn  Rn ? Rn is a bilinear map, the following equation is satisfied

Q uu ab ¼ qða; bÞ þ qðb; aÞ

ðB:1Þ

where Quu is the second derivative of Q with respect to u(also called Hesse tensor), and a, b 2 Rn are any two non-zero vectors. Proof of Proposition 1: The tangent operator Qu 2 L(Rn, Rn) can be expressed as

Q u ¼ qðu; Þ þ qð; uÞ

ðB:2Þ

(B.2) can be easily obtained from

Q u b ¼ lim t!0

ðB:4Þ

Q ðu þ tbÞ  Q ðuÞ qðu þ tb; u þ tbÞ  qðu; uÞ ¼ limt!0 t t ðB:3Þ

The bilinear form in (B.3) is

qðu þ tb; u þ tbÞ ¼ qðu; uÞ þ t½qðu; bÞ þ qðb; uÞ þ t 2 qðb; bÞ We get

Q u b ¼ qðu; bÞ þ qðb; uÞ Quua 2 L(Rn, Rn) is a linear map which is defined by

Q uu ab ¼ qða; bÞ þ qðb; aÞ where Quu is a symmetric bilinear map.

h

References [1] Keller HB. Lectures on numerical methods in bifurcation problems. Berlin: Springer-Verlag; 1987. [2] 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. [3] Chiang HD, Flueck AJ, Shah KS, Balu N. CPFLOW: a practical tool for tracing power system steady state stationary behavior due to load and generation variations. IEEE Trans Power Syst 1995;10(2):623–34. [4] 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. [5] Echavarren FM, Lobato E, Rouco L, Gómez T. Formulation, computation and improvement of steady state security margins in power systems. Part II: Results. Int J Electr Power Energy Syst 2011;33(2):347–58. [6] Flueck AJ, Dondeti JR. A new continuation power flow tool for investigating the nonlinear effects of transmission branch parameter variations. IEEE Trans Power Syst 2000;15(1):223–7. [7] Ejebe GC, Tong J, Waight JG, Frame JG, Wang X, Tinny WF. Available transfer capability calculations. IEEE Trans Power Syst 1998;13(4):1521–7. [8] Mori H, Yamada S. Continuation power flow with the nonlinear predictor of the Lagrange’s polynomial interpolation formula. In: IEEE PES transmission and distribution conference, Yokohama, Japan (vol. 2); October 2002. p. 1133–8. [9] Li SH, Chiang HD. Nonlinear predictors and hybrid corrector for fast continuation power flow. IET Gener Transm Distrib 2008;2(3):341–54. [10] Alves DA, da Silva LCP, Castro CA, da Costa VF. Continuation fast decoupled power flow with secant predictor. IEEE Trans Power Syst 2003;18(3):1078–85. [11] Cochelin B. A path-following technique via an asymptotic–numerical method. Comput Struct 1994;53(5):1181–92. [12] Lahmam H, Cadou JM, Zahrouni H, Damil N, Potier-Ferry M. Higher-order predictor–corrector algorithms. Int J Numer Meth Engng 2002;55(6):685–704. [13] Cochelin B, Damil N, Potier-Ferry M. Méthode Asymptotique Numérique. Paris: Hermès-Lavoisier Publications; 2007. [14] Schwetlick H, Cleve J. Higher order predictors and adaptive steplength control in path following algorithms. SIAM J Numer Anal 1987;24(6):1382–93. [15] Gervais JJ, Sadiky H. A continuation method based on a high order predictor and an adaptive steplength control. Z Angew Math Mech (ZAMM) 2004;84(8):551–63. [16] Kundur P. Power system stability and control. New York: McGraw-Hill; 1994. [17] Yu J, Yan W, Li WY, Wen LL. Quadratic models of AC–DC power flow and optimal reactive power flow with HVDC and UPFC controls. Electr Power Syst Res 2008;78(3):302–10. [18] Dobson I, Lu L. Voltage collapse precipitated by the immediate change in stability when generator reactive power limits are encountered. IEEE Trans Circuits Syst (I) 1992;39(9):762–6. [19] Hiskens IA, Chakrabarti BB. Direct calculation of reactive power limit points. Int J Electr Power Energy Syst 1996;18(2):121–9. [20] Chen K, Hussein A, Bradley ME, Wan HB. A performance index guided continuation method for fast computation of saddle-node bifurcation in power systems. IEEE Trans Power Syst 2003;18(2):753–60. [21] Yang XY, Zhou XX, Du ZC. Efficient solution algorithms for computing fold points of power flow equations. Int J Electr Power Energy Syst 2011;33(2):229–35. [22] Eidiani Mostafa. A reliable and efficient method for assessing voltage stability in transmission and distribution networks. Int J Electr Power Energy Syst 2011;33(3):453–6. [23] Ajjarapu V. Computational techniques for voltage stability assessment and control. New York: Springer; 2006. [24] Power system test case archive. .