IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 25, NO. 2, MAY 2010

929

Computing Rightmost Eigenvalues for Small-Signal Stability Assessment of Large-Scale Power Systems Joost Rommes, Nelson Martins, Fellow, IEEE, and Francisco Damasceno Freitas, Senior Member, IEEE

Abstract—Knowledge of the rightmost eigenvalues of system matrices is essential in power system small-signal stability analysis. Accurate and efficient computation of the rightmost eigenvalues, however, is a challenge, especially for large-scale descriptor systems. In this paper we present an algorithm, based on subspace accelerated Rayleigh quotient iteration (SARQI), for the automatic computation of the rightmost eigenvalues of large-scale (descriptor) system matrices. The effectiveness and robustness of the algorithm is illustrated by numerical experiments with realistic power system models, and we also show how SARQI can be used to compute eigenvalues closest to any damping ratio and repeated eigenvalues. The algorithm can be used for stability analysis in any other field of engineering. Index Terms—Eigenvalues, eigenvectors, large-scale eigenvalue problems, poles, poorly-damped oscillations, power system stability, small-signal stability analysis, sparse systems, specialized eigensolvers, system oscillations, transfer functions.

I. INTRODUCTION

T

HE assessment of the small-signal stability of a mid-sized power system, at various operating points, is obtained by sequential computation of eigenvalue decompositions of the state matrices by the QR method [1]. However, the use of the powerful QR method, and more recently the BR method [2], is not feasible for large-scale power system studies due to their excessively large CPU time and memory requirements. Usually it is not necessary to determine all poles of a power system model, but only those few poles (mostly associated with electromechanical modes) that lie near the imaginary axis and relate to poorly-damped or even negatively damped (unstable) oscillatory modes. Except for very rare instances, unstable real poles are due to very badly adjusted controllers that may be easily found by gross-error data detection routines at a data validation stage [3], [4]. The S-matrix or Cayley transform [5], applied to the state matrix , maps the stable eigenvalues into the unit circle and has been used with moderate success in sparse power system descriptor formulations, mostly in association with Lanczos and Arnoldi methods, for the iterative solution of the rightmost Manuscript received August 05, 2009. First published December 31, 2009; current version published April 21, 2010. This work was supported by EU Marie-Curie project O-MOORE-NICE! Paper no. TPWRS-00419-2009. J. Rommes is with NXP Semiconductors, 5656 AE Eindhoven, The Netherlands (e-mail: [email protected]). N. Martins is with CEPEL, Rio de Janeiro CEP-21944-970, Brazil (e-mail: [email protected]). F. D. Freitas is with the Department of Electrical Engineering, University of Brasilia, Brasilia, DF, CEP:70910-900, 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.2009.2036822

eigenvalues [6]–[8]. A more general variant of the Cayley transform, the Möbius transform [7]–[10], maps the eigenvalues of that have damping ratios above some specified value into the unit circle and has also been moderately used. Despite several advances in both theory and practice [6], [8], [10]–[13], the whole family of Möbius transforms continues to suffer from its fundamental weakness: those poorly-damped eigenvalues of , once spectrally transformed, become quite close to each other in moduli. This fact hampers convergence to the rightmost eigenvalues of the algorithms that use subspace iteration or Arnoldi methods [14] on the Cayley/Möbius transforms of the sparse, descriptor formulation of the power system small-signal stability problem. The authors’ motivation for developing an algorithm for computing rightmost eigenvalues, based on subspace accelerated Rayleigh quotient iteration (SARQI), is the standing difficulty in finding the complete set of unstable poles of the state matrix (located at the right half-plane of the complex plane) for very-large power system models. This difficulty is magnified by the fact that some of these unstable poles may have very weak observability in major variables or nodes of the interconnected system. The proposed technique is also especially interesting in situations where these unstable poles have very small residues (weak controllability and observability) in those mostly evident transfer functions—either single-input single-output (SISO) or multi-input multi-output (MIMO)—and so may not be easily detected by other methods that take into account modal dominance. Though not considered as the ultimate rightmost eigenvalues algorithm, the accumulated experience on SARQI points it as a major improvement to the existing technology and capable to face the challenges of online power system small-signal stability assessment [9]. The proposed technique basically is a subspace accelerated Rayleigh quotient iteration method (SARQI) [15] with a damping ratio selection criterion for the eigenvalue estimates, whose mathematical justification and large-scale system results are all given in this paper. The large-scale matrix data used to generate some of the results are also made available online [16] so as to allow readers to verify the SARQI results reported in this paper and compare them against those produced by other methods. It is also shown that the developed subspace method SARQI is very robust and performs well in the presence of repeated eigenvalues that are unstable, as long as the system analyzed is not defective (i.e., the eigenvectors of the system must form a complete basis). The outline of the paper is as follows. In Section II dynamical systems and the small-signal stability analysis eigenvalue

0885-8950/$26.00 © 2009 IEEE Authorized licensed use limited to: IEEE Xplore. Downloaded on April 21,2010 at 04:00:39 UTC from IEEE Xplore. Restrictions apply.

930

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 25, NO. 2, MAY 2010

problem are described. The intrinsic (numerical) difficulties of computing the rightmost eigenvalues are explained in Section III. The SARQI for computing rightmost eigenvalues is described in Section IV. Applications of SARQI to large-scale practical power system models, showing the effectiveness and robustness of SARQI, are presented in Section V. Section VI concludes. II. DYNAMICAL SYSTEMS AND RELATED EIGENVALUE PROBLEMS The motivation for this paper comes from dynamical systems of the form (1) , singular, is the state vector, is the output vector, and sponding transfer function (matrix) fined as where

, , is the input vector, . The correis de(2)

. If , the system (1) is called a SISO with , system (1) is called a MIMO system. system. If denotes the complex-conjugate Throughout this paper, denotes the inverse of the transpose of a matrix , and complex-conjugate transpose of . of the matrix pencil are The eigenvalues the poles of transfer function (2). The problem of computing and corresponding right and left eigenvecthe eigenvalues tors is also known as the generalized eigenvalue problem [1]:

