Finding points of maximal loadability considering post-contingency corrective controls

Electric Power Systems Research 116 (2014) 187–200

Contents lists available at ScienceDirect

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

Finding points of maximal loadability 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 10 January 2014 Received in revised form 13 May 2014 Accepted 15 June 2014 Keywords: Corrective control Hopf bifurcation Saddle-node bifurcation Switching loadability limit

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 search for maximal loadability while considering post-contingency corrective controls. This objective is different from the usual approach to the problem, where the aim is to include the post-contingency controls in a security-constrained optimal power flow. Our approach gives us an optimal control problem with a variable start point. Optimal control problems are generally very cumbersome to solve in high dimensions. However, under some mild assumptions we find that our infinite dimensional optimization problem can be transformed into a finite dimensional one. More specifically, by assuming that the load recovery is an explicit function of time we can specify a set of constraints that are necessary for optimality. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Contingency analysis is one of the most important issues that a power system operator is faced with. When long term voltage stability is of concern, neglecting post-contingency corrective controls will often lead to an overly conservative operation of the power system. However, introducing post-contingency corrective controls turns the contingency analysis problem into an optimal control problem. The objective, which without corrective controls is to find a feasible post-contingency operating point, is now to find a feasible corrective control that saves the system from losing stability. Several attempts at simplifying this problem appear in the literature. A pioneering paper on corrective control of voltage instability is [1], where a method is suggested that identifies a set of nodes where the load restoration is responsible for collapse. Minimal corrective controls are then determined based a hyperplane approximation of the loadability surface. A similar approach is taken in [2], but here a measure of the severity of unstable cases is introduced and an energy measure method is utilized to estimate the time available for corrective controls.

∗ Tel.: +46 46 222 44 75; fax: +46 46 13 81 18. E-mail addresses: [email protected], [email protected], [email protected] http://dx.doi.org/10.1016/j.epsr.2014.06.008 0378-7796/© 2014 Elsevier B.V. All rights reserved.

In the first instances of corrective security-constrained optimal power flow (CSCOPF) [3] it was proposed that constraints for a feasible post-contingency operating point, allowing for a limited amount of corrective control, should be added to the optimal power flow (OPF) problem. In [4] it was proposed that feasibility immediately following the contingency should be added to the constraints and an example illustrating the importance of this inclusion was given. In [5] a method where a constant control action is applied throughout the load recovery process is proposed. This constant control action is achieved by solving a single OPF problem. Then feasibility of the corresponding path is investigated through quasi-steady state (QSS) simulation. If feasibility is not obtained, parameters (time available for controls) are changed and the process is repeated until a feasible path is rendered or it is clear that no control action that will save the system can be found in this manner. In [6] an approach based on model predictive control (MPC) is suggested to find corrective controls in real-time operation. A number of points in time are chosen and a corrective control that solves an OPF with constraints on feasibility at these times, with an assumed load recovery model that depends explicitly on time, is computed. The first step in this control is then applied and the problem is resolved with a new recovery model based on the latest measurements of the systems evolution. In [7] a technique for maximizing the loadability limit in a specific direction of stress by tuning control parameters was

188

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

developed. The method is designed to be able to handle corner points of the loadability surface, since optimal solutions tend to be located at such points. Although much work has already been done on finding efficient post-contingency corrective controls since the introduction of this problem in [3], no really complete way of tackling the infinitedimensionality of the problem has been suggested yet. The method proposed in [5] has the appealing feature of guaranteeing feasibility of the control by checking each proposed control with a QSS simulation tool. However, it has the limiting factor of only being able to suggest constant controls. The method proposed in [6] seems efficient, but it does not give us any information of additional system loadability with corrective controls. In this article a different approach will be proposed where, instead of discretizing the system trajectory, we find the time where the corrected path meets the stability boundary. To be able to do this an explicit dependence of time will be used in the load recovery model (as in [6]). Furthermore, it will be assumed that, except for in very rare cases, there is an optimal path that only meets the stability surface at one time instance. It will be shown how the loadability in a given direction can be computed when corrective controls are allowed. We use the structure of the feasibility region and approximations of its boundary to estimate a globally optimal corrective control. The method will have its application in the planning stage where the operator of a system seeks to plan operation to be able to withstand contingencies as well as deviations from forecasted levels of uncertain system parameters. The remainder of the article is organized as follows. In the next section the problem of determining maximal loadability with corrective controls is formulated. In Section 3 the solution procedure is outlined. Then, in Section 4 a set of optimality conditions that applies at points where the system trajectory meets the postcontingency stability boundary is derived. Section 5 describes how approximations of the post-contingency stability boundary can be used to predict the optimal trajectory of the system parameters. In Section 6 it is discussed how these predictions can be corrected using the optimality conditions. A small illustrative example is then given in Section 7, before the article is concluded with a discussion of computational aspects and a summary in Section 8. 2. Problem formulation To understand how we can solve the problem of maximizing loadability by post-contingency corrective controls we must first define the stability boundary. 2.1. Loadability limits The following equality–inequality constraint representation of the power system model was introduced in [8]: (x, ) = 0

(1)

f a,i (x) · f b,i (x) = 0,

i = 1, . . ., ns

(2)

f

a,i

(x) ≥ 0,

i = 1, . . ., ns

(3)

f

