Development of bilinear power system representations for small signal stability analysis

Electric Power Systems Research 77 (2007) 1239–1248

Development of bilinear power system representations for small signal stability analysis J. Arroyo a , R. Betancourt b , A.R. Messina a,∗ , E. Barocio c a

Graduate Program in Electrical Engineering, Cinvestav, P.O. Box 31-438, Plaza La Luna, Guadalajara, Jal. 44550, Mexico b School of Electrical Engineering, Universidad de Colima, Manzanillo, Col. 28860, Mexico c School of Electrical and Mechanical Engineering, Universidad Aut´ onoma de Nuevo Le´on, Monterrey, NL. 66450, Mexico Received 4 April 2005; received in revised form 18 September 2006; accepted 19 September 2006 Available online 7 November 2006

Abstract In this paper, a new analysis technique for predicting and characterizing nonlinear behavior of stressed power networks is proposed. Making use of an analytical approximation of the system nonlinear model via the use of a truncated Carleman linearization technique, a bilinear state-space model of the power system is developed in which the second and higher order nonlinear terms are explicitly incorporated in the series expansion representation of the system model. The proposed framework enables the study of the dynamic behavior of nonlinear systems, both analytically and numerically, and can be used to represent a wide class of non-linear systems and oscillatory processes. Analytical criteria are developed based on the structural properties of the bilinear state-space model matrices, which predict the existence and stability character of modal interaction in terms of the eigenstructure of the linear system representation. The properties and behavior of the bilinear model are then investigated, and a number of useful results are derived. The present method is quite general and extends readily to higher-dimensional systems. A simplified 2-area, 4-machine system is used to illustrate the proposed procedure. Detailed nonlinear time-domain simulations are conducted to identify the strength of nonlinear behavior arising from interaction of the fundamental modes of oscillation, as well as to check the validity of the analysis. © 2006 Elsevier B.V. All rights reserved. Keywords: Bilinear dynamical systems; Perturbation theory; Small-signal stability

1. Introduction The analysis of nonlinear power system behavior has been the focus of considerable recent interest due to the possibility to detect a number of important phenomena such as nonlinear coupling and nonlinear interaction between the fundamental modes of oscillation in power systems [1–4]. An understanding of system non-linear characteristics is of both, theoretical and practical importance. This information may be used to design controllers as well as to predict various aspects of nonlinear behavior [5,6]. Over the past several years, interest in the study and characterization of system nonlinear behavior has grown enormously. This has led to the development of several analytical techniques

∗

Corresponding author. Tel.: +52 33 3134 5570; fax: +52 33 3134 5570. E-mail address: [email protected] (A.R. Messina).

0378-7796/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.epsr.2006.09.014

with the capacity to predict and characterize nonlinear processes that cannot be satisfactorily modeled under the assumption of linearity [5–7]. The calculation of such nonlinear behavior, however, is extremely involved and is computationally demanding. Recent studies have revealed that nonlinearities play an important role in determining crucial performance characteristics such as nonlinear coupling involving the electromechanical modes of oscillation in complex power systems [5,7,8]. Most existing approaches quantify nonlinear dynamic behavior using perturbation methods [1,3]. In these formulations, the nonlinear representation is first expanded into a sequence of homogeneous polynomials using a formal power series representation. The power series approximation inherent in perturbation theory can reproduce extremely well the quantitative characteristics of the system behavior but lacks information on the nature of mode interaction. Subsequently, the dynamic behavior of the system is analyzed using techniques from dynamical systems theory,

1240

J. Arroyo et al. / Electric Power Systems Research 77 (2007) 1239–1248

namely center manifold reduction and transformation to normal forms. While the underlying theory is straightforward, the computation of the normal form and the associated nonlinear transformations, however, may pose computational problems especially under heavy stress or heavy loading conditions. Such considerations have stimulated research into parallel methods which maintain the accuracy of nonlinear analysis while overcoming these limitations. In this paper, a bilinear state-space model of the power system is derived. The method improves the usual linearization of the system about an equilibrium point and allows taking advantage of existing linear system theory and software. Further, bilinear analysis methods can be used to understand changes in dynamic behavior and frequency-modulated processes. On the basis of this understanding, some fundamental characteristics of the nonlinear model are deduced rigorously using linear stability theory. Many of the results are of importance in the general field of nonlinear oscillations. The procedures developed are tested on a simplified 2-area, 4-machine test system. Numerical calculations are presented and discussed to check the performance of the method of analysis. The efficiency and accuracy of the method is demonstrated by comparisons to fully nonlinear computations by a commercially available transient stability software. The adoption of a quantitative description to nonlinear system behavior allows direct comparison of nonlinear results to linear analysis. Results generated by both analysis tools are seen to be consistent, and to agree well with the exact, step-by-step solution (SBSS).

2. Problem formulation and assumptions

