1258

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 23, NO. 3, AUGUST 2008

Gramian-Based Reduction Method Applied to Large Sparse Power System Descriptor Models Francisco Damasceno Freitas, Member, IEEE, Joost Rommes, and Nelson Martins, Fellow, IEEE

Abstract—This paper presents an efficient linear system reduction method that computes approximations to the controllability and observability gramians of large sparse power system descriptor models. The method exploits the fact that a Lyapunov equation solution can be decomposed into low-rank Choleski factors, which are determined by the alternating direction implicit (ADI) method. Advantages of the method are that it operates on the sparse descriptor matrices and does not require the computation of spectral projections onto deflating subspaces of finite eigenvalues, which are needed by other techniques applied to descriptor models. The gramians, which are never explicitly formed, are used to compute reduced-order state-space models for large-scale systems. Numerical results for small-signal stability power system descriptor models show that the new method is more efficient than other existing approaches. Index Terms—Control systems, descriptor systems, gramians, large-scale systems, linear systems, low-rank Choleski factor ADI, Lyapunov matrix equations, model reduction, power system, small-signal stability, sparsity.

I. INTRODUCTION

L

INEAR systems and numerical linear algebra are widely used in power system small-signal stability analysis and control. Applications include the design of controllers, eigenanalysis [1], [2] frequency response [3], time response [4], system identification [5], and oscillation monitoring and model reduction [3], [6], [7]. The intense research work on several types of model reduction reflects the fact that reduced systems can approximate well the dynamics of the full-size systems, and can be effectively used for advanced controller design and faster computation of time domain simulations [1], [2], [5], [6]–[8]. Model order reduction techniques have been applied in many areas of engineering and generally are of one of the following types: 1) moment-matching or Krylov subspace-based methods, 2) modal truncation methods, or 3) balanced truncation or gramian-based methods; see, for instance, [9] and [10] for recent overviews. In this paper, a new efficient gramian-based approach for the reduction of large-scale power system descriptor models is presented. Conventional gramian-based Manuscript received November 19, 2007. The work of J. Rommes was supported by O-MOORE-NICE!, Marie Curie Actions FP6 MTKI-CT-2006042477. Paper no. TPWRS-00842-2007. F. D. Freitas is with the Department of Electrical Engineering, University of Brasilia, Brasilia 70910-900, Brazil (e-mail: [email protected]). J. Rommes is with NXP Semiconductors, Corporate I&T/DTF, 5656 AE Eindhoven, The Netherlands (e-mail: [email protected]). N. Martins is with the CEPEL, Rio de Janeiro 20001-970, Brazil (e-mail: [email protected]). Digital Object Identifier 10.1109/TPWRS.2008.926693

methods for model reduction are based on a state-space representation [9]–[15]. Gramian-based approaches using singular or descriptor systems have also been developed [16]–[18], but result expensive for power system applications. The balanced truncation technique requires the computation of the system controllability and observability gramians [9], [12], [18]. These gramians, that are symmetric positive-definite matrices, are computed by solving two matrix equations: the continuous-time algebraic Lyapunov equations (CALEs). Classical methods use the state-space description to determine the controllability and observability gramians. More efficient methods implicitly compute the gramian as a product of two low-rank Choleski factors [19]–[23]. Methods based on lowrank Choleski factors for determining gramian approximations allow exploiting the sparsity of the system model. A method for computing gramians of descriptor systems was described in [17], which produces a reduced descriptor system representation. However, two projected generalized continuous-time Lyapunov equations (GCALEs) must be solved, giving two full rank gramians. A more efficient version of this technique was presented in [16], approximating the gramians by low-rank matrices. A drawback of these techniques is the fact that one must compute spectral projections onto the left and right deflating subspaces corresponding to the finite eigenvalues of the related pencil [16], [24], which is computationally prohibitive for largescale systems. This paper describes a method to compute low-rank Choleski factorizations of the controllability and observability gramians of large-scale descriptor systems, without having to compute spectral projections. The sparse descriptor system representation, based on the unreduced Jacobian matrix, allows for CPU and memory efficient computations. The low-rank Choleski factors are used to construct reduced-order state-space models. The proposed technique is an improvement of the low-rank Choleski factor alternate direction implicit (LRCF-ADI) described in [19]–[21], and will be called sparse LRCF-ADI. The paper is organized as follows. In Section II, system modeling and existing gramian-based techniques are described, and the motivation for this work is stated. Section III presents the improved method to compute low-rank factorizations of the gramians. Sections IV and V describes applications of the new method to large-scale power system models. Section VI concludes. The following notation is used in this paper: the complex plane is denoted by and the open left half-plane by . Symand denote real and complex mabols trices, respectively. The matrix stands for the transpose of , denotes the complex conjugate and transpose of

0885-8950/$25.00 © 2008 IEEE

FREITAS et al.: GRAMIAN-BASED REDUCTION METHOD APPLIED TO LARGE SPARSE POWER SYSTEM DESCRIPTOR MODELS

, matrix of order

, and . The identity is denoted by or simply .

II. POWER SYSTEM DYNAMIC REDUCTION A linear time invariant (LTI) dynamic power system model can be represented by the state-space system

(1) where is the state vector, the input the output vector; , vector and , and are system matrices. The transfer function or frequency response of system (1) is defined . The state-space representation as in (1) is obtained from the unreduced Jacobian matrix or descriptor form [25]

1259

with is the reduced state vector, where is the output, and , , , and are the reduced system state matrices. The reduced-order model (5) must satisfy a number of requirements [9]: the order of the reduced system must be much ; the approxsmaller than that of the original system (in suitable norm) must be small imation error for a given set of inputs; stability must be preserved; the procedure should be computationally efficient; and there should be a stopping criterion based on an error tolerance. Simulation of the reduced-order model should be (much) faster than of the original model: this may introduce additional requirements on the (sparsity) structure of the reduced-order model, since the costs of computing with dense matrices may easily overshadow the advantage obtained with dimension reduction. There exist infinitely many state-space representations that satisfy (5). The general expressions are of the form

