IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

1905

Reduced-Order Transfer Matrices From RLC Network Descriptor Models of Electric Power Grids Francisco Damasceno Freitas, Senior Member, IEEE, Nelson Martins, Fellow, IEEE, Sergio Luis Varricchio, Senior Member, IEEE, Joost Rommes, and Franklin C. Véliz

Abstract—This paper compares the computational performances of four model order reduction methods applied to large-scale electric power RLC networks transfer functions with many resonant peaks. Two of these methods require the state-space or descriptor model of the system, while the third requires only its frequency response data. The fourth method is proposed in this paper, being a combination of two of the previous methods. The methods were assessed for their ability to reduce eight test systems, either of the single-input single-output (SISO) or multiple-input multiple-output (MIMO) type. The results indicate that the reduced models obtained, of much smaller dimension, reproduce the dynamic behaviors of the original test systems over an ample range of frequencies with high accuracy. Index Terms—Descriptor systems, dominant poles, dominant subspaces, eigenvalues, frequency response, index-2 DAE, large-scale systems, low rank Gramians, passivity, reduced models, resonant peaks, RLC circuits, singular systems, transfer function, vector fitting.

I. INTRODUCTION

M

ODEL order reduction (MOR) has been widely used in scalar and multivariable control system applications [1], vibration analysis of large mechanical structures, and VLSI circuit design [2]. Network equivalents are used in power system electromagnetic transient studies [3], as well as in real-time simulators [4] and harmonic distortion analysis [5]. Power system models requiring detailed modeling include those having a large number of resonant peaks in their frequency responses, over a wide range of frequencies [6]. Detailed representations of every component in a large-scale system yield more accurate results but require excessive CPU time, calling therefore for the further development of model order reduction techniques. A good reduced-order model (ROM) should have a reduced state-space vector and be able to reproduce the simulated inputManuscript received December 31, 2009; revised January 13, 2010, October 05, 2010, and January 12, 2011; accepted March 02, 2011. Date of publication May 12, 2011; date of current version October 21, 2011. The work of F. D. Freitas was supported in part by FINATEC. Paper no. TPWRS-01021-2009. F. D. Freitas is with the Department of Electrical Engineering, University of Brasilia, Brasilia, DF, CEP:70910-900, Brazil (e-mail: [email protected]). N. Martins and S. L. Varricchio are with CEPEL, Rio de Janeiro, RJ, CEP:20001-970, Brazil (e-mail: [email protected]; [email protected]). J. Rommes is with NXP Semiconductors, Central R&D, High Tech Campus 46, 5656 AE Eindhoven, The Netherlands (e-mail: [email protected]). F. C. Véliz is with the Pontific Catholic University (PUC), Rio de Janeiro, RJ, Brazil (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TPWRS.2011.2136442

output response of the original system with the desired accuracy while requiring a considerably smaller computational effort. Also, the CPU time needed to generate the reduced-order model should not be excessive, so that the whole ROM utilization process is cost-effective. The main contributions of this paper are: 1) the performances of three existing MOR methods are compared for large power system RLC models; 2) a combination of two methods is proposed, generating a new MOR method, which is also tested against the previous three methods; 3) a simple and effective choice for the alternating direction implicit (ADI) parameters of the low-rank Choleski factor (LRCF) method is presented, when applied to RLC circuits; 4) tests including time and frequency responses as well as spectral (eigenvalue) plots are shown; 5) full data on a 34-bus subtransmission system from practice, together with eight related descriptor system models, are provided in this paper and online, for data exchange in the electromagnetic transients as well as the MOR communities. The first method is the LRCF method, described in [7] and [8]. The second method is based on the computation of dominant subspaces (Subspace Accelerated Dominant Pole Algorithm and Subspace Accelerated MIMO Dominant Pole Algorithm (SADPA [9] and SAMDP [10], respectively)). The third method is vector fitting (VF) [11]–[13], which allows the computation of the reduced model directly from the (finely enough) discretized values of the frequency response of the original system. The fourth method is a new hybrid method, which combines the use of the LRCF method with the VF method, yielding both computational and qualitative gains, allowing transformation of nonpassive ROMs obtained by LRCF into passive ROMs. All transmission lines considered in this paper are modeled by ladder networks, comprised of cascaded RLC PI-circuits, having fixed parameters. Representation of transmission lines by large numbers of RLC PI-circuits results in many resonant peaks in the system’s frequency response, spread over a large frequency range. The formulation and studies of frequency-dependent parameter transmission lines can be found in several classical references [14]–[16] and selected reprints [17]. SISO modal equivalents for -domain models of power systems having long transmission lines with frequency-dependent parameters are described in [18]. Multi-port frequency-dependent equivalents, involving some approximations but given in terms of a passive RLC circuit, which are frequently used in practical EMTP studies, are described in [3]. The terms multi-port and MIMO are used interchangeably in this paper. The frequency dependence of longitudinal parameters is also considered in

0885-8950/$26.00 © 2011 IEEE

1906

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