b,i

(x) ≥ 0,

i = 1, . . ., ns

(4)

where is a smooth function representing the power system. fa,i and fb,i are called switching functions and are also smooth. The vector x ∈ Rn represent the power system state, in general voltage magnitudes and phase angles at all nodes of the system, and generator state variables for the generators of the system. The vector  ∈ Rm contains the system parameters which can represent quantities such as customer demand, or active power production. Here, we split the vector  into a sub-vector y ∈ Rk of controllable

system parameters, and a vector PL ∈ Rl of non-controllable system parameters. The switching constraints are introduced in order to take account of controller limits imposing constraints on the power system control parameters. One such limiting constraint is due to the overexcitation limiters in the generators of the system. In unconstrained operation the generator excitation EMF Ef is in equilibrium given by 0 = −Efi + KAi (Vref − V i ) = f a,i (x).

(5)

However, limits on the generator excitation EMF dictate that −Ef + Eflim = f b,i (x) ≥ 0.

(6)

A point in parameter space where, for some i ∈ {1, . . ., ns }, fa,i (x) = fb,i (x) = 0, is referred to as a breaking point [9] due to the shape that the PV-curve takes at such points, or a constraint switching point [10]. At such points the limit of the control variable, in this case Ef , is reached and the set of active constraints change. From an initial operating point (xp , p ) satisfying the feasibility constraints (1)–(4) the loadability limit in the direction of stress ds ∈ Sm−1 (the unit sphere in Rm ) is the solution to the equality–inequality constrained optimization problem max {r : (1) − (4) holds with  = p + rds }.

x∈Rn ,r∈R+

(7)

This problem can be solved by various methods such as continuation methods [11,12], optimization methods [9], direct methods [13] and quasi steady state (QSS) simulations [14,15]. 2.2. The stability boundary The stability boundary  is the boundary of the domain wherein the system is small-signal stable. The surface  is made up of a number of different smooth manifolds [16]. Due to constraint switching there are two types of loadability limits and we get the following different types of points on the stability boundary. • 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 [7] 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. • TL: A loadability limit corresponding to a TL occurs when the active power transfer over one line reaches the line’s thermal limit. 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 [7] 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 boundary , made up of two SLL-surfaces, is plotted when varying Pg3 and PLoad , while keeping Pg2 fixed at

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

189

and define a feasible corrective control as a piece-wise continuous function u : [0, T] → U with right-limits, 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.  t The corrective control u will result in a trajectory y(t) = y0 + u(s) ds, of the controllable system parameters. 0 2.5. Load recovery model

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

1.2 p.u. In the figure the solid line is the actual surface  and the dotted lines are constraint switching points. As can be seen in the figure, at some level of Pg3 , just below 1 p.u.,  is non-smooth. This is due to the switching of active constraints. The dashed curve corresponds to the set of breaking points where generator 2 reaches the limit on Ef , and generator 3 is already at the limit. The dotted curve corresponds to the set of breaking points where generator 2 reaches its limit on Ef , while both generator 1 and generator 3 are already at the limit. The apparently non-smooth SLL-SLL CP in the figure is thus a point where both generator 1 and generator 2 reach their limit on Ef simultaneously. 2.4. 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 thermostatically controlled loads will try to increase the loads. This leads to a phase of load recovery until, hopefully, a new equilibrium is reached. When a voltage collapse occurs it is often because the system is too heavily loaded when the contingency happens. 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. We assume that there is a specific recovery time, T, available for applying corrective controls,

The influence that corrective controls can have on loadability will highly depend on the characteristics of the load recovery process. One model that has been suggested for load recovery is the exponential recovery model [17]: x˙ p = Ps (V ) − Pd , P˙ d =

xp + Pt (V ), Tp

(8) (9)