(6) (2)

where is a vector of algebraic variables and , , , , , , , , and are sparse matrices with adequate dimensions. For practical power systems and this type is always a nonsingular matrix [25], [26]. of problem, The state-space (1) and descriptor representations (2) are related as follows:

(3)

and and . In where both matrices and so the sequel, we describe techniques to compute that the reduced-order model satisfies the above requirements. B. Balanced Truncation Method Given a stable state-space system, truncated balanced realization (TBR) produces a stable reduced model, that is usually superior to models obtained by moment matching techniques [13], with a global frequency domain error bound [27]. The TBR method, which is computationally highly intensive, consists of three main steps (see [9], [10], and [28] and references therein for more details). The first step is to compute the controllability and observability gramians and dense , respectively, by solving the Lyapunov equations

The set of (2) can be put in a compact form as

(4)

(7) (8)

where , with ; and have adequate dimensions. The and matrices , , transfer function of system (4) (and (1)) is . The main contribution of this paper is to extend the LRCF-ADI method of [19]–[21] to the sparse descriptor system (2) to obtain significant computational gains in power system model reduction.

Since (8) is the dual of (7), solution details are discussed only for (7). It is also assumed that the system (1) is stable and in that case, the Lyapunov (7) and (8) have a unique positive semidefinite solution [14]. and The second step is to factor and as , where and are the Choleski factors. These factors are then multiplied and their product decomposed as follows:

A. Order Reduction of Systems The problem addressed here is to approximate other dynamical system

with an-

(5)

(9) where and are composed of the leading columns of and , respectively. Matrix con, that are also known as the tains the singular values of Hankel singular values of system (1) [9].

1260

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 23, NO. 3, AUGUST 2008

The third step is to construct a reduced-order model of order via (6), where and are given by [27]

(10) , and are chosen such that the largest Here , Hankel singular values are preserved in the reduced-order model [9], [28]

(11) The frequency response of the reduced system satisfies the global error bound [1], [9], [28]

(12) where and are the frequency responses of and , respectively, and is the frequency, in rad/s. Relation (12) is an a priori error bound, so that given an error tolerance, one can decide how many states to truncate without forming the reduced model. The main drawback of TBR in the above form is that the solution of the Lyapunov equations, the computation of the Cholesky factorizations and singular value decomposition, in order to compute all Hankel singular values, involve dense, and sometimes ill-conditioned, matrix computations of com, making TBR of little practical value to systems plexity having more than 500–1000 states [12], [27]. However, the and often have low-rank compared to , and gramians and as well as the in many cases, the eigenvalues of Hankel singular values , decay very rapidly [21], [22]. such that In other words, in many cases there is a and approaches zero. This leads to the idea of approximating the gramians by low-rank matrices [19]–[21] that will be exploited by the algorithms described in this paper. C. ADI Method for Solving Lyapunov Equations The alternating direction implicit (ADI) method [19], [20], [29] iteratively computes the solution by performing the following steps, beginning with and

, and . is a Cayley spectral transform on It is worth noting that [30], [31]. If the matrix is asymptotically stable and the ADI-parameters have strictly negative real parts, the iterative Lyapunov equation solution (15) converges to the exact solution , within an acceptable tolerance [19], [32]. For convergence, in (15) should have spectral radius less than the matrices unity. The ADI parameters affect the maximum spectral radius and hence have significant impact on (the rate of) convergence. Therefore, the ADI-parameters must be chosen with care; see Section II-D and [19]–[21] for more details and proof of the above. . Each paWe must choose ADI-parameters rameter may be real or complex. Complex parameters should and . This enappear in conjugate pairs, i.e., sures that the approximate solutions can be put in a real form. ADI iterations, the ADI-parameters are If there are , used cyclically [19]. In iteration the shift is then where for . For , is the reis zero, mainder after division by . If the remainder .

where

D. Low-Rank ADI Iterative Method The observation that the final solution of (7) has low-rank for a given acceptable tolerance has lead to a more (memory and CPU) efficient low-rank formulation for the ADI iteration [15], in (15) as [21], [27]. The main idea is to rewrite the iterate the following product of Cholesky factors:

(16) Hence, the low-rank ADI algorithm (LR-ADI) [10], [19] for can be formulated from (15) calculating the Choleski factor

(17) (18) , with and is where the maximum number of iterations. A further optimization, where the number of columns to is kept constant by using a clever reordering of modify the ADI-parameters, was presented in [21] as the low-rank Cholesky factor-ADI method (LRCF-ADI)

(13) (14)

(19)

where , with , is a shift or an ADIparameter. Steps (13) and (14) are equivalent to the single step [10]

(20)

(15)

(21) and the ADI-parameters are used in a where cyclic way, as explained in Section II-C.

FREITAS et al.: GRAMIAN-BASED REDUCTION METHOD APPLIED TO LARGE SPARSE POWER SYSTEM DESCRIPTOR MODELS

The choice of a set of nonconstant ADI-parameters will dominate the rate of convergence to a solution depending on the specassocitral radius of the positive semidefinite matrix ated with the error bound in iterate given by [23]

1261

descriptor systems. This LRCF-ADI variant performs the following operations:

(22) where

is the exact solution of (7), and

is defined as

(23) Appendix A presents further details on the ADI-parameter computation problem. The ADI-parameter set of interest is a function of and should be chosen so as to solve the mini-max problem [21], [23], [29]

(24) where

is a region in the open left half-plane and (25)

Since the solution to (24) is not known when is an arbitrary region, or is the full spectrum of , strategies to compute (near-)optimal parameters are used in practice, based on the dominant eigenmodes of . Here we use the technique proposed in [19] to compute a region based on an approximation of the dominant spectrum (in magnitude) of the matrix and its largest and smallest inverse. The set is generated by the (magnitude) stable Ritz-values (approximate eigenvalues of and , computed by the Arnoldi method). Then ADI-parameters are selected among the elements of based on the mini-max optimization criterion (24). At least elements are necessary in to determine ADI-parameters. This produces a suboptimal solution for (24); see also [19]. For other ways to compute the ADI-parameters , see [20], [21] and references therein. The stopping criterion for LRCF-ADI, checking the (relative) and , can be implemented as [19] difference between

