Analysis of Small Oscillations in Complex Electric Power Systems

Allaev Kaxramon Raximovich, Makhmudov Tokhir Farkhadovich

Department “Electric Stations, Networks and Systems”, Tashkent State Technical University Named after Islam Karimov, Tashkent, Uzbekistan.

DOI: 10.4236/eng.2018.105017                

Abstract

In this article the mathematical model of complex regulated electric system in matrix form is developed. This mathematical model makes it possible to study the steady-state stability of a complex electrical system by determining the eigenvalues of the dynamics matrix. The model of an electrical system that reflects transient processes for small deviations is convenient, both algorithmically and computationally, in particular, in cases of their joint solution with steady-state equations—the equations of nodal voltages. The obtained results in the form of the eigenvalues of the matrix spectrum are qualitatively the same as the results of classical studies, which is a consequence of the adequacy of the proposed model and the correct reflection of the dynamic processes occurring in a real electrical system. In addition, the equations obtained are of independent importance for the analysis of various modes, including transient, electrical systems of any complexity.

Keywords

Electric System, Small Oscillations, Automatic Excitation Control, Matrix, Matrix Spectrum

1. Introduction

The present stage in the development of the power industry is characterized by the presence of large concentrated energy systems connected by relatively weak connections, in which the powers of distributed generation are actively included. The change in the composition of generation and the structure of power consumption leads to a decrease in the permanent inertia of the elements of the power systems, increasing the sensitivity of the parameters of the regime of the power system as a whole to small perturbations.

As is known [1] [2] [3] [4] , the dynamic properties of complex electrical systems can differ significantly from those of simple electric power systems (EPS), which is confirmed by numerous field and model experiments and computational and experimental studies. In a multi-machine electrical system, the choice of the parameters of the control devices is much more complicated than in the simplest EPS. Therefore, as a rule, in the case of a multi-machine EPS, one generator or one station is considered to be adjustable and the parameters of their automatic excitation controllers (AECs) are determined based on the task at hand―ensuring equal damping, the required stability factor, etc., and the parameters of AECs of other stations are assumed to be given, with constant emf. for a certain inductive resistance [2] .

In this paper we study the dynamic properties of electrical systems for small deviations (steady-state stability), described by linearized differential equations with constant coefficients.

Verification of the stability of power systems consists in determining the possibility of the existence of a stable regime with small perturbations of the parameters of the regime with given values of the parameters of the power system, the mode of generating sources, the load of node points, and the tuning of automatic mode control devices [3] .

The complication of modern electrical systems, the introduction of digital and logical control devices into their structure requires refined and in-depth studies of the modes of electrical systems. Such a problem can be successfully solved by matrix methods. The article suggests a matrix model of the electric system, resolved based on the absolute angles of the generators, which emphasizes the relevance of the task and the method for solving it [4] [5] [6] .

In the article the mathematical model and equations of multi-machine electric system, resolved concerning absolute angles of load of generators are received. On the basis of the obtained model, the results of steady-state stability analysis will be obtained using the example of a three-generator electric system.

Matrix equations of the elements of the EPS and the whole system were compiled on the basis of the most widely obtained equations of state variables [4] [5] [6] , which are small deviations of the mode parameters―the angles of the rotor load of the synchronous generator, busbar voltages, power and other operating parameters of the EPS. The considered matrix equations are used for analysis of transient processes and steady-state stability of EPS and for the synthesis of optimal parameters of regulators of synchronous machines operating in an electrical system.

2. Mathematical Model of Transients in a Complex Electrical System

This model describes the transient process in the electric system, taking into account the balance of the moments (powers) on the shaft of the i-th aggregate of the electric power system and has the form [5] :

${\text{d}}^2{\delta }_i/\text{d}{t}^2=\left({\omega }_0/T_j{}_i\right)\left[P_{Тi}-P_{Gi}\right],$ (1)

where ω₀ is the synchronous angular frequency; $T_{ji},{\delta }_i,P_{Ti},P_{Gi}$ -is the inertia constant of the i-th aggregate, the load angle of the i-th generator, the mechanical power of the i-th turbine, the electromagnetic power of the i-th synchronous generator, respectively.