If , the problem is called an ordinary or standard eigenvalue problem. One can transform the generalized eigenvalue problem, with singular, into an ordinary eigenvalue problem by applying shift-and-invert [1], [6], [13]:

where for . Another way is to apply a (generalized) Cayley transform [5]:

where . Assuming that the pencil is nondefective, the right and left eigenvectors and corresponding to finite eigenvalues can be scaled so that . Furthermore, right and left eigenvectors corresponding to disfor . The tinct eigenvalues are -orthogonal: can be expressed as a sum of residue matransfer function over finite first-order poles [17]: trices

where the residues

are

and is the number of finite first-order poles, and the contribution of poles at infinity is assumed to be zero (and hence poles at infinity are of no physical importance). In small-signal stability analysis one is interested in assessing whether the system, linearized around an operating point, is dynamically stable or unstable under small disturbances. Since the states of the system in time domain can be written , it is easy to as a linear combination of the terms see that the presence of eigenvalues with will result in states that grow exponentially with time [17], [18]. In this case the operating point, or more generally the system, is called unstable. Hence, the small-signal stability problem boils down to determining whether there are any finite eigenvalues with positive real part. The main goal of this paper is to present a robust algorithm for the small-signal stability problem. 1) Problem 2.1: Given (descriptor) system matrices and , compute the finite eigenvalues (and corresponding eigenvectors) with largest real part (i.e., the rightmost eigenvalues). In general, one is not only interested in the eigenvalues with positive real part, but also in the eigenvalues closest to the imaginary axis in the left half-plane. In small-signal stability analysis, these eigenvalues are characterized as having a small damping is ratio, where the damping ratio of an eigenvalue defined as [19] (3) Roughly speaking, eigenvalues with small damping ratio correspond to frequency modes that are poorly damped, and hence these eigenvalues are of interest for small-signal stability as well [9], [20], [21]. Therefore, we also address this problem. 2) Problem 2.2: Given (descriptor) system matrices and , compute the finite eigenvalues (and corresponding eigenvectors) with smallest damping ratio. Note that Problem 2.1 is in fact a special case of Problem 2.2, since eigenvalues with positive real part have negative damping ratio. However, for clarity of presentation we will consider the two problems separately. III. NUMERICAL DIFFICULTIES OF COMPUTING RIGHTMOST EIGENVALUES The problem of computing the rightmost eigenvalues of matrices or matrix pencils arises in many application areas, including analysis of steady state solutions of Navier-Stokes equations [5], [6], [13] and power system and analogue circuit small-signal stability analysis [22], [23]. For matrices of , one typically uses full space dimensions up to methods such as the QR method for the ordinary eigenproblem, and the QZ method for the generalized eigenproblem [1], to compute the full spectrum. However, for large-scale dynamical and for the (descriptor) systems one easily has corresponding eigenproblems the QR and QZ methods are no longer applicable, due to their memory and CPU require. Hence, one uses iterative methods such as the ments of

Authorized licensed use limited to: IEEE Xplore. Downloaded on April 21,2010 at 04:00:39 UTC from IEEE Xplore. Restrictions apply.

ROMMES et al.: COMPUTING RIGHTMOST EIGENVALUES

Lanczos method [24], [25], the Arnoldi method [26], [27], and the Jacobi-Davidson method [28], [29] that exploit the sparsity of the systems to limit memory and CPU requirements. Since these iterative methods are especially designed for computing a few eigenvalues,1 the major challenge is to use these methods in such a way that they compute all the eigenvalues of interest. In our case the eigenvalues of interest are the rightmost eigenvalues and the eigenvalues with smallest damping ratio. It is clear that in any case all eigenvalues with positive real part must be computed when performing small-signal stability analysis: failing to find at least one eigenvalue with positive real part, in the presence of eigenvalues with positive real part, leads to a wrong conclusion on the stability of the system. and , one For systems with nonsingular matrices usually adopts a strategy to select the rightmost eigenvalue approximations at every iteration, in order to steer the iterative method to converge to the rightmost eigenvalues. This approach has shown to be successful in many cases; see, for instance, [11], [13], [27], [30], and [31]. Here we note that there are no waterproof ways to detect whether all the rightmost eigenvalues are found, although quite robust techniques exist [11]. On the other hand, checking whether a found eigenvalue is indeed an eigenvalue of the matrix at hand is straightforward and hence conclusions on the stability of the system can be drawn with great certainty if eigenvalues with positive real part are found. For descriptor systems with singular matrix , the situation is different. Due to the singularity of there are eigenvalues at infinity. Hence, adopting a strategy to select the rightmost eigenvalue approximations without any precautions by no means prevents approximations of eigenvalues at infinity to be selected. In finite arithmetic, however, one never converges , and the corresponding exact eigenvectors are found to neither. In other words, there is no stable way to determine whether an approximation of an eigenvalue at infinity is indeed accurate, as will be shown by a numerical example in Section V-B. One could argue that transforming the generalized to the ordinary eigenproblem eigenproblem , where would solve the issues with eigenvalues at infinity, since these are now transformed to eigenvalues at zero, while eigenvalues with positive real part are transformed to eigenvalues with positive real part. As is discussed in detail in [5], [6], [13], this transformation only cures the symptoms: in finite arithmetic, the iterative methods that are will consider eigenvalue approximations of sufficiently close to zero as being converged, as soon as the corresponding residuals are below the predefined tolerance , where a typical tolerance is ( ). In finite arithmetic one rarely obtains , instead, with , . If the sign of but is positive, this may lead to the wrong conclusion that the system is unstable, since it is hard to determine whether is (and hence possibly an approximation of an eigenvalue (a so-called an eigenvalue with positive real part) or 1In principle, iterative methods can be used to compute the complete spectrum of a matrix (pencil), but the memory and CPU costs are in that case of the same order as those of the QR and QZ methods.

931

spurious eigenvalue, which corresponds to an eigenvalue at infinity). See also Section V-B for a numerical example. With the preceding discussion in mind, it is not surprising that there have been developments in methods that try to prevent spurious approximations of eigenvalues at infinity (or equivalently, at zero) from entering the iterative process. This technique is also known as purification and has been introduced in [5] and further developed and refined in [6] and [13]. The key idea of purification is to keep the search spaces clean of any components in the direction of eigenvectors that correspond to eigenvalues at infinity. Purification has been shown to be robust for several examples [6], [13], but has as drawbacks that knowledge about the structure of the matrices is required, and that the process can still be disturbed by round off errors. The algorithm presented in the next section does not assume any structure of the system matrices and is conceptually simpler than other methods, while still robust and accurate in the computation of the rightmost eigenvalues. Since the eigenvalues at infinity are of no physical importance, in the following we will no longer mention explicitly that we are interested in the finite rightmost eigenvalues. However, we will continue to stress that eigenvalues at infinity may cause numerical difficulties, and that algorithms should take special care of that. Another way to circumvent problems due to eigenvalues at infinity is to work with the state space matrices instead of with the descriptor matrices, either implicitly or explicitly [32]. Recently this technique has been used for solving Lyapunov equations [33], [34] and also in the context of eigenvalue problems [31]. For power systems, and have the structure (4) where plement [1]

