Approximating the parameter-space stability boundary considering post-contingency corrective controls

Electric Power Systems Research 121 (2015) 313–324

Contents lists available at ScienceDirect

Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr

Approximating the parameter-space stability boundary considering post-contingency corrective controls夽 Magnus Perninge ∗ Division of Automatic Control, Lund University, Lund, Sweden

a r t i c l e

i n f o

Article history: Received 27 June 2014 Received in revised form 12 October 2014 Accepted 9 November 2014 Available online 8 December 2014 Keywords: Corrective control Loadability Optimal control Stability boundary Taylor’s expansions Voltage stability

a b s t r a c t Lately, much work in the area of voltage stability assessment has been focused on finding postcontingency corrective controls. In this article a contribution to this area will be presented where we investigate the surface of maximal loadability while allowing for post-contingency corrective controls. This objective is different from the usual, where the aim is to include the post-contingency controls in a security-constrained optimal power flow. Our aim is rather to find approximations of the postcontingency stability boundary, in pre-contingency parameter space, while including the possibility for post-contingency corrective controls. These approximations can then be used in, for example, a chanceconstrained optimal power flow routine. © 2014 Elsevier B.V. All rights reserved.

1. Introduction With the increasing amount of variable generation that we see in power systems today, knowing the stability margin in a given direction of stress is usually not sufficient for operating the system. To obtain more reliable measures of system security we want to have global information about the stability margins. This global information is expressed through the parameter-space stability boundary , that separates parameters corresponding to stable equilibrium points from those corresponding to unstable or infeasible ones. Work has already been done in approximating the stability boundary without corrective controls. An approach that will often give conservative estimates is to define the stability margin as the (minimal) distance to the stability boundary. Efficient algorithms for the computation of this distance have been proposed [1–3]. In [4,5], sensitivities of the distance to the stability boundary to changes in system parameters are given for small-signal and voltage stability, respectively. The use of the sensitivities can help the system operator take optimal actions to either steer the system away from instability or make the system return to stability following a contingency. In [6], a formula for a unified sensitivity of

夽 The author is a member of the LCCC Linnaeus Center and the eLLIIT Excellence Center at Lund University. This work was supported by the Swedish Foundation for Strategic Research through the project ICT-Psi. ∗ Tel.: +46 46 222 44 75; fax: +46 46 13 81 18. E-mail address: [email protected] http://dx.doi.org/10.1016/j.epsr.2014.11.024 0378-7796/© 2014 Elsevier B.V. All rights reserved.

the loading margin to changes in system parameters for different types of instabilities is given, and shown to give results that are consistent with the existing sensitivities presented in, for example, [5]. In [7], the stability boundary is approximated by hyperplanes from the inside, so that the approximation is conservative whenever the stability region is convex. Examples are given where the approximations are used for assessing security margins. In [8], first- and second-order approximations of the smallsignal stability boundary are presented. The authors of [8] use the implicit function theorem to express the relationship between the parameters on the stability boundary. In [9,10] it is suggested that a number of points on the stability boundary should be computed by moving in different directions in parameter space, starting from a given initial point. From these points the entire boundary is then approximated using a neural network. In [11,12], the normal to and the curvature tensor of the stability boundary are used to express second-order approximations of the voltage stability boundary, thus giving an intuitive geometrical expression of the second-order approximations. In [13] these approximations were extended to small-signal stability as well as thermal stability. A systematic way of finding appropriate points to compute the approximations in was also suggested. A special situation relating to long-term voltage stability was described in [14] and further investigated in [15]. In [16,17] second order approximations of the transient stability boundary were proposed.

314

M. Perninge / Electric Power Systems Research 121 (2015) 313–324

Although many contributions have been made to the subject of approximating parameter-space stability boundaries, none of them consider the important extension of post-contingency corrective controls. If corrective controls are not considered when dispatching the production, the system will in most situations be run in an overly conservative manner. In [18] the loadability limit surface with corrective controls, CC , was computed in a given direction of stress. The approach considers a quasi-static load recovery and assumes that the post contingency parameter trajectory only intersects the post-contingency stability boundary at one single point, called the point of intersection. In a first step, second order approximations of the post-contingency stability boundary are used to get an initial prediction of the intersection point. Then a set of necessary conditions for optimality gives a system of equations that can be solved starting with the prediction obtained in the first step. In this article we will change view to approximating the entire surface CC by a number of quadratic polynomials. These polynomials will be derived from Taylor’s expansions of the surface at points on the different smooth parts that together constitute CC . To find the normal and the curvature of CC we exploit the necessary conditions of optimality proposed in [18]. The approach taken is analogous to the one presented in [13] in the sense that we search for a point on each encountered smooth part of the stability boundary called the most important point of that surface. At these points we make series approximations of CC based on local information. The remainder of this article is organized as follows. First, in Section 2 the notion of a stability boundary is introduced. Then, in Section 3 a summary of stability limits with corrective controls, based on [18], is given. In Section 4, the post-contingency stability boundary with corrective controls, in pre-contingency parameter space, is described and corner points thereof are discussed. In Section 5, second order Taylor’s approximations for all different types of smooth parts of CC are derived. Section 6 is devoted to the more practical issue of finding points to base the approximations on and how to make sure that different approximations do not interfere with each other. In Section 7, a numerical example is given. The paper then ends with a discussion of the results of the numerical example, computational complexity of the algorithm and generalizations of the admissible controls, and a summary in Section 8. 2. The stability boundary without corrective controls To understand the notion of a stability boundary with postcontingency corrective controls we must first define the stability boundary without corrective controls. 2.1. Stability limits We align the system parameters in the vector  ∈ Rm . The system parameters can represent quantities such as load demand parameters, or active power production. Later on, we will split the vector  into a sub-vector y ∈ Rk of controllable system parameters, and a vector PL ∈ Rl of non-controllable system parameters. From an initial parameter vector p , e.g. present production and consumption, the stability limit in the direction of stress ds ∈ Sm−1 (the unit sphere in Rm ) is given by the solution to the optimization problem max{r : stable e.p. exists with  = p + rds } r∈R

Fig. 1. The power system used in the example.

