宽带噪声激励下带有分数阶控制器的强非线性系统的随机平均技术

梁霄，陈林聪，赵珧冰; Liang Xiao; Chen Lincong; Zhao Yaobing

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

摘要

运用随机平均法研究了宽带噪声激励下带有分数阶PI^λD^μ控制器的强非线性系统.首先，应用广义谐波平衡技术，将分数阶PI^λD^μ控制力分解为幅值依赖的等效拟线性阻尼力和拟线性回复力，得到了受控整数阶等效非线性系统.然后，运用基于广义谐和函数的随机平均法得到关于幅值的平均伊藤微分方程.最后，建立并求解相应的简化Fokker-Planck-Kolmogorov（FPK）方程，得到稳态概率密度函数.作为算例，考察了Duffing-van de Pol 振子.数值结果表明随机平均法能够达到较高的精度，分数阶PI^λD^μ控制器能够对系统响应进行有效的控制.此外，宽带噪声参数ξ_i、ω_i及D_i改变时，本文提出的方法仍具有较好的适用性，分数阶控制器仍同样具有非常好的控制效果.

关键词

分数阶控制器; 随机平均法; 宽带噪声; 强非线性系统; 蒙特卡罗模拟

引言

分数阶微积分理论及其应用的研究在过去的几十年中得到了广泛的开展^[

1-6].随着分数阶微积分理论的研究取得长足发展，其应用也逐渐受到重视，特别是与现代控制理论相结合开辟了一个全新的领域——分数阶控制.传统的控制器如：朱位秋等^{[参考文献 7
百度学术}7，8]基于拟哈密顿系统的随机平均法^{[参考文献 9
百度学术}9，10]和随机动态规划原理^{[参考文献 11
百度学术}11]，设计了Bang-Bang控制器，对谐和与随机激励下Duffing系统的随机跳跃分岔进行控制；冯长水等^{[参考文献 12
百度学术}12，13]通过改变时滞参数，研究了基于时滞反馈控制器的窄带噪声激励下Duffing振子的分岔控制；许勇等^{[参考文献 14
百度学术}14]采用正交多项式逼近法，并利用随机反馈控制方法，对一类随机非线性系统的Hopf分岔进行了控制等.虽然，这些控制方法以及控制器已取得了一定的成果，但仍存在一定的问题，例如：Bang-Bang控制的不连续性，可能会导致所得结果存在误差；反馈控制器由于时滞的存在，会影响系统的稳定性，给系统带来更复杂的动力学现象.而分数阶控制器可以在避免这些误差的同时能够获得更好的动力学行为和鲁棒性.

目前，已有的分数阶控制器如：CRONE控制器^[

15]、分数阶PI^λD^μ控制器^{[参考文献 16
百度学术}16]、TID控制器^{[参考文献 17
百度学术}17]、超前滞后校正补偿器^{[参考文献 18
百度学术}18].Chen和Moore^{[参考文献 19
百度学术}19]把分数阶控制的理论应用到实际的系统控制（如分数阶电容）中.应用表明，分数阶控制器比整数阶控制器的控制效果更好.陈林聪等^{[参考文献 20
百度学术}20]利用分数阶PID控制器对有界噪声激励下Duffing系统的随机P分岔进行控制.Badreddine Boudjehem和Djalil Boudjehem^{[参考文献 21
百度学术}21]设计了基于最小性能指标的分数阶PID控制器.Shah和Agashe^{[参考文献 22
百度学术}22]对分数阶PI^λD^μ控制器作了全面的评述.此外，在航空航天、通信等其他领域分数阶控制器也都得到了很广泛的应用^{[参考文献 23-30}23-30].

目前，对分数阶控制器的研究大都集中在白噪声激励情形.如:李伟等^[

28]研究了基于分数阶PID控制器的高斯白噪声激励下随机动力系统的可靠性；最近陈林聪等^{[参考文献 29
百度学术}29]运用随机平均法研究了高斯白噪声激励下带有时滞反馈的分数阶PD控制器的单自由度非线性系统，研究表明时滞会破坏分数阶控制器的控制效果，且时滞参数和分数阶次的改变会引起随机P分岔.然而，白噪声是一个理想的过程，忽略了实际噪声在小时间间隔内相关的事实，并不能代表真实噪声.

随机平均法是用扩散的马尔柯夫过程来近似响应.近年来，随机平均法被扩展到研究具有分数阶控制的非线性随机动力学系统^[

30-32].本文运用随机平均技术研究了宽带噪声激励下带有分数阶PI^λD^μ控制器的强非线性系统.首先,应用广义谐波平衡技术将分数阶PI^λD^μ控制力解耦为幅值相关的拟线性阻尼力和拟线性恢复力,得到受控整数阶系统的等效非线性系统.然后,应用基于广义谐和函数的随机平均法得到平均漂移系数和扩散系数,进而得到关于系统幅值的平均伊藤方程.最后，求解相应的FPK方程得到系统的稳态响应概率密度.作为算例，研究了Duffing-van de Pol 振子，研究表明参数ξ_i、ω_i及D_i发生变化时，本文所提出的方法同样具有较好的适应性，近似解析解与蒙特卡罗模拟结果十分吻合，分数阶控制仍同样具有非常好的控制效果.

