Response of Duffing Oscillator with Time Delay Subjected to Combined Harmonic and Random Excitations

This paper aims to investigate the stationary probability density functions of the Duffing oscillator with time delay subjected to combined harmonic and white noise excitation by the method of stochastic averaging and equivalent linearization. By the transformation based on the fundamental matrix of the degenerate Duffing system, the paper shows that the displacement and the velocity with time delay in the Duffing oscillator can be computed approximately in non-time delay terms. Hence, the stochastic system with time delay is transformed into the corresponding stochastic non-time delay equation in Ito sense. The approximate stationary probability density function of the original system can be found by combining the stochastic averaging method, the equivalent linearization method, and the technique of auxiliary function. The response of Duffing oscillator is investigated. The analytical results are verified by numerical simulation results.


Introduction
It is known that time delay in real active control systems is unavoidable due to the time spent in calculating and executing the control forces, performing online computation, and so on.Hence, time delay causes unsynchronized application of a control force and often leads to poor performance or instability of controlled systems [1].In the theory of random vibration, much attention has been paid to the study of the oscillator system with time delay under random excitation extensively.For the case of systems with time delay excited by Gaussian white noise, some methods/techniques have been developed to treat these such systems [1][2][3][4][5][6].For instance, in [1], Di Paola and Pirrotta used the Taylor expansion of the control force and another approach to finding exact stationary solution to study the effects of time delay on the controlled linear systems.In 2009, Li et al. [2] studied effects of time delay in feedback control on the first-passage failure of controlled systems under stochastic excitation by using the stochastic averaging method for quasi-integrable Hamiltonian systems.In 2012, Liu and Zhu [3] proposed a procedure based on the stochastic averaging method for the time delay stochastic optimal control and stabilization of quasi-integrable Hamiltonian systems subject to Gaussian white noise excitation to study the response and stability of systems.They converted the problem of time delay stochastic optimal control of quasi-integrable Hamiltonian systems into the problem of stochastic optimal control without time delay and the result problem is solved by applying the stochastic averaging method for quasi-integrable Hamiltonian systems and the stochastic dynamical programming principle.In Feng and Liu's work [6], they used stochastic averaging method with assumptions that state variables are slowly varying processes and then solved the Fokker-Planck (FP) equation by finite difference method to yield the stationary joint probability density.To authors' knowledge, this assumption is based on no appropriate reason/proof from now on.Maybe the assumption is based on authors' experiences.In our work, we showed that this assumption is acceptable.
In the present work we investigate Duffing oscillator under time delay state feedback, involving both displacement and velocity feedback, and under combined harmonic and random excitation.Through this study, a new approximate procedure for vibration systems with time delay under harmonic and random excitation is proposed.The system is expressed equivalently in terms without time delay.Then the equivalent system is transformed into Ito stochastic equations which are treated by using the stochastic averaging method.The solution of the FP equation associated with the Ito stochastic system is a very difficult problem.To overcome the arising problem, the technique of combining the linear equivalent method and the auxiliary function technique proposed in [7][8][9] is utilized.The paper is organized as follows: in Section 2, the approximate technique for Duffing oscillator with time delay is discussed.In Section 3, effects of system's parameters on the response are investigated and compared with Monte-Carlo simulation (MCS).Finally, in Section 4, summary and conclusions are given.

