Analysis of robust transient stability of power systems using sum of squares programming

Electrical Power and Energy Systems 115 (2020) 105401

Contents lists available at ScienceDirect

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

Analysis of robust transient stability of power systems using sum of squares programming

T

Shinsaku Izumia, , Hiroki Somekawab, Xin Xina, Taiga Yamasakia ⁎

a b

Faculty of Computer Science and Systems Engineering, Okayama Prefectural University, 111 Kuboki, Soja, Okayama 719-1197, Japan MITSUBISHI MOTORS CORPORATION, 1-1 Mizushimakaigandori, Kurashiki, Okayama 712-8501, Japan

ARTICLE INFO

ABSTRACT

Keywords: Power systems Transient stability Parameter uncertainties Sum of squares programming

This paper analyzes the transient stability of power systems with uncertain parameters. More precisely, we suppose that the inertia constants and damping coefficients of generators are uncertain, and consider the problem of checking the stability of an equilibrium state for all possible values of these uncertain parameters, and estimating the intersection of the region of attraction for those parameter values. By using a polytopic representation for the uncertain parameters and a concept from the field of robust control, we present a method for solving this problem. In the presented method, we solve a sum of squares programming problem with constraints imposed on systems where the values of the uncertain parameters correspond to the vertices of a polytope. We prove that, if the sum of squares programming problem is feasible, our method solves the analysis problem of robust transient stability. A numerical example demonstrates that our method correctly analyzes transient stability despite the parameter uncertainties, whereas an existing method gives an incorrect analysis result owing to the uncertainties.

1. Introduction The transient stability [1] of power systems refers to the ability of power systems to maintain their synchronism despite transient disturbances, e.g., faults, or loss of generators/loads. Analysis of transient stability is important in order to ensure highly reliable power supplies. In fact, it enables us to estimate whether a power supply will be maintained in the event of an accident. A key concept for analysis of transient stability is the region of attraction (ROA), i.e., the set of all system states converging to an equilibrium state. An example of the ROA is depicted in Fig. 1. If the state after the disturbances is in the ROA, it again converges to the equilibrium state; otherwise, it does not converge. Hence, a larger ROA implies higher transient stability. A typical method for investigating the ROA is the Lyapunov function approach in which a sublevel set of a Lyapunov function for a target system is used as an estimate of the ROA. This approach has been actively studied by a number of researchers. For instance, studies based on energy functions [2–4], decomposition into subsystems [5], an operational transform technique [6], extended Lyapunov functions [7], and numerical optimization [8] have been reported. Among these works, [8] has proposed a method to construct Lyapunov functions for power systems based on sum of squares (SOS) ⁎

programming [9,10]. In this method, a Lyapunov function is obtained by solving an SOS programming problem where the conditions to be satisfied by the Lyapunov function are described as constraints that polynomials must be SOS polynomials. This SOS-based method is promising for the following two reasons. First, the method can handle power systems with transfer conductances. Transfer conductances should be considered in practical situations because they correspond to losses in power systems. However, existing works have assumed that transfer conductances are zero [2,4] or sufficiently small [3,7]. By contrast, the SOS-based method can be used without restrictions on transfer conductances. Second, the SOS-based method achieves algorithmic construction of Lyapunov functions by using numerical optimization techniques. This reduces the time and effort needed to find Lyapunov functions. Recently, we have simplified the above method for practical use, and have provided a theoretical guarantee for the simplified method by showing the convergence of its iterative process under certain conditions [11]. Meanwhile, our previous work [11], as well as [8], has assumed that the system parameter values are completely known. Thus, we consider systems with uncertain parameters as the next research step. The motivation is as follows. First, relaxation of the above assumption is important to allow practical applications of the results of [8,11]. In fact, since power systems are generally large-scale complex

Corresponding author. E-mail addresses: [email protected] (S. Izumi), [email protected] (X. Xin), [email protected] (T. Yamasaki).

https://doi.org/10.1016/j.ijepes.2019.105401 Received 30 October 2018; Received in revised form 29 March 2019; Accepted 3 July 2019 0142-0615/ © 2019 Elsevier Ltd. All rights reserved.

Electrical Power and Energy Systems 115 (2020) 105401

S. Izumi, et al.

