Specialized eigenvalue methods for large-scale model order reduction problems Joost Rommes∗ NXP Semiconductors Corporate I&T/DTF High Tech Campus 37, WY4.042 5656AE, Eindhoven, The Netherlands [email protected]

Abstract Physical structures and processes are modeled by dynamical systems in many application areas, such as the design of very large-scale integration chips or large power systems. Since these dynamical systems can become very large, the essential simulation before production may consume hours or days of computing time. Hence there is need for efficient approaches that limit the computing time while preserving the accuracy. In this paper it will be shown how specialized eigenvalue methods and model order reduction techniques can be used to perform fast and accurate simulations of large dynamical systems. Results will be illustrated by numerical experiments with realistic examples.

1. Introduction Dynamical systems are used to model physical structures and processes in a wide range of application areas. The modeling of very large-scale integration (VLSI) chips, for instance, results in large systems of differential-algebraic equations [7]. The small signal stability assessment and controller design of large electrical power systems relies on sparse descriptor models linearized in various operating scenarios [8]. In structural engineering, large constructions such as buildings and bridges are modeled by dynamical systems to estimate, among others, the natural (resonance) frequencies [33]. It is often cheaper and safer to simulate the system before actually realizing it, because simulation gives insight in the behavior of the system and may help to check if design requirements are met. During the design of the system, several types of simulation and analysis may be required: in stability analysis one is interested in the right most eigenvalues [21], while ∗ This

author was supported by O-MOORE-NICE!, Marie Curie Actions FP6 MTKI-CT-2006-042477.

Nelson Martins CEPEL CEP-21944-970 Rio de Janeiro, Brazil [email protected]

in modal analysis one is interested in the resonance frequencies or dominant poles [19]. The increasing demand for complex components and large structures makes the models larger and more complicated. To be able to analyze and simulate these large-scale systems, there is need for efficient eigenvalue methods and reduced-order models of much smaller size, that approximate the behavior of the original model and preserve the important characteristics. In this paper methods for the solution of several specific large-scale eigenvalue problems are presented. It is shown how specialized eigenvalue methods can be used to compute the eigenvalues of interest in an efficient way. These eigenvalue methods can be used for several types of analysis, such as stability and sensitivity analysis, and for the tuning of (control) system behavior. Secondly, it is described how these eigenvalue methods can be used to improve existing model order reduction techniques for the efficient construction of accurate reduced order models. The operation and performance of the methods are illustrated by numerical experiments with real-life applications. The paper is organized as follows. In Section 2 dynamical systems, transfer functions, and related eigenvalue and model order reduction problems are described. In Section 3 several specialized eigenvalue methods are presented. Model order reduction methods are described in Section 4, where also a new approach for improving reduced order models is explained. In Section 5 numerical experiments with real-life applications are shown. Section 6 concludes.

2. Dynamical systems and model reduction 2.1

Dynamical systems

This paper is concerned with dynamical systems (E, A, B, C, D) of the form  ˙ E x(t) = Ax(t) + Bu(t) (1) y(t) = C ∗ x(t) + Du(t),

where A, E ∈ Rn×n , E can be singular, B ∈ Rn×m , C ∈ Rn×p , x(t) ∈ Rn is the state vector, u(t) ∈ Rm is the input vector, y(t) ∈ Rp is the output vector, and D ∈ Rp×m . The corresponding transfer function H : C −→ Cp×m is defined as

0 Response Dominant poles Dominant zeros

−10

−20

−1

H(s) = C (sE − A)

B + D,

(2)

