New direct lyapunov-type method for studying synchronization problems

Automatica, Vol. 13, pp. 139-15l. Pergamon Press, 1977. Printed in Great Britain

New Direct Lyapunov-Type Method for Studying Synchronization Problems* E R I K J. N O L D U S t

When studying synchronization which occurs in synchronous machines or in phase lock circuits, special stability problems, requiring a new method of solution, arise if the phenomenon of synchronization after skipping a number oj cycles is considered. Key Word Index Stability; nonlinear systems; power system control; synchros; Lyapunov methods. Summary--A new approach of the Lyapunov type is presented for the stability analysis of a class of time invariant, nonlinear dynamical systems, which are characterized by a periodic nonlinearity, and the existence of an infinite number of stationary states. In particular, the method can be applied for studying transient stability problems in power systems, and the synchronization of phase lock circuits. Frequency criteria are developed for the convergence of all solutions to an equilibrium state, and a procedure, based on these criteria, is outlined for constructing stability regions for a collection of neighboring equilibria. When applied to power--or phase lock systems, these regions can be interpreted as synchronization regions within an arbitrary, given number of cycles. Also, some new results are obtained with respect to synchronization conditions within the first cycle.

The system (1) has considerable practical relevance. Under the usual simplifying a s s u m p t i o n s [ l ] , a power system consisting of a single synchronous machine, connected to an infinite bus, can be represented by an equation

1. INTRODUCTION

D26+aDiJ+ f(6)+c=O

p(s) and q(s) are polynomials in s, with degree p(s) > degree q (s), and q ( 0 ) ~ 0 . f ( . ) is a time invariant, nonlinear, periodic function with period

f(u+z)=f(u),

IN THE PAST decade L y a p u n o v techniques have been extensively applied to the study of synchronization problems. In particular, they have been used in the area of transient stability studies in p o w e r s y s t e m s [ l ] , and to investigate the synchronization of phase lock circuits. Interesting results have been obtained by representing such systems in the form of a nonlinear feedback loop, and then applying systematic procedures, available for constructing L y a p u n o v functions for such loops[2]. Because of the existence of an infinite n u m b e r of equilibrium states, the stability p r o b l e m arising in this context is not one of global stability, but one of evaluating the d o m a i n of attraction of a single equilibrium state. The i m p o r t a n c e of developing accurate L y a p u n o v m e t h o d s for synchronization studies, stems from the fact that these methods seem to constitute the only available, general analytical tool, for estimating finite stability regions. The present paper is concerned with the stability analysis of a dynamical system of the form

Dp(D)z+f[q(D)z]=O,

Da=d/dt

allu.

(2)

(3)

6 is the angle, in electrical degrees, between the rotor shaft and a shaft running at synchronous speed, while a, f ( . ) and c represent the normalized d a m p i n g coefficient, synchronous torque characteristic and speed governor output. F o r the unregulated machine, with c - 0 , (3) clearly belongs to the class (1). If governor action is represented by a transfer function c(s)/s6(s)= m(s)/n(s), then, defining 6 a__n (D)z, (3) becomes

D[Dn(D) + an(D) + m(D)]z + f[n(D)z] = O, which again has the form (1). The dynamic equation of phase lock circuits[3], [4], can also be written as (1). Then

H o(s) a=q(s)/p(s) represents the transfer function of the feedback filter, while f ( . ) is the phase c o m p a r a t o r characteristic. In the classical techniques for deriving stability regions for nonlinear systems[5],[6], and in a n u m b e r of papers dealing with their application to transient stability and synchronization problems[7-9], the d o m a i n of attraction of an equilibrium state is generally b o u n d e d by a surface of the form V ( x ) = V,,, where V(x) is a suitable L y a p u n o v function. This m a y be the L y a p u n o v function constructed for the p r o o f of the well k n o w n P o p o v criterion, or some of its generalizations. The scalar V,, is obtained as follows : Let { - a < u < b} be

(1

*Received 10 October 1975; revised 24 August 1976. The original version of this paper was not presented at any IFAC meeting. It was recommended for publication in revised form by associate editor I. Landau. tPresent address: Automatic Control Laboratory, University of Ghent, Grote Steenweg Noord 12, B-9710 Zwijnaarde, Belgium. 139

140

E . J . NOLDUS

the range of the argument u, for which the nonlinearity f ( u ) satisfies the sector inequality derived from the considered stability criterion. Then V,~ is selected as the maximum value for which the corresponding stability domain lies inside the region { - a < q ( D ) z < b } . Often substantial improvements can be obtained, using Willems' method of open Lyapunov surfaces[10]. Here the stability region is bounded by parts of the hyperplanes q(D)z =-a and q(D)z=b, and by a surface V(x) =constant. An obvious disadvantage of these methods is that the region { - a < q (D)z < b} cannot possibly contain other equilibria than the origin. As a result these methods are only suitable for studying synchronization within the first cycle. In the present paper an alternative approach is developed towards the stability problem of the system (1). The method proceeds in basically two distinct steps: First, sufficient conditions are derived, which guarantee that the system is nonoscillatory. By 'nonoscillatory' is meant that, in state space, every bounded trajectory converges to an equilibrium state, as time increases. This property implies that sustained, bounded oscillatory motions cannot exist in the system. Next, suppose that a closed, bounded region f~(x) in state space can be found, such that any trajectory x(t, Xo), starting at t = 0 at some point Xoef~(x ), remains in {)(x) for all subsequent time. Then, assuming that the system is nonoscillatory, x (t, x 0) must converge to an equilibrium state in ~(x). If f~(x) contains a single equilibrium, say the origin, then f~(x) is a stability region for this equilibrium. If f~(x) contains several equilibria, then any trajectory starting in f~(x) converges to one of them. In other words, ~(x) is a stability region for the collection of all equilibria contained in it. The concept proposed here has two specific advantages: (a) The stability region f~(x) may be constructed completely independent from the Lyapunov function that was used for proving the system's nonoscillatory behavior. In particular, analytically 'simple' stability regions may be looked for, such as regions bounded by straight line segments, in the case of second-order systems. (b) Due to the periodicity of the function j{. ), there exists a spatial periodicity in the pattern of the system's trajectories in state space. Partially as a result of this periodicity, it is equally simple to find a stability region for a collection of neighboring equilibria, as a stability region for a single equilibrium. Hence regions of initial states may be found, from which the system synchronizes, eventually after skipping m cycles or less. In other words, synchronization regions within m cycles can be obtained. The structure of the paper is as follows. In section 2 some frequency criteria are reviewed for the convergence of all bounded solutions of equation (1)

to an equilibrium state. These criteria have been derived in previous papers by the author[11], [12] for certain classes of systems which are somewhat more general than (1). They have been repeated here for easy reference, and for use in the examples. Sections 3 and 4 present criteria, also of the frequency type, which ensure synchronization from all initial states, eventually after skipping a finite number of cycles. Procedures are outlined for constructing synchronization regions within the first cycle, and within the mth cycle. First. secondorder systems are discussed, for which stability regions can conveniently be visualized in the state plane. Subsequently the extension is presented to systems of arbitrary order, and the merits and shortcomings of the method are pointed out. In the final sections some examples are worked, dealing with power systems and phase lock circuits. 2. C R I T E R I A F O R N O N O S C I L L A T O R Y

BEHAVIOR

Let 2 = Ax - bJ(c'x)

t4)

be a vector representation of the system (1), with e ' ( s l - A ) - l b = H ( s ) ~ ( l ~ s ) H ~ { s ) . The state vector x is related to z by x = A d j ( D I - A ) , hz. The basic principle for proving the convergence of all bounded solutions to an equilibrium state, proceeds as follows: Let V(x) be a scalar function of x, such that, along the solutions of equation (4), l?(x) does never change sign. Furthermore, assume that the largest invariant set where l?(xl= 0, consists of the points where 2 = 0 , i.e. the collection of the system's stationary states. Then, according to LaSalle[13], every bounded trajectory approaches an equilibrium. The following theorem is based on this principle Theorem l. Assume that either I

Im HUgo)>0

