- Email: [email protected]
An optimal modal approximation method for model reduction of linear power system models
Electrical Power and Energy Systems 44 (2013) 293–300
Contents lists available at SciVerse ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
An optimal modal approximation method for model reduction of linear power system models A. López Ríos, A.R. Messina ⇑ Department of Electrical Engineering, The Center for Research and Advanced Studies (Cinvestav), Av. del Bosque 1145, colonia el Bajío, Zapopan, CP 45019 Jalisco, Mexico
a r t i c l e
i n f o
Article history: Received 1 February 2011 Received in revised form 1 August 2012 Accepted 1 August 2012 Available online 26 September 2012 Keywords: Modal approximation Proper orthogonal decomposition Optimized modal approximation
a b s t r a c t A number of recent studies have examined the problem of modal reduction of large-scale power system linear models. In this paper, an optimal model reduction strategy for the analysis and control of inter-area oscillations in power systems that uses concepts from proper orthogonal decomposition and modal analysis and preserves exactly the input–output properties of the original system is presented. The technique uses a new projection matrix that minimizes the error between the states of the original system and the states of its reduced-order model, and is suitable for analyzing large power system models described by differential–algebraic equation (DAE) systems. Based on this theoretical framework, an efficient algorithm for model reduction is proposed which combines the ease of use of linear modal analysis, with the flexibility of extracting relevant behavior that statistical multivariate methods can offer. Methods for analysis of reduced system representations using linear approximations are critically reviewed and compared. The proposed technique is tested on a 50-machine, 145-bus, and 453-line test system. The effects of system dimension, accuracy of the approximations to capture inter-system oscillations, and the generality of the techniques are discussed in detail in the parametric study. Accuracy of the solution and weaknesses of the model reduction are also investigated. Direct numerical simulation of the power system dynamic model is performed to investigate the accuracy and robustness of the method. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Recent years have seen the emergence of a number of model reduction techniques based on linear analysis methods. The basic motivation for reduced-order system approximation is the need for simplified models of dynamical systems, which capture the main features of the original complex model. This need arises from limited computational, accuracy, and storage capabilities [1]. The simplified model is then used in place of the original complex model, for either simulation or control [2,3]. A significant number of methods of modal reduction are based on the retention of the dominant modes of the system in the reducedorder model. Reduced-order models are essential for extracting specific system performance and achieving real-time control of dynamic processes in power systems. Algorithms such as optimal Hankel model reduction, Krylov-based methods [4], balanced truncation [5], and proper orthogonal decomposition (POD), among others, have been widely used throughout the control and power system communities to generate reduced linear models with strong guarantees of accuracy and efficacy. These techniques have been
⇑ Corresponding author. E-mail address: [email protected] (A.R. Messina). 0142-0615/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijepes.2012.08.001
broadly used for reduced order modeling of linear and nonlinear models in a wide-range of disciplines, including signal analysis and data compression. Considerable effort has been applied towards development of algorithms that extend balanced truncation to large-scale, linear time-invariant systems [6,7]. Foremost among these techniques, the POD method [8] has emerged as a popular alternative for dimension reduction of large, stable DAE models. Proper orthogonal decomposition [9] has been used in a wide variety of disciplines such as aeroelasticity [10], and flow identification and control [11]. More recently, the technique has been successfully applied to power system data [12]. In this approach, reduced modes selected adaptively from the dataset are used to approximate system response. Despite considerable success in the development of reduced-order models (ROMs), efficient algorithms for very large systems remain a challenge, and new methods continue to be developed. Models described by DAE systems [13–15], are often of high order resulting in difficulties for online control due to the extensive computational effort. Reducing the size of the model, while retaining important system properties for controller design, is a major research issue in control-relevant model reduction techniques [16]. To be of practical use, the technique should preserve the outputs of interest and must include the effects of actuation.
294
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
In this paper, one new model reduction technique for the analysis and control of inter-area oscillations which preserve the inputs and outputs of interest is presented. The method use concepts from linear system theory and proper orthogonal decomposition, and provides a framework for the determination of efficient reducedorder dynamic models for large-scale systems, in which critical system features and structure are retained. This technique can be used to analyze large power system models described by DAE models and takes account of inputs and outputs of interest. Results are presented that demonstrate improved performance over alternative techniques that use conventional modal approximations. 2. Mathematical preliminaries Model reduction plays a crucial role in the model-based prediction of dynamical systems. This section highlights those aspects of linear systems theory that are relevant to the study of reducedorder system representations. Notation is also introduced. 2.1. Background Consider a stable multi-input multi-output linear dynamical system described by [17]
R:
_ xðtÞ ¼ AxðtÞ þ BuðtÞ yðtÞ ¼ CxðtÞ þ DuðtÞ
$ R :¼
A
B
C
D
;
ð1Þ
where A 2 Mn,n, B 2 Mn,m, C 2 Mp,n, D 2 Mp,m are matrices describing system behavior, u 2 Rm is the input vector, y 2 Rp is the output vector, and x 2 Rn is the vector of state variables. The solution of Eq. (1) with initial conditions x0 and inputs u can be written as
xðtÞ ¼ Hðx0 Þ þ LðuÞ;
ð2Þ
where H(x0) and L(u) are lineal operators defined by
Hðx0 Þ ¼ eAðtt0 Þ x0 ; LðuÞ ¼ eAt
Z
t
eAs BuðsÞds:
t0
Assume now that the system in Eq. (1) is diagonalizable. Let the matrix A have an eigenvalue set k ¼ fk1 ; k2 ; . . . ; kn g with associated eigenvectors u and reciprocal eigenvectors w. Use of the transformations x = uz, and z = u1x = wx in Eq. (1) transforms the linear dynamical system given by Eq. (1) into its Jordan canonical form
_ zðtÞ ¼ Az zðtÞ þ Mc uðtÞ;
ð3Þ
yðtÞ ¼ Mo zðtÞ þ DuðtÞ;
ð4Þ
where z is the modal state vector, Az = wAu and Mc = wB and Mo = Cu are the controllability and observability matrices, respectively. The projection matrix u can be computed from the eigenvectors of the system matrix. As discussed below, model reduction involves finding lowdimensional models that approximate the full high-dimensional dynamics in Eq. (2) and preserve modal properties. 2.2. The model reduction problem Given the linear dynamical system R in Eq. (1), the fundamental concept of model reduction is to find a reduced-order linear system Rr of the same form
Rr :
x_ r ðtÞ ¼ Ar xr ðtÞ þ Br uðtÞ yðtÞ ¼ Cr xr ðtÞ þ DuðtÞ
$ Rr :¼
Ar
Br
Cr
D
ð5Þ
where xr 2 Mk,1 is the state vector of the retained components, k is the order of the reduced-order model (ROM) with k n, and
Ar 2 Mk,k, Br 2 Mk,m, Cr 2 Mp,k, are matrices describing the reducedorder state-space system. It is assumed that the ROM preserves the relevant input–output behavior and that the reduced system Rr accurately describes the desired dynamics of the original system R with many fewer states, k n. The corresponding controllability and observability matrices of the ROM are given by:
z_ r ðtÞ ¼ Azr zr ðtÞ þ Mcr uðtÞ;
ð6Þ
yðtÞ ¼ Mor zr ðtÞ þ DuðtÞ;
ð7Þ
where zr is the dominant modal state vector, Azr contains the dominant modes (i.e., the k most energetic modes) of A, Mcr = wrB is the controllability matrix corresponding to dominant modes, Mor = Cur is the observability matrix corresponding to dominant modes, ur and wr are the right and left eigenvectors corresponding to dominant modes of A, respectively. To be of practical interest, the following properties must be satisfied [17]: (a) the approximation error (x xr) must be small; (b) system properties such as stability and modal controllability and observability are preserved; and (c) the procedure is computationally efficient. Methods which preserve exactly a limited number of parameters of the full-order model are of particular interest here. In what follows, an overview of basic model reduction methods is presented. Then, a technique is proposed that improves system characterization.
3. Projection-based ROM techniques The fundamental concept of model reduction is to identify a low-order, dominant subspace. By projecting the high-order system onto this subspace, a low-order model is obtained which accurately captures the dynamics of interest. Consider a dynamic linear system described by Eq. (1). The objective of the reduction process is to determine a kth-order reduced space basis onto which the state vector can be projected, that is:
xr ; T ¼ ½ T1 T2 ; x ¼ Txproj ; xproj ¼ ~ x Ik 0 T1i T1 T1i T2 Ti T ¼ ¼ ; 0 Ink T2i T1 T2i T2
Ti ¼ T1 ¼
T1i T2i
;
where T1 2 Mn,k, T2 2 Mn,n-k, T1i 2 Mk,n, T1 2 Mn-k,n are components of the T and Ti projection matrices, xr 2 Mk,1 contains the states of the ~ 2 M nk;1 repreROM; that is, the states with the highest energy, x sents the states with lower energy. Substituting x = Txproj into Eq. (1), results in:
~Þ þ T1i Bu; x_ r ¼ T1i AðT1 xr þ T2 x
ð8Þ
~_ ¼ T2i AðT1 xr þ T2 x ~Þ þ T2i Bu; x
ð9Þ
~Þ þ Du: y ¼ CðT1 xr þ T2 x
ð10Þ
In practice, the dimension of the linear system given by Eqs. ~ (the (8)–(10) can be reduced by neglecting the terms involving x states with lower energy); this results in a dynamical system that evolves in a lower dimensional state space. The projected state and output equations now become
x_ r ¼ T1i AT1 xr þ T1i Bu; y ¼ CT1 xr þ Du; and the associated approximation Rr of R is
295
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
Rr ¼
T1i AT1 CT1
Ar T1i B ¼ D Cr
Br ; D
ð11Þ
where xr = T1ix and x = T1xr. Several different approaches exist for systematically obtaining the basis transformation matrix T. We next offer a critical review of some existing techniques for reducing the order of large-scale systems. We then compare the performance of these methods in their ability to extract modal information in simulated data. 4. Approximation methods In this section, reduced modeling techniques based on modal approximation (MA) and a proper orthogonal decomposition method are reviewed. Then an optimal projection-based model reduction method is proposed that combines concepts from the linear modal analysis and proper orthogonal decomposition methods. This framework is constituted by three techniques of approximation; the first is the modal approximation (MA), the second is the proper orthogonal decomposition (POD) and, the third is the optimized modal approximation (OMA). 4.1. Modal approximation One widespread application of the truncation approach to model reduction, besides balanced truncation, is modal truncation. In this case, however, modal approximation is not derived from balancing. It is derived by looking at the eigenvalue decomposition (EVD) not of the product of Gramians PQ but simply of A [17]. The problem of selecting the significant modes from those undistinguishable from R is of considerable interest and has been studied extensively. Assuming that A is diagonalizable and that the system is stable, in addition the eigenvalues of A can be ordered so that Reðkkþ1 Þ 6 Reðkk Þ 6 6 Reðk1 Þ 0; the reduced system Rr obtained by truncation now preserves the k dominant modes. Now we can define the transformation matrices T1 and T1i as follows:
T1 ¼ u1 ; T1i ¼ w1 ;
ð12Þ
where u1 and w1 are the first k right and left eigenvectors corresponding to dominant modes of A respectively. The limitation of this techniques is that the modes with the most energy are often left out of the reduced-order representation and, consequently, that the reduced system Rr, obtained by truncation may not be optimal in sense of the eigenvalues of A. 4.2. Proper orthogonal decomposition Proper orthogonal decomposition (POD) has become popular as a means of extracting and using dominant (ideally optimal) energycontaining structures from measured or simulated data as basis functions for generating low-order dynamic models. In the context of this research, consider the output vector of a MIMO system of the form in Eq. (1). The observation matrix X is defined as [17],
2
X ¼ ½ x1
x2
x1 ðt 1 Þ 6 x2 ðt 1 Þ 6 xN ¼ 6 6 .. 4.
x1 ðt2 Þ x1 ðt N Þ
3
x2 ðt2 Þ x2 ðt N Þ 7 7 7; .. .. .. 7 5 . . .
X 2 Rn;N ;
xn ðt 1 Þ xn ðt2 Þ xn ðt N Þ where xi(tj) is the instantaneous system state or snapshot at time, tj, obtained from the numerical simulation of the system linear model and N is the number of samples. The objective is to find a set of orthonormal basis vectors (basis functions), {U1, U2, . . . , Uk}.
Let X = UrW⁄ denote the singular value decomposition of the observation matrix, where U and W are orthogonal matrices containing the right and left singular vectors, and r is a diagonal matrix whose diagonal elements consist of non-negative numbers ri (the singular values). It can be shown that the k basis functions, {U1, U2, . . . , Uk}, are given by the first k left singular vectors; the singular vectors represent the amount of energy captured by each basis function. The asterisk denotes complex conjugation. A reduced-order model can then be derived by assuming that the output vector xi, can be expressed approximately by a truncated expansion involving the first few basis functions, that is
xi
k X
cji Uj ; cji ¼ rj Wij ; i ¼ 1; 2; . . . ; N;
j¼1
where rj is jth singular value, Wij is the (i, j) element from W, the coefficients cji are to be determined such that xi will most resemble X in the sense that it maximizes the projection onto the allowed subspace. In practical applications, an energy criterion is used to determine the number of basis functions Uj, needed to accurately approximate the transient response. The transformation matrices T1 and T1i, can be written in the form
T1 ¼ ½ U1
U2
Uk ;
T1i ¼ T1 :
ð13Þ
The disadvantage of the POD technique is that Rr is not always stable even if the system R is stable [17], this is observes in the eigenvalues of Rr. Further, the accuracy of the POD-based ROMs is heavily dependent upon the sampling frequency and the time window in which the snapshots are taken. While the accuracy of the POD method can be enhanced by increasing the number of snapshots and POD bases, the method becomes impractical for large system models. Another drawback of these techniques is an inability to account effectively for the system inputs and outputs of interest. A good overview of POD-based methods is given in [18]. Next, a reduced-order modeling framework for approximating the transient response of large-scale linear system models that combines the advantages of modal methods and the POD technique is proposed. 4.3. Optimized modal approximation Consider the dynamic system in Eq. (1). Let now the linear transformations be defined by
x ¼ uz;
u 2 Mn;n ;
z ¼ wx;
w 2 Mn;n ;
where u and w are the right and left eigenvectors of A. Using the above coordinate transformations, the linear dynamic system in Eq. (1) can be expressed as
z_ ¼ wAuz þ wBu; z_ ¼ Az z þ Bz u; where, z is the modal state vector. _ k Þ be defined Let now the unit vector z_ u ðtk Þ in the direction of zðt as
z_ u ðt k Þ ¼
_ kÞ zðt ; _ k Þk2 kzðt
kz_ u ðtk Þk2 ¼ kz_ u ðtk Þk22 ¼ 1;
ð14Þ
It then follows that
kz_ u ðtk Þk22 ¼
Pn
_
j¼1 ½zj ðt k Þ
½z_ j ðt k Þ
_ k Þk22 kzðt
Pn ¼
_ _ j¼1 ½zj ðt k Þ ½zj ðt k Þ _ k Þ _ k Þ ½zðt ½zðt
¼ 1;
296
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
or
4.4. Numerical algorithm r X
n X
j¼1
j¼rþ1
ð15Þ
A step-by-step description of the algorithm used for optimal projection-based model reduction is described here.
in which it is assumed that E1(tk) P E2(tk) P P En(tk) P 0, and the terms Ej(tk), j = 1, . . . , r are associated to the states with the highest energy. It then follows that the fraction of total energy captured by a set of relevant modes, j = 1, . . . , r is given by
(1) Given an initial operating condition, compute the right u an w and left eigenvectors of the state-space matrix A. Compute _ 1 Þ with z(t1) = wx(t1). linear solutions for zðt (2) Estimate the associated energies using Eqs. (15)–(17). (3) Determine the k most energetic modes using the Eq. (16). and determine the dynamic system Rr ; T1i ¼ w (4) Set T1 ¼ u using Eq. (11). (5) Compute xr ðt1 Þ ¼ A1 r T1i Axðt 1 Þ and determine the solution in time domain of the dynamic system Rr. (6) The final step of the method is to project the dynamic behavior of the Rr using the Eq. (22), i.e., x(t) = A1(T1Ar) xr(t) + A1(T1T1i In)Bu(t).
1¼
Ej ðtk Þ þ
Ecaptured
Ej ðtk Þ;
! r X 6 100 Ej ðt k Þ ;
ð16Þ
j¼1
where, in practice, Ecaptured is a predefined energy percentage, i.e. 99.9% and
Ej ðtk Þ ¼
n X ½z_ j ðtk Þ ½z_ j ðtk Þ wji x_ i ðtk Þ; ; 8z_ j ðt k Þ ¼ _ _ ½zðtk Þ ½zðtk Þ i¼1
ð17Þ
Substituting Eq. (17) into Eq. (16) we obtain the first r modes that satisfy the energy balance Eq. (16). We can define the projec as and T1i ðt k Þ ¼ w tion matrices T1 ðtk Þ ¼ u 2 3 2 3 u11 u12 u1r w11 w12 w1n 6u 7 6w 7 6 21 u22 u2r 7 6 21 w22 w2n 7 7 6 7 T1 ðtk Þ ¼ 6 .. .. .. 7: .. .. .. 7; T1i ðtk Þ ¼ 6 .. 6 .. 4. 4. . . . 5 . . . 5 wr1 wr2 wrn un1 un2 unr With this approach, it becomes possible to study model reduction in an analytical setting where the inputs and outputs of interest are preserved. These transformation matrices, however, are not optimized to reduce the approximation error kx_ T1 x_ r k22 . 4.3.1. Optimization based on the error kx_ T1 x_ r k22 Based on the above formulation, the error kx_ T1 x_ r k22 can be defined as
e ¼ kx_ T1 x_ r k22 ¼ ðx_ T1 x_ r Þ ðx_ T1 x_ r Þ:
ð18Þ
Using the properties of the L2-norm, we can rewrite the Eq. (18) as
e ¼ ½Ax þ Bu T1 ðAr xr þ Br uÞ ½Ax þ Bu T1 ðAr xr þ Br uÞ:
ð19Þ
An optimal reduced order model of the system can be obtained by choosing the solution with minimum L2-norm. Differentiating in Eq. (19) with respect to x and xr, we obtain
@ e ¼ 2ðAxÞ A 2ðT1 Ar xr Þ A þ 2ðBuÞ A 2ðT1 Br uÞ A ¼ 0; @x @ e ¼ 2ðAxÞ ðT1 Ar Þ þ 2ðT1 Ar xr Þ ðT1 Ar Þ 2ðBuÞ ðT1 Ar Þ @xr þ 2ðT1 Br uÞ ðT1 Ar Þ ¼ 0;
ð20Þ
1
xðtÞ ¼ A ðT1 Ar Þxr ðtÞ þ A ðT1 T1i In ÞBuðtÞ;
The methodology described applies in the general case of nonlinear dynamic models described by an implicit mixed set of differential and algebraic equations (DAEs) arise when modeling power system dynamic processes. To motivate our analysis, consider a general power system represented by the state-space realization
x 2 Rn ;
_ xðtÞ ¼ fðx; y; u; tÞ;
u 2 Rm ;
y 2 Rp ;
ð23Þ
0 ¼ gðx; yÞ;
where x is the vector of state variables, y is the vector of algebraic variables, and u is the input vector. For the ease of discussion, each synchronous generator is represented using the two-axis model and a static excitation system. Machine saturation is neglected [19]. This approach has been adopted here as it offers simplicity and can be easily extended to the case of more detailed system representations. The differential and algebraic equations describing the steady state and dynamic behavior for the ith generator and the excitation system are given by [19]. 5.1. Linear system representation Assuming that f(x, y, u, t) is continuous and can be expanded, the Taylor power series expansion up to first order of Eq. (23) about a stable equilibrium point results in
Dx_ ¼ Ax Dx þ AG DVG þ Au Du ð24Þ
0 ¼ Bx Dx þ BG DVG þ BL DVL
ð21Þ 1
5. Modeling framework
0 ¼ C G DV G þ C L DV L ;
Finally, solving Eq. (20) for x and xr, results in
xr ðtÞ ¼ A1 r T1i AxðtÞ;
The next sections describe in more detail that nature of the approximations used here to construct ROMs of large-scale systems.
ð22Þ
where In is the n n identity matrix. This approach results in the best possible energy transformations in the mean energy sense. Eqs. (21) and (22) are the new transformation matrices that minimize the error e. These equations are used by the OMA technique. It is noted that this procedure is more general than the conventional approaches proposed in previous sections.
where
Dx ¼ ½DE0q1 DR f 1
DE0qng
DRfng
DE0d1 DV r1
DE0dng DV rng
DEfd1 Dd1
Dxrng T ;
DVG ¼ ½Dh1
¼ ½Dhngþ1
Du ¼ ½DPm1
Dhng
Dh n DPmng
DV 1
DV ngþ1 DV ref 1
DV ng T ; DVL
DV n T ; DV refng T ;
DEfdng
Ddng
Dxr1
297
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
Fig. 1. Single-line diagram of the IEEE 50-machine test system.
Table 1 Comparison of inter-area modes of ROM given by the reduction techniques and FSM. Scenario 1. Number mode
FSM (350 states)
MA-ROM (230 states)
OMA-ROM (62 states)
POD-ROM (70 states)
1 2 3 4
1.9508 ± j2.6669 2.6942 ± j3.5691 0.9017 ± j6.3225 2.6566 ± j6.8713
1.9508 ± j2.6669 2.6942 ± j3.5691 0.9017 ± j6.3225 2.6566 ± j6.8713
– 2.6942 ± j3.5691 – 2.6566 ± j6.8713
2.0724 ± j3.1062 0.9628 ± j3.4098 0.9076 ± j6.3141 3.0325 ± j6.8082
Fig. 2. Percentage energy content for the modes in Table 1 (OMA-ROM). Scenario 1. Fig. 3. Rotor speed deviation of machine G27. Scenario 1.
,
Ax ¼
¼
@ f1 @x
@ g @VG ngþ1
@ fng @x
T
; Bx ¼
@ g @VG n
T ;
@ g @x 1
@ g @x ng
T
; CG
T @ @ @ f1 fng ; BG ¼ g @VG @VG @VG 1 T @ @ ¼ g g ; @VL ngþ1 @VL n
AG ¼
@ g @VG ng
T ; CL
298
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
Fig. 4. Percentage energy content for the modes in Table 2 (OMA-ROM). Scenario 2.
Au ¼
@ f1 @u
@ fng @u
T
; BL ¼
@ g @VL 1
@ g @VL ng
Fig. 5. Rotor speed deviation of machine G7. Scenario 2.
T :
Further, the constraint functions given in Eq. (24) are related as follows
DVL ¼ CGL DVG
ð25Þ
DVG ¼ BxG Dx; where 1 CGL ¼ C1 L CG ; BxG ¼ ðBG þ BL CGL Þ Bx ;
Substituting (25) into (24), leads to the linear DAE model
Dx_ ¼ Axc Dx þ Auc Du; "
DVG DVL
#
¼
BxG CGL BxG
For the present study, a sixth-order machine model with a fast static exciter was used for all machines. The full system model (FSM) comprises 350 differential equations, and 200 algebraic equations. The aim of this analysis is to seek a low-dimensional dynamic description of system behavior which preserves the most energetic modes and is suitable for analysis and control purposes. The performance of the various reduced-order modeling techniques is evaluated in terms of numerical accuracy, order reduction, and computational burden under different operating conditions. Simulation results are benchmarked against direct numerical simulation of the original (unreduced) system model. Three representations were considered in the studies, namely:
(1) A reduced order model obtained from the proper orthogonal decomposition method (POD-ROM). (2) A reduced order model obtained from the modal approximation method (MA-ROM). (3) The new optimized modal approximation method based ROM.
ð26Þ
Dx;
where Axc = Ax + AGBxG. Eq. (26) constitutes the analytical model used in this document. The derived model is consistent with large-scale system representations and permits consideration of general dynamical systems. Once an efficient method for computing the linear system representation is available, it is possible to use any of the above models to identify a ROM. The main consideration is then how efficiently this reduction may be performed. In the following, results will be presented that compare the above methods, and the ROMs that they produce.
The solutions of the OMA-ROM technique were compared with the full model solution for two case studies: Case 1. Simultaneous tripping without fault the transmission lines 1–6 and 2–6 for 0.2 s. This scenario is known from observability studies to excite inter-area modes 2 and 4. Case 2. Simultaneous 0.10 p.u. step change in the mechanical power of machines G7, G13, G15 and G16. These machines participate strongly in inter-area modes 1 and 3.
6. Study results The ability of the reduced-order modeling techniques to capture the relevant inter-area dynamics of multi-machine systems has been studied on the IEEE 50-machine test system [20,21]. Fig. 1 shows a single-line diagram of the test system showing the location of selected generators and test sites.
6.1. Case study 1 Using the procedures described in Section 4, a ROM was determined using the proposed procedures. Table 1 gives the dominant
Table 2 Comparison of inter-area modes of ROM given by the reduction techniques and FSM. Scenario 2. Number mode
FSM (350 states)
MA-ROM (350 states)
OMA-ROM (12 states)
POD-ROM (16 states)
1 2 3 4
1.9508 ± j2.6669 2.6942 ± j3.5691 0.9017 ± j6.3225 2.6566 ± j6.8713
1.9508 ± j2.6669 2.6942 ± j3.5691 0.9017 ± j6.3225 2.6566 ± j6.8713
1.9508 ± j2.6669 – 0.9017 ± j6.3225 –
2.3028 ± j2.4852 – 0.9343 ± j6.1538 –
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
(most energetic modes) captured by the MA, OMA and POD-based ROMs. Column 2 in Table 1 shows the eigenvalues of the FSM (unreduced system). Fig. 2 gives the percentage of energy captured by the dominant modes. In all cases, the size of the resulting ROM representations is significantly lower than that of the original system (1/6–1/5 the size of the original system representation). The MA-ROM has 230 states, the OMA-ROM has 62 states and the POD-ROM has 70 states. For the results presented in this study, the observation matrix was derived from the normalized ensemble of state snapshots X (refer to Section 4.2) were obtained from a detailed time domain simulation of the linear model in Eq. (26). 1750 snapshots of data were employed in order to construct the reduced-order models. These snapshots were then used to construct the projection matrices T1, T1i in Eq. (13). Several general observations can be drawn from this analysis. As shown in Table 1, the MA-ROM tends to capture several major inter-area modes at the expense of a high reduced-order model. The POD-ROM, on the other hand, fails to properly characterize inter-area modes 1, 2 and 4. In order to provide a basis of assessing computational efficiency, the unreduced system was integrated numerically and the computational time was compared with that of the reduced order models. The CPU times are: 1.25 s for the MA-ROM, 0.40 s for the OMAROM and 0.44 s for the POD-ROM. As a reference, computation of the full system response using a transient stability program takes about 1.90 s. To further quantify the accuracy of the optimized and PODbased ROMs, a relative global error (mean square error) between the exact solution and the numerical approximations was evaluated. Noting that the largest singular value of a matrix X is equal to its induced L2-norm, i.e. r1 = kXk2, we define the relative error (RE) between the full system response XFSM and the approximate ROM response XROM, as
RE ¼
kXFSM XROM k2 rFSMROM ¼ ; kXFSM k2 rFSM
where rFSM-ROM is the largest singular value of XFSM XROM, rFSM is the largest singular value of XFSM; the relative error ranges from 0 to 1. The relative errors (REs) of the various approximations relative to the FSM are: 0.067 for the MA-ROM, 0.055 for the OMA-ROM and 0.0005 for the POD-ROM. Fig. 3 provides a comparison of the speed deviation of machine G27 obtained from direct numerical integration of the full system models with the solution of the various reduced-order representations. As shown in this plot, the agreement is very good over the entire study period.
6.2. Case study 2 For case 2, Fig. 4 gives the captured energy for the modes retained in the OMA-ROM. A summary of the results obtained by the various methods is given in Table 2. The MA-ROM has 350 states, the OMA-ROM has 12 states and the ROM-POD has 16 states. The computational times needed to compute the reduced order models are: 0.755 s for the MA-ROM, 0.062 s for the OMA-ROM and 0.098 s for the POD-ROM. Direct numerical simulation of the unreduced system model takes about 0.755 s. Fig. 5 shows the variation of the rotor speed deviation of machine G7. The relative error (RE) of each reduction method with respect to the full system model is: 1.271 1016 for the MA-ROM, 0.0030 for the OMA-ROM technique and 0.0001 for the POD-ROM.
299
7. Conclusions Reduced-order modeling represents an attractive approach for the analysis and control of inter-area oscillations in large power system models. In this research, we introduce a new linear projection for large systems of differential–algebraic equations. This approach reduces the quadratic error between the exact solutions and the reduced-order model approximation and results in an optimum low-dimensional description of the system dynamics. The approach is particularly well suited to treat control or design problems of large-scale linear systems governed by differential–algebraic equations. Experience with the analysis of various test systems suggest that the proposed ROM approach is well suited for analysis of large-scale dynamical systems and that the number of basis vectors needs to be increased only linearly with the number of retained states to maintain a particular ROM performance. Preliminary results show that the optimized modal approximation method is capable of representing the linear response over a range of temporal scales and can be to identify highest energy modes. The key advantages of the proposed model over more formal approaches are simplicity and robustness. Future work should focus on model reduction of large-scale dynamical systems involving multiple FACTS controllers and the application of model-based control design. The promising results for linear systems motivate the extension of the proposed projection methods to nonlinear systems. References [1] Ghosh Sudipta, Senroy Nilanjan. The localness of electromechanical oscillations in power systems. Electr Power Energy Syst 2012;42(1):306–13. [2] de Campos VAF, da Cruz JJ, Zanetta Jr LC. Robust and optimal adjustment of power system stabilizers through linear matrix inequalities. Electr Power Energy Syst 2012;42:478–86. [3] Hashmani Ashfaque A, Erlich István. Mode selective damping of power system electromechanical oscillations for large power systems using supplementary remote signals. Electr Power Energy Syst 2012;42:605–13. [4] Young JS, Wie LF. Super-optimal model reduction in sense of Hankel-norm. In: Proceedings of the 2004 IEEE international conference on networking, sensing and control, Taipei, Taiwan, 21–23; March 2004. [5] Sandberg H. Model reduction of linear systems using extended balanced truncation. In: American control conference, Seattle, WA, 11–13; June 2008. [6] Benner P, Quintana-Orti E, Quintana-Orti G. State-space truncation methods for parallel model reduction of large-scale systems. In: Proceedings of 43rd IEEE conference on decision and control, vol. 23; 2004. p. 3078–83. [7] Gugercin S, Antoulas A. A survey of model reduction by balanced truncation and some new results. Int J Control 2004;77:748–66. [8] Lin WZ, Zhang YJ, Li EP. Proper orthogonal decomposition in the generation of reduced order models for interconnects. IEEE Trans Adv Packag 2008;31(3): 627–636. [9] Kunisch K, Volkwein S. Proper orthogonal decomposition for optimality systems. ESAIM: Math Model Numer Anal 2008;42:1–23. [10] Romanowski M. Reduced order unsteady aerodynamic and aeroelastic models using Karhunen–Loéve eigenmodes. AIAA Paper 1996:96–194. [11] Kunisch K, Volkwein S. Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics. SIAM J Numer Anal 2002;40:492–515. [12] Messina AR, Vittal V. Extraction of dynamics pattern from wide area measurements using empirical orthogonal functions. IEEE Trans Power Syst 2005;22(2):682–92. [13] Chowdhry S, Linninger AA. Automatic structure analysis of large scale differential algebraic systems. In: IEEE instrumentation and measurement technology conference, 21–23; May 2001. [14] Byrne GD, Ponzi PR. Differential–algebraic systems, their application and solutions. Comput Chem Eng 1998;12(5):377–82. [15] Li YH, Yuan WP, Chan KW, Liu MB. Coordinated preventive control of transient stability with multi-contingency in power systems using trajectory sensitivities. Electr Power Energy Syst 2011;33:147–53. [16] Sun C, Hahn J. Reduction of stable differential–algebraic equation systems via projections and system identification. J Process Control 2005;15:639–50. [17] López Ríos A, Messina AR. An energy-based modal approximation method for reduced-order modeling analysis of power system models. Electric Power Compon Syst 2011;39(14):1523–41. [18] Kerschen G, Galinval JC, Vakakis AF, Bergman LA. The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: an overview. Nonlinear Dynam (Special issue on Reduced Order Models: Methods and Applications) 2005;41:147–69.
300
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
[19] Sauer PW, Pai MA. Power system dynamics and stability. New Jersey: Prentice Hall; 1998. p. 161–9. [20] IEEE Stability Test Systems Task Force of the Dynamic System Performance Subcommittee. Transient stability test systems for direct stability methods. IEEE Trans Power Syst 1992; 7(1): 37–43.
[21] Power System Toolbox Ver. 2.0. Dynamic tutorial and functions. Colborne, ON, Canada: Cherry Tree Scientific Software; 1999.
Contents lists available at SciVerse ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
An optimal modal approximation method for model reduction of linear power system models A. López Ríos, A.R. Messina ⇑ Department of Electrical Engineering, The Center for Research and Advanced Studies (Cinvestav), Av. del Bosque 1145, colonia el Bajío, Zapopan, CP 45019 Jalisco, Mexico
a r t i c l e
i n f o
Article history: Received 1 February 2011 Received in revised form 1 August 2012 Accepted 1 August 2012 Available online 26 September 2012 Keywords: Modal approximation Proper orthogonal decomposition Optimized modal approximation
a b s t r a c t A number of recent studies have examined the problem of modal reduction of large-scale power system linear models. In this paper, an optimal model reduction strategy for the analysis and control of inter-area oscillations in power systems that uses concepts from proper orthogonal decomposition and modal analysis and preserves exactly the input–output properties of the original system is presented. The technique uses a new projection matrix that minimizes the error between the states of the original system and the states of its reduced-order model, and is suitable for analyzing large power system models described by differential–algebraic equation (DAE) systems. Based on this theoretical framework, an efficient algorithm for model reduction is proposed which combines the ease of use of linear modal analysis, with the flexibility of extracting relevant behavior that statistical multivariate methods can offer. Methods for analysis of reduced system representations using linear approximations are critically reviewed and compared. The proposed technique is tested on a 50-machine, 145-bus, and 453-line test system. The effects of system dimension, accuracy of the approximations to capture inter-system oscillations, and the generality of the techniques are discussed in detail in the parametric study. Accuracy of the solution and weaknesses of the model reduction are also investigated. Direct numerical simulation of the power system dynamic model is performed to investigate the accuracy and robustness of the method. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Recent years have seen the emergence of a number of model reduction techniques based on linear analysis methods. The basic motivation for reduced-order system approximation is the need for simplified models of dynamical systems, which capture the main features of the original complex model. This need arises from limited computational, accuracy, and storage capabilities [1]. The simplified model is then used in place of the original complex model, for either simulation or control [2,3]. A significant number of methods of modal reduction are based on the retention of the dominant modes of the system in the reducedorder model. Reduced-order models are essential for extracting specific system performance and achieving real-time control of dynamic processes in power systems. Algorithms such as optimal Hankel model reduction, Krylov-based methods [4], balanced truncation [5], and proper orthogonal decomposition (POD), among others, have been widely used throughout the control and power system communities to generate reduced linear models with strong guarantees of accuracy and efficacy. These techniques have been
⇑ Corresponding author. E-mail address: [email protected] (A.R. Messina). 0142-0615/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijepes.2012.08.001
broadly used for reduced order modeling of linear and nonlinear models in a wide-range of disciplines, including signal analysis and data compression. Considerable effort has been applied towards development of algorithms that extend balanced truncation to large-scale, linear time-invariant systems [6,7]. Foremost among these techniques, the POD method [8] has emerged as a popular alternative for dimension reduction of large, stable DAE models. Proper orthogonal decomposition [9] has been used in a wide variety of disciplines such as aeroelasticity [10], and flow identification and control [11]. More recently, the technique has been successfully applied to power system data [12]. In this approach, reduced modes selected adaptively from the dataset are used to approximate system response. Despite considerable success in the development of reduced-order models (ROMs), efficient algorithms for very large systems remain a challenge, and new methods continue to be developed. Models described by DAE systems [13–15], are often of high order resulting in difficulties for online control due to the extensive computational effort. Reducing the size of the model, while retaining important system properties for controller design, is a major research issue in control-relevant model reduction techniques [16]. To be of practical use, the technique should preserve the outputs of interest and must include the effects of actuation.
294
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
In this paper, one new model reduction technique for the analysis and control of inter-area oscillations which preserve the inputs and outputs of interest is presented. The method use concepts from linear system theory and proper orthogonal decomposition, and provides a framework for the determination of efficient reducedorder dynamic models for large-scale systems, in which critical system features and structure are retained. This technique can be used to analyze large power system models described by DAE models and takes account of inputs and outputs of interest. Results are presented that demonstrate improved performance over alternative techniques that use conventional modal approximations. 2. Mathematical preliminaries Model reduction plays a crucial role in the model-based prediction of dynamical systems. This section highlights those aspects of linear systems theory that are relevant to the study of reducedorder system representations. Notation is also introduced. 2.1. Background Consider a stable multi-input multi-output linear dynamical system described by [17]
R:
_ xðtÞ ¼ AxðtÞ þ BuðtÞ yðtÞ ¼ CxðtÞ þ DuðtÞ
$ R :¼
A
B
C
D
;
ð1Þ
where A 2 Mn,n, B 2 Mn,m, C 2 Mp,n, D 2 Mp,m are matrices describing system behavior, u 2 Rm is the input vector, y 2 Rp is the output vector, and x 2 Rn is the vector of state variables. The solution of Eq. (1) with initial conditions x0 and inputs u can be written as
xðtÞ ¼ Hðx0 Þ þ LðuÞ;
ð2Þ
where H(x0) and L(u) are lineal operators defined by
Hðx0 Þ ¼ eAðtt0 Þ x0 ; LðuÞ ¼ eAt
Z
t
eAs BuðsÞds:
t0
Assume now that the system in Eq. (1) is diagonalizable. Let the matrix A have an eigenvalue set k ¼ fk1 ; k2 ; . . . ; kn g with associated eigenvectors u and reciprocal eigenvectors w. Use of the transformations x = uz, and z = u1x = wx in Eq. (1) transforms the linear dynamical system given by Eq. (1) into its Jordan canonical form
_ zðtÞ ¼ Az zðtÞ þ Mc uðtÞ;
ð3Þ
yðtÞ ¼ Mo zðtÞ þ DuðtÞ;
ð4Þ
where z is the modal state vector, Az = wAu and Mc = wB and Mo = Cu are the controllability and observability matrices, respectively. The projection matrix u can be computed from the eigenvectors of the system matrix. As discussed below, model reduction involves finding lowdimensional models that approximate the full high-dimensional dynamics in Eq. (2) and preserve modal properties. 2.2. The model reduction problem Given the linear dynamical system R in Eq. (1), the fundamental concept of model reduction is to find a reduced-order linear system Rr of the same form
Rr :
x_ r ðtÞ ¼ Ar xr ðtÞ þ Br uðtÞ yðtÞ ¼ Cr xr ðtÞ þ DuðtÞ
$ Rr :¼
Ar
Br
Cr
D
ð5Þ
where xr 2 Mk,1 is the state vector of the retained components, k is the order of the reduced-order model (ROM) with k n, and
Ar 2 Mk,k, Br 2 Mk,m, Cr 2 Mp,k, are matrices describing the reducedorder state-space system. It is assumed that the ROM preserves the relevant input–output behavior and that the reduced system Rr accurately describes the desired dynamics of the original system R with many fewer states, k n. The corresponding controllability and observability matrices of the ROM are given by:
z_ r ðtÞ ¼ Azr zr ðtÞ þ Mcr uðtÞ;
ð6Þ
yðtÞ ¼ Mor zr ðtÞ þ DuðtÞ;
ð7Þ
where zr is the dominant modal state vector, Azr contains the dominant modes (i.e., the k most energetic modes) of A, Mcr = wrB is the controllability matrix corresponding to dominant modes, Mor = Cur is the observability matrix corresponding to dominant modes, ur and wr are the right and left eigenvectors corresponding to dominant modes of A, respectively. To be of practical interest, the following properties must be satisfied [17]: (a) the approximation error (x xr) must be small; (b) system properties such as stability and modal controllability and observability are preserved; and (c) the procedure is computationally efficient. Methods which preserve exactly a limited number of parameters of the full-order model are of particular interest here. In what follows, an overview of basic model reduction methods is presented. Then, a technique is proposed that improves system characterization.
3. Projection-based ROM techniques The fundamental concept of model reduction is to identify a low-order, dominant subspace. By projecting the high-order system onto this subspace, a low-order model is obtained which accurately captures the dynamics of interest. Consider a dynamic linear system described by Eq. (1). The objective of the reduction process is to determine a kth-order reduced space basis onto which the state vector can be projected, that is:
xr ; T ¼ ½ T1 T2 ; x ¼ Txproj ; xproj ¼ ~ x Ik 0 T1i T1 T1i T2 Ti T ¼ ¼ ; 0 Ink T2i T1 T2i T2
Ti ¼ T1 ¼
T1i T2i
;
where T1 2 Mn,k, T2 2 Mn,n-k, T1i 2 Mk,n, T1 2 Mn-k,n are components of the T and Ti projection matrices, xr 2 Mk,1 contains the states of the ~ 2 M nk;1 repreROM; that is, the states with the highest energy, x sents the states with lower energy. Substituting x = Txproj into Eq. (1), results in:
~Þ þ T1i Bu; x_ r ¼ T1i AðT1 xr þ T2 x
ð8Þ
~_ ¼ T2i AðT1 xr þ T2 x ~Þ þ T2i Bu; x
ð9Þ
~Þ þ Du: y ¼ CðT1 xr þ T2 x
ð10Þ
In practice, the dimension of the linear system given by Eqs. ~ (the (8)–(10) can be reduced by neglecting the terms involving x states with lower energy); this results in a dynamical system that evolves in a lower dimensional state space. The projected state and output equations now become
x_ r ¼ T1i AT1 xr þ T1i Bu; y ¼ CT1 xr þ Du; and the associated approximation Rr of R is
295
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
Rr ¼
T1i AT1 CT1
Ar T1i B ¼ D Cr
Br ; D
ð11Þ
where xr = T1ix and x = T1xr. Several different approaches exist for systematically obtaining the basis transformation matrix T. We next offer a critical review of some existing techniques for reducing the order of large-scale systems. We then compare the performance of these methods in their ability to extract modal information in simulated data. 4. Approximation methods In this section, reduced modeling techniques based on modal approximation (MA) and a proper orthogonal decomposition method are reviewed. Then an optimal projection-based model reduction method is proposed that combines concepts from the linear modal analysis and proper orthogonal decomposition methods. This framework is constituted by three techniques of approximation; the first is the modal approximation (MA), the second is the proper orthogonal decomposition (POD) and, the third is the optimized modal approximation (OMA). 4.1. Modal approximation One widespread application of the truncation approach to model reduction, besides balanced truncation, is modal truncation. In this case, however, modal approximation is not derived from balancing. It is derived by looking at the eigenvalue decomposition (EVD) not of the product of Gramians PQ but simply of A [17]. The problem of selecting the significant modes from those undistinguishable from R is of considerable interest and has been studied extensively. Assuming that A is diagonalizable and that the system is stable, in addition the eigenvalues of A can be ordered so that Reðkkþ1 Þ 6 Reðkk Þ 6 6 Reðk1 Þ 0; the reduced system Rr obtained by truncation now preserves the k dominant modes. Now we can define the transformation matrices T1 and T1i as follows:
T1 ¼ u1 ; T1i ¼ w1 ;
ð12Þ
where u1 and w1 are the first k right and left eigenvectors corresponding to dominant modes of A respectively. The limitation of this techniques is that the modes with the most energy are often left out of the reduced-order representation and, consequently, that the reduced system Rr, obtained by truncation may not be optimal in sense of the eigenvalues of A. 4.2. Proper orthogonal decomposition Proper orthogonal decomposition (POD) has become popular as a means of extracting and using dominant (ideally optimal) energycontaining structures from measured or simulated data as basis functions for generating low-order dynamic models. In the context of this research, consider the output vector of a MIMO system of the form in Eq. (1). The observation matrix X is defined as [17],
2
X ¼ ½ x1
x2
x1 ðt 1 Þ 6 x2 ðt 1 Þ 6 xN ¼ 6 6 .. 4.
x1 ðt2 Þ x1 ðt N Þ
3
x2 ðt2 Þ x2 ðt N Þ 7 7 7; .. .. .. 7 5 . . .
X 2 Rn;N ;
xn ðt 1 Þ xn ðt2 Þ xn ðt N Þ where xi(tj) is the instantaneous system state or snapshot at time, tj, obtained from the numerical simulation of the system linear model and N is the number of samples. The objective is to find a set of orthonormal basis vectors (basis functions), {U1, U2, . . . , Uk}.
Let X = UrW⁄ denote the singular value decomposition of the observation matrix, where U and W are orthogonal matrices containing the right and left singular vectors, and r is a diagonal matrix whose diagonal elements consist of non-negative numbers ri (the singular values). It can be shown that the k basis functions, {U1, U2, . . . , Uk}, are given by the first k left singular vectors; the singular vectors represent the amount of energy captured by each basis function. The asterisk denotes complex conjugation. A reduced-order model can then be derived by assuming that the output vector xi, can be expressed approximately by a truncated expansion involving the first few basis functions, that is
xi
k X
cji Uj ; cji ¼ rj Wij ; i ¼ 1; 2; . . . ; N;
j¼1
where rj is jth singular value, Wij is the (i, j) element from W, the coefficients cji are to be determined such that xi will most resemble X in the sense that it maximizes the projection onto the allowed subspace. In practical applications, an energy criterion is used to determine the number of basis functions Uj, needed to accurately approximate the transient response. The transformation matrices T1 and T1i, can be written in the form
T1 ¼ ½ U1
U2
Uk ;
T1i ¼ T1 :
ð13Þ
The disadvantage of the POD technique is that Rr is not always stable even if the system R is stable [17], this is observes in the eigenvalues of Rr. Further, the accuracy of the POD-based ROMs is heavily dependent upon the sampling frequency and the time window in which the snapshots are taken. While the accuracy of the POD method can be enhanced by increasing the number of snapshots and POD bases, the method becomes impractical for large system models. Another drawback of these techniques is an inability to account effectively for the system inputs and outputs of interest. A good overview of POD-based methods is given in [18]. Next, a reduced-order modeling framework for approximating the transient response of large-scale linear system models that combines the advantages of modal methods and the POD technique is proposed. 4.3. Optimized modal approximation Consider the dynamic system in Eq. (1). Let now the linear transformations be defined by
x ¼ uz;
u 2 Mn;n ;
z ¼ wx;
w 2 Mn;n ;
where u and w are the right and left eigenvectors of A. Using the above coordinate transformations, the linear dynamic system in Eq. (1) can be expressed as
z_ ¼ wAuz þ wBu; z_ ¼ Az z þ Bz u; where, z is the modal state vector. _ k Þ be defined Let now the unit vector z_ u ðtk Þ in the direction of zðt as
z_ u ðt k Þ ¼
_ kÞ zðt ; _ k Þk2 kzðt
kz_ u ðtk Þk2 ¼ kz_ u ðtk Þk22 ¼ 1;
ð14Þ
It then follows that
kz_ u ðtk Þk22 ¼
Pn
_
j¼1 ½zj ðt k Þ
½z_ j ðt k Þ
_ k Þk22 kzðt
Pn ¼
_ _ j¼1 ½zj ðt k Þ ½zj ðt k Þ _ k Þ _ k Þ ½zðt ½zðt
¼ 1;
296
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
or
4.4. Numerical algorithm r X
n X
j¼1
j¼rþ1
ð15Þ
A step-by-step description of the algorithm used for optimal projection-based model reduction is described here.
in which it is assumed that E1(tk) P E2(tk) P P En(tk) P 0, and the terms Ej(tk), j = 1, . . . , r are associated to the states with the highest energy. It then follows that the fraction of total energy captured by a set of relevant modes, j = 1, . . . , r is given by
(1) Given an initial operating condition, compute the right u an w and left eigenvectors of the state-space matrix A. Compute _ 1 Þ with z(t1) = wx(t1). linear solutions for zðt (2) Estimate the associated energies using Eqs. (15)–(17). (3) Determine the k most energetic modes using the Eq. (16). and determine the dynamic system Rr ; T1i ¼ w (4) Set T1 ¼ u using Eq. (11). (5) Compute xr ðt1 Þ ¼ A1 r T1i Axðt 1 Þ and determine the solution in time domain of the dynamic system Rr. (6) The final step of the method is to project the dynamic behavior of the Rr using the Eq. (22), i.e., x(t) = A1(T1Ar) xr(t) + A1(T1T1i In)Bu(t).
1¼
Ej ðtk Þ þ
Ecaptured
Ej ðtk Þ;
! r X 6 100 Ej ðt k Þ ;
ð16Þ
j¼1
where, in practice, Ecaptured is a predefined energy percentage, i.e. 99.9% and
Ej ðtk Þ ¼
n X ½z_ j ðtk Þ ½z_ j ðtk Þ wji x_ i ðtk Þ; ; 8z_ j ðt k Þ ¼ _ _ ½zðtk Þ ½zðtk Þ i¼1
ð17Þ
Substituting Eq. (17) into Eq. (16) we obtain the first r modes that satisfy the energy balance Eq. (16). We can define the projec as and T1i ðt k Þ ¼ w tion matrices T1 ðtk Þ ¼ u 2 3 2 3 u11 u12 u1r w11 w12 w1n 6u 7 6w 7 6 21 u22 u2r 7 6 21 w22 w2n 7 7 6 7 T1 ðtk Þ ¼ 6 .. .. .. 7: .. .. .. 7; T1i ðtk Þ ¼ 6 .. 6 .. 4. 4. . . . 5 . . . 5 wr1 wr2 wrn un1 un2 unr With this approach, it becomes possible to study model reduction in an analytical setting where the inputs and outputs of interest are preserved. These transformation matrices, however, are not optimized to reduce the approximation error kx_ T1 x_ r k22 . 4.3.1. Optimization based on the error kx_ T1 x_ r k22 Based on the above formulation, the error kx_ T1 x_ r k22 can be defined as
e ¼ kx_ T1 x_ r k22 ¼ ðx_ T1 x_ r Þ ðx_ T1 x_ r Þ:
ð18Þ
Using the properties of the L2-norm, we can rewrite the Eq. (18) as
e ¼ ½Ax þ Bu T1 ðAr xr þ Br uÞ ½Ax þ Bu T1 ðAr xr þ Br uÞ:
ð19Þ
An optimal reduced order model of the system can be obtained by choosing the solution with minimum L2-norm. Differentiating in Eq. (19) with respect to x and xr, we obtain
@ e ¼ 2ðAxÞ A 2ðT1 Ar xr Þ A þ 2ðBuÞ A 2ðT1 Br uÞ A ¼ 0; @x @ e ¼ 2ðAxÞ ðT1 Ar Þ þ 2ðT1 Ar xr Þ ðT1 Ar Þ 2ðBuÞ ðT1 Ar Þ @xr þ 2ðT1 Br uÞ ðT1 Ar Þ ¼ 0;
ð20Þ
1
xðtÞ ¼ A ðT1 Ar Þxr ðtÞ þ A ðT1 T1i In ÞBuðtÞ;
The methodology described applies in the general case of nonlinear dynamic models described by an implicit mixed set of differential and algebraic equations (DAEs) arise when modeling power system dynamic processes. To motivate our analysis, consider a general power system represented by the state-space realization
x 2 Rn ;
_ xðtÞ ¼ fðx; y; u; tÞ;
u 2 Rm ;
y 2 Rp ;
ð23Þ
0 ¼ gðx; yÞ;
where x is the vector of state variables, y is the vector of algebraic variables, and u is the input vector. For the ease of discussion, each synchronous generator is represented using the two-axis model and a static excitation system. Machine saturation is neglected [19]. This approach has been adopted here as it offers simplicity and can be easily extended to the case of more detailed system representations. The differential and algebraic equations describing the steady state and dynamic behavior for the ith generator and the excitation system are given by [19]. 5.1. Linear system representation Assuming that f(x, y, u, t) is continuous and can be expanded, the Taylor power series expansion up to first order of Eq. (23) about a stable equilibrium point results in
Dx_ ¼ Ax Dx þ AG DVG þ Au Du ð24Þ
0 ¼ Bx Dx þ BG DVG þ BL DVL
ð21Þ 1
5. Modeling framework
0 ¼ C G DV G þ C L DV L ;
Finally, solving Eq. (20) for x and xr, results in
xr ðtÞ ¼ A1 r T1i AxðtÞ;
The next sections describe in more detail that nature of the approximations used here to construct ROMs of large-scale systems.
ð22Þ
where In is the n n identity matrix. This approach results in the best possible energy transformations in the mean energy sense. Eqs. (21) and (22) are the new transformation matrices that minimize the error e. These equations are used by the OMA technique. It is noted that this procedure is more general than the conventional approaches proposed in previous sections.
where
Dx ¼ ½DE0q1 DR f 1
DE0qng
DRfng
DE0d1 DV r1
DE0dng DV rng
DEfd1 Dd1
Dxrng T ;
DVG ¼ ½Dh1
¼ ½Dhngþ1
Du ¼ ½DPm1
Dhng
Dh n DPmng
DV 1
DV ngþ1 DV ref 1
DV ng T ; DVL
DV n T ; DV refng T ;
DEfdng
Ddng
Dxr1
297
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
Fig. 1. Single-line diagram of the IEEE 50-machine test system.
Table 1 Comparison of inter-area modes of ROM given by the reduction techniques and FSM. Scenario 1. Number mode
FSM (350 states)
MA-ROM (230 states)
OMA-ROM (62 states)
POD-ROM (70 states)
1 2 3 4
1.9508 ± j2.6669 2.6942 ± j3.5691 0.9017 ± j6.3225 2.6566 ± j6.8713
1.9508 ± j2.6669 2.6942 ± j3.5691 0.9017 ± j6.3225 2.6566 ± j6.8713
– 2.6942 ± j3.5691 – 2.6566 ± j6.8713
2.0724 ± j3.1062 0.9628 ± j3.4098 0.9076 ± j6.3141 3.0325 ± j6.8082
Fig. 2. Percentage energy content for the modes in Table 1 (OMA-ROM). Scenario 1. Fig. 3. Rotor speed deviation of machine G27. Scenario 1.
,
Ax ¼
¼
@ f1 @x
@ g @VG ngþ1
@ fng @x
T
; Bx ¼
@ g @VG n
T ;
@ g @x 1
@ g @x ng
T
; CG
T @ @ @ f1 fng ; BG ¼ g @VG @VG @VG 1 T @ @ ¼ g g ; @VL ngþ1 @VL n
AG ¼
@ g @VG ng
T ; CL
298
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
Fig. 4. Percentage energy content for the modes in Table 2 (OMA-ROM). Scenario 2.
Au ¼
@ f1 @u
@ fng @u
T
; BL ¼
@ g @VL 1
@ g @VL ng
Fig. 5. Rotor speed deviation of machine G7. Scenario 2.
T :
Further, the constraint functions given in Eq. (24) are related as follows
DVL ¼ CGL DVG
ð25Þ
DVG ¼ BxG Dx; where 1 CGL ¼ C1 L CG ; BxG ¼ ðBG þ BL CGL Þ Bx ;
Substituting (25) into (24), leads to the linear DAE model
Dx_ ¼ Axc Dx þ Auc Du; "
DVG DVL
#
¼
BxG CGL BxG
For the present study, a sixth-order machine model with a fast static exciter was used for all machines. The full system model (FSM) comprises 350 differential equations, and 200 algebraic equations. The aim of this analysis is to seek a low-dimensional dynamic description of system behavior which preserves the most energetic modes and is suitable for analysis and control purposes. The performance of the various reduced-order modeling techniques is evaluated in terms of numerical accuracy, order reduction, and computational burden under different operating conditions. Simulation results are benchmarked against direct numerical simulation of the original (unreduced) system model. Three representations were considered in the studies, namely:
(1) A reduced order model obtained from the proper orthogonal decomposition method (POD-ROM). (2) A reduced order model obtained from the modal approximation method (MA-ROM). (3) The new optimized modal approximation method based ROM.
ð26Þ
Dx;
where Axc = Ax + AGBxG. Eq. (26) constitutes the analytical model used in this document. The derived model is consistent with large-scale system representations and permits consideration of general dynamical systems. Once an efficient method for computing the linear system representation is available, it is possible to use any of the above models to identify a ROM. The main consideration is then how efficiently this reduction may be performed. In the following, results will be presented that compare the above methods, and the ROMs that they produce.
The solutions of the OMA-ROM technique were compared with the full model solution for two case studies: Case 1. Simultaneous tripping without fault the transmission lines 1–6 and 2–6 for 0.2 s. This scenario is known from observability studies to excite inter-area modes 2 and 4. Case 2. Simultaneous 0.10 p.u. step change in the mechanical power of machines G7, G13, G15 and G16. These machines participate strongly in inter-area modes 1 and 3.
6. Study results The ability of the reduced-order modeling techniques to capture the relevant inter-area dynamics of multi-machine systems has been studied on the IEEE 50-machine test system [20,21]. Fig. 1 shows a single-line diagram of the test system showing the location of selected generators and test sites.
6.1. Case study 1 Using the procedures described in Section 4, a ROM was determined using the proposed procedures. Table 1 gives the dominant
Table 2 Comparison of inter-area modes of ROM given by the reduction techniques and FSM. Scenario 2. Number mode
FSM (350 states)
MA-ROM (350 states)
OMA-ROM (12 states)
POD-ROM (16 states)
1 2 3 4
1.9508 ± j2.6669 2.6942 ± j3.5691 0.9017 ± j6.3225 2.6566 ± j6.8713
1.9508 ± j2.6669 2.6942 ± j3.5691 0.9017 ± j6.3225 2.6566 ± j6.8713
1.9508 ± j2.6669 – 0.9017 ± j6.3225 –
2.3028 ± j2.4852 – 0.9343 ± j6.1538 –
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
(most energetic modes) captured by the MA, OMA and POD-based ROMs. Column 2 in Table 1 shows the eigenvalues of the FSM (unreduced system). Fig. 2 gives the percentage of energy captured by the dominant modes. In all cases, the size of the resulting ROM representations is significantly lower than that of the original system (1/6–1/5 the size of the original system representation). The MA-ROM has 230 states, the OMA-ROM has 62 states and the POD-ROM has 70 states. For the results presented in this study, the observation matrix was derived from the normalized ensemble of state snapshots X (refer to Section 4.2) were obtained from a detailed time domain simulation of the linear model in Eq. (26). 1750 snapshots of data were employed in order to construct the reduced-order models. These snapshots were then used to construct the projection matrices T1, T1i in Eq. (13). Several general observations can be drawn from this analysis. As shown in Table 1, the MA-ROM tends to capture several major inter-area modes at the expense of a high reduced-order model. The POD-ROM, on the other hand, fails to properly characterize inter-area modes 1, 2 and 4. In order to provide a basis of assessing computational efficiency, the unreduced system was integrated numerically and the computational time was compared with that of the reduced order models. The CPU times are: 1.25 s for the MA-ROM, 0.40 s for the OMAROM and 0.44 s for the POD-ROM. As a reference, computation of the full system response using a transient stability program takes about 1.90 s. To further quantify the accuracy of the optimized and PODbased ROMs, a relative global error (mean square error) between the exact solution and the numerical approximations was evaluated. Noting that the largest singular value of a matrix X is equal to its induced L2-norm, i.e. r1 = kXk2, we define the relative error (RE) between the full system response XFSM and the approximate ROM response XROM, as
RE ¼
kXFSM XROM k2 rFSMROM ¼ ; kXFSM k2 rFSM
where rFSM-ROM is the largest singular value of XFSM XROM, rFSM is the largest singular value of XFSM; the relative error ranges from 0 to 1. The relative errors (REs) of the various approximations relative to the FSM are: 0.067 for the MA-ROM, 0.055 for the OMA-ROM and 0.0005 for the POD-ROM. Fig. 3 provides a comparison of the speed deviation of machine G27 obtained from direct numerical integration of the full system models with the solution of the various reduced-order representations. As shown in this plot, the agreement is very good over the entire study period.
6.2. Case study 2 For case 2, Fig. 4 gives the captured energy for the modes retained in the OMA-ROM. A summary of the results obtained by the various methods is given in Table 2. The MA-ROM has 350 states, the OMA-ROM has 12 states and the ROM-POD has 16 states. The computational times needed to compute the reduced order models are: 0.755 s for the MA-ROM, 0.062 s for the OMA-ROM and 0.098 s for the POD-ROM. Direct numerical simulation of the unreduced system model takes about 0.755 s. Fig. 5 shows the variation of the rotor speed deviation of machine G7. The relative error (RE) of each reduction method with respect to the full system model is: 1.271 1016 for the MA-ROM, 0.0030 for the OMA-ROM technique and 0.0001 for the POD-ROM.
299
7. Conclusions Reduced-order modeling represents an attractive approach for the analysis and control of inter-area oscillations in large power system models. In this research, we introduce a new linear projection for large systems of differential–algebraic equations. This approach reduces the quadratic error between the exact solutions and the reduced-order model approximation and results in an optimum low-dimensional description of the system dynamics. The approach is particularly well suited to treat control or design problems of large-scale linear systems governed by differential–algebraic equations. Experience with the analysis of various test systems suggest that the proposed ROM approach is well suited for analysis of large-scale dynamical systems and that the number of basis vectors needs to be increased only linearly with the number of retained states to maintain a particular ROM performance. Preliminary results show that the optimized modal approximation method is capable of representing the linear response over a range of temporal scales and can be to identify highest energy modes. The key advantages of the proposed model over more formal approaches are simplicity and robustness. Future work should focus on model reduction of large-scale dynamical systems involving multiple FACTS controllers and the application of model-based control design. The promising results for linear systems motivate the extension of the proposed projection methods to nonlinear systems. References [1] Ghosh Sudipta, Senroy Nilanjan. The localness of electromechanical oscillations in power systems. Electr Power Energy Syst 2012;42(1):306–13. [2] de Campos VAF, da Cruz JJ, Zanetta Jr LC. Robust and optimal adjustment of power system stabilizers through linear matrix inequalities. Electr Power Energy Syst 2012;42:478–86. [3] Hashmani Ashfaque A, Erlich István. Mode selective damping of power system electromechanical oscillations for large power systems using supplementary remote signals. Electr Power Energy Syst 2012;42:605–13. [4] Young JS, Wie LF. Super-optimal model reduction in sense of Hankel-norm. In: Proceedings of the 2004 IEEE international conference on networking, sensing and control, Taipei, Taiwan, 21–23; March 2004. [5] Sandberg H. Model reduction of linear systems using extended balanced truncation. In: American control conference, Seattle, WA, 11–13; June 2008. [6] Benner P, Quintana-Orti E, Quintana-Orti G. State-space truncation methods for parallel model reduction of large-scale systems. In: Proceedings of 43rd IEEE conference on decision and control, vol. 23; 2004. p. 3078–83. [7] Gugercin S, Antoulas A. A survey of model reduction by balanced truncation and some new results. Int J Control 2004;77:748–66. [8] Lin WZ, Zhang YJ, Li EP. Proper orthogonal decomposition in the generation of reduced order models for interconnects. IEEE Trans Adv Packag 2008;31(3): 627–636. [9] Kunisch K, Volkwein S. Proper orthogonal decomposition for optimality systems. ESAIM: Math Model Numer Anal 2008;42:1–23. [10] Romanowski M. Reduced order unsteady aerodynamic and aeroelastic models using Karhunen–Loéve eigenmodes. AIAA Paper 1996:96–194. [11] Kunisch K, Volkwein S. Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics. SIAM J Numer Anal 2002;40:492–515. [12] Messina AR, Vittal V. Extraction of dynamics pattern from wide area measurements using empirical orthogonal functions. IEEE Trans Power Syst 2005;22(2):682–92. [13] Chowdhry S, Linninger AA. Automatic structure analysis of large scale differential algebraic systems. In: IEEE instrumentation and measurement technology conference, 21–23; May 2001. [14] Byrne GD, Ponzi PR. Differential–algebraic systems, their application and solutions. Comput Chem Eng 1998;12(5):377–82. [15] Li YH, Yuan WP, Chan KW, Liu MB. Coordinated preventive control of transient stability with multi-contingency in power systems using trajectory sensitivities. Electr Power Energy Syst 2011;33:147–53. [16] Sun C, Hahn J. Reduction of stable differential–algebraic equation systems via projections and system identification. J Process Control 2005;15:639–50. [17] López Ríos A, Messina AR. An energy-based modal approximation method for reduced-order modeling analysis of power system models. Electric Power Compon Syst 2011;39(14):1523–41. [18] Kerschen G, Galinval JC, Vakakis AF, Bergman LA. The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: an overview. Nonlinear Dynam (Special issue on Reduced Order Models: Methods and Applications) 2005;41:147–69.
300
A. López Ríos, A.R. Messina / Electrical Power and Energy Systems 44 (2013) 293–300
[19] Sauer PW, Pai MA. Power system dynamics and stability. New Jersey: Prentice Hall; 1998. p. 161–9. [20] IEEE Stability Test Systems Task Force of the Dynamic System Performance Subcommittee. Transient stability test systems for direct stability methods. IEEE Trans Power Syst 1992; 7(1): 37–43.
[21] Power System Toolbox Ver. 2.0. Dynamic tutorial and functions. Colborne, ON, Canada: Cherry Tree Scientific Software; 1999.