[19]–[21]. In [19], the frequency dependence is taken into account through the use of synthetic circuits, with a cascade of PI-circuits for each mode of interest. In [20], rational approximation of frequency-dependent admittance matrices has been proposed. A descriptor system representation is presented in [21]. The frequency responses of these detailed frequency-dependent line models can be matched by synthesized RLC circuits of higher dimension and containing only fixed RLC elements. Therefore, high-fidelity transfer function models of frequency-dependent transmission lines can be described by high-dimension linear descriptor systems. This fact is of practical importance and was the major drive for the development and testing of MOR techniques, by the authors, for large-scale linear descriptor systems. This work is organized as follows. Section II reviews the three existing MOR methods and the proposed Hybrid method, whose performances are compared in this paper. Section III describes the test systems, while Sections IV and V present MOR results for the single-input single-output (SISO) and multiple-input multiple-output (MIMO) test systems, respectively. Section VI concludes.

with the structures

. The blocks

,

,

,

, and

have

Regarding the structure of , the variables in can be written as , where and , with . When modeling the RLC test systems of Section III, as well as the two-node RLC circuit of Appendix A, with the formuis singular, despite being lation of this paper, matrix nonsingular. This means that just part of the algebraic variables, part, can undergo Gaussian elimination: the (4) is given by (4), Since algebraic equations can be eliminated. The remaining braic equations and the differential equations are

II. METHODOLOGY

alge-

(5)

A. RLC Descriptor Systems of Index-2 and Their Reduction to Index-0

where the reduced matrix blocks are given by

The MOR methods compared in this work either make integral use of the descriptor formulations (generalized state-space) of the test systems or are based on the identification of reduced models from discretized frequency response data of the test systems. These MOR methods apply to linear, time invariant systems comprising differential and algebraic equations (DAEs) of the form (1) is the generalized state vector, where is the input vector and is the output vector, , with a nonzero matrix, , , and . The generalized state vector comand algebraic variables , prises states respectively. With this notation, matrices , , , and can be written as

(2)

and . equations In the RLC application of this paper, the last in the DAE (1) do not contain any input , which implies in . Also, the network equations in the DAE are organized and are such that the products zero, although each matrix individually is not necessarily zero (see an example of two-node circuit in Appendix A illustrating how all matrices are formed). As a result, the matrix blocks and in (5) are both equal to zero. and , an Considering (5) and the characteristics of and index-2 DAE system [22] can be formed, with states remaining algebraic variables , as follows:

(3)

(6)

where , possibly singular, contains at least one nonzero element per line. In this work, we assume that is square and nonsingular and matrix consists of blocks , , , and ,

In (6), the block is as described in (2) and has . The equations in (6), therefore, constitute an generalized states, index-2 DAE system having states and algebraic variables. If all the algebraic being

FREITAS et al.: REDUCED-ORDER TRANSFER MATRICES

1907

equations in (1) could undergo Gaussian elimination, and no redundant states existed in the DAE, then the descriptor system would be index-1, allowing the use of the sparse low-rank Choleski factor (SLRCF) algorithm [23]. However, considering the index-2 nature of the original systems, and the modest number of algebraic equations and variables in (6), these can all be eliminated by employing a symbolic math procedure (see Appendix C for further details) to form an index-0 (just differential equations) system. This system representation permits the use of the LRCF method [7]. The transformation is possible states in the set of equations can be expressed as a because states (the effective states). linear combination of Appendix C describes the elimination of the algebraic equa. Also, the reduntions, corresponding to the block matrix dant differential equations are eliminated, which then allows the use of LRCF [7]. Eliminations invariably lead to some matrix fill-in, which is however acceptable for the sparse RLC sysis much tems being studied here, since . Symbolic math elimination of the smaller than the rank of algebraic variables and redundant states in (6) yields the index-0 descriptor system model: (7) where . Matrices , , , , and have appropriate dimensions and is now nonsingular, allowing therefore the use of LRCF [7], [8]. , Frequency response results are obtained by setting over a discretized frequency where is the frequency in interval, and evaluating (8) where is the transfer function matrix, whose elements are scalar transfer functions (TF). The symbols and denote the Laplace transforms of variables and , respectively. As shown in [23], frequency response as well as SADPA or SAMDP computations can be performed either on the index-2 descriptor model or on the state-space model in (7). B. Model Order Reduction The problem addressed here is to approximate (1) or (7) with another dynamical system (9) where vector,

with is the reduced order of the state is the output, and , , , and are the reduced model state ma-

trices. In this section, the four methods for model order reduction are briefly described: low-rank Choleski factor (LRCF) ADI [7], [8], Subspace Accelerated Dominant Pole Algorithm (SADPA) [9]/Subspace Accelerated MIMO Dominant Pole Algorithm (SAMDP) [10], vector fitting [11], and a Hybrid method that combines vector fitting with LRCF. 1) Low-Rank Choleski Factor (LRCF) ADI: The LRCF method has a direct analogy with the conventional state-space

