分数阶微积分已有300多年的历史，自然界的物理现象大多以分数阶的形式存在，整数阶微分方程正好是分数阶微分方程的特例.与整数阶模型相比，分数阶模型更接近真实的世界，具有更诱人的发展前景，近年来已得到了越来越多的关注^[1-2].值得注意的是，与只有固定翼混沌吸引子的分数阶混沌系统相比，具有多种多翼混沌吸引子共存的分数阶混沌系统显示出更复杂的动力学行为和更好的性能^[3-4].在安全通信^[5]和图像加密^[6]中，此类分数阶混沌系统具有更高的序列复杂度以及更大的密钥空间，提高了系统安全性能.因此，发现和构造具有多种多翼混沌吸引子共存的低维分数阶混沌系统具有更大的价值.Zhou等人基于四翼整数阶忆阻混沌系统，构造了相应的分数阶忆阻系统，出现三翼与三翼、三翼与四翼混沌吸引子的共存^[7]；Borah和Roy设计了一种新的分数阶混沌系统，此系统具有三翼和四翼混沌吸引子共存^[8]；Xian等人构造了一个双翼与四翼混沌吸引子共存的分数阶混沌系统^[9].目前，构造具有更多种多翼混沌吸引子共存的分数阶混沌系统仍然存在一定挑战.在已有文献中，能够产生双翼到四翼混沌吸引子共存的三维分数阶混沌系统还比较少见.

本文设计了一个新型三维分数阶混沌系统，对其进行了动力学分析.在q=0.98时，出现双翼、双翼、四翼等混沌吸引子的共存；q=0.83时，出现双翼、三翼、四翼等混沌吸引子的共存，表明系统具有丰富的混沌特性.与文献[10]和文献[11]相比，此系统不存在不光滑的非线性函数，更易于用硬件电路实现.对系统进行了Multisim模拟电路设计和仿真，仿真结果与数值分析相符，进一步验证了系统的混沌行为.基于分数阶Lyapunov稳定理论，设计了一个自适应同步控制器，仿真结果表明了所设计控制器的有效性.

1 系统模型及其特性 1.1 分数阶系统模型

增广Lü系统^[12]状态方程为