Duffing with Time Delay Subjected to Combined Periodic and Random Excitations
Let us consider the Duffing system subjected to a harmonic force  cos ] and the white noise (): where , ℎ, , , ], ,  1 ,  2 ,  > 0 are constants,  is a small positive parameter,   = ( − ), ẋ  = ẋ ( − ), dot denotes differentiation with respect to time, and function () is a Gaussian white noise process of unit intensity with the correlation function   () = [()( + )] = (), where () is the Dirac delta function, and notation (⋅) denotes the mathematical expectation operator.Consider the primary resonant domain where there is a relation between  and ] as where Δ is a detuning parameter.Substituting (2) into (1) one obtains To the authors' knowledge, the stochastic averaging technique is developed for non-time delay systems which are in standard form and derivatives of their variables are proportional to the small parameter.There is no version for systems with the time delay.Thus, so as to apply this technique, the system with time delay should be replaced appropriately by the non-time delay one.
In order to employ the theory of stochastic differential equations, in many works, time delay terms were approximated by non-time delay ones under some assumptions such as "small time delay," for example, Bilello et al. (2002) [10], Feng et al. (2009) [5], Feng et al. (2011) [11], and Liu and Zhu (2012) [3], or "small parameter," for example, Li et al. ( 2009) [2] and Li and Feng (2010) [4].Under such assumptions, the amplitude () of the system's response can be treated as a slow varying term.Thus, one obtains the approximation: ( − ) ≈ ().Then, the time delay terms are approximated by non-time delay ones.Maybe these assumptions are based on the authors' experiences.In the next paragraph, it is showed that these assumptions are acceptable because they come from the fact that the new variables are, as shown in (8) and (10), varying slowly processes since their derivatives are proportional to the small parameter .To do that, (3) is rewritten in matrix form as follows: where Make the following transformation in (4): where is the fundamental matrix of ordinary differential equation (7), obtained from (6) for  = 0: Substituting ( 6) into (4) one obtains where It is observed from (8) that () can be assumed to be a slowly varying random process since  is the small parameter; that is, one has the following approximation: Using ( 6) and (10) with noting that ( − ) = (−)() gives where Substituting ( 11) into (4) gives the equivalent system with (4) without time delay as follows: Rewriting obtained (13), with noting ( 12), in second-order stochastic differential equation yields Let us denote Then ( 14) becomes It is noted that ( 16) is non-time delay one which is an approximation of original system (1).
We seek the solution of ( 16) in the form of where  1 and  2 are random processes satisfying an additional condition Substituting ( 17) into ( 16) and then solving the resulting equation and (18) with respect to the derivatives ȧ 1 and ȧ 1 yield This pair of stochastic differential equations can be simplified by using the stochastic averaging method [9,12,13]: Here  1 () and  2 () are independent white noises with unit intensity, and the drift coefficients  1 ( 1 ,  2 ) and  2 ( 1 ,  2 ) are determined as follows: where ⟨⋅⟩ is a time-averaging operator over one period T defined by From ( 21), one obtains the drift coefficients of system (20): The FP equation written for the stationary PDF ( 1 ,  2 ) associated with system (20) has the form Mathematical Problems in Engineering Solution (24) is still a difficult problem because functions  1 ( 1 ,  2 ) and  2 ( 1 ,  2 ) are nonlinear functions in  1 ,  2 and they do not fulfill the potential condition as showed in [12,14].To overcome this difficulty, the equivalent linearization method [15][16][17] is employed.Following this method, the nonlinear functions  1 ,  2 are replaced by linear ones.Denote According to the stochastic equivalent linearization method, nonlinear terms (25) are replaced by where equivalent coefficients   ,  = 1, 2,  = 1, 2, 3, are to be determined by an optimization criterion.Thus, the functions   ,  = 1, 2, in (23) are replaced by linear functions So far, it is seen that ( 24), where  1 ( 1 ,  2 ) and  2 ( 1 ,  2 ) are replaced by (27), is still hard to solve because it does not fulfill the potential condition [12].In order to integrate (24), where  1 ( 1 ,  2 ) and  2 ( 1 ,  2 ) are replaced by (27), following the technique of auxiliary function [8,9], we introduce an auxiliary function ( 1 ,  2 ) =  0 = const as follows: We will choose the function ( 1 ,  2 ) =  0 so that the equalities below are fulfilled: Substituting (30) into system (29) one obtains Solving system (31) in Φ/ 1 and Φ/ 2 yields where Eliminating Φ( 1 ,  2 ) from system (32) gives the equation for the auxiliary function After finding the function ( 1 ,  2 ) =  0 , the stationary PDF ( 1 ,  2 ) can be found from (30), (32), and (33) by the quadrature where  is a normalization constant [9,18].Hence, the corresponding FP equation to (24), where drift coefficients are the linear functions (27), has the following solution: where  is a normalization constant and coefficients   ,  = 1, 5 are determined as follows: It is noted that the joint PDF ( 1 ,  2 ) determined by (36) has finite integral if coefficients  1 and  2 are positive.This condition is always fulfilled because (24), where drift coefficients are linear functions (27), is associated with a linear system under Gaussian white noise in the form of (20).Therefore, the approximate stationary PDF of ( 24) is determined by (36) whose coefficients are given in (37).It is seen from ( 36) that random variables  1 and  2 are jointly Gaussian [19].Thus, from (36), one obtains where  2  1 and  2  2 are variance of  1 and  2 , respectively, and   1  2 is covariance of  1 and  2 .It is seen from ( 38) that necessary statistics of processes  1 and  2 can be computed algebraically based on coefficients of joint PDF ( 1 ,  2 ).
From transformation (17), the mean response of the oscillator can be rewritten in the form Mathematical Problems in Engineering where tan  = −( 2 )/( 1 ).Thus it is periodic with amplitude , where Also from transformation (17), the mean square response of (1) can be determined as follows: Taking averaging with respect to time (46) gives Substituting ( 38) into (47) and reducing the obtained result yield the time-averaging of mean square response to be where   ,  = 1, 5, are given by (37).It is noted from (48) that the approximate time-averaging value of mean square response of the Duffing oscillator is calculated algebraically.
Moreover, from transformation (17), one obtains It follows from transformation (49) and PDF (36) that the joint PDF of  and ẋ takes the form of sin ]) It is seen from ( 36), (50), and (51) that the joint PDF of  and ẋ and the marginal PDF of  depend on time , although two variables  1 and  2 are described in a stationary joint PDF.In Tables 1 and 2, effects of the delay coefficient  2 and the time delay  on time-averaging of mean square response of the system are considered.The analytical approximate results obtained by the proposed technique are compared to simulation ⟨( In Table 1, the time-averaging mean square response of the system is performed for computation with various values of the parameter  2 while the system parameters are chosen to be  = 0.2,  = 1,  = 1, ℎ = 4,  = 2,  2 = 1, ] = 1.01,  1 = 0.1, and  = 1.It is seen that the time-averaging mean square response of the system slightly increases as the delay coefficient  2 decreases from 0.01 to 0.2.Table 2 presents time-averaging values of the mean square response of the system evaluated versus the parameter  in the primary resonant region with the system parameters chosen to be  = 0.2,  = 1,  = 1, ℎ = 5,  = 2,  2 = 1, ] = 1.01, and  1 =  2 = 0.1.In Table 2, one expects that the error will decrease as the time delay decreases.However, sometimes this does not often happen in analysis of a stochastic system with time delay; for example, see Figure 1 in Feng et al. (2009) [5].In the paper, the Duffing oscillator under studying is subjected to both periodic and random excitation.The random excitation would make the diversity of the system's response unpredictable.On the other hand, the error in the paper consists of 3 errors from replacing the time delay system with the non-time delay one, applying the stochastic averaging method, and applying the equivalent linearization method.Thus, in a particular range of the time delay, the values of error may be fluctuated.However, the errors in Table 2 are quite small (from 0.34% to 3.86%) as the time