all real ~o

(5.1)

- I m H(jco)<0

all real co

(5.2)

SO or

1

(O

Then all bounded solutions of the system (4} approach an equilibrium state. A short proof is outlined in Appendix A, mainly for easy reference, since the associated Lyapunov function will be needed in the examples. The conditions (5) may be relaxed, at the expense of imposing some restrictions on the nonlinearity f ( . ). The following theorem can be found in [12]: Theorem 2. Assume that (1) For some scalars k 1 and k z, the slope of the nonlinearity satisfies kl < f~l~(u)<=k 2

allu.

New direct Lyapunov-type method for studying synchronization problems ^

141

ADz

~(×,)

WlXl = -~

FIG. 1. Synchronization region within the first cycle for a second-order system.

(2) Scalars e and fl, fl > O, can be found such that R e - - (~+je)fl) jco

x (J)J

> 0 all real e). (6)

Then all bounded solutions of the system (4) approach an equilibrium state. Theorem 1 is a special case of theorem 2, as can easily be shown by letting fl = 0 and ~ = + 1. 3. S Y N C H R O N I Z A T I O N REGIONS FOR S E C O N D - O R D E R SYSTEMS

It will be helpful to consider first the case where the system (1) is of second order. The nonlinearity f ( . ) satisfies (2), such that f(u) = 0 for u = kz and for u='c 1 + k r , 0
Dp(D)y+ f~[q(D)y] =0, where

(7)

