Designing a Robust Power System Stabilizer via Probabilistic Analysis of System Eigenvalues

Power system stability has always been a serious issue for system operators. Power System Stabilizer (PSS) is a common equipment to improve stability conditions in generators. PSS improve system stability by shifting eigenvalues of system. A common method in traditional PSS is to presume grid parameters, as well as load characteristics as constant values. The fact is, however, in real operating conditions, system and load parameters are not constant. Therefore, different variables and uncertainties such as load value, network design and changes in controller parameters have to be considered. Proper statistic distribution is a common method for modeling uncertainties of parameters. This study focuses on designing PSS via analysis of probabilistic eigenvalues. Then, to prove design authenticity, the method is applied on a test 8-machie system.


INTRODUCTION
The common eigenvalue analysis is based on stable operation and constant parameters of system (Kundur, 1994;Kundur et al., 1989;Lin et al., 1996;Mori et al., 1993;Tse and Tso, 1998;Wang and Semlyen, 1990;Yu, 1983).However, operating the power system is based on continuous load dispatch and random variable factors.This implicates that more uncertainty sources and more variables have to be taken into account in system calculation and system planning.This study is based on analysis of probabilistic eigenvalues and it is assumed that variables have a unique distribution.
Probabilistic theory was first introduced in 1978 as a tool to study dynamic stability of power system (Burchett and Heydt, 1978a,b).Except widespread studies within 1979 to 1982, which led to different papers, there is not much work in this field of research.
In early probabilistic theory, which was suggested for analyzing system dynamic behavior, a probabilistic dispatch was suggested for the real part of eigenvalues.The suggested dispatch was based on probabilistic nature of random parameters of system.Random variables may have either a unique distribution or any other kind of dispatch (Burchett and Heydt, 1978a,b).Then, using normal dispatch of sensitive eigenvalues, stability probability of system is calculated.This method is widely used in different studies for calculating random load dispatch (Tse and Tso, 1993).Measurement inaccuracy, inaccuracy due to estimation of controller parameters and inaccuracy of load predictions are uncertainties which are consideredin this study (Burchett and Heydt, 1978a,b;Sauer and Heydt, 1978).In probabilistic dynamic studies, presented in (Burchett and Heydt, 1978a,b), uncertainties in parameters of generator, as well asuncertainties in controller parameters are considered.Considering the complexity in calculation of eigenvalues for parameter dispatch, it seems that rotor angle and damping factor are the only parameters which have been regarded for twomachine system presented in (Burchett and Heydt, 1978a,b).Parameters such as the effect of load uncertainty on the current of generator seem to be ignored for simplicity.
In (Loparo and Blankenship, 1979) second order statistic indices have been used to evaluate dynamic stability of system.In (Brucoli et al., 1981) "dynamics stability limits" curves were used to investigate the stability of a grid-connected machine.In 1988, a comprehensive study on different operating conditions was performed using "Gram series" (Maslennikov and Ustinow, 1997).Dynamic model of induction motors was introduced in 1960, in which intelligent methods are used to determine system equivalent matrix, Aeq, via eigenvalues and eigenvectors (Gibbard, 1988).
Distortion sources in a multi-machine system can be listed as: C Changes in system operating parameters C Uncertainty of controller parameters C Changes in structure and parameters of system and output of generator With respect to statistic inherent of distortion sources, the first source can be modeled as continuous random variables with respect to existing data and prediction inaccuracy.The second distortion can also be modeled via continuous random variables, which these variables express measurement inaccuracy and inaccuracy in estimation of controller parameters.The best way to explain the third distortion sources is to model it with discrete random variables.
This study only focuses on the effects of uncertainties caused by the first source.

PROBABILISTIC EQUATIONS FOR RANDOM VARIABLES
Numerical characteristic of random variables: In most practical issues, random variables can be modeled with mean value and variance of values.

