Comparison of Numerical Methods and Open-Source Libraries for Eigenvalue Analysis of Large-Scale Power Systems

: This paper discusses the numerical solution of the generalized non-Hermitian eigenvalue problem. It provides a comprehensive comparison of existing algorithms, as well as of available free and open-source software tools, which are suitable for the solution of the eigenvalue problems that arise in the stability analysis of electric power systems. The paper focuses, in particular, on methods and software libraries that are able to handle the large-scale, non-symmetric matrices that arise in power system eigenvalue problems. These kinds of eigenvalue problems are particularly difﬁcult for most numerical methods to handle. Thus, a review and fair comparison of existing algorithms and software tools is a valuable contribution for researchers and practitioners that are interested in power system dynamic analysis. The scalability and performance of the algorithms and libraries are duly discussed through case studies based on real-world electrical power networks. These are a model of the All-Island Irish Transmission System with 8640 variables; and, a model of the European Network of Transmission System Operators for Electricity, with 146,164 variables.


Introduction
Eigenvalue analysis is a fundamental component of electric power system Small-Signal Stability Analysis (SSSA). It is utilized in order to determine the stability of the operating points of the system [1,2], design controllers [3], determine machine clusters and reduced equivalent models [4], and evaluate the properties of the transmission grid [5]. On the other hand, there is rich literature on numerical algorithms that determine the full or a partial solution of a given eigenvalue problem. Relevant monographs on the topic are, for example [6,7]. However, not all available algorithms are suitable for the SSSA of power systems.
The vast majority of numerical algorithms, in fact, solve exclusively symmetric, Hermitian, or tridiagonal Hermitian eigenvalue problems, which arise in many engineering, physics, and computer applications. Such algorithms are, for example, the Davidson [8], the implicitly restarted Lanczos [9], and the locally optimal block preconditioned conjugate gradient [10] methods. However, the matrices that describe a linearized power system model are typically non-symmetric. When compared to symmetric problems, non-symmetric eigenvalue problems are more difficult and computationally demanding to solve. This is also due to the fact that symmetric eigenvalue problems are generally well-conditioned, while non-symmetric eigenvalue problems are not. Consequently, the available free and open-source software libraries that solve non-symmetric eigenvalue problems are a small subset of all existing libraries.
The scalability of the numerical solution of eigenvalue problems is also very important, since real-world electrical power networks are large-scale dynamic systems. Unfortunately, the most reliable methods for finding the full spectrum of an eigenvalue problem are dense-matrix methods, and their computational complexity and memory requirements increase more than quadratically (in some cases even cubically) as the size of the matrix increases. As a matter of fact, the largest ever eigenvalue analysis with a dense algorithm to date was the solution of a 10 6 × 10 6 problem in about 1 h, and it was carried out in 2014 by the Japanese K computer in Riken. To be able to obtain this result, the K computer includes 88,000 processors that draw a peak power of 12.6 MW, while its operation costs annually US$10 million. The solution of large-scale power system eigenvalue problems is challenging even when using sparse matrices and limiting the search to a subset of the spectrum.
A coarse taxonomy of existing algorithms for the solution of non-symmetric eigenvalue problems is as follows: • Vector iteration methods, which, in turn, are separated to single and simultaneous vector iteration methods. Single vector iteration methods include the power method [11] and its variants, such as the inverse power and Rayleigh quotient iteration. Simultaneous vector iteration methods include the subspace iteration method [12] and its variants, such as the inverse subspace method.
Regarding the application of vector iteration methods in power systems, the inverse power iteration was first discussed in [13]. A Rayleigh quotient iteration based algorithm was proposed in [14,15], while recent studies have provided subspace accelerated and deflated versions of the same algorithm [16][17][18]. Finally, simultaneous vector iteration methods were studied in [19,20]. • Schur decomposition methods, which mainly include the QR algorithm [21], the QZ algorithm [22], and their variants, such as the QR algorithm with shifts. Schur decomposition based methods have been the standard methods employed for the eigenvalue analysis of small to medium size power systems [23,24]. • Krylov subspace methods, which basically include the Arnoldi iteration [25] and its variants, such as the implicitly restarted Arnoldi [26] and the Krylov-Schur method [27]. In this category also belong preconditioned extensions of the Lanczos algorithm, such as the non-symmetric versions of the Generalized Davidson and the Jacobi-Davidson method. In power systems, different versions of the Arnoldi method were proposed in [19,28,29], while a parallel version of the Krylov-Schur method was implemented very recently, in [30]. The Jacobi-Davidson method and an inexact version of the same method were discussed in [31,32], respectively. • Contour integration methods, which basically include a moment-based Hankel method [33] and a Rayleigh-Ritz-based projection method [34] proposed by Sakurai; and, the FEAST algorithm [35]. Contour integration for power system eigenvalue analysis was only very recently discussed [36,37].
The main objective of this paper is to provide a state-of-art reference of eigenvalue numerical algorithms that can be actually utilized for the analysis of large power system models. To this aim, we have carefully selected only methods for non-Hermitian matrices for which open-source numerical libraries are available and that we were able to test and that were successful to solve relatively "large" eigenvalue problems. In particular, the paper provides a fair comparison of algorithms that have been proposed for power system analysis only very recently, such as versions of the parallel Krylov-Schur and contour integration methods discussed in [30,37], as well as of state-of-art implementations of standard algorithms for unsymmetric dense and sparse matrices, such as the QR and Arnoldi method. The paper provides an overview and comparison of the state-of-art open-source libraries that solve the non-symmetric eigenvalue problem. These are Linear Algebra PACKage (LAPACK) [38] and its GPU-based version, namely Matrix Algebra for GPU and Multicore Architectures (MAGMA) [39]; ARnoldi PACKage (ARPACK) [40]; Anasazi [41]; Scalable Library for Eigenvalue Problem computations (SLEPc) [42]; FEAST [43]; and, z-PARES [44].
The remainder of the paper is organized, as follows. Section 2 provides the background and problem formulation of power system eigenvalue analysis and a description of matrix pencil spectral transforms. Section 3 outlines the most relevant numerical algorithms for the eigenvalue analysis of power systems. Open-source libraries that implement state-of-art implementations of these methods are presented in Section 4. The examined algorithms are comprehensively tested through two real-world power system models in the case studies discussed in Section 5. Finally, conclusions are drawn in Section 6.

