Iterative time-scale separation of synchronous machines

Iterative time-scale separation of synchronous machines

Iterative time-scale separation of synchronous machines M. S. Mahmoud Electronics and Communications Department, Faculty of Engineering, Giza, Egypt ...

679KB Sizes 1 Downloads 68 Views

Iterative time-scale separation of synchronous machines M. S. Mahmoud Electronics and Communications Department, Faculty of Engineering, Giza, Egypt

Cairo University,

A. S. Al-Fuhaid Electrical and Computer Engineering Department, Kuwait University, PO Box 5969, Kuwait

F. A. Saleh Computer Centre, Kuwait University, PO Box5969,

Kuwait

(Received May 1986; revised August 1986)

A scheme for iterative separation of time-scales in linear singularlyperturbed systems is developed and applied to synchronous machines, Both the current and flux-linkage dynamic models are considered and it is shown that they are linearly related for small disturbances. Simulation studies on a typical system have indicated the usefulness of time-scale separation schemes in anaiysing the dynamic behaviour of synchronous machines. Keywords: tion

iteration,

Mathematical models of electric power systems are nonlinear and generally display several interacting phenomena of various dynamic responses. For the case of a small disturbance, it is often desirable to seek linearized models. To analyse the dynamic behaviour of such models, it would seem beneficial to describe the linearized power models in terms of reduced-order models.’ Recent results have indicated the significance of multitime-scale theory* in developing reduced models operating at different time-scales which provide a good approximation to the full-order models. Previous work3 has focused on the analysis of the effects of neglecting stator transients in machine performance. The technique of asymptotic modal decomposition has been utilized in output feedback control design with high gains.4 Subsequently, others5 have used singular perturbation theory with a single parameter on two time-scales to further clarify the meaning of several classical results. Their work has centred around the electrical part of the synchronous machine model. 0307-904X/87/02132-08/$03.00 0 1987 Butterworth & Co. (Publishers)

Ltd

synchronous

machines,

flux-linkage,

simula-

In this paper, an iterative scheme of time-scale separation is developed and applied to linearized models of synchronous machines. Both the flux-linkage and current models are considered and it is shown that the two are linearly related. Therefore the flux-linkage model is focused on. Using typical data, it is concluded that the rotor angle, rotor angular speed and the field flux linkage are slow variables, sufficiently separated from the dq components of the stator and rotor flux linkages as fast variables. Simulation results have indicated the usefulness of the scheme developed in analysing the models of synchronous machines for small disturbances. Development

of dynamic models

The starting point in developing dynamic models of a synchronous generator connected to an infinite bus is the generator’s equivalent circuit. Figure I depicts the generator’s equivalent circuit using the so-called dq reference frame.6 It should be noted that this circuit is applicable when using physical quantities such as

Appl. Math. Modelling,

1987, Vol. 11, April

133

Iterative

time-scale

separation

of synchronous

machines:

M. S. Mahmoud [;I]

et al.

= (3)‘12V, [~si~~~]

(unit-‘)

Observe that 6, the generator’s with respect to a synchronously (i.e. infinite bus), is given by:

(4)

rotor-angle, as measured rotating reference frame

8 = q$ + S + II/2 (rad)

(5)

It is also worth noting that, in power systems, the equivalent stator voltage E, corresponding to the field voltage “r is more commonly used than “r. It can easily be shown6 that:

1

KMF Efd = uvF =

