CN113642119B - 快速确定不连续系统Lyapunov指数谱的方法 - Google Patents

Abstract

本发明公开了一种快速确定不连续系统Lyapunov指数谱的方法，该方法以不连续系统状态参数初始值和沿相空间每个基本矢量的扰动量为初始条件，确定相轨迹；其次，采取差商近似导数法，获得Jacobi矩阵的近似矩阵。然后，对Jacobi矩阵进行特征值提取，得到不连续系统的Lyapunov指数谱；最后，将新算法应用到二自由度干摩擦冲击振荡器系统实例中，并将计算结果与同步方法的计算结果进行对比，对新方法的有效性进行了验证。本方法不仅适用于离散系统和连续时间系统，而且能够快速确定复杂不连续系统的Lyapunov指数，为确定复杂不连续系统的动力学行为提供了新思路。

Description

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

技术领域

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

背景技术

Lyapunov指数是在局部和全局稳定性准则下确定系统稳定性的数值，它表征了系统在相空间中相邻轨道间收敛或发散的平均指数率。从数学角度，Lyapunov指数是状态转移矩阵的特征值，决定了系统吸引子上轨迹无穷小扰动的演化进程。因此，Lyapunov指数被认为是衡量系统对初始条件敏感性的重要指标。所有Lyapunov指数值的和决定了系统在相空间中体积的发散度，所有Lyapunov指数值的和为负值是存在全局稳定吸引子的充要条件，并且表明系统是耗散的；若所有Lyapunov指数均为负值，则吸引子为稳定点；对于n维系统，如果前m个Lyapunov指数等于零而另外n-m个Lyapunov指数均为负值，则吸引子为m维环面；若最大Lyapunov指数为正值，则为混沌运动，若存在不止一个Lyapunov指数为正值，则为超混沌运动。对于保守系统，所有Lyapunov指数值的和为零，若所有Lyapunov指数值的和为正值，则表明系统是全局不稳定的，相轨迹呈指数发散。

Benettin和Shimada等人首先提出了一个基于Oseledets定理计算Lyapunov指数谱的有效算法，Benettin等人对算法进行了改进。近年来，提出了基于系统扰动的数量积和导数计算Lyapunov指数的方法，上述算法可以计算连续系统的Lyapunov指数谱，但是如果系统中存在不连续性则会出现较大误差或数值溢出现象。

为了确定具有冲击或干摩擦的机械系统、分段线性特性电振荡器等不连续系统的Lyapunov指数，可以采用等效的连续函数近似不连续分量、最小二乘阴影法和伪轨道法消除状态扰动矢量的不连续性等方法，虽然不需要重构相空间和Jacobi矩阵，但容易引起系统动力学特性的质变。Takens和Wolf等人根据时间序列信号重构系统吸引子，提出了计算不连续系统最大Lyapunov指数的算法，多年以来，Rosenstein等人针对该类方法的数据量、稳定性和计算速率等方面进行了评估改进，该类方法的优点是物理意义明显，便于理解，计算结果不易受到拓扑复杂性的影响，有一定的抗噪能力。近年来，部分学者又提出了映射法、补偿矩阵法、完全同步法、克隆动力学法等，这些方法对一些特性类型的不连续系统具有较高的准确性，但算法通用性不强且计算流程较为复杂。

发明内容

本发明的目的就是要提供一种快速确定不连续系统Lyapunov指数谱的方法，该方法以一般性的理论为依据，更具有通用性；另一方面，以经典的数值理论为基础，对广泛应用的Jacobi法进行改进，保证了准确性。

为实现此目的，本发明所设计的一种快速确定不连续系统Lyapunov指数谱的方法，包括如下步骤：

步骤1，设定n维不连续系统的初始条件，初始条件包括不连续系统初始运动状态的初始时间t₀、不连续系统在相空间中的初始位移矢量x₀和初始扰动量i＝1,…,n；j＝1,2,…,T，i表示不连续系统的维数，n为不连续系统的最大维数，j表示不连续系统的运行时刻，T表示不连续系统运行结束的时刻；

