THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING A STUDY OF NONLINEAR SYSTEMS WITH RANDOM INPUTS 'Kuei Chuang.:.Research Associate Louis F. Kazda Professor of Electrical Associate Engineering December, 1958 IP- 341

A 1 ' I I i - r

TABLE OF CONTENTS LIST OF FIGURES.....o..........o..................... I. INTRODUCTION. o oo.. o..oo o o e................... o II. STATISTICAL IDEA FROM PHASE SPACE POINT OF VIEW.....o III. CONDITIONS IMPOSED ON THE FUNCTIONS AND EFFECT OF FEEDBACK......o....o...,........ooo. IV. MARKOFF RANDOM PROCESS AND FOKKER-PLANCK EQUATIONo.... V. A SECOND ORDER NONLINEAR SERVO SYSTEM AND ITS RELATION WITH FOKKER-PLANCK EQUATION "..o........ VI. HIGHER-ORDER SYSTEMo.o o o.o o o o...o o. o........... o o VII. COMPUTER STUDY.....,o....o..o.....aooo.....o o o...o.. SUMMARYo..oo o o oo o 0..~o e.. ~.. ~ o o ~ e e * o REFERENCES.... oo........................................ Page iii 1 1 2 7 10 14 15 24 25 ii

LIST OF FIGURES Figure Page 1o Power Spectrum of True and Effective Inputs. o...... 5 2. Closed Loop Representation of an nth Order Nonlinear Servomechanism.................. 8 3. Open Loop Representation of an nth Order Nonlinear Servomechanism....o..o.,,.....O 9 4. Functional Diagram of Amplitude Distribution Analyzer o.o o o. o. o....o o o o o.. o. o o o.. o o... 1.6 5o Error Probability Distribution of a Second Order Servomechanism with Saturation...o................. 18 6. Error Rate Probability Distribution of a Second Order Servomechanism with Saturation.... o....... 19 7. Error Probability Distribution of a Second Order Relay Servomechanism,.......o.......... o..o... 21 8. Error Probability Distributions of a Second Order Relay Servomechanism Plotted on Gaussian Paper..... 22 9. Error Rate Probability Distributions of a Second Order Relay Servomechanism.................*..o.. 23 iii

I, INTRODUCTION The response of a nonlinear system to a random input usually cannot be obtained using the phase space or the describing function methods of analysis, since these methods are applicable only to a nonlinear system which contains only prescribed functions for its input. On the other hand, when a system contains a random input, this input can be characterized only by a set of statistical properties. To overcome this difficulty, Booton first introduced a linearization technique(l), which was based upon the assumption that the probability distribution of the response of a nonlinear closed-loop system containing a zero memory nonlinearity when excited by a Gaussian random input, is approximately Gaussian. This statistical linearization method has been checked experimentally and the reqsults agree with the theoretical calculations within the limits of experimental measurement. This paper determines the theoretical probability distributions of Booton's type of nonlinear system when some special conditions are imposed on the random input function. It simultaneously explains, therefore, the validity of Booton's assumption under these special conditions. II, STATISTICAL IDEA FROM PHASE SPACE POINT OF VIEW It is well known by most engineers that the response of any dynamical system can be represented by a point in phase space. For a second order system the phase space reduces to a phase plane. In this plane a point represents the position and the velocity of the dynamical system at a certain instant of time. With increasing time the representative point, which is continuously moving in this plane, describes a phase trajectory. Its path is analogous to the stream line of flow in

hydrodynamics, and it represents the history of the dynamical system. If the input to the system is a random function of time which has a wide uniform power spectrum, then under this condition the representative point of the system, when described in the phase space, undergoes a random motion. This random motion if observed on an oscilloscope reminds one of the Brownian movement of a particle on the surface of liquid. It is this analogy that led the authors to consider a Markoff random process technique for this motion. It is interesting to point out that the analogy between the phase space flow and the nonlinear differential equation was first considered by Kaplano (2) The present idea of the random motion of a representative point in phase space could be considered as an extension of the ideas of phase space flow. A moment's reflection would have revealed that both of these ideas could have been derived from Boltzmann's equation (3) III. CONDITIONS IMPOSED ON THE FUNCTIONS AND EFFECT OF FEEDBACK In order to apply the Markoff random process technique to a dynamical system, some conditions must be imposed on the random input, F(t). These conditions(4) are: <F(t)>AV = (1) <F(t) F(t + T)>AV = 2D6(T) (2) F(t) is Gaussian distributed (3) Condition (2) implies that the input function F(t) has a white spectrum. By imposing conditions (1) and (2) on the random input function, one might assume that the Markoff random process technique is not useful

