Optimal Bounded Control for Stationary Response of Strongly Nonlinear Oscillators under Combined Harmonic and Wide-Band Noise Excitations

Abstract

We study the stochastic optimal bounded control for minimizing the stationary response of strongly nonlinear oscillators under combined harmonic and wide-band noise excitations. The stochastic averaging method and the dynamical programming principle are combined to obtain the fully averaged Itô stochastic differential equations which describe the original controlled strongly nonlinear system approximately. The stationary joint probability density of the amplitude and phase difference of the optimally controlled systems is obtained from solving the corresponding reduced Fokker-Planck-Kolmogorov (FPK) equation. An example is given to illustrate the proposed procedure, and the theoretical results are verified by Monte Carlo simulation.

1. Introduction

The well-known tool to solve the problem of stochastic optimal control is the dynamical programming principle, which was proposed by Bellman [1]. According to this principle, a stochastic optimal control problem may be transformed into the problem of finding a solution to the so-called Hamilton-Jacobi-Bellman (HJB) equation. However, in most cases, the HJB equation cannot be solved analytically [2], especially in the multidimensional case. A powerful technique to solve stochastic optimal control problems is to combine the so-called stochastic averaging method [3] with Bellman’s dynamic programming principle. The idea is to replace the original stochastic system by the averaged one and then to apply the dynamic programming principle to the averaged system. It has been justified that the optimal control for the averaged system is nearly optimal for the original system [4]. This combination strategy has two notable advantages. Firstly, the dimension of the averaged system can be reduced remarkably, and the corresponding HJB equation is low-dimensional. Secondly, the diffusion matrix of the original stochastic system is usually singular, while the averaged one is usually nonsingular. This unique characteristic enables the HJB equation for the averaged system to have classical rather than viscous solution. In recent years, a notable nonlinear stochastic optimal control strategy has been proposed for the control of quasi-Hamiltonian systems under random excitations by Zhu based on the stochastic averaging method for quasi-Hamiltonian systems and the stochastic dynamical programming principle [3]. Examples have shown that this strategy is very effective and efficient.

In practice, the excitations of dynamical systems can be classified as either deterministic or random. The random excitation is often modeled as Gaussian white noise, wide-band colored noise, or bounded noise. On the other hand, the magnitudes of the control forces are usually limited due to the saturation in actuators; that is, the control forces are bounded. The optimal bounded control of linear or nonlinear systems under random excitations has been studied by many researchers [5–8]. In all these studies, the excitations of the systems are random excitations alone. However, many physical or mechanical systems are subjected to both random and deterministic harmonic excitations. For example, such combined excitations arise in the study of stochastic resonance [9] or uncoupled flapping motion of rotor blades of a helicopter in forward flight under the effect of atmospheric turbulence [10]. Due to the existence of the deterministic harmonic excitation, the response of the stochastic system is not time homogeneous and the stationary behavior is not easy to capture. Stochastic averaging method is such a method that can approximate the original system by time-homogeneous stochastic processes. By using stochastic averaging method, many researchers have studied the linear or strongly nonlinear systems under combined harmonic and wide-band random excitations [11–16]. On the other hand, optimal bounded control of a linear or nonlinear oscillator subject to combined harmonic and Gaussian white noise excitations has been studied [17–19]. However, Gaussian white noise is an ideal model. In most cases, the random excitation should be modeled as wide-band colored noise. So far, little work has been done on the optimal bounded control of a strongly nonlinear oscillator subject to combined harmonic and wide-band noise excitations [20].

In the present paper, a procedure for designing optimal bounded control to minimize the stationary response of a strongly nonlinear oscillator under combined harmonic and wide-band colored noise excitations is proposed. In Section 2, based on the stochastic averaging method, the equation of motion of a weakly controlled strongly nonlinear oscillator under combined harmonic and wide-band noise excitations is reduced to partially averaged Itô stochastic differential equations. In Section 3, a dynamical programming equation for the control problem of minimizing response of the system is formulated from the partially averaged Itô stochastic differential equations by applying the dynamical programming principle. The optimal control law is determined by the dynamical programming equation and the control constraint. In Section 4, the reduced FPK equation governing the stationary joint probability density of the amplitude and phase difference of the optimally controlled system is established. In Section 5, a nonlinearly damped Duffing oscillator under combined harmonic and wide-band noise excitations is taken as an example to illustrate the application of the proposed procedure. All theoretical results are verified by Monte Carlo simulation.

2. Partially Averaged Itô Stochastic Differential Equations

Consider a weakly controlled strongly nonlinear oscillator subject to weak linear or nonlinear damping and weak external or parametric excitation by a combination of harmonic functions and wide-band colored noises. The equation of motion of the system has the following form𝑋̈+𝑔(𝑋)=𝜀𝑓𝑋,𝑋̇,Ω𝑡+𝜀𝑢𝑋,𝑋̇+𝜀^1/2ℎ_𝑘𝑋,𝑋̇𝜉_𝑘(𝑡),𝑘=1,2,…,𝑟, (2.1) where the term 𝑔(𝑋) is stiffness, which is assumed as an arbitrary nonlinear function of 𝑋. 𝜀 is a small parameter; 𝜀𝑓(𝑋,𝑋̇,Ω𝑡) accounts for light damping and weak harmonic excitation with frequency Ω; 𝜀^1/2ℎ_𝑘𝜉_𝑘(𝑡) represent weak external and/or parametric random excitations; 𝜉_𝑘(𝑡) are wide-band stationary and ergodic noises with zero mean and correlation functions 𝑅_𝑘𝑙(𝜏) or spectral densities 𝑆_𝑘𝑙(𝜔); 𝜀𝑢 denotes a weakly feedback control force. The repeated subscripts indicate summation.

When 𝑢=0, (2.1) describes a large class of physical or structural systems. Under the conditions specified by Xu and Cheung [21], the solutions of system (2.1) can be expressed as the following form:𝑋(𝑡)=𝐴cosΦ(𝑡)+𝐵(𝐴),𝑋̇(𝑡)=−𝐴𝜈(𝐴,Φ)sinΦ(𝑡), (2.2) whereΦ(t)=𝜇(t)+Θ(t), (2.3)𝜈(𝐴,Φ)=𝑑𝜇=𝑑𝑡2[𝑃(𝐴+𝐵)−𝑃(𝐴cosΦ+𝐵)]𝐴²sin²Φ, (2.4)𝑃(𝑥)=^𝑥₀𝑔𝑦𝑑𝑦. (2.5)𝐴, Φ, 𝜇, Θ, and 𝜈   are all stochastic processes. 𝑃(𝑥) is the potential energy of the system (2.1). The functions cosΦ(𝑡) and sinΦ(𝑡) are called generalized harmonic functions. Obviously𝜈(𝐴,Φ) is the instantaneous frequency of the system (2.1).