fl(U)a=f(unt-'Cl). Let (7) have the vector

representation

:¢x = Ax,

-- bf l (C'Xl ),

(8)

where xl ~=Adj(DI-A).by. Now assume that an A p indefinite* quadratic form P(xl)=xlPxl,P=P', can be found such that, along the solutions of (8), and for a scalar r > 0,

P(xl)+2rP(xl)
all xt #0.

(9)

In Fig. 1 the region O(xl)a={xl;P(xl) 0 , then by (9), P ( x i ) < 0 , and x t cannot be an equilibrium. Furthermore, since P(xl ) = 0 implies P ( x i ) < 0 for xl:p0, the trajectories point to the interior of@(x 1) on (1 and (2. We now invoke the periodicity of the function f(. ), and the resulting spatial periodicity in the pattern of the trajectories in state space. The substitution z = y + z l - r again transforms (1) in (7), such that, proceeding as before, a translated region ~ ( x t ) is *By indefinite is meant neither positive (or negative) definite, nor semidefinite.

142

E . J . NOLDUS

obtained, with identical properties as (I)(x 1 ). (I)r (X I ) is generated by translating (I)(Xl) in the direction of the negative y-axis (or z-axis) over a distance r. It follows that on the boundary 8O(x), all trajectories point to the interior of O(x). Hence any trajectory starting in ~(x) converges to an equilibrium in O(x), i.e. to either I ~, So or I o. If a stability region for so is desired, the unstable equilibria I_ ~ and I o must be cut out of O(x). This is done by introducing an additional set A A )2 __ ~,2 LIJ(x1 } = {.X'I ; W (x 1 ) = (I4/x 1 , ~0},

having the properties that (a) The straight line w'x~ = 0 lies outside the sector (l)(x ~). (b) If W(x~)=O, P l X l ) < 0 and [eI is sufficiently small, then W (x ~) > 0. In other words, on that part of the boundary 8~(x~ ) which lies in (I)(Xl) the trajectories point to the interior of ~(x~). Hence no trajectory can leave (l)(x~) c~ Ud(x~). Since ]eI may be chosen arbitrarily small, and since the argument can be repeated for the equilibrium I_ ~, it follows that any trajectory starting in ~(x) converges to S o, except the equilibrium solutions I ~ and I o. With this restriction,* f~(x) is a stability region for So. The following theorem yields a frequency criterion, which guarantees the existence of the quadratic forms x]Px~ and (w'x~) z with the desired properties. Theorem 3. Assume that for the second-order system (7) with vector representation (8) (1)

k3tt 2 <=uJ1(U)~ k4u 2

(10)

for all u and for some scalars k 3 and k 4 (2) The linearized equation 2~ = A t x ~, A~ ~=A -/~l~(O)bc',+ has one characteristic value in the open right half plane, and the other one in [Re s < -d<0',. (3) For a scalar r, 0 < r < d , [1 + k,,H(jo) - r)] Re L ~ + ~ / - i ~ - _ r)j > 0 all real ~).

(11)

Then: (1) There exists an indefinite quadratic form P(xl)A=x'~Px~, P=P', such that (9) holds. (2) There exists a quadratic function W(x~) a= (w,xl)2 _ c 2 such that P(xl)l,,x~ 0 > 0

for

x~v~0,

*A similar restriction applies to other stability regions, discussed below. fFhis is the system's equation, when linearized at the equilibrium I,~.

and W ( x t ) > 0 Jor W(x1)=0, P(xl)<-'O and ici sufficiently small. The proof will be deferred to section 4, where the extension to nth order systems is discussed. Let us first examine the implications of the present result. First, as already explained, the conditions of theorem 3 guarantee the existence of the stability region O(x) for the origin S o. ~(x) can be determined by calculating the symmetric matrix P, which, as shown in Appendix B, satisfies a quadratic equation of the Riccati type. The resulting stability region ~(x ) has the advantage of simplicity, since it is bounded by straight line segments. This may be useful for systems analysis at the design level, where simple, albeit conservative expressions of the stability boundary may be desirable. However, further exploiting the spatial periodicity of the system's trajectory pattern, substantially more information can be drawn from the obtained results. Indeed, for any integer hi, let ~,~(xl ) be the sector, derived from ~(x~ ) by a translation in the direction of the z-axis, over a distance ni~ (in the negative sense for n~>0, and in the positive sense for n~<0). Then any region of the form O,~(x: ) c~(1),~(x~), n~, n2 integers, has the following properties" (a) No trajectory starting in this region can ever leave it, as time increases. (b) The region contains a bounded, closed subset ~.i(x), delimited by straight line segments, and separated from the other subsets by two unstable equilibrium points, located at the vertices of the sectors cb,,~(x t ) and (l),~(Xl). The subset fl~,j(x) is a stability region for the collection of all equilibria contained in it. In this collection, the two unstable equilibria on the boundary of f~.j(x) need not be included, for the reasons explained before. For example, in Fig. 2, {I-~Q411Q 41 1} is a stability region for the set of equilibria [So, I o. S1}. Since I 0 is locally unstable, this means practically that, starting from any initial state in this region, the system will settle down in either S o or S~. Indeed, ifa trajectory in the stability region approaches I0, then the slightest disturbance transfers the system state from this trajectory to a different one, which converges to either S o or S~. Hence a synchronization region within two cycles is obtained. Similarly, {I_ l Q612Q_6 I_ 1} is a stability region for the set of equilibria {So, I0, Sx,ll, $2], i.e. a synchronization region within three cycles, etc. Furthermore, for any point x 0 of the state plane, a pair of sectors (I),,~(x x ) and q),~(x~) can be found, such that the corresponding synchronization region f~,j(x) contains Xo. It follows that, starting from any initial state, the system will eventually synchronize, after skipping a finite number of cycles. For any Xo, a collection of equilibria can be found, in one of which the corresponding trajectory will settle down. This means that an upper bound can be found for the

143

New direct Lyapunov-type method for studying synchronization problems

Q6

Oz

~

h Z

FIG. 2. Synchronization regions within m cycles for a second-order system.

number of cycles the system will skip before synchronization. Finally, the z-range from which the system synchronizes within the first cycle, in the case where Dz = 0, consists of the entire segment of the z-axis, between the two unstable equilibria, neighboring the origin. This is the largest segment that can possibly exist. In transient phenomena in power systems, initial states of this type occur, for example, in the case of step changes in mechanical power input. The property implies that, for arbitrary such step changes, within the static stability limits, the system always synchronizes within the first cycle. For phase lock circuits the implication is that, for arbitrary step changes in the frequency of the input signal, again within static stability limits, acquisition occurs within the first cycle. One condition for a system to possess these stability properties, is that the relation (11) can be satisfied for a scalar r within the range 0 < r < d. The condition (11) has in general the physical significance, that the system must be sufficiently damped. This is the main drawback of the method: If (11) cannot be satisfied, no synchronization regions can be obtained. However, for weakly damped systems, it may still be possible to construct synchronization regions by combining the sectors ~(xl) with a suitable Lyapunov function. Some results of this type will be derived in the examples.

4. SYNCHRONIZATION IN HIGHERORDER SYSTEMS

Systems of arbitrary order can be studied using the following theorem, which generalizes theorem 3. Theorem 4. Assume that for the nth order system (7) with vector equation (8)

(1)

(12)

k3u 2 ~ ttfl (u) =< kau 2

for all u and for some scalars k 3 and k4. (2) The linearized equation ~1=Alxl, Ala=A -f(l)(O)bc', has one characteristic value in the open right half plane, and (n - 1) characteristic values in {Res< - d < 0 } . (3) For a scalar r,0< r
+ k3H(jc°-

>0

all real o.

(13)

Then: (1) There exists an indefinite quadratic form P(Xl ) a=x'lPx 1, P = P', such that

[a(xl)+2rP(xl)
all x1¢0.

(14)

(2) There exists a quadratic function W(x~) A(W'X1)2--~ 2, such that P(x~)lW,Xl=o is positive definite in the remaining ( n - 1 ) independent state variables, (15) and, fl/(xl)>O for W(x,)=O, P(x,)<__O and Is[ sufficiently small. (16) A proof is given in Appendix B. Theorem 3 is a s__~ecial case of the present result. The set ~(x 1) {x 1;P(x I)__<0} represents a quadratic cone in state space, which is the natural generalization of the sector O(xl) in the second-order case. By (15), • (x 1) can be written as (~(xl)=fx1,P(Xl)=

- - a 2(w,xl

"~i z/Z~ 0},

"[-1"= 1

)2 a:~0

(17)

144

E . J . NOLDUS

where w'x~ and z~,i= 1..... n - 1 , are linear independent combinations in the components of xp It follows that qb(x~ ) can be partitioned in the subsets

{ (I)

(X1)=

nl

X1, •

-

a 2 (w'x~)2

AV E Z250

}

W'xI
i=1

(1)+(X1) =

X l ; - - a 2 ( w ' x l ) 2AV E

2/250,

w'xl>0

(a) The condition that equation 2 : - A l X : has exactly one unstable characteristic value, somewhat restricts the applicability of the result. This condition is necessary however, in order that the region fl(x) is bounded. If A 1 has n i > 1 unstable characteristic values, then the corresponding quadratic form P(x 1) takes the form:

,

nI

i=1

and the unstable equilibrium state xl = 0. Furthermore, by (14) and (16), no trajectory can leave the set

(1)-(x1)O{Xl;W'Xl ~__-g.},

(18.1)

O+(xl) m {xl;w'xl >~),

(18.2)

or the set

for e > 0 sufficiently small. Now proceeding as in section 3, the set (18.2) is translated in the direction of the negative y-axis, over a distance ~, yielding

O+(xl)~{xl;w'xl+w(O)z~e},

(18.3)

where we have defined w'x I ~=w(D)y. Then (18.3) is intersected with (18.1), to obtain a region t2~(x) with the following properties: (a) fl~(x) is closed and bounded. To see the boundedness of tl~(x), note that xeric(x) implies w'x~ < - ~ , and w'x I + w ( 0 ) : > e . Hence (W'X1) 2 is bounded, and since

E z/2

n n1

- E w/2 + E z/2,

):,

i-i

Ix~l is bounded. (b) Any trajectory starting in tl~(x) remains m it for all subsequent time. (c) fl~(x) contains no equilibrium except the origin

So={z=Dz . . . . =D" lz=()l. Indeed, all equilibria are located on the z-axis, and the segment of the z-axis in t2~(x) {T 1 - - T ~'- 0(~,,)---~ Z 5 Z'l - - 0 ( ~ 3 ) } ,

i=1

i-1

and no bounded region tl(x) can be constructed from it. (b) For third-order systems, any intersection of the stability region fl(x) with a plane passing through the z-axis, is again bounded by straight line segments. Intersections with other planes are bounded by segments of quadratic curves. (c) All other remarks following theorem 3 can immediately be extended to the nth order case. In particular, the conditions of theorem 4 ensure boundedness of all solutions, and, together with a criterion for nonoscillatory behavior, yield sufficient conditions for synchronization from all initial states. Synchronization regions with rn cycles, for a given m, can be found by calculating the quadratic form x'lPXl, using equations (49.1, 49.2) in Appendix B. (d) For given system parameters, and for a given r within the permissible range, such that (13) holds, the system (49.1, 49.2) in Appendix B has no unique solution. Any solution yields a different stability region fl(x). In addition, for a given set of system parameters, (13) may be satisfied for a range of permissible r-values, r I < r < r 2, and for any r in this range, corresponding stability regions fl(x) may be found. The union of these regions is again a stability region, possibly larger than any fl(x) taken separately. Presently no general computer algorithms seem to be available for detecting all solutions of a system of equations such as (49.1, 49.2) in Appendix B. Therefore some further research in this matter is desirable.

5. EXAMPLE 1: TRANSIENT STABILITY OF A SYNCHRt3;',OUS MACHINE

Consider a system described by the equation contains no equilibrium except So. It follows that fl~.(x) is a stability region for S o. Finally, since lel may be chosen arbitrarily small, define fl(x) lim fl~(x) for ~.--+0. Then t2(x) is a stability region for So in the sense that any trajectory starting in tl(x) converges to So, except the unstable equilibria at the vertices of O - ( x i ) and q~+(x~ ). Relative to theorem 4 some remarks are in order:

DXz+aDz+JIbz)=O,

a,b>O

(19)

J(u) is periodic with period 2n, and for 0 < r 1 <2n and for k an integer:

f ( 2 k n ) = O,

f~U(2kn)>O,

.f(:l + 2kn)=O, f(1)(T 1 + 2kn)
New direct Lyapunov-type method for studying synchronization problems

145 ~rm

The system has the stable equilibria Sk = {bz = 2kn, Dz=0}, and the unstable equilibria Ik={bz=zx +2krt, Dz=O}. Equation (19) may represent a phase lock circuit, having a feedback filter with transfer function Ho(s)=b/(s+a), and a phase comparator characteristic

j(u) = K [ s i n

P Re

(u + ~bo) - sin ~bo],

K>0,

0 < 4~o< re/2

(20)

~bo is the phase error at synchronization, which depends on the frequency of the input signal. Equation (19) may also describe an unregulated synchronous generator, connected to an infinite bus, and subject to the usual simplifying assumptions of a linear damping torque, and machine representation by a constant E M F behind the transient reactance[l]. For the system (19), H(s) = b/s(s + a), such that (1/09) Im H (j~) = - ba/co 2 (a 2 + 602) < 0 for all real 09. Hence, by theorem l, the system is nonoscillatory. The substitution bz = by + zl transforms (19) in

DZy+aDy+ f~(by)=O,

FIG. 3. E x a m p l e 1: frequency c u r v e s H~(i.;-r) for v a r y i n g values of r; 0 < r I < a/2 < r 2.

first cycle can be obtained by the procedure outlined before. For the present example it is instructive however, to compute directly the largest stability region, bounded by straight line segments. To this end let

X l l ~'~-X 12

(21)

Yq2 = - a x 1 2 - fl (bx11) with fl (u)= f (u + z 1). The linearization of (21),

D2y+aDy+knby=O,

ko~f~l)(0)<0,

(22)

has one characteristic value, 21 A ( _ a + (a 2 _ 4kob )1/2 )/2, in the open right half plane, while the other one, ).2a=(-a-(a2-4kob)l/2)/2, is stable. Assume that fl (u) satisfies inequality (10), with k a < 0 and k4>0, and rewrite (11) as --+ReHa(jog-r)>0 k4 -

allreale~,

(23)

be a vector representation of equation (21), with statevectorxl=[xll x12]=[y Dy],andlet

v(xl)a=x12+px11=O,

#>0,

be a straight line passing through I o as shown in Fig. 4. Then ~(X1 ) =3~12 -1-/t/X11 = (F/-a)xt2 --ft (bxl 1)

k3

= - # ( / / - a)xl i - j l ( b x l l ) ,

where

H1 (s) a=H(s)/[1 + k3H(s)] = b/(s 2 + as+ kab ). Plots of Hl(jo3-r), for different values of r in the permissible interval 0 < r < - 2 2 , are displayed in Fig. 3. T h e diagrams reveal that, with respect to condition (23), the optimal choice for r is r=a/2. Then (23) holds if l / ( k 4 - k 3 ) + H l ( - a / 2 ) > O , which, after some calculations, simplifies to

a > 2x/k4b.

(24)

By theorem 3, (24) ensures synchronization from all initial states within a finite number of cycles. An estimate of the synchronization region within the

ifv(x~ ) =0. Some more manipulations show that, on the line V(Xl ) =0, x11~(x1)= - ( p 2 - a # + k a b ) x ~ l - X l l [ f l ( b X l l -k3bx11] <0

for

) (25)

x 1~0,

p2 - a p + k3b>0 ,

(26)

using (10). Selecting # as small as possible, such that (26) holds, yields, with k3 < 0: a a + (a 2 -4k3b) 1/2 #=#1 -

2

~-e,

(27)

146

E . J . NOLDUS Dz=xl z

t~2 +/./.iXll =0

•0. r.o~( x

) bz" bxII+rL X._-I-/J... X. = 0

"W'R I =E" W'RI=-E

FIG. 4. Synchronization region fl,,.~(x ), bounded by straight line segments.

e a small positive number. Next, rewrite (25) as

(b) Ifx 1 e ~)max(X), W( X 1 )~---0 and ]g] is sufficiently small, then

Xl i t~(xl ) = - (/t 2 - alt + k.~b)x21

(xl) = 2 ( 0 - a)(w'x 1)2 + 2w,xlO(x21 )

- x l 1[.[~( b x l l ) - k4bxl ~]

= 2 ( O - a ) e z -~ 0(g3) > 0.

>0

for

xl@0,

if

~2-a#+k4b
(28)

again using (10). The condition (28) implies that (24) must be satisfied. Then, selecting /~ as large as possible, such that (28) holds, results in

#=J/2

A a + (a 2 -4kgb) 1/2 2

=

e,.

(29)

It follows that, if (10) and (24) hold, then, as time increases; no trajectory can leave the set flma~(X) {I_ 1Q'floQ'- 1I - i} defined in Fig. 4. Finally let A

W'X1 ~ X12 ~ OXI 1,

such that

The same argument applies to I_ 1. As before this implies that any trajectory, starting in f~madX), converges to the origin So, except the unstable equilibria I_ 1 and I o on its boundary. ~m~x(X) is the largest stability region, bounded by straight line segments, that can be obtained by the proposed method. The condition (24) means that the system (19) must be sufficiently damped. The restriction may be rather severe, such as in the case of power systems. To be specific, consider a round-rotor synchronous machine, described by its normalized equation (19), where b = 1. flu) is given by (20), with K = 1, and q5o the equilibrium torque angle. (24) becomes: a > 2x/k ~, where k4 depends on ~bo. For example, for ~bo=0.2, a numerical calculation reveals that k~ =0.27, hence amin= 1.04. As a second example, consider a machine of the salient type, with

W'..~1 -----(0 -- a)x12 -- bf(ll)(0)x11 jr_ 0(x12l ),

flu) = 0.935[sin (u + 4~o)- 0.2 sin 2(u + q~o) where we have substituted ./'1(u) by f]x)(0)u + 0(u2). Select 0 > 0 such that 02 - aO+ bf]l)(0) = 0, or a + [a 2 - 4bf ]1)(0)] 1/2 0-

2

Then w'21 = (O-a)w'xl +0(x21). Now, as in theorem 3, the quadratic function W(xa)= (W'X1)2--g 2 has the properties that (a) The line w'xl=O lies outside the region flm,x(x), since k 3 ~ fill)(0 ) =
-

sin 4)0 + 0.2 sin 2~o].