where is the Frobenius norm, and . The tolerance depends on the characteristics of each problem, since the dominant invariant subspace of the gramian solutions and related decay of Hankel singular values may change with the systems and their inputs and outputs. For other stopping criteria, see for instance [19], [21], and [29]. E. Previous Work on LRCF-ADI for Descriptor Systems In [16] a low-rank iterative method is presented for solving projected generalized Lyapunov equations associated with

where can be a singular or nonsingular matrix, , and is the spectral projection onto the left deflating subspace corresponding to the finite eigenvalues of the pencil (assumed to be regular) [16]. For the computation of the observability gramian, also the spectral projection onto the right deflating subspace, , i.e., the eigenspace corresponding to the , is needed. These projecfinite eigenvalues of the pencil tions separate the finite modes from the infinite modes and hence can be seen as a deflation process. This strategy is used to solve projected (generalized) Lyapunov equations [16] for descriptor systems. In principle, the projections could be restricted to the dominant modes only, but this is subject of further research. In some cases these projections are available explicitly (but expensive), due to the structure of the underlying problem [16]. and A better strategy would be to avoid the computation of and directly operate on the descriptor model. Therefore, in Section III we propose a modified LRCF-ADI method for descriptor systems that does not require the computation of and stays as close as possible to LRCF-ADI (19)–(21), while fully exploiting matrix sparsity. III. SPARSE LRCF-ADI METHOD (SLRCF-ADI) The most expensive operations in LRCF-ADI (19)–(21) are the linear system solutions with shifted matrices . For large systems with poor sparsity characteristics, such as state-space representations of power system small-signal stability models [4], [25], the huge computational load represents a severe bottleneck. The Sparse LRCF-ADI method presented in Section III-A solves the same linear system equation, but in a more efficient way by operating on the sparse descriptor system representation. The advantage is that the original LRCF-ADI framework (19)–(21) and its results remain intact, without having to resort to the expensive generalized variants developed in [16]: only the way the linear systems are solved is changed. Section III-B shows how a reduced-order system model can be built from the computed low-rank approximations. A. Sparse LRCF-ADI Via Unreduced Jacobian Model The sparse LRCF-ADI (SLRCF-ADI) performs operations such as and implicitly by using the sparse descriptor matrices of the unreduced Jacobian representation (2) and the relations in (3). is performed as the The operation solution of the linear system

(26)

1262

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 23, NO. 3, AUGUST 2008

Using (3), expression (26) is transformed to

, replace by , In order to compute by , and by a null matrix of appropriate dimensions and perform steps 3–7 of Algorithm 1. (27) B. Reduced-Order System Via Unreduced Jacobian Model

With

, it follows that

(28)

and of the Once the low-rank Choleski factors ADI controllability and observability gramians are obtained, the reduced-order system can be determined by using the low-rank square root method [20]. This technique is a slight modification of the classical square root method, which in turn is a numerically advantageous version of the balanced truncation technique [20]. The complete model order reduction procedure is summarized in Algorithm 2.

Although matrix (28) has a larger dimension , it is highly sparse and can efficiently be solved by any suitable sparse is needed. LU factorization method [4], [25]. Note that only See steps 1–2 of Algorithm 1. . In this A similar approach is followed for by , and by case, solve (28) by replacing by , a null matrix of appropriate dimensions. The complete sparse LRCF-ADI method is shown in Algorithm 1.

Algorithm 2 Reduced-Order System Via SLRCF-ADI.

Algorithm 1 Sparse LRCF-ADI.

OUTPUT: Reduced-order state-space system ( ,

INPUT: Unreduced Jacobian system matrices ( , , , , , ) and ADI-parameters , (real and/or in complex conjugate pairs).

1) Compute the ADI-parameters

OUTPUT: Low-rank Cholesky factors 1) Solve

of

.

INPUT: Unreduced Jacobian system matrices ( , , , , , , , , ) and ADI-parameters , . ,

,

).

.

2) Compute low-rank Choleski factors and of the observability and controllability gramians, respectively, via SLRCF-ADI (Algorithm 1) 3) Compute thin [33] singular value decomposition , with singular values sorted in descending order.

from

4) Truncate after the th largest singular value 2) Set

(cf. (19)),

.

3) for 4) Solve

5) Compute real matrix transformations

where and respectively. 5) Set 6) Set

and

as

from

(cf. (20))

6) Compute

,

are the real forms of ,

and

,

and

(cf. (21))

7) end for When solving for the low-rank approximation of the observis used instead of and operations are perability gramian, formed with where (29) Using (3),

can be solved from

(30)

,

,

, and

. In step 5 of Algorithm 2, a conversion from complex to (or ) real form is applied to each pair of columns of with complex-valued entries, where one column is the con(or ) has jugate of the other (by construction). When columns, a block-diagonal transformation matrix is used,

FREITAS et al.: GRAMIAN-BASED REDUCTION METHOD APPLIED TO LARGE SPARSE POWER SYSTEM DESCRIPTOR MODELS

TABLE I TEST-SYSTEMS DERIVED FROM BIPS98

TABLE II TEST-SYSTEMS DERIVED FROM BIPS07

, with blocks defined as follows: if a column has only real elements; in case of two complex conjugate columns and , the related diagonal block has the form in (31), with

(31) Applying this transformation to

, we have (32)