is nonsingular. We can use the Schur com-

but in general it is not practical (nor needed) to construct explicitly, since it will be dense. Instead, it is more economical to implement product and solves with implicitly, by using the Schur complement [32]. This approach is efficient if per iteration either only a solve or a product is needed. However, for algorithms such as SADPA [35], one needs both a solve and a product per iteration, which may make this approach more expensive than working with descriptor matrices, as is also shown in Section V. For more details on (implicit) state space computations, see [31]. IV. ALGORITHM FOR RIGHTMOST EIGENVALUES In Section IV-A we give a brief summary of iterative methods with subspace acceleration. In Section IV-B we describe how the subspace accelerated dominant pole algorithm (SADPA) [35], with a special selection strategy, can be used to compute rightmost eigenvalues. In Section IV-C we show that in fact a SARQI [36] with a special selection strategy can be used to compute eigenvalues of any damping ratio, and we explain why this is numerically robust. We revisit computation of rightmost eigenvalues with positive real part, this time using SARQI, in Section IV-D.

Authorized licensed use limited to: IEEE Xplore. Downloaded on April 21,2010 at 04:00:39 UTC from IEEE Xplore. Restrictions apply.

932

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 25, NO. 2, MAY 2010

A. Iterative Methods With Subspace Acceleration Single vector iterations such as the dominant pole algorithm [37] and Rayleigh quotient iteration [36] compute approximate eigenvectors based on the current shift only. The idea behind subspace acceleration is that we take previously computed eigenvector approximations into account as well: in the th and are kept iteration, the approximate eigenvectors and , respectively. in search spaces In the th iteration, this leads to a projected eigenproblem of order . This small problem can be solved exactly using the QZ method [1], leading to approx. Of these triplets, the one that imate eigentriplets of satisfies a specific criterion, e.g., rightmost, best is selected, and the corresponding eigenvalue approximation is used as shift for the next iteration. satisfies the converIf an approximate eigentriplet gence criterion (with ), deflation is used to avoid convergence to this eigentriplet again: the search are kept orthogonal to and , respecspaces and tively, during the computation of other eigenvalues. The general framework for an iterative method with subspace acceleration is shown in Algorithm 1. We will refer to this framework in the following subsections. Algorithm 1: Framework for Iterative Method With Subspace Acceleration INPUT: Unreduced Jacobian system matrices and selection criterion, number of wanted eigenvalues OUTPUT: Eigentriplets of 1) for

,

satisfying selection criterion

until convergence do

2) Compute new search space vectors (using current approximate eigentriplet 3) Expand the search spaces 4) Project:

and

See [15], [35], and [38] for more details, including full code for Matlab, on subspace acceleration in the context of computing dominant poles. B. Computing Rightmost Eigenvalues With SADPA Usual dominance indices for poles are based on the (magnitude of the) residues (sometimes scaled by the inverse of the magnitude of the real part of the corresponding poles); see [35] and [39]. Since the residue of a pole is defined as , the left and right eigenvectors and play an important role. Hence, in SADPA approximate eigenvectors are used to select the most dominant pole approximations; see [35] for more details. However, in general there is no way to recognize rightmost poles by their corresponding eigenvectors. Consequently, the selection criterion in subspace based methods typically depends on the pole approximations only—select the rightmost pole approximation as shift for the next iteration. In the presence of poles at infinity this leads to numerical difficulties, as explained in detail in the previous section. Therefore, we propose to use a combined selection criterion where the approximate pole with largest

is selected. This criterion consists of two components. Firstly, . we define as input and output vectors The reasoning is as follows: rightmost poles can not be recognized by their eigenvectors, so all poles are equally dominant only (in other when looking at the residue words, we do not discriminate with respect to eigenvectors). to make poles at infinity the least Secondly, we divide by dominant:

and

with

) and

and

5) Select approximate eigentriplet from according to selection criterion 6) Check convergence and deflate if converged 7) endfor The nature of the method is determined by Step 2 and 5. In SADPA [35], Step 2 is

Note that since the sign of is taken into account, approximate poles with positive real part are the most dominant in this index. The advantage of this approach is that SADPA can be used directly: only the selection criterion has to be adapted. As is described in Section V, this approach is successful if there are indeed poles in the right half-plane. As soon as all finite poles in the right half-plane have been found, or if there are no finite poles in the right half-plane at all, the method stagnates. Faced with this practical inconvenience, we developed more robust selection criteria, as described in the next subsections. C. Computing Eigenvalues With Smallest Damping Ratio

