Steady-state analysis of nonlinear dynamic systems with periodic excitation based on linearization in harmonic space

Analyse en regime permanent des systemes dynamiques a excitation periodique par linearisation dans I'espace des harmoniques

By A. Semlyen, Department of Electrical Engineering, University of Toronto, Toronto, Ontario.

This paper considers a very general formulation of the differential equations of a dynamic system in periodic steady state. These are linearized around an operating point, in algebraic form in terms of incremental harmonic phasor components. The solution is iterative, either Newton-type with quadratic convergence in the neighbourhood of the solution, or it has linear convergence if the Jacobian is not updated at each iteration.

On considere pour un systeme dynamique la formulation generalisee des equations differentielles en regime permanent. On linearise ces equations autour d'un point de fonctionnement sous forme algebrique par moyen des composantes harmoniques incrementales. La solution est obtenue par une methode d'iteration, soit de convergence quadratique de type Newton, ou de convergence lineaire si la matrice Jacobien n'est pas recalculee pour chaque iteration.

Introduction

The analysis of linear systems in sinusoidal steady state is conve­niently performed by means of traditional phasor methods and it will be shown that , by generalization, the analysis of non-linear dynamic systems in periodic steady state becomes an efficiently tractable algebraic problem in the space of harmonic components . References 1 to 3 describe different solution methods and reference 4 gives a review of non-linear analysis procedures of electric cir­cuits. This paper presents a solution method based on linearization in the harmonic space. Contrary to the well-known equivalent linearization and describing function approach in reference 5, where a global transfer function is established by retaining only the fundamental, the method in this paper is based on a local, in­cremental linearization with all significant harmonics being retain­ed.

The motivation for this analysis came from a study of the non­linear magnetizing branch of a transformer as an important source of harmonics in power systems. Since the power system voltage is very nearly sinusoidal, the same is true for the flux in the transformer. Therefore, the usually accepted simplifying assump­tion in such an analysis is that of a base sinusoidal state with har­monic perturbations. However, there are many situations when the flux is strongly distorted, possibly even because of the voltage drop in the winding due to the non-sinusoidal magnetizing current. This clearly indicates the need for a harmonic inpu t /ou tpu t model , which is algebraic in this particular application. Such a model could be produced by linearization of the input /ou tpu t characteristic, and it was found that it is remarkably simple if phasors are used for all harmonics. However, linearization permits solving not only algebraic problems, but also the problem of dynamic systems described by ordinary differential equations with periodic inputs , outputs and state variables. This, of course, is a problem with many applications and it was felt that , at this stage, no particular ones should be examined but only the solution of the general problem should be presented which, while conceptually simple, is very broad in scope and applicability if used with appropriate numerical techniques.

The end result of the procedure presented in the paper is a linearized form of the system differential (or algebraic) equations

in the harmonic space. It contains a harmonic Jacobian matrix and a mismatch vector at each iteration point . The latter has to be recalculated at each step for accurate solution. If the Jacobian is also updated at each step, we have the Newton-Raphson method and we expect quadrat ic convergence in the neighbourhood of the solution point; otherwise the convergence will be only linear.

The procedure

The system is described by n equations

AY,x,u°) = O ( l )

in n unknowns , y representing a redundant set of state descriptors of which the state variables x are a subset. In the case of purely algebraic problems, the derivatives x will not appear in Eq . (1). The input

u° = « (0 (2)

is periodic. The superscript 0 indicates that the function is known and it will be used also with other variables. It is further assumed that the steady state solution y(t) and other related variables are also periodic.

We linearize Eq . (1)

f° + FyAy + I* Ax = 0 (3)

where Fy and Fx denote Jacobian matrices.

Any element of the vectors f°, Ay, Ax is periodic and will be replaced by a vector of harmonic phasor components : f°, Ay, Ax. It will be shown in the next section that Eq . (3) becomes

f° + F°yAy + FiAx = 0 (4)

where F°y and F*i are Toeplitz 6 matrices. We note that

Ax = Z>AJC (5a)

Can . E l e c . Eng . J., Vol. 11 No . 3 , 1986

114

SEMLYEN: NONLINEAR DYNAMIC SYSTEMS 115

where

(5b)

(6)

Equat ion (6) is the end result of the linearization in harmonic space. It is a simple algebraic relation between the harmonic vectors Ay and - / ° - , used to update the state descriptors^, if FJ? and Fs are also updated at each step, the method is a Newton-Raphson pro­cedure; if FJand Fj are kept constant then we have a linear iterative approach.

In the case of purely alegbraic problems, Eq . (6) takes on the simpler form

F°yAy = -f°

Linear increments in the harmonic space

Given the scalar function of a scalar argument

its linearized form is

Av=f'(!t)AZ

(6a)

(7)

(8)

We assume that all variables, £, rj, A£, Arj, and / ' ( f ) are periodic functions of t ime. Then,

(9a)

(9b)

(9c)

We substitute Eqs . (9) into Eq . (8). To simplify the calculations, we consider only one term of Eq . (9c) at a time

(9d)

Then, equating exponentials of the same order, we obtain

AVk = ck-hA£h (10)

This equation can be put in matrix form

Arj = FA£

where

A{ = [ . . . , A£- 2 , A£_!, A£ 0 , A £ 2 , . . . ] T

Ar) = [ .. ^ A ^ . A r y . ^ A r / o , Ar?!,Ar/2, . . - ] r

(11)

(12a)

(12b)

F =

C2 Ci C 0 C-i . . .

. . . C2 Cx C0 C_, . .

. . . C2 Cx Co C-j

(12c)

The matrix F is of Toeplitz form (equal elements in a diagonal) and corresponds to a Jacobian. In the case of a linear function (7), only c0 is non-zero and , in general, with only a limited number of harmonics being considered, F is band-diagonal Toeplitz.

Eq . (11) explains the transition from Eq . (3) to E q . (4).

Example 1 The general procedure culminating in Eq . (6) can be more easily

unders tood through an example. For the system Eqs . (1), we choose, with x = yx and.y = \yx, y 2 ] r , the following two equations:

fx = 2yty2 - sin2r = 0

fi = 2yxy2 - cos2r = 0

(13a)

(13b)

where the input is of second harmonic only (sin2r and cos2r). Because of the coupling between harmonics , due to the nonlineari-ty of Eqs . (13), we may expect that other harmonics will also be pre­sent in the state d e s c r i p t o r s ^ a n d ^ 2 . For simplicity of illustration, we assume that , at the following stage of calculation, yx consists of only the first a n d ^ 2 of the third sinusoidal harmonic . Therefore, we linearize Eqs . (13) a round

y\> = sinr

y% = sin3f

Then Eqs . (13) become

2y °xy°2 - sin2f + 2y?Aj>2 4- 2y°2Ayx = 0

2y\y% - cos2f 4- lylAyx + 2>>?A.y2 = 0

or

2cos/sin3r - sin2r + 2cosf A.y2 + 2sin3/A^j = 0

2sinrsin3r - cos2r + 2s in3 /A^ + 2sin/A>'2 = 0

(14a)

(14b)

(15a)

(15b)

(16a)

(16b)

It is not intended to perform calculations analytically, since this is possible only in contrived problems. Instead, the periodic func­tions are evaluated numerically and then a Fast Fourier Transform (FFT) is used to obtain the harmonics . Then, in Eqs . (16) we obtain sin4r for the expression 2cosrsin3f — sin2r, and — cos4r for 2sinr sin3r — cos2r. Since

sin4r = 2j

(17a)

we have only two harmonics , of order —4 and -f 4 , and of magnitudes j/2 and —y/2, respectively.

In Eq . (16a) the coefficient of A.y2 is

2cosf = eJt + e~Jt (17b)

so that only the harmonics of order — 1 and + 1 exist, bo th of magnitude equal to 1. The Toeplitz matrix will therefore have two

00

— 00

oo

A£ = A£hej{hut)

— OO

00

Ar)=J^ ArjkeJlkut)

— 00