When 𝜀=0, (2.1) degenerates to the following nonlinear conservative oscillator: 𝑋̈+𝑔(𝑋)=0. (2.6) The average frequency of the conservative oscillator (2.6) can be obtained by the following formula:𝜔(𝐴)=2𝜋∫₀^2𝜋𝑑Φ/𝜈(𝐴,Φ). (2.7) Then the following approximate relation exists:Φ(𝑡)≈𝜔(𝐴)𝑡+Θ(𝑡). (2.8)

Treating (2.2) as generalized Van der Pol transformation from (𝑋,𝑋̇) to (𝐴,Θ), one can obtain the following equations for 𝐴 and Θ:𝑑𝐴𝑑𝑡=𝜀𝐹₁⁽¹⁾(𝐴,Φ,Ω𝑡)+𝜀𝐹₁⁽²⁾(𝐴,Φ,𝑢)+𝜀^1/2𝑈_1𝑘(𝐴,Φ)𝜉_𝑘(𝑡),𝑑Θ𝑑𝑡=𝜀𝐹₂⁽¹⁾(𝐴,Φ,Ω𝑡)+𝜀𝐹₂⁽²⁾(𝐴,Φ,𝑢)+𝜀^1/2𝑈_2𝑘(𝐴,Φ)𝜉_𝑘(𝑡), (2.9) where𝐹₁⁽¹⁾=−𝐴𝑔()𝑓(())𝜈()𝐹𝐴+𝐵)(1+ℎ𝐴cosΦ+𝐵,−𝐴𝜈𝐴,ΦsinΦ,Ω𝑡𝐴,ΦsinΦ,₁₍₂₎=−𝐴𝑢𝑔()𝜈()𝐹𝐴+𝐵)(1+ℎ𝐴,ΦsinΦ,₂₍₁₎=−1𝑔()𝑓(())𝜈(),𝐹𝐴+𝐵)(1+ℎ𝐴cosΦ+𝐵,−𝐴𝜈𝐴,ΦsinΦ,Ω𝑡𝐴,Φ)(cosΦ+ℎ₂₍₂₎=−𝑢𝑔()𝜈(),𝑈𝐴+𝐵)(1+ℎ𝐴,Φ)(cosΦ+ℎ_1𝑘=−𝐴𝑔()ℎ𝐴+𝐵)(1+ℎ^𝑘(())𝜈()𝑈𝐴cosΦ+𝐵,−𝐴𝜈𝐴,ΦsinΦ𝐴,ΦsinΦ,_2𝑘=−1𝑔()ℎ𝐴+𝐵)(1+ℎ^𝑘(𝐴cosΦ+𝐵,−𝐴𝜈(𝐴,Φ)sinΦ)𝜈(𝐴,Φ)(cosΦ+ℎ),ℎ=𝑑𝐵=𝑔()()𝑑𝐴−𝐴+𝐵+𝑔𝐴+𝐵𝑔()().−𝐴+𝐵−𝑔𝐴+𝐵 (2.10)

According to the Stratonovich-Khasminskii limit theorem [22, 23], 𝐴 and Θ converge weakly to 2-dimensional diffusive Markov processes in a time interval of 𝜀⁻¹ order as 𝜀→0, which can be represented by the following Itô stochastic differential equations:𝑑𝐴=𝜀𝑆₁(𝐴,Φ,Ω𝑡)+𝐹₁⁽²⁾(𝐴,Φ,𝑢)𝑑𝑡+𝜀^1/2𝐺_1𝑘(𝐴,Φ)𝑑𝐵_𝑘(𝑡),𝑑Θ=𝜀𝑆₂(𝐴,Φ,Ω𝑡)+𝐹₂⁽²⁾(𝐴,Φ,𝑢)𝑑𝑡+𝜀^1/2𝐺_2𝑘(𝐴,Φ)𝑑𝐵_𝑘(𝑡), (2.11) where 𝐵_𝑘(𝑡) are independent unit Wiener processes,𝑆_𝑖=𝐹_𝑖⁽¹⁾+⁰_−∞𝜕𝑈_𝑖𝑘𝜕𝐴_𝑡⋅𝑈_{1𝑙𝑡+𝜏}𝑅_𝑘𝑙(𝜏)+𝜕𝑈_𝑖𝑘𝜕Φ_𝑡⋅𝑈_{2𝑙𝑡+𝜏}𝑅_𝑘𝑙(𝜏𝑏)𝑑𝜏,_𝑖𝑗=𝜀𝐺_𝑖𝑘𝐺_𝑗𝑘=^∞_−∞𝑈_𝑖𝑘𝑡⋅𝑈_{𝑗𝑙𝑡+𝜏}𝑅_𝑘𝑙(𝜏)𝑑𝜏,𝑖,𝑗=1,2,𝑘,𝑙=1,…,𝑟. (2.12)

System (2.1) has harmonic excitation and two cases can be classified: resonant case and nonresonant case. In the nonresonant case, the harmonic excitation has no effect on the first approximation of the response. Thus, we are interested in the resonant case, namely,Ω𝜔(𝐴)=𝑚𝑛+𝜀𝜎, (2.13) where 𝑚 and 𝑛 are relatively prime positive small integers and 𝜀𝜎 is a small detuning parameter. In this case, multiplying (2.13) by 𝑡 and utilizing the approximate relation (2.8) yield𝑚Ω𝑡=𝑛𝑚Φ+𝜀𝜎𝜇−𝑛Θ. (2.14)

Introduce a new angle variable Γ such that 𝑚Γ=𝜀𝜎𝜇−𝑛Θ (2.15) which is a measure of the phase difference between the response and the harmonic excitation. Then (2.14) can be rewritten as 𝑚Ω𝑡=𝑛Φ+Γ. (2.16)

Using the Itô differential formula, one can obtain the following Itô stochastic differential equations for 𝐴, Γ, and Φ:𝑑𝐴=𝜀𝑆₁(𝐴,Φ,Γ)+𝐹₁⁽²⁾(𝐴,Φ,𝑢)𝑑𝑡+𝜀^1/2𝐺_1𝑘(𝐴,Φ)𝑑𝐵_𝑘(𝑡),(𝐴)−𝑚𝑑Γ=𝜀𝜎𝜔𝑛𝑆₂(𝐴,Φ,Γ)+𝐹₂⁽²⁾(𝐴,Φ,𝑢)𝑑𝑡−𝜀^1/2𝑚𝑛𝐺_2𝑘(𝐴,Φ)𝑑𝐵_𝑘(𝑡),𝑑Φ=𝜔(𝐴)+𝜀𝑆₂(𝐴,Φ,Γ)+𝐹₂⁽²⁾(𝐴,Φ,𝑢)𝑑𝑡+𝜀^1/2𝐺_2𝑘(𝐴,Φ)𝑑𝐵_𝑘(𝑡). (2.17)

Obviously, 𝐴 and Γ are slowly varying processes, while Φ is a rapidly varying process. Averaging the drift and diffusion coefficients in (2.17) with respect to Φ yields the following partially averaged Itô stochastic differential equations: 𝑑𝐴=𝜀𝑚₁⁽¹⁾(𝐴,Γ)+𝐹₁⁽²⁾(𝐴,Φ,𝑢)_Φ𝑑𝑡+𝜀^1/2𝜎_1𝑟(𝐴)𝑑𝐵_𝑟(𝑡),𝑑Γ=𝜀𝑚₂⁽¹⁾()−𝑚𝐴,Γ𝑛𝐹₂⁽²⁾(𝐴,Φ,𝑢)_Φ𝑑𝑡+𝜀^1/2𝜎_2𝑟(𝐴)𝑑𝐵_𝑟(𝑡), (2.18) where𝑚₁⁽¹⁾=𝑆₁(𝐴,Φ,Γ)_Φ,𝑚₂⁽¹⁾(𝐴)−𝑚=𝜎𝜔𝑛𝑆₂(𝐴,Φ,Γ)_Φ,𝑏₁₁(𝐴)=𝜀𝜎_1𝑟𝜎_1𝑟=𝑏₁₁(𝐴,Φ)_Φ,𝑏₂₂(𝐴)=𝜀𝜎_2𝑟𝜎_2𝑟=𝑚²𝑛²𝑏₂₂(𝐴,Φ)_Φ,𝑏₁₂(𝐴)=𝑏₂₁(𝐴)=𝜀𝜎_1𝑟𝜎_2𝑟=𝑚𝑛𝑏₁₂(𝐴,Φ)_Φ,⟨•⟩_Φ=12𝜋₀^2𝜋(•)𝑑Φ. (2.19)

Herein, ⟨•⟩_Φ denotes the averaging with respect to Φ from 0 to 2𝜋⋅𝑏_𝑖𝑗 are diffusion coefficients.

Note that there are two procedures of averaging in this section. One is stochastic averaging, and the other is deterministic time averaging. The procedures for obtaining 𝑆_𝑖 and 𝑏_𝑖𝑗 in (2.12) are called stochastic averaging. The procedures for obtaining 𝑚_𝑖⁽¹⁾ and 𝑏_𝑖𝑗 are called deterministic time averaging. To complete averaging and obtain the explicit expressions of 𝑆_𝑖 and 𝑏_𝑖𝑗, the functions 𝐹_𝑖⁽¹⁾ and 𝑈_𝑖𝑘 in (2.12) are expanded into Fourier series with respect to Φ and the approximate relation (2.8) is used.

3. Dynamical Programming Equation and Optimal Control Law

For mechanical or structural systems, 𝐴(𝑡) is the response amplitude and Γ(𝑡) is the phase difference between the response and the harmonic excitation. They represent the averaged state of system (2.1). Usually, only the response amplitude is concerned, so it is meaningful to control 𝐴(𝑡) in a semi-infinite or finite time interval. Consider the ergodic control problem of system (2.18) in a semi-infinite time interval with the following performance index: 𝐽=lim_𝑇→∞1𝑇^𝑇₀𝑓(𝐴(𝑡))𝑑𝑡. (3.1) Herein, 𝑓 is a cost function. Based on the stochastic dynamical programming principle [24], the following simplified dynamical programming equation can be established from the first equation of (2.18):𝜂=min_𝑢𝑓(𝐴)+𝑚₁⁽¹⁾+𝐹₁⁽²⁾_Φ𝑑𝑉+1𝑑𝐴2𝑏₁₁𝑑²𝑉𝑑𝐴², (3.2) where 𝑉(𝐴) is called value function; min_𝑢∈𝑈[•]denotes the minimum value of [•] with respect to 𝑢; 𝜂 is optimal performance index; 𝑢∈𝑈 indicates the control constraint.

The optimal control law can be determined from minimizing the right-hand side of (3.2) with respect to 𝑢 under control constraint. Suppose that the control constraint is of the form |𝑢(𝐴,Φ)|≤𝑢₀, (3.3) where 𝑢₀ is a positive constant representing the maximum control force. The control 𝑢 enters the performance index (3.2) only through the term𝐹₁⁽²⁾_Φ𝑑𝑉=−𝑢𝑑𝐴𝑔(𝐴+𝐵)(1+ℎ)𝐴𝜈(𝐴,Φ)sinΦ_Φ𝑑𝑉𝑑𝐴. (3.4) This expression attains its minimum value for 𝑢=±𝑢₀, depending on the sign of the coefficient of 𝑢 in (3.4). So the optimal control is𝑢^∗=𝑢₀()sgn𝐴𝜈𝐴,ΦsinΦ𝑔()⋅𝐴+𝐵)(1+ℎ𝑑𝑉𝑑𝐴, (3.5) where sgn denotes sign function.

It is reasonable to assume that the cost function 𝑓 is a monotonously increasing function of 𝐴(𝑡) because the controller must do more work to suppress larger amplitudes 𝐴(𝑡). Then 𝑑𝑉(𝐴)/𝑑𝐴 is positive. Under the specified conditions [21], 𝑔(𝐴+𝐵) and (1+ℎ) are both positive. Equation (3.5) is reduced to 𝑢^∗=−𝑢₀sgn(−𝐴𝜈(𝐴,Φ)sinΦ)=−𝑢₀sgn𝑋̇(𝑡). (3.6)

Equation (3.6) implies that the optimal control is a bang-bang control. 𝑢^∗ has a constant magnitude 𝑢₀. It is in the opposite direction of 𝑋̇(𝑡) and changes its direction at𝑋̇(𝑡)=0.

For the optimal bounded control of system (2.18) in a finite time interval, the following performance index is taken:𝐽=𝐸^𝑡𝑓₀𝑓(𝐴(𝑡))𝑑𝑡+𝜒𝐴𝑡_𝑓, (3.7) where 𝑡_𝑓 is the final control time and 𝜒 is the final cost function. 𝐸[•] denotes the expectation operator. The following simplified dynamical programming equation can be established from the first equation of (2.18):𝜕𝑉𝜕𝑡=−min_𝑢𝑓(𝐴)+𝑚₁⁽¹⁾+𝐹₁⁽²⁾_Φ𝜕𝑉+1𝜕𝐴2𝑏₁₁𝜕²𝑉𝜕𝐴², (3.8) where 𝑉=𝑉(𝐴,𝑡) is the value function. Equation (3.8) is subject to the following final time condition: 𝑉𝐴,𝑡_𝑓=𝜒𝐴𝑡_𝑓. (3.9)

It is obvious that the same optimal control law as that in (3.6) can be derived if the control constraint is of the form of (3.3).

Note that the optimal control for the averaged system (2.18) is nearly optimal for the original system (2.1). For simplicity, here it is called optimal control for both original and averaged systems.

4. Stationary Response of Optimally Controlled System

Inserting 𝑢^∗from (3.6) into (2.18) to replace 𝑢 and averaging𝐹_𝑖⁽²⁾, the following fully averaged Itô stochastic differential equations for 𝐴 and Γ can be obtained:𝑑𝐴=𝜀𝑚₁(𝐴,Γ)𝑑𝑡+𝜀^1/2𝜎_1𝑘(𝐴)𝑑𝐵_𝑘(𝑡),𝑑Γ=𝜀𝑚₂(𝐴,Γ)𝑑𝑡+𝜀^1/2𝜎_2𝑘(𝐴)𝑑𝐵_𝑘(𝑡), (4.1) where𝑚₁(𝐴,Γ)=𝑚₁⁽¹⁾+−𝑢^∗𝐴𝜈(𝐴,Φ)sinΦ𝑔(𝐴+𝐵)(1+ℎ)_Φ,𝑚₂(𝐴,Γ)=𝑚₂⁽¹⁾+𝑚𝑛𝑢^∗𝜈(𝐴,Φ)(cosΦ+ℎ)𝑔(𝐴+𝐵)(1+ℎ)_Φ. (4.2)

The reduced FPK equation for the optimally controlled system is of the following form𝜕0=−𝜕𝑎𝑚₁𝑝−𝜕𝜕𝛾𝑚₂𝑝+12𝜕²𝜕𝑎²𝑏₁₁𝑝+𝜕²𝜕𝑎𝜕𝛾𝑏₁₂𝑝+12𝜕²𝜕𝛾²𝑏₂₂𝑝, (4.3) where 𝑝=𝑝(𝑎,𝛾) is the stationary joint probability density of the amplitude 𝐴 and the phase difference Γ. Since 𝑝 is a periodic function of 𝛾, it satisfies the following periodic boundary condition with respect to 𝛾: 𝑝𝑎,𝛾=𝑝𝑎,𝛾+2𝜋. (4.4)

The boundary condition with respect to 𝑎=0 is 𝑝=finiteat𝑎=0 (4.5) which implies that 𝑎=0 is a reflecting boundary. The other boundary condition is 𝑝,𝜕𝑝𝜕𝑎⟶0as𝑎⟶∞. (4.6)

In addition to the boundary conditions, the stationary joint probability density 𝑝(𝑎,𝛾) satisfies the following normalization condition:₀^2𝜋∞₀𝑝𝑎,𝛾𝑑𝑎𝑑𝛾=1. (4.7)

Usually, the partial differential equation (4.3) can be solved only numerically.

5. Example

To illustrate the proposed strategy in the previous sections, take the following controlled nonlinearly damped Duffing oscillator as an example. Duffing oscillator is a typical model in nonlinear analysis. The equation of motion of the system is of the form𝑋̈+𝛽₁+𝛽₂𝑋²𝑋̇+𝜔²₀𝑋+𝛼𝑋³=𝐸cosΩ𝑡+𝜉₁(𝑡)+𝑋𝜉₂(𝑡)+𝑢, (5.1) where 𝛽₁, 𝛽₂, 𝜔₀, 𝛼, 𝐸, and Ω are positive constants; 𝑢 is the feedback control with the constraint defined by (3.3); 𝜉_𝑘(𝑡)(𝑘=1,2) are independent stationary and ergodic wide-band noises with zero mean and rational spectral densities𝑆_𝑖(𝜔)=𝐷_𝑖𝜋1𝜔²+𝜔²_𝑖,𝑖=1,2. (5.2)

𝜉_𝑖(𝑡) can be regarded as the output of the following first-order linear filter:𝜉̇_𝑖+𝜔_𝑖𝜉_𝑖=𝑊_𝑖(𝑡),𝑖=1,2, (5.3) where 𝑊_𝑖(𝑡) are Gaussian white noises in the sense of Stratonovich with intensities 2𝐷_𝑖. Note that the maximum value of 𝑆_𝑖(𝜔) is D_𝑖/(𝜋𝜔²_𝑖). It is assumed that 𝛽₁, 𝐷_𝑖/(𝜋𝜔²_𝑖), and 𝐸 are all small.

For the system (5.1), the instantaneous frequency defined by (2.4) has the following form:𝜈()=𝜔𝐴,Φ²⁰+3𝛼𝐴²4(1+𝜆cos2Φ)^1/2,𝜆=𝛼𝐴²4𝜔²₀+3𝛼𝐴²./4 (5.4)𝜈(𝐴,Φ) can be approximated by the following finite sum with a relative error less than 0.03%:𝜈(A,Φ)=𝑏₀(A)+𝑏₂(A)cos2Φ+𝑏₄(A)cos4Φ+𝑏₆(A)cos6Φ, (5.5) where𝑏₀(A)=𝜔²₀+3𝛼A²4^1/2𝜆1−²,𝑏16₂(A)=𝜔²₀+3𝛼A²4^1/2𝜆2+3𝜆³,𝑏64₄(A)=𝜔²₀+3𝛼A²4^1/2−𝜆²,𝑏16₆(A)=𝜔²₀+3𝛼A²4^1/2𝜆³.64 (5.6)

Then the averaged frequency 𝜔(𝐴) of the system (5.1) can be approximated by 𝑏₀(𝐴). In the case of primary external resonance,Ω𝜔(A)=1+𝜀𝜎. (5.7)

Introduce the new angel variable defined by (2.15), and complete the procedures shown in Sections 2–4 we obtain the following fully averaged Itô stochastic differential equations for 𝐴 and Γ:𝑑𝐴=𝑚₁(𝐴,Γ)𝑑𝑡+𝑏₁₁(𝐴)𝑑𝐵₁(𝑡),𝑑Γ=𝑚₂(𝐴,Γ)𝑑𝑡+𝑏₂₂(𝐴)𝑑𝐵₂(𝑡), (5.8) where 𝑚_𝑖 and 𝑏_𝑖𝑖(𝑖=1,2) are drift and diffusion coefficients, respectively. 𝑚_𝑖 and 𝑏_𝑖𝑖 are given in the appendix.

The stationary joint probability density 𝑝(𝑎,𝛾) of the optimally controlled system (5.1) is governed by the following reduced FPK equation:𝜕0=−𝜕𝑎𝑚₁𝑝−𝜕𝜕𝛾𝑚₂𝑝+12𝜕²𝜕𝑎²𝑏₁₁𝑝+12𝜕²𝜕𝛾²𝑏₂₂𝑝. (5.9)

Solving the FPK equation (5.9) by finite difference method under the conditions (4.4)–(4.7), the stationary joint probability density of the amplitude and the phase difference of the optimally controlled system (5.1) can be obtained. Furthermore, the stationary mean amplitude of the optimally controlled system (5.1) can be obtained as follows:𝐸[𝑎]=^∞₀₀^2𝜋𝑎𝑝𝑎,𝛾𝑑𝑎𝑑𝛾. (5.10)

To check the accuracy of the proposed method, Monte Carlo digital simulation of the original system (5.1) is performed. The sample functions of independent wide-band noise 𝜉_𝑖(𝑡) were generated by inputting Gaussian white noises to the linear filter (5.3). The response of system (5.1) was obtained numerically by using the fourth-order Runge-Kutta method with time step 0.02. The long-time solution after 1500,000 steps was regarded as the stationary ergodic response. 100 samples are used. For every sample, the amplitude 𝐴 and angle variable Γ are calculated from step 1500,001 to step 2000000 to obtain the statistical probability density of 𝑝(𝑎,𝛾). Figures 1 and 2 show the typical sample function of wide-band noise 𝜉₁(𝑡) and control force 𝑢, respectively. Figure 3 shows 𝑝(𝑎,𝛾) of the optimal controlled system (5.1). It is seen that the theoretical result agrees very well with that from digital simulation. Figures 4(a) and 4(b) show the sample functions of displacement and velocity of the system (5.1), respectively. It is obvious that the bounded control can reduce the displacement and velocity. Also, the amplitude of the original system (5.1) is reduced by control, which is verified by Figure 5.

Figure 1 

Sample function of wide-band noise 𝜉1(𝑡)⋅𝐷1=10, 𝜔1=30.

Figure 2 

Sample function of the control force 𝑢, 𝑢0=0.1.

(a)

(b)

(c)

(d)

Figure 3 

Stationary joint probability density 𝑝(𝑎,𝛾) of the optimally controlled system (5.1) in primary external resonance case. 𝜔0=1.0, Ω=1.1,  𝛼=1.0, 𝐸=0.3, 𝛽𝑖=0.05, 𝜔𝑖=30, 𝐷1=10, 𝐷2=5, 𝑢0=0.1 (𝑖=1,2). (a) Theoretical result; (b) Monte Carlo simulation of the original system (5.1); (c) stationary marginal probability density 𝑝(𝑎); (d) stationary marginal probability density 𝑝(𝛾). — theoretical results; ● results from Monte Carlo simulation.

(a)

(b)

Figure 4 

Sample functions of displacement and velocity of the system (5.1), 𝑢0=0.4. Other parameters are the same as those in Figure 3. Red line: with control; blue line: without control.

Figure 5 

Mean amplitude of the response of the optimally controlled system (5.1) in primary external resonance. The parameters are the same as those in Figure 3 except that 𝑢0 is a variable. — theoretical results; ●  ▲ results from Monte Carlo simulation.

Note that the system (5.1) is strongly nonlinear. Increasing the nonlinearity coefficients 𝛽₁ or 𝛼 in (5.1), one can see that the agreement between theoretical results and Monte Carlo digital simulation is still acceptable (see Figures 6 and 7). This demonstrates that the proposed method is powerful to deal with strongly nonlinear problems, even though the nonlinearity is extremely strong.

Figure 6 

Stationary marginal probability density 𝑝(𝑎) of the optimally controlled system (5.1) in primary external resonance case. Other parameters are the same as those in Figure 3. — theoretical results; ●  ▲  ◆ results from Monte Carlo simulation.

Figure 7 

Stationary marginal probability density 𝑝(𝑎) of the optimally controlled system (5.1) in primary external resonance case. Other parameters are the same as those in Figure 3. — theoretical results; ●  ▲  ◆ results from Monte Carlo simulation.

6. Conclusions

In the present paper, a combination procedure of the stochastic averaging method and Bellman’s dynamic programming for designing the optimal bounded control to minimize the response of strongly nonlinear systems under combined harmonic and wide-band noise excitations has been proposed. The procedure consists of applying the stochastic averaging method for weakly controlled strongly nonlinear systems under combined harmonic and wide-band noise excitations, establishing the dynamical programming equation for the control problem of minimizing the response based on the partially averaged Itô stochastic differential equations and the dynamical programming principle, determining the optimal control from the dynamical programming equation and the control constraint without solving the dynamical programming equation. Then the stationary joint probability density and mean amplitude of the optimally controlled averaged system are obtained from solving the reduced FPK equation associated with the fully averaged Itô stochastic differential equations. A nonlinearly damped Duffing oscillator with hardening stiffness has been taken as an example to illustrate the application of the proposed procedure. The comparison between the theoretical results and those from Monte Carlo simulation shows that the proposed procedure works quite well even though the nonlinearity is extremely strong. The results show that the response amplitude of the system can be reduced remarkably by the feedback control.

The advantages of the proposed method are obvious. Note that; in the example, the wide-band noise is generated by a first-order linear filter. In principle, one could apply the stochastic dynamical programming method to the extended system [𝑋,𝑋̇,𝜉₁,𝜉₂]^𝑇 to study the optimal control problem. However, the corresponding HJB equation is 5-dimensional or 4 dimensional. The corresponding reduced FPK equation is 4 dimensional, which is very difficult to solve. After stochastic averaging, the original system is represented by two-dimensional time-homogeneous diffusion Markov processes of amplitude and phase difference with nondegenerate diffusion matrix. The dynamical programming equation derived from the averaged equations is two dimensional or one dimensional. The corresponding reduced FPK equation is two dimensional, which is easy to solve. The other advantage of the proposed procedure is that it is not necessary to solve the dynamical programming equation for obtaining the optimal control law. Furthermore, the proposed method can be extended to multi-degrees-of-freedom (MDOF) systems easily. This will be our future work.

Appendix

The drift coefficients in (5.8) are as follows: 𝑚₁=𝐻₁(𝐴)+𝐹₁₀(𝐴,Γ)+2𝑢₀−105𝑏₀(𝐴)+35𝑏₂(𝐴)+7𝑏₄(𝐴)+3𝑏₆(𝐴)105𝜋𝜔²₀+𝛼𝐴²,𝑚₂=𝐻₂(𝐴)+Ω−𝑏₀(𝐴)+𝐸cosΓ2𝑏₀(𝐴)+𝑏₂(𝐴)4𝐴𝛼𝐴²+𝜔²₀,𝐹₁₀(𝐴,Γ)=𝐸sinΓ2𝑏₀(𝐴)−𝑏₂(𝐴)4𝛼𝐴²+𝜔²₀−𝐴𝛽₁16𝜔²₀+10𝛼𝐴²+𝐴²𝛽₂4𝜔²₀+3𝛼𝐴²32𝛼𝐴²+𝜔²₀,𝐻₁(𝐴)=𝑚₁₁+𝑚₁₂+𝑚₁₃+𝑚₁₄,𝐻₂(𝐴)=𝑚₂₁+𝑚₂₂+𝑚₂₃+𝑚₂₄,𝑚₁₁=𝑚₁₁₁𝑆₁(𝜔(𝐴))+𝑚₁₁₃𝑆₁(3𝜔(𝐴))+𝑚₁₁₅𝑆₁(5𝜔(𝐴))+𝑚₁₁₇𝑆₁(7𝜔(𝐴)),𝑚₁₃=𝑚₁₃₁𝑆₁(𝜔(𝐴))+𝑚₁₃₃𝑆₁(3𝜔(𝐴))+𝑚₁₃₅𝑆₁(5𝜔(𝐴))+𝑚₁₃₇𝑆₁(7𝜔(𝐴)),𝑚₁₂=𝑚₁₂₂𝑆₂(2𝜔(𝐴))+𝑚₁₂₄𝑆₂(4𝜔(𝐴))+𝑚₁₂₆𝑆₂(6𝜔(𝐴))+𝑚₁₂₈𝑆₂(8𝜔(𝐴)),𝑚₁₄=𝑚₁₄₂𝑆₂(2𝜔(𝐴))+𝑚₁₄₄𝑆₂(4𝜔(𝐴))+𝑚₁₄₆𝑆₂(6𝜔(𝐴))+𝑚₁₄₈𝑆₂(8𝜔(𝐴)),𝑚₂₁=𝑚₂₁₁𝐼₁(𝜔(𝐴))+𝑚₂₁₃𝐼₁(3𝜔(𝐴))+𝑚₂₁₅𝐼₁(5𝜔(𝐴))+𝑚₂₁₇𝐼₁(7𝜔(𝐴)),𝑚₂₃=𝑚₂₃₁𝐼₁(𝜔(𝐴))+𝑚₂₃₃𝐼₁(3𝜔(𝐴))+𝑚₂₃₅𝐼₁(5𝜔(𝐴))+𝑚₂₃₇𝐼₁(7𝜔(𝐴)),𝑚₂₂=𝑚₂₂₂𝐼₂(2𝜔(𝐴))+𝑚₂₂₄𝐼₂(4𝜔(𝐴))+𝑚₂₂₆𝐼₂(6𝜔(𝐴))+𝑚₂₂₈𝐼₂(8𝜔(𝐴)),𝑚₂₄=𝑚₂₄₂𝐼₂(2𝜔(𝐴))+𝑚₂₄₄𝐼₂(4𝜔(𝐴))+𝑚₂₄₆𝐼₂(6𝜔(𝐴))+𝑚₂₄₈𝐼₂((𝐴,𝐼8𝜔))_𝑖(𝜔)=𝐷_𝑖𝜔_𝑖𝜔𝜔²+𝜔²_𝑖(𝑖=1,2),𝑚₁₁₁=𝜋𝑏₂(𝐴)−2𝑏₀(𝐴×)2𝛼𝐴2𝑏₀(𝐴)−𝑏₂(𝐴)+𝐴²𝛼+𝜔²₀𝑑𝑏₂(𝐴)/𝑑𝐴−2𝑑𝑏₀(𝐴)/𝑑𝐴8𝐴²𝛼+𝜔²₀³,𝑚₁₁₃=𝜋𝑏₂(𝐴)−𝑏₄(𝐴×𝑏)2𝛼𝐴₄(𝐴)−𝑏₂(𝐴)+𝐴²𝛼+𝜔²₀𝑑𝑏₂(𝐴)/𝑑𝐴−𝑑𝑏₄(𝐴)/𝑑𝐴8𝐴²𝛼+𝜔²₀³,𝑚₁₁₅=𝜋𝑏₄(𝐴)−𝑏₆(𝐴×𝑏)2𝛼𝐴₆(𝐴)−𝑏₄(𝐴)+𝐴²𝛼+𝜔²₀𝑑𝑏₄(𝐴)/𝑑𝐴−𝑑𝑏₆(𝐴)/𝑑𝐴8𝐴²𝛼+𝜔²₀³,𝑚₁₁₇=𝜋𝑏₆(𝐴)−2𝛼𝐴𝑏₆(𝐴)+𝐴²𝛼+𝜔²₀𝑑𝑏₆(𝐴)/𝑑𝐴8𝐴²𝛼+𝜔²₀³,𝑚₁₂₂=𝜋𝐴2𝑏₀(𝐴)−𝑏₄(𝐴×)2𝑏₀(𝐴)−𝑏₄(𝐴)𝐴²𝛼−𝜔²₀−𝐴𝐴²𝛼+𝜔²₀2𝑑𝑏₀(𝐴)/𝑑𝐴−𝑑𝑏₄(𝐴)/𝑑𝐴32𝐴²𝛼+𝜔²₀³,𝑚₁₂₄=𝜋𝐴𝑏₂(𝐴)−𝑏₆(𝐴×𝑏₆)(𝐴)−𝑏₂(𝐴)𝐴²𝛼−𝜔²₀+𝐴𝐴²𝛼+𝜔²₀𝑑𝑏₂(𝐴)/𝑑𝐴−𝑑𝑏₆(𝐴)/𝑑𝐴32𝐴²𝛼+𝜔²₀³,𝑚₁₂₆=𝜋𝐴𝑏₄(𝐴)𝑏₄(𝐴)𝐴²𝛼−𝜔²₀+𝐴𝐴²𝛼+𝜔²₀𝑑𝑏₄(𝐴)/𝑑𝐴32𝐴²𝛼+𝜔²₀³,𝑚₁₂₈=𝜋𝐴𝑏₆(𝐴)𝑏₆(𝐴)𝐴²𝛼−𝜔²₀+𝐴𝐴²𝛼+𝜔²₀𝑑𝑏₆(𝐴)/𝑑𝐴32𝐴²𝛼+𝜔²₀³,𝑚₁₃₁=−𝜋𝑏²₂(𝐴)−4𝑏²₀(𝐴)8𝐴𝐴²𝛼+𝜔²₀²,𝑚₁₃₃=−3𝜋𝑏²₄(𝐴)−𝑏²₂(𝐴)8𝐴𝐴²𝛼+𝜔²₀²,𝑚₁₃₅=−5𝜋𝑏²₆(𝐴)−𝑏²₄(𝐴)8𝐴𝐴²𝛼+𝜔²₀²,𝑚₁₃₇=7𝜋𝑏²₆(𝐴)8𝐴𝐴²𝛼+𝜔²₀²,𝑚₁₄₂=𝜋𝐴2𝑏₀(𝐴)−𝑏₄(𝐴)2𝑏₀(𝐴)+2𝑏₂(𝐴)+𝑏₄(𝐴)16𝐴²𝛼+𝜔²₀²,𝑚₁₄₄=𝜋𝐴𝑏₂(𝐴)−𝑏₆(𝐴)𝑏₂(𝐴)+2𝑏₄(𝐴)+𝑏₆(𝐴)8𝐴²𝛼+𝜔²₀²,𝑚₁₄₆=3𝜋𝐴𝑏₄(𝐴)𝑏₄(𝐴)+2𝑏₆(𝐴)16𝐴²𝛼+𝜔²₀²,𝑚₁₄₈=𝜋𝐴𝑏²₆(𝐴)4𝐴²𝛼+𝜔²₀²,𝑚₂₁₁=2𝑏₀(𝐴)+𝑏₂(𝐴)²8𝐴²𝛼𝐴²+𝜔²₀²,𝑚₂₁₃=3𝑏₂(𝐴)+𝑏₄(𝐴)²8𝐴²𝛼𝐴²+𝜔²₀²,𝑚₂₁₅=5𝑏₄(𝐴)+𝑏₆(𝐴)²8𝐴²𝛼𝐴²+𝜔²₀²,𝑚₂₁₇=7𝑏₆(𝐴)²8𝐴²𝛼𝐴²+𝜔²₀²,𝑚₂₂₂=2𝑏₀(𝐴)+2𝑏₂(𝐴)+𝑏₄(𝐴)²16𝛼𝐴²+𝜔²₀²,𝑚₂₂₄=𝑏₂(𝐴)+2𝑏₄(𝐴)+𝑏₆(𝐴)²8𝛼𝐴²+𝜔²₀²,𝑚₂₂₆=3𝑏₄(𝐴)+2𝑏₆(𝐴)²16𝛼𝐴²+𝜔²₀²,𝑚₂₂₈=𝑏₆(𝐴)²4𝛼𝐴²+𝜔²₀²,𝑚₂₃₁=2𝑏₀(𝐴)−𝑏₂(𝐴×−𝑏₂)(𝐴)+2𝑏₀(𝐴)3𝐴²𝛼+𝜔²₀+𝐴𝛼𝐴²+𝜔²₀2𝑑𝑏₀(𝐴)/𝑑𝐴+𝑑𝑏₂(𝐴)/𝑑𝐴8𝐴²𝛼𝐴²+𝜔²₀³,𝑚₂₃₃=𝑏₂(𝐴)−𝑏₄(𝐴×−𝑏₂)(𝐴)+𝑏₄(𝐴)3𝐴²𝛼+𝜔²₀+𝐴𝛼𝐴²+𝜔²₀𝑑𝑏₂(𝐴)/𝑑𝐴+𝑑𝑏₄(𝐴)/𝑑𝐴8𝐴²𝛼𝐴²+𝜔²₀³,𝑚₂₃₅=𝑏₄(𝐴)−𝑏₆(𝐴)8𝐴²𝛼𝐴²+𝜔²₀³×−𝑏₄(𝐴)+𝑏₆(𝐴)3𝐴²𝛼+𝜔²₀+𝐴𝛼𝐴²+𝜔²₀𝑑𝑏₄(𝐴)/𝑑𝐴+𝑑𝑏₆(𝐴)/𝑑𝐴8𝐴²𝛼𝐴²+𝜔²₀³,𝑚₂₃₇=𝑏₆(𝐴)−𝑏₆(𝐴)3𝐴²𝛼+𝜔²₀+𝐴𝛼𝐴²+𝜔²₀𝑑𝑏₆(𝐴)/𝑑𝐴8𝐴²𝛼𝐴²+𝜔²₀³,𝑚₂₄₂=−2𝐴𝛼2𝑏₀(𝐴)+2𝑏₂(𝐴)+𝑏₄(𝐴)+𝐴²𝛼+𝜔²₀2𝑑𝑏₀(𝐴)+𝑑𝐴2𝑑𝑏₂(𝐴)+𝑑𝐴𝑑𝑏₄(𝐴)×𝐴𝑑𝐴2𝑏₀(𝐴)−𝑏₄(𝐴)32𝛼𝐴²+𝜔²₀³,𝑚₂₄₄=−2𝐴𝛼𝑏₂(𝐴)+2𝑏₄(𝐴)+𝑏₆(𝐴)+𝐴²𝛼+𝜔²₀𝑑𝑏₂(𝐴)+𝑑𝐴2𝑑𝑏₄(𝐴)+𝑑𝐴𝑑𝑏₆(𝐴)×𝐴𝑏𝑑𝐴₂(𝐴)−𝑏₆(𝐴)32𝛼𝐴²+𝜔²₀³,𝑚₂₄₆=𝐴𝑏₄(𝐴)×−2𝐴𝛼𝑏₄(𝐴)+2𝑏₆(𝐴)+𝐴²𝛼+𝜔²₀𝑑𝑏₄(𝐴)/𝑑𝐴+2𝑑𝑏₆(𝐴)/𝑑𝐴32𝛼𝐴²+𝜔²₀³,𝑚₂₄₈=𝐴𝑏₆(𝐴)−2𝐴𝛼𝑏₆(𝐴)+𝐴²𝛼+𝜔²₀𝑑𝑏₆(𝐴)/𝑑𝐴32𝛼𝐴²+𝜔²₀³. (A.1)

The diffusion coefficients in (5.8) are as follows: 𝑏₁₁=𝑏₁₁₁+𝑏₁₁₂,𝑏₂₂=𝑏₂₂₁+𝑏₂₂₂,𝑏_𝑖𝑗=0,𝑖≠𝑗,𝑏₁₁₁=𝑏₁₁₁₁𝑆₁(𝜔(𝐴))+𝑏₁₁₁₃𝑆₁(3𝜔(𝐴))+𝑏₁₁₁₅𝑆₁(5𝜔(𝐴))+𝑏₁₁₁₇𝑆₁(7𝜔(𝐴)),𝑏₂₂₁=𝑏₂₂₁₁𝑆₁(𝜔(𝐴))+𝑏₂₂₁₃𝑆₁(3𝜔(𝐴))+𝑏₂₂₁₅𝑆₁(5𝜔(𝐴))+𝑏₂₂₁₇𝑆₁(7𝜔(𝐴)),𝑏₁₁₂=𝑏₁₁₂₂𝑆₂(2𝜔(𝐴))+𝑏₁₁₂₄𝑆₂(4𝜔(𝐴))+𝑏₁₁₂₆𝑆₂(6𝜔(𝐴))+𝑏₁₁₂₈𝑆₂(8𝜔(𝐴)),𝑏₂₂₂=𝑏₂₂₂₀𝑆₂(0)+𝑏₂₂₂₂𝑆₂(2𝜔(𝐴))+𝑏₂₂₂₄𝑆₂(4𝜔(𝐴))+𝑏₂₂₂₆𝑆₂(6𝜔(𝐴))+𝑏₂₂₂₈𝑆₂(8𝜔(𝐴)),𝑏₁₁₁₁=𝜋𝑏₂(𝐴)−2𝑏₀(𝐴)²4𝛼𝐴²+𝜔²₀²,𝑏₁₁₁₃=𝜋𝑏₂(𝐴)−𝑏₄(𝐴)²4𝛼𝐴²+𝜔²₀²,𝑏₁₁₁₅=𝜋𝑏₄(𝐴)−𝑏₆(𝐴)²4𝛼𝐴²+𝜔²₀²,𝑏₁₁₁₇=𝜋𝑏²₆(𝐴)4𝛼𝐴²+𝜔²₀²,𝑏₁₁₂₂=𝜋𝐴²𝑏₄(𝐴)−2𝑏₀(𝐴)²16𝛼𝐴²+𝜔²₀²,𝑏₁₁₂₄=𝜋𝐴²𝑏₂(𝐴)−𝑏₆(𝐴)²16𝛼𝐴²+𝜔²₀²,𝑏₁₁₂₆=𝜋𝐴²𝑏²₄(𝐴)16𝛼𝐴²+𝜔²₀²,𝑏₁₁₂₈=𝜋𝐴²𝑏²₆(𝐴)16𝛼𝐴²+𝜔²₀²,𝑏₂₂₁₁=𝜋2𝑏₀(𝐴)+𝑏₂(𝐴)²4𝐴²𝛼𝐴²+𝜔²₀²,𝑏₂₂₁₃=𝜋𝑏₂(𝐴)+𝑏₄(𝐴)²4𝐴²𝛼𝐴²+𝜔²₀²,𝑏₂₂₁₅=𝜋𝑏₄(𝐴)+𝑏₆(𝐴)²4𝐴²𝛼𝐴²+𝜔²₀²,𝑏₂₂₁₇=𝜋𝑏²₆(𝐴)4𝐴²𝛼𝐴²+𝜔²₀²,𝑏₂₂₂₀=𝜋2𝑏₀(𝐴)+𝑏₂(𝐴)²8𝛼𝐴²+𝜔²₀²,𝑏₂₂₂₂=𝜋2𝑏₀(𝐴)+2𝑏₂(𝐴)+𝑏₄(𝐴)²16𝛼𝐴²+𝜔²₀²,𝑏₂₂₂₄=𝜋𝑏₂(𝐴)+2𝑏₄(𝐴)+𝑏₆(𝐴)²16𝛼𝐴²+𝜔²₀²,𝑏₂₂₂₆=𝜋𝑏₄(𝐴)+2𝑏₆(𝐴)²16𝛼𝐴²+𝜔²₀²,𝑏₂₂₂₈=𝜋𝑏²₆(𝐴)16𝛼𝐴²+𝜔²₀². (A.2)

Acknowledgments

The work reported in this paper was supported by the National Natural Science Foundation of China under Grant nos. 10802030, 10902096 and Specialized Research Fund for Doctoral Program of Higher Education of China under Grant no. 200802511005.