while in SARQI we have

Note that in both cases , i.e., the approximate eigenvalue that was selected in Step 5 (e.g., the most dominant or rightmost eigenvalue).

It is clear that any selection criterion for computing rightmost poles should consider approximations of poles at infinity as the least dominant. Hence, simply selecting the rightmost pole approximations does not satisfy this requirement. We have also seen that in general the eigenvectors reveal no information on ”how rightmost” poles are. Therefore, one may expect that the selection criterion does not involve the approximate eigenvectors. Indeed, the selection criterion we propose is to select the with approximate pole closest to a certain damping ratio

Authorized licensed use limited to: IEEE Xplore. Downloaded on April 21,2010 at 04:00:39 UTC from IEEE Xplore. Restrictions apply.

ROMMES et al.: COMPUTING RIGHTMOST EIGENVALUES

[c.f (3)]. This means that in iteration we sort the approximate poles in increasing order of distance of their damping ratio

to

In the next iteration of SARQI we take as shift. A typical damping ratio of interest for power system applications is , but this selection strategy works for any damping ratio with . This criterion has been used in combination with other methods (Arnoldi, Jacobi-Davidson) [22], [23] as well, but not in the context of computing rightmost eigenvalues, and moreover, not in the presence of eigenvalues at infinity. Here we explain now why the damping ratio selection strategy can safely be used in the presence of eigenvalues at infinity. Since no active purification is performed to keep the search spaces clean from components in the direction of eigenvectors corresponding to eigenvalues at infinity (as used in [6] and [13]), we can be sure that eventually such components enter the search spaces; with our selection strategy we prevent following these directions. Let a finite approximation of an infinite eigenvalue , with . Note that be given by has not converged to , it will have a real and as long as an imaginary part; since we are computing in finite arithmetic and hence convergence to infinity is not possible, this will be the case throughout the complete process. A reasonable assumption , with a corresponding damping ratio of is that

933

where one typically uses fixed shifts (manually chosen), in SARQI the shift is automatically updated every iteration using the approximations available. Remains the choice of the input and output vectors and : in , as before (note principle one can take that and no longer play a role in the selection criterion). However, since and have no critical influence on the process, nothing prevents us from using SARQI instead of SADPA: in other words, we update the right hand-sides every iteration with the latest eigenvector approximations, cf. Algorithm 1. This increases the speed of convergence in the neighborhood of a pole from quadratic to cubic, and no longer requires a choice for and . D. Computing Rightmost Eigenvalues With SARQI Here we distinguish between complex conjugated pairs of eigenvalues and real eigenvalues. For complex eigenvalues, the damping ratio selection strategy as described in the previous , for instance, we find section suffices: by setting damping ratio line. In printhe eigenvalues closest to the to , which corresponds to ciple, we can decrease targeting positive real poles. This way, we effectively scan the right half-plane on all possible unstable eigenvalues. , we also increase the chance of seOf course, as lecting approximations of eigenvalues at infinity, as described in the previous section. However, their corresponding damping (in finite arithmetic), while ratios will never reach the damping ratios of approximate finite real eigenvalues will , and hence will always be selected quickly converge to during the SARQI process. E. Complete Algorithm The complete algorithm, with SARQI at its core, for computing eigenvalues within a specific range of damping ratios is presented in Algorithm 2. Algorithm 2: Algorithm for Computing Rightmost Eigenvalues, or Eigenvalues Within Damping Ratio Range . INPUT:Unreduced Jacobian system matrices and , , number of wanted damping ratio range eigenvalues per damping ratio OUTPUT:Eigentriplets of

If tends to become (much) smaller than (since the eigenwill become even larger. Hence, values at infinity are real), , we do not exwith small target damping ratios, i.e., pect any influence of approximations of infinite eigenvalues on the convergence to the finite eigenvalues with small damping ratio. This is confirmed by all numerical experiments described in Section V. The advantage of selecting the best (rightmost, smallest damping ratio) eigenvalue approximation as shift for the next iteration is that at every iteration the spectrum is re-transformed in order to enforce convergence to the eigenvalues of interest. Unlike with shift-and-invert, Cayley, and Möbius transforms,

1)

in damping ratio range //matrices for found eigentriplets

2) for all

do

3) //deflate against found eigentriplets 4) 5) 6) 7) end for

Authorized licensed use limited to: IEEE Xplore. Downloaded on April 21,2010 at 04:00:39 UTC from IEEE Xplore. Restrictions apply.

934

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 25, NO. 2, MAY 2010