second and third order, with coefficients  n  2  ∂ fm xi xj Dx2 fmo (x) = ∂xi ∂xj o i,j=1  n   ∂ 3 fm 3 Dx fmo (x) = xi xj xk ∂xi ∂xj ∂xk o i,j,k=1

.. . obtained from the power series expansion associated with f(x). The series approximation becomes precise as k increases; this, however, inevitably leads to considerable complexity and an explosion in computational cost. The stability of the above system can be conveniently studied using Carleman linearization. In order to derive the nonlinear Carleman representation for the system in (1) some background knowledge of this technique is necessary. 2.2. Properties of the Kronecker-product The Kronecker product is a matrix binary operator that maps two arbitrarily dimensioned matrices into a larger matrix with special block structure. Let A = (aij ) and B = (bij ) be arbitrary matrices of dimensions na × ma and nb × mb , respectively. Then their Kronecker product, denoted, A⊗B is an na mb × na mb block-matrix whose (i, j) block, is the p × q matrix aij B, defined as [10,11] ⎡ ⎤ a11 B a12 B · · · a1ma B ⎢ a B a B ··· a B ⎥ 22 2ma ⎢ 21 ⎥ ⎥ = [aij B] A⊗B=⎢ (3) . . .. ⎢ . ⎥ .. . . ⎣ . ⎦ . . ana 1 B

2.1. Basic concepts Consider a nonlinear set of coupled differential equations in a physical phase space of dimension n, x˙ = f(x, u)

(1)

where x is the vector of system states and f(x) is a real vector function; the vector u represents the control inputs. Further, it is assumed that f : n → n is continuous and sufficiently differentiable in the neighborhood of the origin. Expanding f(x, u), in a Taylor series about an equilibrium point (xsep ), and retaining terms up to order k, yields the perturbed model [9] x˙ ≈ Ax +

1 2 D f(xsep )(x − xsep )2 + · · · 2! x

1 + Dxk f(xsep )(x − xsep )k + O(x − xsep )k+1 + Bu k!

(2)

where A fx (x, u) = df/dx|(x=xsep ,u=uo ) is the Jacobian matrix of f(x) at xsep , and O(x − x)k+1 denotes an error of order k + 1; the p terms 1/p!Dx f(xsep ), p = 2, . . ., k represent variational effects of

am1 B

...

ana xma B

The following properties of the Kronecker product are easily proven from Eq. (3). (a) (b) (c) (d) (e)

(A + B)⊗C = (A⊗C) + (B⊗C) (A⊗B)T = AT ⊗BT (A⊗B)(C⊗D) = AC⊗BD An ⊗Bm = (An ⊗Im )(In ⊗Bm ) (A⊗B)−1 = A−1 ⊗B−1 if A−1 and B−1 exist.

In what follows we examine the use of these properties to determine structural properties of system behavior. The main motivation for using Kronecker products is that several properties of linear systems theory exist that allow analyzing the structure of the state representation as discussed in our simulations. 3. Mathematical formulation: bilinear system (BLS) representation 3.1. Nonlinear power system model Let the nonlinear system behavior be given by Eq. (1). A useful alternative representation to the series expansion model

J. Arroyo et al. / Electric Power Systems Research 77 (2007) 1239–1248

in (2) is obtained by introducing the variational matrices ⎡ 2 ∂ f1 ∂ 2 f1 ··· ⎢ ∂x2 ∂x1 ∂x2 ⎢ 1 ⎢ .. .. .. A2 = ⎢ . ⎢ . . ⎢ 2 ⎣ ∂ 2 fn ∂ fn ··· ∂xn ∂x1 ∂xn ∂x2 ⎡ 3 ∂ f1 ∂ 3 f1 ··· ⎢ ∂x3 ∂x12 ∂x2 ⎢ 1 ⎢ .. .. .. A3 = ⎢ . ⎢ . . ⎢ 3 ⎣ ∂ 3 fn ∂ fn ··· 2 ∂xn ∂x1 ∂xn2 ∂x2 .. .

the tensor product. Defining ∂2 f

1

∂x1 ∂xn .. . ∂2 f

n

∂xn2 ∂ 3 f1 ∂x12 ∂xn .. . ∂ 3 fn ∂xn3

x˙ = A1 x

+

r 

Further, noting that dx(1) dx(1) dx(2) = ⊗ x(1) + x(1) ⊗ dt dt dt is given by (6), we obtain

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

dx(2) = (A1 x(1) + A2 x(2) + A3 x(3) + · · · + Bu) ⊗ x(1) + x(1) dt ⊗(A1 x(1) + A2 x(2) + A3 x(3) + · · · + Bu)

⎤

dx(2) = [(A1 ⊗ I)(x(1) ⊗ x(1) ) + · · ·] + (B ⊗ I)x(1) u dt +[(I ⊗ A)(x(1) ⊗ x(1) ) + · · ·] + (I ⊗ B)x(1) u = [(A1 ⊗ I + I ⊗ A1 )x(2) + · · ·] + (B ⊗ I + I ⊗ B)x(1) u (8)