in analyzing a feedback control system. This impression is based upon the fact that the power spectrum of the random input to a feedback control system usually is limited to an extremely low frequency region, while condition (2) requires the random input to have a white power spectrum. However, in the following discussion it is shown that this impression is not true. In reality the feedback action of the system creates an effective input, which is a linear combination of the true random input and its higher order derivatives, rather than a true random input., Consider, for example, the first order servomechanism which contains an integrator and a nonlinear element K(E) in its forward path. The following two equations define the dynamic behavior of the system. E = i -oo, (4) and dao t = K(E), (5) where 0o, 0i and e are respectively the output, input and error Substitution of Equation (4) into Equation (5) results in an equation relating the error, E, and input, Qio aE + K(C) = iT = F(t) (6) dt dt The random function F(t) is called an effective random input and is assumed to have the properties prescribed by Equations (1), (2), and (3). Since F(t) is a linear function of the true random input Qi(t) which is a Gaussian random function, it follows immediately that the conditions (1) and (3) are satisfied by F(t). According to power spectrum analysis, the power spectrum qpF(O) of the effective random input

-4 - F(t) and the power spectrum q) i(W) of the true random input ~i(t) are related by F(W)) = |j 12 Q (n). (7) If the power spectrum cpq (c) is of the following form: cpn (W) (8) pi(m) = (2 + 2)(D2 + 2) where a, P, D are constants, then from Equation (7) the power spectrum <pF(c) of the effective random input, F(t), will have the following form: a)2D q (a) == (9) 'F = (a + 4 c)(a2 + o2) If in the frequency range of inter'set the value of the constant, U2, in the expression (9) is very much smaller than o2, and if the value of the constant p2 is very much larger than u2, then under these conditions the power spectrum q)F() can be approximated by the following equations: cpF(c) )-.2o2 for 2 < 2 (lO.a) D) for c2< C2 < ~2 (lO0b) F P pF() 2 for d2 > 2 (10.c) In Figure 1 will be found the power spectra cpg (.) and F(w) as function of co Inspection of Figure 1 reveals the power spectrum of the effective random input, F(t), to be almost flat over the frequency range a < w < P. If the frequency range a < D < P is much larger than the system bandwidth, and, if in addition, the value a is very small, then the