D = col{diag(/iaj)l

h denoting the order of the harmonics.

Substitution of Eq . (5a) into (4) yields

(F°y + F°xD)Ay = -/°

116

diagonals with entries equal to 1. Similarly, in Eq . (16a) the coeffi­cient of Ayi is

CAN. ELECT. ENG. J., VOL. 11 NO. 3, 1986

(24a)

(17c)

so that the Toeplitz matrix will have two conjugate diagonals with entries j and - y , for the only harmonics of order — 3 and -I- 3 respectively.

After solving Eqs . (6) derived from Eqs . (16) for [AyuAy2], new values are obtained for >>? a n d j > 2 . These now contain a number of harmonics and we certainly have to use them i n / 0 of Eq . (4), i.e. in the first terms of Eqs . (15). The Toeplitz matrices, however, do not have to be recalculated unless we wish to have a Newton-Raphson solution of our nonlinear equations. In any case, the solution will converge to the accurate one except for the effect of harmonics which may have been neglected or have been truncated in the discrete Fourier transform process.

Example 2 Consider a nonlinear inductance supplied from a sinusoidal

voltage source through a constant impedance. We have then

v di d\I/ „ . ^ Ri + -jt =£ m s ina>f (18a)

where is assumed to be of the simple polynomial form

t = ai+bP (18b)

The harmonic phasor solution of Eqs . (18) is Iu I3, Is,I7j... bu t , for the simplicity of illustration, we shall consider only the first three harmonics in Eqs . (18), i.e.

k

(19a)

(19b)

with

h,k= - 5 , - 3 , - 1 , 1 , 3 , 5 (20)

The linearized form of Eq. (18b) is (similarly to Eq . (8) obtained from Eq. (7))

A+ = ^(i)Ai

where

^ ' ( / ) = a + 3bi2

(21a)

(21b)

Substitution of Eq . (19a) into Eq. (21b) yields for f ( / ) an ex­pression of the form (9b)

with

m = -10 , - 8 , - 6 , - 4 , - 2 , 0 , 2 , 4 , 6 , 8 , 1 0

(22)

(23)

(24b)

F rom the known coefficients cm of Eq . (22), we can assemble the Toeplitz m a t r i x F o f Eq . (12b) so tha t , for the incremental variables of Eqs . (24), Eq . (11) becomes

Atf-7 Co C-2 C-4 C-6 C-s C-io" ~A/-s~

c2 Co C-2 C-4 C-6 C-s A/-3

A\p-i c4 c2 Co C-2 C-4 C-6 ALx

Afc c6 c2 Co C-2 C-4 A/i

Afc c 8 c6 c4 c2 Co C-2 A/3

A^s Cio c6 c4 c2 Co A/ 5

(25)

We note that in Eq . (25), the Toeplitz coefficient matrix is t run­cated because of the prescribed number of harmonics contained in the vector A^.

Our basic differential Eq . (18a) (which corresponds to the general Eq . (1)) can be rewritten in terms of / = i° + Ai and \p = \p° -I- A\p, using Eqs . (19) and (24)

Substitution of A\ph from Eq . (25) into Eq . (26) yields an equation containing only the increments AIk. By separating the expressions multiplied by eJ{kut) for the different values of k given in Eq . (20), we obtain a set of linear equations in AT* (which correspond to the general Eq . (6)). These can serve for updating the harmonics Ik.

Computing tasks

After the vector of harmonics y is updated using Eq . (6), we know^y 0 ,* 0 and JC° (using Eq . (5a)). We subdivide a period TintoN equal intervals and evaluate, at the equally spaced discrete values tK, the t ime functions yi, x% , and u%. With these we calculate

f°K=f(yl,x°K,u°K)

and, via an FFT , we o b t a i n / 0 of Eq . (6)

(27)

For Eq . (6), we also need to process the partial derivatives of the Jacobian matrices F ° and I* of Eq . (3). Any such partial derivative can be computed by a central difference formula which, with the scalar notat ions of Eqs . (7) and (8), is

(28)

The incremental form of the variables Eq . (19) is where

\^(0= c - e J l m u , t )

m

k

^MVH + = ^ (**- - e*") (26) h

gjlt _g-j3t

2sin3/ = :

J

Ai = £ AIkej(k^ (AI.k = Alt)

SEMLYEN: NONLINEAR DYNAMIC SYSTEMS 117 (29)

The complex Fourier coefficients, obtained by an F F T program for each partial derivative symbolized by Eq . (28), are then assembled in Toeplitz matrices F as shown in Eq . (12c) and then in the matrices F ? a n d Fj of Eq . (6). The addition in Eq . (6) is a direct sum. Thus Eq . (6) is complete for updating again the state descrip­tors y.

directly from Fourier coefficients. They are thus easy to compute and to manipulate in the solution process. • The result of harmonic linearization makes available phasor pro­cedures to approach the final solution. The possible strategies in­clude a Newton-Raphson algorithm and iterations with a constant iteration matrix. Other iterative algorithms may be more efficient, depending on the particular system structure.

Concluding remarks

The linearization in harmonic space considered in this paper is different from time domain linearization by the fact that in­crements of harmonics are themselves sinusoidal functions of t ime. Therefore, the process of updating in harmonic space corresponds to a closer approximation to the steady state periodic solution of the system.

The following details are important : • All computat ions must be computerized. This includes numerical differentiation and Fourier transforms (FFT). • The Jacobian matrices are band-diagonal Toeplitz, obtained

References

1. M.S . Nakhla and J. V l a c h , " A Piecewise Harmonic Balance Technique for Deter­mination of the Periodic Response o f Nonlinear Systems", IEEE Trans, on Cir­cuits and Systems, Vol. C A S - 2 3 , 1976, pp. 85 -91 .

2. A . Fukuma, "Jump Resonance in Nonlinear Feedback Systems: I. Approximate Analysis by the Describing-Function Method", IEEE Trans, on Automatic Con­trol, Vol. A C - 2 3 , Oct. 1978, pp. 891-896.

3. A .G.J . MacFarlane, Frequency-Response Methods in Control Systems, Wiley, 1979.

4. L .O. Chua and C.-Y. Ng. "Frequency-Domain Analysis of Nonlinear Systems'', IEEE J. Electronics Circuits Syst., Vol. 3 , 1979, pp. 165-185.

5. J .C. Hsu and A . U . Meyer, Modern Control Principles and Applications, McGraw-Hill, 1968.

6. G .H. Golub and C.F. Van Loan, Matrix Computations, The John Hopkins University Press, 1983.