2.2. The stability boundary The stability boundary surface  is the boundary of the domain wherein the system is small-signal stable. The surface  ⊂ Rm is made up of a number of different smooth manifolds [2]. Due to constraint switching there are two types of feasibility limits and we get the following different types of points on the stability boundary [13] • SNB: A saddle-node bifurcation loadability limit is a loadability limit that may occur when the system Jacobian becomes singular. This type of loadability limit is the most commonly addressed loadability limit in voltage stability assessment VSA. • SLL: Switching loadability limits [25] correspond to cases when the power system becomes immediately unstable when a control variable limit is reached. • HB: Hopf Bifurcation points are points in parameter space where the real part of one pair of complex eigenvalues of the dynamic Jacobian becomes positive as the system parameters change so that the system is no longer small-signal stable. The stability boundary is not smooth but rather made up of a number of smooth manifolds which intersect at non-smooth points that are referred to as corner points (CPs). 2.3. Example Consider the system depicted in Fig. 1. This system was analyzed in [25] and consists of three generators and one load. It is assumed that node 1 is the slack node (where all power deviations are balanced) and that the load is of the constant power type with a fixed power factor. The system has three parameters that are allowed to vary; Pg2 , Pg3 , and PLoad . It is also assumed that each generator has a limited Ef with Eflim = 2.5968 p.u. for each generator. In Fig. 2 the stability surface pre corresponding to the original system is plotted together with the stability boundary post with doubled impedance between nodes 5 and 6, when varying Pg3 and PLoad , while keeping Pg2 fixed at 1.2 p.u. As can be seen in the figure, at some level of Pg3 , between 0.5 and 1 p.u., the stability surfaces are non-smooth.

(1)

This problem can be solved by various methods such as continuation methods [19,20], optimization methods [21], direct methods [22] and quasi-steady state (QSS) simulations [23,24].

Fig. 2. The pre- and post-contingency stability surfaces in the example.

M. Perninge / Electric Power Systems Research 121 (2015) 313–324

315

The apparently non-smooth points of the stability surfaces in the figure are points where both Generator 1 and Generator 2 reach their limit on Ef simultaneously. 3. Stability limits with corrective controls In [18] the author proposes a method for deciding postcontingency stability limits, in pre-contingency parameter space, including corrective controls. The proposed method assumes a load recovery model that depends explicitly on time, and hence is independent of the corrective control. The load recovery is assumed to take place in a finite time interval [0, T]. It is also assumed that at almost all points  of CC (all points except a codimension-one subset of CC ) either  ∈ pre or the post contingency parameter trajectory touches the post-contingency stability surface at one single point. This point, denoted sIP is called the intersection point, and the corresponding time after occurrence of the contingency is denoted tIP . In this section a brief description of the proposed problem formulation and the derived necessary conditions of optimality in [18] will be given. These conditions will then serve as basis for approximating the stability boundary with corrective controls. 3.1. Contingencies and corrective controls When contingencies occur the voltage dependent loads will in general instantly drop to a lower level. If the system remains transiently stable following the contingency, load tap changers will start to increase the load-side voltage and thermostatical loads will try to increase their consumption. This leads to a phase of load recovery until, hopefully, a new equilibrium is reached. In this article the assumption made in [18] will be adopted and the following load recovery model will be applied: Pd (t) = Pc + ˇ(t)[P0 − Pc ], where P0 is the pre-contingency loading, Pc is the loading immediately following the contingency and ˇ : [0, T ] → Rl×l is a smooth matrix-valued function with ˇ(0) = 0 and ˇ(T) = I, the identity matrix. When a voltage collapse occurs it is often because the system is too heavily loaded when the contingency occurs. In this case the load recovery phase will drive the voltage past the tip on the nose curve leading to an accelerated decrease in voltage, sometimes referred to as a voltage avalanche, until the collapse is complete. Since the load recovery phase will last for a while, there is some room for applying corrective controls to save the system from collapse following a contingency. Since we assume that there is a specific recovery time, T, available for applying corrective controls, we define an admissible corrective control as a piecewise continuous function u : [0, T] → U, where U is an interval U = {z ∈ Rk : ui ≤ zi ≤ u¯ i }. Here, u can be seen as the ramping in power plants where maximal ramp rates give the upper and lower bounds on u. Furthermore, we define V as the set of all admissible controls. 3.2. Stability limits

 t The corrective control u will result in a trajectory y(t) = y0 +

u(s) ds, of the controllable system parameters. Now, let {(t, u, 0 0 ) : t ∈ [0, T]} be the trajectory corresponding to the control u when starting with the pre-contingency parameter vector 0 = (y0 , P0 ), i.e.

⎡

(t, u, 0 ) = ⎣



y0 +

⎤

t

u(s)ds 0

Pc + ˇ(t)(P0 − Pc )

⎦

(2)

Fig. 3. An optimal trajectory that intersects the stability boundary at one point.

From an initial parameter vector p , e.g. present production and consumption, the stability limit with corrective controls in the direction of stress ds ∈ Sm−1 (the unit sphere in Rm ) is given by the solution to the optimization problem max {r : stable e.p. exist for (t, u, 0 ), ∀t ∈ [0, T ] with 0

u∈V, r∈R

= p + rds }.

(3)

Note here that feasibility implies that there exists a precontingency operating point as well, otherwise the map 0 → Pc is not well defined. We might also want to add constraints limiting us to stable pre-contingency operating points. Not also that the requirement of a stable equilibrium point assumes that the load recovery is quasi-static. We let s(· , ) be an optimal trajectory that intersects post in the least possible number of points, where  point on the stability boundary CC . Fig. 3 shows an optimal trajectory with one intersection point. In the figure the horizontal axis represents the controllable parameters and the vertical axis the non-controllable (load) parameters. The point 0 corresponds to a maximum in the above optimization problem. 3.3. Necessary conditions for optimality The necessary conditions for optimality given in [18] are based on two observations. The first, denoted local feasibility, is that at tIP the distance of the optimal trajectory to the set of unstable parameter values will have a unique (local) minimum and equal zero. The second, called local optimality, is that no admissible perturbation of the control can make the trajectory corresponding to the optimal control lose contact with post . It turns out that there are three fundamentally different cases when it comes to defining necessary conditions for optimality, tIP = 0, 0 < tIP < T and tIP = T. The necessary conditions for optimality in these different cases are given below. 3.3.1. tIP = 0 If tIP = 0 then clearly no control can make the trajectory lose contact with the stability boundary. Local feasibility is fulfilled by any control that makes the initial direction point away from the stability boundary into the stable domain. Hence, the necessary conditions for optimality when tIP = 0 are trivial. 3.3.2. 0 < tIP < T When 0 < tIP < T the situation is a bit different. Now, the trajectory must touch post at one point and then be able to bend of inwards as in Fig. 3.