Ak x(k) + O(r + 1) + Bu

(4)

k=2

where r is the desired order of approximation. Here, A1 is the linear dynamic part of dimension n × n, A2 is the quadratic term of dimension n × n2 , and A3 is the n × n3 tensor; the term x(k) = x ⊗ · · · ⊗ x is the vector that includes the first-order Kroneckerk times

where I is the n × n identity matrix. Carrying out the operations indicated in Eq. (8) yields dx(2) = A21 x(2) + A22 x(2) + · · · + B20 x(1) u dt

(9)

where A21 = (A1 ⊗ I + I ⊗ A1 )

A22 = (A2 ⊗ I + I ⊗ A2 )

A23 = (A3 ⊗ I + I ⊗ A3 )

B20 = (B ⊗ I + I ⊗ B)

In a similar manner, it can be proved that

product vectors of kth order. In particular,

dx(1) d(x(2) ⊗ x(1) ) dx(2) dx(3) = = ⊗ x(1) + x(2) ⊗ dt dt dt dt

x(1) = [x1 x2 . . . xn ]T x(2) = x(1) ⊗ x(1) = [x12 x1 x2 . . . x1 xn . . . xn x1 xn x2 . . . xn2 ]

(7)

Application of the property (c) in Section 2.2 results in

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

and using the basic definitions of the tensor product, the nonlinear system behavior can be expressed in compact form as (1)

1241

T

x(3) = x(1) ⊗ x(1) ⊗ x(1) = [x13 . . . x12 xn x22 x1 . . . x22 xn . . . xn2 x1 . . . xn3 ] .. .

T

This approach has been used in the past to compute secondorder normal form approximations [12]. By extending this idea, a bilinear analysis approach is suggested in this paper for the analysis of more complex systems. The method computationally accounts for nonlinear interactions up to a specified order r, allowing to describe with more accuracy the dynamic behavior of the system.

(10)

From this expression, it follows that dx(3) = A31 x(3) + A32 x(4) + · · · + B30 x(2) u dt

(11)

and correspondingly A31 = (A1 ⊗ I ⊗ I + I ⊗ A1 ⊗ I + I ⊗ I ⊗ A1 ) A32 = (A2 ⊗ I ⊗ I + I ⊗ A2 ⊗ I + I ⊗ I ⊗ A2 ) A33 = (A3 ⊗ I ⊗ I + I ⊗ A3 ⊗ I + I ⊗ I ⊗ A3 )

(12)

B30 = (B ⊗ I ⊗ I + I ⊗ B ⊗ I + I ⊗ I ⊗ B) The procedure can be continued sequentially up to the desired order of approximation. At kth order, application of the above procedure yields

3.2. Bilinear state-space representation In order to determine a bilinear system representation, let now the new state vector be expressed in tensor notation as

dx(k) dx(k−1) dx(k−1) = ⊗x+x⊗ = Ak1 x(k) + Bk1 xu (13) dt dt dt

x⊗ = [ x(1)

where

x(2)

x(3)

···]

T

(5)

Using the properties of the Kronecker product, it can be readily proved that dx(1) = A1 x(1) + A2 x(2) + A3 x(3) + · · · + Bu dt

(6)

Ak1 = A11 ⊗ · · · ⊗ I + I ⊗ A11 ⊗ · · · ⊗ I + · · · + I ⊗ · · · ⊗ I ⊗ A11 Bk1 = A10 ⊗ · · · ⊗ I + I ⊗ B10 ⊗ · · · ⊗ I + · · · + I ⊗ · · · ⊗ I ⊗ B10

Based on the above derivation, an algorithm to compute higher order bilinear representations has been developed. A summary of the main steps of the analysis is as follows.

1242

J. Arroyo et al. / Electric Power Systems Research 77 (2007) 1239–1248

1. Determine high order derivatives as is the new plant matrix of dimension dx(k) = dt

k 

x ⊗ ··· ⊗

x˙ ⊗ · · ·x

k times

0

k  dx(k) x ⊗ · · · ⊗ f(x) ⊗ · · ·x = dt

C⊗ = [ C

i=1

=

k 

x ⊗ ··· ⊗

i=1

∞ 

A1j x(j) ⊗ x· · · ⊗ x

(15)

j=1

j=1

with

i = 2, 3, . . . , k;

p 

Nj xuj + Bu

(17)

j=1

y = Cx

I ⊗ · · · ⊗ I ⊗ A1j ⊗ I ⊗ · · ·I

j=1

A similar procedure is used to obtain Bk0 . 3. Derive the bilinear system representation up to the desired order of approximation k. We remark that the power series terms of order k generate components of order k and above, (k + 1, k + 2). A suitable choice of truncation order will depend upon the nature of the physical system and the objectives of the study; we have found that a third order bilinear model suffices for practical applications as discussed in the numerical evaluation of the method. Combining (6) and (9), (11) and (13), the BLS model can be expressed in the bilinear state-space form d ⊗ x = A⊗ x⊗ + N⊗ x⊗ u + B⊗ u dt

(16)

y = C⊗ x⊗ where A1

⎡ ⎤ B ⎢ ⎥ ⊗ B = ⎣ 0 ⎦, 0

0]