been proposed, respectively; however, those methods are not directly applicable to uncertain (non-polynomial) power systems. This paper is based on our conference paper [23], but differs from it in the following respects. First, this paper includes full explanations and rigorous proofs of the main results, omitted from [23]. Second, we discuss the choice of an initial Lyapunov function which is an important parameter that affects the behavior of the proposed algorithm. Finally, we present a numerical example to demonstrate the effectiveness of our method. Notation: Let , +, and 0 + be the real number field, the set of positive real numbers, and the set of nonnegative real numbers, respectively. We denote by 0 both the zero scalar and the zero vector. For the number c and the function f : n , a sublevel set of f is n | f (x ) c } . We use to expressed by c (f (x )) , i.e., c (f (x )) {x n represent the set of polynomials. Moreover, for x , let {p (x ) | p (0) = 0} , and let + be the set of positive definite 0 n\ {0}} . Finally, x polynomials, i.e., + {p (x ) 0 p (x ) > 0 deg(p) denotes the degree of the polynomial p.

Fig. 1. Region of attraction (ROA).

systems, accurate modeling of such systems is impractical. Furthermore, the system parameter values change depending on various factors such as temperature and humidity, and measuring all of them requires considerable time and effort. Second, although studies have been conducted on the topic of ROA estimation for systems with uncertain parameters [12,13], the methods proposed in those works are not directly applicable to power systems. This is because those methods are for systems described by polynomials, whereas power systems are described by non-polynomial functions. The purpose of the present paper is to extend the method we have proposed in [11] to power systems with uncertain parameters. As a first step toward achieving this purpose, we suppose that the inertia constants and damping coefficients of generators are uncertain; that is, their minimum and maximum values are known, but the actual values are unknown. Then, we consider the problem of checking the stability of an equilibrium state for all possible values of the uncertain parameters, and estimating the intersection of the ROA for those parameter values. To solve this problem, we introduce a polytopic representation [14], i.e., to represent the set of possible parameter values as a polytope whose vertices are specified by the (known) minimum and maximum values of the parameters. Based on this representation and a concept from the field of robust control, we present an SOS programming problem as an extension of the problem formulated in [8]. In our problem, constraints to obtain a Lyapunov function are imposed on systems where the values of the uncertain parameters correspond to the vertices of the polytope. We further present a method to solve this SOS programming problem by modifying the method proposed in [11]. Then, by utilizing the polytopic representation and properties of SOS polynomials, we prove that, if there exists a solution to our SOS programming problem, then it solves the analysis problem of transient stability stated above. As a result, we can conduct transient stability analysis taking into account all possible values of the inertia constants and damping coefficients of generators, merely by solving a single SOS programming problem. As final remarks in this section, we discuss related works. Some previous studies on stability analysis of power systems with uncertain parameters have been conducted [15–19]. In [15–17], methods for performing time-domain simulations that take into account the influence of uncertain parameters have been proposed. The authors of [18] have modeled an uncertain power system as a stochastic hybrid system, and have investigated its stochastic properties. However, ROA estimation has not been considered in those works. Although [19] has proposed a method for evaluating the stability and resiliency of power systems based on Lyapunov functions, that work has assumed that power injections are uncertain, whereas the inertia constants and damping coefficients of the generators are completely known. On the other hand, ROA estimation based on the Lyapunov function approach has been studied in [20–22]. The authors of [20] have proposed an ROA estimation method for power systems based on SOS programming by approximating trigonometric functions describing the system dynamics via the Taylor series expansion. However, they have not considered the parameter uncertainties. In [21,22], estimation methods using rational Lyapunov functions and piecewise polynomial Lyapunov functions have

2. Problem formulation 2.1. System description Consider a power system consisting of n generators. The swing equation describes the dynamics of generator i (i {1,2, …, n}) as

Mi ¨i (t ) = Pmi

Pei ( (t ))

(1)

Di i (t )

where i (t ) is the phase angle of the generator voltage, (t ) is the collective phase angle for all the generators, i.e., (t ) [ 1 (t ) 2 (t ) is n (t )] , Mi + is the inertia constant, Pmi the mechanical input, and Di + is the damping coefficient. The variable Pei ( (t )) is the electrical output defined as n

Pei ( (t ))

Ei Ej (Bijsin( i (t )

j (t ))

+ Cij cos( i (t )

j (t )))

j {1,2, … , n}

(2) where Ei are the suscep+ is the generator voltage and Bij,Cij tance and conductance between generators i and j, respectively. By considering i (t ) j (t ) in (2), let the state variable vector x (t ) 2n 1 . be x (t ) [ 1 (t ) n (t ) 2 (t ) n (t ) n 1 (t ) n (t ) 1 (t ) 2 (t ) n (t )] Then, (1) and (2) give the state equation of the system as

xn (t ) x2n 1 (t ) xn + 1 (t ) x2n 1 (t )

x1 (t ) x2 (t ) xn 1 (t ) x n (t ) xn + 1 (t ) x2n

1 (t )

=

x2n 2 (t ) x2n 1 (t ) (1/ M1)(Pm1 Pe1, x (x (t )) (1/ M2)(Pm2 Pe2, x (x (t ))

D1 xn (t )) D2 x n+ 1 (t ))

(1/ Mn)(Pmn

Dn x2n

Pen, x (x (t ))

1 (t ))

(3)

where x i (t ) (i {1,2, …, 2n 1}) is the i-th element of x (t ) and Pei, x (x (t )) (i {1,2, …, n}) is defined as

Pei, x (x (t )) Ei En (Bin sin x i (t ) + Cin cos x i (t )) + × (Bij sin(x i (t )

j {1,2, … , n 1}

xj (t )) + Cij cos(x i (t ) if i

En2 Cnn

j {1,2, … , n

E E (Bnj sin x j (t ) 1} n j

Ei Ej

xj (t ))) {1,2, …, n

1},

Cnj cos xj (t )) if i = n.

From (3), the equilibrium points of the system [x e1 x e2 x e (n 1) 0 0 0] where x ei (i = 1,2, …, n solution to

Pmi = Pei, x (x ) (i = 1,2, …, n 2

1).

(4)

take the form of 1) are given as a (5)

Electrical Power and Energy Systems 115 (2020) 105401

S. Izumi, et al.

Eq. (5) means that the electrical outputs of generators i (i = 1,2, …, n 1) are balanced with the mechanical inputs to them. Furthermore, Pmn is determined so that Pmn = Pen, x (x ) at the equilibrium point. To simplify the discussion, we assume that has an equilibrium point at x = 0 . Otherwise, we perform a coordinate transformation to shift the equilibrium point to x = 0 . 2.2. Analysis problem of robust transient stability We suppose that Mi and Di in (1) are uncertain. That is, Di [Di , Di+] Mi [Mi , Mi+] and hold for known Mi , Mi+, Di , Di+ , but the actual values of Mi and Di are unknown. + We further assume that the values of these uncertain parameters are time-invariant. Note that, since (5) and Pmn = Pen, x (x ) do not depend on Mi, Di (i = 1,2, …, n) from (4), the locations of the equilibrium points of the system are invariant for every Mi [Mi , Mi+], Di [Di , Di+] (i = 1,2, …, n) . Next, we provide examples to demonstrate that using the existing method given in [11] is risky when the system parameters are uncertain. Consider the power system with n 3 and uncertain M3 and D3. 0.0340 s2 , Pm1 0.900 pu, Based on [24], let M1 0.125 s2, M2 D 0.159 s, D 0.159 s , E1 1.05 pu, Pm2 1.00 pu, Pm3 1.29 pu , 1 2 E2 1.05 pu, E3 1.02 pu , and

B

C

2.99 1.51 1.23

1.51 2.72 1.01

1.23 1.01 , 2.37

0.845 0.287 0.210 0.287 0.420 0.213 0.210 0.213 0.277

3×3 where B,C are matrices whose (i, j) -elements are Bij and Cij in units of pu, respectively. The uncertain parameters M3 and D3 are set as M3 0.8M3*, M3+ 1.25M3*, M3* 0.0160 s2 , D3 0.25D3*, D3+ 4D3*, 0.159 s , where M3 ,D3 and D3 + are the nominal values of M3 and D3 , respectively. Then, has an equilibrium point at x = [ 0.288 0.176 0 0 0] , and this is asymptotically stable for M3 M3 and D3 D3 . We shift this equilibrium point to x = 0 through a coordinate transformation. For the shifted system, we estimate the ROA by using the existing method given in [11] assuming no uncertain parameters, where the D3 , revalues of M3 and D3 are supposed to be M3 M3 and D3 spectively. The resulting estimate is illustrated in Fig. 2 where x i = 0 for every i {3,4,5} , the solid line expresses the boundary of the estimate, D3 ). and the gray area expresses the true ROA (for M3 M3 and D3 We see that the resulting estimate is included in the true ROA. Moreover, Fig. 3 shows the time response of the system for x (0) [ 2.57 0.0400 0 0 0] indicated by the dot in Fig. 2, where each line corresponds to each element of x (t ) . It is apparent that x (t ) converges to the equilibrium point x = 0 . By contrast, Fig. 4 depicts the D3 , M3 and D3 D3 , i.e., M3 M3+ and D3 true ROA for M3 where x i = 0 for every i {3,4,5} and the solid line and dot are the same as those in Fig. 2. We see from (b) that the estimated ROA is not included in the true ROA. Furthermore, Fig. 5 depicts the time response of for x (0) indicated by the dot in Fig. 4 (i.e., the same x (0) as that in Fig. 3). It is apparent that, although x (0) is in the estimated ROA, x (t ) does not converge to x = 0 . These results imply that the existing method given in [11] can misestimate the ROA when the system parameters are uncertain. Operation of power systems based on such incorrect estimates can lead to accidents. To avoid this risk, we address the following problem.

Fig. 2. Estimate of ROA of

Fig. 3. Time response of in Fig. 2.

with M3

with M3

M3 and D3

M3 and D3

D3 .

D3 for x (0) indicated by dot

(b) if the equilibrium point x = 0 is asymptotically stable for all possible values of Mi, Di (i = 1,2, …, n) in (a), then estimate the intersection of 2n 1 such that the solution of starting the ROA (i.e., the set of all x from x converges to the equilibrium point) for every Mi [Mi , Mi+], Di [Di , Di+] (i = 1,2, …, n) . In Problem 1, (a) is for guaranteeing that the equilibrium point x = 0 is asymptotically stable for all possible values of the uncertain parameters Mi, Di (i = 1,2, …, n) . This is necessary in order to address (b). (b) is intended to obtain a region where the convergence of the system state is guaranteed regardless of the uncertainties of Mi, Di (i = 1,2, …, n) . If such a region is obtained, we can avoid misestimating the ROA. Here, we should note that the intersection of the ROA for all possible values of Mi, Di (i = 1,2, …, n) is always non-empty. This is because the locations of the equilibrium points of the system are the same for every Mi [Mi , Mi+], Di [Di , Di+] (i = 1,2, …, n) as explained at the beginning of this section. Meanwhile, it is difficult to solve Problem 1

Problem 1. For the power system , suppose that Mi, Di (i = 1,2, …, n) are uncertain. Then, (a) check the asymptotic stability of the equilibrium point x = 0 for every Mi [Mi , Mi+], Di [Di , Di+] (i = 1,2, …, n) ; 3

Electrical Power and Energy Systems 115 (2020) 105401

S. Izumi, et al.

3.1. Preliminary As a direct consequence of the Lyapunov stability theory (see, e.g., [25]), the following result is obtained. Lemma 1. Consider the power system with the equilibrium point x = 0 2n 1 and no uncertain parameters. Assume that there exist a set containing x = 0 and a continuously differentiable function V: 0+ such that

V (0) = 0 and V (x ) > 0

x

\ {0},

(6)

V (0) = 0 and V (x ) < 0

x

\ {0}.

(7)

Then, the equilibrium point x = 0 is asymptotically stable and the set is included in the ROA, where c + is c (V (x )) satisfying c (V (x )) a positive number. In this lemma, V (x ) is called the Lyapunov function. Lemma 1 provides a sufficient condition for the asymptotic stability of the equilibrium point x = 0 , and shows that c (V (x )) is an estimate of the ROA. Therefore, if there are no uncertain parameters, we can show the asymptotic stability of the equilibrium point and estimate the ROA by finding a pair (V (x ), c ) satisfying (6), (7), and c (V (x )) for a given . 3.2. Problem formulation based on sum of squares programming In [8], the problem of finding such a (V (x ), c ) has been formulated as an SOS programming problem. 3.2.1. Sum of squares programming problems First, the definition of the SOS is introduced. Fig. 4. Estimated ROA in Fig. 2 and true ROA for M3

M3+ and D3

m , the polynomial p (z ) Definition 1. For z is said to be the SOS if there exist polynomials qi (z ) (i = 1,2, …, µ ) satisfying

D3 .

µ

p (z ) = i=1

Fig. 5. Time response of dot in Fig. 4.

with M3

M3+ and D3

qi2 (z ).

(8)

Definition 1 means that SOS polynomials are polynomials that can be decomposed as the sum of the squares of polynomials. An example of an 9z12 + 2z22 + 6z1 z2 + 4z2 + 4 for z [z1 z2] . SOS polynomial is p (z ) 3z1 + z2 and q2 (z ) z2 + 2 satisfy (8). We denote by In fact, q1 (z ) the set of SOS polynomials. SOS programming problems are optimization problems with the constraints that polynomials must be SOS polynomials. This type of problem can be transformed into a semidefinite programming problem [26] under certain conditions, and can be efficiently solved by numerical optimization techniques.

D3 for x (0) indicated by

3.2.2. Coordinate transformation However, the SOS concept cannot be directly used for the system because (3) is not expressed by polynomials from the trigonometric h (x (t )) has functions in (4). Thus, the coordinate transformation z (t ) 3n 2 is the function such that been performed in [8], where h: 2n 1

through direct application of the method given in [11]. In fact, a straightforward approach to Problem 1 is to analyze the stability of the equilibrium point and estimate the ROA for each Mi [Mi , Mi+], Di [Di , Di+] (i = 1,2, …, n) according to [11]; however, this is impractical because the stability analysis and ROA estimation for all possible values of the uncertain parameters would require considerable time and effort.

z i (t ) 3. Transient stability analysis by sum of squares programming To solve Problem 1, we extend the method given in [11] to the case where Mi, Di (i = 1,2, …, n) are uncertain. Since [11] is based on [8], this section briefly introduces the result reported in [8], and then explains the method given in [11].

sin xN (t ) if N {1,2, …, n 1} s.t. i = 3N 2, 1 cos xN (t ) if N {1,2, …, n 1} s.t. i = 3N 1, x n 1 + N (t ) if N {1,2, …, n 1} s.t. i = 3N , x2n 1 (t ) if i = 3n 2

(9)

for the i-th element z i (t ) (i {1,2, …, 3n 2}) of z (t ) . Here, for simplicity of explanation, we exchange the order of the elements of z (t ) . More precisely, we take 4

Electrical Power and Energy Systems 115 (2020) 105401

S. Izumi, et al.

sin x i (t ) if i x i (t ) if i

z i (t )

{1,2, …, n 1}, {n,n + 1, …, 2n

1},

cos xi 2n + 1 (t ) if i {2n,2n + 1, …, 3n

1

2},

(10)

and use the resulting h. From (3), (4), and (10), the transformed system is given by

z (t ) = f (z (t )), g (z (t )) = 0

Fig. 6. Expanding interior.

(11)

3n 2 and g: where f : 3n 2 elements fi (z (t )) (i {1,2, …, 3n are defined as

are functions whose i-th 2}) and gi (z (t )) (i {1,2, …, n 1})

3n 2

n 1

z i + 2n 1 (t ))(z i + n 1 (t ) z2n 1 (t )) if i {1,2, …, n 1}, (1/Mi n + 1)(Pm (i n + 1) Pe (i n + 1), z (z (t )) if i {n,n + 1, …, 2n 1}, z i 2n + 1 (t )(z i n (t ) z 2n 1 (t )) if i {2n,2n + 1, …, 3n 2},

q+ (z ) 0 for q+ (z ) + and applying the Positivstellenstaz theorem (see, e.g., [10]) to the resulting constraints, we obtain the following SOS programming problem:

(1 fi (z (t ))

zi2 (t ) + zi2+ 2n

gi (z (t ))

1 (t )

2z i + 2n

(SOSP) Di

s.t. V (z )

n + 1 z i (t ))

1 (t )

q+ (z )

(15)

,

p+ (z ))