1 问题阐述

考察宽带噪声作用下的强非线性振子，其运动方程为：

\ddot{X} (t) + ε f (X, \dot{X}) \dot{X} (t) + g (X)

= ε^{1 / 2} h_{i} (X, \dot{X}) W_{i} (t) + ε u

(1)

其中， $X (t), \dot{X} (t)$ 分别表示广义位移和广义速度，ε是正小参数； $ε f (X, \dot{X})$ 是弱非线性阻尼系数；g(X)是线性或非线性回复力； $ε^{1 / 2} h (X, \dot{X})$ 是随机力的振幅； $ε u = ε k_{0} X (t) + ε k_{1} I^{λ} X (t) + ε k_{2} D^{μ} X (t)$ 是分数阶PI^λD^μ控制力，其中，k₀, k₁, k₂为常数.I^λX(t)和D^μX(t)分别是Riemann-Liouville定义下的分数阶积分和分数阶微分

D^{μ} X (t) = \frac{1}{Γ (1 - μ)} \frac{d}{d t} \int_{0}^{t} \frac{X (t - τ)}{τ^{μ}} d τ

I^{λ} X (t) = \frac{1}{Γ (λ)} \int_{0}^{t} \frac{X (t - τ)}{τ^{1 - λ}} d τ,

0 < λ, μ < 1

(2)

W_i(t)为宽带噪声激励，其有理谱密度

S_{i} (ω) = \frac{D_{i}}{π} \frac{1}{(ω^{2} - {ω_{i}}^{2})^{2} + 4 ξ_{i}^{2} ω^{2} ω_{i}^{2}}

i = 1,2, 3

(3)

ξ_i,ω_i, D_i为常数.

当ε=0时，系统(1)退化为非线性保守振子

\ddot{x} + g (x) = 0

(4)

系统的总能量H为

H = {\dot{x}}^{2} / 2 + U (x)

(5)

式中

U (x) = \int_{0}^{x} g (u) d u

(6)

表示为系统的势能.假设 $U (x) \geq 0$ 相对于x对称,且在x=0处取得最小值，则系统(4)存在一族周期解

x (t) = a c o s θ

\dot{x} (t) = - a ν (a, θ) s i n θ

(7)

式中

θ = ϕ + γ_{0}

(8)

ν (a, θ) = \sqrt[]{\frac{2 [U (a) - U (a c o s θ)]}{a^{2} s i n^{2} θ}}

(9)

这里, $c o s θ$ 和 $s i n θ$ 称为广义谐和函数； $a$ 为位移的幅值； $ν (a, θ)$ 为振子的瞬时频率； $γ_{0}$ 为振子的初相位角.

将 $ν (a, θ)$ 展开为Fourier级数，得

ν (a, θ) = b_{0} (a) + \sum_{r = 1}^{\infty} b_{r} (a) c o s r θ

(10)

可近似为

ν (a, θ) = \sum_{r = 0}^{3} b_{2 r} (a) c o s (2 r θ)

(11)

式中

b_{0} = [ω_{0}^{2} + 3 α_{0} A^{2} {/ 4]}^{1 / 2} (1 - η^{2} / 16)

b_{2} = [ω_{0}^{2} + 3 α_{0} A^{2} {/ 4]}^{1 / 2} (η / 2 + 3 η^{3} / 64)

b_{4} = [ω_{0}^{2} + 3 α_{0} A^{2} {/ 4]}^{1 / 2} (- η^{2} / 16)

b_{6} = [ω_{0}^{2} + 3 α_{0} A^{2} {/ 4]}^{1 / 2} (η^{3} / 64)

η = α_{0} A^{2} / (4 ω_{0}^{2} + 3 α_{0} A^{2}) \leq 1 / 3

(12)

α₀, ω₀为常数.利用式(12)前两个式子，通过数值迭代可求得b₀和b₂，进而可求得b₄和b₆.

2 平均伊藤方程

由于ε为小量，系统(1)的解在平衡点附近有周期解，表示为

X (t) = A (t) c o s Θ (t)

Y (t) = \dot{X} (t) = - A (t) ν (A (t), Θ (t)) s i n Θ (t)

(13)

式中

Θ (t) = Φ (t) + Γ (t)

ν (A, Θ) = \frac{d Φ}{d t} = \sqrt[]{\frac{2 [U (A) - U (A c o s Θ)]}{A^{2} s i n^{2} Θ}}

(14)

这里的A(t), Θ(t), Φ(t), Γ(t) 都是随机过程.

将方程(13)代入系统方程(1)，可解得关于A, Θ 和Γ 的随机微分方程