Summary and Conclusions
In the present paper, a technique for predicting the response of the Duffing system with time delay under combined harmonic and white noise excitation has been proposed.The technique can be summarized as follows.Firstly, the time delay terms are expressed approximately in terms of the system state variables without time delay and, correspondingly, the original system with time delay is approximately transformed into one without time delay.The paper has shown this transformation in detail by ( 6), ( 7), ( 8), ( 10), ( 11), ( 12), ( 13), ( 14), (15), and (16).Secondly, the state coordinates (, ẋ ) are transformed to Cartesian coordinates ( 1 ,  2 ).In these coordinates, the averaged equations obtained by the stochastic averaging method are nonlinear ones whose solution is still a difficult problem.Thirdly, the stochastic equivalent linearization method and the technique of auxiliary function are employed to solve approximately the Fokker-Planck equation associated with the averaged equations.The linearization coefficients of the equivalent linear system are determined from a closed nonlinear algebraic system as presented.Finally, the response of the system is obtained

Table 1 :
The error between the simulation result and approximate values of the time-averaging of mean square response versus the parameter  2 .

Table 2 :
The error between the simulation result and approximate values of the time-averaging of mean square response versus the parameter .It shows that when the time delay increases from 0.05 to 1, the time-averaging values of the mean square response of the system slightly decrease.The two tables show that the proposed technique gives a good prediction.