-5 - TRUE INPUT POWER SPECTRUM /F (' ) EFFECTIVE INPUT / POWER SPECTRUM 0 11 w Figure 1. Power Spectrum of True and Effective Inputs

power spectrum qF(c() of the effective random input can be considered as a white spectrum,* Under these conditions the effective random input has approximately the properties described by Equation (2). The Markoff random process technique can be applied therefore to the first order system provided that the random input Qi has the power spectrum qi (co) of the form described by Equation (9)0 In general, if the error of a servomechanism or of a feedback system can be described by the following equation: N n N dn9i an d + K(C) a. n (11) n=0 dtn n tn where an's are constants, as previously defined, then the right-hand side of Equation (11) can be lumped together to form an effective input function F(t) which is given by the following equation: N dngi F(t) = an i (12) n=O dt Accordingly, the power spectrum cpF(w) of the effective input F(t) and the power spectrum cp. i() of the true input Qi(t) are related as follows: N nn=O ~ a.nt2 (~ (13) F(Cp ) = | Z (jo)n an qQ (13) n=O ~ By properly choosing the form of the power spectrum 'cpQ (), it * In an actual case, if a power spectrum is uniform over a frequency range which is greater than the frequency range of the system, then the power spectrum can be considered as a white spectrum.

-7 - is always possible to obtain a power spectrum qF(() which is approximately a white spectrum from a physical point of view. The block diagram representation of Equation (11) can be given in two equivalent forms~ (I) As a closed-loop representation as shown in Figure 2; (II) as an open-loop representation as shown in Figure 3. From Figure 3 it is readily seen that the effective input to the system is essentially white noise provided that the input power spectrum is restricted to a special form determined by Equation (13). Because of the above power spectrum1s changing property in feedback system, it is possible to use Markoff random process technique to analyze several commonly encountered nonlinear servomechanisms. IVo MARKOFF RANDOM PROCESS AND FOKKER-PLANCK EQUATION A discrete Markoff process can be interpreted as a process in which the occurrence of an event depends only on the occurrence of an event immediately preceding ito For a continuous case it is usually defined in terms of its conditional probability density function by an important relation given by Equation (14) which is known as the ChapmannKolmogoroff equation f(t2i Y1! Y2' ~~ Yn I S; S1 2 ooo ~n) = f(t1; x x, ~~~ Xn | S; t; 2, oo~j2 n) (14) f(t2> Yln Y2, OOd Ynl tl; X1, x2 Xn) dxl, dx2 o dxn The function f(t; x1, x oo xn IS; 1, 2, n) stands for the probability density function which states mathematically that at time t1 the random variables are at x1, x2, o.. xn on the assumption that at time S they were at 1, 2, ooo ~. no For simplification in writing,

-8 - Figure 2. Closed Loop Representation of an nth Order Nonlinear Servome chanism

-9 - F f 0 W - --- I,, i (t) N dn n-o dtn F(t) N F(t):= X-O dn+K(E) a dtK,, Figure 3. Open Loop Representation of an nth Order Nonlinear Servome chanism

vector notation is often used, Equation (14) when *written in vector form becomes f(t2 yI |S. t = f f(tL; x 5S; ) f(t2; yx tl; dIt (15) Based upon some conditions imposed on the conditional moments of those random variables, Equation (15) can be reduced to a linear par"tial differential equation known as Fokker-Planck equation(4) ^ - - [l f] * 2 Z [Pijf] (16) i- Yi i3,j where Pi At At f (Ayi) f(t + At, y + y | j d (A Dij = O lim I (AY)(Ayj)f(t + At; y- + y It yI d(A4y d (Ay) ijk. o o = 0 for any B whose subscripts have three or more letters. V o A SECOND ORDER NONLINEAR SERVO SYSTEM AND ITS RELATION WITH FOKKER-PLANCK EQUATION The dynamical behavior of a second order servo system can be described by the following second-order differential equation~ d2o d~o. 0 + a + b - K(c) (17) dt dt = ( where K(c) is a function e, a and b are constant, and ~i and ~o are as previously definedo By substituting e = QQi Equation (17) can be written into the following form: d2 d d2i di at b~ + K(C) a (8) dt2 dt dt2 dt i

If the left side of' Equation (l8) is -written as a single function F(t), Equation (i8) will reduce to the following forms d + a + be + K(c) = F(t) (19) dt Since the function F(t) is a linear combination of a Gaussian input function Qi, it follows that F(t) is a Gaussian random functiono Furthermore by properly choosing the power spectrum of Qi as sho'wn in Section III, F(t) can be approximately treated as Gaussian white noise, and the Markoff random process technique will apply under this condition, If we now let de y1, Ec y2 then the above equation becomes dt two first order simultaneous differential equations in phase space y1 - y2~ dy2 d~t-'+ ay1 + by2 + K(Y2) = F(t) (20ca) dt Yl y1 = ^ (20.b) Let t be increased to t + At where At is very small compared with the response time of the system, but is not necessarily small compared with the function F(t)o Putting the above two simultaneous equations into the incremental form they becomeo Y1 + aylAt + [by2 + K(y2)]jAt - f F(t) dt (21.a) t and AY2 - YlAt (21.b) where J F(t) dt indicates the fact that F(t), in this small time intert val At, is a rapidly changing function of timeO By using Equations (21oa) and (2l.b), it is possible to calculate the functions iij, i. jk. o' io ij.om of the Fokker-Planck equation given in Section IV.

-12 - For the case of the second order servomechanism expressed by Equation (l9), the Fokker-Planck equation has the following form~ 62f yf) (y1f) + {[ayl + by2 + K(y2)] f} + D 22 (22) Here we use the assumption that <F(t) F(t + Tr)>AV = 2DD(T) where 6(T) is a delta function and D a constanto The solution of Equation (22) gives the conditional probability density function of the nonlinear servo system. This solution, in general, is a function of yl1 y2, t and initial conditions y10 Y20o to0 However, as t -t>o the solution becomes independent of Y10, yY20 and to. This is what one would expect, since any random variable after a long time interval becomes independent of its initial conditionso Consequently, as t ->oo the conditional probability function becomes a first probability function. (5) Since we, are interested in the first probability density function, this is equivalent to finding the steady state solution or having 0 in Equation (22)o The Fokker-Planck equation (22) reduces, in this case, to the following form. 0 - (f) + [ay + by1 + K(y2) f} + D a (23) dY1 The additional conditions that the stationary solution must satisfy are listed below~ lo f(y2, y, t = )- 0 as Iyl 1-o or as lY2 1 o or both 'y11 and Y21 —> o 20 f(Y2,Yl, t = 00) should always be a positive function 3 f+ fJ f(Y2, y1 t = oo) dy1 y2 1 o00 ^*o00 2

If we assume a stationary solution of the form f[my12 + n f [by2+K(y2) dy2 ] and use the above three conditions it can be shown that -a 2 a f(Y2,y1 t 0) = A exp [D y1 - f [bY2 + K(Y2)J dy2J (24) where A is determined by the normalized condition 30 The steady-state solution or the first probability density function in phase-space is the most general type solution that can be obtained from a second order system 'which has an arbitrary zero memory nonlinearity ip its forward path, when it has been subjected to a Gaussian random inputo From this solution it is readily seen that in spite of a nonlinearity in this type of second order system- the rate of the response of the system, namely y1, still retains the Gaussian propertyo However, the response y2 itself is changed to a random process that is of a nonGaussian type For the case when the nonlinear function K(y2) is approximately linear in the region around y2, such as in the case of saturated amplifier approximated by K(y2) = tanh (cy2), the first probability function f(Y2, yl' t =.0) can almost be considered as approximating a Gaussian probability density functiono This can be seen from the expressiono a2 2 f(y2, y!, t = ) = A exp [j2D- YA 1 exp L[ J tanh (cy2) dy2] -a2 2 -a A exp [-D Y12 exp [ — ln cosh (cy2)] when ley2 < 1 in cosh (cy2))(cy:2)2 (. then a2 2a f(Y2, Yi, t = ) - A exp [2- Y1 ] exp [ (cy2)21 which is evidently the same as a Gaussian probability density functiono

VIo HIGHER-ORDER SYSTEM The above technique which was applied to a class of second order systems in Section V can be extended to higher-order systemso However, in the case of higher order systems, only under very special conditions could the steady state solution be obtained by the authors. To be specific, only those systems that can be decomposed into the following form can the steady state solution of the Fokker-Planck equation be obtainedo i dx N 2 + ai d- + Ki(xi) + Z g1j(x) - F(t) i 1, 2, o N (25) dt2 t j 1 J where ai's are constants, gij (xj) and gji (xj) are of the same functional form, and F (t) s are Gaussian white noise. The first probability density function usually has the following form~ N 2 N N N f = A exp biVi + Z C1J K (xi)dxi + Z d i l gj (xj x (i ivl i=1 ii +1 i Z dxi where V =- 't and bi s, Ci s, d is are constants. The constant, A, is the normalization factoro It should be pointed out that when a higher-order system is decomposed into a system of equations (25) then it is no longer of the form of Equation (11), namelyN dnQ0 N dng Z an tn + K(c) = Zan tn n=O n=O Consequently, the argument used in Section IV is no longer true, because in this case one cannot write the effective input, F(t), as a linear combination of the true input Q,. However, if the input is nonGassn dtr ted an te oeral e t of the nonnear combin on Gaussian. distributed. and. the overall effect of the nonlinear combination

-15 - of the true input Oi still gives an effective input, F(t) which is Gaussian distributed and has white spectrum, then it is possible to obtain the first probability density function of this system. VIIo COMPUTER STUDY In order to verify the result of the analysis presented in the previous sections, an analog computer study was made on a second order nonlinear servomechanism. Before the results of the computer study are discussed, the measuring equipment or the single channel amplitude distribution analyzer will be described. The functional diagram of the analyzer is given in Figure 4o This analyzer measures the probability distribution function P(x) of a random variable x, the relation between the probability distribution function P-(x) and the probability density function f(x) which is of the form: x P(x) = f f(x) dx. (26) -00 The resolution of the analyzer was checked against a known function of time such as a triangular wave. It can be shown the probability distribution function P(x) vs. the amplitude, x, of the triangular wave theoretically is a straight lineo The agreement between the theoretical straight line and the experimental points was found to be excellento The noise generator which was used, was constructed in the Servomechanism 1 Laboratory. It had a flat spectrum from approximately 375 cps to 25 cps and was Gaussian distributed. The equation of the second order servomechanisms under study is of the form: d2E d_ d 2i dQ d t — + — F(t) (27) dt2 dt dt2 dt

INPUT RANDOM SIGNAL VARIABLE TRIGGERING LEVEL J 0 t I 7L FL RANDOM INPUT I I VARIABLE EXPANDING D.C. AMPLIFIER VOLTAGE SCHMITT CIRCUIT I I I I Ill II 118111111 11 1111111 Ill~ll~fl~l 1111111l I I I I L.1I.U.LL LL1.LJ.1L I.....L...LLLI PULSE GEN. I I 4 r I= I L L Ir A I I I I!L I I II Li1 _.|. l............................. SAMPLING RATE = PULSE - REPETITION RATE TIMER =S PULSE COUNTER Figure 4. Functional Diagram of Amplitude Distribution Analyzer

-17 - Two different kinds of function K(E) were used, these are listed below. For the case of the saturated amplifiers 100 e > 1 (1) K(e) = o100e | I < 1 and l = 10 (28) -100 E < 1 For the case of the relay amplifier: ( K, e>0) (2) K(c) = where K = 100 and D = 20. (29) -K C< o0 In the first case, the measurements were intended to show the departure of the distribution of the error, E, from normality in the nonlinear region; therefore, under these conditions only, is the test of normality of the error important~ In order to confirm the belief that the departure from normality was not due to statistical fluctuation of the data, a confidence interval test was performed. The results of the experimental data were plotted on normal probability papers and are given in Figure 5 and Figure 6. In the second case, a detailed study was made because of the nonlinearity occurring for small values of error as well as for large values of error. It presents, therefore, a severe test of the validity of the analysis given in previous sections. For Equation (27) with K(c) and P as given in (29), the first probability density function f(c e), according to Equation (24), is of the form: f(e, e) = A exp [2 2 + ] for e<O (30.a) - A exp [2 D ] for (30b)

-18 - or 0 z n~ b. 0 - CD a. a -=I -T- - - -- - - - -- - -- -99 95% LIMITS OF THE RANDOM VARIATION 98 2 98 ABOUT THE THEORETICAL STRAIGHT LINE 95 10 90 30 -- -70 50 50 70 ---— 30 90- I 10 j95_- LINEAR5 / REGION 98 — 2 99 -~- I -2 -1 1 2 e, ERROR IN VOLTS I Figure 5. Error Probability Distribution Servomechanism with Saturation of a Second Order

-19 - LJ cr 0 ncr w LU. 0 z 0 Icn ICE a(U 0 nQI 4 5 --- 10 -- - 30 50 -70 / 95 98 -- -- 99 I 0 J- 0 n A B A 1a I rr 0 0% A 30 98 90 70 50 30 10 5 2> -'Z -1U -Il -I -4 - U 4 U i i, ERROR RATE IN VOLTS 1b zU z4 Figure 6. Error Rate Probability Distribution of a Second Order Servomechanism with Saturation

where A is a constant. By using the normalization condition -+00 +00 L SL0 f(e e) de dc =, it is possible to determine the constant Ao The first probability density function of E alone is determined from the relationf(c) -= f(,c)e de - exp [D c] for e0 (31.a) -00 2D = 2 exp [ e c] for c>0 (31.b) The probability distribution function P(c), which is defined below, is related to the probability density function f(c) by the following relations: P(E) f f(t) dt = 2 exp [D c] for E< 0 (32.a) -00 P(C) = 1 - 2 exp [ c] for e>0 (32ob) Since the functions P(e) and 1- P(c) are equal to some exponential functions of c, it is evident that the graph of the data plotted on semi-logarithmic paper should be a straight lineo Furthermore, if e =- is substituted in Equation (32.a), P(c) is found to equal 0.1834. PK Based upon this fact, by finding the point on the graph of P(c) vs.o at which the function P(c) is equal to 18034 per cent, it is possible to D determine the value of E which corresponds to the quantity PK. The results of the experimental data were plotted on semi-logarithmic paper and normal probability paper. They are given in Figure 7, Figure 8 and Figure 9. The effects of the parameters on the probability distribution were checked by comparing the change with an arbitrarily chosen standard system. The theoretical calculation and the measured data are listed below, the results were amazingly close:

-21 - or 0 UJ,L. 0 z 0 In I-J Cn 0 0 (0 aw a. (3 z 0" 70 60 50 40 C,,, // \\ \ X 30 20O I 9 l/ D |D/0 2 --— _ --- ~I_11 -5 -4 -3 -2 -1 0 ERROR IN VOLTSe~ 2 3 5 e, Figure 7. Error Probability Distribution of a Relay Servomechanism Second Order

1 __._ __ __ I 1 TT1 1 F 7 I I I I I I1 I 0 I02 02 DD 02 0 0r I'ar 5 --------- --- -- --- 95 100. I -------- 90 so~~~~~~~~~~~~~~~~~- I I I 30 — 0 7- 70 ~ _8... — ~. --- —- - --— ~- —... --- -2 950-. / 95L LL1 I I ~-1 R 0 1 e, ERROR IN VOLTS 2 3 4 Figure 8. Error Probability Distributions of a Second Order Relay Servomechanism Plotted on Gaussian Paper

-23 - w I0 LrJ w U0 z 0 aa I(I) -J 0 cr -a. e, ERROR RATE IN VOLTS Figure 9. Error Rate Probability Distributions of a Second Order Relay Servomechanism

Curve Changed Parameter B DN = Ds C ~M2sPKN = 2 Ps D KN = 1.5 Ks Theoretical Calculation Ds ~ = 2.25 Ps - - 2~0 PN Ks - o 666 KN Measured Data Ds = 2.22 DN s PSN = 2.0 Ks -K = 0.690 KN Remarks DN = changed effective input power spectrum PN = changed damping KN - changed force A STANDARD SYSTEM D5 = effective input power spectrum of the standard system s = damping of the standard system Ks restarting force of the standard system SUMMARY The purpose of this paper is to present a method of determining the first probability density function of a nonlinear system subjected to a Gaussian random input~ The Markoff random process technique was used. The results from the computer study agree with the theoretical calculations within the limits of experimental erroro Although this study is limited to a Gaussian input with some definite power spectrum, this analysis can be applied to a large class of nonlinear systems

REFERENCES 1. Booton, Ro C., Proc. of the Symposium on Nonlinear Circuit Analysis, Polytechnic Instit. of Brooklyn, Volo 2, 1953o 2. Kaplan, Wo, Proc. of the Symposium on Nonlinear Circuit Analysis-, Polytechnic Instit. of Brooklyn, Vol. 2, 19530 3o Wang, Mo C., PhoD Thesis The University of Michigan. 4. Wang) Mo C. and Uhlenbeck, G. Eo, Review of Modern Physics, Volo 17, Noo 2 and No. 3, pp. 323-342, 1945o 50 Kolmogoroff, Ao, Mathematische Annalen, Volo 104, p. 456, 1931.

UNIVERSITY OF MICHIGAN 3 9015 02829 52211111li 3 9015 02829 5221