\frac{d A}{d t} = ε F_{11} (A, Θ) + ε F_{12} (A, Θ, u) + ε^{1 / 2} \sum_{k = 1}^{3} G_{1 k}^{'} (A, Θ) W_{i} (t)

\begin{array}{l} \frac{d Γ}{d t} = ε F_{21} (A, Θ) + ε F_{22} (A, Θ, u) \\ + ε^{1 / 2} \sum_{k = 1}^{3} G_{2 k}^{'} (A, Θ) W_{i} (t) \end{array}

\begin{array}{l} \frac{d Θ}{d t} = ω (A) + ε F_{21} (A, Θ) + ε F_{22} (A, Θ, u) \\ + ε^{1 / 2} \sum_{k = 1}^{3} G_{2 k}^{'} (A, Θ) W_{i} (t) \end{array}

(15)

其中

\begin{array}{l} F_{11} = - \frac{A^{2} ν^{2} s i n^{2} Θ}{g (A)} f (A c o s Θ, - A v s i n Θ) \\ F_{12} = - \frac{- A ν s i n Θ}{g (A)} u (t) \end{array}

F_{21} = - \frac{A ν^{2} c o s Θ s i n Θ}{g (A)} f (A c o s Θ, - A v s i n Θ)

F_{22} = - \frac{ν c o s Θ}{g (A)} u (t)

G_{1 k}^{'} = \frac{- A ν s i n Θ}{g (A)} h (A c o s Θ, - A v s i n Θ)

G_{2 k}^{'} = \frac{- ν c o s Θ}{g (A)} h (A c o s Θ, - A v s i n Θ)

(16)

由方程(15)可知，A(t)和Γ(t)为慢变量，Θ(t)为快变量.当τ 较小时，存在如下近似关系：

A (t - τ) \approx A (t) = A

，

Γ (t - τ) \approx Γ (t) = Γ

，

Θ (t - τ) \approx Θ (t) - ω (A) τ

，

\frac{d Θ (t - τ)}{d t} \approx ω (A)

(17)

应用方程(13)和(17)，可得如下近似关系：

X (t - τ) = A c o s Θ c o s ω (A) τ + A s i n Θ s i n ω (A) τ

\begin{array}{l} \dot{X} (t - τ) = - A ω (A) s i n Θ c o s ω (A) τ \\ + A ω (A) c o s Θ s i n ω (A) τ \end{array}

(18)

根据广义谐波平衡技术，u(t)可解耦为幅值依赖的等效拟线性阻尼力和拟线性回复力

ε u (t) = ε C (A) \dot{X} (t) + ε K (A) X (t)

(19)

其中

\begin{array}{l} ε C (A) = - ε k_{1} {b_{0}}^{- λ - 1} s i n (\frac{λ π}{2}) \\ + ε k_{2} {b_{0}}^{μ - 1} s i n (\frac{μ π}{2}) \end{array}

(20)

\begin{array}{l} ε K (A) = ε k_{0} + ε k_{1} b_{0}^{- λ} c o s (\frac{λ π}{2}) \\ + ε k_{2} b_{0}^{μ} c o s (\frac{μ π}{2}) \end{array}

(21)

εC(A)和εK(A)的详细计算过程见附录A.

将方程(19)代入系统方程(1)，得到与系统方程(1)等价的系统

\begin{array}{l} \ddot{X} (t) + ε [f (X, \dot{X}) - C (A)] \dot{X} (t) + g (X) \\ - ε K (A) \end{array}

= ε^{1 / 2} h (X, \dot{X}) W_{i} (t)

(22)

与系统(22)相应的能量为

H = {\dot{X}}^{2} / 2 - ε K (A) X^{2} / 2 + \int_{0}^{X} g (x) d x

(23)

相应地，方程(15)可表述为

\frac{d A}{d t} = ε F_{1} (A, Θ) + ε^{1 / 2} \sum_{k = 1}^{l} G_{1 k} (A, Θ) W_{k} (t)

\frac{d Θ}{d t} = ω (A) + ε F_{2} (A, Θ) + ε^{1 / 2} \sum_{k = 1}^{l} G_{2 k} (A, Θ) W_{k} (t)

(24)

其中

\begin{array}{l} F_{1} = - \frac{A^{2} ν^{2} s i n^{2} Θ}{ε K (A) A + g (A)} \\ \times ε [f (A c o s Θ, - A v s i n Θ) + C (A)] \end{array}

\begin{array}{l} F_{2} = - \frac{A ν^{2} c o s Θ s i n Θ}{ε K (A) A + g (A)} \\ \times ε [f (A c o s Θ, - A v s i n Θ) + C (A)] \end{array}

G_{1 k} = \frac{- A ν s i n Θ}{ε K (A) A + g (A)} f_{k} (A c o s Θ, - A v s i n Θ)

G_{2 k} = \frac{- ν c o s Θ}{ε K (A) A + g (A)} \times f_{k} (A c o s Θ, - A v s i n Θ)

(25)

基于Khasminskii极限定理^[

33]，当

ε \to 0

时，在

ε^{- 1}

量级时间区间上，A弱收敛于一维扩散过程.对Θ 作0到2π的平均，得

d A = {\bar{m}}_{1} (A) d t + {\bar{σ}}_{1} (A) d B (t)

(26)

其中平均漂移和扩散系数分别为

{\bar{m}}_{1} (A) = ε {<F_{11} (A, Θ) + F_{12} (A, Θ)>}_{Θ}

+ {<\int_{- \infty}^{0} ({\frac{\partial G_{1 k}}{\partial A}|}_{t} G_{1 l} |_{t + τ} R_{k l} (τ) + {\frac{\partial G_{1 k}}{\partial Θ}|}_{t} G_{2 l} |_{t + τ} R_{k l} (τ)) d τ>}_{Θ}

{\bar{σ}}_{1} (A) = ε {<\int_{- \infty}^{0} G_{1 k} |_{t} G_{1 l} |_{t + τ} R_{k l} (τ) d τ>}_{Θ}

(27)

式中 $< • >_{Θ}$ 表示对 $Θ$ 作平均

< • >_{Θ} = \frac{1}{2 π} \int_{0}^{2 π} <•> d Θ

(28)

S_{k l} (ω) = \frac{1}{π} \int_{- \infty}^{0} R_{k l} (τ) c o s ω τ d τ

I_{k l} (ω) = \frac{1}{π} \int_{- \infty}^{0} R_{k l} (τ) s i n ω τ d τ

(29)

$R_{k l} (τ)$ 表示系统激励的相关函数，即

E [W_{k} (t) W_{l} (t + ι)] = R_{k l} (τ)

(30)

现将 $G_{i k}$ 展开成Fourier级数

G_{i k} = G_{i k 0} (A) + \sum_{r = 1}^{\infty} [G_{i r k}^{(c)} (A) c o s r Θ + G_{i r k}^{(s)} (A) s i n r Θ]

(31)

式(31)代入式(27)，完成对τ的积分和对θ的平均后，得

\begin{array}{l} {\bar{m}}_{1} (A) = F_{10} (A) + π \frac{d G_{1 k 0}}{d A} G_{1 l 0} S_{k l} (0) \\ + \frac{π}{2} \sum_{r = 1}^{\infty} {[\frac{d G_{1 k r}^{(c)}}{d A} G_{1 l r}^{(c)} + \frac{d G_{1 k r}^{(s)}}{d A} G_{1 l r}^{(s)} + r (G_{1 k r}^{(s)} G_{2 l r}^{(c)} \\ - G_{1 k r}^{(c)} G_{2 l r}^{(s)})] \times S_{k l} (r ω (A)) + [\frac{d G_{1 k r}^{(c)}}{d A} G_{1 l r}^{(s)} - \frac{d G_{1 k r}^{(s)}}{d A} G_{1 l r}^{(c)} \\ + r (G_{1 k r}^{(s)} G_{2 l r}^{(s)} + G_{1 k r}^{(c)} G_{2 l r}^{(c)})] \times I_{k l} (r ω (A))} \end{array}

\begin{array}{l} {\bar{σ}}_{1} (A) = 2 π G_{1 k 0} G_{1 l 0} S_{k l} (0) \\ + π \sum_{r = 1}^{\infty} [(G_{1 k r}^{(c)} G_{1 l r}^{(c)} + G_{1 k r}^{(s)} G_{1 l r}^{(s)}) \times S_{k l} (r ω (A)) \\ + (G_{1 k r}^{(c)} G_{1 l r}^{(s)} - G_{1 k r}^{(s)} G_{1 l r}^{(c)}) I_{k l} (r ω (A))] \end{array}

(32)

3 稳态概率密度函数

与伊藤随机微分方程相应的FPK方程

\frac{\partial p}{\partial t} = - \frac{\partial ({\bar{m}}_{1} (A) p)}{\partial A} + \frac{1}{2} \frac{\partial^{2} ({\bar{σ}}_{1}^{2} (A) p)}{\partial A^{2}}

(33)

这里 $p = p (A, t | A_{0})$ 是关于幅值的转移概率密度，平均FPK方程的初始条件为

p (A, 0) = δ (A - A_{0})

(34)

相应的边界条件为

p = f i n i t e

，当

A = 0

(35)

p, \partial p / \partial A \to 0

，当

A \to \infty

(36)

在上述边界条件下，方程(32)具有如下形式的平稳解

p (A) = \frac{C_{0}}{{\bar{σ}}_{1}^{2} (A)} e x p [\int_{0}^{a} \frac{2 {\bar{m}}_{1} (u)}{{\bar{σ}}_{1}^{2} (u)} d u]

(37)

其中C₀为归一化常数.

相应地，由 $p (A)$ 可以得到关于广义位移和速度的联合概率密度

p (x, y) = p (x, \dot{x}) = {\frac{p (A)}{g (A)} \frac{ω (A)}{2 π}|}_{A = V^{- 1} (H)}

(38)

这里 $A = V^{- 1} (H)$ 是 $H = H (A)$ 的反函数.

4 算例

考虑宽带噪声激励下带有分数阶PI^λD^μ控制的Duffing-van del Pol振子.其运动方程为

\ddot{X} (t) + (- β_{1} + β_{2} X^{2}) \dot{X} (t) + ω_{0}^{2} X + α_{0} X^{3}

= W_{1} (t) + X W_{2} (t) + \dot{X} W_{3} (t) + ε u

(39)

其中β₁, β₂, ω₀, α₀为常数.

于是，与系统(39)等效的非线性随机系统为

\ddot{X} (t) + (- β_{1} + β_{2} X^{2} - ε C (A)) \dot{X} + (ω_{0}^{2} - ε K (A)) X + α_{0} X^{3} = W_{1} (t) + X W_{2} (t) + \dot{X} W_{3} (t)

+ α_{0} X^{3} = W_{1} (t) + X W_{2} (t) + \dot{X} W_{3} (t)

(40)

系统(39)的能量H为

H = {\dot{X}}^{2} / 2 + (ω_{0}^{2} - ε K (A)) X^{2} / 2 + α_{0} X^{4} / 4

(41)

系统(40)对应的平均伊藤方程形如式(26)，其中平均漂移系数 ${\bar{m}}_{1} (A)$ 和平均扩散系数 ${\bar{σ}}_{1}^{2} (A)$ 见附录B.

方程(40)的FPK方程可由(33)式求得，稳态概率密度函数 $p (x, \dot{x})$ 可由式(38)得到.图1-2分别为未控与受控时关于广义位移与速度的联合概率密度函数，其中x轴表示广义位移，y轴表示广义速度.左图(a)为近似解析结果，右图(b)为模拟结果.由图可知，近似解析解与蒙特卡罗模拟结果都吻合较好，这表明本文所提方法的有效性.另外，系统未控情形时(见图1)，响应表现为一个扩散的极限环.当加上控制时，由图2知，响应都得到显著减小，极限环消失，变为一个单峰，说明了分数阶控制器可以达到较好的控制效果.

图1 未控系统关于广义位移与速度的联合概率密度. 其他参数为β₁=0.01, β₂=0.02, ω₀=1.0, α₀=1.0, ω₁=ω₂=ω₃=5.0, ξ₁ =ξ₂ =ξ₃ =0.5,εk₀=εk₁=εk₂=0.0, D₁=D₂=0.2, D₃=0.1

Fig.1 Joint probability density of generalized displacement and velocity of uncontrolled system. The other parameters are β₁=0.01, β₂=0.02, ω₀=1.0, α₀=1.0,ω₁=ω₂=ω₃=5.0, ξ₁ =ξ₂ =ξ₃ =0.5, εk₀=εk₁=0.0，εk₂=0.0, D₁=D₂= 0.2, D₃=0.1

图2 λ=μ=0.35, εk₀=εk₂=-0.1, εk₁=0.1时受控系统关于广义位移与速度的联合概率密度. 其他参数同图1

Fig.2 Joint probability density of generalized displacement and velocity of controlled system with λ=μ=0.35 and εk₀=εk₂=-0.1, εk₁=0.1. The other parameters are the same as those in Fig.1

另外，本文还分别考察了宽带激励参数ξ_i和ω,_i以及激励强度D_i变化时本文所提方法的适用性,以及分数阶PI^λD^μ控制器对系统响应的控制效果(分别见图3-8).类似图1和图2，图3(a),5(a)与7(a)表示近似解析结果，图3(b), 5(b)与7(b)表示蒙特卡罗模拟结果.由这些图可知，近似解析解与蒙特卡罗模拟结果仍然吻合较好.同时，分数阶PI^λD^μ控制器仍具有较好的适应性，消除了极限环，成功实现了对系统响应的控制.

图3 ξ_i =0.3时未控系统关于广义位移与速度的联合概率密度.其他参数同图1

Fig.3 Joint probability density of generalized displacement and velocity of uncontrolled system with ξ_i=0.3. The other parameters are the same as those in Fig.1

图4 ξ_i =0.3时受控系统关于广义位移与速度的联合概率密度.其他参数为同图2

Fig.4 Joint probability density of generalized displacement and velocity of controlled system with ξ_i=0.3. The other parameters are the same as those in Fig.2

图5 ω_i=2.5时未控系统关于广义位移与速度的联合概率密度.系统其他参数同图1.

Fig.5 Joint probability density of generalized displacement and velocity of uncontrolled system with ω_i=2.5. The other parameters are the same as those in Fig.1

图6 ω_i=2.5受控系统关于广义位移与速度的联合概率密度. 其他参数同图2.

Fig.6 Joint probability density of generalized displacement and velocity of controlled system with ω_i=2.5. The other parameters are the same as those in Fig.2

图7 D₁=D₂=0.4, D₃=0.2时未控系统关于广义位移与速度的联合概率密度. 系统其他参数同图1

Fig.7 Joint probability density of generalized displacement and velocity of uncontrolled system with D₁=D₂=0.4, D₃=0.2. The other parameters are the same as those in Fig.1

图8 D₁=D₂=0.4, D₃=0.2时受控系统关于广义位移与速度的联合概率密度. 系统其他参数同图2

Fig.8 Joint probability density of generalized displacement and velocity of controlled system with D₁=D₂=0.4, D₃=0.2. The other parameters are the same as those in Fig.2

5 结论

本文运用随机平均法研究了基于分数阶PI^λD^μ控制器的宽带噪声激励下的强非线性系统.基于广义谐波平衡技术，将分数阶控制力表示为幅值依赖的等效拟线性阻尼力和拟线性回复力，再用广义谐和函数将等效非线性系统转化为关于幅值的平均伊藤随机微分方程，建立并求解相应的FPK方程.在求解伊藤随机微分方程时，运用随机平均法得到了系统平均漂移系数、扩散系数以及平稳概率密度的闭合表达式，简化了求解过程，尤其是将分数阶控制力等效为整数阶控制力，克服了分数阶微积分模拟耗时长的问题.数值结果表明近似解析解与原方程的蒙特卡罗模拟结果吻合较好，说明随机平均法的适用性和高精度.分数阶PI^λD^μ控制器对系统响应进行了有效控制.另外，宽带随机激励参数发生变化时，文中所提出的方法同样具有较好的适应性，近似解析解与蒙特卡罗模拟结果十分吻合，分数阶控制仍同样具有非常好的控制效果.

附录

附录A

\begin{array}{l} ε u = ε k_{0} X + ε k_{1} I^{λ} X (t) + ε k_{2} D^{μ} X (t) = ε k_{0} X + \frac{ε k_{1}}{Γ (λ)} \int_{0}^{t} \frac{X (t - τ)}{τ^{1 - λ}} d τ + \frac{ε k_{2}}{Γ (1 - μ)} [\frac{X (0)}{t^{μ}} + \int_{0}^{t} \frac{\dot{X} (t - τ)}{τ^{μ}} d τ] \\ = ε k_{0} A (t) c o s Θ (t) + \frac{ε k_{1}}{Γ (λ)} A (t) \times [c o s Θ (t) \int_{0}^{t} \frac{c o s ω τ}{τ^{1 - λ}} d τ + s i n Θ (t) \int_{0}^{t} \frac{s i n ω τ}{τ^{1 - λ}} d τ] + \frac{ε k_{2}}{Γ (1 - μ)} {\frac{X (0)}{t^{μ}} \\ - A (t) ω [s i n Θ (t) \int_{0}^{t} \frac{c o s ω τ}{τ^{μ}} d τ - c o s Θ (t) \int_{0}^{t} \frac{s i n ω τ}{τ^{μ}} d τ]} \end{array}

(A-1)

为进一步简化，引入如下渐进积分

\int_{0}^{t} \frac{c o s (ω τ)}{τ^{α}} d τ = ω^{α - 1} \int_{0}^{s} \frac{c o s (u)}{u^{α}} d u = ω^{α - 1} [Γ (1 - α) s i n (\frac{α π}{2}) + \frac{s i n (s)}{s^{α}} + Ο (s^{- α - 1})]

\int_{0}^{t} \frac{s i n (ω τ)}{τ^{α}} d τ = ω^{α - 1} \int_{0}^{s} \frac{s i n (u)}{u^{α}} d u = ω^{α - 1} [Γ (1 - α) c o s (\frac{α π}{2}) - \frac{c o s (s)}{s^{α}} + Ο (s^{- α - 1})]

(A-2)

其中, $u = ω τ, s = ω t$

可得

ε u = (- ε k_{1} {b_{0}}^{- λ - 1} s i n (\frac{λ π}{2}) + ε k_{2} {b_{0}}^{μ - 1} s i n (\frac{μ π}{2})) \dot{X} (t) + (ε k_{0} + ε k_{1} b_{0}^{- λ} c o s (\frac{λ π}{2}) + ε k_{2} b_{0}^{μ} c o s (\frac{μ π}{2})) X (t)

(A-3)

其中

ε C (A) = = - ε k_{1} {b_{0}}^{- λ - 1} s i n (\frac{λ π}{2}) + ε k_{2} {b_{0}}^{μ - 1} s i n (\frac{μ π}{2})

ε K (A) = ε k_{0} + ε k_{1} b_{0}^{- λ} c o s (\frac{λ π}{2}) + ε k_{2} b_{0}^{μ} c o s (\frac{μ π}{2})

(A-4)

附录B

算例一得到的平均漂移系数和平均扩散系数

{\bar{m}}_{1} (A) = \frac{β_{1} (A ω_{0}^{2} + K (A) A + 0.625 α_{0} A^{3})}{2 [ω_{0}^{2} + K (A) + α_{0} A^{2}]} - \frac{β_{2} (A^{3} ω_{0}^{2} + K (A) A^{3} + 0.75 α_{0} A^{5})}{8 [ω_{0}^{2} + K (A) + α_{0} A^{2}]}

\begin{array}{l} + \frac{(A ω_{0}^{2} + K (A) A + 0.625 α_{0} A^{3})}{2 [ω_{0}^{2} + K (A) + α_{0} A^{2}]} \times (- ε k_{1} {b_{0}}^{- λ - 1} s i n (\frac{λ π}{2}) + ε k_{2} {b_{0}}^{μ - 1} s i n (\frac{μ π}{2})) \\ + b_{6} S_{2} (8 ω) \frac{d}{d A} [\frac{A b_{6}}{ω_{0}^{2} + K (A) + α_{0} A^{2}}]\} + \frac{π A}{16 {[ω_{0}^{2} + K (A) + α_{0} A^{2}]}^{2}} \end{array}

\times [(2 b_{0} - b_{4}) (2 b_{0} + 2 b_{2} + b_{4}) \times S_{2} (2 ω) + 2 (b_{2} - b_{6}) (b_{2} + 2 b_{4} + b_{6}) S_{2} (4 ω)

\begin{array}{l} + 3 b_{4} (b_{4} + 2 b_{6}) S_{2} (6 ω) + 4 b_{6}^{2} S_{2} (8 ω)] + \frac{π A}{8 {[ω_{0}^{2} + K (A) + α_{0} A^{2}]}^{2}} {(2 b_{0} - b_{2}) S_{1} (ω) \\ \times \frac{d}{d A} [\frac{2 b_{0} - b_{2}}{ω_{0}^{2} + K (A) + α_{0} A^{2}}] + (b_{2} - b_{4}) S_{1} (3 ω) \frac{d}{d A} [\frac{b_{2} - b_{4}}{ω_{0}^{2} + K (A) + α_{0} A^{2}}] \\ + (b_{4} - b_{6}) S_{1} (5 ω) \frac{d}{d A} [\frac{b_{4} - b_{6}}{ω_{0}^{2} + K (A) + α_{0} A^{2}}] + b_{6} S_{1} (7 ω) \frac{d}{d A} [\frac{b_{6}}{ω_{0}^{2} + K (A) + α_{0} A^{2}}] \\ + \frac{π}{8 A {[ω_{0}^{2} + K (A) + α_{0} A^{2}]}^{2}} \times [(4 b_{0}^{2} - b_{2}^{2}) S_{1} (ω) + 3 (b_{2}^{2} - b_{4}^{2}) S_{1} (3 ω) + 5 (4 b_{4}^{2} - b_{6}^{2}) S_{1} (5 ω) \\ + 7 b_{6}^{2} S_{1} (7 ω)] + \frac{π A}{256 {[ω_{0}^{2} + K (A) + α_{0} A^{2}]}^{3}} (5 A^{2} α_{0} + 8 ω_{0}^{2}) \times (5 A^{4} {α_{0}}^{2} + 7 A^{2} α_{0} ω_{0}^{2} + 8 ω_{0}^{4}) S_{3} (0) \end{array}

\begin{array}{l} + \frac{π A^{3} (A^{2} α_{0} + 2 ω_{0}^{2}) (7 A^{2} α_{0} + 8 ω_{0}^{2})}{32 {[ω_{0}^{2} + K (A) + α_{0} A^{2}]}^{2}} \times S_{3} (2 ω) + \frac{π A^{7} {α_{0}}^{2} S_{3} (4 ω)}{64 {[ω_{0}^{2} + K (A) + α_{0} A^{2}]}^{2}} \end{array}

(B-1)

\begin{array}{l} {\bar{σ}}_{1}^{2} (A) = \frac{π A^{2}}{16 {[ω_{0}^{2} + K (A) + α_{0} A^{2}]}^{2}} [{(2 b_{0} - b_{4})}^{2} S_{2} (2 ω) + {(b_{2} - b_{6})}^{2} S_{2} (4 ω) + b_{4}^{2} S_{2} (6 ω) + b_{6}^{2} S_{2} (8 ω)] \\ + \frac{π}{4 {[ω_{0}^{2} + K (A) + α_{0} A^{2}]}^{2}} [{(2 b_{0} - b_{2})}^{2} S_{1} (ω) + {(b_{2} - b_{4})}^{2} S_{1} (3 ω) + {(b_{4} - b_{6})}^{2} S_{1} (5 ω) \end{array}

(B-2)

参考文献

Tong D K,Wang R H.Analysis of the flow of non-Newtonian viscoelastic fluids in fractal reservoir with the fractional derivative.Science in China Series G Physics Mechanics & Astronomy, 2004, 47(4): 424~441 [百度学术]

Yu S W, Huang T Z, Liu X Y, et al. Information measures based on fractional calculus. Information Processing Letters, 2012, 112 (23): 916~921 [百度学术]

Ionescu C, Lopes A, Copot D, et al. The role of fractional calculus in modeling biological phenomena: A review. Communications in Nonlinear Science and Numerical Simulation, 2017, 51: 141~159 [百度学术]

Clara M I, James F K. Fractional calculus for respiratory mechanics: Power law impedance, viscoelasticity, and tissue heterogeneity. Chaos, Solitons & Fractals, 2017, 102: 433~440 [百度学术]

Zhao D Z, Pan X Q, Luo M K. A new framework for multivariate general conformable fractional calculus and potential applications. Physica A: Statistical Mechanics and its Applications, 2018, 510: 271~280 [百度学术]

Tarasov v E, Tarasova V V. Macroeconomic models with long dynamic memory: Fractional calculus approach. Applied Mathematics and Computation, 2018, 338: 466~486 [百度学术]

Zhu W Q, Ying Z G. Optimal nonlinear feedback control of quasi Hamiltonian systems. Science in China Series A: Mathematics, 1999, 42: 1213~1219 [百度学术]

Zhu W Q, Ying Z G, Soong T T. An optimal nonlinear feedback control strategy for stochastically excited structural systems. Nonlinear Dynamics, 2001, 24: 31~51 [百度学术]

Zhu W Q, Huang Z L, Yang Y Q. Stochastic averaging of quasi integrable Hamiltonian systems. Journal of Applied Mechanics. 1997, 64: 975~984 [百度学术]

Huang Z L, Zhu W Q. Stochastic averaging of quasi-integrable Hamiltonian systems under bounded noise excitations. Probabilistic Engineering Mechanics, 2004, 19(3): 219~228 [百度学术]

Yong J M, Zhou X Y. Stochastic Control, Hamiltonian Systems and HJB Equations. New York: Springer, 1999 [百度学术]

Feng C S, Wu Y J, Zhu W Q. Response of Duffing system with delayed feedback control under combined harmonic and real noise excitations. Communication in Nonlinear Science and Numerical Simulation, 2009, 14: 2542~2250 [百度学术]

Feng C S, Liu R. Response of Duffing system with delayed feedback control under combined harmonic and real noise excitations. Archive of Applied Mechanics, 2012, 82(12): 1753~1761 [百度学术]

Xu Y, Ma S, Zhang H. Hopf bifurcation control for stochastic dynamical system with nonlinear random feedback method. Nonlinear Dynamics, 2011, 65(1): 77~84 [百度学术]

Oustaloup A, Sabatier J, Lanusse P. From fractal robustness to CRONE control. Fractional Calculus and Applied Analysis, 1999, 2: 1~30 [百度学术]

Podlubny I. Fractional-order systems and PI^λD^μcontroller. IEEE Transactions on Automatic Control, 1999, 44: 208~214 [百度学术]

Lurie B J. Three-parameter tunable tilt-integral-derivative controller. U.S. Patent Documents, 2004, 11(2): 125~142 [百度学术]

汪纪锋, 李元凯. 分数阶P(ID)^u控制器和超前滞后校正器的设计. 电路与系统学报, 2006, 11: 21~25(Wang J F, Li Y K. Design for two types of fractional-order controllers: P(ID) controller and lead-lag compensator. Journal of Circuits and Systems^u, 2006, 11: 21~25(in Chinese)) [百度学术]

Chen Y Q, Moore K L. Discretization schemes for fractional-order differentiators and integrators. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 2002, 49(3): 363~367 [百度学术]

Chen L C, Zhao T L, Li W, et al. Bifurcation control of bounded noise excited Duffing oscillator by a weakly fractional-order PID feedback controller. Nonlinear Dynamics, 2016, 83(1-2): 529~539 [百度学术]

Grigoriu M. Control of time delay linear systems with Gaussian white noise. Probabilistic Engineering Mechanics, 1997, 12: 89~96 [百度学术]

Shah P, Agashe S. Review of fractional PID controller. Mechatronics, 2016, 38: 29~41 [百度学术]

Chen Z H, Yuan X H, Ji B, et al. Design of a fractional order PID controller for hydraulic turbine regulating system using chaotic non-dominated sorting genetic algorithm II. Energy Conversion and Management, 2014, 84: 390~404 [百度学术]

Zamani A, Barakati S M, Yousofi-Darmian S. Design of a fractional order PID controller using GBMO algorithm for load-frequency control with governor saturation consideration. ISA Transactions, 2016, 64:56~66 [百度学术]

Ates A, Yeroglu C. Optimal fractional order PID design via Tabu search based algorithm. ISA Transactions, 2016, 60: 109~118 [百度学术]

Chen N, Chen Y. On fractional control method for four-wheel-steering vehicle. Science in China Series E: Technological Sciences, 2009, 52(3): 603~609 [百度学术]

Xu Y H, Zhou J Z, Xue X M, et al. An adaptively fast fuzzy fractional order PID control for pumped storage hydro unit using improved gravitational search algorithm. Energy Conversion and Management, 2016, 111: 67~78 [百度学术]

Nucera F, McFarland D M, Bergman L A, et al. Application of broadband nonlinear targeted energy transfers for seismic mitigation of a shear frame: Computational results. Journal of Sound and Vibration, 2010, 329:2973~2994 [百度学术]

Ates A，Yeroglu C. Optimal fractional order PID design via Tabu Search based algorithm. ISA Transactions, 2016, 60: 109~118 [百度学术]

Zhong G, Deng H, Li J. Chatting-free variable structure controller design via fractional calculus approach and its application. Nonlinear Dynamics, 2015, 81(1): 679~694 [百度学术]

Li W, Chen L C, Zhao J F, et al. Reliability estimation of stochastic dynamical systems with fractional order PID controller. International Journal of Structural Stability and Dynamics, 2018, (18): 1850083-1 [百度学术]

Chen L C, Liang X, Zhu W Q, et al. Stochastic averaging technique for SDOF strongly nonlinear systems with delayed feedback fractional-order PD controller. Science in China Series E: Technological Sciences, 2019, 62: 287~297 [百度学术]

Chen L C, Zhao T L, Li W, et al. Bifurcation control of bounded noise excited Duffing oscillator by a weakly fractional-order PI^λD^μfeedback controller. Nonlinear Dynamics, 2016, 83: 529~539 [百度学术]

Huang Z L, Zhu W Q, Suzuki Y. Stochastic averaging of strongly nonlinear oscillators under combined harmonic and white noise excitations. Journal of Sound and Vibration, 2000, 238(2): 233~256 [百度学术]

Khasminskii R Z. On the averaging principle for Itô stochastic differential equations. Kibernetka, 1986, 3: 260~279 [百度学术]