Aij = A1j ⊗ Ini−1 + In ⊗ Ai−1,j ,

x˙ = Ax +

Akj =

k=1

The above analytic procedure can be carried out up to any desired order of approximation to find higher-order nonlinear system representations in terms of the basic system eigenvalues. It is also interesting to note that Eq. (16) can be equivalently expressed in the form

∞

k 

B30 .. .

... 0 ... 0⎥ ⎥ ⎥ ⎥, ... 0⎥ ⎦ .. . 0

nk , and

j = 0, 1, . . . , k − i + 1

 dx(k) = Akj x(k+j−1) dt

⎢ ⎢ ⎢ ⎢ ⊗ A =⎢ ⎢ ⎢ ⎣

0

0 0

⎤

r 

with

or, equivalently

⎡

⎡

0 ⎢ B20 ⎢ ⎢ N⊗ = ⎢ .. ⎢ . ⎣

where x(k) = x ⊗ · · · ⊗ x. 2. Rewrite (14) as

nk ×

k=1

(14)

ith position

i=1

r 

A2 A21

A3 A22 A31

⎤

Ar ⎥ . . . A2r ⎥ ⎥ . . . A3r ⎥ ⎥ .. ⎥ .. ⎥ . . ⎦ Arr ...

where x ∈ n is the vector of system states, Nj ∈ n×n and uj denotes the jth input of the p-dimensional input vector u. As it is seen from this representation, bilinear systems are linear in control and linear in state, but not jointly linear in state and linear in control [11] and can be used to describe a wide variety of phenomena in the fields of engineering and physics. In the following sections, we derive several characteristics of nonlinear behavior based on the block-upper triangular characteristics of matrix A⊗ and provide a physical interpretation of the bilinear system model in terms the fundamental modes of oscillation of the system. 3.3. Analytical properties of the BLS An important issue in the analysis of small signal stability is that of determining the nature of system oscillations. In this section the basic properties and behavior of the bilinear model are investigated and a number of useful results are derived using the notion of a Kronecker product. For an n × n square matrix, the following general properties and definitions of the Kronecker operator can be written down. Definition 1. Let A be an arbitrary n × n matrix. Then the determinant of A, denoted by det A or |A|, is given by [10] det A = ai1 Ai1 + ai2 Ai2 + · · · + ain Ain =

n 

aik Aik

(18)

k=1

where Akj is the minor of matrix A obtained by deleting the kth and jth column.

J. Arroyo et al. / Electric Power Systems Research 77 (2007) 1239–1248

Definition 2. The eigenvalues of any arbitrary n × n matrix are the roots of its characteristic polynomial p(λ) = det(A − λI) = 0

(19)

Remark. With these definitions, we can express the eigenvalues of A⊗ in the form ⎧⎡ ⎤⎫ A2 A3 ⎪ ⎪ ⎬ ⎨ A1 − 1 ⎥ ⎢ 0 A21 − 2 A22 ⎦ = 0 det A⊗ = det ⎣ ⎪ ⎪ ⎭ ⎩ 0 0 A31 − 3 (20) where 1 , 2 and 3 are diagonal matrices of eigenvalues of appropriate dimensions. The nature of system behavior becomes evident from the analysis of system eigenvalues. Because matrix A⊗ is block-upper triangular, the system eigenvalues can be readily computed using some basic identities. Let det(z) denote the determinant of a matrix z. It is straightforward to show that the determinant of an upper triangular matrix is equal to the product of the diagonal elements. Therefore det A⊗ = det (A1 − 1 ) det(A21 − 2 ) det(A31 − 3 ) = 0

1243

arise in the output response of nonlinear systems. Some important remarks are in order in the interpretation of these results. The first observation that we can make is that second and higher order terms give rise to modal components with frequencies formed by the sums and differences of those of the primary (linear) modes. In addition, it follows from the foregoing discussion that the nonlinear interaction between different linear modes occurs in the form of sum of second-order and third-order combinations of the linear eigenvalues, i.e. λr + λs and λr + λs + λk that are combinations of the input frequencies. Thus for instance, for electromechanical modes of the form, λ = α ± jω, this interaction gives raise to intermodulation components at frequencies given by ωr ± ωs , and ωr ± ωs ± ωk . Furthermore, as can be seen from the form of Eq. (23) the BLS exhibits eigenvalues with frequencies two and three times higher than the corresponding linear frequencies; this happens when s, r and k are equal. Clearly, the amplitude, and importance, of these components depends upon the level of stress in the system. Further, the dimension of the BLS model increases with the degree of the nonlinear approximation order; the resulting system model is extremely sparse and is amenable for computer analysis. These observations are substantially confirmed in our numerical tests.