where V is the voltage at the load bus, Ps is the static load function, Pt is the transient load function, xp is a state variable introduced in order to be able to write the load recovery model in first-order normal form and Tp is a time constant for the load recovery. With this model the recovery speed will be affected by the corrective control implicitly through the load-bus voltage. This will make optimization difficult. To prevent this, and due to the fact that load model parameters are often not known by the system operator, it was proposed in [6] to use a load model that explicitly evolves with time. In this article the same assumption will be made and the following load 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). 2.6. Stability limits with corrective controls 

] ∈ Starting from the present parameter vector p = [yp PLp m R we want to find the maximal increase in parameters in the direc tion ds = [dy dL ] ∈ S m−1 such that there is a feasible corrective control that will save the system from voltage collapse.

Fig. 2. The loadability surface from the example.

190

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

We consider a specific contingency with system equations denoted by the subscript post while the pre-contingency equations have subscript 0. The problem of maximizing loadability with postcontingency corrective controls can then be formulated as: Problem 1. Let V be the set of all piece-wise continuous functions u : [0, T] → U with right-limits. Given yp , dy , PLp and dL solve: max

r

x0 ∈Rn , x(t)∈Rn , r∈R, u∈V

subj. to 0 (x0 , y0 , P0 )

= 0,

f0a,i (x0 ) · f0b,i (x0 )

= 0,

i = 1, . . ., ns ,

f0a,i (x0 )

≥ 0,

i = 1, . . ., ns ,

f0b,i (x0 )

≥ 0,

i = 1, . . ., ns ,

= 0,

t ∈ [0, T ],

a,i b,i (x(t)) · fpost (x(t)) fpost

= 0,

t ∈ [0, T ], i = 1, . . ., ns ,

a,i (x(t)) fpost

≥ 0,

t ∈ [0, T ], i = 1, . . ., ns ,

b,i fpost (x(t))

≥ 0,

t ∈ [0, T ], i = 1, . . ., ns ,

R{j (t)}

≤ 0,

t ∈ [0, T ], j = 1, . . ., M(t),

y0

= yp + rdy ,

P0

= PLp + rdL ,

Pc

= g(x0 , x(0), P0 ).

post (x(t), y(t), Pd (t))

where  j (t) is the jth eigenvalue of the system Jacobian at equilibrium at time t and M(t) is the number of dynamic state variables t at time t, y(t) = y0 + 0 u(s) ds, g : Rn × Rn × Rm → Rm is a function that is derived from the voltage dependence of the loads and Pd (t) = Pc + ˇ(t)[P0 − Pc ] for a given smooth function ˇ : [0, T ] → Rm×m . The load-flow feasibility constraints in the optimization problem mean that we require that the parameter vector always remains within the domain of parameter vectors corresponding to stable operating points. 3. Solution outline The method for solving Problem 1 can be divided into two steps. In the first step we use approximations of the stability boundary to predict the shape of the corrective control employed in the optimum and the corresponding value of r (Section 5). Then, in the second step, we identify the points of the stability boundary that correspond to the binding constraint using a set of necessary constraints for optimality (Section 4). In the following we drop the subscript for the given contingency and denote by pre the pre-contingency stability boundary and by post the post-contingency stability boundary. Furthermore, we let CC be the set of points (p + r* ds ) in parameter space that correspond to an optimum of Problem 1. Let {s(t, ) : t ∈ [0, T]} be an optimal trajectory corresponding to  ∈ CC . We have that

⎡

s(t, ) = ⎣



⎤

t ∗

y0 +

u (s) ds 0

⎦

(10)

Pc + ˇ(t)(P0 − Pc ) Proposition 1. For any  ∈ CC , either  ∈ pre , or the parameter trajectory, s(· , ), intersects post in at least one point. This follows from the fact that s(t, p + rds ) is continuous in r for a given control u.

Fig. 3. The method proposed in this article.

Now, assume that s(· , ) is an optimal path with least possible number of intersections. Denote the set of times that s(t, ) intersects post by . Hence,  = {t ∈ [0, T] : s(t, ) ∈ post }. We denote such a time, called time of intersection, by tIP . We name the corresponding point of parameter space an intersection point, denoted sIP , i.e. sIP () = s(tIP , ), for tIP ∈ . At certain points  ∈ CC we may get a path that either intersects post at several time 1 , . . ., t N , or at one time instance while we also have instances tIP IP that  ∈ pre . To solve Problem 1 we assume that either  ∈ pre and  =∅, or ∈ / pre and || = 1. The idea is then to find sIP and tIP , which in turn will give r* , by first approximating the optimal corrective control using second order approximations of post and then solving a set of optimality constraints to find the intersection point. A flow chart of the proposed algorithm is given in Fig. 3. 3.1. Example continued Assume now that the active power production in generator 3 can be used for correct control. Fig. 4 shows a trajectory s(· , ) that intersects post at a single point.

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

191

5.2

Σ pre

5.1 5 4.9

Σ post

4.8 4.7 4.6 4.5 4.4 4.3 4.2 0.5

1

1.5

2

P g3

2.5

Fig. 4. The pre- and post-contingency stability surfaces and an optimal trajectory.

In Fig. 5 a trajectory that intersects post on two different occasions is plotted. If we make a small change of the initial point to the left (decreasing Pg3 ) along CC we would move away from the rightmost stability boundary and end up with only one intersection point. This intersection point will be located on the leftmost smooth part of the stability boundary. If we instead make a small change to the right we also end up with one intersection point but this time on the rightmost stability boundary. Hence, the points of CC corresponding to a trajectory that intersects post at two different points in this way will here be a codimension one set in CC . This will later motivate the assumption that when we search for the maximal loadability, with corrective controls, in a given direction we may assume that the trajectory only intersects post once.

domain. In the second case we may end up in one of three different situations: 1. tIP = 0, 2. 0 < tIP < T, or 3. tIP = T. We now formulate local feasibility constraints and derive necessary optimality criteria for the three different situations. 4.1. Local feasibility To define what we mean by local feasibility we let h(t) = (npost (c (t))) (c (t) − s(t, 0 )), where c (t) is the point of post closest to s(t, 0 ) and npost (c (t)) is the outward normal of post at c (t). By this definition h(t) is the signed distance between s(t, 0 ) and post , i.e. it tells us how far inside the stable operation domain s(t, 0 ) is. We say that s(· , 0 ) is locally feasible at sIP whenever h has a local minimum at tIP . For the different situations listed above, local feasibility implies the following.

4. Local feasibility and optimality Assume now that our search direction ds is pointing towards a point 0 where either 0 ∈ pre or s(· , 0 ) intersects post at one specific time instance, tIP . In the first case we only have to validate that the proposed trajectory is inside the post-contingency stable

46 44

2 s IP

42 4.4 38 36 34

1 s IP

32 4.3 28 26 0.55

0.6

0.65

0.7

0.75

0.8

P g3

0.85

0.9

Fig. 5. A trajectory with two intersection points.

0.95

1

1.05

192

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

87

86

85

84

83

82

81

0.03

0.04

0.05

0.06

0.07

P g3

0.08

0.09

0.1

0.11

0.12

Fig. 6. A trajectory for which tIP = 0.

4.1.1. tIP = 0 When tIP = 0 the critical point of the post-contingency trajectory is that immediately following the contingency (see Fig. 6). To be feasible, the post-contingency trajectory will have to point away from post . A local criteria for feasibility is thus that n post st (tIP , 0 ) ≤ 0,

4.1.3. tIP = T When s(· , 0 ) intersects post at time T local feasibility conditions dictate that n post st (tIP , 0 ) ≥ 0

(11)

since the path has to reach the post-contingency stability boundary from within the stable region (see Fig. 8).

where npost is the outward normal to post at s(tIP , 0 ). 4.1.2. 0 < tIP < T If 0 < tIP < T, the path touches post at one point of the trajectory and then bends off inwards again (see Fig. 7). At such a point n post st (tIP , 0 ) = 0.

4.2. Local optimality We start by assuming that post is smooth at the intersection point. Hence, we assume

(12)

Assumption 1. We assume that either s(t, ) intersects post at exactly one point t ∈ [0, T], or  ∈ pre and s(· , ) does not intersect post . Furthermore, we assume that post is smooth at sIP .

Furthermore, the curvature of s(t, 0 ) at tIP in the direction npost cannot be larger than the curvature of post along st (tIP , 0 ), i.e. [18] T n post stt (tIP , 0 ) + st Apost st ≤ 0.  Here, Apost = −Cpost dN post C

(13) post

(14)

By this assumption the function h(t) corresponding to the optimal path will have a global unique minimum at tIP , where h(tIP ) = 0. Our definition of local optimality will be with respect to small changes in tIP and y(tIP ). The local optimality constraints stated

, where Cpost ∈ Rm×m−1 is a

basis for the tangent plane of post at s(tIP , 0 ) and dN post is the corresponding curvature tensor. 4.7

68

66

64

62

4.6

58

56

54 1

1.05

1.1

1.15

1.2

P g3

1.25

Fig. 7. A trajectory for which 0 < tIP < T.

1.3

1.35

1.4

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

193

4.8

78

76

74

72

4.7

68

66 1.88

1.9

1.92

1.94

1.96

P g3

1.98

2

2.02

2.04

2.06

Fig. 8. A trajectory for which tIP = T.

below will, together with the local feasibility constraints, serve as necessary conditions for optimality.

Using the two constraints on the trajectory at tIP , we can compute derivatives of r, with respect to u. Local optimality is then guaranteed whenever (15b)–(15d) hold and

4.2.1. tIP = 0 If tIP = 0 any corrective control leading to a feasible path will by default give an optimal solution to Problem 1.

dr ≤ 0, dui

i ∈ Kl

(16a)

dr ≥ 0, dui

i ∈ Ku

(16b)

dr = 0, dui

i ∈ Kf

(16c)

4.2.2. 0 < tIP < T When 0 < tIP < T there are two important aspects of the trajectory, the position s(tIP , 0 ) and the derivative st (tIP , 0 ), at time tIP . These can be independently controlled when no limits on the controls are reached. We denote by uA the average control in [0, tIP ] and let v be the control action taken at time tIP . We define



A

(tIP , u , r) =

so that t (tIP , v, r) =

yp + rdy + tIP uA Pc + ˇ(tIP ) [P0 − Pc ]



v ˇ (tIP ) [P0 − Pc ]



.

 .

Assume now that uA and r* are the average control and the scalar multiple corresponding to an optimal path, respectively, so that s(tIP , 0 ) = (tIP , uA , r* ). Let Kl = {i : uAi = ui }, Ku = {i : uAi = u¯ i }, Kf = {i : uAi ∈ (ui , u¯ i )} and let Kc = Ku ∪ Kl . When i ∈ Kc we must have

vi = uAi by right-continuity of u. If, on the other hand, i ∈ Kf we A must have that n post ui (tIP , u , r) = 0, otherwise just a small change in ui in the right direction would mean that we lose contact with post thus violating Proposition 1. Hence, in the optimum, v is either decided by uA or does not affect (12). However, it may still affect (13). To see if we have a feasible solution we thus pick vi , i ∈ Kf , to minimize the left hand side of (13). To define local optimality we can thus drop v and end up with:

max

r

tIP ∈(0,T ), r∈R, u∈Rk

(15a)

The derivatives of r with respect to u can be obtained by eliminating tIP using (15b) and (15c). The necessary optimality condition given above is slightly impractical as a Newton based solution method would require us to compute second order derivatives of r with respect to u. To find a more applicable set of necessary conditions we define the perturbation cone U := {u ∈ Rm : ∃ε > 0, ui ≤ uAi + εui ≤ u¯ i }, so that U is the cone containing all vectors u such that we can make a small change of uA in the direction u without violating the restrictions on the control. We also define the feasibility cone Du : = {u ∈ Rm : [u 0]npost ≤ 0}. Hence, Du is the first order approximation of the restriction of the post-contingency stable domain to the u-space. The necessary condition for optimality that no possible (small) change in the control should make the trajectory lose contact with the stability boundary can be reformulated as U∩ int(Du ) = ∅. This holds true if and only if (npost )i = 0,

∀i ∈ Kf ,

(17a)

(npost )i ≥ 0,

∀i ∈ Kl ,

(17b)

(npost )i ≤ 0,

∀i ∈ K .

(17c)

u

Hence, instead of trying to find a vector (t, u, r) that solves (15b)–(15d) and (16), we may search along post for a point such that (15c), (15d) and (17) holds using, for example, a continuation method [19].

n post t (tIP , u, r) = 0,

(15c)

4.2.3. tIP = T When tIP = T, the only limiting factor is the post-contingency point when the load has recovered at time T. We thus get

¯ u ≤ u ≤ u,

(15d)

r∈R,u∈Rk

s.t.

(tIP , u, r) ∈ post ,

(15b)

max

r

(18a)

194

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

s.t.

p + rds +

Tu



0

∈ post

¯ u ≤ u ≤ u,

(18b) (18c)

Since tIP is fixed to T in this case, we only need the one available constraint to relate changes in u to changes in r. The local optimality constraints for this case is thus (18b), (18c) and (16) or (17). 4.3. Intersection at CPs of post

which together with (19b)–(19f) forms a necessary condition for optimality in this case. If N = nf , then the free control parameters are all decided by constraint (19b). If the CP is optimal then there is a j such that f,j f,k npost · npost < 0, for k = 1, . . ., j − 1, j + 1, . . ., N . 4.3.2. tIP = T In analogy to the previous extension, we get, for tIP = T, that max

r

(21a)

r∈R,u∈Rk

As noted in [7], when controlling the system to obtain maximal loadability we often end up in CPs of the stability surface. In Assumption 1 we have assumed that the intersection point is located away from a CP. We now adapt to the case where sIP is located at a CP of post . For tIP = 0 CPs of post will have no effect since we only have to find a feasible path. Hence, we can omit the case with tIP = 0. Assume now that the intersection point is located at a CP  where N ≥ 2 smooth manifolds, 1post , . . ., N post , intersect nontangentially. The smooth manifold consisting of such CPs will thus be of codimension N . 4.3.1. 0 < tIP < T This time we need to include v, since the ith component of the normal is not necessarily zero even when i ∈ Kf , and get: max tIP ∈ (0, T ), r ∈ R, u, v ∈ s.t.

r (19a)

Rk

(tIP , u, r) ∈ ∩ kpost ,

(19b)

k=1



(nkpost ) t (tIP , v, r) = 0, ¯ u ≤ u ≤ u,

k = 1, . . ., N

(19c) (19d)

∀i ∈ Kf ,

ui ≤ vi ≤ u¯ i ,

f,k

⎢ ⎢ ⎣

npost = ⎢

∀i ∈ Kc ,

(nkpost )

N

∈ ∩ kpost

(21b)

k=1

¯ u ≤ u ≤ u,

(21c)

The analysis of the previous case specializes to this case by replacing (19b)–(19f) with (21b) and (21c). 4.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 at its capacity limit cannot ramp up in order to save the system from collapse. To include the constraint yi ≤ (yp + rdy + tu)i ≤ y¯ i ,

(22)

5. Prediction of globally optimal corrective controls To predict the shape of the corrective control in the optimum we use the second order approximations of the stability boundary developed and tested in [20–23]. 5.1. Approximating the stability boundary As mentioned above the stability boundary is smooth at any point 0 ∈ L0 ⊂ , which is not a CP. Here, we define L0 to be the largest smooth subset of  to which 0 belongs. We denote by n = n (0 ) ∈ Sm−1 the normal vector to  at 0 . Here Sm−1 is the unit sphere in Rm . Given the normal vector n (0 ), a basis {c1 , . . ., cm−1 } for the tangent hyperplane T0  can be computed by, for example, Gram-Schmidt orthogonalization. Let C(= C(0 )) = [c1 . . . cm−1 ]. The second order Taylor expansion of  around the point 0 is given by H 0 : Rm−1 → Rm , H 0 (xc ) = 0 + Cxc +

1 II (xc )n (0 ), 2 0

(23)

where II0 is the second fundamental form of  at 0 ,



⎥ ⎥ ⎥, ⎦



II0 (xc ) = − dN 0 (xc ), xc .

in

f

cannot be linearly independent. Assume that there is a subset of the normal vectors f,jM 1 {nf,j post , . . ., npost }, with M ≤ nf that is not linearly independent and let this subset be chosen to give the minimal M. Then any (nf + 1 − M) × nf matrix H−j1 whose rows are orthonormal and f,jM 2 orthogonal to the vectors {nf,j post , . . ., npost } will have 1 H−j1 nf,j post = 0,

0



(19f)

⎤

i1

.. . (nkpost )

Tu

(19e)

where nkpost , for k = 1, . . ., N , are the normals to the smooth parts of post intersecting at the CP. Here, constraint (19b) specifies that the trajectory should intersect the set of CPs, constraints (19c) makes sure that the trajectory can follow the set of CPs and constraint (19f) guarantees that the trajectory can be taken to be smooth at the intersection point. Now, the feasibility cone changes to Du := {u ∈ Rm : [u 0]nkpost ≤ 0, k = 1, . . ., N }. At the optimal intersection point we must thus have U∩ int(Du ) = ∅. If we restrict our attention to the subset of free f,1 f,N , where controls Kf , this means that the vectors npost , . . ., npost

⎡

p + rds +

we add the hypersurfaces yi = yi and yi = y¯ i to post . N

vi = ui ,

s.t.

(20)

(24)

The map dN0 : T0  → T0  is the Weingarten map [18] of  at 0 . The Weingarten map is the differential of the Gauss map N :  → Sm−1 , i.e. the map that takes the point 0 ∈  to the normal vector n (0 ) ∈ Sm−1 . Let J = {1, . . ., M}, where M is the number of manifolds that j



make up the stability boundary. For each j ∈ J, let Hj = H 0 be the approximate parameterization of the jth smooth part of the staj bility boundary, j , with 0 the base point of the approximation of j . For the purpose of checking whether a parameter vector corresponds to a stable operating point we can evaluate the

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

distance to the approximation of  in the direction of the normal vector at 0 , dj ()

j 

= n (0 ) (Hj (C T ( − 0 )) − ) j 

j

= n (0 ) (0 − ) +

(25)

1 j II j (C T ( − 0 )). 2 0

To find a point of post around which to make approximations of the stability surface we solve Problem 1 for T = 0 and ˇ ≡ 0, i.e. we only require stability immediately following the contingency. If the limiting factor is post-contingency stability, the process will give such a point, c . If not, we take c to be the point of post closest to the post-contingency parameter vector obtained as the solution to Problem 1. Having found c we set pre to the value of the corresponding pre-contingency parameter vector. To predict the optimal corrective control and the corresponding maximal r we discretize time by  = {0, t, . . ., T − t, T} where t divides T with N = T/t. We then find an approximate optimal corrective control and a corresponding translation from pre along ds by solving

s.t.

˛

(26a)

dj (˜y(s, u, ˛), P˜ d (s, ˛)) ≥ 0,

∀s ∈  , j = 1, . . ., l,

(26b)

where initially l = 1, d1 is the distance to the approximation of post made at c ,



P˜ d (t, ˛) = Pc + ˇ(t)(P0 − Pc ) + ˛ (1 − ˇ(t))

dP c ds + ˇ(t)dL dpre



ˆ c (t) is the predicted optimal trajectory projected for t ∈  . Hence,  onto post . Now,

for t ∈  . This can be done by following the surface in a continuation manner [19,22]. Now, since tIP should be a minimum of h, we let tˆIP = arg minhˆ  (t). t∈ 

ˆ c (tˆIP ) belongs to a different smooth manifold than the ones we If  have approximated to compute u* in Section 5.2, then we add an ˆ c (tˆIP ) and include the new approximaapproximation of post at  tion when solving (26) (i.e. with l ← l + 1). ˆ c (tˆIP ) belongs to one of the The process is then repeated until  smooth parts of the stability boundary approximated by dj , for j ∈ {1, . . ., l}. When this is the case we expect hˆ  (tˆIP ) to be quite close to zero. We choose r so that sˆ (tˆIP , pre + rds ) ∈ post . We denote this first prediction of sIP with sˆIP . 6. Correction using local information Once a prediction of the optimal path has been computed we should correct the predicted path to obtain a continuous path that is an optimal solution to Problem 1. We also need to make sure that this prediction is in fact feasible. 6.1. Predicting sIP

and y˜ (t, u, ˛) = y0 +

c ∈post

ˆ (t, pre )), hˆ  (t) = n post (c (t))(c (t) − s

5.2. Predicting the optimal corrective control using the approximations of the stability boundary

max

to non-smoothness of post we may have to update this approximation by adding more approximated surfaces in the constraints, thus increasing l. We start by computing the distance from sˆ (tk , pre ) to post , which we denote hˆ  :  → R. To do this we define ˆ c (t) = arg min  sˆ (t, pre ) − c , 

Formulas for calculating the normal vector and the derivative of the Gauss map for SNB, SLL, HB and TL surfaces are collected in [22]. How to adapt the distance functions (25) to the case with tangential intersections of SNB-SLL surfaces is described in [23].

˛∈R,u( · )∈U

195

u(k)t + ˛dy .

To predict sIP we take all points corresponding to local minima of the function hˆ  . These points will then serve as initial guesses in the search for sIP .



6.2. Computing sIP

t/t  k=1

Here, pre = [y0 P0 ] is the initial guess of the pre-contingency loading and Pc = Pc (pre ) is the corresponding loading level immediately following the contingency. Now, if initially ˛* ≥ 0 the discretized trajectory

⎡

⎢ y0 +

sˆ (tk , pre ) = ⎣

k 

⎤

∗

u (j)t

j=1

⎥ ⎦

(27)

Pc + ˇ(tk )(P0 − Pc ) is optimal based on the approximation of post made at c . If ˛* < 0 we cannot find an approximately feasible trajectory corresponding to pre . Instead the solution predicts that we should move ˛* units in the direction ds from pre . We thus update by letting pre ← pre + ˛* ds , compute the corresponding post contingency loading level Pc and solve (26) again. By repeating this process until the magnitude of ˛* is sufficiently small we find a prediction of the optimal corrective control. 5.3. Updating the post-contingency stability boundary approximation Above, l was set to 1 to indicate that one second order surface approximation was used to approximate the optimal control. Due

Depending on the type of intersection point we have, we get different types of optimality criteria as explained in Section 4. Starting with the initial guesses, we solve the system of optimality equations using for example Newton’s method. 6.3. Checking feasibility Once we think that we have found the full set of intersection points we can make a new continuation along the corresponding path to validate that we have found a feasible corrective control. If the validation fails it will give us a new and hopefully improved hˆ  . 7. Numerical example In the numerical example we continue working with the test system from Section 2.3. We will consider a contingency situation where the impedance of the line connecting nodes 5 and 6 is doubled. Since the emphasis is on corrective controls rather than on stability in general we assume that the system is small-signal stable and thus only consider load-flow feasibility. The loadability surface corresponding to the contingency is plotted in Fig. 9. The green lines

196

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

Fig. 9. Maximal loadability in the post-contingency setting without corrective control and maximal loadability curves.

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. (For interpretation of the references to color in this text, the reader is referred to the web version of the article.) 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 , where p0 is a constant determined by the post-contingency loading (p0 = P0 /V02 , with V0 the post-contingency load-node voltage) and 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.

7.1. Unlimited production capacity Fig. 10 shows the result of solving Problem 1 starting with values of Pg2 and Pg3 ranging from 0 to 3 p.u., when the direction of parameter increase, ds , is that of only increasing the load. In all cases plotted in Fig. 10 there was only one intersection point. In Fig. 11 the corresponding intersection times are given. It can be noted that, when starting with y0 = [Pg2 Pg3 ] in the regions where post has a high slope in PLoad , the limiting factor seems to be feasibility immediately following the contingency, and for values of y0 close to where post reaches its maximal PLoad-value then tIP = T. Between these regions is a region where 0 < tIP < T. If neglecting that the path between t = 0 and t = T has to lie within the stable operation domain we would thus overestimate the

Fig. 10. The loadability surface with corrective control for the example system.

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

197

Fig. 11. The intersection times.

Fig. 12. The difference between to solving only for t = 0 and t = T.

loadability in this region. In Fig. 12 a plot of the error made when neglecting the path between times t = 0 and t = T is given. Although the magnitude of the error, for this case, is minor, there is a clear overestimation in a large region of the domain of the starting vector y0 . 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 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. 9) intersect the capacity limit hypersurface. Assume, therefore, that the output of generator 3 is limited with P¯ g3 = 1.5. In Fig. 13 the post-contingency stability boundary for this case is plotted along with the curves of maximal loadability. As can be seen in the figure the curve of maximal loadability when fixing Pg2 is located along the CP where the grey surface of points with Pg3 = P¯ g3 intersects the loadability surface. As a result, an optimal corrective control will attempt to drive the trajectory towards the corner point as this is the subset of the stability boundary that gives highest loadability.