IV. NUMERICAL RESULTS FOR SISO SYSTEMS The performance of the technique was assessed with the use of several test-systems derived from two practical Brazilian Interconnected Power System (BIPS) models. The first one (BIPS98) relates to a 1998 heavy load condition and has 2380 buses, 2536 nonlinear loads, 3450 ac branches, eight HVDC converters, up to 124 synchronous machines and six FACTS devices. The second model, BIPS07, represents an evolved state of the same power grid and topology for the year 2007. BIPS07 has 191 synchronous machines, 3647 buses, 5175 branches, 3639 nonlinear loads, eight HVDC converters and eight FACTS devices. The state-space and descriptor models were obtained by a production-grade software [34], with a procedure for data conversion to Matlab format, for seven test-systems (whose characteristics are summarized in Tables I and II). The numbers given in Tables I and II relate to synchronous machines having dynamic models with subtransient effects fully represented. The 606- and 1142-state systems have 106 generators with dynamic models and another 18 that were converted to negative loads. A. SISO Test-Systems All SISO results in this paper relate to the same scalar transfer function (TF), whose input is the reference exciter voltage, , for the Itaipu generator, and whose output is the terminal power ( is the same output for the reduced-model). deviation,

1263

Although all test-systems are stable, there are eigenvalues near zero [35], since individual rotor angle deviations, rather than relative values, were taken as state variables. In the case of TBR-based methods, this may lead to numerical problems. To circumvent this, a small shift of value is applied to the real axis of the complex plane (see Appendix B for details), so that all eigenvalues of the spectrally transformed system are located in the left half-plane and TBR via SLRCF-ADI can be applied. The maximum number of ADI iterations to compute the (see (21)). The gramians was set to a constant value and the values of and were fixed ADI-parameters heuristically, the latter parameters ranging from 40 to 180 depending on the test-system. The final dimension of the reduced-order model is system dependent and controlled by a tolerance for the global error bound (see (12)) in Algorithm 2. The gramians were computed for all test-systems and used to find reduced-order models. Frequency domain simulations were carried out for each test-system and corresponding reduced-order model, and compared over a wide frequency range, between 0.01 rad/s and 100 rad/s. Time domain simulations were carried out for a step disturbance of 0.01 pu applied to the TF input, while monitoring the TF output. The obtained results were then used to determine the absolute errors between the full and reduced-order model responses. All computations were carried out using Matlab 6.5 on an Intel Pentium 4 processor with a 2.66-GHz clock and 1 GB of RAM. Identical simulations were carried out for all test-systems in Tables I and II. Since the results are very similar, only results for the 606-state system are discussed in detail; results for the other systems are discussed briefly. The 606-state test-system has several machines modeled without their power system stabilizers (PSSs) and automatic voltage regulators (AVRs), while others of smaller capacity were converted into negative static loads. The speed governors were also not considered (see Table I for details). The system has a higher number of poorly damped modes. In the experiments for the 606-state system, a shift was used to move the eigenvalue closest to origin to the complex plane left-side. Algorithm 1 was used with 70 ADI-parameters, , and . SLRCF-ADI needed 90 iterations for the computation of the gramians, yielding a 39th-order reduced model. Table III shows a summary of these data for the 606-state system and six other systems derived from BIPS98 and BIPS07. The step responses and absolute errors for the 606-state system in Figs. 1 and 2, respectively, indicate a very good match between the original and reduced-order model. The 606-state system and its 39th-reduced-order frequency responses are compared in Fig. 3; the corresponding error is shown in Fig. 4. Both time and frequency simulations show that the reduced-order model produces accurate responses, despite the presence of various poorly damped modes. Fig. 5 shows eigenvalue loci of the 606-state system and its 39th-reduced-order model, together with the ADI-parameters. Fig. 5(a) shows the eigenvalues with largest magnitude, while Fig. 5(b) shows those of smallest magnitude (near the origin). It can be seen from Fig. 5(a) that the ADI-parameters accurately

1264

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 23, NO. 3, AUGUST 2008

TABLE III PARAMETERS USED AND ORDER OF THE REDUCED-ORDER MODELS

Fig. 4. Frequency response deviations between the full and reduced model (k = 39) for the 606-state system.

Fig. 1. Time domain responses for the 606-state system (solid) and 39th-order reduced model (dashed).

Fig. 2. Absolute deviations between time domain responses for the full and reduced model (k = 39) for the 606-state system.

Fig. 5. Poles and ADI-parameters for the 606-state test system.

Fig. 3. Frequency responses for the 606-state system (solid) and 39th-order reduced model (dashed).

approximate the largest eigenvalues of the complete model. On the other hand, Fig. 5(b) shows that the ADI parameters do not

approximate the smallest eigenvalues very well: this can be explained by the fact that the smallest eigenvalues are generally more difficult to compute with iterative methods such as the Arnoldi method. Similar tests are performed for the other test-systems, considering their characteristics in Table III. The absolute errors in the frequency response of the 30th-order reduced model for the 1142-state system, and of the 28th-order reduced model for

FREITAS et al.: GRAMIAN-BASED REDUCTION METHOD APPLIED TO LARGE SPARSE POWER SYSTEM DESCRIPTOR MODELS

1265

Fig. 6. Absolute error in frequency response of the 30-state reduced model for the 1142-state system.

Fig. 9. Absolute errors in the frequency responses of the 29-state reduced-order models for the 3078-state system, computed by SLRCF-ADI (solid) and PVL (dashed), respectively.

Fig. 7. Absolute error in frequency response of the 28-state reduced model for the 1450-state system.

Although similar tests were carried out for the 1998-state, 2476-state and 3078-state systems, only the latter is described, since the simulations, parameters, and results for the other systems are similar. In order to compare the SLRCF-ADI with other system reduction methods, we also include results obtained by the Padé-via-Lanczos (PVL) method [36], [37] (see Appendix C for additional details). For reducing the 3078-state system, we used a shift and , 45 ADI-parameters, and . The low-rank gramians had dimension 3078 180, and the reduced-order model obtained by SLRCF-ADI had order 29. The total computational time was 168 s (24 s for computing the ADI parameters, 130 s for computing the gramians, and 14 s for constructing the reduced-order model). A model of the same order was computed using PVL (expansion point rad/s), requiring only 13 s. This model, however, needed to be stabilized by heuristically removing its spurious unstable poles. The frequency response plots in Fig. 9 show that the reducedorder model calculated by SLRCF-ADI is more accurate than the reduced-order model computed by PVL. The quality of the PVL model is known to be very sensitive to the expansion point PVL was not able [38]. Also for different expansion points to compute a stable reduced-order model. SLRCF-ADI, on the other hand, is more robust and preserves stability. B. Comparing the Performances of LRCF-ADI and SLRCF-ADI