truncated balanced realization (TBR) method [1], being however designed to operate on the sparse, large-dimension state-space matrices [7]. This method computes the controllability and observability Gramians, both of which are decomposed into low rank Choleski factors (CF). The product of these low rank CF yields the approximate Hankel matrix, whose singular value decomposition (SVD) will determine the appropriate dimension of the reduced system for a given model error tolerance. The index-2 DAE structure of the test system models does not allow the direct use of the SLRCF method to (1), since this method, as described in [23], can only be applied to models is nonsingular (index-1 DAE where the matrix-block system). 2) Vector Fitting: The VF method identifies the scalar transfer functions from their frequency responses, once some user-specified parameters such as model order and set of initial estimates for the poles are provided. The VF method requires the prior computation of the frequency response [sequence of numerical values for the scalar transfer function (8)], for . The level sufficiently discretized values of frequency of discretization can be varied along the frequency range of interest, so as to allow the identification of each scalar transfer function to the user-specified accuracy. Therefore, the sampling rate should be higher along those frequency ranges showing large numbers of resonant peaks. In this paper, we use an implementation of VF [11] with all the improvements described in [12] and [13]. 3) SADPA and SAMDP: The SADPA method [9] is used for the reduction of SISO systems, and directly operates on the formulation (1), by computing the dominant poles, associated eigenvectors, and transfer function residues. The SAMDP method [10], of the same family of methods, is used for the reduction of MIMO systems. 4) Hybrid Method-VF+LRCF: VF requires initial estimates for the poles, and convergence of VF is improved by providing a fairly accurate set of initial estimates. Hence, we propose the Hybrid method: VF+LRCF, where initial pole estimates are computed from less accurate ROMs obtained by relaxed LRCF solutions. III. TEST SYSTEMS The basic system of this paper is shown in Fig. 1 and relates to a section of a transmission and distribution network from practice. Its full RLC data is given in Appendix B, so results may be reproduced. All transmission lines (TL) of this basic system were then represented by cascaded PI-circuits and no frequency-dependent line parameters were modeled, generating eight large descriptor system models. The performances of four MOR methods were assessed for these eight test systems. The DAEs of the eight test systems were built for both SISO and MIMO cases, by varying the number of PI-sections per TL, which are represented by 10 PIs (S10PI, M10PI), 20 PIs (S20PI, M20PI), 40 PIs (S40PI, M40PI), and 80 PIs (S80PI, M80PI). The dimensions of these DAEs and their spectral radii (moduli of their largest eigenvalues) are shown in Table I. The four SISO systems (S10PI, S20PI, S40PI, S80PI) have the electrical current injected into bus #21 as their input, and the mon-

1908

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

TABLE II METHODS’ PERFORMANCES FOR THE SISO TEST SYSTEMS

Fig. 1. One-line diagram of the electrical power system.

TABLE I NAMES AND STATISTICS OF THE EIGHT TEST SYSTEMS

where is the th frequency response value of the th scalar is the equivalent element of the TF of the original model; and are the numbers of scalar TFs and reduced model; samples for each one of these TFs, respectively. itored voltage at the same bus as their output. The four MIMO systems (M10PI, M20PI, M40PI, M80PI) have the electrical currents simultaneously injected into buses #18, #21, and #23 as their inputs, and the voltages at these same three buses as their outputs.1 The interest here is to determine the ROMs for these test systems described in Table I, which should adequately reproduce, for all practical purposes, their input-output behavior. The complete eigensolutions (eigenvalues and associated right/left eigenvectors) for the test systems can in principle be obtained by either the QR or QZ routines. However, due to the complexity of these methods, where is the number of states, they are not applicable to large-scale systems with many thousands of states. Such large sparse systems can, however, be dealt with by the MOR methods in this paper, that exploit the structure and sparsity of the system matrices. All results of this paper were obtained with Matlab, its toolboxes, and M-code routines. The MOR algorithms’ frequency response, time domain, and pole spectra results are compared with those obtained for the original system. The results of these test systems are described and compared in the next two sections. The adopted definition for model error was the RMS value of the error difference between the frequency response values of the original system and that of the ROM [11]: 1The datafiles for these eight descriptor systems are available for direct download at: http://sites.google.com/site/rommes/software.

IV. RESULTS FOR THE SISO TEST SYSTEMS A. Vector Fitting Method This method is widely used in practice for producing linear ROMs that match the frequency response plots of the original systems [11]. The quality pole-residue results obtained for the four SISO test systems described in Table I, plus its robust convergence fully confirmed the reputation of the method. The results obtained for the S80PI system are described next. The simulations with the other three SISO test systems required similar procedures and the choice of system-specific initial parameters: pole estimates, number of iterations, and the dimension of the reduced models. As already mentioned, a system having many resonant peaks in its frequency response requires a large number of samples. Obviously, when increasing the “resolution” of the frequency response plots, for more accurate VF results, there is an associated rise in CPU time. The CPU time, ROM dimensions, ROM errors, and spectral radii are the four performance indices used in Table II for comparing the various methods. The numbers of DAE equations for the S10PI, S20PI, S40PI, and S80PI systems are given within brackets in this table. All systems have, by construction, the same number of algebraic equations: 119, the remainder being state variables (still containing some redundant states). The QR method is not included in Table II because it is not applicable to large-scale systems in general.

FREITAS et al.: REDUCED-ORDER TRANSFER MATRICES