(21) The eigenvalues of the BLS are related to a linear space of matrices generated from matrix A1 are equal of the state or plant matrix, and the matrices A21 and A31 contain eigenvalues associated with the higher-order approximation. Recalling that A21 = A1 ⊗I + I⊗A1 , and using (21), it can readily be verified from the properties of the Kronecker product that det(A21 − 2 ) = det((A1 ⊗ I + I ⊗ A1 ) − 2 ) = 0

(22)

Definition 3. The eigenvalues of (A⊗Im ) + (In ⊗B) are the m × n numbers λr + λs , r = 1, 2, . . ., m, s = 1, 2, . . ., n. This matrix is often called the Kronecker sum of A and B is defined as ⊕(A + B). K

More formally, the Kronecker sum can be defined as ⊕ Ak = K  k=1

k=1

Ink−1 ⊗ Ak ⊗ InK where Ik is the k × k identity matrix. This 1

k+1

expression enables to study higher-order representations. Definition 3 can be easily extended to compute the eigenvalues of matrix A31 . The details are omitted. After applying this procedure, the eigenvalues of the system in (20) can be expressed in the form ⎤ ⎡ ␭r ⎥ ⎢ ␭r + ␭s eig (A⊗ ) = ⎣ ⎦ r = 1, 2, . . . , n; ␭r + ␭s + ␭k s = 1, 2, . . . , n;

k = 1, 2, . . . , n

(23)

The key fact of importance to us here, is that Eq. (23) can be used to predict how harmonics, and intermodulation terms can

3.4. Time domain closed form solutions Once the system is represented by the bilinear model (16), time domain solutions can be easily derived. Assume to this end, that the system (16) is taken to the Jordan form by a linear transformation x⊗ = U⊗ y⊗ : y˙ ⊗ = (V⊗ A⊗ U⊗ )y⊗ = ⊗ y⊗

(24)

where U⊗ is the matrix of right eigenvectors, V⊗ = (U⊗ )−l is the matrix of left eigenvectors, y⊗ is the new state vector in Jordan coordinates, and ⎤ ⎡ 1 ⎥ ⎢ 2 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⊗ 3  =⎢ ⎥ ⎥ ⎢ .. ⎥ ⎢ . ⎦ ⎣ k From linear system theory, the time evolution of the modal states is yj (t) = yj0 exp(λ⊗ j t),

j = 1, . . . , k −1

(25)

⊗ ⊗ ⊗ ⊗ with initial conditions y⊗ o = (U ) xo = (V )xo , where the ⊗ λj , j = 1, . . . , k are the elements of matrix ⊗ . At any instant t ≥ 0, Eq. (25) gives the modal system response in terms of the natural frequencies λ⊗ j , for any initial conditions yj0 , j = 1, . . . , n. Under proper considerations, the above approach can also be used to obtain the normal form representation of the system [12].

1244

J. Arroyo et al. / Electric Power Systems Research 77 (2007) 1239–1248

Fig. 1. One-line diagram of the test system.

In terms of the modal components, λ⊗ j , the approximate time evolution of the system states are xj (t) =

k 

⊗ ⊗ p⊗ ji xjo exp(λi t)

Under the above simplifications, the differential equations of motion for the classical system model are given by dδi = ωi ; dt

(26)

1 dωi = [Pmi − Pei − Di ωi ], dt 2Hi

i = 1, . . . , n (27)

i=1 ⊗ ⊗ where the p⊗ ji = uji vij represents the nonlinear participation of bilinear mode λ⊗ i on the jth state. Note that this definition of participation factors permits the calculation of both, linear and higher order measures of the interaction between states and modes. An important feature of these participation measures is that higher order effects are explicitly incorporated.

With n = 4, where δi is the angular position of the rotor, Pmi is the mechanical input power, Di is the generator damping coefficient and Hi is the inertia constant of the ith machine. Defining ⎡ ⎢ x⊗ = ⎣δ1 , . . . δn , ω1 , . . . ωn , δ1 δ1 , δ1 δ2 , . . . δn ωn , . . . ,      x(1)

x(2)

⎤T 4. Case study

⎥ δ1 δ1 δ1 , . . . δ2 δn ω1 , . . . ωn ωn ωn ⎦ .   

4.1. Modeling considerations

x(3)