步骤2，利用初始位移矢量x₀和初始时间t₀，采用数值方法得到不连续系统吸引子的轨迹

步骤3，对不连续系统吸引子的轨迹利用差商近似倒数的方法得到不连续系统Jacobi矩阵的近似矩阵U_t(x)；

步骤4，将初始扰动量i＝1,…,n；j＝1,2,…,T和Jacobi矩阵的近似矩阵U_t(x)相乘得到新的初始扰动量/>并对新的初始扰动量/>利用Gram-Schmidt正交归一化获得相互正交的初始扰动量/>i＝1,…,n；j＝1,2,,T；

步骤5，对相互正交的初始扰动量i＝1,…,n；j＝1,2,,T采用Lyapunov指数公式估算法得到不连续系统的整个Lyapunov指数谱λ_i。

本发明的有益效果是提供了一种切实可行的快速确定不连续系统Lyapunov指数谱的方法，仅需要选定初始条件，无需确切知道系统的方程式，便能够确定系统的最终状态，而且对不连续系统和连续系统均适用。丰富了不连续系统Lyapunov指数谱的方法，为确定复杂不连续系统的动力学行为提供了新思路。

附图说明

图1为本发明方法流程图；

图2为本发明实施例的二自由度干摩擦冲击振荡器模型；

图3为本发明实施例的系统状态变量的分叉图和三个最大Lyapunov指数图；

图4为本发明实施例的系统状态变量的分叉图和三个最大Lyapunov指数图；

图5为本发明实施例的系统状态变量的分叉图和三个最大Lyapunov指数图；

图6为矩阵U_t(x)逼近算法示意图。

具体实施方式

以下结合附图和具体实施例对本发明作进一步的详细说明：

二自由度干摩擦冲击振荡器系统由连接到刚度为k₁的弹簧的质量块M和阻尼系数为c₁的粘滞阻尼器组成，如图2所示，另一个质量块m位于质量块M上，其相对于M的位置由缓冲器A和B限制，质量块之间存在摩擦力F_T和阻尼系数为c₂的粘性阻尼力。二自由度冲击振子是一个质量-弹簧-阻尼系统，其中两个质量块可以相互碰撞，并在相对运动过程中存在干摩擦。

二自由度冲击振荡器可由以下方程组描述：

式中，表示x₁的一阶导数，/>表示x₁的二阶导数，/>表示x₂的一阶导数，/>表示x₂的二阶导数，Ω是谐波激励的角频率；t是不连续系统的运行时间。

上述方程组定义了振荡器系统的演化；

上述方程组描述了在振荡器系统发生碰撞时刻t_c处不连续性的状态转换。

式中，M是质量块M的质量；m是质量块m的质量；x₁是质量块M的位移；是质量块M的速度；/>是质量块M的加速度；x₂是质量块m的位移；/>是质量块m的速度；/>是质量块m的加速度；k₁是弹簧刚度；c₁是粘滞阻尼器的阻尼系数；c₂是质量块M和质量块m之间的阻尼系数；F_T是摩擦力；A是强迫振幅；δ₀是内部质量块m与外部质量块M之间的间隙；μ是摩擦系数；sgn()是阶跃函数，当内部参数为正数时返回1，当内部参数为0时返回0，当内部参数为负数时返回-1；R是恢复系数；t_c是发生碰撞的时间；/>是碰撞时间的右极限时间；/>是碰撞时间的左极限时间。

将式(1)、(2)、(3)转换为无量纲表示得到：

其中, y₅＝Ωt；/> 为振荡器系统解的一阶导数。

在这个振荡器系统中，冲击和干摩擦的两种不连续性都存在。假设对于任何初始条件y(0)，使得|y₃-y₁|＜y₀，可以数值估计轨迹S_τ(y(0))。唯一的限制是两个物体在碰撞过程中的速度没有定义。此外，还假设当|y₃(0)-y₁(0)|＜y₀和y₂≠y₄时，轨迹S_τ(y(0))相对于初始条件是可微的。