Fig. 14 shows a sketch of the shape that Du takes along the set of CPs together with the restrictions of the normal vectors to the u-plane. The different situations in the figure can be described as follows: 1. At Point 1 the feasibility cone is pointing towards higher values of Pg2 . If this point is the intersection point for an optimal trajectory, i.e. if U∩ int(Du ) = ∅, we must have U = R × (−∞, 0]. Hence, uA2 = u¯ 2 which means that the ramping of generator 2 must be as large as possible. 2. In the second point the feasibility cone is pointing in the same direction but with a smaller spread. Hence, the same type of control as in the previous case is still optimal. 3. If we continue increasing Pg2 along the set of CPs the spread of the cone will decrease until eventually, at Point 3, Du has empty interior. This means that the set of CPs attains its maximal loadability here and thus any trajectory ending in this point will be optimal. Note that since ˇ is monotonically increasing, tIP = T for all feasible trajectories intersecting the stability boundary at this point. Point 3 corresponds the point in Fig. 13 where the curves of maximal loadability intersect. 4. Below Point 3, the direction of Du will flip and the cone is now pointing towards lower values of Pg2 . Hence, any optimal trajectory intersecting the set of CPs at Point 4 will have uA2 = u2 .

198

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

Fig. 13. Maximal loadability in the post-contingency state when P¯ g3 = 1.5 without corrective control, and maximal loadability curves. Here, the grey surface is the hyperplane with Pg3 = P¯ g3 .