For the S80PI system, the selected frequency response . The initial pole range was from 10 Hz to estimates were obtained by distributing the imaginary part of the pole estimates logarithmically spaced on the frequency range [11]. In our application, it is not appropriate to distribute these points linearly over the frequency range since the resonant peaks are spread over several decades. The scalar TF was represented by a sequence of 800, 1400, 1800, and 1800 discrete values of complex frequency, for the test systems: S10PI, S20PI, S40PI, and S80PI, respectively. All SISO ROMs were assumed to be strictly proper system models, among other constraints. As there are no deriva[terms of this nature would exist, if matrix tive terms in in (5) was nonzero], it was assumed that such term also did not exist in the reduced model. The VF simulations were carried out using the version 3 of the m-files available at http://www.energy.sintef.no/Produkt/VECTFIT/index.asp. All pole estimates were constrained to remain on the left half of the complex plane for the SISO systems of this section as well as the MIMO systems of the following section. The convergence of the VF method, based on a specified RMS error tolerance, for the four SISO test systems of this paper usually occurred between 7 and 17 iterations. The convergence tolto a RMS error depending on erance varied from a the test system. In all computations, the frequency range of interest was 10 Hz to 40 kHz. The descriptor formulation for the S80PI system has 4182 equation (4063) states, 119 algebraic variables; the descriptor matrix is highly sparse, with only 10 261 nonzero elements. The reduced model, computed by VF, the S80PI VF ROM, has 394 states, as shown in Table II. Accurate ROMs were obtained for the systems S10PI, S20PI, and S40PI (of reduced order 152, 206, and 297, respectively), as also shown in Table II. B. LRCF Method The LRCF-ADI method requires the use of ADI parameters [7], [24]. The computation and choice of these parameters are quite challenging for they impact the LRCF convergence as well as the ROM quality [7], [8]. The LRCF computations were carried out utilizing multiple ADI parameters, all of which were chosen to be real numbers, whose values were equal to the frequencies of some of the dominant resonant peaks, in rad/s. These peaks can be readily determined from the visual inspection of the frequency response plots of the RLC electric power networks. The number of adopted LRCF iterations for the systems S10PI, S20PI, S40PI, and S80PI were 300, 300, 450, and 600, respectively. The order of the reduced model is a function of the smallest Hankel singular value to be retained. Considering that the LRCF is an approximation of the TBR method, it is fair to expect its results to be slightly worse than those of the latter [23]. This expectation was confirmed in most cases, since the obtained reduced models usually had spurious unstable poles. These spurious poles always have TF residues that are of insignificant magnitude, which makes the task of discarding them an easy one (one should note that theory points the TBR method [1] cannot be directly applied to unstable systems).

1909

Fig. 2. Relevant spectra of S10PI and ROMs obtained by the LRCF and Hybrid , with most poles methods. The actual S10PI spectrum goes up to . showing real parts roughly equal to

The LRCF method cannot be applied to systems that have poles which are unstable or very close to zero. However, by , as described in [23], this using the spectral shift problem is circumvented for the test systems. The chosen ADI parameters for the test systems S10PI, , , S20PI, S40PI, and S80PI are the elements in vectors , and , respectively, which are all shown in the following:

where and . In summary, the ROMs were obtained by truncating the balanced realization (computed using the low rank factors) so that the Hankel singular values with magnitudes above error tolerance are preserved. The LRCF computations were performed on the index-0 descriptor (7), derived from the index-2 DAE (1) of the test systems S10PI, S20PI, S40PI, and S80PI, which have 682, 1182, 2182, and 4182 variables (states and algebraic), respectively. The number of algebraic variables (119) and redundant states (35) are the same for all systems. The index-0 descriptor matrices (cf. (7)) are obtained through the elimination of the redundant states and the algebraic variables of the S10PI, S20PI, S40PI, and S80PI systems, and have 528, 1028, 2028, and 4028 states, respectively. The partial spectra of S10PI and of its 158-state ROM obtained by LRCF, the S10PI LRCF ROM, are shown in Fig. 2. The poles of this ROM are seen to match many of the system eigenvalues obtained by the QR eigensolution. Fig. 4 shows the magnitude plots of S80PI and ROMs computed by, amongst others, LRCF and VF, and clearly shows their excellent matching, both at low and high frequency ranges. The RMS errors for both LRCF and VF were smaller than 0.001. It (as pointed by the arrow is clearly seen that up to tip in the figure), there are some minor resonant peaks.

1910

Fig. 3. Relevant spectra of S10PI and ROMs obtained by the SADPA and VF , with most poles methods. The actual S10PI spectrum goes up to . showing real parts roughly equal to

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

Fig. 5. Time simulation plots for the S80PI system and its 394-state ROM obtained by the Hybrid method.

The partial spectra of the S10PI and that of its 152-states ROM obtained by SADPA are shown in Fig. 3. The poles computed by SADPA perfectly match 152 system eigenvalues obtained by the QR eigensolution of the S10PI system. D. Hybrid Method (LRCF + VF)

Fig. 4. TF magnitude versus frequency of self-impedance of bus #21 for the S80PI system and ROMs, obtained by LRCF, SADPA, VF, and Hybrid methods.

C. SADPA Method A similar set of ROMs was obtained with SADPA [9], considering the descriptor system in (1). This method requires that several parameters that control the iterative process be specified, but most of these parameters remained fixed for all test systems of this paper. The minimum and maximum dimensions of the search spaces were 5 and 10, respectively, for all test systems. The number of specified poles (maximum restarts) were 80 (50), 103 (100), 150 (100), and 200 (120) for the S10PI, S20PI, S40PI, and S80PI test systems, respectively. The option “turbo deflation” of the SADPA method [25] was used and allowed a significant computational speed-up. The modal dominance index was “LM” (largest residue magnitude), since practically all dominant poles are near the imaginary axis, yielding a response with numerous resonant peaks. The initial , in all runs. estimate was The magnitude plot of the ROM produced by SADPA, the SADPA ROM, was similar to those obtained by the other two methods; see Fig. 4. The CPU time required by SADPA was however larger than that of the LRCF method.