Fig. 8. Absolute error in frequency response of the 32-state reduced model for the 1693-state system.

the 1450-state system, are shown in Figs. 6 and 7, respectively. Fig. 8 shows the absolute errors in the frequency response of the 32th-order reduced model for the 1693-state system. Several other systems were derived from the BIPS07 model by keeping/removing the dynamic representation of some machine controllers, leading to different number of states and algebraic variables for each test-system as described in Table II. The 3078-state system is a complete model of the actual Brazilian interconnected power grid, and the 2476- and 1998-state systems were derived from it by neglecting many of its original controllers (see Table II).

The gramians were computed operating on both the sparse Jacobian and the dense state-space matrices. Fig. 10 shows CPU timings for the computation of the controllability and observability gramians by LRCF-ADI applied to state-space systems (solid) and SLRCF-ADI applied to sparse descriptor systems (dashed). The plots were generated considering 100, 200, and 300 ADI iterations. It can be seen that the SLRCF-ADI approach is faster in all cases; for state-space systems with more than 1500 states, the state-space approach even failed (due to memory limitations). In general, the more ADI iterations, the better the accuracy of the low-rank approximation, as can be seen in Fig. 11. Note that the first 60 singular values of the low-rank approximation, after 300 iterations, become very close to the singular values of the original system. Concerning the computation of

1266

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 23, NO. 3, AUGUST 2008

TABLE IV CPU TIMINGS (S) FOR TIME DOMAIN SIMULATIONS OF TEST-SYSTEMS AND REDUCED-ORDER MODELS OBTAINED BY THE SLRCF-ADI PROCESS AND FOR COMPLETE REDUCTION PROCESS

Fig. 10. CPU timings for computation of the gramians by LRCF-ADI (solid) and SLRCF-ADI (dashed).

TABLE V PARAMETERS USED AND ORDER OF THE REDUCED-ORDER MODELS FOR THE MIMO SYSTEMS

V. NUMERICAL RESULTS FOR MIMO SYSTEMS This section presents results for the 3078- and 606-state MIMO systems. Table V shows parameters used to compute the reduced-order models with SLRCF-ADI. A. Well-Damped MIMO Test-System Fig. 11. Singular values of low-rank approximations of the 606-state system gramians become more accurate as the number of ADI iterations increases.

only the ADI-parameters, LRCF-ADI is faster for systems with fewer than 2000 states. In Fig. 10, a linear increase in CPU time can be observed for SLRCF-ADI (operating on sparse Jacobian matrices), against a quadratic increase for LRCF-ADI (operating on dense statespace matrices). This difference is due to the fact that sparse LU-factorization techniques are used for the Jacobian model, while the state-space model requires techniques for dense ma. trices Table IV shows CPU timings for time domain simulations (CPU sim.) of the complete systems and reduced-order models computed by SLRCF-ADI. The CPU timings for computing the reduced-order models (CPU ROM) are also shown. Clearly, time domain simulation of the reduced-order models is much faster in all cases. For these examples, computation of the reduced-order models is slower than time simulation of the test-systems. However, when time domain simulations over larger time intervals are needed, possibly including smaller time steps as well, CPU timings for reduced-order model computation and simulation will be factors lower than CPU timings for full system simulation. Furthermore, the reduced-order model can be (re)used in applications where full system simulation may be expensive, such as in large-scale power system controller design. [26].

This subsection describes results for the 3078-state system. The reference voltage signals from Itaipu, Xingo, Jacui and Tucurui power stations comprise the input control vector. The terminal power deviations for the same four generators define the output matrix. These four generators are located at different geographic areas of the Brazilian interconnected power system, which are very distant from each other. CPU timings for computing the ADI parameters, the gramians and the reduced model were 24.7 s, 140.21 s and 56.67 s, respectively. A convergence was used for SLRCF-ADI. tolerance of Figs. 12 and 13, respectively, show time domain simulations for two scalar transfer functions that have the best and worst matches between full and reduced models, out of the 16 input–output combinations. The best match is shown in Fig. 12, where the input and output relate to the same power station (Itaipu). In the worst case, the input is at Tucurui, located in the North, and the output at Jacui in the South. These two power plants are more than 3000 km apart, with various other power plants in between. Roughly 1 s delay at the beginning of the step response in Fig. 13 confirms this fact. Despite the huge distance between the sites, the MIMO reduced-order model is very accurate. Fig. 14 shows that the sigma plots [39] of the transfer function of the complete system and the transfer matrix of the reduced-order model function matrix virtually coincide, and confirm that the 106-state reduced-order model is accurate.

FREITAS et al.: GRAMIAN-BASED REDUCTION METHOD APPLIED TO LARGE SPARSE POWER SYSTEM DESCRIPTOR MODELS

Fig. 12. Time domain responses for 3078-state system (solid) and 106-state reduced-order model (dashed), for input and output at Itaipu power station (best match).

1267

Fig. 15. Sigma plots (maximum and minimum singular value) of 606-state system and 92-state reduced-order model.

same type of input/output signals as for the 3078-state MIMO test-system). Fig. 15 shows the sigma plot for the test-system and the 92-state reduced-order model computed by the SLRCF-ADI , where there method. Except for frequencies below (most likely caused by rounding errors, are deviations in considering the small magnitude), the reduced-order model is very accurate. VI. CONCLUSIONS

Fig. 13. Time domain responses for 3078-state system (solid) and 106-state reduced-order model (dashed), for input at Tucurui and output at Jacui power stations (worst match).

Fig. 14. Sigma plot (maximum and minimum singular value) of 3078-state system and 106-state reduced-order model.