5. At the fifth point the spread of the cone has increased to the level that one ray points away from the set of corner point. This corresponds the point in Fig. 13 where the curve of maximal loadability with fixed Pg2 leaves the set of CPs. 6. At Point 6 the vector pointing in the direction of decreasing values of Pg2 lies inside the feasibility cone. Hence, any optimal corrective control would try to steer away from the set of CPs here. Hence the only way that point 6 can be an intersection point is of tIP = 0. In Fig. 15 the loadability surface with corrective controls, CC is plotted. The loadability surface with corrective controls, CC , for the case of limited production capacity coincides with CC from the previous case for low values of Pg3 . However, as we reach the levels of Pg3 where the limits comes in play at the intersection point, the loadability will decrease as a consequence of the limited ability to control the trajectory.

In Fig. 16 the intersection times are plotted. As can be seen from this figure, in the vicinity of Point 3 in Fig. 14, i.e. the point of the set of CPs with maximal loadability, the intersection time equals T as anticipated. 8. Discussion 8.1. Optimality The necessary conditions for optimality stated in Section 4 guarantees that r will have an extremum when the intersection point satisfies the optimality condition. If the predicted optimal path is close enough to the actual optimal path we might expect that the solution we get from the correction step of Section 6 is in fact a global maximum. 8.2. Computational complexity An analysis of the computational time required for approximating the post-contingency stability boundary was made in [22]. Once the approximation is obtained, solving the arising quadratically constrained quadratic program (QCQP) in (26) is a convex problem if the post-contingency stable domain is convex. If this is not the case the problem becomes more challenging (see e.g. [24]). To compute hˆ  we need to search along the post-contingency stability boundary using, for example, a continuation method as proposed in [19]. This will be the most computationally demanding part and may make the method intractable for large systems. An alternative would be to do a regular continuation power flow [12] along the predicted optimal trajectory, sˆ . As this approach would not give us hˆ  we would have to predict the intersection point by observing a different index, such as the load bus voltage. The correction phase requires us to solve a system of equations and in the validation phase yet another continuation power flow is required. 8.3. Summary