The equation of electromagnetic power of the i-th synchronous generator in the positional idealization has the form [6] :

$P_{Gi}=E_1^2{y}_{ii}\sin {\alpha }_{ii}+{\displaystyle \underset{j=1,j\ne i}{\overset{n}{\sum }}{E}_i{E}_j{y}_{ij}}\sin \left({\delta }_{ij}-{\alpha }_{ij}\right),$ (2)

where E_i, E_j―emf. i-th and j-th synchronous generators; y_ii, y_ij―intrinsic and mutual conductivity of the network; α_ii, α_ij are complementary angles.

${\delta }_{ij}={\delta }_i-{\delta }_j,{\delta }_i={\delta }_{i0}+\Delta {\delta }_i,{\delta }_j={\delta }_{j0}+\Delta {\delta }_j,{\delta }_{ij}=-{\delta }_{ij},$ (3)

and beyond

$\begin{array}{c}\sin \left({\delta }_{ij}-{\alpha }_{ij}\right)=\sin \left[\left({\delta }_{i0}+\Delta {\delta }_i\right)-\left({\delta }_{j0}+\Delta {\delta }_j\right)-{\alpha }_{ij}\right]\\ =\sin \left[\left(\Delta {\delta }_i-\Delta {\delta }_j\right)+\left({\delta }_{i0}-{\delta }_{j0}-{\alpha }_{ij}\right)\right]\\ =\Delta {\delta }_i\cos {\beta }_{ij}-\Delta {\delta }_j\cos {\beta }_{ij}+\sin {\beta }_{ij},\end{array}$ (4)

where ${\beta }_{ij}={\delta }_{i0}-{\delta }_{j0}-{\alpha }_{ij}$ .

It should be noted that the derivation of formula (4) uses the obvious relationships:

$\sin \left(\Delta {\delta }_i-\Delta {\delta }_j\right)\cong \left(\Delta {\delta }_i-\Delta {\delta }_j\right)\text{andcos}\left(\Delta {\delta }_i-\Delta {\delta }_j\right)\cong 1,$

valid for small deviations in the load angles of generators.

After transformations (2), taking into account (3), (4), Equation (1) takes the form:

${\text{d}}^2{\delta }_i/\text{d}{t}^2=\left({\omega }_0/T_{ji}\right)\left[P_{Ti}-\left(E_i^2{y}_{ii}\sin {\alpha }_{ii}-{\displaystyle \underset{j=1,j\ne i}{\overset{n}{\sum }}{b}_{ij}\Delta {\delta }_j}+b_{ii}\Delta {\delta }_i+c_{ij}\right)\right],$ (5)

and taking into account the parameters of the initial regime and the relation ${\delta }_i={\delta }_{i0}-\Delta {\delta }_i$ , finally leads to a differential equation in the deviations:

${\text{d}}^2\Delta {\delta }_i/\text{d}{t}^2=\left({\omega }_0/T_{ji}\right)\left[{\displaystyle \underset{j=1,j\ne i}{\overset{n}{\sum }}{b}_{ij}\Delta {\delta }_j}-b_{ii}\Delta {\delta }_i\right],$ (6)

where

$b_{ij}=a_{ij}\cos {\beta }_{ij},a_{ij}=E_i{E}_j{y}_{ij},b_{ii}={\displaystyle \underset{j=1,j\ne i}{\overset{n}{\sum }}{b}_{ij},}{c}_{ij}={\displaystyle \underset{j=1,j\ne i}{\overset{n}{\sum }}{a}_{ij}\sin {\beta }_{ij},}$

$P_{Ti}-\left(E_i^2{y}_{ii}\sin {\alpha }_{ii}+c_{ij}\right)=0.$

In the case of the damper contours of the rotor of the i-th synchronous generator, Equation (6) takes the form:

${\text{d}}^2{\delta }_i/\text{d}{t}^2=\left({\omega }_0/T_{ji}\right)\left[{\displaystyle \underset{j=1,j\ne i}{\overset{n}{\sum }}{b}_{ij}\Delta {\delta }_j}-b_{ii}\Delta {\delta }_i-P_{di}\left(\text{d}\Delta {\delta }_i/\text{d}t\right)\right],$ (7)

