Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 2, 40, 2002 OPTIMAL CONTROL OF A DUFFING OSCILLATOR UNDER PARAMETRIC AND EXTERNAL EXCITATIONS Artur Nowoświat Department of Building Processes, Technical University of Silesia, Gliwice e-mail: now2@wp.pl The present paper deals with investigations of a stochastic model of the Duffing oscillator operating continuously. The standard theory of the optimal stochastic control described in Fleming andRichel (1975) is presented, as well as the application of these concept to the perturba- tion control techniques suggested by Suhardio et al. (1992). Themethod presented in this paper is, therefore, a generalization of the method of the development of nonlinear control units into a series, as described for Duffing’s deterministic oscillator involving stochastic systems deri- ved by the author and applied in his master thesis (Nowoświat, 1999). Detailed analysis and numerical calculations have been done using the Runge-Kutta procedures. Key words: perturbation control techniques, Duffing’s oscillator 1. Introduction Theproblemof optimal control of stochastic parametrically and externally excited nonlinear systems was recently studied by Beaman (1984), Yoshida (1984), Young andChang (1988). By the combined use of Gaussian statistical linearization and linear quadraticGaussian (LQG) theory a suboptimal linear Beaman (1984), Yoshida (1984) or nonlinear Young and Chang (1988) state feedback controller is synthesized. At the same time, a perturbation technique was applied to the determination of optimal control in the deterministicmodel of nonlinearDuffing’s oscillator inSpencer et al. (1966), Suhardio et al. (1992). The effect of higher-order feedback corrections based upon series expansions of the optimal cost function and the optimal control function in theHamilton- Jacobi-Bellman approach were shown. The application of the perturbation 436 A.Nowoświat method to the determination of optimal control of a nonlinear discrete-time systemwas given in Shefer and Breakwell (1987). The studyof optimal control of nonlinear stochastic systemshasbeengiven considerable attention in recent years. In involved stochastic systems which are described by the Itô stochastic differential equations with the Gaussian parametric and external noise excitations. The success obtained due to the applied approach indicates that it is advi- sable to include external excitations in the Duffing deterministic oscillator. It has been achieved by adding a model of the external constraint to the classi- cal formulation of Duffing’s oscillator. In effect, it led to analysis of stochastic differential equations known in the literature as the Itô equations. These equ- ations are the most popular mathematical means used when considering this problem. This work contains only theoretical considerations involving exa- mination of dynamic systems. The mentioned considerations, as well as the obtained results can be applied to technology, e.g. while analysing suspension systems of wheeled vehicles as well as other problems related with Duffing’s oscillator. This paper is concerned with the study of a stochastic continuous time model of Duffing’s oscillator. Using the concepts of standard methods of sto- chastic optimal control, see e.g. Fleming and Richel (1975), and combining them with the concept by Suhardio et al. (1992) on the application of a per- turbation technique to the control problems, we derive results on the optimal control ofDuffing’s oscillator underparametric andexternal excitations.Then, a second-order stochastic parametrically and externally excited Duffing’s ty- pe system is selected to illustrate the application of the ’development into a series’ technique to the optimal control of nonlinear stochastic systems and to compare the procedure with the other methods. 2. Perturbation method for a non linear Duffing’s oscillator Consider a nonlinear stochastic model of a dynamic system described by the Itô vector differential equation dx(t)= [f(x, t)+Bu] dt+ M ∑ k=0 Gk dξk (2.1) where x – state vector, x= [x1, ...,xn]⊤ u – controlling vector, u= [u1, ...,un]⊤ Optimal control of a Duffing oscillator... 437 B – matrix with constant elements with the dimensions n× l, B= [Bij] Gk – deterministic vectors, Gk = [G 1 k , ...,Gn k ]⊤ ξk – independent standardWiener processes. Let us assume nonlinear vector functions f : Rn → Rn in the form of polynomials fi(x)= n ∑ j=1 Aijx j + n ∑ j,k=1 Aijkx jk+ ... (2.2) where Aij,A i jk, ... are constant coefficients and x jk =xjxk, and xjkl =xjxkxl. Assume also that the control u takes the form of a polynomial uj =K j 0 + n ∑ k=1 K j k xk+ n ∑ k,l=1 K j kl xkl+ ... (2.3) where K j 0,K j k ,Ki kl are constant coefficients. According to the model presented by Suharido et al. (1992) and Yoshida (1984) we minimize the following cost function I = E [ V0(t1)+Px(t1)+x ⊤(t1)Mx(t1)+ (2.4) + t1 ∫ t0 ( x ⊤ Qx+u⊤Ru+ n ∑ j,k,l=1 Qjklx jkl+ ... ) dt ] where E[·] – moment vector P – constant vector M,Q – symmetric semi-positively determinedmatrices R – positively determinedmatrix Qjkl, ... – selected symmetric tensors. Theelements K j 0,K j k ,K j kl , ... occurring inequation (2.3), canbecalculated bymeans ofBellmann’s equation (FlemingandRichel, 1975),which in our case takes the form ∂V ∂t +min u [L∗(V )+L(x,u, t)] = 0 (2.5) where 438 A.Nowoświat L(x,u, t)=x⊤Qx+u⊤Ru+ n ∑ j,k,l=1 Qjklx jkl+ ... (2.6) V =V0(t)+Px(t1)+x ⊤(t1)Mx(t1)+ t1 ∫ t L(x,u, t) dt whereas L is an operator, which, according to (2.6)1, can be defined in the following way L∗(·)= ∂(·) ∂t + n ∑ i=1 ( fi+ m ∑ j=1 Biju j )∂(·) ∂xi + 1 2 M ∑ k=1 ( Gk, ∂ ∂x )2 (·) (2.7) or L∗(·) = ∂(·) ∂t + n ∑ i=1 ( fi+ m ∑ j=1 Biju j )∂(·) ∂xi + 1 2 M ∑ i,j=1 ∂2(·) ∂xi∂xj aij (2.8) where a= [aij] = M ∑ k=0 GkG ⊤ k The function V will be considered in the form ofMaclaurin’s series V =V0(t)+ n ∑ i=1 Vi(t)x i+ n ∑ i,j=1 Vij(t)x ij + n ∑ i,j,k=1 Vijk(t)x ijk+ ... (2.9) where Vij, Vijk, ... are tensors symmetric with respect to their indices. In addition, we consider the quadratic terms in the cost function given in equation (2.4), that is Qijk = ... = 0. Furthermore, we set A i jk = Ai jklm = ...=0 and Vi =Vijk = ...=0. Similarly, the elements of the equal order of control (2.3) will be equated to zero Ki0 =K i jk = ...=0 (2.10) The remaining coefficients of Lapunov’s function (2.10) and of control function (2.3) can be found by solving the equations for Vij,Kji, Vijk,Kjkli, ... min u {( V̇0+ n ∑ i,j=1 V̇ijx ij + n ∑ i,j,k,l=1 V̇ijklx ijkl+ ... ) + Optimal control of a Duffing oscillator... 439 + n ∑ i=1 [( 2 n ∑ j=1 Vijx j +4 n ∑ j,k,l=1 Vijklx jkl+ ... ) · · ( n ∑ j=1 Aijx j + n ∑ j,k,l=1 Aijklx jkl+ ...+ m ∑ j=1 Biju j )] + (2.11) + 1 2 n ∑ i,j=1 ( 2Vij +12 n ∑ k,l=1 Vijklx kl+ ... ) aij + + m ∑ i,j=1 uiRiju j + n ∑ i,j=1 xiQijx j + n ∑ i,j,k,l=1 Qijklx ijkl } =0 Taking the differential coefficients an the left-hand side of equation (2.11) in relation to u, and taking into account that the matrix R is symmetric, we get 2 n ∑ i,j=1 Vijx jBip+4 n ∑ i,j,k,l=1 Vijklx jklBip+ ...+2 n ∑ j=1 Rpju j =0 (2.12) for p=1,2, ...,n. Substituting equation (2.3) into equations (2.11) and (2.12) we get an equation for ui, this control being the optimal one ui∗, with quantities V and x corresponding to the quantities V ∗ and x∗ for the control u∗ (for convenience the asterisk will be omitted in the following study). Thus, we get {( V̇0+ n ∑ i,j=1 V̇ijx ij + n ∑ i,j,k,l=1 V̇ijklx ijkl+ ... ) + + n ∑ i=1 [( 2 n ∑ j=1 Vijx j +4 n ∑ j,k,l=1 Vijklx jkl+ ... ) · · ( n ∑ j=1 Aijx j + n ∑ j,k,l=1 Aijklx jkl+ ...+ + m ∑ j=1 Bij ( n ∑ l=1 K j l xl+ n ∑ l,p,q=1 K j lpq xlpq+ ... ))] + (2.13) + 1 2 n ∑ i,j=1 ( 2Vij +12 n ∑ k,l=1 Vijklx kl+ ... ) aij+ 440 A.Nowoświat + n ∑ i,j=1 ( n ∑ l=1 K j l xl+ n ∑ l,p,q=1 K j lpq xlpq+ ... ) Rij · · ( n ∑ l=1 K j l xl+ n ∑ l,p,q=1 K j lpq xlpq+ ... ) + + n ∑ i,j=1 Qijx ij + n ∑ i,j,k,l=1 Qijklx ijkl+ ... } =0 Substituting equation (2.3) into (2.12) we have 2 n ∑ i,j=1 Vijx jBip+4 n ∑ i,j,k,l=1 Vijklx jklBip+ ...+ (2.14) +2 n ∑ j=1 Rpj ( n ∑ l=1 Kilx l+ n ∑ l,p,q=1 Kilpqx lpq+ ... ) =0 for p=1,2, ...,n. When we take into account the symmetric operator sym[·] and when we repeat such a procedure, we can obtain the required conditions to ensure the optimality in equation (2.4).These conditions, expressedbya set of differential equations concerning thecoefficients Vij,Vijk, ... are tobe foundbysubsequent equating of expressions of the zero order, the second order, the fourth order and so on, and by taking into account the fact that xji =xij V̇0+ n ∑ i,j=1 V̇ijaij =0 sym [ V̇ij +2 n ∑ l=1 VliA l j +6 n ∑ k,l=1 Vijklakl+Qij + n ∑ k,l=1 RklK k iK l j + +2 n ∑ k,l=1 VklB k lK l j ] =0 sym [ V̇ijkl+2 n ∑ p=1 VpiA p jkl +4 n ∑ p=1 Vpijk ( A p l + n ∑ q=1 BpqK q l ) + (2.15) +15 n ∑ p,q=1 Vpqijklapq+ n ∑ p,q=1 K p ijRpqK q kl +Qijkl ] =0 Optimal control of a Duffing oscillator... 441 sym [ V̇ij +2 n ∑ l=1 VliA l j +6 n ∑ k,l=1 Vijklakl+Qij + + n ∑ k,l=1 Rkl ( n ∑ p1,q1=1 R−1 kp1 Vq1iB q1 p1 )( n ∑ p2,q2=1 R−1 lp2 Vq2iB q2 p2 ) − −2 n ∑ k,l=1 VklB k l ( n ∑ p,r=1 R−1 lp VrjB r p )] =0 3. Duffing’s oscillator We consider a stochastic model of Duffing’s oscillator described by the Itô vector stochastic differential equation dx1 = x2 dt (3.1) dx2 = (−ω 2 0x1−2hx2−εx 3 1− bu) dt+ + (−ω20x1−εx 3 1)δ1 dξ1+ δ2x2dξ2+ δ0dξ0 where ω20, h, ε, b and σ are constant parameters, ξ1, ξ2, ξ3 are mutually independent standardWiener processes. Assuming that the control is a state feedback in the form of odd order polynomials i.e. u= ∑ q=1,3,5,... q ∑ i=0 ( q i ) kq−i ix q−i 1 x i 2 (3.2) where kij are constant coefficients which are designed tominimize the quadra- tic cost function I = lim T→∞ 1 T E [ ∞ ∫ 0 ( q1x 2 1(s)+q2x 2 2(s)+ru 2(s) ) ds ] (3.3) where q1, q2 and r are weight coefficients. 442 A.Nowoświat 3.1. Main results To obtain the control force in equation (3.2), which minimizes the cost function in equation (3.3) constrained by the system of equations (3.1), one has to solve the following Bellman’s equation (Fleming and Richel, 1975) ∂V ∂t +min u [L∗(V )+L(x,u)] = 0 (3.4) where V =V (x, t) is the Lyapunov function proposed in the form V (x, t)=α00+ ∑ p=2,4,6,... p ∑ i=0 ( p i ) αp−i ix p−i 1 x i 2 (3.5) L(x, t)= q1x 2 1+q2x 2 2+ru 2 L∗(·) is an infinitesimal operator defined by L∗(·) = x2 ∂(·) ∂x1 +(−ω20x1−2hx2−εx 3 1− bu) ∂(·) ∂x2 + (3.6) + 1 2 [ δ20 + δ 2 1(−ω 2 0x1+εx 3 1) 2+ δ22x 2 2 ]∂2(·) ∂x22 For simplicity, we consider three cases: when n=1 and p=2, n=3 and p=4, n=5 and p=6. Substituting quantities (3.5) and (3.6) into equation (3.4) we minimize the obtained equation with respect to u. Then, we obtain the linear algebraic relationships between coefficients kij and αij k10 = a1α11 k01 = a1α02 k30 =2a1α31 k21 =2a1α22 k12 =2a1α13 k03 =2a1α04 k50 =3a1α51 k41 =3a1α42 k32 =3a1α33 k23 =3a1α24 k14 =3a1α15 k05 =3a1α06 (3.7) where a1 = b/r. Substituting quantities (3.5)-(3.7) into equation ∂V ∂t +L∗(V )+L(x,u)= 0 (3.8) and averaging with respect to the stationary measure generated by the out- put process (solution to equation (3.1)) and equating the coefficients of even Optimal control of a Duffing oscillator... 443 order polynomials to zero we find the nonlinear differential equations for the coefficients αij, i+ j= p, i,j =0,1, ...,p, p=2,4,6, ... α̇ij(t)=Fij(α00,α10, ...,α06,ω 2 0,h,ε,b,σ,q1,q2,r) (3.9) In a general case, when both external and parametric excitations are conside- red i.e. δ0 6=0, δ1 6=0 and δ2 6=0 the stationary solutions to equations (3.9) αij(∞) determine the optimal stationary coefficients kij = kij(∞), i+ j = q, i,j = 0,1, ...,q, q = 1,3,5, ... by relations (3.7). To find the optimal costs we have to substitute the feedback control u defined by (3.3), (3.7) and αij(∞) into the cost function I in equation (3.3). After averaging we find that the cost function I depends on the stationary coefficients kij and stationary mo- ments of the output process x= [x1,x2] ⊤ i.e. the quantities E[x p1 1 x p2 2 ] where p1+p2 =2,4,6, ... . These stationary moments can be found (approximately) from moment equations for the original system (3.1) where u is defined by (3.3), (3.7) and αij(∞) where, for instance, the cumulant closure technique is applied. This can be treated as a general procedure. In a particular case when only linear parametric excitations are considered i.e. δ0 6=0, δ1 =0 and δ2 = 0, then the general procedure of the determination of the optimal con- trol andminimal cost function reduces significantly. It has the form described below. Simplified procedure: Step 1. We find the stationary solutions for αij, i+ j=2, i,j=0,1,2 from the system of equations (3.9) which is closed. Step 2. From equalities (3.7) we find the coefficients kij, i+j =1, i,j=0,1 and we substitute them into the moment equations. Step 3. We find the stationary moments E[xixj], i+ j=2, i,j =0,1,2. Step 4. Using the stationary solutions αij, kij and E[xixj] obtained in steps 1 to 3, similarly we find the stationary solutions for αij, i+ j = 4, i,j = 0,1,2,3,4 from system (3.9), kij, i+ j = 3, i,j = 0,1,2,3 from system (3.7) and E[xixjxkxl], i+j+l+k=4, i,j, l,k=0,1,2,3,4 from the moment equations. Step 5. Using the stationary solutions αij, kij,E[xixj] and E[xixjxkxl] ob- tained in steps 1 to 4, similarly we find the stationary solutions for αij, i+j=6, i,j =0,1, ...,6 from system (3.9), kij, i+j=5, i,j=0,1, ...,5 from system (3.7) and E[xixjxlxkxpxq], i+ j + l + k + p+ q = 6, i,j, l,k,p,q =0,1, ...,6 from themoment equations. 444 A.Nowoświat Step 6. We calculate the stationary value of criterion (3.3) i.e. the integral given by (2.4). The advantage of this procedure, when compared with the general one, lies is in finding the stationary solutions to system of equations (3.9) and moment equations. In the general procedure the- se equations are coupled and have to be solved together while in the simplified procedure they can be solved separately. We illustrate this approach by the following example. 3.2. Example We consider the Duffing oscillator for the following parameters ω20 = 1, h=0.01, δ0 6=0, δ1 =0, δ2 =0, q1 =100, q2 =10, r=1. The plot of I versus parameters t, ε and δ0 is given in Fig.1, Fig.2 and Fig.3, respectively. The diagrams enable one to compare the first-order and higher order controls. Fig. 1. Cost function I versus coefficient t 3.3. Final conclusions Theminimization of the cost function I, (3.3), has been analysed for three approximations of the control function, viz. approximation of the 1st order (li- near approximation), 3rd order and 5th order. The obtained numerical results Optimal control of a Duffing oscillator... 445 Fig. 2. Cost function I versus coefficient α Fig. 3. Cost function I versus coefficient δ0 446 A.Nowoświat indicate that the approximation of the 5th order is the closest to the solutions of statistical linearization described in Socha (2000). The method of series involving systems with one degree of freedom pre- sented in this paper provides good and interesting results at a rather low expenditure of labour, but in the case of systems with a larger number of degrees of freedom this method proves to be ineffective due to a considerable number of calculations to be done. Yet, the numerical results remain satis- factory when compared with other methods of the statistical linearization, as indicated in the master thesis by Nowoświat (1999). References 1. Beaman J.J., 1984,Nonlinear quadraticGaussian control, Int. J. Control, 39, 343-361 2. FlemingW.H.,RichelR., 1975,Deterministic and Stochastic Optimal Con- trol, Springer, Berlin 3. Nowoświat A., 1999, Zastosowaniemetody szeregóww analizie stochastycz- nych układów dynamicznych, Pracamagisterska nie publikowana, Gliwice 4. Shefer M., Breakwell J.V., 1987, Estimation and control with cubic non- linearities, J. Opt. Theory Appl., 53, 1- 7 5. Spencer B.F. Jr., Timlin T.L., Sain M.K., Dyke S.J., 1996, Series solu- tion of a class of nonlinear optimal regulatores, J. of Optimization Theory and Applications, 91, 321- 345 6. SuhardioJ., SpencerB.F. Jr., SainM.K., 1992,Nonlinear optimal control of a Duffing system, Int. J. Non-Linear Mechanics, 27, 157-172 7. Socha L., 2000,Application of statistical linearization techniques to design of quasi-optimal active control of nonlinear systems, J. of Theoretical andApplied Mechanics, 38, 3, 591-605 8. YoshidaK., 1984,Amethodof optimal control of nonlinear stochastic systems with non-quadratic criteria, Int. J. Control, 39, 279-291 9. YoungG.E.,ChangR.J., 1988,Optimal control of stochastic parametrically and externally excited nolinear control systems, Trans. ASME J. Dyn. Syst. Meas. Cont., 110, 114-119 Optimal control of a Duffing oscillator... 447 Optymalne sterownie parametrycznym i zewnętrznie wzbudzonym oscylatorem Duffinga Streszczenie W pracy zaproponowanometodę perturbacji dla sterowania układów dynamicz- nych z wymuszeniami o charakterze białych szumów Gaussowskich i kryteriami wprzestrzeni funkcji gęstości prawdopodobieństw.Przywyznaczaniuwspółczynników sterowania, jak również współczynników funkcji Lapunowa, korzysta się z równania Bellmana.Takwyznaczonewspółczynniki niezbędne są dominimalizacji funkcji kosz- tów I. Szczegółową analizę i obliczenia numeryczne wykonano za pomocą procedur Rungego-Kutty. Manuscript received April 2, 2001; accepted for print October 30, 2001