The procedures developed are tested on a 2-area 4-machine test system taken from ref. [13] modified to consider different operating conditions. Fig. 1 gives a single line diagram of the study system. For the sake of mathematical clarity and simplicity, it is assumed that each generator is represented by the classical model; the loads L1 and L2 are modeled as constant impedances and the network is reduced to the generator internal nodes. At the base case condition, all four generators have similar loadings, and 200 MW is exported from Area 1 to Area 2. To stress the system, circuit #1 of the tie line (5–6) is taken out of service. Based on this scenario, several operating conditions are developed representing various levels of stress in the system as discussed in our numerical simulations. Appendix A gives the generator and system parameters, along with the base case operating condition used in the analysis.

and u = [Pm1 . . . Pmn ]T , the system can then be expressed in the form (16). Closed-form expressions for linear and higher-order approximations of system behavior are derived. In addition, three system representations were considered in the studies, namely: (a) a first-order system representation obtained by neglecting nonlinear terms in (6); (b) a second-order system representation obtained by expanding (9) up to order two; and (c) a third-order system representation obtained by expanding (16) up to order three. 4.2. Linear stability analysis Linear analysis studies were conducted to determine the fundamental modes of oscillation of the system. Tables 1 through 3 synthesize the linear system eigenvalues along with the states

Table 1 Linear system eigenvalues Mode

Eigenvalue (λr )

Damping (%)

Frequency (Hz)

Dominant state

1 2 3

−0.096 ± j2.870 −0.096 ± j8.555 −0.096 ± j8.612

3.348 1.124 1.116

0.450 1.362 1.370

␻3 ␻2 ␻3

Base case condition.

J. Arroyo et al. / Electric Power Systems Research 77 (2007) 1239–1248

1245

Table 2 System eigenvalues associated with matrix A21 Mode

Eigenvalue (λr + λs )

Damping (%)

Frequency (Hz)

Dominant state

4 5 6 7 8 9 10 11 12 13 14 15

−0.192 ± −0.192 ± −0.192 ± −0.192 ± −0.192 ± −0.192 ± −0.192 ± −0.192 ± −0.192 ± −0.192 −0.192 −0.192

1.124 1.120 1.116 1.683 1.675 3.381 3.348 3.348 95.83 100.0 100.0 100.0

2.723 2.732 2.741 1.818 1.827 0.905 0.914 0.914 0.009 0.000 0.000 0.000

␦ 2 ␻2 ␦ 2 ␻2 ␦ 3 ␻3 ␦ 3 ␻2 ␦ 3 ␻3 ␦ 4 ␻2 ␦ 4 ␻3 ␦ 3 ␻3 ␦ 4 ␻2 δ22 ␻ 3 ␻4 ␦ 3 ␻4

j17.11 j17.17 j 17.22 j11.42 j11.48 j5.684 j5.741 j5.741 j0.057

Base case condition.

with more contribution to the system modes for the base case condition. For the sake of clarity they are grouped according to the source of system modeling, namely, eigenvalues associated with the linear part A1 and eigenvalues associated with secondorder and third-order effects. Note that both, angle and speed terms are present and that the damping ratio of each modal component is the complex conjugate of the corresponding frequency term. The two-area system exhibits three electromechanical modes of interest (see Table 1): one inter-area mode representing the exchange of oscillating energy between Areas 1 and 2 (mode 1),

and two local modes (modes 2 and 3) associated with the local dynamics between generators in each area. From the results in Tables 2 and 3, we observe that the natural response of the system contains components of the linear matrix (fundamental modes) λ1 , . . ., λn in addition to harmonics of frequencies 2ωr , . . ., 3ωr , . . . Moreover, it is also apparent from these tables, that nonlinear system behavior gives rise to harmonics and intermodulation of the form ωr ± ωs . For instance, it can be seen that ω4 = 2ω2 , ω5 = (ω2 + ω3 ). Similarly, it can be seen from Table 3 that ω41 = 3ω1 , ω20 = (ω1 + 2ω2 ), etc. as suggested from the analysis of Eq. (23).

Table 3 System eigenvalues associated with matrix A31 Mode

Eigenvalue (λr + λs + λk )

Damping (%)

Frequency (Hz)

Dominant states

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

−0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288 −0.288

1.124 1.121 1.116 1.119 1.444 1.439 1.435 2.025 2.017 2.017 2.009 2.009 3.393 3.326 9.805 10.201 10.199 9.997 3.348 9.999 9.999 9.999 3.370 3.370 3.370 3.348 3.348 3.348

4.084 4.094 4.112 4.103 3.180 3.189 3.198 2.266 2.275 2.275 2.284 2.284 1.352 1.380 0.466 0.448 0.448 0.457 1.370 0.457 0.457 0.457 1.361 1.361 1.361 1.371 1.371 1.371

