- Email: [email protected]
Fast determination of transient stability regions for the second Lyapunov method
Fast determination of transient stability regions for the second
Lyapun0vmeth0d J Machowski Institute of Power Systems, Technical University of Warsaw, 00-662, Warsaw, Koszykowa 75, Poland
In the present paper, methods o f finding the criterial value fi)r the Lyapunov function are considered. A new method is described which efiminates generation ofall possible starting points. The criterial starting point has been determined by findhlg the minimum cut set for a certain network representing transfer synchronizing powers o f generators. The proposed method is faster than techniques published so far.
I. Introduction Although numerous publications have dealt with the application of the direct Lyapunov method to power system transient stability studies, the practical use of this method is far from satisfactory. One of the main drawbacks is the long time necessary to determine the criterial value Vc for the Lyapunov function referring to the transient stability region. When the Lyapunov function in the form of total energy is applied, the criterial value corresponds to the least value of energy necessary to divide the system into two generator groups operating in an asynchronous way. So the critedal value is equal to the value of the Lyapunov function in the lowest saddle point, called the criterial point. Methods of finding the V~. value can be divided into four groups. i. In the starting-point method, the criterial point is found by: a. generating the ith starting point, near to which an equilibrium point may be located b. determining the coordinates of the equilibrium point by solving the corresponding nonlinear equations using, for example, the Newtonian method c. calculating the Vi value for the determined equilibrium point and returning to stage a until all the possible starting points are over d. selecting I7,. = rain Vi. This procedure requires (2 n I 1)-fold solutions of (n - 1) nonlinear equations, and the same number of calculations Received: 9 July 1980
Vol 2 No 4 October 1980
referring to the value of the Lyapunov function. These operations are very time-consuming, and for any power system of reasonable size the procedure becomes computationally prohibitive. A simplified version of the method proposed in Reference 1 consists in displacing stage b from the second position to the last, in order to eliminate the multiple solution of the nonlinear equations and thus make all the computations nmch faster. If the reference generator is selected properly (e.g. the generator remarkable for its highest self-synchronizing power is chosen) the results obtained are usually true. However, in spite of this simplification, the computing time is too long, and still far from satisfactory for practical purposes. Some other simplified versions of this method, proposed in References 2 and 3, will be discussed in Section IV. ii. The method proposed in Reference 4 consists of two stages. During the first stage, the operation should begin from a stable point and move in small steps from one point on an equiscalar surface to the next point having a greater V-function value. In this case, every subsequent point may be determined by minimization of a certain function illustrating the path of least ascent. All these calculations are ceased after an attraction area of the saddle point has been reached. At this moment, the second stage of the operation should be started. It consists in strict determination of the saddle point coordinates by solution o f ( n 1) nonlinear equations. iii. The criterial value Vc is formulated as the following constrained optimization problem: find the minimum of tile V-function on the surface ~ = 0, where 4/is a certain function obtained from the Rauth-Hurwitz condition s or simply the directional derivative of the Lyapunov flmction. The authors of Reference 6 applied the Monte Carlo method to solve the above task. iv. The fourth method differs from the previous three in having a quite different approach. It consists in some particular decomposition of the system and the application of the vector Lyapunov function 7. The method proposed at a later stage in the present paper
0142-0615/040207-07 $02.00 © 1980 IPC Business Press
207
belongs, in general, to the first group. However, the simplifications made allow for relatively fast determination of the criterial point, without generation of all the possible starting points.
II. Conservative model of power system For transient stability studies using the direct Lyapunov method, the generators are represented by constant voltage E ' behind reactances X~ as well as by the equations of generator rotor inovement which can be written as follows: d2~i
- 2 = Pmi Mi -dt
Pi
at i = 1, . . . , n
I I. 1 Application o f coherent generator aggregation In the majority of large electrical power systems, after shortcircuit localization, groups o f coherent generators can be found. Each of them can be replaced by an equivalent generator, by means o f aggregation. The equivalent generator is connected to the system by certain fiction lines. The values of conductance for these lines depend to a great extent on the differences in the angles 6 i of the aggregated generators. In the case of there being great differences in these angles and any reasonable aggregation method being applied, the conductance values are too high to be neglected. Such conductances should be inw~lved in the Vp function.
( 1)
where M i is the coefficient o f inertia, 8i is the generator rotor angle, Prni is the mechanical power (constant), and Pi is the electrical power. Load nodes can be eliminated by assuming that reactive powers o f loads can be replaced by susceptances and that active powers of loads can be replaced by constant injection currents. If this assumption is taken into account, then, after elimination of load nodes, the values of transfer conductances are so small that they can be neglected. Thus, the following formula for electrical power can be derived:
Later, it will be proved that when the aggregation is made with the aid of the Zhukov method (see Appendix 1) the conductances introduced by aggregation are antisymmetrical and form potential moments only, so that these momellts can be exactly taken into account in the calculation of t5,. The Zhukov aggregation is remarkable for one of its properties: if the admittance nodal matrix Yc is of the . . . . *T *T antlhermltlan t y Pe , l • e . YRR =. Y RR, Y~4~I * r ..' , = . Y I ~ M , YmR . . . . YRM, then the matrix Yg obtamed alter aggregation is of the antihermitian type, too. It appears so because"the matrix YRR does not change and y~er
=
(yRM e)*T
= e
* r ,t
,,,7' IR,M
) = e
*'r YMR = YeR
The above is true for any number of aggregated generator groups.
/l
Pi = Psi + Z bi] sinSi/
i = 1. . . . . n
(2) It should be noted that the matrix YG obtained after elimination of the load nodes and resulting negligence of the conductance value is of the antihermitian type (only symmetrical imaginary elements exist).
where Psi is the power o f equivalent loads resulting from elimination o f load nodes - so called self-power and
bi] = BijEiE j The mathematical model in equations (1) and (2) is conservative, as the condition OPi/38 i = aPl/O6i referring to the existence o f potential Vp is valid. In this model, the potential is defined as follows: (Pmi Pi) = aVp/O6i and hence after integrating equation (1), the Vp function may be found by the following formula
Thus, it has been stated that after aggregation, the nodal matrix to describe the equivalent system is of the antihermitian type, i.e. the following expressions are true for every fiction line:
=
rT,
Or Vp = -
~ (Pmi i=1
~..
n
- Psi)(6i
-- 8 s )
Bij = Bii
1
y
%(cosa;j - cosS,';)
and (3)
i=l ]=i+1
where 8is ~s • the coordinate o f the postfault equilibrium point. Each stationary point of the Vp function is considered as an equilibrium point for the system of equation (1) and so the equilibrium points are to be classified (according to the shape of the potential field Vp) into the minimum point, saddle points, and maximum points. The saddle point is classified as the point of maximum type with reference to the k variables and as the point o f minimum type with reference to the remaining variables. Such a point will be called the unstable equilibrium point of the k type.
208
Gi/ = _ Gii The electrical power for the equivalent system can be determined by using the following formula:
Pi = Psi + ~, (bij singq + gij cos6ij)
{4}
at g i / = GijEiEI. As the conductances introduced within the aggregation are antisymmetrical, the condition OPi/~8/= OP//OgJi referring to the existence of potential Vp is valid, and after integra-
Electrical Power & Energy Systems
tion of equation (1) the Vp function is as follows: n
Vp =
~
X (Pmi - Psi)(Si - 6~) i
1
n -- 1
X i~l]
x [bd(cos6ii - cos6Si ) -gi/(sin6i/
i-1
sin6Si)]
(5)
II I-. Proposed method If the trajectory of generator rotor movement crosses the crest of the type k saddle point, k generators lose their synchronism in relation to the remaining (n k) generators. In this respect, searching for the criterial saddle point is equivalent to looking for such a division of generators (into two groups) that will require minimum energy. Function (3) or 15) determining the energy in question may be expanded (about stable point 6 s) into a Taylor series. Taking into consideration that the gp function and its gradient in point 6 s are equal to zero, it appears that
where A8 : 6 8 s and H is the square, symmetrical Hessian matrix of potential energy H = [OZVp/aSi OOj] and, because 3 Vp/O8 i = -Pi, it may be considered as a Jacobian matrix for equation system (2) H = [3Pi/36il. Approximation (6) is true in the neighbourhood of the stable point 6 s only and should not be used for an approximate calculation of the Vp value for unstable equilibrium points located at some distance from the stable point 8 s. However, after performing calculations for a number of exemplary test models of electric power systems, some valuable properties of approximation (6) were noted. It shows the least value at the same unstable equilibrium points for whidl the value of function (3) or (5) is the least, too. Apart from that, in most cases, approximation (6) permits proper classification of all the unstable equilibrium points, according to their increasing potential energy (3) or (5). The application of these properties will now be discussed. I I I. 1 Unstable e q u i l i b r i u m p o i n t s k = 1 In this case, the following vector can be assumed: A8 ~ [0 . . . . .
O, A6i, 0 . . . . .
111.2 Unstable equilibrium points k > 1 The situation in which k > 1, i.e. when the group of generators {L} loses its synchronism with respect to the remaining generators {R}, may be considered as the case when the synchronism is lost by one generator e which replaces the whole group {L}. The self-synchronizing power for such a generator results from principles of aggregation and is equal to the sum of the transfer synchronizing powers between the group of generators {L} and the group {R} (see Appendix 2)
Hee:HL/~:
Y. iCL
E
¢8)
( H~i)
]~R
The Hee value will be called the equivalent synchronizing power. Assuming that a set {L} consists of one element only, the fornmla in equation (8) becomes simplified to the formula given in equation (7). The group of generators {Lc} for which the equivalent synchronizing power (8) does not exceed the self-synchronizing power of each generator belonging to the group will be called the compact group. The term compact group is justified by the fact that in order to meet its definition the transfer synchronizing powers inside the group are much higher than the transfer synchronizing powers between the generators of this group and the remaining generators of the system. It is also worth noting that each generator satisfies the above definition and forms a single-element (trivial) compact group. In conclusion, it should be stated that determination of the unstable equilibrium point of the k > 1 type, showing the least Vp value (or, more exactly, a suitable division of generators) requires the establishment of a compact group showing the least equivalent synchronizing power (8). Such a group will be called the criterial compact group, II 1.3 A l g o r i t h m s f o r research o f criterial c o m p a c t group
Let the generator nodes of the examined system be the vertices {W} of the graph G with the lines {U}. Let the real numbers ( Hij ) describe the lines of the graph, and the real n u m b e r s H u describe the vertices. The subset {P} of set{U} can be defined as follows:
{P} = [e{u} :x e { wl},y e{W2}l
O] 7"
where
Approximation (6) has then been simplified to the form : i A , ~ 2T/ Vp ~ 2 ~ i * - i i .
Thus, the energy necessary to cause the loss of synchronism by generator is approximately proportional to its selfsynchronizing power Hii where
{Wl} U {W2} = {W}
{~]1} ("1 {~¥2} = ~)
On removal of the lines representing cut set {P} there is no pail1 from {W1} to {W2} and in this way the cut set {P} may be identified by the pair of subsets {Wl} and {W2}. The real function for the cut set {P}
Hii = ~ ( Hi i)
(7)
Thus. the point of the k = 1 type showing the least Vp value is easily determined by means of finding the generator showing the least self-synchronizing power (see Section Ill.5.1).
Vol 2 No 4 October 1980
H(P)= ~ ~ (-Hxy) x~iw,} y~iw2)
(9)
may be defined as the maximum flow of the cut set. The cut set for which the value of equation (9) is the least possible one is considered as the minimal cut set.
209
The above definition (9) is illustrated by Figures 1 and 2. Formulas (8) and (9) are identical, and hence the determination of the criterial compact group is equivalent to obtaining the minimal cut sets between every pair of vertices. Such a task can be solved very quickly even for large networks through use of the Gomory and Hu theorem s, by nreans of the (n - l )-fold application of the labelling procedure 9, m or by use of the Korzanov procedure n.
I
2
4
3
5
6
7
Anita
~"pi
1 1961<1 317 2693 3 566 1581 366 157}I i ; _ l ~ ~ 2 11315 725011939 2 531 i 122 241 112 [ (L(~90irI4095! 3 267t t933 9 8 9 ~ 3 083 1368 580 270ii oioi211i7328 i 4 3562 2530 307314 54814663 495 23n 0.094!28 198 5 11 577 1 121 1366 4 661 9046[ 220 103{ f/(:igg~-1762~ 6 i 334 240 580 494 219 3 6 5 4 1 1 7 > n.030~7t80! 7 i 155 111 269 ,02,787 =6, t t 704i
I 11.4 Calculation of coordinates for criterial point On determination of the minimal cut set, the criterial starting point is established by providing the angles of generators { W=} with the value ~is, and the angles of generatars{W1} with the value (-+lr --6s). The plus sign is given for zr in the case y~
(e~; - P~) > o
',
i~{w,} i.e. when the generators of the group {W~} export the power outside. The minus sign is used if the generators of the group { W1} import the power so that
}2
O % ~ - G ; ) < o.
iE{W,}'
2
(
%
F i g u r e 3 Example l: system containing single-element compact groups only
The above arrangements of plus and minus signs have been proved and described in detail in Reference 12. The coordinates of the criterial point are obtained by solving equation (1) in the neighbourbood of the starting point while d25i/dt 22 O. For this the Newton-Rap!tson or Davidon-Powell-Fletcher method can be used.
Figure 1 Example of cut set w i t h k = card. {Wt} = 1; maxim u m f l o w for cut set is H(P) = - ( H t 2 + H t 3 + H t 4 + H l s )
"~v
111.5 Examples The data for the examples is given in Figures 3, 4 and 5 where the numbers of the rows and columns correspond to the numbers of particular generators. The values given in the upper triangular matrix refer to the parameters bii included in equation (2). The values given in the lower triangular matrix together with the diagonal elements refer to transfer synchronizing powers as well as to self-synchronizing powers. The vector on the right contains the angles
For better illustration of the proposed method, the dendrites have also been included in all the examples; the dendrites, which are the flow-equivalent trees, enable fast visual detennination of the maximal flows for all possible cut sets. 111.5.1 Example 1 This system contains the single-element compact groups only. Thus, the criterial point is o f tire k = 1 type and corresponds to the loss of synchronism of generator 7. The value of function (3) calculated for this point is equal to Vc = V,~7 = 4 704. hr this exmnple, the values of the Vp functions for all unstable equilibrium points of k = 1 are also given.
F i g u r e 2 Example of cut set w i t h k = c a r d . { W ~ } = 2 ; m a x i m u m f l o w for cut set is H(P) = -- (Hi3 + HI4 + H i s + H23 + H24 + H25)
210
111.5.2 Example 2 This system contains tile doubleelenrent compact group {6, 7} which is the criterial compact
Electrical Power & Energy Systems
1
2
1 ~1276
3
4
3160 3299
5
6
7
Angles
1732
280
230
0.070
1713
902
148
121
0.072
3 3159
16621110713415
1794
608
473
0.043
4 3291 5 1731 6 279
1709 3 4 1 2 1 3 7 8 6 1 4 6 7 3
389
319
0.000
168 203 4 2 9 6 1 2 6 7 9
0.049 0.005
167 2 6 7 t
0.069
2 1264 5797]_1673
7
230
895 1794 4 6 6 7 147 607 389 120
473
318
94571 204
3979
ally, attention should be paid to the fact that the advantage of the proposed method increases as the number of generators increases. Further, it should be pointed out that white preparing the above program, not all the possibilities of computation time reduction were taken into account and so the times given in Figure 6 can be reduced further, e.g. by the application of the algorithm in Reference 1 1 instead of Reference 13.
IV. Other simplified methods
-.~
">--< ,,,jC-J
~c~ 1 ~ . . .
/ J J J
/
/
J
Figure 4 Example 2: system containing criterial compact group{6, 7}
group. In this case, the criterial point is of the k = 2 type and corresponds to the loss of synchronism of generators 6 and 7. The function (3) value of this point is equal to Vc = 5 864. The minimal point of the k = 1 type corresponds to the loss of synchronism of generator 7 which belongs to the criterial compact group. 111.5.3 Example 3 This system contains (as in Example 2) the criterial compact group {6, 7}, provided that Ve = 5 782. In this case, the minimal point of the k = 1 type at Ve = 6 843 corresponds to the loss of synchronism in generator 2, which does not belong to the criterial compact group. Thus, for this example, the Gupta and El-Abiad assumption is not valid. The values Vp in points of the k = 1 type corresponding to the criterial compact group are Ve6 = 7 756 and Vp7 = 7 207.
In accordance with the observations mentioned in Reference 2, it was noted that in typical systems the criterial starting points are of the k = 1 type. On the basis of these observations, the suggestion has been made of limiting the calculations to 'n' starting points of the k = 1 type. This property is true for the majority of electrical power systems, and more so when considering systems after aggregation of coherent generators. However, some cases can be found where the criterial points are of the k > 1 type in spite of aggregation. This results from the fact that the compact group can be coherent only if the shortcircuit is situated at some distance and involves approximately the same accelerations. If the coherency criteria are not satisfied because of nonuniform accelerations or relative distances 14, ls, the program of automatic reduction leaves the compact group in nonaggregated state, and the criterial point can be of the k > 1 type. The other simplified method is given in Reference 3, and its basis may be explained briefly as follows. Let {Gx} be a set of generators losing their synchronism in minimal 1
2
3
4
6
7
280
230
0.139
685 3 821 l 1 039 1 078 760 147 120 3 1 5 7 1 0 3 6 1 0 4 4 8 1 3 4 1 6 1794 577 472 3 4 3 2 9 2 1077 3 4 1 3 1 3 15614674 388 319 757 1794 4 6 6 8 9 3 2 1 ] 204 168 5 1731 279 147 576 388 203 4 262 {2 689 6 230 119 472 318 168 2669 3 9 7 6 7
0.035
1 9 3 7 4 I 689 3159 3 2 9 9
5
1731
Angles
2
I [
0.112 0.069 0.118 0.064 0.138 //
~//
\\\ _ ////
111.6 Computation time On the basis of the above method, the programme has been prepared for a CDC 3170. In the program, the algorithm in Reference 13 and the standard Newton-Raphson procedure were used. The computations have been performed for many different test systems, and compared with the results of calculations using the exact starting points method and with the results of the Prabhakara and EI-Abiad method (the generator with the highest self-synchronizing power being taken as the reference generator). The results of calculations have confirmed the exactness of the proposed method. Average values of computation times are given in Figure 6. The time requirenrent of the computations for the proposed method is much less than for the remaining two methods. Addition-
Vol 2 No 4 October 1980
// /
/
/
/
Figure5 Example 3: system containlngcriterialcompact group{6, 7}, provided that Vc = 5 782
211
4o 3o
Ribbens-Pavella, M and Lemal B 'Fast determination of stability regions for online transient power-system studies' Proc. IEEE Vol 123 No 7 (1976)
I I
3 Gupta, C L and EI-Abiad, A H 'Determination of the clOSest unstable equilibrium state for Lyapunov methods in transient stability studies' IEEE Trans. Power Appar. & Syst. Vol 95 (September/October 1976)
I
I I I
Tavora, C J and Smith, O J M 'Stability analysis of power systems' IEEE Trans. Power Appar. & Syst. Vol 91 (May/June 1972)
I I
E
E
Weiman, M J 'Regions of permissible deviations in determining synchronous transient stability of power systems by Lyapunov's second method' (in Russian) Elektrotekhnika No I (197 I)
I / I
Strakhov, S V and Weiman, M J 'The present state of
f I
(:z
E
o (.2
development and possibilities for practical application of Lyapunov's second method in determining the transient stability of power systems' (in Russian) Elektrotekhnika No 10 (1977)
/
10
//
// O
5
/ /
10 Number
15 of
f
20
1
f
/
f
Darwish, M and Fantin, J 'Stability analysis of large scale power systems using the comparison principle and vector Lyapunov functions'J. AppL Sci. & Eng. No 3 (1978)
25
generators
Christofides, N Graph t h e o r y - an algorithmic approach
Figure6 Average values of computation times of CDC 3170; - - - - exact method, - . . . . . Prabhakara and EI-Abiad's method, - . . . . proposed method
saddle point o f k = x type. According to file Gupta and E1-Abiad assumption, {Gx_l} should be considered as tile subset of the set {Gx}. This assumption has been supported by examples; however, it is not generally true. This is illustrated in Example 3, where the minimal point of the k = 1 type corresponds to the loss of synchronism of generator 2, and the criteri',d point is determined by generators 6 and 7. The error made in Example 3, while using the method in Reference 3 to determine Vc, amounts to about 1g~/~, and can be much larger if the parameters in Example 3 are changed.
V. Conclusions One o f the main disadvantages o f the application of tile direct Lyapunov method to power system transient stability studies is the lack of methods that would enable relatively fast determination of the criterial value. In the presen! paper, a new technique has been proposed, in which tire time-consuming task of finding the criterial point has been replaced by determination of the minimal cut set of a cerlain network. The proposed method has been tried on many test systems and tire results obtained were considered to be correct. The computation times are then significantly shorter than with previously used methods.
References 1 Prabhakara, F S and EI-Abiad, A H 'A simplified
VI.
determination of transient stability regions for Lyapunov methods' IEEE Trans. Power Appar. & Syst. Vol PAS-94 (March/April 1975)
212
Academic Press, USA and UK (1975) Ford, L R and Fulkerson, D R Flows in networks Princeton University Press, USA (1962) 10 Edmonds, J and Karp, M 'Theoretical improvement in algorithmic efficiency for network flow problems' J. Assoc. Comput. Mach. No 19 (1972) 11 Hamacher, H 'Numerical investigation on maximal flow algorithm of Korzanov' Comput. No 22 (1979) 12 Maciejewski, Z 'Direct method of power system stability evaluation' Prace Instytutu Energetyki Zeszyt 6 (1978) 13 Bayer, G 'Algorithm 324, Maxflow' Commun. ACM No 11 (1968) 14
Lee, S T g and Schweppe, F C 'Distance measures and coherency recognition for transient stability equivalents IEEE Trans. Power Appar. & Syst. Vol PAS-92 (September/October 1973)
15 Machowski, J Simplified models for power system stability studies PhD Thesis, Technical University of Warsaw, Poland (1978)
VII.
Bibliography
Bayer, G 'Remark on algorithm 324' Commun. ACM No 16 (1973) Dinic, E A 'Algorithm for solution of a problem of maximal
Electrical Power & Energy Systems
flow in a network with power estimation' Sov. Math. DokL No 11 (1970)
where
OPi
Hq-
Zhukov, L A 'Construction of power system electromechanical equivalents' Izv. Akacl. Nauk SSSR Energ. & Transp. No 2 (1965)
3Pi aii = _
transfer synchronizing power
self-synchronizing power
Appendix 1. Aggregation of coherent generators Prior to aggregation of the selected generator group {M}, the system may be described by nodal equations in the following form :
Linearization of equation (11)results in the following equation :
[APn] = [ H R R APe 3 tHen LYMR
IM
YMM 3
or I(;=YGEG
or VeR
Yee J
lg = YgEg
(11)
~"
where Yne = YRM e, Yen = e*TyMR, Yee = e*TyMM e, and e = Ee~EM;Ee is any voltage at the centre of the angle (the above values are complex). Attention should be paid to the fact that this aggregation is equivalent to the aggregation considered by Mello, Podmore and Stanton 16, assuming that ideal equipotentializing transformers taken into account by the authors have a complex transformation ratio equal to e i = Ei/E e.
Appendix 2. Synchronizing power of equivalent generator Taking into account that E i = const., the following linearization of the equation (t0) can be made:
[a e R ] _ [ . . . APM3
HRM1
ApG=HGA6 G
L H M R HMMJ LAfMJ or (12)
Vol 2 No 4 October 1980
APg=HgAgg
(13) and it can be proved that
Assuming that the power of the equivalent generator is equal to the sum of aggregated generator powers, and assuming that the aggregation does not involve changes in currents and voltages of nodes {R}, it appears that =
or
(10)
where {R} are the remaining nodes of the system.
Ie
HRe] [ A f R ] Hee tASe J
HRe = HRM 1M, HeR = 1THMR, Hee = 1THMM 1M This means that the Zhukov aggregation may be considered as a nonlinear generalization of the linear Marconato-Marzio aggregation 17. As the sum of elements in the row of matrix H G is equal to zero, the synchronizing power nee may be represented in the following form
gee = 1MHMM1M = -ffIH
R1R =
Y. Y. ,~{M} j~{~}
-g,p
The symbols 1M, 1R represent vectors consisting of ones only, with cardinality {M} and {R}, respectively.
16 Mello, R W, Podmore, R and Stanton, K N 'Coherency based dynamic equivalents: Applications in transient stability studies' Proc. PICA Conf. (1975) 17 Marcenato, R and Marzio, L 'A dynamic model of Italian power system based upon area equivalents' Energ. Elet. No 4 (1973)
213
Lyapun0vmeth0d J Machowski Institute of Power Systems, Technical University of Warsaw, 00-662, Warsaw, Koszykowa 75, Poland
In the present paper, methods o f finding the criterial value fi)r the Lyapunov function are considered. A new method is described which efiminates generation ofall possible starting points. The criterial starting point has been determined by findhlg the minimum cut set for a certain network representing transfer synchronizing powers o f generators. The proposed method is faster than techniques published so far.
I. Introduction Although numerous publications have dealt with the application of the direct Lyapunov method to power system transient stability studies, the practical use of this method is far from satisfactory. One of the main drawbacks is the long time necessary to determine the criterial value Vc for the Lyapunov function referring to the transient stability region. When the Lyapunov function in the form of total energy is applied, the criterial value corresponds to the least value of energy necessary to divide the system into two generator groups operating in an asynchronous way. So the critedal value is equal to the value of the Lyapunov function in the lowest saddle point, called the criterial point. Methods of finding the V~. value can be divided into four groups. i. In the starting-point method, the criterial point is found by: a. generating the ith starting point, near to which an equilibrium point may be located b. determining the coordinates of the equilibrium point by solving the corresponding nonlinear equations using, for example, the Newtonian method c. calculating the Vi value for the determined equilibrium point and returning to stage a until all the possible starting points are over d. selecting I7,. = rain Vi. This procedure requires (2 n I 1)-fold solutions of (n - 1) nonlinear equations, and the same number of calculations Received: 9 July 1980
Vol 2 No 4 October 1980
referring to the value of the Lyapunov function. These operations are very time-consuming, and for any power system of reasonable size the procedure becomes computationally prohibitive. A simplified version of the method proposed in Reference 1 consists in displacing stage b from the second position to the last, in order to eliminate the multiple solution of the nonlinear equations and thus make all the computations nmch faster. If the reference generator is selected properly (e.g. the generator remarkable for its highest self-synchronizing power is chosen) the results obtained are usually true. However, in spite of this simplification, the computing time is too long, and still far from satisfactory for practical purposes. Some other simplified versions of this method, proposed in References 2 and 3, will be discussed in Section IV. ii. The method proposed in Reference 4 consists of two stages. During the first stage, the operation should begin from a stable point and move in small steps from one point on an equiscalar surface to the next point having a greater V-function value. In this case, every subsequent point may be determined by minimization of a certain function illustrating the path of least ascent. All these calculations are ceased after an attraction area of the saddle point has been reached. At this moment, the second stage of the operation should be started. It consists in strict determination of the saddle point coordinates by solution o f ( n 1) nonlinear equations. iii. The criterial value Vc is formulated as the following constrained optimization problem: find the minimum of tile V-function on the surface ~ = 0, where 4/is a certain function obtained from the Rauth-Hurwitz condition s or simply the directional derivative of the Lyapunov flmction. The authors of Reference 6 applied the Monte Carlo method to solve the above task. iv. The fourth method differs from the previous three in having a quite different approach. It consists in some particular decomposition of the system and the application of the vector Lyapunov function 7. The method proposed at a later stage in the present paper
0142-0615/040207-07 $02.00 © 1980 IPC Business Press
207
belongs, in general, to the first group. However, the simplifications made allow for relatively fast determination of the criterial point, without generation of all the possible starting points.
II. Conservative model of power system For transient stability studies using the direct Lyapunov method, the generators are represented by constant voltage E ' behind reactances X~ as well as by the equations of generator rotor inovement which can be written as follows: d2~i
- 2 = Pmi Mi -dt
Pi
at i = 1, . . . , n
I I. 1 Application o f coherent generator aggregation In the majority of large electrical power systems, after shortcircuit localization, groups o f coherent generators can be found. Each of them can be replaced by an equivalent generator, by means o f aggregation. The equivalent generator is connected to the system by certain fiction lines. The values of conductance for these lines depend to a great extent on the differences in the angles 6 i of the aggregated generators. In the case of there being great differences in these angles and any reasonable aggregation method being applied, the conductance values are too high to be neglected. Such conductances should be inw~lved in the Vp function.
( 1)
where M i is the coefficient o f inertia, 8i is the generator rotor angle, Prni is the mechanical power (constant), and Pi is the electrical power. Load nodes can be eliminated by assuming that reactive powers o f loads can be replaced by susceptances and that active powers of loads can be replaced by constant injection currents. If this assumption is taken into account, then, after elimination of load nodes, the values of transfer conductances are so small that they can be neglected. Thus, the following formula for electrical power can be derived:
Later, it will be proved that when the aggregation is made with the aid of the Zhukov method (see Appendix 1) the conductances introduced by aggregation are antisymmetrical and form potential moments only, so that these momellts can be exactly taken into account in the calculation of t5,. The Zhukov aggregation is remarkable for one of its properties: if the admittance nodal matrix Yc is of the . . . . *T *T antlhermltlan t y Pe , l • e . YRR =. Y RR, Y~4~I * r ..' , = . Y I ~ M , YmR . . . . YRM, then the matrix Yg obtamed alter aggregation is of the antihermitian type, too. It appears so because"the matrix YRR does not change and y~er
=
(yRM e)*T
= e
* r ,t
,,,7' IR,M
) = e
*'r YMR = YeR
The above is true for any number of aggregated generator groups.
/l
Pi = Psi + Z bi] sinSi/
i = 1. . . . . n
(2) It should be noted that the matrix YG obtained after elimination of the load nodes and resulting negligence of the conductance value is of the antihermitian type (only symmetrical imaginary elements exist).
where Psi is the power o f equivalent loads resulting from elimination o f load nodes - so called self-power and
bi] = BijEiE j The mathematical model in equations (1) and (2) is conservative, as the condition OPi/38 i = aPl/O6i referring to the existence o f potential Vp is valid. In this model, the potential is defined as follows: (Pmi Pi) = aVp/O6i and hence after integrating equation (1), the Vp function may be found by the following formula
Thus, it has been stated that after aggregation, the nodal matrix to describe the equivalent system is of the antihermitian type, i.e. the following expressions are true for every fiction line:
=
rT,
Or Vp = -
~ (Pmi i=1
~..
n
- Psi)(6i
-- 8 s )
Bij = Bii
1
y
%(cosa;j - cosS,';)
and (3)
i=l ]=i+1
where 8is ~s • the coordinate o f the postfault equilibrium point. Each stationary point of the Vp function is considered as an equilibrium point for the system of equation (1) and so the equilibrium points are to be classified (according to the shape of the potential field Vp) into the minimum point, saddle points, and maximum points. The saddle point is classified as the point of maximum type with reference to the k variables and as the point o f minimum type with reference to the remaining variables. Such a point will be called the unstable equilibrium point of the k type.
208
Gi/ = _ Gii The electrical power for the equivalent system can be determined by using the following formula:
Pi = Psi + ~, (bij singq + gij cos6ij)
{4}
at g i / = GijEiEI. As the conductances introduced within the aggregation are antisymmetrical, the condition OPi/~8/= OP//OgJi referring to the existence of potential Vp is valid, and after integra-
Electrical Power & Energy Systems
tion of equation (1) the Vp function is as follows: n
Vp =
~
X (Pmi - Psi)(Si - 6~) i
1
n -- 1
X i~l]
x [bd(cos6ii - cos6Si ) -gi/(sin6i/
i-1
sin6Si)]
(5)
II I-. Proposed method If the trajectory of generator rotor movement crosses the crest of the type k saddle point, k generators lose their synchronism in relation to the remaining (n k) generators. In this respect, searching for the criterial saddle point is equivalent to looking for such a division of generators (into two groups) that will require minimum energy. Function (3) or 15) determining the energy in question may be expanded (about stable point 6 s) into a Taylor series. Taking into consideration that the gp function and its gradient in point 6 s are equal to zero, it appears that
where A8 : 6 8 s and H is the square, symmetrical Hessian matrix of potential energy H = [OZVp/aSi OOj] and, because 3 Vp/O8 i = -Pi, it may be considered as a Jacobian matrix for equation system (2) H = [3Pi/36il. Approximation (6) is true in the neighbourhood of the stable point 6 s only and should not be used for an approximate calculation of the Vp value for unstable equilibrium points located at some distance from the stable point 8 s. However, after performing calculations for a number of exemplary test models of electric power systems, some valuable properties of approximation (6) were noted. It shows the least value at the same unstable equilibrium points for whidl the value of function (3) or (5) is the least, too. Apart from that, in most cases, approximation (6) permits proper classification of all the unstable equilibrium points, according to their increasing potential energy (3) or (5). The application of these properties will now be discussed. I I I. 1 Unstable e q u i l i b r i u m p o i n t s k = 1 In this case, the following vector can be assumed: A8 ~ [0 . . . . .
O, A6i, 0 . . . . .
111.2 Unstable equilibrium points k > 1 The situation in which k > 1, i.e. when the group of generators {L} loses its synchronism with respect to the remaining generators {R}, may be considered as the case when the synchronism is lost by one generator e which replaces the whole group {L}. The self-synchronizing power for such a generator results from principles of aggregation and is equal to the sum of the transfer synchronizing powers between the group of generators {L} and the group {R} (see Appendix 2)
Hee:HL/~:
Y. iCL
E
¢8)
( H~i)
]~R
The Hee value will be called the equivalent synchronizing power. Assuming that a set {L} consists of one element only, the fornmla in equation (8) becomes simplified to the formula given in equation (7). The group of generators {Lc} for which the equivalent synchronizing power (8) does not exceed the self-synchronizing power of each generator belonging to the group will be called the compact group. The term compact group is justified by the fact that in order to meet its definition the transfer synchronizing powers inside the group are much higher than the transfer synchronizing powers between the generators of this group and the remaining generators of the system. It is also worth noting that each generator satisfies the above definition and forms a single-element (trivial) compact group. In conclusion, it should be stated that determination of the unstable equilibrium point of the k > 1 type, showing the least Vp value (or, more exactly, a suitable division of generators) requires the establishment of a compact group showing the least equivalent synchronizing power (8). Such a group will be called the criterial compact group, II 1.3 A l g o r i t h m s f o r research o f criterial c o m p a c t group
Let the generator nodes of the examined system be the vertices {W} of the graph G with the lines {U}. Let the real numbers ( Hij ) describe the lines of the graph, and the real n u m b e r s H u describe the vertices. The subset {P} of set{U} can be defined as follows:
{P} = [
O] 7"
where
Approximation (6) has then been simplified to the form : i A , ~ 2T/ Vp ~ 2 ~ i * - i i .
Thus, the energy necessary to cause the loss of synchronism by generator is approximately proportional to its selfsynchronizing power Hii where
{Wl} U {W2} = {W}
{~]1} ("1 {~¥2} = ~)
On removal of the lines representing cut set {P} there is no pail1 from {W1} to {W2} and in this way the cut set {P} may be identified by the pair of subsets {Wl} and {W2}. The real function for the cut set {P}
Hii = ~ ( Hi i)
(7)
Thus. the point of the k = 1 type showing the least Vp value is easily determined by means of finding the generator showing the least self-synchronizing power (see Section Ill.5.1).
Vol 2 No 4 October 1980
H(P)= ~ ~ (-Hxy) x~iw,} y~iw2)
(9)
may be defined as the maximum flow of the cut set. The cut set for which the value of equation (9) is the least possible one is considered as the minimal cut set.
209
The above definition (9) is illustrated by Figures 1 and 2. Formulas (8) and (9) are identical, and hence the determination of the criterial compact group is equivalent to obtaining the minimal cut sets between every pair of vertices. Such a task can be solved very quickly even for large networks through use of the Gomory and Hu theorem s, by nreans of the (n - l )-fold application of the labelling procedure 9, m or by use of the Korzanov procedure n.
I
2
4
3
5
6
7
Anita
~"pi
1 1961<1 317 2693 3 566 1581 366 157}I i ; _ l ~ ~ 2 11315 725011939 2 531 i 122 241 112 [ (L(~90irI4095! 3 267t t933 9 8 9 ~ 3 083 1368 580 270ii oioi211i7328 i 4 3562 2530 307314 54814663 495 23n 0.094!28 198 5 11 577 1 121 1366 4 661 9046[ 220 103{ f/(:igg~-1762~ 6 i 334 240 580 494 219 3 6 5 4 1 1 7 > n.030~7t80! 7 i 155 111 269 ,02,787 =6, t t 704i
I 11.4 Calculation of coordinates for criterial point On determination of the minimal cut set, the criterial starting point is established by providing the angles of generators { W=} with the value ~is, and the angles of generatars{W1} with the value (-+lr --6s). The plus sign is given for zr in the case y~
(e~; - P~) > o
',
i~{w,} i.e. when the generators of the group {W~} export the power outside. The minus sign is used if the generators of the group { W1} import the power so that
}2
O % ~ - G ; ) < o.
iE{W,}'
2
(
%
F i g u r e 3 Example l: system containing single-element compact groups only
The above arrangements of plus and minus signs have been proved and described in detail in Reference 12. The coordinates of the criterial point are obtained by solving equation (1) in the neighbourbood of the starting point while d25i/dt 22 O. For this the Newton-Rap!tson or Davidon-Powell-Fletcher method can be used.
Figure 1 Example of cut set w i t h k = card. {Wt} = 1; maxim u m f l o w for cut set is H(P) = - ( H t 2 + H t 3 + H t 4 + H l s )
"~v
111.5 Examples The data for the examples is given in Figures 3, 4 and 5 where the numbers of the rows and columns correspond to the numbers of particular generators. The values given in the upper triangular matrix refer to the parameters bii included in equation (2). The values given in the lower triangular matrix together with the diagonal elements refer to transfer synchronizing powers as well as to self-synchronizing powers. The vector on the right contains the angles
For better illustration of the proposed method, the dendrites have also been included in all the examples; the dendrites, which are the flow-equivalent trees, enable fast visual detennination of the maximal flows for all possible cut sets. 111.5.1 Example 1 This system contains the single-element compact groups only. Thus, the criterial point is o f tire k = 1 type and corresponds to the loss of synchronism of generator 7. The value of function (3) calculated for this point is equal to Vc = V,~7 = 4 704. hr this exmnple, the values of the Vp functions for all unstable equilibrium points of k = 1 are also given.
F i g u r e 2 Example of cut set w i t h k = c a r d . { W ~ } = 2 ; m a x i m u m f l o w for cut set is H(P) = -- (Hi3 + HI4 + H i s + H23 + H24 + H25)
210
111.5.2 Example 2 This system contains tile doubleelenrent compact group {6, 7} which is the criterial compact
Electrical Power & Energy Systems
1
2
1 ~1276
3
4
3160 3299
5
6
7
Angles
1732
280
230
0.070
1713
902
148
121
0.072
3 3159
16621110713415
1794
608
473
0.043
4 3291 5 1731 6 279
1709 3 4 1 2 1 3 7 8 6 1 4 6 7 3
389
319
0.000
168 203 4 2 9 6 1 2 6 7 9
0.049 0.005
167 2 6 7 t
0.069
2 1264 5797]_1673
7
230
895 1794 4 6 6 7 147 607 389 120
473
318
94571 204
3979
ally, attention should be paid to the fact that the advantage of the proposed method increases as the number of generators increases. Further, it should be pointed out that white preparing the above program, not all the possibilities of computation time reduction were taken into account and so the times given in Figure 6 can be reduced further, e.g. by the application of the algorithm in Reference 1 1 instead of Reference 13.
IV. Other simplified methods
-.~
">--< ,,,jC-J
~c~ 1 ~ . . .
/ J J J
/
/
J
Figure 4 Example 2: system containing criterial compact group{6, 7}
group. In this case, the criterial point is of the k = 2 type and corresponds to the loss of synchronism of generators 6 and 7. The function (3) value of this point is equal to Vc = 5 864. The minimal point of the k = 1 type corresponds to the loss of synchronism of generator 7 which belongs to the criterial compact group. 111.5.3 Example 3 This system contains (as in Example 2) the criterial compact group {6, 7}, provided that Ve = 5 782. In this case, the minimal point of the k = 1 type at Ve = 6 843 corresponds to the loss of synchronism in generator 2, which does not belong to the criterial compact group. Thus, for this example, the Gupta and El-Abiad assumption is not valid. The values Vp in points of the k = 1 type corresponding to the criterial compact group are Ve6 = 7 756 and Vp7 = 7 207.
In accordance with the observations mentioned in Reference 2, it was noted that in typical systems the criterial starting points are of the k = 1 type. On the basis of these observations, the suggestion has been made of limiting the calculations to 'n' starting points of the k = 1 type. This property is true for the majority of electrical power systems, and more so when considering systems after aggregation of coherent generators. However, some cases can be found where the criterial points are of the k > 1 type in spite of aggregation. This results from the fact that the compact group can be coherent only if the shortcircuit is situated at some distance and involves approximately the same accelerations. If the coherency criteria are not satisfied because of nonuniform accelerations or relative distances 14, ls, the program of automatic reduction leaves the compact group in nonaggregated state, and the criterial point can be of the k > 1 type. The other simplified method is given in Reference 3, and its basis may be explained briefly as follows. Let {Gx} be a set of generators losing their synchronism in minimal 1
2
3
4
6
7
280
230
0.139
685 3 821 l 1 039 1 078 760 147 120 3 1 5 7 1 0 3 6 1 0 4 4 8 1 3 4 1 6 1794 577 472 3 4 3 2 9 2 1077 3 4 1 3 1 3 15614674 388 319 757 1794 4 6 6 8 9 3 2 1 ] 204 168 5 1731 279 147 576 388 203 4 262 {2 689 6 230 119 472 318 168 2669 3 9 7 6 7
0.035
1 9 3 7 4 I 689 3159 3 2 9 9
5
1731
Angles
2
I [
0.112 0.069 0.118 0.064 0.138 //
~//
\\\ _ ////
111.6 Computation time On the basis of the above method, the programme has been prepared for a CDC 3170. In the program, the algorithm in Reference 13 and the standard Newton-Raphson procedure were used. The computations have been performed for many different test systems, and compared with the results of calculations using the exact starting points method and with the results of the Prabhakara and EI-Abiad method (the generator with the highest self-synchronizing power being taken as the reference generator). The results of calculations have confirmed the exactness of the proposed method. Average values of computation times are given in Figure 6. The time requirenrent of the computations for the proposed method is much less than for the remaining two methods. Addition-
Vol 2 No 4 October 1980
// /
/
/
/
Figure5 Example 3: system containlngcriterialcompact group{6, 7}, provided that Vc = 5 782
211
4o 3o
Ribbens-Pavella, M and Lemal B 'Fast determination of stability regions for online transient power-system studies' Proc. IEEE Vol 123 No 7 (1976)
I I
3 Gupta, C L and EI-Abiad, A H 'Determination of the clOSest unstable equilibrium state for Lyapunov methods in transient stability studies' IEEE Trans. Power Appar. & Syst. Vol 95 (September/October 1976)
I
I I I
Tavora, C J and Smith, O J M 'Stability analysis of power systems' IEEE Trans. Power Appar. & Syst. Vol 91 (May/June 1972)
I I
E
E
Weiman, M J 'Regions of permissible deviations in determining synchronous transient stability of power systems by Lyapunov's second method' (in Russian) Elektrotekhnika No I (197 I)
I / I
Strakhov, S V and Weiman, M J 'The present state of
f I
(:z
E
o (.2
development and possibilities for practical application of Lyapunov's second method in determining the transient stability of power systems' (in Russian) Elektrotekhnika No 10 (1977)
/
10
//
// O
5
/ /
10 Number
15 of
f
20
1
f
/
f
Darwish, M and Fantin, J 'Stability analysis of large scale power systems using the comparison principle and vector Lyapunov functions'J. AppL Sci. & Eng. No 3 (1978)
25
generators
Christofides, N Graph t h e o r y - an algorithmic approach
Figure6 Average values of computation times of CDC 3170; - - - - exact method, - . . . . . Prabhakara and EI-Abiad's method, - . . . . proposed method
saddle point o f k = x type. According to file Gupta and E1-Abiad assumption, {Gx_l} should be considered as tile subset of the set {Gx}. This assumption has been supported by examples; however, it is not generally true. This is illustrated in Example 3, where the minimal point of the k = 1 type corresponds to the loss of synchronism of generator 2, and the criteri',d point is determined by generators 6 and 7. The error made in Example 3, while using the method in Reference 3 to determine Vc, amounts to about 1g~/~, and can be much larger if the parameters in Example 3 are changed.
V. Conclusions One o f the main disadvantages o f the application of tile direct Lyapunov method to power system transient stability studies is the lack of methods that would enable relatively fast determination of the criterial value. In the presen! paper, a new technique has been proposed, in which tire time-consuming task of finding the criterial point has been replaced by determination of the minimal cut set of a cerlain network. The proposed method has been tried on many test systems and tire results obtained were considered to be correct. The computation times are then significantly shorter than with previously used methods.
References 1 Prabhakara, F S and EI-Abiad, A H 'A simplified
VI.
determination of transient stability regions for Lyapunov methods' IEEE Trans. Power Appar. & Syst. Vol PAS-94 (March/April 1975)
212
Academic Press, USA and UK (1975) Ford, L R and Fulkerson, D R Flows in networks Princeton University Press, USA (1962) 10 Edmonds, J and Karp, M 'Theoretical improvement in algorithmic efficiency for network flow problems' J. Assoc. Comput. Mach. No 19 (1972) 11 Hamacher, H 'Numerical investigation on maximal flow algorithm of Korzanov' Comput. No 22 (1979) 12 Maciejewski, Z 'Direct method of power system stability evaluation' Prace Instytutu Energetyki Zeszyt 6 (1978) 13 Bayer, G 'Algorithm 324, Maxflow' Commun. ACM No 11 (1968) 14
Lee, S T g and Schweppe, F C 'Distance measures and coherency recognition for transient stability equivalents IEEE Trans. Power Appar. & Syst. Vol PAS-92 (September/October 1973)
15 Machowski, J Simplified models for power system stability studies PhD Thesis, Technical University of Warsaw, Poland (1978)
VII.
Bibliography
Bayer, G 'Remark on algorithm 324' Commun. ACM No 16 (1973) Dinic, E A 'Algorithm for solution of a problem of maximal
Electrical Power & Energy Systems
flow in a network with power estimation' Sov. Math. DokL No 11 (1970)
where
OPi
Hq-
Zhukov, L A 'Construction of power system electromechanical equivalents' Izv. Akacl. Nauk SSSR Energ. & Transp. No 2 (1965)
3Pi aii = _
transfer synchronizing power
self-synchronizing power
Appendix 1. Aggregation of coherent generators Prior to aggregation of the selected generator group {M}, the system may be described by nodal equations in the following form :
Linearization of equation (11)results in the following equation :
[APn] = [ H R R APe 3 tHen LYMR
IM
YMM 3
or I(;=YGEG
or VeR
Yee J
lg = YgEg
(11)
~"
where Yne = YRM e, Yen = e*TyMR, Yee = e*TyMM e, and e = Ee~EM;Ee is any voltage at the centre of the angle (the above values are complex). Attention should be paid to the fact that this aggregation is equivalent to the aggregation considered by Mello, Podmore and Stanton 16, assuming that ideal equipotentializing transformers taken into account by the authors have a complex transformation ratio equal to e i = Ei/E e.
Appendix 2. Synchronizing power of equivalent generator Taking into account that E i = const., the following linearization of the equation (t0) can be made:
[a e R ] _ [ . . . APM3
HRM1
ApG=HGA6 G
L H M R HMMJ LAfMJ or (12)
Vol 2 No 4 October 1980
APg=HgAgg
(13) and it can be proved that
Assuming that the power of the equivalent generator is equal to the sum of aggregated generator powers, and assuming that the aggregation does not involve changes in currents and voltages of nodes {R}, it appears that =
or
(10)
where {R} are the remaining nodes of the system.
Ie
HRe] [ A f R ] Hee tASe J
HRe = HRM 1M, HeR = 1THMR, Hee = 1THMM 1M This means that the Zhukov aggregation may be considered as a nonlinear generalization of the linear Marconato-Marzio aggregation 17. As the sum of elements in the row of matrix H G is equal to zero, the synchronizing power nee may be represented in the following form
gee = 1MHMM1M = -ffIH
R1R =
Y. Y. ,~{M} j~{~}
-g,p
The symbols 1M, 1R represent vectors consisting of ones only, with cardinality {M} and {R}, respectively.
16 Mello, R W, Podmore, R and Stanton, K N 'Coherency based dynamic equivalents: Applications in transient stability studies' Proc. PICA Conf. (1975) 17 Marcenato, R and Marzio, L 'A dynamic model of Italian power system based upon area equivalents' Energ. Elet. No 4 (1973)
213