(30)

The coefficient 0.935 has been so selected that the maximum torque is the same as for the round-rotor machine. Now k 4=0.23, and am~,=0.96, the decreased value of amin being due to the stabilizing effect of pole saliency. Figure 5 shows the stability domain ~)max(X), for a system with a = 1, 0o=0.2, and f ( . ) defined by (30). The minimum values for a may be decreased, if the power systems are endowed with fast, stabilizing speedgovernors.

New direct Lyapunov-type method for studying synchronization problems f=(u)

I

147

j

I I

I

I

',xl ^ '~,~,

~A2

,.-~2 I\-~,

l Dy

.1-

IN,k3(r/) I

I t

FIG. 5. ~max(X) for a synchronous generator of the salient type, with a = 1 and ~bo = 0.2.

0'3

If the conditions of theorem 3 cannot be satisfied, it may still be possible to construct synchronization regions, by combining the sectors ~(Xx) with a suitable Lyapunov function, for example the function associated with the proof of the system's nonoscillatory behavior. As a simple example how this can be done, consider an unloaded synchronous machine, with equation D2z + aDz + sin z=0.

(31)

Constructing the Lyapunov function V(x) as explained in Appendix A, where we let 5 = 0 , produces V (x) =½(Dz) z - cos z + 1,

(32)