This section presents the results obtained with the Hybrid method, which combines the fast pole estimation of the LRCF with the quality fitting of the VF, yielding reduced models that preserve the passivity and symmetry of the original system. An additional motivation for the Hybrid method is the undesirable sensitivity of the VF to poor initial pole estimates and to the discretization level of the frequency response. The LRCF initializing step for the Hybrid method used 200, 206, 297, and 394 iterations, to produce the ROMs for the SISO systems S10PI, S20PI, S40PI, and S80PI, respectively. The number of iterations used in all cases was about 30% smaller than the number required when the LRCF method is utilized alone. Two iterations of the VF method were applied next, improving on the poles and residues, so as to produce a good matching between the responses of the original system and the computed ROM. Figs. 2 and 3 show the partial spectra for the full S10PI system and its 152-state ROMs obtained by Hybrid and VF methods, respectively. The performance of the MOR methods were also assessed by the quality of the transient time responses from their ROMs. For the SISO case, a step disturbance of electrical current was injected into bus #21 and the voltage at this bus was monitored. was used in all simulations An integration time step of of Section IV. The plot in Fig. 5 shows the good matching between the time simulation results computed for the full system and for its 394th-order Hybrid ROM. Similar time simulation results (not shown due to space limitations) were obtained for the ROMs produced by the other three methods studied in this paper. E. Discussion on SISO System Results Table II shows the CPU times for the methods being compared, to obtain the ROMs for the SISO systems. The order of

FREITAS et al.: REDUCED-ORDER TRANSFER MATRICES

1911

these ROMs and their RMS errors are also given in this table. The LRCF method presented the lowest CPU time, regardless the test system analyzed. For the test system S10PI the variation in CPU time among the four methods was not significant. The Hybrid method had the second best performance, being faster than SADPA and VF for the three largest test systems, as can be verified in Table II. It should be noted that despite the high number of resonant peaks of the test systems, the ROMs do not contain much of their higher frequency modes, as shown in the pole spectra plots of Figs. 2 and 3. This can also be seen by comparing the spectral radii of the ROMs obtained by the various methods (in Table II) with those of the test systems (in Table I). V. RESULTS FOR THE MIMO TEST SYSTEMS

Fig. 6. Sigma plots for the M80PI system and ROMs, obtained by LRCF, SAMDP, VF, and the Hybrid methods.

The four test systems considered in this section are M10PI, M20PI, M40PI, and M80PI. The input (output) variables for these systems are the injected currents (voltages) at buses #18, #21, and #23 (cf. Fig. 1). Again, the ROMs were produced using the four methods: LRCF, VF, Hybrid, and SAMDP (the MIMO version of SADPA). The model errors were measured over the entire frequency range of interest, and simultaneously considered all the element TFs located in the upper triangular part of the 3 3 . transfer matrix A. MIMO Results Obtained by LRCF The ADI parameters and eigenvalue shift used for systems M10PI, M20PI, M40PI, and M80PI were identical to those used for the LRCF SISO systems. Systems M10PI, M20PI, and M40PI, required 300 iterations, while 400 iterations were needed for M80PI. At every LRCF iteration, a new set of vectors is added to the low-rank Choleski factor of the controllability matrix [23]. The number of vectors in this set is equal to the number of system inputs. Similarly, a set of vectors is also added to the low-rank observability matrix. A large number of inputs and outputs of a MIMO system, therefore, implies a high CPU time as well as large memory requirements for all MOR methods. The quality of the MIMO ROMs can be assessed from their Sigma plots, a control systems denomination for the plots of both the maximum and minimum singular values (SV) of their transfer matrices as a function of frequency, and , respectively [10]. Fig. 6 presents the Sigma plots for the M80PI system and its ROMs obtained by all methods. Fig. 7 shows the absolute errors between the for the full test system and its ROMs. Time domain simulations were computed for the M80PI system and its LRCF ROM. These simulations involved simultaneously injecting into buses #18, #21, and #23 the current signals:

Fig. 7. Absolute deviations between of system M80PI and ROMs, obtained by LRCF, SAMDP, VF, and Hybrid methods.

where

is the step function, , , and . Values , , and are time delays associated with these input signals. The integration time step used in the . numerical simulations was The monitored variables were the voltages at the three buses #18, #21, and #23, but only the voltage at bus #18 is shown. Fig. 8 compares the time simulation plots for the full M80PI system and its LRCF ROM, indicating perfect matching for practical purposes. B. MIMO Results Obtained With SAMDP MOR computations similar to those described in Section V-A were carried out with the SAMDP method. The dimensions of the search spaces were kept between 5 and 10 vectors. The number of specified poles (maximum restarts) were 100 (50), 110 (100), 150 (100), and 230 (130) for the M10PI, M20PI, M40PI, and M80PI, respectively. The “turbo deflation” option and the modal dominance index “LM” were adopted. The initial for all test systems. pole estimate was The computed SAMDP ROMs were of equivalent quality to those ROMs obtained with LRCF for the MIMO test systems. The difference, once again, was in the computational speeds

1912

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

frequency range of interest (see Figs. 6 and 7 and Table III for VF). D. MIMO Results Obtained With the Hybrid Method

Fig. 8. Voltage at bus #18 computed for the full M80PI system and its 457-state ROM obtained by the LRCF method.