B. Poorly Damped MIMO Test-System As a final example we discuss the performance of the SLRCF-ADI technique for the 606-state system with five inputs and five outputs at the Itaipu, Jacui, I. Solteira, G. B. Munhoz and Promissao generators of BIPS98 (using the

The sparse LRCF-ADI method is an improvement on the low-rank Cholesky factorization—ADI method (LRCF-ADI), designed for the fast computation of the gramians of large-scale descriptor systems of power system dynamic models. As the method exploits the sparsity of the descriptor systems, it is applicable to large-scale systems while keeping memory and CPU requirements at modest levels. Since the improvement proposed in this paper only involves the way linear systems are solved, all results of the original LRCF-ADI method are preserved. The gramians were always used implicitly to compute reduced-order models for power systems, so they never had to be explicitly factored or stored. Results included time simulations, frequency responses, eigenvalue loci and the absolute errors between the full and reduced system responses. They also showed superior performance on CPU times and memory requirements for the SLRCF-ADI approach, operating on the unreduced Jacobian (descriptor) systems, over the conventional LRCF-ADI approach, operating on the dense state-space systems, while maintaining the same accuracy. The method showed high accuracy and effectiveness when applied to large-scale MIMO systems as well. APPENDIX A ADI-PARAMETER COMPUTATION PROBLEM In this section it is assumed that all eigenvalues of have negative real part. For more details on the results in this section, see [19], [21], and [23]. At any iteration , when (15) is used as an approximated solution of (7), the error bound related to the exact solution is given by (22).

1268

From (15), when

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 23, NO. 3, AUGUST 2008

, we have the error bound [23] (33)

Once , then , same way, at iteration Therefore, in iteration , we have

in (33). In the .

where matrices between parentheses commute with each other. Then it is possible to establish an upper bound for based on (40). Suppose that is diagonal, where is a diagonal izable in such way that and is the related right matrix with the eigenvalues of eigenvector matrix. Then, (40) can be put as

(41) with . We can find an upper bound for as (34) Defining the rational functions

(35)

where , with tion number of . Equation (42) is equivalent to

(42) the condi-

we can associate them as mappings of the left half-plane into the unit disk. Note that (34) has this characteristic. From (15) (43) (36) where . Then, . An upper bound for the decay-rate of the eigenvalues of can be evaluated by the Schmidt–Mirsky theorem [23]. For the th eigenvalue of , this relation is given by

(37) Then, for the eigenvalue equality can be obtained [23]

of

the following in-

(38)

which is valid for any choice of . in such a The idea is now to choose the ADI parameters way that the bound in (43) is as tight as possible, since this would imply an accurate low rank approximation of . This leads to the mini-max problem

(44) Since the exact eigenvalues are usually not known for large scale state- or descriptor matrices, one has to use approximate eigenvalues instead, as is described in Section II-D.

THE

APPENDIX B -SHIFT APPROACH

The -shift approach used in this paper can be interpreted as a change in the Laplace domain variable . Let us consider , where is a Laplace variable and is a strictly positive constant. Then, for the system (1), we have

The result (38) can be written as

(39) The expression is equal to the square root of the largest eigenvalue of the positive semidefinite matrix [23]

(40)

Applying the change of variable , assigning new variables, and inverse Laplace transforms leads to

(45)

FREITAS et al.: GRAMIAN-BASED REDUCTION METHOD APPLIED TO LARGE SPARSE POWER SYSTEM DESCRIPTOR MODELS

where , with and the Laplace transform operator and its inverse, respectively; ; and . Using low-rank gramians of (45) computed by (S)LRCFADI, a reduced-order model can be constructed following (6)

(46) where

. It follows that . Since , where is the order of the reduced system, the reduced-order state-space matrix of the original system is (47) Back-substitution of of the original system (44)

gives the reduced-order model

(48) Note that matrices , and puted for the -transformed system.

are equal to those com-

APPENDIX C KRYLOV SUBSPACE-BASED REDUCED-ORDER MODEL The transfer function of system (4) is (49) Let

be an expansion point such that the LU factorization is nonsingular. Then can be written as (50)

where , and . of can be obtained by An th Padé approximant running steps of the Lanczos process. This procedure, called Padé-via-Lanczos (PVL) [36], [37], is applied to the matrix , with and as the right and left vector, respectively. As result, a is obtained and a reduced-order tridiagonal Lanczos matrix model is computed as follows.

(51) where is chosen such that some appropriate stopping criterion is an unitary vector whose first element is equal is satisfied, is a normalizing scalar. to one, and ACKNOWLEDGMENT The authors would like to thank the Brazilian Electrical Energy Research Center (CEPEL) for the power system software and data as well as the reviewers for their valuable suggestions.

1269

REFERENCES [1] J. Sanchez-Gasca and J. H. Chow, “Power system reduction to simplify the design of damping controllers for interarea oscillations,” IEEE Trans. Power Syst., vol. 11, no. 3, pp. 1342–1349, Aug. 1996. [2] N. Martins, F. G. Silva, and P. C. Pellanda, “Utilizing transfer function modal equivalents of low-order for the design of power oscillation damping controllers in large power systems,” in Proc. IEEE Power Eng. Soc. Meeting, 2005, pp. 2642–2648. [3] J. Sanchez-Gasca and J. H. Chow, “Computation of power system low order models from time domain simulation using a Hankel matrix,” IEEE Trans. Power Syst., vol. 12, no. 4, pp. 1461–1467, Nov. 1997. [4] N. Martins, L. T. G. Lima, and H. J. C. P. Pinto, “Computing dominant poles of power system transfer functions,” IEEE Trans. Power Syst., vol. 11, no. 1, pp. 162–170, Feb. 1996. [5] J. R. Smith, J. F. Hauer, D. J. Trudnowski, F. Fatehi, and C. S. Woods, “Transfer function identification in power system application,” IEEE Trans. Power Syst., vol. 8, no. 3, pp. 1282–1290, Aug. 1993. [6] G. Troullinos, J. Dorsey, H. Wong, and J. Myers, “Reducing the order of very large power system models,” IEEE Trans. Power Syst., vol. 3, no. 1, pp. 127–133, Feb. 1988. [7] D. Chaniotis and M. A. Pai, “Model reduction in power systems using Krylov subspace methods,” IEEE Trans. Power Syst., vol. 20, no. 2, pp. 888–894, May 2005. [8] J. Rommes and N. Martins, “Efficient computation of transfer function dominant poles using subspace acceleration,” IEEE Trans. Power Syst., vol. 21, no. 3, pp. 1216–1226, Aug. 2006. [9] A. C. Antoulas, Approximation of Large-Scale Systems. Philadelphia, PA: SIAM, 2005. [10] P. Benner, Numerical Linear Algebra for Model Reduction in Control and Simulation, 2005. [Online]. Available: http://www-user.tu-chemnitz.de/benner/pub/. [11] U. Baur, “Control-oriented model reduction for parabolic systems,” Ph.D. dissertation, Tech. Univ. Berlin, Berlin, Germany, 2008. [12] A. C. Antoulas and D. C. Sorensen, “Approximation of large-scale systems: an overview,” Int. J. Appl. Math. Comput. Sci, vol. 11, no. 5, pp. 1093–1121, 2001. [13] A. C. Antoulas, D. C. Sorensen, and S. Gugercin, “A survey of model reduction methods for large scale systems,” in Contemporary Mathematics. Boston, MA: AMS, 2005. [14] A. C. Antoulas, D. C. Sorensen, and Y. Zhou, “Solving large-scale control problems,” IEEE Control Syst. Mag., vol. 14, no. 1, pp. 44–59, Feb. 2004. [15] D. Lukashevich, F. Coccetti, and P. Russer, “Identification and model order reduction for TLM analysis,” in Numerical Modelling. New York: Wiley InterScience, 2006, pp. 75–92. [16] T. Stykel, Low-Rank Iterative Methods for Projected Generalized Lyapunov Equations, DFG Research Center MATHEON, Tech. Univ. Berlin, Berlin, Germany, 2006, Tech. Rep. Preprint 198. [17] T. Stykel, “Gramian based model reduction for descriptor systems,” Math. Control, Signals, Syst., vol. 16, pp. 297–319, 2004. [18] T. Stykel, “Analysis and numerical solution of generalized Lyapunov equations,” Ph.D. dissertation, Tech. Univ. Berlin, Berlin, Germany, 2002. [19] T. Penzl, “A cyclic low-rank Smith method for large-scale Lyapunov equations,” SIAM J. Sci. Comp., vol. 2, no. 4, pp. 1401–1418, 2000. [20] T. Penzl, Lyapack—A Matlab Toolbox for Large Lyapunov and Riccati Equations, Model Reduction Problems, and Linear-quadratic Optimal Control Problems. Users’ Guide (Version 1.0), 1999. [Online]. Available: http://www.tu-chemnitz.de/sfb393/lyapack. [21] J. R. Li and J. White, “Low-rank solution of Lyapunov equations,” SIAM J. Matrix Anal. Appl., vol. 24, pp. 260–280, 2002. [22] A. C. Antoulas, D. C. Sorensen, and Y. Zhou, “On the decay rate of Hankel singular values and related issues,” Syst. Control Lett., vol. 46, no. 5, pp. 323–342, 2002. [23] D. C. Sorensen and Y. Zhou, Bounds on Eigenvalue Decay Rates and Sensitivity of Solutions to Lyapunov Equations, Rice Univ., Houston, TX, 2002, Tech. Rep. TR02-07. [24] Y. Saad, Numerical Methods for Large Eigenvalue Problems. Manchester, U. K.: Manchester Univ. Press, 1992. [25] N. Martins, “Efficient eigenvalue and frequency response methods applied to power system small-signal stability studies,” IEEE Trans. Power Syst., vol. 1, no. 1, pp. 217–226, Feb. 1986. [26] F. D. Freitas and A. J. A. Simoes-Costa, “Computationally efficient optimal control methods applied to power systems,” IEEE Trans. Power Syst., vol. 14, no. 3, pp. 1036–1045, Aug. 1999.

1270

[27] J. R. Li, “Model reduction of large linear systems via low rank system gramians,” Ph.D. dissertation, Massachusetts Inst. Technol., Dept. Math., Cambridge, 2000. [28] K. Glover, “All optimal Hankel-norm approximations of linear multivariable systems and their L -error bounds,” Int. J. Control, vol. 39, pp. 1115–1193, 1984. [29] E. L. Wachspress, “Iterative solution of the Lyapunov matrix equation,” Appl. Math. Lett, vol. 107, no. 1, pp. 87–90, 1988. [30] L. T. G. Lima, L. H. Bezerra, C. Tomei, and N. Martins, “New methods for fast small-signal stability assessment of large scale power systems,” IEEE Trans. Power Syst., vol. 10, no. 4, pp. 1979–1985, Nov. 1995. [31] G. Angelidis and A. Semlyen, “Efficient calculation of critical eigenvalue clusters in the small-signal stability analysis of large power systems,” IEEE Trans. Power Syst., vol. 10, no. 1, pp. 427–432, Feb. 1995. [32] Z. Gajic and M. T. J. Qureshi, Lyapunov Matrix Equation in System Stability and Control. San Diego, CA: Academic, 1995. [33] G. H. Golub and C. F. van Loan, Matrix Computations , 3rd ed. Baltimore, MD: John Hopkins Univ. Press, 1996. [34] PacDyn® CEPEL, PacDyn User’s Manual, CEPEL—Brazilian Electrical Energy Research Center, 2008, Tech. Rep. Version 7.0. [35] P. Kundur, Power System Control and Stability . New York: McGraw-Hill, 1994. [36] R. W. Freund and P. Feldmann, “Small-signal analysis and sensitivity computations with the PVL algorithm,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 43, no. 8, pp. 577–585, Aug. 1996. [37] R. W. Freund, Krylov-Subspace Methods for Reduced-Order Modeling in Circuit Simulation, Bell Lab., Murray Hill, NJ, 1999, Tech. Rep. Preprint-99-3-20. [38] E. J. Grimme, D. C. Sorensen, and P. M. van Dooren, “Model reduction of state space systems via implicitly restarted Lanczos method,” Numer. Algor., vol. 12, pp. 1–31, 1996.

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 23, NO. 3, AUGUST 2008

[39] J. M. Maciejowski, Multivariable Feedback Design. Reading, MA: Addison-Wesley, 1989.

Francisco Damasceno Freitas (M’94) received the B.Sc. and M.Sc. degrees in electrical engineering from University of Brasilia, Brasilia, Brazil, in 1985 and 1987, respectively. In 1995, he received the Ph.D. degree in electrical engineering from the Federal University of Santa Catarina, Florianopolis, Brazil. Since 1986, he is with University of Brasilia, where he is an Assistant Professor. His area of interest is power system dynamic performance analysis and control and numerical techniques applied to power systems.

Joost Rommes received the M.Sc. degree in computational science, and M.Sc. degree in computer science, and the Ph.D. degree in mathematics from Utrecht University, Utrecht, The Netherlands, in 2002, 2003, and 2007, respectively. Dr. Rommes is currently a researcher at NXP Semiconductors, Eindhoven, The Netherlands, and works on model order reduction.

Nelson Martins (SM’91–F’98) received the B.Sc. degree in electrical engineering from University of Brasilia, Brasilia, Brazil, in 1972 and the M.Sc and Ph.D. degrees from the University of Manchester Institute of Science and Technology, Manchester, U.K., in 1974 and 1978, respectively. He has been with CEPEL, the Brazilian Electric Power Research Center, Rio de Janeiro, since 1978, developing methods and computer tools for power system dynamics and control.

Gramian-Based Reduction Method Applied to Large ...

and observability gramians of large sparse power system de- scriptor models. The method .... transfer function or frequency response of system (1) is defined as.

358KB Sizes 3 Downloads 167 Views

Recommend Documents

Artificial Intelligence Techniques Applied to Reduction ...
Sep 22, 2007 - uncertainty in decision analysis and improve the decision system's .... of the decision situation that have some predictive significance. Since the ...

pdf-1458\the-montessori-method-scientific-pedagogy-as-applied-to ...
... apps below to open or edit this item. pdf-1458\the-montessori-method-scientific-pedagogy-as ... ation-in-the-childrens-houses-by-maria-montessori.pdf.

Conductance method applied to high mobility ...
3. 4. 5. 6. 0.5. 0.6. 0.7. 0.8. 0.9. 1. 1. /f. D. [ ] [ ] p p f p f p. G. G ω ω. 5. * For details see Nicollian-Brews. A: Capacitor area ω: Pulsation (2πf) q: Electron charge f p. 5f p ... Weak inversion responses, due to interactions of minority

A Network Traffic Reduction Method for Cooperative ...
[11], our work in [9] operates across the physical layer and the medium access layer. In this paper, we extend our previous work by developing a network traffic ...

A Network Traffic Reduction Method for Cooperative ...
Wireless positioning has been providing location-based ser- vices in ... Let us consider a wireless network with two types of ..... Cambridge University Press,.

A POD Projection Method for Large-Scale Algebraic ...
based method in a projection framework. ...... continuous-time algebraic Riccati equation. Electron. Trans. Numer. Anal., 33:53–62, 2009. 1,. 3, 3, 6, 6, 7.

A Re-quantization Noise Reduction Method in MPEG-2 ...
the 4x4 integer DCT of H.264 has also been reported[5]. We have focused our attention on the fact that the com- plete re-use of encoding information, such as picture type, motion vector, and macro-block type, makes it possible to suppress re-quantiza

17-08-022. Disaster Risk Reduction Reduction and Management ...
17-08-022. Disaster Risk Reduction Reduction and Management Program.pdf. 17-08-022. Disaster Risk Reduction Reduction and Management Program.pdf.

Transferred Dimensionality Reduction
propose an algorithm named Transferred Discriminative Analysis to tackle this problem. It uses clustering ... cannot work, as the labeled and unlabeled data are from different classes. This is a more ... Projection Direction. (g) PCA+K-means(bigger s

Reduction of Bit Errors Due to Intertrack Interference ...
sity data storage beyond 10 Tera bits/inch . Practical imple- mentations require ... A MIMO channel model is commonly used for wireless digital communications.

Approximate Reduction from AUC Maximization to 1 ...
Our hypothesis class F is the set of convex combination of ranking functions in H, that is, ... Since the problem is a linear program, by duality, we have ρ. ∗ − 1.

An Overview of Peak-to-Average Power Ratio Reduction ... - IJRIT
IJRIT International Journal of Research in Information Technology, Volume 3, Issue 3, ... BER reduction, and their advantages and disadvantages in detail. ... in communications which can be used in both wired and wireless environments.

Photo reduction of CO2 to methanol via TiO2 ...
Page 2 ..... [8] G. R. R. A. Kumara, F. M. Sultanbawa, V. P. S. Perera,. I. R. M. Kottegoda, and K. Tennakone, Continuous flow photochemical reactor for solar ...

Reduction of time to first fix in an SATPS receiver
Sep 15, 2000 - Use of faster-than-real-time signal correlators, together. With GPS .... signal tracking system that rapidly recaptures a lost signal by repeatedly ...

Counseling and Harm Reduction Centers for Vulnerable Women to ...
Dept. of Public Health, School of Public Health, Shahid Beheshti University of Medical ... have been considered as important bridging population for driving HIV.

Issues related to combining risk factor reduction and ...
Address correspondence to C. Barr Taylor, MD, Professor, Department of ..... Winzelberg AJ, Taylor CB, Altman TM, Eldredge KL, Dev P, Constantinou PS.

Reduction of time to first fix in an SATPS receiver
Sep 15, 2000 - adjustable frequency and is ?ltered to produce an audio beat frequency ..... frequency generator (“clock”) is used to provide a time base for an.

A CONTINUATION METHOD TO SOLVE POLYNOMIAL SYSTEMS ...
the path of pairs (/t,7t), where /t,t ∈ [0,T] is a polynomial system and /t(7t) = 0. He proved ... namely H(d) is the vector space of systems of n homogeneous polyno-.