如图3～5所示为二自由度干摩擦冲击振荡器系统状态变量的分叉图和三个最大Lyapunov指数图，图中，y₁表示振荡器系统的位移x₁的无量纲解，y₂表示振荡器系统的位移x₂的无量纲解，y₃表示振荡器系统的位移x₃的无量纲解，λ₁表示振荡器系统1维方向的Lyapunov指数值，λ₂表示振荡器系统2维方向的Lyapunov指数值，λ₃表示振荡器系统3维方向的Lyapunov指数值。

将本发明方法应用到二自由度干摩擦冲击振荡器系统实例中，并将计算结果与同步方法的计算结果进行对比，对新算法的有效性进行验证。

为了将数值仿真结果与不同算法获得的结果进行比较验证，控制参数η在[1.2:1.45]范围内变化，并设定以下参数值：γ＝0.693，σ＝0.50，μ＝0.02，β＝0.05，y₀＝0.80，p＝1.00，R＝0.60。

一种快速确定不连续系统Lyapunov指数谱的方法，它包括如下步骤，如图1所示：

步骤1，设定n维不连续系统的初始条件，初始条件包括不连续系统初始运动状态的初始时间t₀、不连续系统在相空间中的初始位移矢量x₀和初始扰动量i＝1,…,n；j＝1,2,…,T，i表示不连续系统的维数，n为不连续系统的最大维数，j表示不连续系统的运行时刻，T表示不连续系统运行结束的时刻。在本实施例中，设定初始时间t₀和系统在相空间中的初始状态x₀均为零，初始扰动量Δx₀为10^-3。

步骤2，利用初始位移矢量x₀和初始时间t₀，采用数值方法得到不连续系统吸引子的轨迹在本实施例中，利用自适应步长的龙格-库塔(Runge-Kutta)方法对二自由度冲击振荡器系统进行了数值求解，得到不连续系统吸引子的轨迹/>在开始计算Lyapunov指数和将轨迹点保存到分岔图之前，每组参数针对/>运行该振荡器系统，以减少瞬态影响。

步骤3，对不连续系统吸引子的轨迹利用差商近似倒数的方法得到不连续系统Jacobi矩阵的近似矩阵U_t(x)。在本实施例中，将上述求得的二自由度干摩擦冲击振荡器系统吸引子的轨迹/>带入以下方程式估算矩阵U_t(x)：

其中，Δ值取为10^-6，e₁,…,e_n是实数集Rⁿ中的标准基。

步骤4，将初始位移矢量i＝1,…,n；j＝1,2,…,T和Jacobi矩阵的近似矩阵U_t(x)相乘得到新的初始扰动量/>并对新的初始扰动量/>利用Gram-Schmidt正交归一化获得相互正交的初始扰动量/>i＝1,…,n；j＝1,2,…,T；

步骤5，对相互正交的初始扰动量i＝1,…,n；j＝1,2,…,T采用Lyapunov指数公式估算法得到不连续系统的整个Lyapunov指数谱λ_i。在本实施例中，将不同的正交扰动量/>用于求解不连续系统相对应的Lyapunov指数，对于不同的扰动量/>i＝1,…,n，j＝1,2,…,T，使用以下公式估算不连续系统的整个Lyapunov指数谱λ_i：

其中，t_T为不同初始扰动量对应的时间。

最后，基于Python软件平台，编写程序实现上述计算流程。

