快速确定不连续系统Lyapunov指数谱的方法与流程

快速确定不连续系统lyapunov指数谱的方法
技术领域
1.本发明涉及混沌系统识别技术领域，具体地指一种快速确定不连续系统lyapunov指数谱的方法。


背景技术：

2.lyapunov指数是在局部和全局稳定性准则下确定系统稳定性的数值，它表征了系统在相空间中相邻轨道间收敛或发散的平均指数率。从数学角度，lyapunov指数是状态转移矩阵的特征值，决定了系统吸引子上轨迹无穷小扰动的演化进程。因此，lyapunov指数被认为是衡量系统对初始条件敏感性的重要指标。所有lyapunov指数值的和决定了系统在相空间中体积的发散度，所有lyapunov指数值的和为负值是存在全局稳定吸引子的充要条件，并且表明系统是耗散的；若所有lyapunov指数均为负值，则吸引子为稳定点；对于n维系统，如果前m个lyapunov指数等于零而另外n
‑
m个lyapunov指数均为负值，则吸引子为m维环面；若最大lyapunov指数为正值，则为混沌运动，若存在不止一个lyapunov指数为正值，则为超混沌运动。对于保守系统，所有lyapunov指数值的和为零，若所有lyapunov指数值的和为正值，则表明系统是全局不稳定的，相轨迹呈指数发散。
3.benettin和shimada等人首先提出了一个基于oseledets定理计算lyapunov指数谱的有效算法，benettin等人对算法进行了改进。近年来，提出了基于系统扰动的数量积和导数计算lyapunov指数的方法，上述算法可以计算连续系统的lyapunov指数谱，但是如果系统中存在不连续性则会出现较大误差或数值溢出现象。
4.为了确定具有冲击或干摩擦的机械系统、分段线性特性电振荡器等不连续系统的lyapunov指数，可以采用等效的连续函数近似不连续分量、最小二乘阴影法和伪轨道法消除状态扰动矢量的不连续性等方法，虽然不需要重构相空间和jacobi矩阵，但容易引起系统动力学特性的质变。takens和wolf等人根据时间序列信号重构系统吸引子，提出了计算不连续系统最大lyapunov指数的算法，多年以来，rosenstein等人针对该类方法的数据量、稳定性和计算速率等方面进行了评估改进，该类方法的优点是物理意义明显，便于理解，计算结果不易受到拓扑复杂性的影响，有一定的抗噪能力。近年来，部分学者又提出了映射法、补偿矩阵法、完全同步法、克隆动力学法等，这些方法对一些特性类型的不连续系统具有较高的准确性，但算法通用性不强且计算流程较为复杂。


技术实现要素：