(/(x) = - a(Dz) 2 < O,

(33)

so that

and l?(x)=0 implies 2 = 0 . Hence the system is nonoscillatory. Letting y = z + re, transforms (31) in D Z y + a D y + fl(y)=O, For any scalar r/, let

f 3 (t])

fx(u)=-sinu.

and 14(r/) be such that

f3(rl)uZ
"~L

FIG. 6. Synchronization regions for weakly damped systems.

within the first cycle. Similarly, ~'-~3 is a stability region for $1, while a trajectory starting in ~2 or f14 converges to either So or Io or $1. Hence the union of f~x, ~2, f~3 and f~4 is a synchronization region within two cycles. Larger synchronization regions can be found by taking the union of these sets of all r/. Figure 7.1 displays synchronization regions for a =0.5 obtained in this way. In Fig. 7.2 synchronization regions within the first cycle are shown for a = 0.3 and for a---,0. The latter region approaches the set {½(Dy)2 + cos y =<1 ; - 2n < y < 0}. It is easy to show that this is the exact stability region for the unloaded, undamped generator. 6. E X A M P L E 2: SYNCHRONIZATION O F A PHASE LOCK CIRCUIT

This example deals with the acquisition problem of a phase lock circuit, whose feedback filter has the transfer function 1 +as Ho(s)=~+bs,

