A Critical Eigenvalues Tracing Method for the Small Signal Stability Analysis of Power Systems

The continuation power flow method combined with the Jacobi-Davidson method is presented to trace the critical eigenvalues for power system small signal stability analysis. The continuation power flow based on a predictorcorrector technique is applied to evaluate a continuum of steady state power flow solutions as system parameters change; meanwhile, the critical eigenvalues are found by the Jacobi-Davidson method, and thereby the trajectories of the critical eigenvalues, Hopf bifurcation and saddle node bifurcation points can also be found by the proposed method. The numerical simulations are studied in the IEEE 30-bus test system.


Introduction
The stability analyses of the rapid increase in complexity of modern power systems require the development of efficient and robust numerical methods.In terms of small signal stability analysis, power systems are usually modeled by systems of differential-algebraic equations with large nonlinear dynamics.Simulation in frequency domain of the resulting dynamical systems heavily relies on numerical methods for eigenvalue problems and systems of linear equations.Direct application of conventional methods (e.g., QR method) for power system eigenvalue problems is computationally not feasible or inefficient.Thus, many special intensive researches for the small signal stability analysis of power systems have been proposed in the last two decades, which devoted to evaluating a selected (critical) subset of eigenvalues associated with the complete system response.
In [1,2], only the electromechanical modes; i.e. the rotor angle ones, are considered to be the critical eigenvalues and computed by a frequency response approach without formulating the system state matrix.A more general modal model reduction method, Dominant Pole Algorithm (DPA), has been proposed in [3] to compute the dominant poles in any specified Single Input Single Output (SISO) transfer function.The selective modal analysis approach [4] uses a reduced order model to compute the desired critical eigenvalues relevant to the selected modes.The invariant subspace iteration methods [5][6][7][8][9][10][11][12][13][14][15] are based on the use of the augmented system state equations, exploiting the Jacobian matrix sparsity, and most of them incorporate spectral transformations such as the shift-invert transformation [5,6], the S-matrix method [7] and the Cayley transformation [8,9] for the identification of the set of dominant eigenvalues (i.e., the largest modulus), which are the mapping of the critical (i.e., rightmost) eigenvalues of the system state matrix.Among them, the Jacobi-Davidson (JD) method [14,15] is a recently developed subspace iteration method.Based on its characteristic on iterative construction of a partial Schur form, the deflation technique, and the effective restart, the JD method is suitable for the efficient computation of a selected partial eigenvalues and is highly effective in detecting the clustered eigenvalues.
All of the above mentioned critical eigenvalue solvers are normally applied to calculate critical eigenvalues under a given system operating scenario.As the system equilibrium point often vary with respect to any control and load variations, critical eigenvalues should also be recalculated.In such case, the critical eigenvalue trace methods are of interest to the researchers.In [16], an integration based eigenvalue trajectory tracing method was proposed to identify both oscillatory stability margin and damping margin, and indexes were derived to identify Hopf bifurcation.A critical eigenvalue tracing method based on the continuation of the invariant subspace algorithm combined with the projected Arnoldi method is proposed in [17].
In this paper, the continuation power flow combined with the Jacobi-Davidson method is proposed to trace the critical eigenvalues for power system small signal stability analysis.

Small Signal Stability
The dynamic behavior of power systems can be described by a set of nonlinear differential equations together with a set of algebraic equations (DAEs) [18] ( , , ) 0 ( , , ) where x is the vector of the state variables, y is the vector of the algebraic variables, and u is the vector of system parameters, e.g., the load level of the whole system.In the case of small signal stability analysis, (1) is linearized around a system operating point to yield the linearized DAEs model where A total is the total system Jacobian matrix,  denotes an incremental change in steady state value.By eliminating the algebraic variables y in (2), the equivalent system state matrix can be formulated as 1 ( ) ( ) ( ) ( ) According to the Lyapunov's first method [19], the small signal stability of the nonlinear system (1) in the neighborhood of the operating point can be analyzed by inspecting the eigenvalues of the linearized system (2); e.g., all of the eigenvalues of the system state matrix A should have negative real parts to maintain the stability.Accordingly, eigenvalues close to the imaginary axis of the complex plane are the critical eigenvalues of interest to analyze the stability of power system oscillations.Equation ( 3) describes that the critical eigenvalues and corresponding eigenvectors construct a dominant invariant subspace with respect to the entire eigenspace of the system state matrix A.
where the column vectors in V(u) are the basis in the invariant subspace, and the eigenvalues of (u) are a fraction of the corresponding eigenvalues in the full matrix A(u).Additionally, when the operating point varies with respect to the system parameters changes, the critical eigenvalues and corresponding eigenvectors may also change.A Hopf bifurcation arises when there is a complex conjugate eigenpair cross the imaginary axis first.In such case, the system becomes oscillatory unstable.At the loadability limit, the traditional Newton-Raphson method of obtaining load flow solution will break down.In this case, the system Jacobian will become singular, thus so-called saddle node bifurcation.