Two hundred iterations of the LRCF method were used to generate the initial pole estimates for the VF step of the Hybrid method. Once again, the Hybrid ROMs were produced using only two iterations of the VF method, and the quality of the results compared in a similar manner. The first step of the Hybrid method (200 LRCF iterations) produced a 457-state ROM for the M80PI system, in 39.6 s of CPU time. The two VF iterations carried out next, to improve the pole estimates computed by LRCF, required another 43.3 s of CPU. Therefore, the total time needed to obtain the Hybrid ROM for the M80PI system was equal to 82.9 s (see Figs. 6 and 7 and Table III for Hybrid). E. Discussion on MIMO System Results

TABLE III METHODS’ PERFORMANCES FOR THE MIMO TEST SYSTEMS

Table III presents a summary on the performances of the various methods. Considering the two smallest test systems (M10PI and M20PI), the SAMDP method showed the best performance (CPU time) closely followed by the LRCF and Hybrid methods. On the other hand, for the two largest test systems (M40PI and M80PI), the LRCF method showed the best performance closely followed by the Hybrid method. It is important to point out that for the largest system (M80PI), the LRCF and Hybrid methods showed practically the same performance. Another important fact is that the VF method presented the worst performance regardless of the test system analyzed. The results in Table III show that all ROMs had RMS errors for the frequency range of interest. Once again, below the authors attempted to define a fair basis of comparison for the performances of the four methods, in generating MIMO ROMs of good quality. Different stopping criteria could be used, but the authors believe that any other consistent set of criteria would yield approximately the same ranking of the methods. VI. CONCLUSIONS

of the MOR methods rather than in the quality of the obtained ROMs (see Figs. 6 and 7 and Table III for SAMDP). C. MIMO Results Obtained With Vector Fitting The frequency response discretization adopted for the SISO systems was used again for the MIMO systems. Six scalar TFs and , which correwere generated spond to the elements in the upper triangular part of the transfer , which is symmetric. These six scalar TFs are then matrix imposed to share the same poles, which is a problem equivalent to that of a transfer matrix having one input and six different outputs:

Once the six scalar VF ROMs are identified, in accordance with the previously described constraints, the Sigma plots for the transfer matrices of the obtained MIMO VF ROMs are compared against the original test system models, over the whole

Four MOR methods were compared for the computation of reduced-order transfer function models of descriptor systems having many resonant peaks, to be used for studying high-frequency phenomena in power systems. Four SISO and four MIMO test systems were used in the performance assessment of the MOR methods. The largest system had 4028 effective states, 35 redundant states, and 119 algebraic equations. All four methods yielded quality results, but, in general, the LRCF required less memory and CPU time. The VF was expensive for MIMO systems, but the combination (LRCF + VF), here proposed as the Hybrid method, proved practically as efficient as LRCF and has the advantage that passivity, an important system property, is preserved in the reduced model. The choice of the ADI parameter set for the LRCF method, as the values of the most dominant resonant frequencies, proved simple and effective and is believed to be another contribution of this paper. The spectral radii of the obtained ROMs were seen to be, in most cases, around 60 times smaller than those of the test systems, which is computationally advantageous to both frequency and time domain simulations.

FREITAS et al.: REDUCED-ORDER TRANSFER MATRICES

1913

capacitor. The output variables are the voltages at the two and . It follows that nodes

Fig. 9. Two-node RLC circuit with a redundant state and a current source as input.

In future MOR work, the authors will focus on the preservation of system passivity in the ROMs, for example [26]. The adoption of a quadratic polynomial model and preservation of passivity [25] will be investigated, so as to ensure that the ROMs will retain passivity as well as the RLC structure of the original test system. Approximating the frequency dependence of transmission line parameters by adding RL networks into the original linear descriptor system models will rise their dimension by, at least, an order of magnitude. The new descriptor models could then easily reach dimensions of about 100 000 but the MOR methods of this paper would be adequate to deal with such problems mainly if equipped with robust and reliable stopping criteria and adaptive error bounds. With such high orders, these multi-port network equivalents should be automatically computed by powerful MOR algorithms since their computation would be very laborious if done interactively by the engineer. APPENDIX A REPRESENTATION OF A SINGLE RLC CIRCUIT The RLC circuit in Fig. 9 is used to illustrate the descriptor representation of the system as described in Section II-A. The circuit is represented by six equations instead of the eight equations in [5]. Using the notation of Section II-A, the following matrices are obtained by applying Kirchhoff’s laws to the circuit:

where , and

, . The input variable is is the voltage across the

Since and are zero, the term is zero in (5). Simiis zero, is null. larly, as This circuit has a redundant state since an inductor is connected in series with another inductor. Hence, there are effectively two states, one due to the capacitor and the other due to the two inductors. This physical arrangement, with a redundant state, is correctly represented by the descriptor system model, and . The rank of the latter since matrix is equal to the number of redundant states. The elimination of the algebraic variables and redundant states produces the system representation