b>a>O.

Then in equation (1), p(s)= 1 +bs and q(s)= 1 +as. f ( u ) is given by (20). For this system, synchronization regions within the first cycle have been obtained by Holtzman and Rue[9], using a classical

148

E . J . NOLDtJS ,Dy

(1)

f

"O

2

Dy

(2)

FiG. 7. Synchronization regions for an unloaded r o u n d - r o t o r s y n c h r o n o u s generator. I. Case a = 0 . 5 : ( ~ ) = b o u n d a r y of synchronization region within the first cycle; (/~) = b o u n d a r y of synchronization region within two cycles. 2. Boundary of synchronization regions within the first cycle: (:~)= case a = 0.3; (/~)= case a--* 0.

L y a p u n o v approach. Let us derive the conditions which ensure synchronization from all initial states. It is easy to verify that H ( s ) = (1 +as)/s(1 +bs) satisfies (5.2), hence that the system is nonoscillatory. Letting z = y + n - 24)o yields equation (7), where

and k 4 > 0 depend on q~0, and rewrite ( 11 ) as ....

(34)

+ Re Hl ( j c o - r ) > O `

k 4 -- k 3

ft (u) --K[sin (u + n - ~bo } - sin ~bo].

k =0

k=k o

X~

In Fig. 8.l the root locus plot is shown for the equation . . . . . .

D(bD+l)y+k(aD+l)y=O,

-~.
:k
: k>O

lm

i.e. for the linearization of'(7) with f l ( u ) = ku. The plot reveals that the system, when tinearized at the equilibrium Io, with k = k o = f ] l ) ( O ) < O , has one characteristic value, 21, in the open right half plane, while the other one, 22, is stable and smaller than --

Re lw,

S O .

"

qL"\

b/J

is the break point in the root locus plot, which is closest to zero. Let f l ( u ) satisfy (10), where k3 < 0

r~O FIG. 8. Example 2: I. Root locus plot for the linearized equation. 2. Frequency curves H1 (J~'~- ~') for varying values of r :0 < r 1 <.% < r 2 < I"h.

New direct Lyapunov-type method for studying synchronization problems

149

with

k-%

n~ (s)= (1 + as )/[s(1 + bs) + k 3 (1 +as)]. Figure 8.2 displays the curves H1 ( j o g - r ) for varying values of r in the permissible interval 0 < r < - 2 , . A straightforward calculation reveals that ReHl(je)-So)>Hl(-So) for all real e). Hence, with respect to condition (34), the optimal choice for r equals r = So, and (34) becomes:

.

.

.

.

//I k=O

.

k=kn

k=k°

1/(k4 - k3)+Hl(-So)>O. This is equivalent to Im

k/ k , ga< N/ ai - ~ / a/1

b1

(35,

By theorem 3, (35) ensures that (a) All solutions converge to an equilibrium state, i.e. the system synchronizes from all initial states, after skipping a finite number of cycles. (b) The system synchronizes within the first cycle for all initial states on the z-axis, between the two neighboring unstable equilibria.

Re

r=t" I

A THIRD-ORDER SYSTEM Finally consider a third-order case, where 7. E X A M P L E 3:

FIG. 9. Example3 : 1. Root locus plot for the linearizedequation. 2. Frequencycurves HI (io~- r) for varyingvalues of r; 0 < r l < p
3.

(1 + as)

H(S)=s(l+bs+cs2 ), a,b,c>O. This system may represent a synchronous generator, equipped with a fast speed governor with transfer function 1/(1 +as). Theorem 1 yields, as a sufficient condition for nonoscillatory behavior:

a > c/b.

(36)

If (36) is not satisfied, it may still be possible to prove nonoscillatory behavior by means of theorem 2, at the expense of imposing some condition on the nonlinearity. Let us assume that (36) holds. When linearized at an unstable equilibrium, the system's equation

consider the case where the plot has three real break points. Assume that (12) holds, and again rewrite (13) as (34), where now H 1(s)= (1 +as)/[s(1 + bs+cs2)+ k3(l + as)]. The curves H l(jo~- r) have been sketched in Fig. 9.2. Note that the sum of the characteristic values of (37) is independent of k. This implies that s o < d, such that s o is a permissible choice for r. It follows that if

ReHl(flO-So)>Hl(-So)

allrealeg,

(38)

then r=s o is the optimal choice for r, and (13) is equivalent with

D(1 + bD +cD2)y + ko(1 + aDy)y=O, ko L J~ll)(O) < 0 ,

has exactly one characteristic value in the open right half plane, as can be verified using the R o u t h Hurwitz criterion. The other two characteristic values are stable, i.e. they satisfy a condition {Re s < - d ] for some d > O. Suppose that the root locus of the system

1/(k 4 - k3)+ Hl(-So)>O.

(39)

If (38) does not hold, the optimal r is smaller than s o. An analytical solution cannot be found, since the break points of (37) satisfy a third-order equation. For the numerical case a = 0.75, b --2, c = 1, one has So - 1 - (x/~/3), (38) holds and (39) becomes

D(1 +bD+cD2)y+ k(1 +aD)y=O, -~:
+oc,

(37)

has the general form sketched in Fig. 9.1, i.e. we

k4<~[2x/3-3]

(40)

(40)is analogous to condition (35)in example 2.

150

E.J.

8. CONCLUSION A new technique of the L y a p u n o v type has been p r o p o s e d for the stability analysis of a class of d y n a m i c systems which can be r e p r e s e n t e d in the form of a n o n l i n e a r feedback l o o p with the following characteristics: The l o o p c o n t a i n s a linear element with a characteristic value at the origin, and a time i n v a r i a n t p e r i o d i c nonlinearity. Such systems have an infinite n u m b e r of s t a t i o n a r y states, a n d are a d e q u a t e for s t u d y i n g i m p o r t a n t classes of technical p r o b l e m s , including the transient b e h a v i o r of a s y n c h r o n o u s generator, a n d the s y n c h r o n i z a t i o n of phase lock systems. The a p p r o a c h p r o c e e d s in two distinct steps, the first of which is the p r o o f of convergence of all b o u n d e d s o l u t i o n s to an equilibrium state. The m a i n c o n t r i b u t i o n of the p a p e r consists in d e v e l o p i n g frequency criteria for such n o n o s c i l l a t i n g systems, which (a) ensure that all s o l u t i o n s converge to an e q u i l i b r i u m state, i.e. the system synchronizes from all initial states, eventually after s k i p p i n g a n u m b e r of cycles. (b) permit the c o n s t r u c t i o n of s y n c h r o n i z a t i o n regions within an a r b i t r a r y , given n u m b e r of cycles. T h e c o n s t r u c t i o n of these regions involves the c o m p u t a t i o n of indefinite s o l u t i o n s to an algebraic Riccati e q u a t i o n , in which the system m a t r i x is unstable. Systematic c o m p u t e r a l g o r i t h m s for this p r o b l e m seem as yet to be unavailable. Therefore this p o i n t needs some further research, in view of the a p p l i c a t i o n s to high o r d e r systems. T h e o b t a i n e d criteria also ensure s y n c h r o n i z a t i o n within the first cycle, from all initial states on the z-axis, between the two n e i g h b o r i n g unstable equilibria. This p r o p e r t y is of c o n s i d e r a b l e technical significance, a n d m a y be verified by checking a simple frequency c o n d i t i o n . The m a i n weakness of the m e t h o d is that s y n c h r o n i z a t i o n regions can only be easily o b t a i n e d for systems which are sufficiently d a m p e d . The e x a m p l e s that have been presented are of low order. C o n s i d e r a t i o n of h i g h e r - o r d e r a p p l i c a t i o n s does not i n t r o d u c e a n y new c o n c e p t u a l difficulties, but c o m p u t a t i o n a l p r o b l e m s m a y arise. The optimal choice for the p a r a m e t e r r in t h e o r e m s 3 a n d 4 must then be o b t a i n e d by a c o m p u t e r search a l g o r i t h m , a n d the required c o m p u t a t i o n a l effort increases with increasing system order. The same r e m a r k holds with respect to the c o n s t r u c t i o n of s y n c h r o n i z a t i o n regions, where the a l g e b r a i c Riccati e q u a t i o n to be solved grows in size with the o r d e r of the syslem. The theory has been e x p o s e d for systems c o n t a i n i n g a single nonlinearity. H o w e v e r the extension to certain classes of m u l t i v a r i a b l e systems, c o n t a i n i n g several nonlinearities, should cause no essential difficulties. Specifically it m a y be inleresting to look for extensions to p o w e r systems with nonlinear d a m p i n g torques, and to multim a c h i n e p o w e r systems.

NOLDUS REFERENCES [1] J. L. WILLEMS: Direct methods for transient stability studies in power system analysis. IEEE Trans. Aut. Control AC-16, 332 341 (1971). [2] R. W. BROCKETT: The status of stability theory for deterministic systems. IEEE Trans. Aut. Control AC-II, 596--606 (1966). [3] W.J. GRt:t:N:TheoryofAE('synchl*mi/ati~m I'r,c. IRE 41, 1043-1048 11953). [4] G. W. PRESTON and J. C. TELLIER: The lock-in performance of an AFC circuit. Proc. I R E 41,249--251 (1953). [5] S. WEISSENBERGER: Application of results from the absolute stability problem to the computation of finite stability domains. IEEE Trans. Aut. Control AC-13, 124125 {1968). [6] J. A. WALKERand H. N. MCCLAMROCH:Finite regions of attraction for the problem of Lur'e. Int. d. Control 6, 331 336 (1967). [7] N. DHARMA RAO: Routh-Hurwitz conditions and Lyap-

[8] [9]

[10] [11] [12] [13]

[14]

unov methods for the transient stability problem. Proc. Inst. Elec. Engr. 116,539-547 (1969). M. A. PAt, M. A. MOHAN and J. GOPALARAO: Power system transient stability regions using Popov's method. IEEE Tran.v Power App. Srst. PAS-89, 788 794 (1970L J. M. HOITZM'~N and A. K. Rt~t-: Regions t~f asymptotic stability for phase-lock loops. IEEE Trans. Power App. Svst. PAS-10, 45 46 (1964). J. L. WILLEMS:The computation of finite stability regions by means of open Liapunov surfaces. Int. J. Control 10, 537-544 (1969). E. NOLDUS: On the stability of systems having several equilibrium states. Appl. Sci. Res. 21,218 233 (1969). E. NOLDUS,A. GALLEand L. JOSSON:The computation of stability regions for systems with many singular points. Int. J. Control 17,641 652 (1973). J. P. LASALLE:An invariance principle in the theory of stability. In Differential Equations and Dynamical Systems, Eds. J. K. Hale and J. P. LaSalle, pp. 277 286. Academic Press, New York (1967). R. BELLMAN:Methods oJNonlinear Analysis, Vot. 1. p. 125 Academic Press, New York (1970).

APPENDIX A: PROOF OF THEOREM 1 Some difficulty is raised by the fact that in equation (4) the matrix A is singular. Therefore, rewrite (4) as .,~= Aox - I?1~,(c'x ),

where A o ~ = A - k o b c ' and Jo(U)a=J~u)-kou, k0 a scalar. For k0 4=0, Ao is nonsingular, since q(0) @0. Define V(x)ax'vx-~t~xfolO)dO,

V=V'

(41)

Then f/(x) = x' Vie + 5c' Vx - :~fo(c'x )c'Y~, or, after substituting x by means of x ~ Ao 1[~ + bYo Ic'x)], (/(x)=Yc'[VA O-' + (Ao 1)'V]~ + ~'(2 VAo-~b - ~c)fo(c'x).

This can be written as 12(x)= _ (q,.i)2_e,~,M,,~, aM=~M'>0

(42)

A'oW + W A o = - q q ' - e M } 2Wb_~(Aol),c=O j >

(43)

if

where we have defined W a__(A ° 1 ), V A o i

(44)

(42) shows that 12(x) < 0, and that I2(xt = 0 implies ~ = 0, which is the desired property. By the sufficiency part of the KalmanYakubovich lemma, the system (43) can be solved for W = W'

New direct Lyapunov-type method for studying synchronization problems

By the assumptions of the theorem, A ~ + rl has one characteristic value in { R e s > r > 0 } , and ( n - 1) characteristic values in {Res < r - d < 0}. But (13) implies that all matrices of the form

and q, if [el is small enough and if

(a)

I ® A o + Ao ® l is nonsingular

(45) A+rI-kbc',

(b)

c~Rec'Agl(jo)I-Ao)

151

lb>O

allrealw

(45) is satisfied if all sums of the form 2, + 2j ~ 0, where 2, are the characteristic values of A 0 [ 14]. This condition may be assumed to hold, since in the definition of Ao, the scalar ko m a y be chosen freely. Furthermore, using the identity

have the same number of stable characteristic values, and the same number of unstable characteristic values. Hence, since 0 = k 3 ~ k o ~ k 4,

(53)

A + rl also has one characteristic value in {Re s > 0}, and (n - 1 ) characteristic values in {Re s < 0}. This in its turn, together with (49.1), implies that P is indefinite. Next split up At, P, b and x 1 accordingto

c,Ao l (sl - A o )(sI - A o )- 1b =-c'A o i b, or,

sc'Aol(sI-Ao)-lb=_c'(sl-Ao)-lb+c'Aolb, (46) can be written as

)-t

AI=

1

ot R e _ _ c ' ( j o g I - A o ) - l b jo9 all real to.

j(o kl + koH(jm)J

k3
(46)

I

0 .

0...0

At,.

0

P'3

h=

P

Xl=

This is equivalent to (5) by letting ct = 1 and ct = - 1.

A P P E N D I X B: P R O O F O F T H E O R E M 4 There is no loss of generality in assuming that k 3 = 0. Deriving P(x I ) and substituting xt by means of(8) yields P(Xl )+ 2rP(xl )= x'l[P(A + r l ) + (A + rl)'P]xl -(2Pb-c)'xlfl(c'xt)-(l/k4)f2(c'xl)-g(c'xl)

(47)

where g(u)a=fl(u)[u - (1/k4)fl(u)] > 0 for all u. Writing this in the form 15(Xl ) + 2rP(x 1) = - [q'x I + k21/2fl (c'x I )]2 -ex'lMxl-g(c'xt),

eM=eM'>O,

-- xrlSXl

(54)

where S = S' a=(q + kok21/2c)(q + kok21/%),

P (A + rI) + (A + rI)'P = - qq' - eM

(49.1)

2Pb - c = 2 k~ 1/2q

(49.2)

(48) shows that (14) holds. By the sufficiency part of the K a l m a n Yakubovich lemma, equations (49) can be solved for P = P' and q, if lel is sufficiently small, and if I ® (A + rl) + (A + r l ) ® / i s nonsingular

x'I[P(A1 + r l ) + (A 1 + r l ) ' P ] x I =

(48)

and identifying the right hand sides of (47) and (48) produces the system of algebraic equations

(a)

21 > 0 is the unstable characteristic value of A1, and A t 2 has the (n - 1 ) stable characteristic values of A t. A 1 may be assumed to have the indicated form, by adopting a suitable choice of the state vector x 1. Substituting (52) in (49.1) and using (49.2), yields

+ ko[l - {ko/k, )]cc' + ~:M > 0 because of (53). Letting x~ ~ = 0 in (54) shows that P2tA12+rl)+(At2+rl)'P 2 =-$2,

$2=S~>0

A 12 + rl is stable, having all characteristic values in {Re s < r - d <0). Hence P2 >0. This implies that (15) holds, and that O(xl) can be written in the form (17), where, for the adopted choice of the state vector x 1 W'X I A X I I "

(50) Let

(b)

k2t+Rec'(flol-A-rl)-Xb>O

all real ~o

(51)

(50) means that all sums of the form 0~+ 0~ + 2r ~ 0, where 0~ are the characteristic values of All4]. This condition may be assumed to hold, since r is a free parameter. (If necessary, modify r by an infinitesimal amount.) (51) is equivalent to (13). Next let f l ( u ) A=kou+ f t(u),

ko a=f~l)(O).

Wtxl)~-x~l - & If x 1 e O ( x 1) and W ( x 1) = 0, then x~ 1 = e2, and n-I

E z,~--
such that Ix1[ has the order of magnitude e. Now

Then (8) becomes

W ( x , ) = 2xll [21Xl 1 -- b l f l (c'xl)] J

= 221x~1-2blXllfl(C'Xl).

21 = A I x I - b f l (c'x I ),

where A1 A A - kobc'.

(52)

(55)

In the right hand side of (55) the first term equals 2).,& > 0 , while the second term has the order of magnitude ga, since ] ; " q 0 ) = 0 . So W ( x 1) > 0 if [eI is sufficiently small.