where P_di is the coefficient of the generalized damper moment of the i-th generator.

If the deviation of the emf is taken into account. i-th synchronous generator, Equation (7) takes the form [7] :

${\text{d}}^2{\delta }_i/\text{d}{t}^2=\left({\omega }_0/T_{ji}\right)\left[{\displaystyle \underset{j=1,j\ne i}{\overset{n}{\sum }}{b}_{ij}\Delta {\delta }_j}-b_{ii}\Delta {\delta }_i-P_{di}\left(\text{d}\Delta {\delta }_i/\text{d}t\right)-\left(\text{d}{P}_i/\text{d}{E}_{qi}\right)\Delta E_{qi}\right].$ (8)

The peculiarity of Equation (8) is that it is allowed with respect to the absolute angles of the system generators and, for example, for the three-generator electric system has the form [6] [7] [8] :

$\begin{array}{l}{\text{d}}^2\Delta {\delta }_1/\text{d}{t}^2=\left({\omega }_0/T_{j1}\right)\left[-b_{11}\Delta {\delta }_1+b_{12}\Delta {\delta }_2+b_{13}\Delta {\delta }_3-P_{d1}\left(\text{d}\Delta {\delta }_1/\text{d}t\right)-\left(\text{d}{P}_1/\text{d}{E}_{q1}\right)\Delta E_{q1}\right],\\ {\text{d}}^2\Delta {\delta }_2/\text{d}{t}^2=\left({\omega }_0/T_{j2}\right)\left[-b_{21}\Delta {\delta }_1+b_{22}\Delta {\delta }_2+b_{23}\Delta {\delta }_3-P_{d2}\left(\text{d}\Delta {\delta }_2/\text{d}t\right)-\left(\text{d}{P}_2/\text{d}{E}_{q2}\right)\Delta E_{q2}\right],\\ {\text{d}}^2\Delta {\delta }_3/\text{d}{t}^2=\left({\omega }_0/T_{j3}\right)\left[-b_{31}\Delta {\delta }_1+b_{32}\Delta {\delta }_2+b_{33}\Delta {\delta }_3-P_{d3}\left(\text{d}\Delta {\delta }_3/\text{d}t\right)-\left(\text{d}{P}_3/\text{d}{E}_{q3}\right)\Delta E_{q3}\right].\end{array}$ (9)

The equations of electromagnetic transient processes in the excitation circuit of the i-th synchronous machine in the deviations have the form [1] [2] [4] :

$T_{di}\left(\text{d}\Delta E_{qi}/\text{d}t\right)=\Delta E_{qi}-\Delta E_{qei},$ (10)

$T_{ei}\left(\text{d}\Delta E_{qei}/\text{d}t\right)=\Delta U_{АECi}-\Delta E_{eqi},$ (11)

$T_{pi}\left(\text{d}\Delta U_{AECi}/\text{d}t\right)=\Delta e_i-\Delta U_{АECi},$ (12)

where $T_{di},T_{ei},T_{pi}$ ―the transition time constant of the excitation winding, the exciter time constant, the automatic excitation controller, respectively; $\Delta E_{qi},\Delta E_{qei},\Delta U_{АECi}$ ―deviations of the synchronous, forced emf. and the voltage at the output of the automatic excitation controller, respectively. The formation of signals via the AEC Dе_i channels in an idealized form (provided that the constant times of the differentiating elements of the AEC are considered to be zero) can be represented in the form [9] [10] :

$\Delta e={\displaystyle \underset{1}{\overset{k}{\sum }}{k}_{0Pk}\Delta P_k+k_{1Pk}\left(\text{d}\Delta P_k/\text{d}t\right)+k_{2Pk}\left({\text{d}}^2\Delta P_k/\text{d}{t}^2\right)},$ (13)

where k_0Pk, k_1Pk, k_2Pk are the gain factors of the AEC on the deflection channels, the first and second derivatives of the regime parameters ΔP_k, respectively, k is the number of adjustable mode parameters.

The advantage of Equations (7) and (8) is their dependence on the deviations of the absolute load angles of the generators (Dd_i) rather than the relative angles (Dd_ij), which provides computational convenience, since these equations can be joined to the equations of node voltages, whose solutions give absolute angles [11] [12] [13] [14] [15] .