[

(3)1,23F

(6)

(unit-‘)

“F

Combining equations (1) and (2) to eliminate vq, and using equation (6), then:

“d and

I I

wxd

Rotor

Stator

Figure 7 dq equivalent circuit for synchronous generator operating under balanced conditions volts, amps, etc. Moreover,

the same circuit when per-unit or normalized variables are provided that the base or reference values properly.6 In this work, per-unit quantities exclusively for all variables, including time. The dynamic equations of a synchronous can be readily obtained from the equivalent yield: -“d “F

0 ---

=

0

I wL,

0

rF

0

I

0

0

0

rD

1

0

- _

0 I

Ld KMF KMDl

ld LF

KMF

ID

KMD

LF MR

+ ------------------

_-

0

% .iQ_

0 0

0 0

I

aKM,

IF iD -t

0

Id

MR j 0

0

1F

LD;

0

iD

I o 1 L, KM, 0 ~KM, L,

_-

0

I

0

KMD M, LD j 0 ----_____--_---__

0

1 i,

0

iQ

0

0

0,

IKM,

0

O

?-q _

KM, L,

i=brfR,

i,=L,+L,

i, = L, -I L,

Define: i = [id iF i,

R

=

[-“cad

: iq iQ

&do

]’

: -Vmq

01’

diag [i arr rn i ro]

0

(1)

jg=

1 i, IO

0

jo

0

iQ

r

I Ld

KM,

KM,;

aKMF

aLF

aMR 1

KMD _- ----

1

I I

MR LD -----------------

; :

i,

jKM,

v,,(t)

134

Appl.

of phase a of the infinite

bus is given then equation

= (2) iI2 v, cos (wRt + a) (unit-‘)

Math.

Modelling,

KM,

%

For a generator connected to an infinite bus via an impedance of R, + jwLe, the generator’s terminal voltage dq components can be shown6 to be related to the dq components of the voltage of the infinite bus according to:

where the voltage by:

yq

where:

V =

(Unit-‘)

0

0’

aLF aM, 1 0

0

.. 0

0

0

I 0

J

KM, KM,:

0

r

I

A

Ld

ld

6

____

-uKMDf

0

wKMo

0

I

0

I

--

0

- h,Ld -wKMF

0

generator circuit to

r

---------------~___--

-“q 0

can be used considered, are chosen are utilized

1987, Vol. 11, April

(3)

v=Z&+

(7) reduces wfii+ ii

to:

KM, L,

J

Iterative time-scale separation of synchronous machines: M. S. Mahmoud et al. Rearranging gives: . i = -i-l{R + wti}i

+ i-‘v

The nonlinear current and flux-linkage machine models should be used when large disturbances such as shortcircuit faults are considered. However, when only small disturbances such as incremental load shedding or additions are of concern, then the linearized versions of the current and flux-linkage models developed here are sufficient to describe the performance of a synchronous generator. Thus, given an initial operating point (P or A” and o”, 8’) as specified by a pre-disturbance steadystate load flow program, the linearized flux-linkage model is given by:

(10)

which describes the electrodynamics of the stator’s side of the generator. For the rotor’s dynamics, Newton’s second law as applied to rotational motion yields: 2H%dw/dt

= T, - T,

(unit-‘)

(11)

where:

T, = f(i& - i&)

(unit-‘)

(12)

T, is the mechanical

torque input (unit-‘), T, the electrical torque output (unit-‘), H the generator inertia constant (s), w = d8/dt the angular velocity of the rotor (rad(unit time)-‘), Ad and A, the dq components of the stator flux linkages (unit-‘), and t the time (unit-‘). The Aux linkages, $, and currents i, are related via the inductance matrix L:

X=Li

i=AA,x+B,u

x = [Ad - A!1 Ar - A$1 An - AtI A, - Ail Ao-

(13)

equation

j=Ajy+Biu

(15)

i=-{li+bJfi}L-‘X+v

(unit-‘)

dF

d LMQ

1 i

qQ

4

(unit-‘)

LMD

-LMD

-LMD ldlF

-LMD ldlD

ldlF

IF -

LMD

l2F -LMD lFID

(unit-‘)

u = [AEfdAT,,,lf

(22)

where & and T,,, can be varied according tation and speed-governing systems.

to the exci-

(17)

A case study

(18)

In order to analyse the performance of the dynamic models of synchronous machines, typical practical data of a single-generator infinite-bus system6 is used. The various parameters have the following values (unit-‘):

where the mutual inductances LMQand LMD, and the leakage inductances lq, IF, I,, lD and 1, are used to describe the inverse of the inductance matrix.

12 d

(20) and

Bi are of appropriate

L 2HoR; = T, - 4 +&Aql

=

(21)

(21), the matrices A,, B,, Ai, dimensions and can be obtained by using partial differentiation with respect to the proper state and input variables followed by evaluating the resultant partial derivatives at the chosen operating point. Finally, u is given by: In equations

(15) reduces

(16) Now, the complete state-space model of a synchronous generator connected to an infinite bus is described by equations (lo), (14) and (16). Since currents are used in conjunction with S and o, these equations form the so-called generator’s current model. Equivalently, a flux-linkage model can be obtained by using equation (13) to eliminate i. Thus:

L-1

(unit-‘)

(14)

(rad-‘)

Id -

model:

(5) becomes:

With %ase = OR and tbase= l/ma, equation to 6=W-1 (rad(unit time))‘)

dD

current

8’1’

where:

w = oR + dS/dt

LMD Ilh,+hdl-+hd

A&l W- ~“1 6-

For the linearized

= T,,,

-f[LdiqIKMFiqIKM~iqI-Lqidi-KMOid]i Since o = de/dt,

(20)

where:

where A = [A,hrAnA,Ao] and L is the second matrix shown in equation (1). Using equation (13) to substitute for Ad and A, in equations (12) and (ll), leads to: 2Hw,dw/dt

(unit-‘)

,!,d = 1.7

1, = 0.15

lq = 0.15

L, = 1.64

IF = 0.101

/n = 0.055

-LMD ldlD

-LMD lFID

ID -

LMD 2 1,

------I,-

LMQ

12 q -LMQ j

L,LQ

(unit-‘)

(19)

-LMQ lqlQ

l,

-

LMQ

12 Q

Appt. Math. Modelling,

1987, Vol. 11, April

135

Iterative

time-scale

separation

of synchronous

Ld = 1.7

1,=0.15

l,=O.lS

L, = 1.64

lr = 0.101

1D= 0.055

machines:

-3.48831672 1.20223566 2.20774184 -0.03605909 -0.03520841 -0.00079928816 0

KM, = KM, = MR = 1.55 lo = 0.036 LF = 1.651

LD = 1.605

L, = 1.526

KM, = 1.49

L MD= 0.02825

LMQ = 0.02848

r = 0.001096

rr = 0.00074

rn = 0.0131

r, = 0.054

R, = 0.02

L, = 0.4

Bi=

id = -1.591

Ad= 1.676

vmq= 1.025

iF= 2.826

A, = 2.2

V,d = 1.397

i, = 0

A, = 1.914

vd = -1.148

iq = 0.701

A, = 1.15

vq = 1.675

i, = 0

A, = 1.045

6= 53.74

cr=o Using the foregoing system parameters, the nonlinear machine models are first computed and then they are linearized about the operating point, to yield models (20) and (22), where:

B,=

l-

-26.31867 -1000.122 3.78058 0 -115.28011 0 -435.15385 152.3191 0 283.59947 -0.73795 - 1.99696 ‘0 0 324 0 0 441.62334 I x 1O-3 0 0 (23a)

J

0

:: 0 0 0

I

x 10-j

0.55962 0

-2.54783918 0.87810349 1.61251733 0.09010684 -0.12336775 -0.0004422066 0

-0.48866 5.54861 -4.88656

For the operating point, the infinite bus voltage is specified to be V, = l/O” and the machine is assumed to be delivering its rata power at a power factor lag of 0.85. From the specified conditions, the following initial values can be readily computed:

0.14927 0.82917 0 0 0 0 0

et al.

-2.44531 1.751887 0.842767 -0.603781 1.547627 - 1.108762 1.775953 2.386655 -1.734056 -2.330351 0 0 1 0 1 (24a)

H=2.37s f,=60Hz 0, = WR = 377 rad/s

-16.12939 40.37268 1.38622 -5.26798 45.06393 66.92679 A, = 1 999.8833 -236.9350 0 0 1.02376 -0.40185 0 0 664.3501 -455.89655 0 0 0 0 - 176.42042 321.76923 0 -318.3382 1.64193 0 0 1000

M. S. Mahmoud

I

0 0 0

x 10-3

0

0

0 0 0

0 0.55962 0 1

describe the current model. On examining the properties of both models, it is found that Ai and Ah have the same eigenvalues (-0.35948 f j0.99827, -0.0014388 f jO.029199, -0.12854, -0.098768, -0.000609). This confirms the result stated by Anderson and Fouad.6 Thus, A, and A, are related by a similarity transformation and, as such, it is sufficient to consider either the flux-linkage or the current model. For the remainder of this work, and due to its appeal to power engineers, the flux-linkage machine model is used exclusively. For the purpose of simulation, it turns out that the time, t, would be more meaningful if expressed in seconds. This implies multiplying (A, and BJ in equation (23) by a factor of 2tio = 377. On examining the dynamic behaviour of the flux-linkage model it is found that the variables 6, w and A, represent one group sufficiently separated from the other group, A,, Ad, An and A,. More importantly, the first group, 6, w and Ar, has a slower response (corresponding to smaller-eigenvalues) than the second group, A,, Ad, An and ho. This suggests the use of multi-time scale theory* in analysing the performance of the synchronous machine, on a decoupled basis. Iterative separation of time-scales In order to capture the mode-separation property, the flux-linkage model (23) is rewritten in the standard form of singular perturbation:* ~I=AxI+B~2+F~ &i2

=

C.q + Dx, + Eu

XI(O) =

x10

W-9

x*(O)

x20

(25b)

0.00043611 -0.0049519 0.00436105 2.64889018 -2.58639998 -0.00020268354 0

A

136

Appl. Math. Modelling,

=

where: (23b)

xi = [S~hr]

x2

=

[A,,

Ad,

0.01417717 0.07720239 -0.09641016 2.64889018 -2.58639998 -0.00020268354 0

1987, Vol. 11, April

AD,

&I

and E>O is a small parameter expressing the ratio of speeds of fast and slow modes. The matrices in equation (25) have the values (after permutation and scaling):

describe the flux-linkage model, and: -0.03608014 0.01243489 0.02283498 3.58881896 -3.50415481 -0.00000784582 0

(24b)

0 A=0

377 0

[0

0

0 -0.1514 -1.986 0 0.3839 0.522

1

0 0.2782 1.425

0 0.619 0 I

Iterative time-scaleseparationof synchronousmachines: M. S. Mahmoud 166.492 122.148 0 0

121.3 -171.87 0 0

L

C=

57.42 -377.04 o

D= L

-89.324 15.22 25.23 0 I

376.95 -164.05 -6.08 -9.922 16.98 -43.4 0 0

106.94

1

-66.51 250.46 0 -120.012I

k2 = 0, (26)

It is well-known that nl and x2 will differ mainly from xls and x2s by their fast parts. To determine this part, ql is introduced as the difference between x2 and x2s, i.e.: + EU,]

(27)

The use of equation (27) converts equation (25) to: ;i = Alxl+

Bql + Flus

(284

E$ = Clxl + Dlql + F2u,

(28b)

where: A,=A-BD-‘C

Fl = F - BD-‘F

(294

Cl = ED-‘CA,

D, = D + ED-‘CB

(29b)

F2 = SD-‘CF,

(294

Note from equation (29b), that Ci is of order E, i.e. x is weakly coupled to the state vi. More importantly, the QSS of ql, obtained (equation (28b)) as: 0 = Dr’Cl~ls

+ qls + Di’F2U.s

+

Di’[CA

(304

is introduced as the error due to the QSS assumption &qls = 0 and equation (30a) is substituted into equations (28) and (29). Up on repeating this step K-times with:

&ZI[CK-IXII

rl~ = W-I+

Wb)

rlo=x,

(3Oc)

one arrives at the system: ii = Ak_q + Bqk + FlkU

(314

&tik

WI

=

ck+

+

Dk’lk

+

qj-1)

= T)k - x2

(33)

0~~=

xl -

&BD;‘qk

(344

%%i~i

+ Bkirlk + ]F,k - BDk’F,,]u

= [Ak - BDklCk]o!i + [&AklBDkl]qk + [Fik - BD,1F2k]u

P4b)

It is clear from equation

(34b) that Bkl is of order E which is small; this implies that the fast input has been reduced. By comparison with equation (30), equation (34a) is iteratedj times to give: oLj+i= OLi- EBk,Dkj’y)k

(354

q)=x1

(35b)

which when used in equations (31), reduces them to: &j = Akjol, + [F~k - BkjDGlFzk]U

(364

&k = Dkj’lk + [FZk]U

(36b)

It should be emphasized that equations desired decoupled forms, where:

(36) are the

Akj+i=Akj-BI&~: A,=Ak

is only of order E, which implies that ql is predominantly fast. Continuing this process with constant us: ~2 = w

-

is the slow part of x1 and can be expressed as:

+ EUJ

qi =x2 + D-‘[Cxl,

(qj

j=l

j=l

The quasi-steady state (QSS) condition, applied to equation (25b) gives:’ x2s= -D-‘[Cq,

5

In the limit of equation (32) (k = ~0) A, contains all the slow modes and Da/& all the fast modes; both depend on the original subsystem matrices without modal transformation. Note from equations (31b) and (32b) that Ck is of order .sk which is small. Careful examination of equations (30b) and (31b) suggests that:

0 0.02109 0 I

0 F= 0 [0.03129

It should be noted from equation (31) that after k iterations, the flux-linkage model still contains the full fast input Bqk in the slow subsystem (31a) which has to be reduced. To recover x2 from qk and x1 it is observed from equation (30) that:

=(~DjI\Cj-*)X*

0 0 0 0

0 EC 0.0074 0 L0

eta/.

F2kU

B k]+1 =

&A,+

lBkjD&’

BkO=B

Dkj+l = Dkj + ECkBkjDkj’ Dk,,=Dk

(37)

Note from equations (37), that Bkj + 0 as i + 00, i.e. the weakening of the fast input has been accomplished. In general, the slow and fast subsystems of the resulting system (36) are only weakly coupled, and the matrices A, and Dkj are of order Ek’j apprOXimatiOIIS of A, and DOZ. Finally, to recover xi from aj and qk it is observed from equation (33) that:

where: A, = Ak-i - B&?iCk_i

Ao=A

Ck = D-lCk_lAk

co = c

Flk = E - BDi_!,F Dk = Dk - 1 + &D;!lCk-lB Fzk = ED;! 1CkF,,

2

(324 Do=D

(O!s-cWS-l)=~j-X~

S=l

(324 (3W

= -@ S=l

B,,D,-,)q,

(38)

and that ~j has the same physical meaning as xi.

Appl. Math. Modelling, 1987, Vol. 11, April 137

Iterative

time-scale

separation

of synchronous

machines:

Using the values of A, B, C, D, E and F, the results from simulating the foregoing scheme of time-scales show, after lOk-iterations and two j-iterations and an accuracy of lo+ (the tolerance of the numerical computation was initially selected to be lo-$ the k and j iterations recorded were obtained when this tolerance was achieved; indeed, for tolerances less than 10m6e.g. 1O-8 or lo-lo, more iterations will be required) that the slow submodel is described by: xr =

0 -0.322 [ -0.503

377 -1.038 3.782

[

-0.0001 -0.002 + 0 0.02 -0.001 0.03 and the fast submodel is given by: 303.1 -5.619 12.28 2.51

&&=

M. S. Mahmoud

et al.

It is generally regarded that the results obtained for case (3) are the superposition of the results of cases (1) and (2). A sample of the results is displayed in Figures 2-8. 0.08r

0

1.6

I

3.2

I

Time (s) 4.9 6.4

I

I

9.6

8.0

I

1

11.2

I

1 1 (394 0 -0.16 -0.476

x1

u

- 132.4 -10.476 -35.2 -1.81

Simulation results Both the full model (23) and the decoupled model (38) have been simulated using three different input excitations: O]l (1) u = [O.l (2) u = [O

O.l]l

(3) u = [O.l

0.11t

-0.72l 01: (---_)

Figure 3 Behaviour of angular velocity for input [O.l exact model; (- - -1 QSS model

Time (s 1 o”60i 0.144-

0.128 -

0.112-

_ 0.096; .z 3 g O.O%O= z L 0.064-

0.032 -

0.016

Y

OO

J Figure2 Behaviour of rotor angle for input IO.1 model; (- - -) QSS model

138

Appl.

Math. Modelling,

01: (-_)

exact

1987, Vol. 11, April

I

1.6

I

3.2

I

I

4.8 6.4 Time (s)

I

I

I

8.0

9.6

Il.2

01: (--_)

exact

Figure 4 Behaviour of field flux for input [O.l model; (- - -) QSS model

Iterative time-scale c[

$6

Time (s) 4;8 6;4

3;2

8;O

separation

9;s

of synchronous

machines:

et al.

M. S. Mahmoud

“;

T .:

0.096-

2 5 ;

o.oao-

8e 4 0.064

_//___

::

A 0.048

0.032

0.016 I/

OOV Time (s)

-0.18’ Figure 5 Behaviour of d-axis statorfluxfor exact model; (- - -) QSS model

input [0

Figure

0.11: (-1

(-)

7 Behaviour of d-axis damper exact model; (- - -) QSS model

0.144

0 128-

flux for input [O

-

0.128-

0 112-

i

0096-

9 i

0.000-

e : g

I I I

I I

0.064--

b

I

I 0.040.-

I I I

0.032 .-

0.016 -

0’

0

1

I

I .6

3.2

Figure 6 Behaviour (-) exact model;

I

I

4.8 6.4 Time (s)

I

0.0

I

9.6

1

II.,

of q-axis stator flux for input LO.1 0.11: (- - -) QSS model

0

I

1.6

I

I

I

3.2

I

4.8

I

6.4

8.0

9.6

I

II

Trme (s) Figure 8 Behaviour of q-axis damper (-) exact model; (- - -) QSS model

Appl.

Math. Modelling,

flux for input [O 0.11:

1987, Vol. 11, April

139

Iterative

time-scale

separation

Table 1 Steady-stateerror

of synchronous

machines:

M. S. Mahmoud

et al.

fordifferentinputexcitations u= [O.l 01' Model(38)

u= [O 0.11' Model(38)

u= IO.1 0.11' Model(38)

State

Model(23)

6 t

-0.06555414 0.13132711 0

-0.06561475 0.13236202 0.26525x10-7

0.14848186 -0.16169987 0

0.13290154 -0.14249876 0.5305x10~5

0.08292772 -0.03037276 0

0.067286783 -0.01013674 0.5332x 1O-5

A,

-0.05117291 0.11206552 0.12004046 -0.45599033

-0.0003133 -0.0002287 -0.0006698 -0.0004735

0.11917576 -0.17370065 -0.16172927 0.10619485

0.0001408 -0.0000767 -0.0003226 0.0001082

-0.06800285 -0.06164513 -0.04168881 0.0605958

-0.0001725 -0.0003053 -0.0009924 -0.0003653

hd

Jb &I

The following conclusions simulation studies:

can be drawn

Model(23)

from

the

The discrepancy in the trajectories of the slow variables, particularly for large t, depends, to some extent, on the input excitation. As shown in Table 1, the steady-state error is almost negligible. There is close agreement between the trajectories of the original and approximate fast variables during the initial transients. The steady-state errors, however, are quite large (Table I). Both points verify the fundamental results of singular perturb&ion theory.* In turn, the slow model can be used fairly well for long-term studies, while the fast model can be used in analysing transient phenomena.

Conclusions This paper has been concerned with mode-separation in linearized models of synchronous machines. In view of available results,s5 the present work contributes to the further development of multi-time-scale theory and its application to power system engineering. Specifically, the present work extends and builds upon earlier work in several ways: (1) It treats the electro-mechanical model of the synchronous machine as an integrated model and analyses its dynamic behaviour in terms of slow-fast separation properties.

140

Appl.

Math. Modelling,

1987, Vol. 11, April

Model(23)

(2) It employs an iterative system of time-scales to yield better decoupling of the slow and fast subsystems and to considerably improve the accuracies of these subsystems. (3) It applies the developed scheme to an actual synchronous machine to obtain realistic results. Thus, the proposed technique has emphasized the usefulness of singular perturbation for long-time studies of power system problems. For such problems, a large saving in the computational costs of system simulation is generally expected.

References 1 Mahmoud, M.S. and Singh, M.G. Large-scale systems modelling, Pergamon Press, Oxford, 1981, pp 265-325 2 Saksena, V.R., O’Reilly, J. and Kokotovic, P.V. Singular perturbations and time-scale methods in control theory: survey 1976-1983, Aufom., 1984, 20(3), 273-293 3 Krause, P.C., Thomas, C. H. and Riaz, M. The theory of neglecting stator transients, IEEE Trans. Power Appl. Syst., 1973, PAS-92, 417-422 4 Inan, M.K. Modal analysis of synchronous machine dynamics, Electronic Res. Lab., Berkeley, 1981, UCBIERL MS1 165 5 Zaid, S.A., Sauer, P.W., Pai, M.A. and Sariogh, M.K. Reducedorder modelling of synchronous machines using singular perturbation, IEEE Trans. Circuits and Sys., 1982, CAS-29,782-786 6 Anderson, P.M. and Fouad, A.A. Power system control and stability, Iowa State University Press, Ames, Iowa, 1977, pp 83-307 7 Kokotovic, P.V., Allemong, J.J., Winkelman, J.R. and Chow, J.B. Singular perturbation and iterative separation of time scales, Autom., 1980, I&23-32