Fig. 14. A schematic sketch of the feasibility cones along the CP surface.

In this article we investigate how we can reduce the infinite dimensional problem of proposing post-contingency corrective controls that maximize loadability in a given direction, to a finite

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

199

Fig. 15. The loadability surface with corrective control when P¯ g3 = 1.5.

Fig. 16. The intersection times when P¯ g3 = 1.5.

dimension problem. This is effectively done by introducing socalled intersection points, where the post contingency parameter trajectory intersects the stability boundary. Local feasibility and optimality conditions are posed at the intersection points. To solve the problem a global optimum is first approximated by discretizing time and approximating the stability boundary by polynomial expansions. This initial guess is then corrected using the local optimality constraints. In a numerical example on a small example system the results concur with the assumptions made.

References [1] T. Van Cutsem, An approach to corrective control of voltage instability using simulation and sensitivity, IEEE Trans. Power Syst. 10 (1995) 616–622. [2] T.J. Overbye, Computation of a practical method to restore power flow solvability, IEEE Trans. Power Syst. 10 (1995) 280–287. [3] A. Monticelli, M.V.F. Pereira, S. Granville, Security-constrained optimal power flow with post-contingency corrective rescheduling, IEEE Trans. Power Syst. 2 (1) (1987) 175–180. [4] F. Capitanescu, L. Wehenkel, Improving the statement of the corrective security-constrained optimal power-flow problem, IEEE Trans. Power Syst. 22 (2) (2007) 887–889. [5] F. Capitanescu, T. Van Cutsem, L. Wehenkel, Coupling optimization and dynamic simulation for preventive-corrective control of voltage instability, IEEE Trans. Power Syst. 24 (2) (2009) 796–805.