v2 (z ) g (z )

(V (z )

c)

(13)

s2 (z )(c

V (z ))

s3 (z ) V (z )

v3 (z ) g (z )

,

(16)

q+ (z )

(17)

v1 (z ), v2 (z ), v3 (z ) where are polynomial vectors, s1 (z ), s2 (z ), s3 (z ) are SOS polynomials, and is assumed to be given. The constraints (15)–(17) correspond to sufficient conditions for satisfying the three constraints of the OP. n 1

Ej En (Bjn zj (t ) + Cjn (1 × (Bjk (zj (t )(1 + Cjk ((1

zj + 2n

z k + 2n

z j + 2n 1 (t ))

1 (t ))(1

1 (t )))

+

k {1,2, … , n 1}

z k (t )(1

z k + 2n

zj + 2n

1 (t ))

k {1,2, … , n

+ Cnk (1

z k + 2n

Ej Ek

1 (t )))

3.3. Algorithm for solving SOSP

+ zj (t ) z k (t )))

if j En2 Cnn +

{1,2,…, n

1},

Since the SOSP has bilinear terms of the variables such as s2 (z ) V (z ) , it cannot be transformed into a semidefinite programming problem and cannot be efficiently solved. To address this issue, an iterative algorithm to solve the SOSP has been proposed in [8]; we have improved upon this algorithm, as reported in [11]. The improved algorithm is shown in Algorithm 1. Here, k + and + are the iteration indices, and superscripts are used to represent the variables at each iteration. For example, V k (z ) represents V (z ) at iteration k.

E E ( Bnk z k (t ) 1} n k

1 (t )))