with s ∈ C. If m = p = 1, the system (1) is called a single-input single-output (SISO) system, and b = B, c = C ∈ Rn are vectors and d = D ∈ R is a scalar. If m, p > 1, system (1) is called a multi-input multi-output (MIMO) system. M ∗ denotes the complex-conjugate transpose of a matrix M . The eigenvalues λi ∈ C of the matrix pencil (A, E) are the poles of transfer function (2). Assuming that the pencil is nondefective, the right and left eigenvectors xi and yi corresponding to finite eigenvalues λi can be scaled so that yi∗ Exi = 1. Furthermore, right and left eigenvectors corresponding to distinct eigenvalues are E-orthogonal: yi∗ Exj = 0 for i 6= j. The transfer function H(s) can be expressed as a sum of residue matrices Ri ∈ Cp×m over finite first-order poles [16]: H(s) =

n ˜ X i=1

Ri + D, s − λi

where the residues Ri are Ri = (C ∗ xi )(yi∗ B), n ˜ ≤ n is the number of finite first-order poles, and the contribution of poles at infinity is assumed to be zero.

2.2

Dominant poles

The following definition of modal dominance, related to modal model reduction [11], is used: Definition 2.1 A pole λi of H(s) with corresponding right and left eigenvectors xi and yi (yi∗ Exi = 1) is called the bi = kRi k2 /|Re(λi )| > R bj for all most dominant pole if R j 6= i.

Gain (dB)

−30

∗

−40

−50

−60

−70

−80 0

2

4

6

8

10

12

14

16

18

20

Frequency (rad/s)

Figure 1. Bode magnitude plot of test system (66 states, [19]). Stars (squares) denote dominant poles (zeros), with imaginary part on horizontal axis and real part (×100) on vertical axis.

2.3

Dominant zeros

In this paper only transmission zeros [32] are considered: Definition 2.2 A number z0 ∈ C is called a transmission zero if it satisfies rank(H(z0 )) < normalrank(H(s)), where normalrank(H(s)) [36] is the maximally possible rank of H(s) for at least one s ∈ C. In the SISO case, z0 ∈ C is called a transmission zero if H(z0 ) = cT (z0 E − A)−1 b + d = 0. A dominant zero causes a dip in the Bode plot of H(s), see also Fig. 1. The zeros are poles of the inverse transfer function H −1 (s), if it exists (Lemma 3.38 of [36]). Definition 2.3 A zero zi of an invertible transfer function H(s) is called dominant if it is a dominant pole of H −1 (s).

(c∗ xi )(yi∗ b) Ri = Re(λi ) Re(λi )

The following theorem defines a system Σz = (Ez , Az , Bz , Cz , Dz ), with transfer function Hz (s), such that Hz (s)H(s) = H(s)Hz (s) = I, where H(s) is the transfer function of the square system Σ = (E, A, B, C, D).

is dominant, i.e. a pole that is well observable and controllable in the transfer function. In a Bode plot of H(s), peaks occur at frequencies close to the imaginary parts of the dominant poles of H(s) (cf. Fig. 1).

Theorem 2.4 [20, Thm. 3.1] Let H(s) = C T (sE − A)−1 B + D be the transfer function of the descriptor realization Σ = (E, A, B, C, D), and let Hz (s) = CzT (sEz − Az )−1 Bz + Dz be the transfer function of the augmented

In the SISO case, a pole λi with large

descriptor realization Σz = (Ez , Az , Bz , Cz , Dz ), where     A B E 0 , , E = Az = z 0 0 CT D     0 0 , Cz = , Dz = 0, Bz = Im −Im and Im is the m × m identity matrix. Then Hz (s) = H −1 (s). For SISO systems, B and C are vectors b and c, respectively, and D is a scalar d. Consequently, Bz and Cz become bz = −en+1 and cz = en+1 , where en+1 is (n + 1)th elementary vector. For more details, see [20].

2.4

Model order reduction