Power System Model
Electric power system models for transient stability analysis can be formulated as a set of non-linear Differential Algebraic Equations (DAEs), as follows [45]: where f : R n+m → R n , g : R n+m → R m ; x = x(t), x ∈ R n , are the state variables, y = y(t), y ∈ R m , are the algebraic variables; F ∈ R n×n , G ∈ R m×n , are assumed to be constant matrices; and t ∈ [0, ∞) is the simulation time. Finally, 0 n,m denotes the zero matrix of dimensions n × m. For the purpose of analysis and, for sufficiently small disturbances, (1) can be linearized around an equilibrium (x o , y o ), as follows: where ∆x = x − x o , ∆y = y − y o ; f x , f y , g x , g y , are the Jacobian matrices calculated at (x o , y o ). This system can be rewritten in the following form: where z = ∆x ∆y , E I = F 0 n,m G 0 m,m , A I = f x f y g x g y .
the linear system (2) is often reduced to a system of Ordinary Differential Equations (ODEs), by elimination of algebraic variables. The linearized ODE power system model can be described, as follows: where A S = (F − f y g −1 y G) −1 ( f x − f y g −1 y g x ) is the state matrix, and where we have assumed that g y , F − f y g −1 y G are non-singular.

Eigenvalue Problem
In this paper, we are concerned with finding the eigenvalues of the following system: where E, A ∈ R r×r , x : [0, +∞] → C r×1 . Matrix E in (5) can be: • non-singular, i.e., det(E) = 0. This is the case of the ODE power system model (4). In particular, (4) can be obtained from (5) for r = n, x ≡ ∆x, A ≡ A S , E ≡ I r . • singular, i.e., det(E) = 0. This is the case of the DAE power system model (2). In particular, (2) can be obtained from (5) for r = n + m, x ≡ z, A ≡ A I , E ≡ E I .
The stability of (5) can be assessed by calculating its eigenvalues, which are defined as the roots of the characteristic equation: where s ∈ C denotes the complex Laplace variable. In (6), sE − A is called the matrix pencil and det(sE − A) the characteristic polynomial of system (5). In this paper, we assume that det(sE − A = Π(s) ≡ 0, where Π(s) is polynomial of order equal or less than r, in which case sE − A is called a regular pencil [46]. In general, analytical solution of (6) is possible only if r ≤ 4. For higher degrees, general formulas do not exist and only the application of a numerical method is possible. In addition, algorithms that explicitly determine the characteristic polynomial det(sE − A and then numerically calculate its roots, may be extremely slow, even for small problems. Alternatively, the eigenvalues of sE − A can be found from the solution of the Generalized Eigenvalue Problem (GEP): where u ∈ C r×1 and w ∈ C 1×r . Every value of s that satisfies (7) is an eigenvalue of the pencil sE − A, with the vectors u, w being the corresponding right and left eigenvectors, respectively. Thus, the solution of the GEP consists in calculating the eigenpairs, i.e., eigenvalues and eigenvectors, which satisfy (7). It may be required that only right (or left) or both right and left eigenvectors are calculated, depending on the analysis that needs to be carried out. In general, the pencil sE − A has rank(sE − A) finite eigenvalues and the infinite eigenvalue with multiplicity r − rank(sE − A). Note that, if E is singular, then the pencil will have the infinite eigenvalue with multiplicity of at least one.
In the special case that the left-hand side matrix of (5) is the identity matrix, e.g., (4), the general problem (7) is reduced to the following Linear Eigenvalue Problem (LEP): the solution of the LEP consists in calculating the r finite eigenvalues and eigenvectors of sI r − A.

Spectral Transforms
In general, the solution of the eigenvalue problem involves finding the full or partial spectrum of the pencil sE − A. However, depending on the applied numerical method as well as on the structure of the system matrices, it is common that the eigenvalues are not found by directly using sE − A, but through the pencil that arises from the application of a proper spectral transform. Spectral transforms are utilized by eigenvalue numerical methods usually for one of the following reasons: We describe here the Möbius transformation, which is a general variable transformation that includes as special cases all spectral transforms used in practice by eigenvalue algorithms. The formulation of the Möbius transformation is: The restriction in (9) is necessary, because, if ad = bc, then s is constant which is not possible. Applying the transform (9) in (6) we obtain or, equivalently, by using determinant properties which is the Characteristic equation characteristic equation of a linear dynamical system with pencil z(aE − cA) − (dA − bE). System (5) will be referred as the prime system, and the family of systems (10) will be defined as the proper "M-systems". An important property is that the solutions and stability properties of system (5) can be studied through (10) without resorting to any further computations, see [47]. The utilities of the family of systems of type (10) have been further emphasized by the features of some particular special cases. Table 1 summarizes the most commonly employed Möbius transforms and the corresponding matrix pencils for the GEP. The values of the parameters a, b, c, d that lead to each of these transforms are given in Table 2. In the case that σ > 0, the Cayley transform is equivalent to the bilinear transform z := ( T 2 s + 1)/( T 2 s − 1), where T = 2 σ . The choice of the best transform for a specific system and eigenvalue problem is a challenging task to solve, since the selection of shift values is always more or less heuristic. Finally, note that the importance of such spectral transforms for the eigenvalue analysis has also been identified in power system literature, see [48,49].

Description of Numerical Algorithms
This section provides an overview of the following classes of eigenvalue numerical methods: vector iteration methods, Schur decomposition methods, Krylov subspace methods, and contour integration methods. These are, in fact, the methods implemented by the open-source software libraries that are compared in Section 4 and in the case studies of Section 5 of this paper.

Single Vector Iteration Methods
The vector iteration or power method is the oldest and probably the simplest and most intuitive numerical method for solving an eigenvalue problem. Consider the LEP (8) with matrix pencil sI r − A. The main idea of the power method is that, if one starts with a vector b and repeatedly multiplies by A, then the subsequence converges to a multiple of the right eigenvector that corresponds to the eigenvalue of A with LM. More strictly put, under the assumption that the magnitude of the largest eigenvalue λ m is larger than that of all other eigenvalues, (11) linearly converges to a multiple of the corresponding right eigenvector u m . In order to guarantee the convergence, it is also required that the initial vector is selected with a non-zero component in the direction of λ m . Applying the above idea, the method starts with an initial vector b (0) , which is updated at each step of the iteration through multiplication with A and normalization. The k-th step of the iteration is as follows: Finally, under the assumptions discussed above, b (k) converges to the vector u m , which corresponds to the eigenvalue of LM λ m . Given u m , an estimation of the eigenvalue λ m is given by the Rayleigh quotient. In general, given a right eigenvector u i , the Rayleigh quotient provides an estimate of the corresponding eigenvalue λ i , as follows:λ where u H i is the conjugate transpose of u i . The power method finds an eigenvector that corresponds to the eigenvalue of LM. However, the eigenvalue of LM is typically not of interest for SSSA. The inverse power method addresses this issue by using (A − σI r ) −1 instead of A in (12), where the spectral shift σ represents an initial guess of an eigenvalue λ i . The k-th step of the inverse power method is as follows: Provided that σ is a good initial guess of λ i , b (k) converges to the corresponding right eigenvector u i . Hence, the inverse power method is able to calculate the eigenvector that corresponds to any single eigenvalue, by computing the dominant eigenvector of (A − σI r ) −1 . Then, the eigenvalue λ i is approximated by the Rayleigh quotient that is given by (13).
A bad selection of the spectral shift σ leads to very slow convergence, or convergence to a different eigenvector than the one desired, despite the ability of the inverse power method to compute any eigenpair (λ i , u i ). The Rayleigh quotient iteration method is an extension of the inverse power method that updates the applied spectral shift at each step of the iteration. Starting from an initial eigenpair guess (σ (0) , b (0) ), the k-th step, k = 1, 2, 3, . . ., of the Rayleigh quotient iteration is: The method is able to converge to any eigenpair (λ i , u i ), depending on the initial guess.
For the sake of simplicity, numerical methods were presented in (12), (14), (15) for the solution of the LEP; however, they can be modified to handle also the GEP. For example, considering the GEP (7), the inverse power method for finding an eigenvector of the pencil sE − A can be described as: The methods that are described above calculate a single eigenvector in one iteration and, for this reason, they are referred in the literature as single vector iteration methods. Calculating multiple eigenpairs with a single vector iteration method is possible by applying a deflation technique and repeating the iteration.

Simultaneous Vector Iteration Methods
Simultaneous vector iteration or subspace method is a vector iteration method that calculates multiple eigenvectors simultaneously in one iteration. Suppose that we want to compute the p LM eigenvalues of the LEP (8) with pencil sI r − A. Extending the idea of the power method, consider the subsequence (11), where now b is not a vector, but b ∈ C r×p . As it is, the subsequence does not converge, as we would desire, to a matrix with columns the p eigenvectors of sI r − A. However, convergence to these eigenvectors can be achieved by ensuring that the columns of the k-th element A k b are orthonormal vectors. Subsequently, we say that the columns of A k b span the subspace: The steps of the subspace iteration algorithm are the following.
1. The iteration starts with an initial matrix b (0) , b (0) ∈ C r×p , with orthonormal columns, i.e., (b (0) ) H b = I p . 2. At the k-th step of the iteration, k = 1, 2, 3, . . .: (a) Matrix C, C ∈ C r×p , is formed as: is updated by maintaining column-orthonormality, through the QR decomposition of matrix C: where Q is a unitary matrix and R is an upper triangular matrix.
The columns of b (k) eventually converge to p right eigenvectors that correspond to the p LM eigenvalues.
Note that the LM eigenvalues may not be the important ones for the needs of SSSA. Similarly to the discussion of the inverse power method, finding the eigenvalues of smallest magnitude is possible by forming an inverse subspace iteration.
The convergence of the subspace iteration is, in general, slow and, thus, its use is avoided for complex problems. Still, acceleration of the speed of convergence is possible by combining the method with the Rayleigh-Ritz procedure. Because the utility of the Rayleigh-Ritz procedure goes far beyond the subspace iteration, we describe it separately in Appendix A.

Schur Decomposition Methods
Schur decomposition-based methods take advantage of the fact that, for any r × r matrix A, there always exists a unitary matrix Q ∈ C r×r , such that: where T, T ∈ C r×r , is an upper triangular matrix with diagonal elements all the eigenvalues of the pencil sI r − A. Decomposition (20) is known as the Schur decomposition of A, and matrix T as the Schur form of A. It appears that the most efficient algorithm for computing the Schur decomposition of a given matrix is the QR algorithm. In an analogy to the previously described methods, the QR algorithm can be understood as a nested subspace iteration or as a nested sequence of r power iterations [50]. Starting with matrix A (0) = A, the main steps of the QR algorithm at k-th iteration, k = 1, 2, 3, . . . are: where Q (k) is orthogonal and R (k) is upper triangular. The QR decomposition (21) is computed by applying Householder transformations, see [51]. If A (k) converges, say, after κ steps, then A (κ) = T. The diagonal elements of A (κ) are the eigenvalues of A. The corresponding right eigenvectors are found as the columns of the matrix: where Q ∈ C r×r . Notice that, throughout the algorithm, matrices A (k) , A (k−1) are similar and, thus, they have the same eigenvalues with the same multiplicities. The QR algorithm uses a complete basis of vectors and computes all of the eigenvalues of the matrix pencil sI r − A. In practice, the QR algorithm is not typically used as it is, but multiple shifts are applied to accelerate convergence, which leads to the implicitly shifted QR algorithm [52]. The analog of the QR algorithm for the GEP is the QZ algorithm. The QZ algorithm takes advantage of the fact that for any r × r matrices A, E, there always exist unitary matrices Q ∈ C r×r , Z ∈ C r×r , such that: are upper triangular, and the eigenvalues of the pencil sE − A are the ratios w 1,ii /w 2,ii of the diagonal elements of T 1 , T 2 , respectively. The decomposition (24) is known as generalized Schur decomposition. From the implementation viewpoint, the algorithm computes all of the eigenvalues of the matrix pencil sE − A by first reducing A, E, to upper Hessenberg and upper triangular form, respectively, through Householder transformations. The method can be perceived as equivalent to applying the implicitly shifted QR algorithm to AE −1 , without, nonetheless, explicitly computing E −1 [52].

Krylov Subspace Methods
The power method consists in the calculation of the products Ab, A 2 b, A 3 b, and so on, until convergence. The method utilizes the converged vector, let us say A κ b, which corresponds to the LM eigenvalue. However, the power iteration discards the information of all vectors calculated in the meantime. The main idea of Krylov subspace eigenvalue methods is to utilize this intermediate information in order to calculate the p LM eigenvalues of sI r − A. The Krylov subspace that corresponds to A and vector b is defined as: Subspace (25) is spanned by the columns of the Krylov matrix, which is defined as: The matrix-vector multiplications typically render the columns of K p linearly dependent. Thus, these columns require orthonormalization. The orthonormalized vectors can be then employed in order to provide approximations of the eigenvectors that correspond to the p eigenvalues of LM. The eigenvalues are extracted by projecting A onto the Krylov subspace, typically through the Rayleigh-Ritz procedure (see the Appendix A).

Arnoldi Iteration
The Arnoldi iteration finds the p LM eigenvalues of sI r − A. In particular, it employs the Gram-Schmidt process in order to find an orthonormal basis of the Krylov subspace. Subsequently, the eigenvalues are extracted by applying the Rayleigh-Ritz procedure. Starting from an initial normalized vector q 1 , the vector q k+1 at the k-th step is calculated, as follows: where h i,k = q H i Aq k , h k+1,k = ||q k ||. At the k-th step, the Arnoldi relation is formulated, as follows: where Q k , Q k ∈ C r×k , is the matrix with columns q k ; H k , H k ∈ R k×k , is in Hessenberg form and it has elements h i,j ; e k denotes the k-th column of the identity matrix I k . The algorithm stops when h k+1,k = 0, suppose when k = p. Then, relation (28) becomes: where Q p = [q 1 , q 2 , . . . , q p ] is an orthonormal basis of the Krylov subspace. Following the Rayleigh-Ritz procedure, the eigenvalue The Arnoldi iteration converges quickly if the initial vector q 1 has larger entries at the direction of desired eigenvalues. This is usually not the case and, thus, it is common that many iterations are required. On the other hand, the ability to compute the columns of Q p is constrained, because of its high memory requirements. For this reason, in practice, the Arnoldi iteration is typically restarted after a number of steps, while using a new, improved initial vector. The restarts can be explicit or implicit. The idea of explicit restarts is to utilize the information of the most recent factorization in order to produce a better initial vector. This is usually combined with a locking technique, i.e., a technique to ensure that converged vectors do not change in successive runs of the algorithm. On the other hand, the idea of the Implicitly Restarted (IR) Arnoldi iteration is to combine the Arnoldi iteration with the implicitly shifted QR algorithm.

Krylov-Schur Method
A more recently proposed idea that achieves the effect of the implicit restart technique in a simpler way is the Krylov-Schur method, which is a generalization of the Lanczos thick restart [53] for non-Hermitian matrices. The method considers the Arnoldi relation (28) and applies the QR algorithm to bring matrix H k to its Schur form, i.e., T k = W H k H k W k , where T k is an upper triangular matrix with diagonal elements the eigenvalues of H k (Ritz values). Subsequently, (28) becomes: or equivalently, where (31) is reordered in order to separate undesired Ritz values from the desired ones. The part that corresponds to the desired Ritz values is kept and the rest is discarded: where D k,d , T k,d , and w H k,d represent the parts of D k , T k and w H k that correspond to the desired Ritz values after the reordering, respectively. Subsequently, the reduced decomposition (32) is expanded to order k. The above process is repeated until convergence to an invariant subspace is achieved.

Generalized Eigenvalue Problem
Krylov subspace methods above were described for the solution of the LEP. Using the Krylov subspace methods for the solution of the GEP is usually done by employing a spectral transform. For example, applying the shift & invert transform to the pencil sE − A, the Arnoldi relation becomes: where, in this case, it is Q H k EQ k = I k . Each eigenvalue µ obtained after the convergence of (33) to an invariant subspace, is then transformed in order to find the corresponding eigenvalue λ of the original problem as λ = 1/µ + σ.

Contour Integration Methods
Contour integration methods find the eigenvalues of the pencil sE − A that are inside a given, user-defined domain of the complex plane. The method that was proposed by Sakurai and Sugiura [33] and its variants [34,54], are the most characteristic examples of this class.
Suppose that p distinct eigenvalues λ 1 , λ 2 , . . . , λ p are located inside positively oriented simple closed curve Γ in C. The contour integration method with Hankel matrices (CI-Hankel) applies a moment-based approach in order to reduce the problem of finding the p eigenvalues of sE − A to finding the eigenvalues of a p × p matrix pencil. For a non-zero vector b, consider the moments: where k = 0, 1, . . . , p − 1; γ is a user-defined point that lies in Γ. A standard case is to define Γ as a circle with center γ and radius ρ.
In the case that some eigenvalues in the defined region are very close, the accuracy of the CI-Hankel method decreases, as the associated Hankel matrices become ill-conditioned. An alternative approach is to use contour integration to construct a subspace that is associated to the eigenvalues in Γ, and then extract the eigenpairs from such subspace while using the Rayleigh-Ritz procedure. The resulting method is the CI-RR (Contour Integration with Rayleigh-Ritz).
Based on (35), it can be proven that, if the column vectors {q 1 , q 2 , . . . , q p } form an orthonormal basis of span{ψ o , ψ 1 , . . . , ψ p−1 }, then the eigenvalues λ i and the corresponding eigenvectors u i , i = 1, 2, . . . , p can be extracted while using the basis Q p = [q 1 , q 2 , . . . , q p ] and applying the Rayleigh-Ritz procedure for the pencil sE − A. In order to construct the orthonormal basis Q p , the contour integral (35) is first approximated through the N-point trapezoidal rule. In practice, an orthonormal basis Q p with size larger than p may be used to improve accuracy. When compared to the Hankel method, the CI-RR method is typically more accurate and, thus, is preferred most of the time.
Finally, another promising contour integration algorithm is FEAST. FEAST was first proposed in [35] for the solution of symmetric eigenvalue problems, but it was later extended to include non-symmetric eigenvalue problems. The algorithm can be understood as an accelerated subspace iteration method combined with the Rayleigh-Ritz procedure [55] and, from this point of view, it has some similarities with the CI-RR method. The unique characteristic of the FEAST algorithm is that it implements an acceleration through a rational matrix function that approximates the spectral projector onto the subspace.

Open-Source Libraries
This section provides an overview of-to the best of our knowledge-all of the open-source numerical eigensolvers that implement state-of-art numerical algorithms for non-symmetric eigenvalue problems. These are LAPACK, ARPACK, Anasazi, SLEPc, FEAST, and z-PARES.  [59]. The first version of LAPACK was released in 1992. When compared to EISPACK, LAPACK was restructured to include the block forms of QR and QZ algorithms, which allows for exploiting Level 3 BLAS and leads to improved efficiency [60]. The latest version of LAPACK is 3.7 and it was released in 2016.

ARPACK
• Summary: ARPACK [40] is a library developed for solving large eigenvalue problems with the IR-Arnoldi method. • Eigenvalue Methods: IR-Arnoldi iteration. • Matrix Formats: ARPACK supports the Reverse Communication Interface (RCI), which provides to the user the freedom to customize the matrix data format as desired. In particular, with RCI, whenever a matrix operation has to take place, control is returned to the calling program with an indication of the task required and the user can, in principle, choose the solver for the specific task independently from the library.

•
Returned Eigenvectors: only right eigenvectors are calculated. • Dependencies: ARPACK depends on a number of subroutines from LAPACK/BLAS. Moreover, ARPACK requires to be linked to a library that factorizes matrices. This can be either dense or sparse. In the simulations that are described in this paper, we linked ARPACK to the efficient library KLU, which is part of SuiteSparse [61], and that is particularly suited for sparse matrices whose structure originates from an electrical circuit. • GPU-based version: to the best of the authors' knowledge, a library that provides a functional GPU-based implementation of ARPACK is not available to date.

FEAST
• Summary: FEAST [43] is the eigensolver that implements the FEAST algorithm, which was first proposed in [35]. Among other characteristics, the package includes the option to switch to IFEAST, which uses an inexact iterative solver to avoid direct matrix factorizations. This feature is particularly useful if the sparse matrices are very large and carrying out direct factorization is very expensive.

Summary of Library Features
Tables 3 and 4 provide a synoptic summary of the methods and relevant features of open-source libraries that solve non-symmetric eigenvalue problems. All of the libraries can handle both real and complex arithmetic types.
As it can be seen from Table 4, not all libraries provide algorithms that allow for calculating both left and right eigenvectors at once. However, in order to calculate participation factors, which are an important tool of power system SSSA [66][67][68], both right and left eigenvectors are required. With this regard, a formula to directly calculate left eigenvectors from right ones was proposed in [69]. In general, which is the most efficient solution for the calculation of the left eigenvalues depends on the library and the eigenvalue problem.

Case Studies
This section presents simulation results based on two real-world size power system models. The first system is a detailed model of the All-Island Irish Transmission System (AIITS), which includes 1443 state variables and 7197 algebraic variables. The second system is a dynamic model of the European Network of Transmission System Operators for Electricity (ENTSO-E) system, which includes 49,396 state variables and 96,770 algebraic variables. Table 5 summarizes the versions and dependencies of the open-source libraries considered in this section. Note that we only consider the open-source libraries that we managed to compile and install on Linux and Mac OS X operating systems and that worked for relatively "large" eigenvalue problems. All of the simulations are obtained using the Python-based software tool Dome [70]. The Dome version utilized for this paper is based on Fedora Linux 28, Python 3.6.8, CVXOPT 1.1.9, and KLU 1.3.9. Regarding the computing times reported in both examples, we have two comments. First, all of the simulations were executed on a server mounting two quad-core Intel Xeon 3.50 GHz CPUs, 1 GB NVidia Quadro 2000 GPU, 12 GB of RAM, and running a 64-bit Linux OS. Second, since not all method implementations include two-sided versions and, in order to provide as a fair comparison as possible, all of the eigensolvers are called so as to return only the calculated eigenvalues and not eigenvectors.

All-Island Irish Transmission System
In this case study, we consider a real-world model of the AIITS. The topology and the steady-state operation data of the system have been provided by the Irish transmission system operator, EirGrid Group, whereas the dynamic data have been defined based on our knowledge about the technology of the generators and the controllers. The dynamic model has been also validated while using frequency data from a severe event that occurred in the real system in 2018, see [71]. The system consists of 1479 buses, 796 lines, 1055 transformers, 245 loads, 22 synchronous machines, with Automatic Voltage Regulators (AVRs) and Turbine Governors (TGs), 6 Power System Stabilizers (PSSs), and 176 wind generators. In total, the dynamic model has n = 1443 state variables and m = 7197 algebraic variables. The results of the eigenvalue analysis of the AIITS are discussed for both LEP and GEP and for a variety of different numerical methods, namely, QR and QZ algorithms by LAPACK, GPU-based QR algorithm by MAGMA, subspace iteration, ERD-Arnoldi, and Krylov-Schur methods by SLEPc, IR-Arnoldi by ARPACK; and, CI-RR by z-PARES. In particular, we provide, for each method, the rightmost eigenvalues, as well as the time that is required to complete the solution of the eigenvalue problem. The dimensions of the LEP and GEP for the AIITS are shown in Table 6. Table 7 presents the results obtained with Schur decomposition methods. Both QR and QZ algorithms find all 1443 finite eigenvalues of the system. For the GEP, the QZ algorithm also finds the additional infinite eigenvalue with its algebraic multiplicity. The obtained rightmost eigenvalues are the same for both LEP and GEP. Since LAPACK is the most mature software tool among those considered in this paper, the accuracy of the eigenvalues found with all other libraries is evaluated by comparing them with the reference solution that was computed with LAPACK. Figure 1 shows the system root loci plot.

Dense Algorithms
Regarding the computational time, we see that, for the LEP, both LAPACK and the GPU-based MAGMA are very efficient at this scale, with MAGMA only providing a marginal speedup. On the other hand, when it comes to solving the GEP with LAPACK's QZ method, scalability becomes a serious issue, since the problem is solved in 3669.77 s, which is computationally cumbersome.   We show the image of the spectrum of the AIITS system for a couple of common special Möbius transforms, in particular for the shift & invert and the Cayley transform. Figures 2 and 3 present the results, in whichλ denotes an eigenvalue of the transformed pencil. These results refer to the LEP, and they are obtained using LAPACK. In each figure, the stable region is shaded, while the stability boundary is indicated with a solid line. The 5% damping boundary is indicated with a dash-dotted line.
For the shift & invert transform, the stability boundary is defined by the circle with center γ = (1/2σ, 0) and radius ρ = 1/2σ. If σ < 0, wgucg is the case of Figure 2, stable eigenvalues are mapped outside the circle. On the other hand, if σ > 0, stable eigenvalues are mapped inside the circle. If σ = 0, we obtain the dual pencil with the corresponding invert transform, and the stable region is the full negative right have plane. Finally, Figure 3 shows the image of the Cayley transform of the AIITS for σ = 1.2. All of the stable eigenvalues are located inside the unit circle with center the origin.

Sparse Algorithms
The implementation of the subspace iteration by SLEPc only finds the desired number of LM eigenvalues. However, in the s-domain, the relevant eigenvalues from the stability point of view are not the LM ones, but the ones with LRP or SM. Especially for the GEP, the LM eigenvalue is infinite and, hence, does not provide any meaningful information on the system dynamics. For this reason and for the needs of power system SSSA, the subspace method and, in general, any method that looks for LM eigenvalues must always be combined with a spectral transform. For the needs of this example, we apply the invert transform and pass to SLEPc the pencil of the dual system, i.e., zA − E. Subsequently, the method looks for the 50 LM eigenvalues of the dual system, which correspond to the 50 SM eigenvalues of the prime system. With this setup, the eigenvalues found by the subspace iteration for the GEP are shown in Table 8. As can be seen, the pair −0.1322 ± 0.4353 is not captured, since its magnitude is larger than the magnitudes of the 50 SM eigenvalues. In order to also obtain this pair, one can customize the spectral transform or simply increase the number of the eigenvalues to be returned. However, the best setup is not known a priori and, thus, some heuristic parameter tuning is required. Finally, the method does not scale well, since the solution of the GEP is completed in 6807.24 s. The rightmost eigenvalues found with Krylov subspace methods for the LEP and GEP are shown in Tables 9 and 10, respectively. For the LEP, ARPACK is set up in order to find the 50 LRP eigenvalues. Although all of the eigenvalues shown in Table 9 for ARPACK are actual eigenvalues of the system, some of the LRP ones are missed. Furthermore, no correct eigenvalues were found for the GEP, which we attribute to the fact that a non-symmetric E is not supported. In SLEPc methods, both for LEP and GEP and in order to obtain the eigenvalues with good accuracy, we use the option "Target Real Part" (TRP), which allows targeting eigenvalues with specified real part. In particular, we set the TRP parameter to −0.01, and apply a shift & invert transform with σ = −0.01. ERD-Arnoldi and Krylov-Schur methods are both able to accurately capture all rightmost eigenvalues. For the sake of completeness, we mention that the eigenvalues obtained with SLEPc, when compared to the ones that were found by LAPACK, appeared to be shifted by a constant offset −σ, i.e., 0.01 was returned instead of 0, and so on. The results shown in Tables 9 and 10 take into account such a shift by adding σ to all output values that were returned by SLEPc. Finally, the Krylov subspace methods by SLEPc appear to be more efficient than ARPACK's IR-Arnoldi. When compared to Schur decomposition methods, at this scale, Krylov methods, although they require some tuning, appear to be by far more efficient for the GEP, but less efficient for the LEP. Table 11 presents the results that are produced by z-PARES' CI-RR method for the LEP and GEP. The method is set to look for solutions in the circle with center the point γ = (−0.01, 4) and radius ρ = 8. In both cases, the eigenvalues found by z-PARES are actual eigenvalues of the system, although the eigenvalues found for the GEP include noticeable errors, when compared to the results that were obtained with LAPACK.
The most relevant issue is that the eigenvalues obtained with z-PARES are not the most important ones for the stability of the system, which means that critical eigenvalues are missed. This issue occurs despite the defined search contour being reasonable. Of course, there may be some region for which the critical eigenvalues are captured, but this can not be known a priori. Regarding the simulation time, the method for the AIITS is faster than SLEPc's Krylov subspace methods for the LEP, but slower for the GEP. Figure 4 shows the search contour and the location of the characteristic roots that are found by z-PARES for the LEP.

European Network of Transmission System Operators for Electricity
This example presents the simulation results for a dynamic model of the ENTSO-E. The system includes 21,177 buses (1212 off-line); 30,968 transmission lines and transformers (2352 off-line); 1144 zero-impedance connections (420 off-line); 4828 power plants represented by 6-th order and 2-nd order synchronous machine models; and, 15, 756 loads (364 off-line), modeled as constant active and reactive power consumption. Synchronous machines that are represented by 6-th order models are also equipped with dynamic AVR and TG models. The system also includes 364 PSSs.
The system has, in total, n = 49, 396 state variables and m = 96,770 algebraic variables, as summarized in Table 12. The pencil sE − A has dimensions 146,166 × 146,166 and the matrix A has 654,950 non-zero elements, which represent the 0.003% of the total number of elements of the matrix. Neither the LEP nor GEP could be solved while using Schur decomposition methods. At this scale, the dense matrix representation that is required by LAPACK and MAGMA libraries leads to massive memory requirements, and a segmentation fault error is returned by the CPU. Among the algorithms that support sparse matrices, here we only test the contour-integration based methods, which, in fact, were the ones that were able to tackle this large eigenvalue problem on the available hardware.
The effect of changing the search region of the z-PARES CI-RR method on the eigenvalue analysis of the ENTSO-E system is shown in Table 13. Interestingly, simulations showed that shrinking the defined contour may lead to a marginal increase of the computation time. Although not intuitive, this result indicates that the large size of the ENTSO-E system mainly determines the mass of the computational burden, and that, at this scale, smaller subspaces are not necessarily constructed faster by the CI-RR algorithm. Regarding the number of eigenvalues obtained, using a region that is too small leads, as expected, to missing an important number of critical eigenvalues. We test the impact of applying spectral transforms to the matrix pencil sE − A, on the eigenvalue analysis of the ENTSO-E with z-PARES' CI-RR. In particular, we test the invert transform that yields the dual pencil zA − E; and, the inverted Cayley transform, i.e., s = (z + 1)/(σz − σ), which yields the pencil z(E − σA) − (−σA − E). Table 14 shows the results. Passing the transformed matrices to z-PARES provides a marginal speedup to the eigenvalue computation. In addition, when considering either the prime system or the inverted Cayley transform with σ = −1, results in finding the same number of eigenvalues, whereas, when the dual system is considered, a number of eigenvalues is missed.

Remarks
The following remarks are relevant: • With regard to dense matrix methods, their main disadvantage is that they are computationally expensive. In addition, they generate complete fill-in in general sparse matrices and, therefore, cannot be applied to large sparse matrices simply because of massive memory requirements. Even so, LAPACK is the most mature among all computer-based eigensolvers and, as opposed to basically all sparse solvers, requires practically no parameter tuning. For small to medium size problems, the QR algorithm with LAPACK remains the standard and most reliable algorithm for finding the full spectrum for the conventional LEP.

•
As for sparse matrix methods, convergence of vector iteration methods can be very slow and, thus, in practice, if not completely avoided, these algorithms should only be used for the solution of simple eigenvalue problems. An application where vector iteration methods may be more relevant, is in correcting eigenvalues and eigenvectors that have been computed with low accuracy, i.e., as an alternative to Newton's method [69]. With regard to Krylov subspace methods, the main shortcoming of ARPACK's implementation is the lack of support for general, non-symmetric left-hand side coefficient matrices, which is the form that commonly appears when dealing with the GEP of large power system models. On the other hand, the implementations of ERD-Arnoldi and Krylov-Schur by SLEPc do not have this limitation and exploit parallelism while providing good accuracy, although some parameter tuning effort is required. In addition, for the scale of the AIITS system and the GEP, these methods appear to be by far more efficient than LAPACK. Moreover, the implementation of contour integration by z-PARES is very efficient and can handle systems at the scale of the ENTSO-E. An interesting result is that, at this scale, smaller contour search regions, even if properly setup, do not necessarily imply faster convergence or, equivalently, smaller subspaces are not necessarily constructed faster by the CI-RR algorithm. The most relevant issue for z-PARES is that, depending on the problem, it may miss some critical eigenvalues, despite the defined search contour being reasonable. Although there may be some parameter settings for which this problem does not occur, those can not be known a priori. • Finally, the presented comparison of numerical methods and libraries is not specific of power system problems, but it is relevant to any system with unsymmetric matrices. A relevant problem of the eigenvalue analysis for any dynamical system is to utilize its physical background in order to obtain an efficient approximation to some (components of) the eigenvectors. A computer algorithm that has successfully addressed this problem for conventional, synchronous machine dominated power systems is AESOPS [23,72], a heuristic quasi Newton-Raphson based method that aimed at finding dominant electromechanical oscillatory modes. On the other hand, how to effectively exploit the physical structure of modern, high-granularity power systems to shape efficient numerical algorithms is a challenging problem and an open research question.

Conclusions
The paper provides a comprehensive comparison of algorithms and open-source software tools for the numerical solution of the non-Hermitian eigenvalue problems that arise in power systems.
In particular, the paper considers state-of-art implementations of methods that are based on Schur decomposition, vector iteration, Krylov subspace, and contour integration.
The most robust and reliable method for the solution of the LEP for small to medium size problems is the QR algorithm with LAPACK. For large-scale problems, which is the case of real-world power systems, the only option is to use software tools that exploit matrix sparsity. While different options exist, such methods are less mature than the QR algorithm. The simulation results discussed in this paper recommend that, among the available options, the most promising for power systems are the parallel versions of the Krylov-Schur method by SLEPc and CI-RR method by z-PARES. In particular, SLEPc's Krylov-Schur provides the best accuracy, while z-PARES's CI-RR was the only one that was able to tackle the 146,166 × 146,166 GEP of the European Network of Transmission System Operators for Electricity (ENTSO-E).