For small perturbations of the parameters of the regime, after the corresponding transformations (1) - (13), it is possible to obtain a generalized block matrix А_S of size (4n × 4n) for the dynamics of an electrical system with n generators having automatic strong excitation regulators in the form:

$A_{\Sigma }=\left[\begin{array}{cccc}{0}_{n\times n}& I_{n\times n}& 0_{n\times n}& 0_{n\times n}\\ A_{21\left(n\times n\right)}& A_{22\left(n\times n\right)}& A_{23\left(n\times n\right)}& 0_{n\times n}\\ 0_{n\times n}& 0_{n\times n}& A_{33\left(n\times n\right)}& A_{34\left(n\times n\right)}\\ A_{41\left(n\times n\right)}& A_{42\left(n\times n\right)}& 0_{n\times n}& A_{42\left(n\times n\right)}\end{array}\right].$ (14)

The components of the matrix A_Σ are defined in [6] .

In this case, the vector-column of the state parameters containing the parameters of the electric system mode has the form:

$x={\left[\Delta {\delta }_1\cdots \Delta {\delta }_n⋮\Delta {\stackrel{\dot{}}{\delta }}_1\cdots \Delta {\stackrel{\dot{}}{\delta }}_n⋮\Delta {\stackrel{\dot{}}{E}}_{q1}\cdots \Delta {\stackrel{\dot{}}{E}}_{qn}⋮\Delta {\stackrel{\dot{}}{E}}_{qe1}\cdots \Delta {\stackrel{\dot{}}{E}}_{qen}\right]}^{\text{T}}.$ (15)

For example, for a three-generator EPS (Figure 1), assuming that the automatic excitation regulators react to voltage and load angle deviation of the generators (Dd_i, DU_Gi), as well as their first derivatives ( $\Delta {\stackrel{\dot{}}{\delta }}_i,\Delta {\stackrel{\dot{}}{U}}_{Гi}$ ).

The equation of the output of the automatic excitation controller for the i-th generator is:

$\Delta U_{АECi}=k_{0\delta Gi}\Delta {\delta }_{Gi}+k_{1\delta Gi}\left(\text{d}\Delta {\delta }_{Gi}/\text{d}t\right)+k_{0UGi}\Delta U_{Gi}+k_{1UGi}\left(\text{d}\Delta U_{Gi}/\text{d}t\right),$ (16)

In this case, the matrix AΣ takes the form [6] :