if j = n.

(14)

The equality g (z (t )) = 0 in (11) is a constraint imposed as the property of the trigonometric functions, i.e., sin2 x i + cos2 x i = 1 for every x i . Eqs. (11)–(14) show that the transformed system is expressed by polynomials. It is remarked that, since z (t ) contains sin x i (t ) and cos x i (t ) (i = 1,2,…, n 1) from (10), a unique relationship between z (t ) < x i (t ) < (i = 1,2,…, n 1) by the fourand x (t ) is established for quadrant arctangent function. For this transformed system, if we find a pair (V (z ), c ) satisfying (6), , where x corresponds to z , by SOS program(7), and c (V (z )) ming, then from Lemma 1, the equilibrium point x = 0 of the original system (3) is asymptotically stable and c (V (h (x ))) is an estimate of the ROA. This is because V (h (x )) is a Lyapunov function for the system (3).

Algorithm 1 [11]: Transient stability analysis based on SOS programming. Step 0

0

Step 1 Step 2

3n

| V (z ) 2 | p (z ) +

Step 3

3n 2

| V (z )

c, g (z ) = 0, V (z )

0, z

c} =

p+

1

Set V (z ) problem: (SOSP1)

Vk

1 (z )

p

+

in Steps 4 and 5. Set

and

max

n 1, s , s , s 1 2 3

k 1,

and solve the SOS programming

c

Save the resulting c, s2 (z ) , and s3 (z ) as c k , s2k (z ) , and s3k (z ) , respectively.

s2k (z ) , and s3 (z ) Set c c k , s2 (z ) ming problem: (SOSP2) max V

Step 4

k

p+

1

k 1

Step 6

0} =

where p+ (z ) and c are assumed to be given. The first, second, and third (p+ (z )) constraints correspond to (6), c (V (z )) , and (7), respec0 with tively, where c (V (z )) corresponds to . By replacing z

and V (z ) as

k

and V k (z ) , respectively.

, then go to Step 5. Otherwise, set k

<

k + 1 and go to Step

2. If > 1 and the largest absolute value of the coefficients of

(z )

V 0 (z )

,

(17).

Save the resulting If

s3k (z ) , and solve the SOS program-

n 1, s 1

0, v1, v 2, v3

s.t. (15)

Step 5

0, g (z ) = 0, z 0} = , , g (z ) = 0, V (z ) c, V (z )

,

1.

(z ) .

Set p+ (z )

v2, v3

0

3n 2

1, and

0, k

s.t. (16) and (17).

max

V