Critical Eigenvalues Tracing Method
The idea of proposed method is based on combing continuation power flow with the Jacobi-Davidson method to identify and trace the critical eigenvalue trajectories.The continuation power flow calculates a series of power system steady state solutions with respect to a given power injection variation scenario.During the continuation power flow process, the Jacobi-Davidson method with deflation technique is applied to effectively calculate the critical eigenvalues, and thereby the critical eigenvalue trajectories, Hopf bifurcation (HB) and saddle node bifurcation (SNB) points can also be found.

Continuation Power Flow
The concept of the continuation power flow [20] is based on a predictor-corrector technique to evaluate a continuum of steady state power flow solutions starting at some base load and leading to the steady state stability limit (SNB) of the system.A salient feature of the continuation power flow is that it remains well-conditioned around the SNB point.As a consequence, divergence due to ill-conditioning is not encountered in the vicinity of SNB point, even when single-precision computation is used.
The continuation power flow based on a predictorcorrector technique is described as follows, and is shown in Figure 1.
 Predictor step: The function of the predictor is to find an approximation point for the next solution.-Suppose the continuation process is at the ith step with the solution ( , , ).
x y u The predictor tries to find an approximation point( , , ) x y u for next solution )  , and can be formulated as Copyright © 2013 SciRes.
where u is the loading parameter which will vary from base case to the point of maximum loadability,  is the step-size for the next prediction, and where e k is a column vector of zeros with a single 1 at the position of the unknown that is chosen to be the continuation parameter.
 Corrector step: Using the predicted value as the initial condition ) , , for the nonlinear iteration, the augmented power flow equations ( 5)-( 6) are then solved by the Newton iterative method to achieve the solution ) , , (

Jacobi-Davidson Method
The Jacobi-Davidson (JD) subspace iteration method [14] is applicable to a selected eigenvalue and corresponding eigenvector approximations of general unsymmetric matrix A, in which each iteration combines the idea of Davidson's and Jacobi's methods.In this approach, a search subspace span{V} is generated onto which the given eigenproblem, , is projected, where V is a complex n matrix and its columns constitute an orthonormal basis v 1 ,v 2 ,…,v j , j<<n.The 'Davidson' part, based on the Ritz-Galerkin condition: , is to select an approximate eigenpair of A with the reduced system, which can be represented as Then the projected eigenproblem ( 7) is solved and a solution is selected.The Ritz pair is defined as to form an approximate eigenpair of A, with the residual vector . The 'Jacobi' part forms an orthogonal correction for the current eigenvector approximation q ~, by the solution of q v  in the following correction equation Equation ( 8) defines an optimal expansion (orthogonal basis) of the search subspace, span{V, v}; the search subspace accumulates valuable information for the desired eigenvector approximation for every iteration.
When the search subspace V reaches a certain maximum dimension n j  max , the JD method adopts an implicit restart strategy to reduce the dimension of the search subspace from j max to j min (j min ≤ j ≤ j max ) by discarding the columns , and continues the next JD iteration with min .Here, the JD algorithm is obviously easy to implement because is already orthogonal, and the JD algorithm is repeated until the norm of the residual ||r|| is smaller than a given tolerance.A JD-style method, called JDQR [15], is designed to find a number of desired eigenpairs by constructing a partial Schur form of A iteratively.A deflation technique has been successfully incorporated into the JDQR method, which would avoid repeated computation of the detected eigenvalues.The JDQR method is described as follows.
Suppose that k-1 Schur pairs have been detected and formed with the partial Schur form , where A suitable new Schur pair will be used to expand and , so that where , which satisfies Thus, the new Schur pair is also an eigenpair of the deflated matrix , and the deflated matrix has the same eigenvalues as A except the detected eigenvalues that are transformed to zero; the deflation part effectively avoids repeated computation of the detected eigenvalues.Then, the search subspace span {V} is expanded by the orthogonal complement of v to V, where v is the solution of the following deflated correction equation (11) Copyright © 2013 SciRes.EPE where and .

Implementation of Algorithm
The flowchart of the solution algorithm is shown in Figure 2, and the step-by-step procedure is described as follows: a) Solve base case power flow solution: Perform a specified base case Newton-Raphson power flow for evaluating the initial conditions of state variables.
b) Find Jacobian matrix: Compute the linearized DAEs model ( 2) by both the initial conditions and the power flow solution.
c) Find critical eigenvalues: Compute the critical eigenvalues of the Jacobian matrix by the JD method with a given sorting criterion, i.e. the rightmost eigenvalues or the least damping ratio eigenvalues.d) Identify Hopf bifurcation: Detect the real part of the critical eigenvalues for Hop bifurcation, i.e. the complex conjugate eigenvalue crosses through the imaginary axis of the complex plane first.
e) Find the next equilibrium Perform the continuation power flow with a given control strategy (a specified control variable and increment factor) for the next operating equilibrium states.