(10) (11) where . Numerical values for the elements of the circuit shown in , , , and Fig. 9 are: . When the QZ method is applied to the descriptor system (10), the pair of finite eigenvalues are obtained. The same result was presented in [5]. Moreover, if the index-2 representation with three states and three algebraic variables is used, besides the same finite eigenvalues, four infinite eigenvalues are obtained. APPENDIX B POWER SYSTEM DATA This appendix describes the basic data for the power system in Fig. 1. Table IV shows the nominal voltage , in kV, for each bus. Also, load and shunt equipment data are presented. Each load is modeled by a resistive branch connected to an inductor or a capacitor. For the inductor, the connection is in series and for the capacitor, the connection is in parallel. Table V shows the transmission line data for a single PI section, including the sending-end and receiving-end buses. The parameters of each TL are the resistance, , the inductance, , and the capacitance, . The parameters and form the series branch of a PI section circuit. Half of the capacitance is placed at the sending and receiving-ends of the PI-section. Table VI shows the buses and data related to existing transformers in the power system. The parameters are the resistance, , the reactance, , and the transformer tap. The circuit eleare in percent of an impedance base which is ments and calculated considering a 100-MVA power base and the nominal

1914

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

TABLE IV LOAD AND SHUNT EQUIPMENT DATA

TABLE V TRANSMISSION LINE DATA

TABLE VI TRANSFORMER DATA

voltage of a bus. The frequency of the system is 60 Hz. The data in Tables IV and VI allow to compute the resistance in and the inductance in H for the transformer equivalent circuit. APPENDIX C REDUNDANT STATES AND ALGEBRAIC VARIABLES A procedure for symbolic math elimination of algebraic variables and redundant state variables of an index-2 DAE system [22] is described in this appendix. Assume that is the operator of the Laplace transform and that the initial conditions of the system are zero. In the Laplace domain, system (1) can be written as (12) From the structure of matrix and (2), one observes algebraic equations in (1). If (see, e.g., [23]), no redundancy exists and the algebraic equations can be algebraic variables. In this paper, however, used to eliminate

. Hence, the difference gives the number of redundant states. Equation (12) including the output equation and definitions from (2) and (3) can be written as the augmented system (13) where

.

FREITAS et al.: REDUCED-ORDER TRANSFER MATRICES

1915

Since , (13) has algebraic variables and redundant states, both of which could be eliminated symboliis singular, a Gaussian elimination process can cally. As be carried out through Kron’s reduction technique [27]. In this process, the operator is handled as a symbolic variable. At the end of the Gaussian elimination process, it is used to transform the model of the resultant system to the time domain. must be submitted to Assume that a matrix . The new a Kron reduction by using as pivot the element of a matrix are elements . Since the algebraic equations will be eliminated, each pivot and . Let the pivot should be chosen from the matrices be chosen from a nonzero row of this latter matrix. The reduction process must be repeated until the elements in the modified become all equal zero. This transforms the set of matrix (13) into (14) ; and , where , , , , comes from the initial Kron reduction. Note that as the pivot elements were chosen from matrix , there was no change in . , as exWithout loss of generality, we assumed that plained in (6). Therefore, the equations and variables related to are eliminated next. This process the system leads to the elimination of redundant states. Choosing the pivots results in from (15) The redundant state equations are then eliminated. The pivots are chosen from and the variables eliminated, yielding (16) which in the time domain, assuming that given by

, is

(17) ACKNOWLEDGMENT The authors would like to thank the Brazilian Electrical Energy Research Center—CEPEL—for permission to use its power system software to generate the DAE models. REFERENCES [1] A. C. Antoulas, Approximation of Large-Scale Systems. Philadelphia, PA: SIAM, 2005.