316

M. Perninge / Electric Power Systems Research 121 (2015) 313–324

Let uA ∈ Rk be the average of an optimal control u* on the interval [0, tIP ], i.e. 1 u = tIP



tIP

A

u∗ (t)dt,

(4)

0

and define

changes in 0 . This means that the optimality conditions are different compared to in the above case. Assume that sIP is located at a  CP of post where the smooth surfaces 1post , . . ., N post intersect. A necessary condition for optimality is that the free coordinates 1,f N ,f , are linearly dependent. Assume of the normals, npost , . . ., npost 1,f



y + tuA

A

(t, u , ) =

Pc + ˇ(t)(P − Pc )

2,f

N ,f }. The necessary conditherefore that nIP ∈ span{nIP , . . ., nIP



1,f

(5)

.

tion for optimality can then be restated as H nIP = 0, where H = [h1 . . . hnf −N +1 ] ∈ Rnf ×(nf −N +1) is a matrix whose orthogonal 2,f

The introduction of the average of the optimal control up to the intersection point gives us a way of reducing the original infinite dimensional optimization problem (3) to a finite dimensional optimization problem, where the feasibility of the entire trajectory is replaced with local feasibility at the intersection point. Furthermore, we let Kl = {i : uAi = ui }, Ku = {i : uAi = u¯ i } and Kf = {i : uAi ∈ (ui , u¯ i )}. Then, for all optimal controls u* giving the same intersection point, we have u∗i (t) = ui for all t ∈ [0, tIP ) if i ∈ Kl ,

u∗i (t) = u¯ i for all t ∈ [0, tIP ) if i ∈ Ku , and neither holds if i ∈ Kf .

N ,f . Optimality columns are orthogonal to the vectors nIP , . . ., nIP is thus given by

H nIP = 0 1,f

(11)



(nkIP ) st = 0,

k = 1, . . ., N ,

(12)

N

sIP ∈ ∪ kpost .

(13)

k=1

With vf undecided we have 2nf + 2 unknown variables and nf − N + 1 +2N = nf + N + 1 constraints, generally leaving this sys1,f

3.3.2.1. Intersection at a smooth part of post . Assume first that the intersection is at a smooth point of post . To have local optimality we must have that at sIP the normal vector, nIP , of post is perpendic (t , uA , ) ≤ ular to ui (tIP , uA , ) for all i ∈ Kf . Furthermore, n IP ui IP l u A 0, for all i ∈ K , and nIP ui (tIP , u , ) ≥ 0, for all i ∈ K . If not, we can add a small perturbation to uAi such that the result is still admissible and the trajectory loses contact with post . As this contradicts local optimality, and ui (t, uA , ) = tei , we get the conditions (nIP )i = 0,

∀i ∈ Kf ,

(6a)

(nIP )i ≥ 0,

∀i ∈ Kl ,

(6b)

(nIP )i ≤ 0,

∀i ∈ Ku .

(6c)

Local feasibility implies that the tangential intersection cons = 0, must hold. Now, the time-derivative of the straints, n IP t optimal trajectory at the intersection point, which we denote st = st (tIP , ), does not have to relate to the average control in the same way that the position does. Let v denote the control at the intersection time (v = u∗ (tIP )), so that



st =

v ˇ (t)(P0 − Pc )



.

(7)

Let Kc = Kl ∪ Ku ∪ {k + 1, k + 2, . . ., m} and define, for any vector a ∈ Rm , af to be the vector (ai )i∈Kf and let ac be the vector (ai )i∈Kc . By continuity of s(· , ) at the intersection point we have

N ,f } only span tem under-determined. However, as {nIP , . . ., nIP a N − 1 dimensional subspace, a suitable change of coordinate leaves us with nf + N + 1 unknown relevant variables and nf − N + 1 irrelevant free variables. In the same manner as the entire vector vf was irrelevant when intersection was at a smooth part of post . To define the above mentioned change of coordinates we let

⎡ ⎢ ⎢ ⎣

f

nIP = 0 (ncIP ) stc

=0

sIP ∈ post to find the nf + 2 unknowns r, tIP and

(10) uA,f .

3.3.2.2. Intersection at a CP of post . As noted in [25], when controlling the system to obtain maximal loadability we often end up in CPs of the stability surface. Above it was assumed that the intersection point is located away from a CP. We also have to give optimality conditions when sIP is located at a CP of post . When the intersection point corresponding to a 0 on a smooth part of CC is at a CP of post , sIP follows the set of CPs for small



f

If st = −Ic stc ,

(15)

where the subscript I is used to clarify that all rows are included, f indicates the columns with indexes in Kf and c the columns with indexes in Kc .   (N −1)×n f such that its Define the matrix Q = q1 . . . qN −1 ∈R  2,f

N ,f }, i.e. the rows form an orthonormal basis for span{nIP . . ., nIP columns of [Q H] form an orthonormal basis for Rnf . Then there exists a matrix  such that Q =  I!f . Hence, f

Qst = −I!c stc .

(16)

The change of coordinates



f

(9)

(14)

Separating the columns of  into free and constrained parameters we get

st → (8)

⎤ ⎥ ⎥ st = st = 0. ⎦

.. . N

f

= 0. gives This gives us the nf + 2 equations (where nf = |Kf |)



(nIP )

that vi = uAi , ∀ i ∈ Kl ∪ Ku [18]. Then, since nIP = 0, local feasibility (ncIP ) stc

(n2IP )



Q

f

st

H

is thus the sought one that lets us separate the important variables from the irrelevant ones. The matrix Q can be found by Gram-Schmidt orthogonalization 2,f N ,f which gives of the vectors nIP , . . ., nIP



  

q1

= nIP / nIP

q2

= nIP − q1 q1 , nIP

2,f



.. . qN −1

2,f

3,f

 =

N ,f nIP

3,f

  3,f   3,f / nIP − q1 q1 , nIP 

 

N −2

−

qj

j=1

N ,f qj , nIP



   N −2    N ,f    N ,f  / nIP − qj qj , nIP    j=1

M. Perninge / Electric Power Systems Research 121 (2015) 313–324

From this we get H as

 =

h1

1 −

N −2  

qj qj , 1

j=1

. . .

 =

hnf −N +1



   N −2       / 1 − qj qj , j    j=1

 

nf −N

nf −N +1 −

hj hj , nf −N +1

N −2   

−

j=1

qj qj , nf −N +1

 

/ · · · ,

j=1

where 1 , . . ., nf −N +1 ∈ Rnf are linearly independent vectors 2,f N ,f }. Now from Q =  , not in span{nIP , . . ., nIP If T 2,f N ,f  If k = nIP · · · nIP k , we find that



 



which gives qk =

.. .

 =

N −1

eN −1 −

N −2   

j

N ,f qj , nIP

j=1



   N −2    N ,f   . N −1,f / nIP − qj qj , nIP    





j=1

3.3.3. tIP = T When tIP = T, local optimality conditions are the same but local feasibility is reduced to

(nkIP ) st > 0.

(17)

Now, that tIP is fixed we have one less variable to determine. Also, we only set vf to guarantee that (17) is fulfilled. We thus have nf + 1 unknown variables, that have to be computed, in both the case of smooth intersection and when intersection is at a CP. When intersection is at a smooth point of post the conditions of optimality give f

nIP = 0

(18)

sIP ∈ post

(19)

and when intersection is at a CP of post the conditions give H nIP = 0 1,f

N

sIP ∈ ∪

k=1

kpost .

(20) (21)

3.4. Adjusting to limited production capacities In optimal power flow it is important that we are able to obtain solutions that do not violate the production capacity of any single unit. This is also the case when working with corrective controls, a generator producing its capacity limit cannot ramp up in order to save the system from collapse. To include the constraint yi ≤ (yp + rdy + tuA )i ≤ y¯ i ,

Fig. 4. The pre- and post-contingency stability surfaces and the stability surface with corrective controls.

4.1. Example continued

 2,f  = e1 / nIP 

1

317

(22)

we add the hypersurfaces yi = yi and yi = y¯ i to post . 4. The stability boundary with corrective controls As in the case of the stability boundary without corrective controls, , the post-contingency stability surface with corrective , is made up of a number of different controls for contingency i, CC i smooth manifolds. From now on we drop the subscript i, but keep in mind that we focus on a specific contingency.

In Fig. 4, pre and post are plotted together with CC for a given load recovery model and an assumed maximal ramp rate (see Section 7). In the figure the plotted trajectories are optimal trajectories for a selection of different starting points and the dashed lines show the load drop immediately following the contingency due to the decreased voltage at the load node. The non-smoothness that we can detect at several points of CC are CPs of CC . 4.2. Corner points of CC Corner points (CPs) of CC are points where CC is non-smooth. Clearly, a CP will occur whenever sIP moves from one smooth part of post to another. However, there are more ways in which a non-smoothness may arise on CC . At some places tIP may be discontinuous due to higher curvature of post than what we may obtain by controlling the trajectory. Other types of CPs will arise when tIP switches between tIP = 0, 0 < tIP < T and tIP = T. A CP will also occur whenever a control variable goes from free to constrained or vice versa. The free versus constrained variable setting CP is analogous to that of an AVR acting under voltage control mode versus under excitation limiter mode, which renders CPs of post . To summarize we have the following types of codimension-one CPs (i) Shifting between pre- and post-contingency constraints on stability. (ii) CP induced when the intersection point shifts between different smooth parts of post . (iii) Switching between different types of intersection time, tIP = 0, 0 < tIP < T and tIP = T. (iv) A control parameter switches between free and constrained. (v) Discontinuity of sIP due to high curvature of post . CPs of type i, ii, iii and iv are all easy to detect since a fundamental change in the property of the intersection point occurs when passing through the CP while moving along CC . A simple comparison of two different points of CC will in these cases reveal if they do, in fact, belong to the same smooth part of post . The same thing holds true for CPs of post , which helped the authors identify the different smooth parts of the stability boundary in [13]. CPs of type v, are not as easy to detect since the only indication we have is that tIP makes a jump as we pass through the CP. Indications of such CPs may be that a series expansion of tIP at two different points give incoherent estimates of the intersection times at the other points. Also, when moving along CC in a continuation manner we may detect such corner points by tracking the changes in tIP .

318

M. Perninge / Electric Power Systems Research 121 (2015) 313–324

5. Polynomial approximation

5.2. 0 < tIP < T

The shape of CC at a point 0 , where CC is smooth, will depend on how sIP (0 ) is affected by small changes in 0 , as well as on the shape of post at sIP (0 ). In this section we will derive formulas for the normal vector and the Weingarten map to CC , thus yielding a second order expansion. Since sIP ∈ post , we must have that for a small change  around 0 ,

What separates this case from the other two is that the intersection time changes when we move along CC . To calculate dtIP /d* we use the tangentiality constraints saying that

dsIP  ∈ TsIP post , d

−st A

(23)

where TsIP post is the tangent plane to post at sIP . Hence, dsIP  n IP d

= 0,

(24)

for all vectors  ∈ T0 CC , where nIP is the normal vector of post at sIP . We thus arrive at the conclusion that the normal vector of CC at 0 is

 ds 

nCC =˛  0

IP

d

nIP ,

(25)

where ˛ is a normalization constant whose sign makes the normal point in the direction away from the feasible domain. Now, since sIP (0 ) = s(tIP , 0 ) we get that dsIP dtIP = s∗ + st , d∗ d∗

(26)

(27)

0

Differentiating (27) w.r.t. * we get d∗



d˛ dtIP s nIP + ˛nIP s∗ + st d∗  d∗



+ ˛s

dnIP d∗

(28)

To obtain the Weingarten map dNCC we project dnCC /d* onto the basis vectors B = [b1 . . . bm−1 ], and get dN CC

= B

dst ds + n = 0. IP d∗ d∗

(31)

Combining this result with (26) while using the conditions of optimality for the two different cases will then give us dtIP /d* . 5.2.1. Intersection at a smooth part We start with the case when sIP is located at a smooth part of post . Since we are also on one of the smooth parts of CC the sets Kc and Kf will remain unchanged following a small perturbation in 0 . Partial derivatives of sc (tIP , 0 ) can therefore be computed straightforwardly. Earlier we saw how local optimality implied that f

nIP = 0.

(32)

Hence, nCC = ˛(sc ) ncIP ,  0

giving us the Weingarten map



dN CC = ˛ncIP sc

∗ ∗

+ sc

dtIP ∗ t d ∗



+ ˛(sc ) ∗

c dnIP

d∗ f

nCC = ˛s nIP . 

=

Differentiating this constraint war’s. * we get



where d/d* (and ∂/∂* ) represents differentiation with respect to local coordinates of CC at 0 giving rise to the basis {b1 . . . bm−1 } of the tangent plane. From now in we drop the argument (tIP , 0 ) of s, but keep in mind that we always evaluate s at (tIP , 0 ). If tIP = 0 or tIP = T, then dtIP /d* = 0, and when 0 < tIP < T then n s = 0. Hence, IP t we can simplify expression (25) to

dnCC

n IP st = 0.

dnCC d∗



= ˛nIP s∗ ∗ + s∗ t

 dtIP d∗

+ ˛s

∗

dnIP . d∗

(29)

Let A = − CdNIP C , where the columns of C = [c1 . . . cm−1 ] form an orthonormal basis for TsIP post and dNIP is the corresponding Weingarten map of post at sIP . Then dnIP ds = −A . d∗ d∗

AfI

ds = 0, d∗

(33)

where AfI is obtained by taking the rows of A with indices in Kf . If we further divide this matrix into Aff and Afc with columns in Kf and Kc , respectively, we get dsf dsc = −Af Afc , d∗ d∗ where Af is the matrix-inverse of Aff . Combined with (31) this gives us

(stc )





Acf Af Afc − Acc

(30)

To complete the calculation of the normal vector and the Weingarten map of CC , we need to compute s∗ , st∗ , s∗ ∗ and dtIP /d* . Since sIP is given as the solution to a set of local optimality conditions it is not surprising that also derivatives of sIP has to be divided into the three different cases tIP = 0, 0 < tIP < T and tIP = T that are treated separately. 5.1. tIP = 0 When tIP = 0, sIP (0 ) does not depend on the control, and dtIP /d* = 0. In this case dsIP /d and d2 sIP /d2 can be calculated directly from the power system equations in the pre- and postcontingency states.

f

So, for this case we do not need to compute either s  or s t . ∗ ∗ ∗ If we take the derivative of (32) w.r.t. * we get (see (30))



sc + stc ∗

dtIP d∗



+ (ncIP )





c c st + stt ∗

dtIP d∗



= 0. (34)

Now, we can compute dtIP /d* as c (stc ) [Acf Af Afc − Acc ]sc + (ncIP ) st dtIP ∗ ∗ =− . c c c d∗ (s ) [Acf Af Afc − Acc ]s + (n ) sc t

t

IP

(35)

tt

5.2.2. Intersection at a CP of post Since sIP follows the set of CPs for small changes in 0 the same reasoning as when setting st in Section 3.3 applies to ds/d, Q

dsf dsc = −Ic . d d

(36)

M. Perninge / Electric Power Systems Research 121 (2015) 313–324

Dropping the superscript 1 in n1IP we get the normal vector = ˛nTIP

nCC



5.3. tIP = T

c f ds T ds f T ds = ˛(ncIP ) + ˛(QnIP ) Q d d d

= ˛ (ncIP ) − (QnIP ) Ic f











= ˛ (ncIP ) − (QnIP ) Ic f



When tIP = T, the directional derivative st is not important, since dtIP /d* = 0. This means that

 dsc d

dtIP d

sc + stc



ds = s∗ , d∗



= ˛ (ncIP ) − (QnIP ) Ic sc + ˛ (ncIP ) st c + (QnIP ) Qst f











f

= ˛ (ncIP ) − (QnIP ) Ic sc + n s IP t f



f

 dt

IP

d

dtIP d

= ˛ (ncIP ) − (QnIP ) Ic sc , f

where the step from the first line to the second follows from (36) and the step from line three to line four-five uses (16). Differentiation and projection onto the basis {b1 , . . ., bm−1 } now yields



dN

CC

=˛

dncIP

 −

d∗

 f − QnIP



 d  d∗

dIc d Ic +  d∗ d∗

f



f



QnIP

Ic



+ ˛ (ncIP ) − (QnIP ) Ic

sc

∗ ∗

+ sc

∗t

dtIP d∗

 .

f



 f

nIP G + H AfI

ds ds =L = 0, d∗ d∗

(37)

dsf dsc dsc = −(LIf )−1 LIc = −Lf LIc . d∗ d∗ d∗



0 = (nIP ) st = (ncIP ) stc + (QnIP ) Qst . f

f

Differentiating w.r.t. * we get c dnIP

d∗



+ (ncIP )



f

f

dstc f d(QnIP ) f d(Qst ) + (Qst ) + (QnIP ) . d∗ d∗ d∗

d(Qst ) d dIc =− Ic +  d∗ d∗ d∗

max

()

(38a)

 stc + Ic

s.t.  ∈ CC . ij

(38b)

By definition, the position of the most important points will highly depend on our choice of the importance function. A suitable choice, which was proposed in [13] and further developed in [26], is to use a prediction of the probability distribution of the parameters fZ and set () = fZ (). 6.2. Moving along CC

To compute dtIP /d* we use the tangential intersection constraints which for this case can be written

, 0 is thus the solution to the For one given smooth part CC ij following optimization problem: 

where G is a tensor with local derivatives of the basis vectors h1 , . . ., hnf −N +1 . It is straightforward to compute these derivatives from the Gram-Schmidt formulas in an iterative manner. Similar to the previous case we can now obtain dsf /d* through

f

to find an adequate approximation point 0 , called in the following . the most important point, on each CC ij In order to get the most accurate approximation of the surface, ij the most important points 0 will be chosen as being the points on the smooth parts of the stability boundary that maximize a given . importance function  for each CC ij ij

Differentiating this constraint with respect to * and including the fact that we have to follow the set of CPs we get

Here,

In Section 5, analytical expressions for series approximations of each type of smooth part of CC were given. These series approximations are all computed around a given point 0 on the original stability boundary. In this section we will define the notion of most important points, which gives an intuitive way of choosing 0 to get as accurate approximations as possible. It will also be explained how to move along the surface in order to find the most important point of each smooth part of CC . One important concept is the dividing surfaces which are surfaces that separate approximations of different smooth parts. It will also show how to compute these surfaces.

ij

H nIP = 0.

0 = (stc )

6. Building the approximations

Let CC , . . ., CC be all smooth parts of CC , for contingency i. i1 iM To approximate the entire surface as accurately as possible we want , and thus, we need to get a polynomial approximation of each CC ij

∗

Since the conditions of optimality are different from the smooth intersection case we must find an equivalent of (33) to be able to find a relation between dsf /d* and dsc /d* . We have the conditions of optimality



which will simplify the equations for the normal vector and the Weingarten map of Section 5.2.

6.1. Most important points sc



319

dstc . d∗

Putting these equations together will allow us to, as in the previous case, compute dtIP /d* , and to finally retrieve dNCC .

To find the most important point of each smooth part of CC we need to be able to move along the different smooth parts and also along their edges of intersection. 6.2.1. Continuation methods Continuation methods have previously been proposed as a means of exploring the power system stability boundary. In [27], a predictor-corrector method was proposed to explore a SNB surface, where the prediction direction was chosen arbitrarily in the tangent hyperplane. The corrector step then uses the Newton–Raphson method to solve equations describing the SNB surface and project the predicted point onto the surface orthogonally to the initial tangent hyperplane. With some minor adjustments, the method from [27] can be used to explore CC as well. First we predict whereto we should

320

M. Perninge / Electric Power Systems Research 121 (2015) 313–324

Fig. 5. During the continuation search, a non-feasible solution may pass by undetected.

move with the next step by, for example, projecting the gradient of the importance function  onto the tangent hyperplane (see [28] for a similar application). p =

m−1  

6.3. Dividing surfaces



d bi ,  bi . d

Fig. 7. A highly curved part of CC destroys the approximation at a distant location.

(39)

i=1

Then in the corrector step the conditions of optimality from Section 3.3 are used to characterize the surface. 6.2.2. Detecting corner points As we move along the surface on a smooth part with the continuation method we will only recognize how the present intersection point moves along post . This may lead to unfeasible solutions. An example is shown in Fig. 5, the intersection point follows one smooth part of post and we fail to recognize that after a few steps of the continuation method, another smooth part of the boundary is breached by the optimal path. In this case we get a solution that is locally feasible, but not globally feasible. As we progress the search we will either find a candidate for the solution of (6.1) or the intersection point will reach the new smooth part of the surface. The second case is depicted in Fig. 6. In the second case we thus realize that we have passed a CP. In the first case a validation of the feasibility has to be performed. 6.2.3. Moving along corner point sets of CC In many cases the solution to (6.1) will be located at a CP of CC . Therefore, it is important that we adapt the continuation method to be able to follow sets of CPs as well. In [12,13] it was shown how to adapt the continuation method proposed in [27] to follow sets of CPs. The proposed method relies on first projecting the projected gradient in (39) onto the tangent hyperplane of the second smooth manifold and then using the characteristic equations of both surfaces to correct the guess onto the set of CPs. This method is directly applicable also when we have corrective controls, by replacing the characteristic equations of the intersecting surfaces with the optimality conditions from [18].

When there are many smooth parts of CC it may be necessary to separate the approximations from each other. Otherwise, a smaller smooth piece of CC that has a high principal curvature may disrupt the approximations at a point far away from the smooth surface (see Fig. 7). In [15] it was proposed that a dividing surface DS should be jk defined for intersections of SNB and SLL surfaces. We can extend this notion to include all types of intersecting surfaces. The idea behind the dividing surface is to, for each codimension-two set of CPs, define an additional surface that divides CC along the set of breaking points. With this praxis, the approximation of the yellow surface in Fig. 7 would only be valid to the right of the intersection of the blue and the yellow surfaces. Hence, avoiding the encircled situation. When the surfaces do not intersect tangentially we can define the normal to the dividing surface as



nDS = ˛(B(2) (B(2) ) nCC,1 − B(1) (B(1) ) nCC,2 ).

(40)

Hence, the dividing surface will be at equal angles to both the intersecting smooth parts. One advantage with this formulation is that the higher order derivatives follow directly from the equations in Section 5. However, (40) does no longer apply when the intersection is tangential, as is the case when sIP reaches a tangential intersection point (TIP) [15] of post , since then the projection of one normal onto the tangent space of the other surface is zero. In this case we have to follow the procedure described in [15] to find the normal vector and the Weingarten map for the dividing surface. 7. Numerical example In the numerical example we continue working with the system introduced in the small example of Section 2, where a contingency at time t = 0 gives a doubled impedance of the line connecting nodes 5 and 6. To define the corrective control problem we need a model of the load voltage characteristics and a model of the load recovery, in addition to the power system model. For the voltage characteristics model we assume that in the transient stage, immediately following the contingency, the active power consumed is Pc = p0 V 2 ,

Fig. 6. Eventually, the non-feasibility is detected.

where p0 is a constant determined by the post-contingency loading (p0 = P0 /V02 , with V0 the post-contingency load-node voltage) and

M. Perninge / Electric Power Systems Research 121 (2015) 313–324

321

Fig. 8. The different smooth parts of CC when the production capacity is unlimited.

V is the load-node voltage. We assume a constant power factor so that also the reactive power follows a constant impedance characteristic. To get a realistic load recovery model we assume that ˇ(t) = (1 − 0 .5t )/(1 − 0 .5T ), where the total recovery time, T, is set to 5. To limit the control it is assumed that the maximal ramp rates, in both directions, are 0.04 and 0.1 p.u. per time unit for generators 2 and 3, respectively.

Assume now that we set

⎡

0.5

⎤

c = ⎣ 0.5 ⎦ , 3 and let () = −   − c  .

(41) 4.58]

