Abstract
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.
1. 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–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–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.
2. Duffing with Time Delay Subjected to Combined Periodic and Random Excitations
Let us consider the Duffing system subjected to a harmonic force and the white noise :where , , , , , , , , 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 :
Substituting (6) into (4) one obtainswhere
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 and are random processes satisfying an additional condition
Substituting (17) into (16) and then solving the resulting equation and (18) with respect to the derivatives and yield
This pair of stochastic differential equations can be simplified by using the stochastic averaging method [9, 12, 13]:
Here and are independent white noises with unit intensity, and the drift coefficients and 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 associated with system (20) has the form
Solution of (24) is still a difficult problem because functions and are nonlinear functions in , and they do not fulfill the potential condition as showed in [12, 14]. To overcome this difficulty, the equivalent linearization method [15–17] is employed. Following this method, the nonlinear functions , are replaced by linear ones. Denote
According to the stochastic equivalent linearization method, nonlinear terms (25) are replaced by where equivalent coefficients , , , are to be determined by an optimization criterion. Thus, the functions , , in (23) are replaced by linear functions
So far, it is seen that (24), where and are replaced by (27), is still hard to solve because it does not fulfill the potential condition [12]. In order to integrate (24), where and are replaced by (27), following the technique of auxiliary function [8, 9], we introduce an auxiliary function as follows:We will choose the function so that the equalities below are fulfilled:DenoteSubstituting (30) into system (29) one obtainsSolving system (31) in and yieldswhere Eliminating from system (32) gives the equation for the auxiliary function After finding the function , the stationary PDF can be found from (30), (32), and (33) by the quadraturewhere 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 , are determined as follows:
It is noted that the joint PDF determined by (36) has finite integral if coefficients and 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 and are jointly Gaussian [19]. Thus, from (36), one obtains where and are variance of and , respectively, and is covariance of and . It is seen from (38) that necessary statistics of processes and can be computed algebraically based on coefficients of joint PDF .
Thus, the approximate solution (36) of (1) is completely determined when the linearization coefficients , ; are found. There are some criteria for determining the coefficients , , . The most extensively used criterion is the mean square error criterion which requires that the mean square of the following errors be minimum [17]. From (23), (25), (26), and (27), the errors in this problem will be
So, the mean square error criterion leads to
From it follows that where , are given by (25). Solving system (42) in , , , with noting that higher moments of and can be expressed in the first and second moments because and are jointly Gaussian [16], gives
Thus, , , , are determined from the system of nonlinear equations obtained by combining (37), (38), and (43). After being found, the values of coefficients , , , are substituted into (36) to obtain the approximate stationary PDF in and of (1).
From transformation (17), the mean response of the oscillator can be rewritten in the formwhere . 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 , , 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 From (50), one gets the marginal PDF of asIt is seen from (36), (50), and (51) that the joint PDF of and and the marginal PDF of depend on time , although two variables and are described in a stationary joint PDF.
3. Numerical Results
In the paper, the obtained results are compared to ones obtained from the Monte-Carlo simulation which is “often the sole tool available for assessing the accuracy of random vibration solutions generated by approximate methods of analysis,” as pointed out in Robert and Spanos (2003) [17]. The various values of response of Duffing equation (1) are compared to the numerical simulation results versus the particular parameter. The numerical simulation of the mean square response is obtained by 10,000-realization Monte-Carlo simulation (MCS) with the time from zero to 300 seconds. The time-averaged mean square response of the Duffing oscillator obtained by formula (48) is compared to a numerical result obtained by MCS in Figure 1.

To obtain theoretical results below, algebraic system (43) is solved with initial values , , , , , . With the input parameters , , , , , , , in Figure 1, we present the effect of periodic force’s amplitude on time-averaging mean square response. It may be seen that the theoretical predictions and the simulations compare very well. Moreover, it can be observed from Figure 1 that the time-averaging mean square response increases when harmonic excitation increases and decreases when the damping coefficient increases.
In Figure 2, we present graphs of PDF of x plotting the time (sec). In order to get the simulation of PDF of by MCS, 20,000 sample paths have been used and counted at the time (sec). Figure 2 shows that the proposed technique gives a good prediction.

In Tables 1 and 2, effects of the delay coefficient 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 results by MCS. The errors in the tables are defined as
In Table 1, the time-averaging mean square response of the system is performed for computation with various values of the parameter while the system parameters are chosen to be , , , , , , , , and . It is seen that the time-averaging mean square response of the system slightly increases as the delay coefficient 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 , , , , , , , and . 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 delay in . 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.
4. 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 . 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 from several algebraic and integral expressions (44)–(51). It has been shown that the results obtained by using the proposed procedure agree well with those from the numerical simulation of the original system. For the given parameters, the time-averaging of mean square response of the system is increasing with the periodic force’s amplitude and the delay coefficient of velocity ; and it decreases when the time delay increases from 0.05 to 1.
A significant contribution of the present paper is that the technique gained through it can be helpful for other nonlinear systems with time delay under harmonic and random excitation. Furthermore, the paper shows that the new variables are, as shown in (8), varying slowly processes since their derivatives are proportional to the small parameter . This approach allows simplifying the research by averaging procedure and getting analytical expression for solution.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research is funded by Vietnam National University, HCM City, under Grant no. C2016-26-04.