The model order reduction problem is to find, given an n-th order (descriptor) dynamical system (E, A, B, C, D) e A, e B, e C, e D): (1), a k-th order system (E, ( ex ex(t) + Bu(t) e ˜˙ (t) = A˜ E (3) ∗ e ˜ (t) ˜ (t) + Du(t), y = C x

e A e ∈ Rk×k , B e ∈ Rk×m , C e ∈ Rk×p , where k < n, and E, k m p ˜ (t) ∈ R , u(t) ∈ R , y ˜ (t) ∈ R , and D ∈ Rp×m . The x number of inputs and outputs are the same as for the original system, and the corresponding transfer function becomes e e ∗ (sE e − A) e −1 B e + D. H(s) =C

The reduced-order model (3) should satisfy some or all of the following requirements [2]:

3

Specialized eigenvalue methods

In this section the Dominant Pole Algorithm (DPA) [19, 26] for computing dominant poles will be described and it will be shown how DPA can be adapted to compute dominant zeros [20], poles most sensitive to parameter changes [28], and how it can be generalized to MIMO [25] and higher order dynamical systems [27]. For brevity, all algorithms are described for SISO systems, but generalizing remarks can be found in Section 3.4.

3.1

Dominant poles

The poles of transfer function (2) are the λ ∈ C for which lims→λ |H(s)| = ∞, or, equivalently, the roots of G(s) = 1/H(s). The Dominant Pole Algorithm (DPA) [19] uses Newton’s method to find the dominant poles as the roots of G(s). This leads to the following Newton scheme: sk+1

1 H 2 (sk ) G(sk ) = s + k G′ (sk ) H(sk ) H ′ (sk ) ∗ c (sk E − A)−1 b . (4) = sk − ∗ c (sk E − A)−1 E(sk E − A)−1 b = sk −

Formula (4) was originally derived in [5]. Using vk = (sk E − A)−1 b and wk = (sk E − A)−∗ c, the Newton update (4) can also be written as the generalized two-sided Rayleigh quotient, provided wk∗ Evk 6= 0: sk+1

=

wk∗ Avk . wk∗ Evk

2. The order of the reduced system is much smaller than the order of the original system: k ≪ n.

An implementation of this Newton scheme, also known as the Dominant Pole Algorithm [19], is shown in Alg. 1. The two linear systems that need to be solved in step 3 and 4 of Alg. 1 can be efficiently solved using one LU factorization LU = sk E − A, by noting that U ∗ L∗ = (sk E − A)∗ . It will be assumed that an exact LU factorization is available. If an exact LU -factorization is not available, one has to use inexact Newton schemes [15, 24]. A detailed convergence analysis of DPA can be found in [29].

3. Preservation of (physical and numerical) properties such as stability and passivity.

3.2

1. The approximation error must be small, and a global error bound should exist. Usually this means that the ˜ (t)k should be minimized for output error ky(t) − y some or even all inputs u(t) in an appropriate norm.

4. The procedure must be stable and efficient. 5. The procedure must be based on some error tolerance. Simulation of the reduced order model should be (much) faster than of the original model: this may introduce additional requirements on the (sparsity) structure of the reduced order model. For an overview of model order reduction methods, see [2].

Dominant zeros

The dominant zeros of a transfer function H(s) can be computed by applying DPA (Alg. 1) to the inverse transfer function, as defined in Section 2.3. The procedure can be summarized as follows: 1. Given the system Σ = (E, A, b, c, d) with transfer function H(s), define the realization Σz = (Ez , Az , bz , cz , dz ) of the inverse transfer function [H(s)]−1 as described in Theorem 2.4.

Algorithm 1 The Dominant Pole Algorithm (DPA). INPUT: System (A, E, b, c), pole estimate s0 , tol ǫ ≪ 1 OUTPUT: Dominant pole λ and corresponding right and left eigenvectors x and y 1: Set k = 0 2: while not converged do 3: Solve vk ∈ Cn from (sk E − A)xk = b 4: Solve wk ∈ Cn from (sk E − A)∗ yk = c 5: Compute the new pole estimate sk+1 = sk − 6:

c∗ vk w∗ Avk = ∗k ∗ wk Evk wk Evk

The pole λ = sk+1 with x = vk /kvk k2 and y = wk /kwk k2 has converged if kAx − sk+1 Exk2 < ǫ

7: 8:

Set k = k + 1 end while

where b, c ∈ Rn are vectors. Then the sensitivity of an eigenvalue λ with left and right eigenvectors y and x (with y∗ Ex = 1) becomes ∂A ∂λ = y∗ x = (y∗ b)(c∗ x) = (c∗ x)(y∗ b). ∂p ∂p

In the right-hand side of (8) one recognizes the residues of the transfer function H(s) = c∗ (sE − A)−1 b. Consequently, the most sensitive eigenvalues of the pencil (A(p), E) can be computed by applying DPA to (E, A, b, c), with b and c defined by (7). If ∂A ∂p has rank higher than 1, one can change Alg. 1 as follows to compute the most  sensitiveeigenvalues: replace ∗

∂A b and c by ∂A , respectively. For ∂p vk−1 and ∂p wk−1 more details and generalizations to higher rank derivatives and multiparameter systems, see [28].

3.4

2. Apply DPA (Alg. 1) to Σz = (Ez , Az , bz , cz , dz ) to compute dominant zeros of H(s).

(8)

MIMO and higher order systems

The main idea of using Newton’s method to find dominant poles can be generalized to higher order systems [27]. For the second-order transfer function H(s) = c∗ (s2 M + sC + K)−1 b, for instance, the scheme becomes

For a more detailed explanation, see [20].

3.3

sk+1

Eigenvalue sensitivity

Let p ∈ R be a system parameter (e.g., a resistor value in a electric circuit), and let A(p) and E(p) be matrices that depend on p. The derivative of an eigenvalue λ of the pencil (A(p), E(p)), with left and right eigenvectors y ≡ y(p) and x ≡ x(p), to a parameter p is given by [22, 13] ∂E y∗ ( ∂A ∂λ ∂p − λ ∂p )x = . ∂p y∗ Ex

(5)

The derivative (5) is often called the sensitivity (coefficient) of λ. Assuming that ∂E ∂p = 0, with y and x scaled so that y∗ Ex = 1, the eigenvalue derivative (5) becomes ∂A ∂λ = y∗ x. ∂p ∂p

(6)

The larger the magnitude of the derivative (6), the more sensitive eigenvalue λ is to changes in parameter p. In practical applications such information is useful when, for instance, a system needs to be stabilized by moving poles from the right half-plane to the left half-plane [28, 35]. Suppose that the derivative of A to parameter p has rank 1 and hence can be written as ∂A = bc∗ , ∂p

(7)

= sk −

c∗ v , w∗ (2sk M + C)v

where v = (s2k M +sk C +K)−1 b and w = (s2k M +sk C + K)−∗ c. For generalizations to MIMO systems, see [25].

4

Model order reduction methods

4.1

Projection based methods

Projection based model order reduction methods construct a reduced order model via Petrov-Galerkin projection e A, e B, e C, e D) ≡ (W ∗ EV, W ∗ AV, W ∗ B, V ∗ C, D), (E,

where V, W ∈ Rn×k are matrices whose k ≪ n columns form bases for relevant subspaces of the state-space. There are three types of projection methods that differ in the way the matrices V and W are chosen: • Modal truncation: the columns of V and W are eigenvectors of (A, E), see for instance [9, 19, 24]. A general framework for modal approximation is: 1. Compute the poles λi and corresponding left and right eigenvectors yi and xi 2. Sort (λi , Ri ) in decreasing |Ri |/Re(λi ) order

3. Truncate at |Ri |/Re(λi ) < Rmin

4. Project with Wk = [y1 , . . . , yk ] and Vk = [x1 , . . . , xk ], with Wk∗ EVk = I and Wk∗ AVk = diag(λ1 , . . . , λk ) ≡ Λk to obtain  ˜ ˜˙ ˜ (t) + bu(t) x = Λk x ∗ ˜ x ˜ (t) y(t) = c

−35

−45

Gain (dB)

−50

with transfer function Hk (s) =

k X i=1

• Moment matching: the columns of V and W are bases for Krylov spaces, see for instance [12, 23]. The series expansion of H(s) = c∗ (sE − A)−1 b around s0 is H(s)

=

i=0

∗

i

−55 −60 −65

Ri s − λi

For large-scale systems step one of this framework is in general not feasible, since full space eigenmethods such as QR and QZ [10] have time complexity O(n3 ). In practical situations an alternative for step 1 is the Subspace Accelerated DPA (SADPA) [26]: given a k ≪ n, SADPA computes k dominant poles and corresponding eigenvectors in an iterative way.

∞ X

SADPA (k=12) PRIMA (k=30) ORIG (n=66)

−40

mi (s − s0 )i , −1

where mi = c G (s0 E − A) b (with G = (s0 E − A)−1 E) are called the moments of H(s). The idea of model order reduction via moment matching is to match only 2k ≪ n moments. One way to do this is by constructing specific Krylov subspaces. A Krylov subspace is defined as Kk (A, v) = span(v, Av, . . . , Ak−1 v). Computing the basis vectors directly is numerically unstable. There are, however, several stable numerical schemes, for instance the Arnoldi method [3] (see [31] for more methods). The general framework for moment matching via Krylov methods is 1. Compute bases V ∈ Rn×k and W ∈ Rn×k for Kk ((s0 E − A)−1 E, (s0 E − A)−1 b) and k −∗ ∗ −∗ K ((s0 E − A) E , (s0 E − A) c), for instance using dual Arnoldi method [12]. 2. Use Petrov-Galerkin projection to obtain k-th order reduced order model:  ˜˙ = (W ∗ AV )˜ x(t)  (W ∗ EV )x + (W ∗ b)u(t)  y˜(t) = (c∗ V )˜ x(t)

• Balanced truncation: V and W are part of a balancing transformation, see for instance [11].

−70 −75 −80 0

2

4

6

8

10

12

14

16

18

20

Frequency (rad/s)

Figure 2. Frequency response of complete system (n = 66), modal approximation (k = 12), and dual Arnoldi model (k = 30).

4.2

Using dominant poles to improve Krylov models

It is well known that for some examples moment matching works well, while reduced order models computed by modal approximation are of low quality, and the other way around [24, 2]. Generally speaking, modal approximation performs best if there are k ≪ n dominant poles with residues much larger than the residues of the non-dominant poles. In other words, there is a k ≪ n for which one has |R1 | ≥ |R2 | ≥ . . . ≥ |Rk | ≫ |Rk+1 | ≥ |Rn−1 | ≥ |Rn |, so that truncation at the kth pole does not give a large error [11]. Moment matching based approaches, on the other hand, perform best if the moments show a similar steep decay. There is, however, one additional complication for Krylov based moment matching approaches, that is best explained by an example. In Figure 2 the Bode magnitude plots of an exact transfer function and two reduced order models are shown: one modal approximation and a moment matching approximation. While the modal approximation captures the dominant dynamics, the moment matching model deviates for ω > 4 rad/s. The modal approximation matches the original transfer function well because it is built from the 7 most dominant poles. The moment matching Arnoldi model (k = 30) was built by building left and right Krylov subspaces with shift s0 = 0. Therefore, the match for frequencies up to ω = 4 rad/s is good. For higher frequencies, however, this approach suffers from a well known property of Arnoldi methods, that were originally developed for the computation of eigenvalues: the eigenvalue approximations, or Ritz values, tend to approximate the eigenvalues at the outside of the spectrum [34]. This can also be seen in Fig. 3, where the cir-

0

Exact

exact poles SADPA (k=12) Dual Arnoldi (k=30)

0.5 0.4

−50

Reduced model (k=89) Error (orig − reduced)

−100

0.3

Gain (dB)

0.2

imag

0.1 0

−150

−200

−250

−0.1 −0.2

−300 −0.3 −0.4

−350

−0.5 −1

−0.9

−0.8

−0.7

−0.6

−0.5

−0.4 real

−0.3

−0.2

−0.1

0

−400 −1 10

0.1

2

10

3

10

4

10

5

10

Figure 4. Bode plot of modal equivalent, complete model and error for transfer function of the PEEC model.

5

1. Compute k ≪ n dominant poles λi = αi ± βi j, with √ i = 1, . . . k and j = −1.

1

10

Frequency (rad/s)

Figure 3. Relevant part of pole spectrum of complete system (n = 66), modal approximation (k = 12), and dual Arnoldi model (k = 30).

cles denote the poles of the moment matching model (note the inverses of the poles are shown): they match the eigenvalues at the outside. The dominant poles, however, may be located anywhere in the spectrum, as can also be seen in Fig. 3 (squares). This explains why the Arnoldi model fails to capture the peaks. Based on the above observations and theory in [12], the idea is to use the imaginary parts of dominant poles as shifts for the rational Krylov approach, so that resonance peaks located well within the system frequency bandwidth can also be captured by Krylov methods. The combined dominant pole – Krylov approach can be summarized as follows:

0

10

5.1

Numerical experiments Circuit simulation

This descriptor model arises from a partial element equivalent circuit (PEEC) of a patch antenna structure and the dimension of the matrices A and E is n = 480 (see [6] for more details and the model data). The system is known as a difficult problem, because it has pole clusters in the interior of the spectrum, and Krylov based methods may suffer from this [14]. Figure 4 shows the Bode plot of the complete model and the reduced model (45 poles, 89 states) computed by SADPA. The approximation is almost exact for the frequency range [0.1, . . . , 102 ] rad/sec. It was not possible to obtain reduced order models of the same order and quality with moment matching based methods.

2. Choose interpolation points si = βi j.

5.2

3. Construct Vi , Wi ∈ Cn×ki such that their columns are bases for the rational Krylov [30] spaces

The Brazilian Interconnected Power System (BIPS) is a year 1999 planning model that has been used in practice (see [26] for more technical details). The size of the sparse matrices A and E is n = 13, 251. To study the sensitivity of the dominant poles and system stability, the gain (Kpss) of one of the generators is varied between 0 and 30, with increments of 0.5. Figure 5 shows the traces for the most sensitive poles as computed by SASPA (Section 3.3, see also [28]). The CPU time for the 60 runs was 1450 s. A root-locus plot for all poles, computed using the QR method, confirmed that the most sensitive poles were found, but needed 33600 s. Hence, for large-scale systems, SASPA is a very effective way to produce relevant root-locus plots.

colspan(Vi ) = Kki ((si E − A)−1 E, (si E − A)−1 b) and colspan(Wi ) = Kki ((si E − A)−∗ E ∗ , (si E − A)−∗ c), respectively. 4. Project with V = orth([V1 , . . . , Vk ]) and W = orth([W1 , . . . , Wk ]), where orth constructs an orthogonal basis for the spaces. The size of the reduced model Pk is at most K = i=1 ki , matching 2K moments.

Power systems

−110

15

k=30 (QDPA) k=35 (QDPA) k=40 (PRIMA) Exact

−120

−130

10 Gain (dB)

imag (rad/s)

−140

−150

−160 5 −170

−180

−190 4 10 0 −1.4

−1.2

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

5

10 Frequency (Hz)

6

10

real (1/s)

Figure 6. Bode magnitude plots of exact transfer function (solid), 40th order Krylov model (dash), and 35th (dash-dot) and 30th (dot) order modal equivalents.

Figure 5. Root locus plot of sensitive poles computed by SASPA. As the gain increases, the critical rightmost pole crosses the imaginary axis and the 5% damping ratio boundary. Squares denote initial pole locations.

5.3

Micro Electro Mechanical Systems

The Butterfly gyro is a vibrating micro-mechanical gyro that is used in inertial navigation applications (see [1, 18] for more details). Model order reduction not only reduces the simulation times significantly, it also helps in developing signal processing algorithms for the gyro. The full system has 17361 degrees of freedom. In Figure 6 the Bode magnitude plot of a 40th order rational Krylov model is plotted together with 30th and 35th order modal equivalents. In the eye norm it appears as if the 35th and 40th order model are qualitatively the same, except near 106 Hz, where the 35th order modal equivalent is more accurate. Closer inspection shows that the 40th order model is more accurate in the lower frequency range, while the accuracy decreases after 105 Hz. This can be explained by the fact that the expansion point was σ0 = 0. The modal equivalents are accurate in the neighborhood of frequencies near the imaginary parts of the dominant poles.

5.4

The breathing sphere

Fig. 7 shows the frequency response of a 70th order Second-Order Arnoldi [4] reduced model of vibrating body from sound radiation analysis (n = 17611 degrees of freedom, see [17]), that was computed using the complex part iβ of dominant poles λ = α + iβ (computed by Quadratic DPA [27]) as interpolation points. This model is more accurate than reduced order models based on standard Krylov methods and matches the peaks up to ω = 1 rad/s, because of use of shifts near the resonance frequencies. This model

is a good example of the combined dominant pole – rational Krylov approach, since modal approximations of similar quality require too much CPU time, while Krylov models with uniformly spaced shifts do not capture the peaks.

6

Conclusions

A family of specialized eigenvalue methods was presented, that can be used to solve a number of different largescale eigenvalue problems arising in real-life applications and simulation of dynamical systems: computation of dominant poles, computation of eigenvalues most sensitive to parameter changes, and construction of reduced-order models. All methods are based on the Dominant Pole Algorithm [19] and can be implemented and generalized easily. Numerical experiments with real-life examples show good performance. Furthermore, the algorithms can be used to improve existing model order reduction techniques.

References [1] Oberwolfach Model Reduction Benchmark Collection. http://www.imtek.uni-freiburg.de/simulation/benchmark/. [2] A. C. Antoulas. Approximation of Large-Scale Dynamical Systems. SIAM, 2005. [3] W. E. Arnoldi. The principle of minimized iteration in the solution of the matrix eigenproblem. Quart. Appl. Math., 9(1):17–29, 1951. [4] Z. Bai and Y. Su. SOAR: a second-order Arnoldi method for the solution of the quadratic eigenvalue problem. SIAM J. Matrix Anal. Appl., 26(3):640–659, 2005. [5] L. H. Bezerra. An eigenvalue method for calculating dominant poles of a transfer function. Applied Mathematics Letters, 21:244–247, Feb 2008.

50

0

Gain (dB)

−50

−100

−150

−200

−250 −1 10

k=70 (RKA) Exact Rel Error 0

10

1

10

Frequency (rad/sec)

Figure 7. Exact and reduced system transfer functions for a vibrating body, computed by a rational Krylov method with resonance frequencies as complex interpolation points.

[6] Y. Chahlaoui and P. van Dooren. A collection of Benchmark examples for model reduction of linear time invariant dynamical systems. SLICOT Working Note 2002-2, 2002. [7] P. G. Ciarlet, W. H. A. Schilders, and E. J. W. ter Maten, editors. Numerical Methods in Electromagnetics, volume 13 of Handbook of Numerical Analysis. 2005. [8] CIGRE TF 38.01.17. Analysis and control of power system oscillations. Technical Report 111, CIGRE, December 1996. [9] E. Davison. A method for simplifying linear dynamic systems. IEEE Trans. Aut. Control, 11(1):93–101, 1966. [10] G. H. Golub and C. F. van Loan. Matrix Computations. John Hopkins University Press, third edition, 1996. [11] M. Green and D. J. N. Limebeer. Linear Robust Control. Prentice-Hall, 1995. [12] E. J. Grimme. Krylov projection methods for model reduction. PhD thesis, University of Illinois, 1997. [13] S. B. Haley. The generalized eigenproblem: pole-zero computation. Proc. IEEE, 76(2):103–120, February 1988. [14] P. J. Heres. Robust and efficient Krylov subspace methods for model order reduction. PhD thesis, Eindhoven University of Technology, 2005. [15] M. E. Hochstenbach and G. L. G. Sleijpen. Two-sided and alternating Jacobi-Davidson. Lin. Alg. Appl., 358(1-3):145– 172, 2003. [16] T. Kailath. Linear Systems. Prentice-Hall, 1980. [17] J. Lampe and H. Voss. Second order Arnoldi reduction: application to some engineering problems. Report 93, Hamburg University of Technology, 2006. [18] J. Lienemann, D. Billger, E. B. Rudnyi, A. Greiner, and J. G. Korvink. MEMS compact modeling meets model order reduction: examples of the application of Arnoldi method to microsystem devices. In Technical Proceedings of the 2004 Nanotechnology Conference and Trade Show, 2004. [19] N. Martins, L. T. G. Lima, and H. J. C. P. Pinto. Computing dominant poles of power system transfer functions. IEEE Trans. Power Syst., 11(1):162–170, Feb 1996.

[20] 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., 22(4):1218–1226, Nov. 2007. [21] K. Meerbergen and A. Spence. Implicitly restarted Arnoldi with purification for the shift-invert transformation. Math. Comp., 66(218):667–689, 1997. [22] D. V. Murthy and R. T. Haftka. Derivatives of eigenvalues and eigenvectors of a general complex matrix. Int. J. Num. Meth. Eng., 26:293–311, January 1988. [23] A. Odabasioglu and M. Celik. PRIMA: Passive Reducedorder Interconnect Macromodeling Algorithm. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 17(8):645–654, 1998. [24] J. Rommes. Methods for eigenvalue problems with applications in model order reduction. PhD thesis, Utrecht University, 2007. [25] J. Rommes and N. Martins. Efficient computation of multivariable transfer function dominant poles using subspace acceleration. IEEE Trans. Power Syst., 21(4):1471–1483, Nov 2006. [26] J. Rommes and N. Martins. Efficient computation of transfer function dominant poles using subspace acceleration. IEEE Trans. Power Syst., 21(3):1218–1226, Aug. 2006. [27] J. Rommes and N. Martins. Efficient computation of transfer function dominant poles of large second-order dynamical systems. Preprint 1360, Utrecht University, 2007. Accepted for publication in SISC. [28] 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., 2008. Accepted for publication. [29] J. Rommes and G. L. G. Sleijpen. Convergence of the dominant pole algorithm and Rayleigh quotient iteration. SIAM J. Matrix Anal. Appl., 30(1):346–363, 2008. [30] A. Ruhe. Rational Krylov: A practical algorithm for large sparse nonsymmetric matrix pencils. SIAM J. Sc. Comp., 19(5):1535–1551, 1998. [31] Y. Saad. Numerical Methods for Large Eigenvalue Problems: Theory and Algorithms. Manchester University Press, 1992. [32] C. B. Schrader and M. K. Sain. Research on system zeros: a survey. Int. J. Control, 50(4):1407–1433, 1989. [33] F. Tisseur and K. Meerbergen. The quadratic eigenvalue problem. SIAM Review, 43(2):235–286, 2001. [34] H. A. van der Vorst. Computational Methods for Large Eigenvalue Problems. In P.G. Ciarlet and J.L. Lions, editor, Handbook of Numerical Analysis, volume VIII, pages 3–179. North-Holland (Elsevier), 2001. [35] 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., pages 246– 252, May 21–24 2001. [36] K. Zhou, J. Doyle, and K. Glover. Robust and Optimal Control. Prenctice Hall, 1996.