7.1. Unlimited production capacity We start by looking at the case of unlimited production capacity. The different smooth stability surfaces with corrective controls for this case, that where computed in [18], are plotted in Fig. 8. In the figure the different colors are used to indicate different smooth parts of the boundary. As we can see there are now 16 different smooth manifolds, compared to the three smooth manifolds of post plotted in Fig. 9. We thus see that the number of smooth surfaces may increase drastically as we allow corrective controls. This is not surprising since we have introduced an additional variable (time) and several new operational limits in the maximal ramp rates. This means that additional efforts have to be put on computing approximations of all the smooth parts. The reward for this is that we should, reasonably, get smaller smooth parts to approximate so that the approximations better fit the entire surface.

we search for the Starting from the point [0.95 1.4 point of the light-brown surface that maximize (41), using the continuation method first proposed in [27] and then refined in [13]. Fig. 10 shows the trajectory that our search takes. First we move along the smooth surface until we reach a set of CPs. Then we follow this set of CPs until we find a maximum of the importance function. At this point we calculate the equations of the dividing surface using (40) and get the dividing surface plotted in Fig. 11. We then make approximations of all the different smooth parts of the surface. The errors in these approximations are given in Fig. 12. 7.2. Limited production capacity In the first example we never obtained any intersection points located at CPs of post . If we add upper limits on the production capacities this is likely to happen. In any case where the capacity

Fig. 9. The different smooth parts of post .

322

M. Perninge / Electric Power Systems Research 121 (2015) 313–324

Fig. 10. Following a smooth part of the surface in search for the most important point.

Fig. 13. The post-contingency stability boundary, post , and curves of maximal loadability.

Fig. 11. Part of the surface together with part of the dividing surface. Fig. 14. The different smooth parts of CC with limited production capacity.

limit has an effect on the total loadability, the intersection point is located at a CP where one of the smooth manifolds of post (in Fig. 8) intersect the capacity limit hypersurface. We assume that P¯ g2 = 2.5 and that P¯ g3 = 1.5. In Fig. 13 the post-contingency stability boundary is plotted. The green lines in the figure are the curves of maximal loadability obtained by varying the output of one generator while keeping the generation in the other generator fixed. What is evident from the figure is that below a certain level of production in generator 2, the maximal loadability that can be reached by varying Pg3 is obtained at Pg3 = P¯ g3 . This means that, for certain starting conditions, the post-contingency corrective control will increase production in generator 3 until the production capacity is reached. The intersection point will in these cases be located along the intersection of the post-contingency stability boundary

Fig. 12. Error in the approximation.

surface and the surface Pg3 = P¯ g3 . Since this can be viewed as a corner point we thus have intersections at a CP when we limit production capacity in this way. The stability boundary with corrective controls is plotted in Fig. 14. The number of different smooth surfaces is now 13. The most interesting new surfaces are the row of surfaces closest to Pg3 = P¯ g3 (except the closest, dark blue, one). When the contingency occurs at these injected powers the intersection point is located along Pg3 = P¯ g3 . The error in the second order approximations is plotted in Fig. 15.

Fig. 15. Error in the approximation.

M. Perninge / Electric Power Systems Research 121 (2015) 313–324

8. Discussion 8.1. Results of the numerical example In the numerical example we investigated two different cases in a 7-node example system. In the first case we do not limit the production in the generators. This means that, compared to the second case in which we have a limited production capacity, we get larger smooth parts of the stability surface CC . With larger smooth surfaces to approximate we also get a worse accuracy, something that can be seen by comparing the error plotted in Fig. 12 to the error plotted in Fig. 15. Although the absolute error is larger in the first case it is only above 0.03 p.u. in very remote locations. Compared to the total loadability this is less than 0.75%.

323

It may be possible to use convex relaxations to remove the arising integer constraints, but any such analysis is out of the scope of this article. 8.4. Summary In this article we changed point of view from the local view, that was taken in [18], where the stability boundary, with corrective controls, in a given direction was sought, to a global view where the aim was to find global approximations of the stability boundary. Second order approximations of the stability boundary were derived for all the different types of intersection points. To get an adequate approximation of the entire surface we developed the notion of most important points, first proposed in [13]. In a numerical example the accuracy of the approximations was investigated for two cases in a small 7-node system.

8.2. Computational complexity References An important notion to any method that deals with power system security assessment is scalability. Since we compute approximations of the stability boundary, which is a global entity, we have the advantage of being able to do the computations off-line and replace the systems stability boundaries with the approximations when optimizing operation in real-time (or close to real-time). Nevertheless, it is important that the requirements, both in terms of computation time and memory, do not blow up with the dimension of the system. When searching along the CC for the most important point we will repeatedly need to find points on post that satisfies the optimality conditions (e.g. Section 3.3.2.1 . If a Newton algorithm is used to this purpose the main computational resources will be occupied by the repeated computation of the Weingarten map. The most demanding situation is when we are dealing with a Hopfbifurcation. In this case we need to compute the second derivative of the system Jacobian which will be a tensor of size O(n4x ), where nx is the number system variables (generator state variables + nodal voltage magnitudes and phase-angles) [13]. Fortunately, these tensors are very sparse. One example that is given in [13] is a computation in the IEEE 39-bus system where the largest tensor with 130 321 entries for a potential memory requirement of about 1 megabyte, only had 1511 nonzero elements and thus only uses around 60 kilobytes. This sparsity can also be exploited when multiplying these large tensors to possibly avoid the seemingly threatening computational burden. Once the most important point of a surface is found we yet again need the Weingarten map of post which is the most computationally expensive step in computing our approximation of CC . As this was already performed upon computing the most important point, this will not induce any additional burden. The number of smooth parts of CC may be substantially higher than the number of smooth parts of post . The approximation of each smooth part, including the search for the most important point, can, however, be done separately, so that parallelization can be used to improve the computational efficiency.

8.3. More general controls This article focuses on applying continuous controls to maintain stability following a contingency. To system operator often has additional discrete control actions to adhere to in emergency states. Such include load tap changer (LTC) blocking, shunt switching, etc. Including such controls in the analysis would give us a combinatorial optimization problem which would typically increase the computational complexity of the algorithms involved substantially.

[1] F. Alvarado, I. Dobson, Y. Hu, Computation of closest bifurcations in power systems, IEEE Trans. Power Syst. 9 (2) (1994) 918–928. [2] I. Dobson, Computing a closest bifurcation in multidimensional parameter space, J. Nonlinear Sci. 3 (1993) 307–327. [3] Y.V. Makarov, D.J. Hill, Z.Y. Dong, Computation of bifurcation boundaries for power systems: a new -plane method, IEEE Trans. Circuits Syst. 47 (2000) 536–544. [4] I. Dobson, F. Alvarado, C. DeMarco, Sensitivity of Hopf bifurcations to power system parameters, in: Proceedings of the 31st IEEE Conference on Decision and Control, IEEE, 1992, pp. 2928–2933. [5] S. Green, I. Dobson, F.L. Alvarado, Sensitivity of transfer capability margins with a fast formula, IEEE Trans. Power Syst. 17 (2002) 34–40. [6] F. Capitanescu, T. VanCutsem, Unified sensitivity analysis of unstable or low voltages caused by load increases or contingencies, IEEE Trans. Power Syst. 20 (1) (2005) 321–329. [7] Y.V. Makarov, P. Du, S. Lu, T.B. Nguyen, X. Guo, J.W. Burns, J.F. Gronquist, M.A. Pai, PMU-based wide-area security assessment: concept, method, and implementation, IEEE Trans. Smart Grid (2012) 1–8. [8] S. Yang, F. Liu, D. Zhang, S. Mei, Polynomial approximation of the small-signal stability region boundaries and its credible region in high-dimensional parameter space, Eur. Trans. Electr. Power (2012) 784–801. ˜ [9] X. Gu, C.A. Canizares, Fast prediction of loadability margins using neural networks to approximate security boundaries of power systems, IET Gener. Transm. Distrib. 3 (1) (2007) 466–475. ˜ [10] V. Gutierrez-Martinez, C.A. Canizares, C.R. Fuerte-Esquivel, A. Pizano-Martinez, X. Gu, Neural-network security-boundary constrained optimal power flow, IEEE Trans. Power Syst. 26 (1) (2011) 63–72. [11] M. Perninge, L. Söder, Risk estimation of the distance to voltage instability using a second order approximation of the saddle-node bifurcation surface, Electr. Power Syst. Res. 81 (2) (2011) 625–635. [12] M. Perninge, L. Söder, On the validity of local approximations of the power system loadability surface, IEEE Trans. Power Syst. 26 (4) (2011) 2143–2153. [13] C. Hamon, M. Perninge, L. Söder, A stochastic optimal power flow problem with stability constraints. Part I: Approximating the stability boundary, IEEE Trans. Power Syst. 28 (2) (2013) 1839–1848. [14] M. Perninge, L. Söder, Geometric properties of the loadability surface at SNB-SLL intersections and tangential intersection points, in: 16th International Conference on Intelligent System Application to Power Systems (ISAP), 2011. [15] M. Perninge, Approximating the loadability surface in the presence of SNB-SLL corner points, Electr. Power Syst. Res. 96 (2013) 64–74. [16] B. Jayasekara, U.D. Annakkage, Derivation of an accurate polynomial representation of the transient stability boundary, IEEE Trans. Power Syst. 21 (4) (2006) 1856–1863. [17] M. Perninge, J. Lavenius, L. Vanfretti, Approximating a post-contingency stable operation region in parameter space through time-domain simulation, in: Bulk Power Systems Dynamics and Control Symposium IX (IREP 2013), 2013. [18] M. Perninge, Finding points of maximal loadability considering postcontingency corrective controls, Electr. Power Syst. Res. 116 (2014) 187–200. [19] F. Milano, Power System Modelling and Scripting, Springer Verlag, London, 2010. [20] V. Ajjarapu, C. Christy, The continuation power flow: a tool for steady state voltage stability analysis, IEEE Trans. Power Syst. 7 (1) (1992) 416–423. [21] C. Vournas, M. Karystianos, N. Maratos, Bifurcation points and loadability limits as solutions of constrained optimization problems, in: Power Engineering Society Summer Meeting, IEEE, 2000. ˜ [22] C.A. Canizares, F.L. Alvarado, Point of collapse and continuation methods for large AC/DC systems, IEEE Trans. Power Syst. 8 (1993) 1. [23] T. Van Cutsem, C.D. Vournas, Voltage stability analysis in transient and midterm time scales, IEEE Trans. Power Syst. 11 (1) (1996) 146–154.

324

M. Perninge / Electric Power Systems Research 121 (2015) 313–324

[24] C.D. Vournas, N. Sakellaridis, M. Karystianos, N. Maratos, Investigating power system stability limits, in: IEEE International Symposium on Circuits and Systems (ISCAS), IEEE, Island of Kos, Greece, 2006, pp. 726–729. [25] M.E. Karystianos, N.G. Maratos, C.D. Vournas, Maximizing power-system loadability in the presence of multiple binding complementarity constraints, IEEE Trans. Circuits Syst. 54 (2007) 1775–1787.

[26] M. Perninge, Stochastic optimal power flow by multi-variate edgeworth expansions, Electr. Power Syst. Res. 109 (2014) 90–100. [27] I.A. Hiskens, R.J. Davy, Exploring the power flow solution space boundary, IEEE Trans. Power Syst. 16 (2001) 389–395. [28] M. Perninge, L. Söder, Optimal distribution of primary control participation with respect to voltage stability, Electr. Power Syst. Res. 80 (11) (2010) 1357–1363.