Exponent (Average value):
If the discrete random variable, X, consists of x 1 , x 2 ,…. and the probability for each variable is defined as P(X = x) = p i , then exponent of X, shown by E[X] or , would be: X (1) If X is a continuous random variable, with f(x) as its dispatch function, then exponent of X would be: (2) Variance: Variance of a discrete random variable can be obtained as: (3) in which F is deviation.Variance represents the distribution of probabilistic values around the average value.Also, for continuous random variable, X: Covariance: Covariance for X and Y, as two discrete random variables, can be expressed as: (5 in which is the exponent of Y.

Y
Positive covariance is an index stating that variations of X and Y are in the same direction.
Usually, when X is big in value, Y would also be big and vice versa.Therefore, is usually positive.
Negative value for covariance represents dissimilarity in changes of X and Y.
Normal dispatch: If probabilistic dispatch for continuous variable, X, is defined as: Then X obeys its normal dispatch.Determining normal dispatch is based on exponent and variance.The value of depicts the function position in align with X horizontal axis.It also determines the deviation of curve, F. Dispatch function for a normal dispatch is presented by: (7) P X x e du u X ( ) Considering and F = 1, shown by N (0, 1), standard normal dispatch will be obtained.For , we can write: In which is achieved via N (0, 1), as: X


And dispatch function would be: If X and Y are continuous random variables, probabilistic dispatch function can be presented by: (12) Correlation and independency: Correlation and independency affect the behavior of random variables significantly.Correlation factor for X and Y, D, can be considered as standard covariance of X and Y: (14) Therefore, it can be mentioned that: C Covariance is dependent to measurement scale, while correlation factor is independent from measurement.C Correlation factor is similar in sign with covariance; thus, it changes between -1 and 1. C The values 1 and -1 for correlation factor represent a linear relation of the two variables.C If D = 0, variables X and Y have no correlation.
X and Y will be independent, if and only if: In which f x (x) and f y (y) are probabilistic density functions for X and Y.It is clear that the data of each variable does not affect the probability of other variable.It can be concluded that, if X and Y are independent random variables, they will be non-correlated.
Behavior of random variables: Considering X, Y and Z as random variables and a and b as real numbers, the following equations can be written based on exponent and covariance.These can also be used for random vectors.If Z = Ax + b, then: Considering a = 1 and b = in Eq. ( 18) reveals  X that has the same variance as X does, which If Z=X+Y, then: (20) For Z=XY: T and if Z can be explained as a nonlinear function of X, as: Then, exponent and covariance for Z can be achieved via Tailor series, as: D represents first-order derivative for matrix C x and stands for element (i, j) in matrix C x .Equation ( 24) states that exponent of Z can be corrected by covariance.

PROBABILISTIC LOAD DISPATCH
Changes in power of generators, changes of load power and distortions in voltages of PV buses are distortion factors in probabilistic eigenvalues analysis.Primary operating condition of generators and primary voltage of buses are determined by probabilistic load dispatch.In a normal dispatch, the necessary parameters are exponent and covariance.
In the N-bus system, the injected power vector, S, for voltage vector of V, which is the second-order function of voltage and can be written as (Stankovic and Lesieutre, 1989): Proper distribution of Eq. ( 25) around the exponent of would make the second-order term similar to G (V): ,..., ,..., ,..., ,..., In which JV is Jacobean matrix.Assuming and , exponent of S can be expressed as: (28) C Vi,j represents covariance of V i and V j .The effect of voltage covariance can be considered from Eq. ( 28) by adding covariance term.Therefore, mismatch power can be explained as: (29)  S S S   0 is the exponent vector of bus injected power.
Implementing linearization will lead to: Covariance matrix, C V , will then be achieved: Consequently, the correction equations for probabilistic load dispatch are: In probabilistic load dispatch calculations, due to different effects of exponent of voltage and covariance of injected power, Eq. (32-a) and (32-b) have different roles in calculation of C V from C S .While Eq. (32-a) is the inner loop in the calculation, Eq. (32-b) acts as the outer loop in the calculation.

DAMPING RATIO
For a specific eigenvalue 8 i = " i +j$, damping ratio, > i , can be calculated as: Exponent and variance for > i will be calculated using and C 8 .If > C remains in defined range (> i for this report ), probabilistic dispatch of damping ratio would be: (34) Determination of probability of stability: Equation ( 35) changes normal dispatch of x to a standard normal dispatch, which : x N X  ( , ) It should be mentioned that z 0 N(0,1).Dispatch function for z will be: (36)


In which z c is calculated by: (37) Probability in Eq. ( 36) can be achieved via probabilistic dispatch table.

NUMERICAL RESULTS
The 8-machine system presented in Fig. 1 is a part of the 36-bus system of PSASP software, in which DC links are eliminated.System parameters and power of buses are presented in Appendix.Buses 1 to 16 represent load buses, while buses 17 to 24 are generator buses.The remaining bus is considered as reference bus.Modeling


all generators with the fifth-order model of synchronous machine will lead to 84 eigenvalues, from which are 40 real eigenvalues and 44 complex eigenvalues.The real part is negative for all eigenvalues.Furthermore, four pairs of complex eigenvalues, shown in Table 1, have the minimum damping ratio.Therefore, only the results of these 4 modes are presented.
In analysis of probabilistic eigenvalues, changes in injected power of buses affect all operating parameters, such as current, voltage and rotor angle.On the other hand, these parameters are involved in forming matrix A. Therefore, they will affect eigenvalues.Using system admittance matrix, known by Y, for producing matrix A and for calculating sensitivity of eigenvalues, enable us to calculate eigenvalues via state matrix of A. With respect to the fact that eigenvalues are random variables, complex eigenvalues, e.g., 8 i = " i +j$ i , distribute within a specific area of s-plane.Discussing the distribution obedience can be done via probabilistic density function (p.d.f).With respect to its random values, the stability or instability of mode i th cannot be determined by sample values.
To demonstrate the point, p.d.f of "53 and "55, mentioned in Table 2, are shown in Fig. 2, which is based on normal dispatch.

Despite the lower exponent of in comparison with
 53 the exponent of , mode 55 is a rigid stable mode, which is due to the fact that all samples of "55 have negative values.At the same time, mode 53 is an instable mode.
Stability order for i can be evaluated by: (38) In which f(0) represents p.d.f for "i, where "i is a random variable.The system stability can be evaluated via probabilistic stability of critical eigenvalues.
For 8-machine system of Fig. 1, it is possible to assume different daily load curves for buses (SL1-SL9,PG1-PG7 and QG3-QG6).The voltage of reference bus changes with respect to A2 curve.All generators are modeled by third-order model of synchronous machine.The results reveal that, within 74 eigenvalues are 30 real eigenvalues and 22 complex eigenvalues.Some eigenvalues are presented in Table 3.The 74 th eigenvalue has a positive real part.

CONCLUSION
Probabilistic dispatch of an eigenvalue, as well as damping ratio can be determined via estimation of exponent and covariance.Figure 2 presents p.d.f for 53 rd and 55 th damping ratios which are shown in Table 3.
Damping factor probability and damping ratio are shown in the last two columns of Table 3, 4 and 5. "*stands for standard exponent, -"/F " .
For matrix A of 8-machine system, a comparison of eigenvalues has been implemented in two different conditions: with covariance correction and without covariance correction.The results are shown in Table 3 and 5.In our test system, without covariance correction exponent has a negative real value for all eigenvalues.The biggest value in this case is -0.085.
Implementing covariance correction will affect the results.For instance, the biggest value for exponent of eigenvalues will be 0.18 with covariance correction, as shown in Table 4.
To reduce the calculation task and time, in systems with numerous machines and, consequently, numerous eigenvalues, only a modicum number of eigenvalues will be investigated.

Table 3 :
Selected real eigenvalues for 8-machine system

Table 4 :
Selected complex eigenvalues for 8-machine system

Table A5 :
Data of transfer blocks for excitation system and governer in 8-machine system