TABLE I RESULTS FOR SEVERAL LARGE-SCALE POWER SYSTEM MODELS AND METHODS, FOR COMPUTING THE 15 RIGHTMOST EIGENVALUES WITH SMALLEST DAMPING POLES ARE FOUND, THIS IS RATIO. THE NUMBER OF STATES IS DENOTED BY . FOR EACH METHOD THE NUMBER OF FOUND RHPS IS SHOWN. IF ONLY . NOTATION (IN THE COLUMNS #RHPS) MEANS INSTANCES ARE FOUND OF POLE WITH MULTIPLICITY DENOTED WITH

Several remarks are in place for Algorithm 2: • Although presented here for clarity as a sequence of calls to SARQI, in our actual implementation we used a single call to SARQI where we provided the complete range of damping ratios at once, and looped over the ratios inside SARQI. Numerically speaking, there is no difference. • In several numerical experiments we have observed that when providing sufficiently spaced damping ratios (e.g., at least 0.05 difference), there was practically no risk of recomputing already found eigenvalues, even without deflation. Nevertheless, we advise to always include deflation against found eigentriplets (cf. Step 3, where the found eigenvectors in and are used for this purpose) to reduce the risk of recomputing found eigentriplets to zero. See [35] for more details on deflation. eigen• It might be the case that there are fewer than values near a specific damping ratio . In this case, SARQI will stagnate or show slower convergence (i.e., needing more iterations per found eigenvalue) as the remaining eigenvalues have damping ratios further away from . In the left half-plane, where generally most of the eigenvalues are located, stagnation is usually very limited. In the right ) half-plane, however, where typically only a few (say eigenvalues are located, stagnation will start as soon as all right half-plane eigenvalues have been found. Based on numerous experiments with the present method, and on years of experience with specialized iterative eigenvalue methods such as SADPA [35] and Jacobi-Davidson [13], [28], we have observed that as long as there are remaining eigenvalues close to the target damping ratio, convergence typically happens within much fewer than 25 iterations. Hence, we have built in logics to ensure that after more than 25 iterations of no convergence, SARQI automatically switches to the next target damping ratio. Similarly, it may eigenalso be the case that there are more than values near a specific damping ratio. In this case we use the same check on stagnation, but now SARQI continues with searching for eigenvalues until stagnation is detected. Note that these settings can be specified by the user. V. NUMERICAL RESULTS This section describes numerical results related to the smallsignal stability analysis of power systems. First, numerical difficulties due to eigenvalues at infinity are illustrated by results obtained for a 41-state model of a five-machine power system.

Then, the ability of SARQI to find the rightmost finite eigenvalues in the presence of eigenvalues at infinity, is illustrated by results obtained for the same system. Finally, SARQI is used to compute rightmost eigenvalues and eigenvalues with small damping ratio for large-scale practical power system models. It is shown that SARQI is more efficient and more robust compared to existing methods. We also show that the method is able to compute repeated eigenvalues, and provide comparisons between descriptor and state space computations. and 10% The dashed lines in some figures denote the damping ratio borders. Unless stated otherwise, initial vectors , and for SARQI were chosen as was used. For the experiments with iminitial shift plicitly restarted Arnoldi and Cayley shifts, we used the shift selection strategy described in [13]. All experiments were carried out in Matlab 7 on a MacBook Pro (2-GHz Intel Core Duo, 1 GB of RAM). Unless stated otherwise, all computations were done with sparse descriptor realizations. We will refer to SARQI in the following, both for the complete process of computing rightmost/poorly damped eigenvalues, as well as to discuss specific convergence properties. Since the descriptor models for all test systems used in this paper only involve real matrices, their eigenvalue spectra are symmetric w.r.t. the real axis. A. Large-Scale Test Systems The performance of the proposed algorithm was assessed for several large-scale test-systems derived from Brazilian Interconnected Power System (BIPS07) models related to a 2007 heavy load condition. BIPS07 has 173 synchronous machines, 3584 buses, 5056 branches, 3469 nonlinear loads, eight HVDC converters, and eight FACTS devices. The descriptor system model has 2898 states, and 20 383 variables in total. Its operating point is stable. In order to evaluate the performances of the proposed algorithm on systems with more stringent conditions, four other unstable test systems were derived by modifying generator controller parameters and topology of the BIPS07, as described in the following. Some main characteristics and results for the systems can also be found in Table I. 1) The 3012-State System (Xingo3012): This system was generated by modeling the Xingo power station with each of its six identical units connected individually to a single bus. Changes were carried out in the automatic voltage regulator (AVR) gains and power system stabilizers (PSS) of the power stations of Itaipu, Angra I, Tucurui and Jacui. The resulting

Authorized licensed use limited to: IEEE Xplore. Downloaded on April 21,2010 at 04:00:39 UTC from IEEE Xplore. Restrictions apply.

ROMMES et al.: COMPUTING RIGHTMOST EIGENVALUES

system is unstable with the following right half-plane (RHP) , , and 2.94. poles: 2) The 2940-State System (Xingo6u): The 3012-state system has five (stable) repeated poles related to the intra-plant modes of the six-unit Xingo power station. To destabilize the intra-plant poles, the PSSs of the six Xingo generator units were removed and negative damping torques were introduced by adding mechanical damping terms to their rotor-speed deviations, which have equal magnitudes but the “wrong” polarities. The resulting system has the same unstable poles as the Xingo3012 system: , , 2.94, which remained practically unchanged (being unrelated to the Xingo power plant dynamics), plus the five new repeated (unstable) Xingo intra-plant poles at . 3) The 5727-State System (Bauru5727): A large-scale test system was formed by linking two BIPS07 systems through a fictitious tie-line between the Bauru (440 kV) and Ivaiporã (765 kV) substations, belonging to the first and second BIPS07 subsystems, respectively. The first subsystem has all features of the original BIPS07, except for the PSS of Itaipu which was switched off. In the second subsystem, the HVDC link of Itaipu was replaced by an equivalent load at the Ibiuna (345 kV) substation. The resulting descriptor system has 40 366 variables, and . 5727 states, and two RHP poles 4) The 5723-State System (Juba5723): To highlight poles with small controllability and observability, we modified the AVR of the Juba power station, which has three units dispatched, totaling 30 MW. The AVR of the aggregated power station was replaced by a fast type AVR with very high gain. Moreover, the related PSS was given a low-frequency gain of high magnitude and wrong polarity. Due to the change in type of AVR the modified system has 5723 states. The descriptor system has 40 337 and variables. Besides the two previous RHP poles , there now exists a RHP pole which can only be observed in a minor subtransmission region of this huge test system, which has more than 130 000 MW of generation. B. Numerical Problems Due to Eigenvalues at Infinity To illustrate possible numerical difficulties due to eigenvalues at infinity, we apply SARQI to a 41-state power system model. The descriptor realization of this 41-state system has dimen. The power system stabilizer gain (Kpss) for sion one generator in this five-machine power system can be conof the descriptor matrix.2 Here we trolled by element , which causes a assigned a low value to this gain, right half-plane pole (see [40] for more details on this system). Since we are interested in the rightmost eigenvalues, we apply SARQI to the descriptor matrices with selection criterion “RightMost”. Fig. 1 shows the computed eigenvalues when the number of wanted eigenvalues is set to ten. Clearly, SARQI manages to find the rightmost eigenvalues. However, suppose we want to compute additional rightmost eigenvalues and set the number of wanted to 15. Due to presence 2Row 116 contains the algebraic equation , with variables “in” (for input to a gain block) and “out” (for output of the gain block) ordered as 115 and 116, respectively, in the augmented state vector.