δ2 ω22 δ2 ω22 ␦3 ␻3 ␻4 ␦2 ␦3 ␻2 ␦2 ␦3 ␻2 ␦3 ␻2 ␻3 δ3 ω32 ␦2 ␦4 ␻2 ␦3 ␻2 ␻3 ␦4 ␻2 ␻3 δ23 ω3 δ4 ω32 δ4 ω22 ␦3 ␦4 ␻2 ␦3 ␻2 ␻3 ␦4 ␻2 ␻3 ␦4 ␻2 ␻3 δ4 ω32 ␦3 ␻3 ␻4 δ23 δ4 δ2 δ24 δ2 δ24 δ32 ␦3 ␦4 ␻2 δ23 ω2 δ32 δ23 δ4 δ33

Base case condition.

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

j25.664 j25.721 j25.836 j25.778 j19.980 j20.037 j20.094 j14.239 j14.295 j14.296 j14.353 j14.353 j8.497 j8.669 j2.928 j2.813 j2.814 j2.871 j8.611 j2.870 j2.870 j2.870 j8.554 j8.554 j8.554 j8.612 j8.612 j8.612

1246

J. Arroyo et al. / Electric Power Systems Research 77 (2007) 1239–1248

A deeper insight into the effect of higher order terms on system behavior is next obtained from the analysis of nonlinear participation factors. 4.3. Nonlinear participation factors To examine the significance of nonlinear effects in the system structure, and determine the states having a greater influence on a given mode of oscillation, we computed nonlinear participation factors using the procedure in Section 3.4. For this study, we consider three operation conditions with different power transfers from Area 1 to Area 2, namely: 200 MW (Base case condition, Case 1), 300 MW (Case 2), and 375 MW (Case 3). In each case, the system was stressed by increasing load L2. The load L1 was adjusted to match the desired level of power transference between the areas. Fig. 2 shows linear and higher order nonlinear participation factors for cases 1, 2 and 3. In the present context, participation factors provide the relative contribution of a give state, or state combination, to modes as explained above. From Fig. 2 we observe that nonlinear modes contribute significantly to the total response; strongly nonlinear dynamics takes places when the linear and nonlinear terms are of the same order of magnitude. Clearly, the coefficients of the nonlinear terms (second and third terms) are very significant compared with the coefficients of the linear terms. This is especially true for the highly stressed condition (Case 3) in Fig. 2(c). This also indicates that higher-order terms will have a significant effect on the output response of the system if the modes are excited by a given contingency. In order to further interpret these results, we computed nonlinear participation factors for the three operating conditions above. Of particular significance and practical importance, simulation results in Fig. 3 show that as stress is increased, the magnitude of the higher order terms relative to linear terms increases. This further motivates the need for including more detailed system representations in the analysis. The proposed procedure provides a precise characterization system behavior and permits the study of other system characteristics such as nonlinear controllability and observability.

Fig. 2. Comparison of participation factors for different modes: (a) Linear terms, (b) Second order terms, (c) Third order terms.

4.4. Model validation A systematic series of studies was finally undertaken in an effort to determine the validity and applicability and effectiveness of the proposed procedure. Following a similar approach to that in ref. [12], approximate time domain solutions were obtained from the bilinear representation in (16). The procedure consists of three main steps: 1. For a given disturbance, determine the post-disturbance stable equilibrium point xSEP . Move the origin to the postdisturbance SEP such that xo = xcl − xSEP , where xcl is the system condition at the end of the disturbance determined from a transient stability program.

Fig. 3. Comparison of participation factors for different operating conditions.

J. Arroyo et al. / Electric Power Systems Research 77 (2007) 1239–1248

1247

with the corresponding estimates for the linear and higher (second and third order) bilinear approximations. The results show that the third-order BLS representation is able to approximate well system behavior even for high stress conditions. The comparisons in all four cases are very satisfactory, both in amplitude and phase. Similar results are obtained for other machines. For low stress conditions (Case 1) simulation results suggest that the response is given mainly by the first order term; the agreement between the linear, second-order and third order approximation and the SBSS is good over the entire study period. As the initial loading is increased, second and thirdorder terms become significant as shown in Figs. 4(b) and (c). Higher-order approximations provide a more accurate approximation for the SBSS over longer times that the approximations given by linear analysis and enable the loss of synchronism to be verified.

5. Conclusion

Fig. 4. Comparison of relative rotor angle swings for different operating conditions: (a) Case study 1 - 200 MW, (b) Case study 2 - 300 MW, (c) Case study 3 - 375 MW.

2. Obtain the second and third order power system representation of the system about xSEP following the procedure in Section 3. 3. Compute approximate bilinear solutions using (26). In order to assess the effect of higher-order approximations on system behavior, three system representations of the system were considered in the studies, namely: (i) a linear approximation obtained by ignoring nonlinear terms in (6), (ii) a second order bilinear representation obtained by expanding (14) up to order 2, and (iii) a third order bilinear representation obtained from the proposed approach in this article (refer to Eq. (16)). For the sake of completeness, the solutions were compared with those obtained by exact step-by-step simulation (SBSS) using a commercial transient stability program. Fig. 4 shows the variation of the internal angle of generator 3 (relative to generator 1) for each case study along