由于最小扰动Δy₄的收缩速度过快，导致在计算自然对数时出现了数值问题。因此，仅给出了该振荡器系统的三个最大的Lyapunov指数。为了观察Lyapunov指数值的稳定性，在估算程序的每次迭代之后，再利用估算不连续系统的整个Lyapunov指数谱λ_i的公式计算结果。在每100次迭代后，估计最后100个λ₁,λ₂,λ₃值的标准差。当所有标准偏差低于阈值10^-4时，终止计算，并取最后100次迭代中λ₁,λ₂,λ₃的平均值作为最终值。图3～图5给出了三个最大Lyapunov指数图和振荡器系统状态分岔图。可以观察到，Lyapunov指数值恰当地反映了振荡器系统在分岔图中的行为。因此，通过新算法在二自由度干摩擦冲击振荡器系统(4)、(5)、(6)中的应用，证明了该Lyapunov指数新算法的有效性。最后，通过对使用本文算法获得的Lyapunov指数计算结果与同步方法得到的计算结果进行比较，可以得到，使用两种算法获得的最大Lyapunov指数图高度一致。从而，两种方法也相互验证了算法的正确性。

上述技术方案中，所述步骤2中，首先使用j-1时刻的位移矢量x^(j-1)在不连续系统中获得j时刻的位移矢量x^(j)，再采用数值方法求解获得不连续系统的轨迹所述数值方法为欧拉方法、自适应步长的龙格-库塔方法、有限差分法或雅可比法。

上述技术方案中，所述步骤3中求解不连续系统Jacobi矩阵的近似矩阵U_t(x)的具体方法如下：

对于不连续系统方程中的矢量场f在相空间轨迹的初始设定状态x处并不是连续可微分的，不能直接从变分方程获得Jacobi矩阵U_t(x_j)，不能够利用经典算法直接得出不连续系统的Lyapunov指数谱。首先考虑连续时间不连续系统，不连续系统的解S_t(x)实际上是将初始设定状态x转换为时间t之后不连续系统状态的映射。U_t(x)为在初始设定状态x处求得轨迹S_t(x)的Jacobi矩阵的近似矩阵：

其中，S_t(x+Δe₁)为在不连续系统在x+Δe₁状态下的解，S_t(x+Δe_n)为不连续系统在x+Δe_n状态下的解，表示不同扰动状态下不连续系统的解，e₁,…,e_n是实数集Rⁿ中的标准基，U_t(x)的列是S_t(x)关于后续状态变量的偏导数，根据导数定义，所述偏导数可以表示为轨迹扰动量Δ趋近于0时差商的极限，则可通过以下方式估算矩阵U_t(x)：

其中，Δ为轨迹扰动量，在选择Δ值时，需要在与Δ成比例的截断误差和取决于实现中使用数据类型的近似误差之间取得平衡，矩阵U_t(x)的近似算法如图6所示，x+Δe₁和x+Δe₂分别为初始设定状态x受Δe₁和Δe₂扰动的初始设定状态，S_t(x+Δe₁)和S_t(x+Δe₂)分别为不连续系统在x+Δe₁和x+Δe₂初始设定状态下的解。

利用未扰动的初始设定状态x以及沿着相空间的每个受Δ值扰动的初始设定状态x+Δe_i，i＝1,…,n，确定映射S_t(x)的值，进而估计Jacobi矩阵的近似矩阵，其中，Δe_i为正交化的扰动量。

上述技术方案中，所述步骤4中利用Gram-Schmidt正交归一化获得相互正交的初始扰动量的方法如下：

其中，表示1维方向j时刻的扰动量，/>表示1维方向j时刻的正交化扰动量，新的初始扰动量/>是由相互正交的初始扰动量/>引起轨迹/>的扰动，即

上述技术方案中，所述步骤5的具体实现方法为：

将相互正交的扰动量用于求解不连续系统相对应的Lyapunov指数，使用以下公式估算不连续系统的整个Lyapunov指数谱λ_i：

其中，T表示不连续系统运行结束的时刻，t_T为不同初始扰动量对应的时间。

从上述实施例可以看出，本发明可以有效的解决不连续系统Lyapunov指数的计算问题，并且该方法具有通用性，其应用是成功的。

上述实施例只为说明本发明的技术构思及特点，其目的在于让熟悉此项技术的人士能够了解本发明的内容并据以实施，并不能以此限制本发明的保护范围。凡根据本发明精神实质所作的等效变化或修饰，都应涵盖在本发明的保护范围之内。