935

Fig. 1. Selection of full spectrum of the 41-state power system matrices , where , computed by the QZ method , and the eigenvalues computed by SARQI (circles).

of eigenvalues of infinity, SARQI starts to compute approximations of these infinite eigenvalues:

As a result, the computational process becomes unstable and may result in unreliable approximate eigenvalues. Although for this example the problems become apparent only after ten eigenvalues have been found, in general one cannot predict when they occur: they might even occur during the first iterations. One way to avoid these problems is to apply SARQI to the state space matrix (instead of to the descriptor matrices), since the state space matrix has no eigenvalues at infinity. Indeed, as also described in detail in [31], this works well in practice. However, as we will see for larger systems in the following examples, the required implicit state space computations can become more expensive than the descriptor computation [note that computing the state space matrix explicitly is not attractive for large sys, since it is dense in general]. Furthermore, we tems show how our method can be used to compute eigenvalues with specific damping ratios. C. Large-Scale Power System Model Results Table I shows results for several large-scale power system models [3], [35], [40]–[42] and methods, for computing fifteen rightmost eigenvalues and eigenvalues with smallest damping ratio. We have considered SARQI with selection of eigenvalues with damping ratio (DR) closest to (equivalent to computing rightmost eigenvalues), for both descriptor and (implicit) state space computations, with initial shift . Furthermore, results for implicitly restarted Arnoldi with shift-and-invert and Cayley shifts are provided (using the ARPACK [43] in Matlab (

Authorized licensed use limited to: IEEE Xplore. Downloaded on April 21,2010 at 04:00:39 UTC from IEEE Xplore. Restrictions apply.

936

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 25, NO. 2, MAY 2010

eigs ), with shift strategy as described in [13]). The most important observations that can be made from Table I are as follows. • SARQI with damping ratio selection outperforms Arnoldi based methods on all fronts: SARQI is faster and computes more poles, including the rightmost ones. • SARQI is also able to deal with eigenvalues with multiplicity greater than one. Although, when computing 15 eigenvalues, not all occurrences of repeated eigenvalues are found, the result is better than for Arnoldi based for these methods, which fail with ARPACK error examples, due to (near) singularity of the system matrices. See Section V-D for more details on repeated eigenvalues. • From these results it is not clear upfront whether state space computations are to be preferred over descriptor state computations: apparently it depends on the system at hand, for instance, if fill-in minimizing reorderings can be used effectively. From a numerical point of view, state space computations are safer since all eigenvalues at infinity are eliminated explicitly. The potential additional costs are not significant and probably worth the increase of robustness. • The sequence in which eigenvalues are computed is not necessarily from right to left in the complex plane; see, for example, Figs. 2 and 3. However, it is clear that the rightmost poles are found. • We do not include here a comparison with Jacobi-Davidson methods [13], [22], [23], [28] because it can be shown that SARQI is mathematically equivalent to two-sided JacobiDavidson [15], [44], but simpler and faster since one does not need to solve the Jacobi-Davidson correction equations. Furthermore, in all our experiments we use exact factorizations) since these can be solves (directly via done efficiently for the sparse power system models and hence there is no need to resort to iterative linear system solvers. If exact solves are not feasible, however, it is advised to consider the Jacobi-Davidson method instead, because it can be shown to be more robust than Rayleigh quotient iteration when exact solves are not available (see, e.g., [13], [22], and [45]). • The numerical problems of implicitly restarted Arnoldi for the larger problems are caused by the fact that the system matrices are nearly singular. Hence, a rather large (real) shift is needed to execute the algorithms, which, however, prevents from convergence to the rightmost eigenvalues. Employing complex shifts, as in SARQI, solves the problems of singularity, but when using complex shifts for Cayley and shift-and-invert transformations, one also loses the desired transformation of the spectrum. SARQI with automatic shift selection does not suffer from this. Furthermore, while choosing the right Cayley shift is crucial for the convergence of the Arnoldi method to the rightmost eigenvalues, the influence of the initial shift for SARQI is typically only reflected in small changes to the order in which (rightmost) eigenvalues are found. D. Repeated Poles The proposed algorithm can also be applied to systems with repeated poles. To illustrate performance, we apply SARQI to

Fig. 2. Rightmost part of spectrum of the 5727-state power system state space ), computed by the QR matrix of Bauru5727 (descriptor order , and the eigenvalues computed by SARQI with damping ratio method (circles). The labels denote the order in which eigenvalues are computed by SARQI. All right half-plane eigenvalues have been found, as well as those with small damping ratio.

Fig. 3. Rightmost part of spectrum of the 5723-state power system state space ), computed by the QR method matrix of Juba5723 (descriptor order , and the eigenvalues computed by SARQI with damping ratio (circles). The labels denote the order in which eigenvalues are computed by SARQI. All right half-plane eigenvalues have been found, as well as those with small damping ratio.

the 2940-state system (Xingo6u). Note that usually, the six identical generator units, considered individually in this system, are agreggated as one equivalent unit. The system has 2940 states; the dimension of the descriptor matrices is . It can be shown that the corresponding system matrices have five nonde. fective repeated (complex conjugate) poles at to compute the We apply SARQI with damping ratio 25 rightmost eigenvalues (CPU times: 320 s for SARQI versus 471 s for the QR method).3 Fig. 4 shows the relevant part of the 3Note that SARQI was implemented in Matlab-code, while for QR the built-in and compiled eig routine was used.

Authorized licensed use limited to: IEEE Xplore. Downloaded on April 21,2010 at 04:00:39 UTC from IEEE Xplore. Restrictions apply.

ROMMES et al.: COMPUTING RIGHTMOST EIGENVALUES

937

ACKNOWLEDGMENT The authors would like to thank the anonymous referees for useful suggestions that helped to improve the readability of the paper. The authors also would like to thank J. ter Maten and T. Beelen (both from NXP Semiconductors) for comments on an earlier version of this paper. REFERENCES

Fig. 4. Rightmost part of spectrum of the 2940-state power system Xingo6u , and the eigenvalues comstate space matrix, computed by the QR method (circles). The labels puted by SARQI with damping ratio denote the order in which eigenvalues are computed by SARQI. All five occurare found. rences of the repeated (complex conjugate pair) pole at

spectrum, including the order in which the eigenvalues are computed by SARQI. The three most important observations here are that 1) all right half-plane eigenvalues have been found; 2) indeed all five repeated occurrences have been found; but 3) they are not found one after another, but as 4th, 9th, 16th, 19th, and 21st, respectively: other eigenvalues are found in between. Although the order is not of crucial importance in practice, future research will focus on specialized (block) algorithms for computing (both defective and nondefective) repeated poles and the corresponding eigenspaces, including ways to verify if indeed all repeated poles have been found. VI. CONCLUSIONS SARQI with the proposed selection strategy, is an effective, efficient, and robust method to compute the rightmost and poorly damped eigenvalues of large-scale power system models. Because SARQI operates on the sparse descriptor matrices of the system and computes the rightmost eigenvalues automatically, it is very useful in small-signal stability analysis of large-scale systems. Numerical experiments with realistic power system models confirmed that SARQI is a reliable and robust method for computing the rightmost and poorly damped eigenvalues. SARQI is faster than the QR/QZ method, and, moreover, is applicable to large-scale sparse descriptor systems as well. Compared to another popular method for computing a few specific eigenvalues, implicitly restarted Arnoldi with shift-and-invert/Cayley shifts, SARQI has the advantage that shift selection is done automatically, while for the former method one has to carefully provide shifts in order to find the desired eigenvalues. Without any adaptation, SARQI can be used for computing rightmost and poorly damped eigenvalues in several other engineering fields as well.