[2] W. H. A. Schilders, H. A. van der Vorst, and J. Rommes, Model Order Reduction—Theory, Research Aspects and Applications. New York: Springer, 2008. [3] A. S. Morched, J. H. Ottevangers, and L. Marti, “Multi-port frequency dependent network equivalents for the EMTP,” IEEE Trans. Power Del., vol. 8, no. 3, pp. 1402–1412, Jul. 1993. [4] G. Jackson, U. D. Annakkage, A. M. Gole, D. Lowe, and M. P. McShane, “A real-time platform for teaching power system control design,” in Proc. Int. Conf. Power System Transients, Montreal, QC, Canada, Jun. 2005. [5] S. L. Varricchio, N. Martins, and L. T. G. Lima, “A Newton-Raphson method based on eigenvalue sensitivities to improve harmonic voltage performance,” IEEE Trans. Power Del., vol. 18, no. 1, pp. 334–342, Jan. 2003. [6] M. A. Rahman, A. Semlyen, and M. R. Iravani, “Two-layer network equivalent for electromagnetic transients,” IEEE Trans. Power Del., vol. 18, no. 4, pp. 1328–1335, Oct. 2003. [7] T. Penzl, “A cyclic low-rank Smith method for large-scale Lyapunov equations,” SIAM J. Sci. Comput., vol. 2, no. 4, pp. 1401–1418, 2000. [8] T. Penzl, “Algorithms for model reduction of large dynamical systems,” Lin. Alg. Appl., vol. 415, no. 2–3, pp. 322–3431, Jun. 2006. [9] 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. [10] J. Rommes and N. Martins, “Efficient computation of multivariable transfer function dominant poles using subspace acceleration,” IEEE Trans. Power Syst., vol. 21, no. 4, pp. 1471–1483, Nov. 2006. [11] B. Gustavsen and A. Semlyen, “Rational approximation of frequency domain responses by vector fitting,” IEEE Trans. Power Del., vol. 14, no. 3, pp. 1052–1061, Jul. 1999. [12] B. Gustavsen, “Improving the pole relocating properties of vector fitting,” IEEE Trans. Power Del., vol. 21, no. 3, pp. 1587–1592, Jul. 2006. [13] D. Deschrijver, M. Mrozowski, T. Dhaene, and D. D. Zutter, “Macromodeling of multiport systems using a fast implementation of the vector fitting method,” IEEE Microw. Wireless Compon. Lett., vol. 18, no. 6, pp. 383–385, Jun. 2008. [14] W. S. Meyer and H. W. Dommel, “Numerical modelling of frequency-dependent transmission line parameters in an electromagnetic transients program,” in Proc. IEEE PES Winter Meeting, New York, Feb. 1, 1974, pp. 1401–1409. [15] A. Semlyen and A. Dabuleanu, “Fast and accurate switching transient calculation on transmission lines with ground return using recursive convolutions,” IEEE Trans. Power App. Syst., vol. PAS-94, no. 2, pp. 361–371, Mar./Apr. 1975. [16] J. R. Marti, “Accurate modelling of frequency-dependent transmission lines in electromagnetic transient simulations,” IEEE Trans. Power App. Syst., vol. PAS-101, no. 1, pp. 147–157, Jan. 1982. [17] J. A. Martinez-Velasco, Computer Analysis of Electrical Power System Transients. New York: IEEE Press, 1997. [18] S. Gomes, Jr, N. Martins, and C. Portela, “Sequential computation of dominant poles of s-domain system models,” IEEE Trans. Power Syst., vol. 24, no. 2, pp. 776–784, May 2009. [19] M. C. Tavares, J. Pissolato, and C. M. Portela, “Mode domain multiphase transmission line model—use in transient studies,” IEEE Trans. Power Del., vol. 14, no. 4, pp. 1533–1544, Oct. 1999. [20] B. Gustavsen, “Computer code for rational approximation of frequency dependent admittance matrices,” IEEE Trans. Power Del., vol. 17, no. 4, pp. 1093–2002, Oct. 2002. [21] L. T. G. Lima, N. Martins, and S. Carneiro, Jr, “Augmented state-space formulation for the study of electrical networks including distributed transmission line models,” in Proc. Int. Conf. Power Systems Transients (IPST’99), Budapest, Hungary, Jun. 24, 1999, pp. 87–92. [22] Y. Ren, “The approximate regularization of index-2 differential-algebraic problems,” Elsevier, J. Comput. Appl. Math., vol. 78, pp. 301–310, 1997. [23] F. D. Freitas, J. Rommes, and N. Martins, “Gramian-based reduction method applied to large sparse power system descriptor models,” IEEE Trans. Power Syst., vol. 23, no. 3, pp. 1258–1270, Aug. 2008. [24] J. R. Li and J. White, “Low-rank solution of Lyapunov equations,” SIAM J. Matrix Anal. Appl., vol. 24, pp. 260–280, 2002. [25] J. Rommes and N. Martins, “Computing transfer function dominant poles of large-scale second order dynamical systems,” SIAM J. Sci. Comput., vol. 30, no. 4, pp. 2137–2157, 2008. [26] B. Porkar, M. Vakilian, R. Iravani, and S. M. Shahrtash, “Passivity enforcement using an unfeasible-interior-point primal-dual method,” IEEE Trans. Power Syst., vol. 23, no. 3, pp. 966–974, Aug. 2008. [27] J. Arrillaga, B. C. Smith, N. R. Watson, and A. R. Wood, Power System Harmonic Analysis. New York: Wiley, 1997.

1916

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

Francisco Damasceno Freitas (SM’08) received the B.Sc. and M.Sc. degrees from the University of Brasilia, Brasilia, DF, Brazil, in 1985 and 1987, respectively, and the Ph.D. degree from the Federal University of Santa Catarina, Florianopolis, SC, Brazil, in 1995, all in electrical engineering. Since 1986, he has worked for the University of Brasilia, where he is an Assistant Professor. His area of interest is power system dynamics.

Joost Rommes received the M.Sc. degree in computational science, the M.Sc. degree in computer science, and the Ph.D. degree in mathematics from Utrecht University, Utrecht, The Netherlands, 2002, 2003, and 2007, respectively. He is currently a researcher at NXP Semiconductors, Eindhoven, The Netherlands, and works on model reduction, specialized eigenvalue methods, and algorithms for problems related to circuit design and simulation.

Nelson Martins (F’98) received the B.Sc. degree from the University of Brasilia, Brasilia, DF, Brazil, in 1972 and the M.Sc. and Ph.D. degrees from the University of Manchester, Manchester, U. K., in 1974 and 1978, respectively, all in electrical engineering. He has worked for CEPEL, the Brazilian Electrical Energy Research Center, Rio de Janeiro, Brazil, since 1978, developing methods and tools for power system dynamics and control.

Franklin C. Véliz received the B.Sc. and M.Sc. degrees in electrical engineering from the Federal University of Rio de Janeiro (UFRJ), Rio de Janeiro, Brazil, in 2001 and 2005, respectively. He works for Pontific Catholic University (PUC), Rio de Janeiro, developing methods and tools for power system analysis.

Sergio Luis Varricchio (SM’06) received the B.Sc. degree from the Catholic University of Petrópolis (UCP), Petrópolis, Brazil, in 1987 and the M.Sc. degree from the Federal University of Rio de Janeiro (UFRJ), Rio de Janeiro, Brazil, in 1994, both in electrical engineering. Since 1989, he has worked for CEPEL, Rio de Janeiro, in modal analysis, power quality, and electromagnetic transients.