5.本发明的目的就是要提供一种快速确定不连续系统lyapunov指数谱的方法，该方法以一般性的理论为依据，更具有通用性；另一方面，以经典的数值理论为基础，对广泛应用的jacobi法进行改进，保证了准确性。
6.为实现此目的，本发明所设计的一种快速确定不连续系统lyapunov指数谱的方法，包括如下步骤：
7.步骤1，设定n维不连续系统的初始条件，初始条件包括不连续系统初始运动状态
的二阶导数，ω是谐波激励的角频率；t是不连续系统的运行时间。
24.上述方程组定义了振荡器系统的演化；
[0025][0026][0027]
上述方程组描述了在振荡器系统发生碰撞时刻t
c
处不连续性的状态转换。
[0028]
式中，m是质量块m的质量；m是质量块m的质量；x1是质量块m的位移；是质量块m的速度；是质量块m的加速度；x2是质量块m的位移；是质量块m的速度；是质量块m的加速度；k1是弹簧刚度；c1是粘滞阻尼器的阻尼系数；c2是质量块m和质量块m之间的阻尼系数；f
t
是摩擦力；a是强迫振幅；δ0是内部质量块m与外部质量块m之间的间隙；μ是摩擦系数；sgn()是阶跃函数，当内部参数为正数时返回1，当内部参数为0时返回0，当内部参数为负数时返回
‑
1；r是恢复系数；t
c
是发生碰撞的时间；是碰撞时间的右极限时间；是碰撞时间的左极限时间。
[0029]
将式(1)、(2)、(3)转换为无量纲表示得到：
[0030][0031][0032][0033]
其中,其中,y5＝ωt；＝ωt；为振荡器
系统解的一阶导数。
[0034]
在这个振荡器系统中，冲击和干摩擦的两种不连续性都存在。假设对于任何初始条件y(0)，使得|y3‑
y1|＜y0，可以数值估计轨迹s
τ
(y(0))。唯一的限制是两个物体在碰撞过程中的速度没有定义。此外，还假设当|y3(0)
‑
y1(0)|＜y0和y2≠y4时，轨迹s
τ
(y(0))相对于初始条件是可微的。
[0035]
如图3～5所示为二自由度干摩擦冲击振荡器系统状态变量的分叉图和三个最大lyapunov指数图，图中，y1表示振荡器系统的位移x1的无量纲解，y2表示振荡器系统的位移x2的无量纲解，y3表示振荡器系统的位移x3的无量纲解，λ1表示振荡器系统1维方向的lyapunov指数值，λ2表示振荡器系统2维方向的lyapunov指数值，λ3表示振荡器系统3维方向的lyapunov指数值。
[0036]
将本发明方法应用到二自由度干摩擦冲击振荡器系统实例中，并将计算结果与同步方法的计算结果进行对比，对新算法的有效性进行验证。
[0037]
为了将数值仿真结果与不同算法获得的结果进行比较验证，控制参数η在[1.2:1.45]范围内变化，并设定以下参数值：γ＝0.693，σ＝0.50，μ＝0.02，β＝0.05，y0＝0.80，p＝1.00，r＝0.60。
[0038]
一种快速确定不连续系统lyapunov指数谱的方法，它包括如下步骤，如图1所示：
[0039]
步骤1，设定n维不连续系统的初始条件，初始条件包括不连续系统初始运动状态的初始时间t0、不连续系统在相空间中的初始位移矢量x0和初始扰动量i＝1,
…
,n；j＝1,2,
…
,t，i表示不连续系统的维数，n为不连续系统的最大维数，j表示不连续系统的运行时刻，t表示不连续系统运行结束的时刻。在本实施例中，设定初始时间t0和系统在相空间中的初始状态x0均为零，初始扰动量δx0为10
‑3。
[0040]
步骤2，利用初始位移矢量x0和初始时间t0，采用数值方法得到不连续系统吸引子的轨迹在本实施例中，利用自适应步长的龙格
‑
库塔(runge
‑
kutta)方法对二自由度冲击振荡器系统进行了数值求解，得到不连续系统吸引子的轨迹在开始计算lyapunov指数和将轨迹点保存到分岔图之前，每组参数针对运行该振荡器系统，以减少瞬态影响。
[0041]
步骤3，对不连续系统吸引子的轨迹利用差商近似倒数的方法得到不连续系统jacobi矩阵的近似矩阵u
t
(x)。在本实施例中，将上述求得的二自由度干摩擦冲击振荡器系统吸引子的轨迹带入以下方程式估算矩阵u
t
(x)：
[0042][0043]
其中，δ值取为10
‑6，e1,
…
,e
n
是实数集r
n
中的标准基。
[0044]
步骤4，将初始位移矢量i＝1,
…
,n；j＝1,2,
…
,t和jacobi矩阵的近似矩阵u
t
(x)相乘得到新的初始扰动量并对新的初始扰动量利用gram
‑
schmidt正交归
一化获得相互正交的初始扰动量i＝1,
…
,n；j＝1,2,
…
,t；
[0045]
步骤5，对相互正交的初始扰动量i＝1,
…
,n；j＝1,2,
…
,t采用lyapunov指数公式估算法得到不连续系统的整个lyapunov指数谱λ
i
。在本实施例中，将不同的正交扰动量用于求解不连续系统相对应的lyapunov指数，对于不同的扰动量i＝1,
…
,n，j＝1,2,
…
,t，使用以下公式估算不连续系统的整个lyapunov指数谱λ
i
：
[0046][0047]
其中，t
t
为不同初始扰动量对应的时间。
[0048]
最后，基于python软件平台，编写程序实现上述计算流程。
[0049]
由于最小扰动δy4的收缩速度过快，导致在计算自然对数时出现了数值问题。因此，仅给出了该振荡器系统的三个最大的lyapunov指数。为了观察lyapunov指数值的稳定性，在估算程序的每次迭代之后，再利用估算不连续系统的整个lyapunov指数谱λ
i
的公式计算结果。在每100次迭代后，估计最后100个λ1,λ2,λ3值的标准差。当所有标准偏差低于阈值10
‑4时，终止计算，并取最后100次迭代中λ1,λ2,λ3的平均值作为最终值。图3～图5给出了三个最大lyapunov指数图和振荡器系统状态分岔图。可以观察到，lyapunov指数值恰当地反映了振荡器系统在分岔图中的行为。因此，通过新算法在二自由度干摩擦冲击振荡器系统(4)、(5)、(6)中的应用，证明了该lyapunov指数新算法的有效性。最后，通过对使用本文算法获得的lyapunov指数计算结果与同步方法得到的计算结果进行比较，可以得到，使用两种算法获得的最大lyapunov指数图高度一致。从而，两种方法也相互验证了算法的正确性。
[0050]
上述技术方案中，所述步骤2中，首先使用j
‑
1时刻的位移矢量x
(j
‑
1)
在不连续系统中获得j时刻的位移矢量x
(j)
，再采用数值方法求解获得不连续系统的轨迹所述数值方法为欧拉方法、自适应步长的龙格
‑
库塔方法、有限差分法或雅可比法。
[0051]
上述技术方案中，所述步骤3中求解不连续系统jacobi矩阵的近似矩阵u
t
(x)的具体方法如下：
[0052]
对于不连续系统方程中的矢量场f在相空间轨迹的初始设定状态x处并不是连续可微分的，不能直接从变分方程获得jacobi矩阵u
t
(x
j
)，不能够利用经典算法直接得出不连续系统的lyapunov指数谱。首先考虑连续时间不连续系统，不连续系统的解s
t
(x)实际上是将初始设定状态x转换为时间t之后不连续系统状态的映射。u
t
(x)为在初始设定状态x处求得轨迹s
t
(x)的jacobi矩阵的近似矩阵：
[0053][0054]
其中，s
t
(x δe1)为在不连续系统在x δe1状态下的解，s
t
(x δe
n
)为不连续系统在x δe
n
状态下的解，表示不同扰动状态下不连续系统的解，e1,
…
,e
n
是实数集r
n
中的标准基，u
t
(x)的列是s
t
(x)关于后续状态变量的偏导数，根据导数定义，所述偏导数可以表示为轨迹扰动量δ趋近于0时差商的极限，则可通过以下方式估算矩阵u
t
(x)：
[0055][0056]
其中，δ为轨迹扰动量，在选择δ值时，需要在与δ成比例的截断误差和取决于实现中使用数据类型的近似误差之间取得平衡，矩阵u
t
(x)的近似算法如图6所示，x δe1和x δe2分别为初始设定状态x受δe1和δe2扰动的初始设定状态，s
t
(x δe1)和s
t
(x δe2)分别为不连续系统在x δe1和x δe2初始设定状态下的解。
[0057]
利用未扰动的初始设定状态x以及沿着相空间的每个受δ值扰动的初始设定状态x δe
i
，i＝1,
…
,n，确定映射s
t
(x)的值，进而估计jacobi矩阵的近似矩阵，其中，δe
i
为正交化的扰动量。
[0058]
上述技术方案中，所述步骤4中利用gram
‑
schmidt正交归一化获得相互正交的初始扰动量的方法如下：
[0059][0060]
其中，表示1维方向j时刻的扰动量，表示1维方向j时刻的正交化扰动量，新的初始扰动量是由相互正交的初始扰动量引起轨迹的扰动，即
[0061]
上述技术方案中，所述步骤5的具体实现方法为：
[0062]
将相互正交的扰动量用于求解不连续系统相对应的lyapunov指数，使用以下公式估算不连续系统的整个lyapunov指数谱λ
i
：
[0063][0064]
其中，t表示不连续系统运行结束的时刻，t
t
为不同初始扰动量对应的时间。
[0065]
从上述实施例可以看出，本发明可以有效的解决不连续系统lyapunov指数的计算问题，并且该方法具有通用性，其应用是成功的。
[0066]
上述实施例只为说明本发明的技术构思及特点，其目的在于让熟悉此项技术的人士能够了解本发明的内容并据以实施，并不能以此限制本发明的保护范围。凡根据本发明精神实质所作的等效变化或修饰，都应涵盖在本发明的保护范围之内。
[0067]
本说明书未作详细描述的内容属于本领域专业技术人员公知的现有技术。