[6] M. Glavic, M. Hajian, W. Rosehart, T. Van Cutsem, Receding-horizon multi-step optimization to correct nonviable or unstable transmission voltages, IEEE Trans. Power Syst. 26 (3) (2011) 1641–1650. [7] 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. [8] T. Van Cutsem, C. Vournas, Voltage Stability of Electric Power Systems, Kluwer Academic Publishers, New York, 1998. [9] C.D. Vournas, M. Karystianos, N.G. Maratos, Bifurcation points and loadability limits as solutions of constrained optimization problems, in: Power Engineering Society Summer Meeting, 2000, IEEE, 2000. [10] Y. Kataoka, Y. Shinoda, Voltage stability limit of electric power systems with generator reactive power constraints considered, IEEE Trans. Power Syst. 20 (2005) 951–962. [11] F. Milano, Power System Modelling and Scripting, Springer Verlag, London, 2010. [12] 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. [13] C.A. Canizares, F.L. Alvarado, Point of collapse and continuation methods for large AC/DC systems, IEEE Trans. Power Syst. 8 (1) (1993) 1–8. [14] T. Van Cutsem, Costas D. Vournas, Voltage stability analysis in transient and mid-term time scales, IEEE Trans. Power Syst. 11 (1) (1996) 146–154. [15] C.D. Vournas, N. Sakellaridis, M. Karystianos, N.G. Maratos, Investigating power system stability limits, in: IEEE International Symposium on Circuits and Systems (ISCAS), IEEE, Island of Kos, Greece, 2006, pp. 726–729. [16] I. Dobson, Computing a closest bifurcation in multidimensional parameter space, J. Nonlinear Sci. 3 (1993) 307–327. [17] D.J. Hill, Nonlinear dynamic load models with recovery for voltage stability studies, IEEE Trans. Power Syst. 8 (1) (1993) 166–176. [18] L. Conlon, Differentiable Manifolds, 2nd ed., Birkhäuser, Boston, 2001.

200

M. Perninge / Electric Power Systems Research 116 (2014) 187–200

[19] I.A. Hiskens, R.J. Davy, Exploring the power flow solution space boundary, IEEE Trans. Power Syst. 16 (Aug 2001) 389–395. [20] 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. [21] 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.

[22] 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. [23] M. Perninge, Approximating the loadability surface in the presence of SNB-SLL corner points, Electr. Power Syst. Res. 96 (2013) 64–74. [24] S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, Cambridge, 2003.