[1] G. H. Golub and C. F. van Loan, Matrix Computations, 3rd ed. Baltimore, MD: John Hopkins Univ. Press, 1996. [2] J. Ma, Z. Y. Dong, and P. Zhang, “Comparison of BR and QR eigenvalue algorithms for power system small signal stability analysis,” IEEE Trans. Power Syst., vol. 21, no. 4, pp. 1848–1855, Nov. 2006. [3] CIGRE TF 38.02.16, Impact of the Interactions Among Power System Controls, CIGRE, 2000, Tech. Rep. 166. [4] CEPEL, PacDyn User’s Manual, 8th ed. Rio de Janeiro: CEPEL—Brazilian Electrical Energy Research Center, 2008. [5] K. A. Cliffe, T. J. Garratt, and A. Spence, “Eigenvalues of block matrices arising from problems in fluid mechanics,” SIAM J. Matrix Anal. Appl., vol. 15, no. 4, pp. 1310–1318, Oct. 1994. [6] K. Meerbergen and A. Spence, “Implicitly restarted Arnoldi with purification for the shift-invert transformation,” Math. Comp., vol. 66, no. 218, pp. 667–689, 1997. [7] L. H. Bezerra and C. Tomei, “Spectral transformation algorithms for computing unstable modes,” Comp. Appl. Math, vol. 18, no. 1, pp. 1–15, 1999. [8] 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. [9] L. Wang, F. Howell, P. Kundur, C. Y. Chung, and W. Xu, “A tool for small-signal security assessment of power systems,” in Proc. Int. Conf. Power Ind. Comput. Appl., May 21–24, 2001, pp. 246–252. [10] G. Angelidis and A. Semlyen, “Improved methodologies for the calculation of critical eigenvalues in small signal stability analysis,” IEEE Trans. Power Syst., vol. 11, no. 3, pp. 1209–1217, Aug. 1996. [11] K. Meerbergen, A. Spence, and D. Roose, “Shift-invert and Cayley transforms for detection of rightmost eigenvalues of nonsymmetric matrices,” BIT, vol. 34, no. 3, pp. 409–423, 1994. [12] F. Ishikawa, H. Sasaki, J. Kubokawa, and H. Terasako, “An efficient algorithm for small signal stability analysis using the Arnoldi and S-matrix methods,” Elect. Eng. Jpn., vol. 121, no. 4, pp. 38–47, 1997. [13] J. Rommes, “Arnoldi and Jacobi-Davidson methods for generalized with singular ,” Math. Comput., eigenvalue problems vol. 77, no. 262, pp. 995–1015, Apr. 2008. [14] D. C. Sorensen, “Implicit application of polynomial filters in a -step Arnoldi method,” SIAM J. Matrix Anal. Appl., vol. 13, pp. 357–385, 1992. [15] J. Rommes, “Methods for eigenvalue problems with applications in model order reduction,” Ph.D. dissertation, Utrecht Univ., Utrecht, The Netherlands, 2007. [16] [Online]. Available: http://sites.google.com/site/rommes. [17] T. Kailath, Linear Systems. Englewood Cliffs, NJ: Prentice-Hall, 1980. [18] M. Green and D. J. N. Limebeer, Linear Robust Control. Englewood Cliffs, NJ: Prentice-Hall, 1995. [19] B. Pal and B. Choudhuri, Robust Control in Power Systems. New York: Springer, 2005. [20] S. Gomes, Jr., N. Martins, and C. Portela, “Computing small-signal stability boundaries for large-scale power systems,” IEEE Trans. Power Syst., vol. 18, no. 2, pp. 747–752, May 2003. [21] N. Martins and L. T. G. Lima, “Eigenvalue and frequency domain analysis of small-signal stability problems,” Eigenanalysis and frequency domain methods for system dynamic performance, no. TH0292–3-PWR, IEEE, pp. 17–33, 1990. [22] Z. Du, W. Liu, and W. Fang, “Calculation of rightmost eigenvalues in power systems using the Jacobi-Davidson method,” IEEE Trans. Power Syst., vol. 21, no. 1, pp. 234–239, Feb. 2006. [23] Z. Du, W. Liu, and W. Fang, “Calculation of electromechanical oscillation modes in large power systems using Jacobi-Davidson method,” IET Gen. Transm. Distrib., vol. 152, no. 6, pp. 913–918, 2005. [24] C. Lanczos, “An iteration method for the solution of the eigenvalue problem of linear differential and integral operators,” J. Res. Nat. Bur. Stand., vol. 45, no. 4, pp. 225–280, 1950.