$ \left\{ {\begin{array}{*{20}{l}} {\dot x = - (bc/(b + c))x - yz + d,}\\ {\dot y = by + xz,}\\ {\dot z = cz + xy.} \end{array}} \right. $

(1)

将式(1)中第一个方程右边添加非线性项-kx-9yz-d，第3个方程右边添加非线性项-9xy+y².再根据Caputo分数阶积分定义，构建了一个新型三维分数阶混沌系统，系统状态方程为

$ \left\{ {\begin{array}{*{20}{l}} {{D^q}x = - ax - 10yz,}\\ {{D^q}y = bx + xz,}\\ {{D^q}z = cz - 8xy + {y^2}.} \end{array}} \right. $

(2)

式中:x，y，z为系统状态变量，a，b为系统参数，q为系统阶数.

固定参数a=2，b=1，c=-0.78，取q=0.98，其最大Lyapunov指数(MAXLE)为1.065 96.利用预估-校正数值计算方法^[13]，结合matlab软件，得到系统的仿真相图.初始值为[2, 1, 1]、[5, 1, 1]和[6, 1, 1]时，分别得到两个孤立的双翼混沌吸引子和一个四翼混沌吸引子，其x-y面和y-z面如图 1所示.其中，图 1(a)和(b)中的蓝色部分为初始值[2, 1, 1]诱发的双翼混沌吸引子，红色部分为初始值[5, 1, 1]诱发的四翼混沌吸引子.

取q=0.83，其MAXLE为1.494 3.初始值为[2, 1, 1]、[5, 2, 1]和[11, 1, 1]时，分别得到双翼、三翼和四翼混沌吸引子，其x-y面和y-z面如图 2所示.

由此可见，系统阶数q=0.98时，系统有双翼、双翼、四翼等混沌吸引子共存；q=0.83时，系统有双翼、三翼、四翼等混沌吸引子共存，表明了阶数q对系统的混沌特性具有较大影响.

1.2 平衡点及其稳定性

参数不变，令系统(2)方程式的左边等于0，即

$ \left\{ {\begin{array}{*{20}{l}} {0 = - 2x - 10yz,}\\ {0 = y + xz,}\\ {0 = - 0.78z - 8xy + {y^2}.} \end{array}} \right. $

(3)

求得平衡点，如表 1所示.

在平衡点$\left( {\begin{array}{*{20}{l}} {\hat x}&, &{\hat y, }&{\hat z} \end{array}} \right)$处线性化系统(2)，得Jacobian矩阵

$ \mathit{\boldsymbol{J}} = \left[ {\begin{array}{*{20}{c}} { - 2}&{ - 10\hat z}&{ - 10\hat y}\\ {\hat z}&1&{\hat x}\\ { - 8\hat y}&{2\hat y - 8\hat x}&{ - 0.78} \end{array}} \right]. $

令det(λE－J)=0，得到其特征多项式为

$ f\left( \lambda \right) = {\lambda ^3} + {A_2}{\lambda ^2} + {A_1}\lambda + {A_0}. $

(4)

其中

$ {A_2} = 1.78, $

$ {A_1} = 8{{\hat x}^2} + 10{{\hat z}^2} - 80{{\hat y}^2} - 2\hat x\hat y - 1.22, $

$ {A_0} = 16{{\hat x}^2} + 80{{\hat y}^2} + 7.8{{\hat z}^2} + 20{{\hat y}^2}\hat z - 4\hat x\hat y - 160\hat x\hat y\hat z - 1.56. $

若平衡点稳定，其特征值满足Re[λ] < 0.依据Routh-Hurwitz判据，当且仅当A₂>0，A₀>0，A₂A₁-A₀>0时，平衡点稳定^[14-15].求得每个平衡点所对应的特征值及平衡点的稳定性如表 2所示.

表 2表明，平衡点S₁是第一类鞍点，S₀，S₂，S₃和S₄是第二类鞍点.混沌吸引子的环围绕第二类鞍点产生，第一类鞍点起到连接环的作用^[16-17].

设λ是第二类鞍点的共轭复根，为使分数阶系统存在混沌吸引子，λ必须保持在不稳定区域.由文献[18]可得，需满足

$ \left| {\arg \left( \lambda \right)} \right| < 0.5{\rm{ \mathit{ π} }}q. $

(5)

其中arg(λ)是特征值λ的辐角，所以有

$ q > \frac{2}{{\rm{ \mathit{ π} }}}\left| {\arg \left( \lambda \right)} \right|. $

(6)

由此，a=2，b=1，c=-0.78时，可得分数阶系统(2)存在混沌吸引子的必要条件，即q>0.822 4.

1.3 时序图

参数不变，取q=0.83, 初始值为[5, 2, 1]，得到相应的时序图，如图 3(a)所示，从图中能直观地看出系统的混沌行为.

1.4 最大Lyapunov指数谱和分岔图

系统阶数的变化会改变平衡点的稳定性, 引起系统状态的变化.参数不变，选择阶数q为控制变量.采用Adams-Bashforth-Moulton数值算法的改进版本^[19]，初始值为[5, 2, 1]时，得到系统变量z随系统阶数q变化的最大Lyapunov指数谱和分岔图，如图 3(b)和图 3(c)所示.最大Lyapunov指数谱和分岔图可以比较直观地反映非线性动力学系统随阶数变化的动态特性.当最大Lyapunov指数>0的时候，系统处于混沌状态.由图 3(c)可知，当q∈ [0.5，0.819]时，系统处于周期状态；当q∈(0.819，1]时，系统处于混沌状态.q从0.5向1增加的过程中，结合图 3(c)可知，系统通向混沌的道路为倍周期分岔道路，进一步验证了q对系统具有较大影响.

2 系统的电路仿真

参数不变，取初始值为[5, 2, 1]，根据文献[20]，当q=0.83时，分数阶算子逼近

$ H\left( s \right) = \frac{{3.1869\left( {s + 0.1981} \right)(s + 5.181)(s + 87.72)}}{{(s + 0.01324)(s + 0.299)(s + 6.792)(s + 153.8)}}. $

(7)

文献[19]中的分数阶等效串并联RC模块电路单元如图 4所示.

其中R_a=44.310 6 MΩ，R_b=1.307 8 MΩ，R_c=0.085 4 MΩ，R_d=0.005 9 MΩ，C_a=1.709 7 μF，C_b=2.215 7 μF，C_c=1.298 μF，C_d=0.724 5 μF.

采用线性电阻、电容、LM2924N型运算放大器、MULTIPLIER模拟乘法器，其中乘法器的输出增益为1，从而设计出了系统的模拟电子电路，如图 5所示.运算放大器U1A相当于加法器，其输出端电压为v_U1A=－(R₃/R₂)x－(R₃/R₁)yz，U2A相当于积分器，利用分数阶算子单元模块，得到

$ {D^{0.83}}x = {v_{{\rm{U1A}}}} = - \left( {{R_3}/{R_2}} \right)x - \left( {{R_3}/{R_1}} \right)yz, $

(8)

同理

$ {D^{0.83}}x = {v_{{\rm{U}}4{\rm{A}}}} = \left( {{R_9}/{R_7}} \right)xz + \left( {{R_9}/{R_8}} \right)y, $

(9)

$ \begin{array}{*{20}{c}} {{D^{0.83}}x = {v_{{\rm{U}}6{\rm{A}}}} = - \left( {{R_{14}}/{R_{13}}} \right)z - \left( {{R_{14}}/{R_{12}}} \right)xy + }\\ {\left( {{R_{14}}/{R_{11}}} \right){y^2}.} \end{array} $

(10)

根据式(2)中的系数，a=2，b=1，c=-0.78，令R₁=100 Ω，R₂=500 Ω，R₃=R₇=R₈=R₉=1 kΩ，R₅=R₆=R₁₆=R₁₇=10 kΩ，R₄=R₁₀=R₁₅=100 kΩ，R₁₁=R₁₄=1.56 kΩ，R₁₂=195 Ω，R₁₃=2 kΩ，R_a1=R_a2=R_a3=44.310 6 MΩ，R_b1=R_b2=R_b3=1.307 8 MΩ，R_c1=R_c2=R_c3=0.085 4 MΩ，R_d1=R_d2=R_d3=0.005 9 MΩ，C_a1=C_a2=C_a3=1.709 7 μF，C_b1=C_b2=C_b3=2.215 7 μF，C_c1=C_c2=C_c3=1.298 μF，C_d1=C_d2=C_d3=0.724 5 μF.

图 6为示波器上观察到的结果，可以看出实验结果与数值仿真结果完全相符，进一步验证了系统的混沌行为.

3 系统的自适应同步

采用自适应控制方法^[21-23]实现分数阶混沌系统的同步.令驱动系统为

$ \left\{ {\begin{array}{*{20}{l}} {{D^{\rm{q}}}{x_1} = - a{x_1} - 10{x_2}{x_3},}\\ {{D^{\rm{q}}}{x_2} = b{x_2} + {x_1}{x_3},}\\ {{D^{\rm{q}}}{x_3} = c{x_3} - 8{x_1}{x_2} + x_2^2.} \end{array}} \right. $

(11)

其中x₁，x₂，x₃为系统变量.响应系统为

$ \left\{ {\begin{array}{*{20}{l}} {{D^{\rm{q}}}{y_1} = - a{y_1} - 10{y_2}{y_3} + {u_1},}\\ {{D^{\rm{q}}}{y_2} = b{y_2} + {y_1}{y_3} + {u_2},}\\ {{D^{\rm{q}}}{y_3} = c{y_3} - 8{y_1}{y_2} + y_2^2 + {u_3}.} \end{array}} \right. $

(12)

其中y₁，y₂，y₃是状态变量，u₁，u₂，u₃是自适应控制器.

令驱动系统(11)和响应系统(12)的同步误差为e_i=y_i－x_i(i=1, 2, 3)，则同步误差系统为

$ \left\{ \begin{array}{l} {D^{\rm{q}}}{e_1} = - a{e_1} - 10{e_2}{e_3} - 10{e_2}{x_3} - \\ \;\;\;\;\;\;\;\;\;10{e_3}{x_2} + {u_1},\\ {D^{\rm{q}}}{e_2} = b{e_2} + {e_1}{e_3} + {e_1}{x_3} + \\ \;\;\;\;\;\;\;\;\;{e_3}{x_1} + {u_2},\\ {D^{\rm{q}}}{e_3} = c{e_3} - 8\left( {{e_1}{x_2} + {e_2}{x_1} + {e_1}{e_2}} \right) + \\ \;\;\;\;\;\;\;\;\;\;e_2^2 + 2{e_2}{x_2} + {u_3}. \end{array} \right. $

(13)

参数估计误差定义为

$ \left\{ {\begin{array}{*{20}{l}} {{e_{\rm{a}}}(t) = a - \tilde a(t),}\\ {{e_{\rm{b}}}(t) = b - \tilde b(t),}\\ {{e_{\rm{c}}}(t) = c - \tilde c(t).} \end{array}} \right. $

(14)

对参数估计误差(14)求q阶导，得

$ \left\{ {\begin{array}{*{20}{l}} {{D^{\rm{q}}}{e_{\rm{a}}}(t) = - {D^{\rm{q}}}\tilde a(t),}\\ {{D^{\rm{q}}}{e_{\rm{b}}}(t) = - {D^{\rm{q}}}\tilde b(t),}\\ {{D^{\rm{q}}}{e_{\rm{c}}}(t) = - {D^{\rm{q}}}\tilde c(t).} \end{array}} \right. $

(15)

设计系统的自适应控制器为

$ \left\{ \begin{array}{l} {u_1} = \tilde a\left( t \right){e_1} + 18{e_3}{x_2} + 17{e_2}{x_3} + \\ \;\;\;\;\;\;9{e_2}{x_3} - {k_1}{e_1},\\ {u_2} = - \tilde b\left( t \right){e_2} + 7{e_3}{x_1} - {e_2}{e_3} - \\ \;\;\;\;\;\;2{e_3}{x_2} - {k_2}{e_2},\\ {u_3} = - \tilde c\left( t \right){e_3} - {k_3}{e_3}. \end{array} \right. $

(16)

式中k₁, k₂, k₃为正的增益常数.

定义参数更新定律为

$ \left\{ {\begin{array}{*{20}{l}} {{D^{\rm{q}}}\tilde a\left( t \right) = {\lambda _1}\left( { - e_1^2 + {e_{\rm{a}}}} \right),}\\ {{D^{\rm{q}}}\tilde b\left( t \right) = {\lambda _2}\left( {e_2^2 + {e_{\rm{b}}}} \right),}\\ {{D^{\rm{q}}}\tilde c\left( t \right) = {\lambda _3}\left( {e_3^2 + {e_{\rm{c}}}} \right).} \end{array}} \right. $

(17)

式中λ₁, λ₂, λ₃为正的增益常数.

对于分数阶混沌系统，由参考文献[24]给出的分数阶Lyapunov稳定性理论，得出引理1.

引理1  对于分数阶系统${D^{\rm{q}}}\mathit{\boldsymbol{X}}(t) = f(\mathit{\boldsymbol{X}})$，其中X=(x₁, x₂, …, x_n)^T为系统状态变量.当系统阶数q∈(0, 1)时，如果存在实对称正定矩阵P，使对任意状态变量矩阵X，恒有$J = \mathit{\boldsymbol{XP}}{D^9}\mathit{\boldsymbol{X}} \le 0$成立，则分数阶系统渐近稳定.

定理1  在自适应控制器(16)和参数更新定律(17)的作用下，同步误差系统(13)渐近稳定.

证明  构造J函数如下:

$ J = \mathit{\boldsymbol{eP}}{D^{\rm{q}}}e = \left[ {e,{e_{\rm{a}}},{e_{\rm{b}}},{e_{\rm{c}}}} \right]\mathit{\boldsymbol{P}}{\left[ {{D^{\rm{q}}}e,{D^{\rm{q}}}{e_{\rm{a}}},{D^{\rm{q}}}{e_{\rm{b}}},{D^{\rm{q}}}{e_{\rm{c}}}} \right]^{\rm{T}}} $

其中，e为同步误差，e_a, e_b, e_c为参数估计误差，P为实对称正定矩阵，选择正定矩阵P=diag(1, 1, 1, 1/λ₁, 1/λ₂, 1/λ₃)，λ₁, λ₂, λ₃为特征值.

将同步误差系统(13)、参数估计误差(14)、自适应控制器(16)和参数更新定律(17)代入上式，简化可得

$ \begin{array}{l} J = \left[ {e,{e_{\rm{a}}},{e_{\rm{b}}},{e_{\rm{c}}}} \right]\mathit{\boldsymbol{P}}{\left[ {{D^{\rm{q}}}e,{D^{\rm{q}}}{e_{\rm{a}}},{D^{\rm{q}}}{e_{\rm{b}}},{D^{\rm{q}}}{e_{\rm{c}}}} \right]^{\rm{T}}} = \\ \;\;\;\;\;{e_1}{D^{\rm{q}}}{e_1} + {e_2}{D^{\rm{q}}}{e_2} + {e_3}{D^{\rm{q}}}{e_3} + \left( {1/{\lambda _1}} \right){e_{\rm{a}}}{D^{\rm{q}}}{e_{\rm{a}}} + \\ \;\;\;\;\;\left( {1/{\lambda _2}} \right){e_{\rm{b}}}{D^{\rm{q}}}{e_{\rm{b}}} + \left( {1/{\lambda _3}} \right){e_{\rm{c}}}{D^{\rm{q}}}{e_{\rm{c}}} = \\ \;\;\;\;\; - ae_1^2 - 10{e_1}{e_2}{e_3} - 10{e_1}{e_2}{x_3} - 10{e_1}{e_3}{x_2} + \tilde a\left( t \right)e_1^2 - \\ \;\;\;\;\;{k_1}e_1^2 + 17{e_1}{e_2}{e_3} + 9{e_1}{e_2}{x_3} + 18{e_1}{e_3}{x_3} + be_2^2 + \\ \;\;\;\;\;{e_1}{e_2}{e_3} + {e_1}{e_2}{x_3} + {e_2}{e_3}{x_1} - \tilde b\left( t \right)e_2^2 - {e_2}{e_3}{x_1} + \\ \;\;\;\;\;8{e_2}{e_3}{x_1} - e_2^2{e_3} - 2{e_2}{e_3}{x_2} - {k_2}e_2^2 - ce_3^2 - \\ \;\;\;\;\;8{e_1}{e_2}{e_3} - 8{e_1}{e_3}{x_2} - 8{e_2}{e_3}{x_1} + e_2^2{e_3} + \\ \;\;\;\;\;2{e_2}{x_2}{e_3} - \tilde c\left( t \right)e_3^2 - {k_3}e_3^2 + \\ \;\;\;\;\;\left( {1/{\lambda _1}} \right){e_{\rm{a}}}{D^{\rm{q}}}{e_{\rm{a}}} + \left( {1/{\lambda _2}} \right){e_{\rm{b}}}{D^{\rm{q}}}{e_{\rm{b}}} + \left( {1/{\lambda _3}} \right){e_{\rm{c}}}{D^{\rm{q}}}{e_{\rm{c}}} = \\ \;\;\;\;\; - {e_{\rm{a}}}e_1^2 - {k_1}e_1^2 + {e_{\rm{b}}}e_2^2 - {k_2}e_2^2 + {e_{\rm{c}}}e_3^2 - {k_3}e_3^2 + \\ \;\;\;\;\;\left( {1/{\lambda _1}} \right){e_{\rm{a}}}{D^{\rm{q}}}{e_{\rm{a}}} + \left( {1/{\lambda _2}} \right){e_{\rm{b}}}{D^{\rm{q}}}{e_{\rm{b}}} + \left( {1/{\lambda _3}} \right){e_{\rm{c}}}{D^{\rm{q}}}{e_{\rm{c}}} = \\ \;\;\;\;\; - {k_1}e_1^2 - {k_2}e_2^2 - {k_3}e_3^2 - e_a^2 - e_b^2 - e_c^2. \end{array} $

显然，J≤0，据引理1误差系统渐近稳定.

证毕.

采用四阶Runge-Kutta算法，仿真步长为0.001.取q=0.83，a=2，b=1，c=-7.8，k₁=k₂=k₃=20, λ₁=λ₂=λ₃=20.x₁(0)=1，x₂(0)=2, x₃(0)=3作为驱动系统(11)的初始条件；y₁(0)=2, y₂(0)=3, y₃(0)=4作为响应系统(12)的初始条件；$\tilde a(0) = 1, \tilde b(0) = 2, \tilde c(0) = 3$作为参数估计的初始条件.

图 7描述了驱动系统(11)和响应系统(12)的同步误差e₁, e₂, e₃和对系统未知参数的识别.结果表明响应系统与驱动系统在0.2 s内达到同步，在0.2 s内完成对未知参数$\tilde{a}, \tilde{b}, \tilde{c}$的识别.

4 结论

本文设计了一个新型三维分数阶混沌系统，分析了此系统的动力学特性，系统具有复杂的动态行为，能够出现双翼、三翼、四翼等不同类型多翼混沌吸引子共存.对系统进行了Multisim模拟电路设计与仿真，仿真结果与数值分析相符，进一步验证了系统的混沌行为.基于分数阶Lyapunov稳定理论及定理1，设计了一个自适应同步控制器，仿真结果表明了所设计控制器的有效性.