Choose the initial Lyapunov function V 0 (z ) and polynomial p+0 (z ) , the polynomial q+ (z ) , and the thresholds

3.2.3. Problem formulation To find an appropriate (V (z ), c ) , [8] has employed the concept of expanding interior [9,10], which is explained as follows. For the positive definite polynomial p+ (z ) +, + and the positive number 3n 2 . If (p+ (z )) consider the set (p+ (z )) c (V (z )) , expanding (p+ (z )) yields the expansion of c (V (z )) , as shown in Fig. 6, which will improve the ROA estimate. Thus, for a given p+ (z ) , we find the pair (V (z ), c ) maximizing under (p+ (z )) c (V (z )) . Based on this idea, we consider the following optimization problem:

{z

v1 (z ) g (z )

s1 (z )(

Pej, z (z (t ))

s.t. {z {z

1, s , s , s 1 2 3

0, v1, v 2, v3

(12)

for

(OP)

maxn

V

Step 1.

p+

2

(z ) is smaller than

V k (z ), p+ (z )

V k (z ),

0

p,

then go to Step 6. Otherwise, set

c k, k

1, and

+ 1, and go to

Output V k (h (x )) and c k .

In Algorithm 1, we iteratively solve the two subproblems, i.e., the SOSP1 and the SOSP2, to obtain a solution to the SOSP. The SOSP1 is 5

Electrical Power and Energy Systems 115 (2020) 105401

S. Izumi, et al.

for determining c, s2 (z ) , and s3 (z ) given in the SOSP2, where V (z ) and are fixed. This problem can be solved by a linear search for c, i.e., a search for the optimal solution by increasing (or decreasing) the value 0.1,0.2, …. In fact, the bilinear terms do of c at a certain rate, e.g., c not appear in (16) and (17) when V (z ), , and c are fixed; thus, the problem of finding v2 (z ), v3 (z ), s1 (z ), s2 (z ) , and s3 (z ) satisfying (16) and (17) can be solved as a semidefinite feasibility problem. The SOSP2 is for finding a better , where c, s2 (z ) , and s3 (z ) are fixed. By a similar discussion to the above, we can show that this problem can be solved by a linear search for . In summary, we overcome the difficulty caused by the bilinear terms in the constraints by fixing some variables. Subsequently, in Step 5, p+ (z ) and 0 are updated based on V (z ) and c obtained as a solution to the SOSP. This is because setting (p+ (z )) c (V (z )) and solving the SOSP again will further improve the ROA estimate (see Fig. 6). We have proven the convergence of this algorithm in [11].

fi (z (t ))= (a /M + a +/ M +)((b + b +) × (Pm

+ a b + (1/M )(Pm

D +z i (t ))

(23)

a1+b1

an bn fp1 (z (t ))

a 2 b2

an bn fp2 (z (t ))

+ a1+b1+a2+b2+ an+bn+f p22n (z (t ))

(24)

Lemma 2 shows that, when Mi, Di (i = 1,2, …, n) are uncertain, the function f (z (t )) describing the system dynamics can be expressed by fpj (z (t )) (j = 1,2, …, 22n) , i.e., f (z (t )) for the combinations of known Mi , Mi+, Di , Di+ (i = 1,2, …, n) . Remark 1. Although one may think that Lemma 2 directly follows from the fact that f (z (t )) is linear in (1/ Mi ) and Di (i = 1,2, …, n) , this is not true. In fact, since the product of (1/ Mi ) and Di appears in f (z (t )) for each i {1,2, …, n} from (12), the numbers ai , ai+, bi , bi+ in (20) and (21) do not directly correspond to j in (19). In addition, Mi and Di are different in each element fi (z (t )) of f (z (t )); thus, we have to separately consider all elements of f (z (t )) and then represent the resulting elements as a vector. These facts render the proof difficult.

(18)

Based on Lemma 2, we extend the SOSP to

= 1 and

(SOSP )

V

max

0, v1, v 2, v3

n 1, s , s , s 1 2 3

s.t. (15), (16), and

j f pj (z (t ))

s2 (z )(c

(19)

j=1

(20) (21)

Di = bi Di + bi+Di+ for every i {1,2, …, n} , where ai , numbers satisfying ai + ai+ = 1 and and (21) for fi (z (t )) (i {n,n + 1,…, 2n

bi , bi+ bi + bi+ =

0+

s2 (z )(

fi (z (t ))=(a / M + a +/ M +) +

+

(b D + b D ) z i (t ))

q+ (z )

j

{1,2, …, 22n}

(25)

(SOSP0 ) find (V 0, v1, v3, s2) s.t. V 0 (z ) v1 (z ) g (z ) q+ (z ) ,

are nonnegative 1. By substituting (20) 1}) in (12), we obtain

ai+,

s3 (z ) Vj (z )

where Vj (z ) ( V (z )/ z ) fpj (z ) . The extended SOSP, i.e., the SOSP , is constructed by replacing the constraint (17) of the SOSP with those for fpj (z ) (j = 1,2, …, 22n) . As a method for solving Problem 1, we solve this SOSP . To obtain a solution to the SOSP , we use Algorithm 1 where (17) is replaced by (25). Here, the initial Lyapunov function V 0 (z ) in Step 0 is determined by solving the following problem:

Proof. Similar to (18), we have

1 1 1 = ai + ai+ + , Mi Mi Mi

V (z ))

v3 (z ) g (z )

where fpj (z (t )) is f (z (t )) where Mi, Di (i = 1,2, …, n) correspond to each vertex of a polytope, specified by each combination of Mi , Mi+, Di , Di+ (i = 1,2, …, n) .

i

D z i (t ))

, z (z (t ))

fi (z (t )) (i {n,n + 1,…, 2n 1}) , where for each (a + a +)(b + b +) = 1 is used. The multiplier n + + ( a + a )( b + b ) is equal to one, and the sum of the i i i i i=1 coefficients in (23) is also equal to one as described above. This implies that the sum of the coefficients in (24) is equal to one. In addition, according to the definitions of ai , ai+, bi , bi+ (i = 1,2, …, n) , the coefficients in (24) are nonnegative. Thus, by letting the coefficients in (24) be j (j = 1,2, …, 22n) , we obtain (19), which completes the proof. □

22n

where

, z (z (t ))

Pe

+

Lemma 2. For f (z (t )) given by (12), suppose that Mi, Di (i = 1,2, …, n) are uncertain. Then, for each Mi [Mi , Mi+], Di [Di , Di+] (i = 1,2, …, n) , 2n there exist nonnegative numbers j 0 + (j = 1,2, …, 2 ) such that

, z (z (t ))

Pe

+ a +b + (1/ M +)(Pm

+

+ for a1 , a1+ 0 + satisfying a1 + a1 = 1. This representation yields the following result.

Pe

+ a b (1/M )(Pm

where the first equality follows from b + b = 1. Note that the sum of the coefficients a b , a b +, a +b , and a +b + in (23) is equal to one from a + a + = 1 and b + b + = 1. We substitute (23) for the corresponding element of the vector f (z (t )) , and then multiply f (z (t )) n n (ai + ai+)(bi + bi+) = 1 is (ai + ai+)(bi + bi+) , where by i=1 i=1 noted. By expanding the resulting vector, we obtain

Our idea for solving Problem 1 is to introduce a polytopic representation, i.e., to represent the set of possible values of the uncertain parameters as a polytope whose vertices are specified by the (known) minimum and maximum values of the parameters. For example, from M1 [M1 , M1+], we can express 1/ M1 in (3) as

× (Pm

D +z i (t ))

, z (z (t ))

+

+

4.1. Proposed method

f (z (t )) =

, z (z (t ))

Pe

+

Now, we present a method for solving Problem 1. As explained in Section 2.2, a straightforward approach to Problem 1 is to use the existing method introduced in the previous section for each Mi [Mi , Mi+], Di [Di , Di+] (i = 1,2, …, n) . By contrast, the method presented here determines the stability of the equilibrium point for all possible values of the uncertain parameters at one time, and directly gives a single inner approximation of the intersection of the ROA for those parameter values.

j

Pe

D z i (t )) + b + (Pm Pe , z (z (t )) D +z i (t ))) = a b (1/ M )(Pm Pe , z (z (t )) D z i (t ))

f (z (t )) = a1 b1 a2 b2

1 1 1 = a1 + a1+ + M1 M1 M1

(b D + b +D +) z i (t ))

, z (z (t ))) + +

= (a /M + a / M )(b (Pm

4. Extension to case of uncertain parameters

22n j =1

Pe

(22)

r+ (z ))

v3 (z ) g (z ) where

n + 1. Eq. (22) yields 6

q+ (z ) +

(26)

0

V j (z ) j

{1,2, …, 22n}

is a positive number, r+ (z )

(27) +

is a positive definite

Electrical Power and Energy Systems 115 (2020) 105401

S. Izumi, et al.

polynomial, and q+ (z ), , and r+ (z ) are assumed to be given. The SOSP0 is an extension of the problem that [11] has addressed to determine V 0 (z ) . More precisely, the constraint of the problem in [11] is imposed on fpj (z ) (j = 1,2, …, 22n) , as shown in (27). Since the SOSP0 is an SOS feasibility problem and does not include the bilinear terms of the variables, this problem can be efficiently solved by numerical optimization techniques. Moreover, when the resulting V 0 (z ) is used, the convergence of Algorithm 1 is guaranteed; that is, k and p+ (z ) converge to a certain value and polynomial, respectively. In fact, by noting that the variables s2 (z ), V (z ), s3 (z ) , and v3 (z ) in (25) are common for every j {1,2, …, 22n} and that (25) takes the same form as that of (17), we can prove the convergence in a similar way to that in [11]. Note here that the optimality of the obtained solution is not guaranteed because we iteratively solve the SOSP by fixing some variables as apparent from Algorithm 1. Here, we comment on the choice of and r+ (z ) in the SOSP0 . By noting that and r+ (z ) correspond to c and V (z ) in (25), respectively, and that (25) is based on the third constraint of the OP, it is apparent 0 that (r+ (z )) represents the region where V j (z ) is negative definite for 2 n every j {1,2, …, 2 } . Hence, to obtain a better estimate of the ROA, we should choose and r+ (z ) so that (r+ (z )) becomes larger, under the constraint that the SOSP0 is feasible.

Proof. Since f (z (t )) for each Mi [Mi , Mi+], Di [Di , Di+] (i = 1,2, …, n) is described by (19), we prove the theorem by showing that a solution to the SOSP satisfies the constraints (15)–(17) of the 22n

2n = 1. SOSP for every j 0 + (j = 1,2, …, 2 ) satisfying j =1 j The constraints (15) and (16) are common to the SOSP and SOSP ; thus, we consider only (17). Since j (j = 1,2, …, 22n) are nonnegative, (25) and (S1) imply that the polynomial j(

2 (z )

V (z ))

v3 (z ) g (z )

2n

2 j =1

This, together with

q+ (z ))

22n

s3 (z )

j

j Vj (z )

.

j =1

j

= 1 and the definition of Vj (z ), shows 22n

s2 (z )(c

V (z ))

v3 (z ) g (z )

q+ (z )

s3 (z )

j j=1

V (z ) f (z ) z pj (30)

which yields

s2 (z )(c

V (z ))

v3 (z ) g (z )

q+ (z )

s3 (z )

V (z ) z

22n j f pj (z ) j=1

(31) because j (j = are scalar. It follows from (19) that (31) is the same as (17). Hence, if (25) holds, (17) holds for every j 0+

1,2, …, 22n)

(j = 1,2, …, 22n) satisfying

22n j =1

j

= 1. This completes the proof. □

Theorem 1 shows that, if the SOSP is feasible, its solution solves Problem 1. 4.3. Example Consider again the power system discussed in Section 2.2. Similar 0.176 0 0 0] , to Section 2.2, the equilibrium point x = [ 0.288 D3 , is shifted to which is asymptotically stable for M3 M3 and D3 x = 0 through a coordinate transformation. For this system, we solve the SOSP by using Algorithm 1, where the SOS problems are handled by the free MATLAB toolboxes: SOSTOOLS [29] version 3.00 and SeDuMi [30] version 1.3. As the inputs to SOSTOOLS, we select the degrees of the polynomials as deg(V ) 2, deg(v1) 2, deg(v2) 0, deg(v3) 2 , deg(s1) 0, deg(s2) 2 , 7 10 3 i = 1 zi2 , 3, and and deg(s3) 0 . Solving the SOSP0 for q+ (z )

and the nonnegative number

+

(28)

q+ (z ))

(29)

Next, we prove that our method provides a solution to Problem 1. To this end, we focus on the following two properties of the SOS.

1 (z )

v3 (z ) g (z )

j=1

4.2. Main result

,

s3 (z ) Vj (z )

22n

( s2 (z )(c

Remark 3. The computational complexity of the proposed method increases as the numbers of the generators and uncertain parameters increase. More precisely, the size of the semidefinite programming problem to be solved depends on 3n 2 and 2 where is the number of the uncertain parameters, because these express the dimension of z (t ) and the number of the polytope vertices, respectively. Hence, it is difficult to apply the proposed method to large-scale systems with many generators and uncertain parameters. Promising approaches for overcomming this issue can be found in [13,28]. The work [13] has proposed an approach to reduce computational complexity by focusing on only a few samples of the polytope vertices. Furthermore, [28] has proposed an approach to reduce the sizes of the systems used in computation by regarding the target system as interconnected (smallscale) subsystems.

2 (z )

V (z ))

is the SOS for every j {1,2, …, 22n} . Hence, by taking the summation of (28) over j {1,2, …, 22n} and using (S2), we obtain

Remark 2. The SOSP and the SOSP0 are based on the idea of finding a common Lyapunov for all vertices of the polytope; however, this concept is not new (see, e.g., [13,27]). Our contribution here is to extend the SOSP proposed in [8] to the SOSP presented herein, and to present a method to solve the SOSP by introducing the SOSP0 and combining it with Algorithm 1 in [11].

(S1) For the SOS polynomial (z ) (z ) is the SOS. 0+, (S2) For the SOS polynomials 1 (z ),

s2 (z )(c

is the SOS.

r+ (z )

These properties directly follow from Definition 1. (S1) and (S2) show that the SOS is preserved under nonnegative scalar multiples and sums, respectively. By utilizing (S1) and (S2), we obtain the following result.

7 i=1

z i2 gives the initial Lyapunov function as

V 0 (z ) = 11.3z12

11.4z1 z2 + 0.287z1 z 3 + 0.143z1 z 4

0.192z1 z5 + 0.691z1 z 6 + 1.23z1 z7 + 11.2z22 0.142z2 z 3 + 0.237z2 z 4 + 0.0889z2 z5 + 0.703z2 z 6 0.0433z2 z7 + 0.627z32 + 0.105z 3 z 4 + 0.0328z 3 z 5 0.393z3 z 6

Theorem 1. For the power system , suppose that Mi, Di (i = 1,2, …, n) are uncertain. If there exists a solution to the SOSP , then the following statements hold.

0.909z 3 z7 + 0.148z 42 0.0607z 4 z 6 0.104z 5 z7 +

(a) The equilibrium point x = 0 is asymptotically stable for every Mi [Mi , Mi+], Di [Di , Di+] (i = 1,2, …, n) . (b) The resulting c (V (h (x ))) is included in the intersection of the ROA for every Mi [Mi , Mi+], Di [Di , Di+] (i = 1,2, …, n) .

0.000220z 4 z5

0.182z 4 z7 + 0.0811z52 11.1z 62 7

0.102z 5 z 6

10.0z 6 z7 + 11.0z 72.

z 2, 0.05, and Furthermore, let p+0 (z ) i=1 i Then, by solving the SOSP , we obtain 7

p

0.2 .

Electrical Power and Energy Systems 115 (2020) 105401

S. Izumi, et al.

method is not. This demonstrates that the proposed method correctly estimates the ROA despite the parameter uncertainties but the existing method does not. This difference directly affects power supply reliability. In fact, power system operation based on incorrect ROA estimates can cause accidents and disrupt power supplies. Remark 4. Similar to [8–10], the proposed method is based on the expanding interior concept explained in Section 3.2.3; however, the resulting estimate in Fig. 7(a) is conservative. The reason for this is that the proposed method gives an estimate of the ROA in consideration of all possible values of the uncertain parameters, whereas the methods reported in [8–10] consider specific values of the system parameters only. 5. Conclusion This paper has addressed an analysis problem pertaining to the transient stability of power systems where the inertia constants and damping coefficients of the generators are uncertain. By using a polytopic representation for these uncertain parameters, we have extended an existing method based on SOS programming. In our method, we solve an SOS programming problem where constraints to obtain a Lyapunov function are imposed on systems characterized by the vertices of a polytope. Then, by appropriately combining the polytopic representation and properties of SOS polynomials, we have proven that our method solves the analysis problem if the SOS programming problem is feasible. With these results, we can perform the transient stability analysis in consideration of all possible values of the uncertain parameters, merely by solving a single SOS programming problem. In the future, our method should be extended to cases where other parameters (e.g., susceptances and conductances between generators) are uncertain. Furthermore, to reduce the conservativeness of the resulting estimate, we should consider more advanced Lyapunov functions, e.g., rational Lyapunov functions [21] and piecewise polynomial Lyapunov functions [22]. On the other hand, we expect to extend our method to controller design for improving robust transient stability.

Fig. 7. Estimate of intersection of ROA for some values of M3 and D3 .

Acknowledgements

V (h (x )) = 1.16sin2 x1 0.412sin x1sin x2 + 1.20sin2 x2 + 1.15cos2 x1 0.432cos x1cos x2 + 1.20cos2 x2 0.181sin x1cos x1 0.199sin x1cos x2 0.163sin x2cos x1 0.157sin x2cos x2 + 0.0400x3sin x1 + 0.00373x3cos x1 0.0583x3sin x2 0.00322x3cos x2 + 0.00219x 4sin x1 0.00373x 4 cos x1 + 0.0710x 4sin x2 0.00923x 4cos x2 0.0546x5sin x1 0.0116x5cos x1 + 0.0204x5sin x2 + 0.00252x5cos x2 + 0.0386x 32 + 0.00346x3 x 4

Funding: This work was supported by JSPS KAKENHI [Grant Nos. 26420425, 16K18124]. References [1] Pavella M, Murthy PG. Transient Stability of Power Systems. Wiley; 1994. [2] Varaiya PP, Wu FF, Chen R-L. Direct methods for transient stability analysis of power systems: Recent results. Proc IEEE 1985;73(12):1703–15. [3] Chiang H-D, Chu C-C. Theoretical foundation of the BCU method for direct stability analysis of network-reduction power system. Models with small transfer conductances. IEEE Trans Circuits Syst I: Fundam Theory Appl 1995;42(5):252–65. [4] Kojima C, Susuki Y, Tsumura K, Hara S. Decomposition of energy function and hierarchical transient stability diagnosis for power networks. Proceedings of the 54th IEEE Conference on Decision and Control. 2015. p. 3266–71. [5] Pai MA, Vittal V. Multimachine stability analysis using vector Lyapunov functions with inertial-centre decomposition. Int J Electric Power Energy Syst 1983;5(3):139–44. [6] Rau VG, Krishnamoorthy N. Construction of improved Lyapunov functions and extension of transient stability regions for power systems. Int J Electric Power Energy Syst 1989;11(1):65–9. [7] Bretas NG, Alberto LFC. Lyapunov function for power systems with transfer conductances: Extension of the invariance principle. IEEE Trans Power Syst 2003;18(2):769–77. [8] Anghel M, Milano F, Papachristodoulou A. Algorithmic construction of Lyapunov functions for power system stability analysis. IEEE Trans Circuits Syst I Regul Pap 2013;60(9):2533–46. [9] Jarvis-Wloszek Z. Lyapunov based analysis and controller synthesis for polynomial systems using sum-of-squares optimization Ph.D. Thesis Berkeley: University of California; 2003. [10] Jarvis-Wloszek Z, Feeley R, Tan W, Sun K, Packard A. Control applications of sum of squares programming. Positive Polynomials in Control. Springer; 2005. p. 3–22. [11] Izumi S, Somekawa H, Xin X, Yamasaki T. Estimation of regions of attraction of power systems by using sum of squares programming. Electr Eng

0.00724x3 x5

+ 0.0120x42 0.00645x 4 x5 + 0.0133x52 + 0.380sin x1 1.87cos x1 + 0.320sin x2 1.98cos x2 0.000504x3 + 0.0130x 4 + 0.00909x5 + 1.92 and c = 2.73. Thus, from Theorem 1, the equilibrium point x = 0 is asymptotically stable for every M3 [M3 , M3+] and D3 [D3 , D3+]. We further show the set c (V (h (x ))) in Fig. 7 where x i = 0 for every i {3,4,5} and the thick and thin lines respectively represent the boundaries of c (V (h (x ))) and an estimate resulting from the existing method given in [11] (which is the same as those in Figs. 2 and 4). The gray area represents the intersection of the true ROA for the pairs of the minimum and maximum values of M3 and D3 , i.e., (M3 , D3 ), (M3 , D3+), (M3+, D3 ) , and (M3+, D3+) , and for 20 pairs of values randomly chosen from the uniform probability distributions on [M3 , M3+] and [D3 , D3+]. We see that c (V (h (x ))) is included in the intersection of the true ROA, but the estimate resulting from the existing 8

Electrical Power and Energy Systems 115 (2020) 105401

S. Izumi, et al. 2018;100(4):2205–16. [12] Chesi G. Estimating the domain of attraction for uncertain polynomial systems. Automatica 2004;40(11):1981–6. [13] Topcu U, Packard AK, Seiler P, Balas GJ. Robust region-of-attraction estimation. IEEE Trans Autom Control 2010;55(1):137–42. [14] Boyd S, Ghaoui LE, Feron E, Balakrishnan V. Linear Matrix Inequalities in System and Control Theory. Soc Industr Appl Math (SIAM) 1994. [15] Wang S, Zheng Z, Wang C. Power system transient stability simulation under uncertainty based on interval method. Proceedings of the 2006 International Conference on Power System Technology. 2006. p. 437–42. [16] Dong ZY, Zhao JH, Hill DJ. Numerical simulation for stochastic transient stability assessment. IEEE Trans Power Syst 2012;27(4):1741–9. [17] Liao X, Liu K, Niu H, Luo J, Li Y, Qin L. An interval Taylor-based method for transient stability assessment of power systems with uncertainties. Int J Electric Power Energy Syst 2018;98:108–17. [18] Dhople SV, Chen YC, DeVille L, Domínguez-García AD. Analysis of power system dynamics subject to stochastic power injections. IEEE Trans Circuits Syst I Regul Pap 2013;60(12):3341–53. [19] Vu TL, Turitsyn K. A framework for robust assessment of power grid stability and resiliency. IEEE Trans Autom Control 2017;62(3):1165–77. [20] Mazumder SK, Fuente EP. Dynamic stability analysis of power network. Proceedings of the 2015 IEEE Energy Conversion Congress and Exposition. 2015. p.

5808–15. [21] Chesi G. Rational Lyapunov functions for estimating and controlling the robust domain of attraction. Automatica 2013;49(4):1051–7. [22] Chen Y-J, Tanaka M, Tanaka K, Wang HO. Stability analysis and region-of-attraction estimation using piecewise polynomial Lyapunov functions: Polynomial fuzzy model approach. IEEE Trans Fuzzy Syst 2015;23(4):1314–22. [23] Izumi S, Somekawa H, Xin X, Yamasaki T. Analysis of robust transient stability of power systems by using sum of squares programming. Proceedings of the 57th IEEE Conference on Decision and Control. 2018. p. 806–10. [24] Sauer PW, Pai MA. Power System Dynamics and Stability. Prentice Hall; 1998. [25] Khalil HK. Nonlinear Systems. Prentice Hall; 2002. [26] Boyd S, Vandenberghe L. Convex Optimization. Cambridge University Press; 2004. [27] Kalemba L, Uhlen K, Hovd M. Stability assessment of power systems based on a robust sum-of-squares optimization approach. Proceedings of the 2018 Power Systems Computation Conference. 2018. p. 1–8. [28] Kundu S, Anghel M. Stability and control of power systems using vector Lyapunov functions and sum-of-squares methods. Proceedings of the 2015 European Control Conference. 2015. p. 253–9. [29] Papachristodoulou A, Anderson J, Valmorbida G, Prajna S, Seiler P, Parrilo PA. SOSTOOLS: Sum of squares optimization toolbox for MATLAB. ; 2013. [30] Sturm JF. SeDuMi. ; 1998.

9