$A_3=\left[\begin{array}{cccccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ -{\omega }_{11}& {\omega }_{12}& {\omega }_{13}& 0& 0& 0& -\frac{\text{d}{P}_1}{\text{d}{E}_{q1}}\frac{{\omega }_0}{T_{j1}}& 0& 0& 0& 0& 0\\ {\omega }_{21}& -{\omega }_{22}& {\omega }_{23}& 0& 0& 0& 0& -\frac{\text{d}{P}_2}{\text{d}{E}_{q2}}\frac{{\omega }_0}{T_{j2}}& 0& 0& 0& 0\\ {\omega }_{31}& {\omega }_{32}& -{\omega }_{33}& 0& 0& 0& 0& 0& -\frac{\text{d}{P}_3}{\text{d}{E}_{q3}}\frac{{\omega }_0}{T_{j3}}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& \frac{1}{T_{d1}}& 0& 0& -\frac{1}{T_{d1}}& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& \frac{1}{T_{d2}}& 0& 0& -\frac{1}{T_{d2}}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{T_{d3}}& 0& 0& -\frac{1}{T_{d3}}\\ \frac{k_{0\delta 1}}{T_{e1}}& 0& 0& \frac{k_{1\delta 1}}{T_{e1}}& 0& 0& 0& 0& 0& -\frac{1}{T_{e1}}& 0& 0\\ 0& \frac{k_{0\delta 2}}{T_{e2}}& 0& 0& \frac{k_{1\delta 2}}{T_{e2}}& 0& 0& 0& 0& 0& -\frac{1}{T_{e2}}& 0\\ 0& 0& \frac{k_{0\delta 3}}{T_{e3}}& 0& 0& \frac{k_{1\delta 3}}{T_{e3}}& 0& 0& 0& 0& 0& -\frac{1}{T_{e3}}\end{array}\right],$ (17)

Vector-column of the space of states of parameters of the EPS regime:

$x={\left[\Delta {\delta }_1\Delta {\delta }_2\Delta {\delta }_3\Delta {\stackrel{\dot{}}{\delta }}_1\Delta {\stackrel{\dot{}}{\delta }}_2\Delta {\stackrel{\dot{}}{\delta }}_3\Delta {\stackrel{\dot{}}{E}}_{q1}\Delta {\stackrel{\dot{}}{E}}_{q2}\Delta {\stackrel{\dot{}}{E}}_{q3}\Delta {\stackrel{\dot{}}{E}}_{qe1}\Delta {\stackrel{\dot{}}{E}}_{qe2}\Delta {\stackrel{\dot{}}{E}}_{qe3}\right]}^{\text{T}}.$

As can be seen, the generalized matrix A3 of the dynamics of the electrical system, consisting of 3 generators, is formed from the parameters of the system

Figure 1. Diagram of a three-generator electrical system.

mode and the automatic regulation of the excitation of machines, and therefore fully characterizes the transient processes in this EPS. The matrix A3 is rather sparse, which is typical for a complex system containing n generators, so this fact determines the computational advantages of the proposed mathematical model in the calculation and experimental studies of EPS.

3. Example

As an example, consider the matrix (17) of the intrinsic dynamics of the three-generator electric system A3 (Figure 1). Table 1 shows the technical parameters of the electrical system under consideration.

Where X_d―synchronous inductive resistance of the generator along the longitudinal axis;

X_c―resistance of branches (inductive resistance of power lines).

The result of the calculation is shown below.

$A_3=\left[\begin{array}{cccccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ -\text{235}\text{.3565}& \text{155}\text{.8297}& \text{99}\text{.8200}& 0& 0& 0& -\text{0}\text{.4167}& 0& 0& 0& 0& 0\\ \text{104}\text{.5232}& -\text{120}\text{.3483}& \text{51}\text{.6890}& 0& 0& 0& 0& -\text{0}\text{.5214}& 0& 0& 0& 0\\ \text{50}\text{.6447}& \text{39}\text{.2440}& -\text{94}\text{.1103}& 0& 0& 0& 0& 0& -\text{0}\text{.5844}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1.25& 0& 0& -1.25& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& \text{1}\text{.4286}& 0& 0& -\text{1}\text{.4286}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1.6667& 0& 0& -1.6667\\ 20& 0& 0& 2& 0& 0& 0& 0& 0& -2& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& -2.5& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& -\text{2}\text{.2222}\end{array}\right].$

The spectrum of the matrix of the intrinsic dynamics of the three-generator EPS A3 at the selected regime and system parameters is equal to: −0.0012 ± 17.9944i; −0.0001 ± 11.9462i; 4.108; −4.0995; −1.9571; 1.2009; 1.4286; 1.6667; −2.5; −2.2222. The electrical system is not stable, since the model under study has positive poles. The inclusion of AECs of other generators by deviation and the first derivative of the angle and voltage of the generators will ensure the stability of the system under study.

Table 1. Parameters of elements of the three-generator power system.

Figure 2. 3D-visualization of the spectrum of a three-generator electrical system with a Hurwitz matrix A₃.

Similar studies of the stability of complex electrical systems were carried out in [2] [5] [8] , but they used mathematical models derived from the relative angles of the load of synchronous generators. In this study, the results were obtained with respect to the absolute angles of loads, which allows us to further determine and investigate a particular generator, the first approaching the stability limit.

Figure 2 shows a 3D visualization of the matrix a pseudo-spectrum obtained with the MATLAB software module EigTool. Note that the EigTool module is a software product developed by Oxford University [15] . Its methodological basis is the method of Arnoldi computation of subspaces of A.N. Krylov [16] . The horizontal axes in Figure 2 correspond to the axes of the complex plane. Logarithm of the norm of the resolvent function was postponed along the vertical axis. Peaks localize the eigenvalues of the matrix.

4. Conclusion

The dynamic properties of complex electrical systems can differ significantly from the properties of simple EPS, which is confirmed by numerous full-scale and model experiments and computational and experimental studies [12] [13] [14] . In a multi-machine electrical system, the choice of the parameters of the control devices is much more complicated than in the simplest EPS. Therefore, as a rule, in the case of a multi-machine EPS, one generator or one station is considered to be adjustable and their AEC parameters are determined proceeding from the task in hand―providing equal damping, the required stability factor, etc., and the parameters of the AEC generators of other stations are selected from the need to provide stability of the entire system and damping of possible oscillations of the regime parameters [17] . Therefore, the introduction of AECs from other generators requires additional studies on the choice of regulatory parameters (synthesis), which is the subject of further research.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Abdellatif, B.M. (2018) Stability with Respect to Part of the Variables of Nonlinear Caputo Fractional Differential Equations. Mathematical Communications, 23, 119-126. http://www.mathos.unios.hr/mc/index.php/mc/article/view/2326/526
[2]	Klos, A. (2017) Mathematical Models of Electrical Network Systems: Theory and Applications—An Introduction. Springer International Publishing AG, Berlin. https://doi.org/10.1007/978-3-319-52178-7
[3]	Holali, K.D., Efimov, D. and Richard, J.-P. (2016) Interval Observers for Linear Impulsive Systems. 10th IFAC Symposium on Nonlinear Control Systems (NOLCOS 2016), Monterey, 23-25 August 2016, 867-872.
[4]	Kothari, D.P. and Nagrath, I.J. (2003) Modern Power System Analysis. McGraw Hill Education, New York.
[5]	Misrihanov, M.Sh. (2004) Klassicheskie i novye metody analiza mnogomernyh dinamicheskih system [Classical and New Methods of Analysis of Multidimensional Dynamic Systems]. Energoatom Publishing, Moskow (In Russian).
[6]	Allaev, K.R. and Mirzabaev, A.M. (2016) Matrichnye metody analiza malyh kolebaniy elektricheskih system [Matrix Methods for the Analysis of Small Oscillations of Electrical Systems]. Fan va texnologiya Publishing, Tashkent (In Russian).
[7]	Gotman, V.I. (2007) Common Algorithm of Static Stability Estimation and Computation of Steady States of Power Systems, Power Engineering, 311, 127-130. http://www.lib.tpu.ru/fulltext/v/Bulletin_TPU/2007/v311eng/i4/30.pdf
[8]	Fazylov, H.F. and Nasyrov, T.H. (1999) Ustanovivshiesya rezhimi elektroenergeticheskih sistem i ih optimizaciya [Established Regimes of Electric Power Systems and Their Optimization]. Moliya Publishing, Tashkent (In Russian).
[9]	Kovalenko, S., Sauhats, A., Zicmane, I. and Utans, A. (2016) New Methods and Approaches for Monitoring and Control of Complex Electrical Power Systems Stability. IEEE 16th International Conference on Environment and Electrical Engineering (EEEIC 2016), Florence, 7-10 June 2016, 270-275.
[10]	Allaev, K.R., Mirzabaev, A.M., Makhmudov, T.F. and Makhkamov, T.A. (2015) Matrix Analysis of Steady-State Stability of Electric Power Systems. AASCIT Communications, 2, 74-81.
[11]	Irwanto, M., et al. (2015) Improvement of Dynamic Electrical Power System Stability Using Riccati Matrix Method, Applied Mechanics and Materials, 793, 29-33.
[12]	Anderson, P.M. and Fouad, A.A. (2003) Power System Control and Stability. 2nd Edition, Willey-Interscience A John Wiley & Sons Inc., Hoboken.
[13]	Kunder, P. (1993) Power System Stability and Control. McGraw-Hill, Inc., New York.
[14]	Pal, M.K. (2007) Power System Stability. Edison, New Jersey.
[15]	MATLAB (2001) User’s Guide. Reference Guide. The Math Works, Inc., Natick.
[16]	Wilkinson, J.H. (1988) The Algebraic Eigenvalue Problem. Clarendon Press, Oxford.
[17]	Albertos, A.P. and Sala, A. (2007) Multivariable Control Systems. Springer, London.