Authorized licensed use limited to: IEEE Xplore. Downloaded on April 21,2010 at 04:00:39 UTC from IEEE Xplore. Restrictions apply.

938

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 25, NO. 2, MAY 2010

[25] B. N. Parlett, The Symmetric Eigenvalue Problem, Ser. Classics in Applied Mathematics. Philadelphia, PA: SIAM, 1998. [26] W. E. Arnoldi, “The principle of minimized iteration in the solution of the matrix eigenproblem,” Quart. Appl. Math., vol. 9, no. 1, pp. 17–29, 1951. [27] R. B. Lehoucq and D. C. Sorensen, “Deflation techniques within an implicitly restarted Arnoldi iteration,” SIAM J. Matrix Anal. Appl., vol. 17, pp. 789–821, 1996. [28] G. L. G. Sleijpen and H. A. van der Vorst, “A Jacobi-Davidson iteration method for linear eigenvalue problems,” SIAM J. Matrix Anal. Appl., vol. 17, no. 2, pp. 401–425, 1996. [29] [Online]. Available: http://www.win.tue.nl/casa/research/topics/jd/. [30] D. R. Fokkema, G. L. G. Sleijpen, and H. A. van der Vorst, “JacobiDavidson style QR and QZ algorithms for the reduction of matrix pencils,” SIAM J. Sci. Comp., vol. 20, no. 1, pp. 94–125, 1998. [31] J. Rommes and N. Martins, “Exploiting structure in large-scale electrical circuit and power system problems,” Lin. Alg. Appl., vol. 431, no. 3–4, pp. 318–333, 2009. [32] N. Martins, “Efficient eigenvalue and frequency response methods applied to power system small-signal stability studies,” IEEE Trans. Power Syst., vol. 1, no. 3, pp. 217–226, Aug. 1986. [33] M. Heinkenschloss, D. C. Sorensen, and K. Sun, “Balanced truncation model reduction for a class of descriptor systems with application to the Oseen equations,” SIAM J. Sci. Comput., vol. 30, no. 2, pp. 1038–1063, 2008. [34] 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. [35] J. Rommes and N. Martins, “Efficient computation of transfer function dominant poles using subspace acceleration,” IEEE Trans. Power Syst., vol. 21, no. 3, pp. 1218–1226, Aug. 2006. [36] B. N. Parlett, “The Rayleigh quotient iteration and some generalizations for nonnormal matrices,” Math. Comp., vol. 28, no. 127, pp. 679–693, Jul. 1974. [37] 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. [38] 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. [39] A. M. A. Hamdan and A. H. Nayfeh, “Measures of modal controllability and observability for first- and second-order linear systems,” J. Guid. Control Dynam., vol. 12, no. 3, pp. 421–428, May 1989. [40] J. Rommes and N. Martins, “Computing large-scale system eigenvalues most sensitive to parameter changes, with applications to power system small-signal stability,” IEEE Trans. Power Syst., vol. 23, no. 2, pp. 434–442, May 2008.

[41] N. Martins, A. de Andrade Barbosa, J. C. R. Ferraz, M. G. dos Santos, A. L. B. Bergamo, C. S. Yung, V. R. de Oliveira, and N. J. P. de Macedo, “Retuning stabilizers for the North-South Brazilian interconnection,” in Proc. IEEE Power Eng. Soc. Summer Meeting, Jul. 1999. [42] N. Martins, P. C. Pellanda, and J. Rommes, “Computation of transfer function dominant zeros with applications to oscillation damping control of large power systems,” IEEE Trans. Power Syst., vol. 22, no. 4, pp. 1218–1226, Nov. 2007. [43] R. H. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK SOFTWARE. [Online]. Available: http://www.caam.rice.edu/software/ARPACK/. [44] M. E. Hochstenbach and G. L. G. Sleijpen, “Two-sided and alternating Jacobi-Davidson,” Lin. Alg. Appl., vol. 358, no. 1–3, pp. 145–172, 2003. [45] M. E. Hochstenbach, “Subspace methods for eigenvalue problems,” Ph.D. dissertation, Utrecht Univ., Utrecht, The Netherlands, 2003.

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, in 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 (SM’91–F’98) received the B.Sc. 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 worked at CEPEL, Rio de Janeiro, Brazil, since 1978 in the development of computer tools for power system dynamics and control and specialized eigensolution methods.

Francisco Damasceno Freitas (M’94–SM’08) received the B.Sc. and M.Sc. degrees in electrical engineering from University of Brasilia, Brasilia, DF, Brazil, in 1985 and 1987, respectively, and the Ph.D. degree in electrical engineering from Federal University of Santa Catarina, Florianopolis, SC, Brazil, in 1995. Since 1986, he has been with the University of Brasilia, Brasilia, Brazil, 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.

Authorized licensed use limited to: IEEE Xplore. Downloaded on April 21,2010 at 04:00:39 UTC from IEEE Xplore. Restrictions apply.