Numerical Results
In this section, numerical examples were examined through the IEEE-30 bus test system.The system has 6 machines and 59 state variables, and the dynamic model is linearized around a base operation point with a total load of 275.2 MW. Figure 3 shows the distribution of the partial rightmost eigenvalues at the based load.The software package MATLAB v7.1 is used to implement the proposed method.In our work, the damping ratio and the real part of the eigenvalues were used as the selection criteria to trace the critical eigenvalues.The convergence tolerance of the JD method is set to 10 -8 .The load increase pattern is mainly performed on the PQ buses starting at the base load and leading to the steady state stability limit of the system.

Tracing of the Rightmost Eigenvalues
For this case, the critical eigenvalues with largest real part are the desired eigenvalues.There are 6 rightmost eigenvalues are traced by the proposed method as shown in Figure 4.The Hopf bifurcation occurs on the 1.859 p.u. load level with respect to the base case, i.e., there is an eigenpair passing through the imaginary axis of the complex plane.When the system has a load level of 2.275 p.u., the system reaches to stability limit, i.e., the SNB occur.In such case, the system operates at point where the Jacobian matrix has a zero eigenvalue.When system power demands go beyond the maximum power limit, the power system is divergent and becomes unsolvable.The continuation power flow remains wellconditioned around the SNB point.

Tracing of the Eigenvalues with Least Damping Ratio
In this case, the load increase pattern is the same as the case in section 4.1 and therefore has predicable same results.The critical eigenvalues with least damping ratios are the desired eigenvalues.There are 6 eigenvalues with least damping ratios are traced by the proposed method as shown in Figure 5.The Hopf bifurcation occurs on the 1.859 p.u. load level with respect to the base case, i.e., there is an eigenpair whose damping ratio is zero.When the system has a load level of 2.275 p.u., the system reaches to stability limit, which is not shown in Figure 5.
The above mentioned simulation results have shown that the critical eigenvalues can be effectively detected by the proposed method.

Conclusions
This paper has presented an effective approach to trace the critical eigenvalue trajectories for the small signal  stability analysis of power systems.Both the rightmost and least damping ratio eigenvalues tracing were effecttively carried out by the proposed method combining the continuation power flow with the JD method.The simulation results have shown that the Hopf bifurcation and saddle node bifurcation points can also be detected by the proposed algorithm and provide useful information for the small signal stability analysis.
lar matrix.The diagonal elements of R k-1 are eigenvalues of A, and the first column of Q k-1 is an eigenvector of A.

Figure 3 .
Figure 3. Distribution of the partial rightmost eigenvalues for the IEEE-30 bus test system.
part of c ritic al eigenv alue