In this paper, a new analysis technique in the study and characterization of the high-dimensional nonlinearity of power system models has been proposed. The proposed procedure introduced herein is thought to have potentially important applications for the dynamic analysis of nonlinear systems. At the same time, the method leads to better insight into the physical phenomena of nonlinear oscillations, such as mode interaction and results in models which can be considerably more accurate than those obtained from conventional methods. Study results indicate that bilinear models based on Carleman linearization can accurately characterize nonlinear behavior even for the most stringent operating conditions. The method could be used to detect operating conditions under which more advanced analytical techniques, such as normal form analyzes could be applied. Future work including improved analytical models and numerical validation under more general operating conditions is needed. Additional aspects, such as the analysis of nonlinear controllability and observability of the developed models are to be further investigated.

Appendix A. Data for the Case study 1 in Section 4 Tables A.1 through A.4 summarize the network and operating conditions used in the base case.

Table A.1 Generator data (p.u. on machine base) Generator

x˙ d

H

D

BMVA

1 2 3 4

0.3 0.3 0.3 0.3

6.5 6.5 6.5 6.5

2.5 2.5 2.5 2.5

900 900 900 900

1248

J. Arroyo et al. / Electric Power Systems Research 77 (2007) 1239–1248

Table A.2 Line data (p.u. on a 100 MVA base) Bus#

Bus#

R

X

1 2 5 5 6 4

2 5 6 (out) 6 4 3

0.0025 0.0010 0.0220 0.0220 0.0010 0.0025

0.025 0.010 0.220 0.220 0.010 0.025

Table A.3 Shunt susceptance (p.u. on a 100 MVA base) Bus#

Shunt susceptance

5 6

2.49 2.51

Table A.4 Base case operating condition (p.u. on a 100 MVA base) Bus#

Voltage

Real power

Reactive power

Real load

Reactive load

1 2 3 4 5 6

1.03 ∠ 26.0 1.01 ∠ 16.4 1.03 ∠ 0.00 1.01 ∠ −9.6 1.00 ∠ 8.55 0.99 ∠ −17.5

7.0 7.0 7.046 7.0 0 0

0.68 0.90 0.68 1.29 0 0

0 0 0 0 11.67 15.67

0 0 0 0 1.0 1.0

References [1] J. Thapar, V. Vittal, W. Kliemann, A.A. Fouad, Application of the normal form of vector fields to predict interarea separation in power systems, IEEE Trans. Power Syst. 12 (2) (1997) 844–850.

[2] A.R. Messina, V. Vittal, Assessment of nonlinear interaction between nonlinearly coupled modes using higher order spectra, IEEE Trans. Power Syst. 20 (1) (2005) 375–383. [3] A.R. Messina, E. Barocio, Assessment of nonlinear modal interaction in stressed power networks using the method of normal forms, Electrical Power & Energy Syst. 25 (1) (2003) 21–29. [4] V. Vittal, N. Bhatia, A.A. Fouad, Analysis of the inter-area mode phenomenon in power systems following large disturbances, IEEE Trans. Power Syst. 6 (4) (1991) 1515–1521. [5] J.J. Sanchez-Gasca et al., Inclusion of higher order terms for small-signal (modal) analysis: Committee report-Task force on assessing the need to include higher order terms for small-signal (modal) analysis, IEEE Trans. on Power Syst. 20 (4) 1886–1903. [6] L. Chih-Ming, V. Vittal, W. Kliemann, A.A. Fouad, Investigation of modal interaction and its effects on control performance in stressed power systems using normal forms of vector fields, IEEE Trans. Power Syst. 11 (2) (1996) 781–787. [7] S. Liu, A.R. Messina, V. Vittal, Assessing placement of controllers and nonlinear behavior using normal form analysis, IEEE Trans. on Power Syst. 20 (3) 1486–1495. [8] N. Pariz, H.M. Shanechi, E. Vaahedi, Explaining and validating stressed power systems behavior using modal series, IEEE Trans. Power Syst. 18 (2) (2003) 778–785. [9] C. Robinson, Dynamical systems—Stability, Symbolic Dynamics and Chaos, CRC Press, Boca Raton, Florida, 1999. [10] P. Lancaster, M. Tismenetsky, The Theory of Matrices: With Applications, Academic Press, San Diego, California, 1985. [11] W.J. Rugh, Nonlinear System Theory: The Volterra/Wiener Approach, The John Hopkins University Press, 1981. [12] A.R. Messina, J. Arroyo, E. Barocio, Analysis of modal interaction in power systems with FACTS controllers using normal forms, IEEE Power Engineering Society General Meeting, July 2003. [13] M. Klein, G.J. Rogers, P. Kundur, A fundamental study of inter-area oscillations in power systems, IEEE Trans. Power Syst. 6 (3) (1991) 914– 921.