本说明书未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (6)

1.一种快速确定不连续系统Lyapunov指数谱的方法，其特征在于，它包括如下步骤：

步骤1，设定n维不连续系统的初始条件，初始条件包括不连续系统初始运动状态的初始时间t₀、不连续系统在相空间中的初始位移矢量x₀和初始扰动量i表示不连续系统的维数，n为不连续系统的最大维数，j表示不连续系统的运行时刻，T表示不连续系统运行结束的时刻；所述不连续系统为二自由度干摩擦冲击振荡器系统；

步骤2，利用初始位移矢量x₀和初始时间t₀，采用数值方法得到不连续系统吸引子的轨迹

步骤3，对不连续系统吸引子的轨迹利用差商近似倒数的方法得到不连续系统Jacobi矩阵的近似矩阵U_t(x)；

步骤4，将初始扰动量和Jacobi矩阵的近似矩阵U_t(x)相乘得到新的初始扰动量/>并对新的初始扰动量/>利用Gram-Schmidt正交归一化获得相互正交的初始扰动量/>

步骤5，对相互正交的初始扰动量采用Lyapunov指数公式估算法得到不连续系统的整个Lyapunov指数谱λ_i。

2.根据权利要求1所述的快速确定不连续系统Lyapunov指数谱的方法，其特征在于：所述步骤2中，首先使用j-1时刻的位移矢量x^(j-1)在不连续系统中获得j时刻的位移矢量x^(j)，再采用数值方法求解获得不连续系统的轨迹

3.根据权利要求2所述的快速确定不连续系统Lyapunov指数谱的方法，其特征在于：

所述数值方法为欧拉方法、自适应步长的龙格-库塔方法、有限差分法或雅可比法。

4.根据权利要求1所述的快速确定不连续系统Lyapunov指数谱的方法，其特征在于：所述步骤3中求解不连续系统Jacobi矩阵的近似矩阵U_t(x)的具体方法如下：

首先考虑连续时间不连续系统，所述不连续系统的解S_t(x)是将初始设定状态x转换为时间t之后不连续系统状态的映射，U_t(x)为在初始设定状态x处求得轨迹S_t(x)的Jacobi矩阵的近似矩阵：

其中，S_t(x+Δe₁)为在不连续系统在x+Δe₁状态下的解，S_t(x+Δe_n)为不连续系统在x+Δe_n状态下的解，表述不同扰动状态下不连续系统的解，e₁,…,e_n是实数集Rⁿ中的标准基，U_t(x)的列是S_t(x)关于后续状态变量的偏导数，根据导数定义，所述偏导数可以表示为轨迹扰动量Δ趋近于0时差商的极限，则可通过以下方式估算矩阵U_t(x)：

其中，Δ为轨迹扰动量，在选择Δ值时，需要在与Δ成比例的截断误差和取决于实现中使用数据类型的近似误差之间取得平衡；

利用未扰动的初始设定状态x以及沿着相空间的每个受Δ值扰动的某一状态x+Δe_i，i＝1,…,n，确定映射S_t(x)的值，进而估计Jacobi矩阵的近似矩阵，其中，Δe_i为正交化的扰动量。

5.根据权利要求1所述的快速确定不连续系统Lyapunov指数谱的方法，其特征在于：所述步骤4中利用Gram-Schmidt正交归一化获得相互正交的初始扰动量的方法如下：

其中，表示1维方向j时刻的扰动量，/>表示1维方向j时刻的正交化扰动量，新的初始扰动量/>是由相互正交的初始扰动量/>引起轨迹S_t(x^(j-1))的扰动，即

6.根据权利要求1所述的快速确定不连续系统Lyapunov指数谱的方法，其特征在于：所述步骤5的具体实现方法为：

将相互正交的扰动量用于求解不连续系统相对应的Lyapunov指数，使用以下公式估算不连续系统的整个Lyapunov指数谱λ_i：

其中，t_T为不同初始扰动量对应的时间。