专利名称:控制系统的李雅普诺夫指数谱的计算方法
技术领域：
本发明属于非线性动力学系统理论应用技术，涉别涉及由测量所得的实验数据序列计算系统的李雅普诺夫指数谱，从而确定系统的混沌特性的方法。
背景技术：
由故障诱发而成的故障混沌控制系统广泛地存在于航空、航天、化工、发电等领域的控制系统中，它的存在严重影响控制系统的正常工作。如果对这种现象有深入研究，知道产生的内在规律，并能进一步寻求将处于无序运动状态的故障控制系统改变成在平衡状态的较大临域内稳定的系统，或者使平衡状态成为相空间中对初始条件不敏感的吸引子，将大大提高控制系统的可靠性。
李雅普诺夫(Lyapunov)指数定量的描述了系统相空间中相邻轨道呈指数发散或收敛的性质。若此系统的李雅普诺夫指数大于零，表示相邻轨道按指数规律发散，长时间行为对初始值敏感，我们称此系统是混沌动力学系统。若控制系统存在正的李雅普诺夫指数，则此系统是非正常工作，不可控的。若常见的控制系统产生非线性故障就有可能使系统输出呈现混沌运动。
对高维系统测量得到的一般是一维的数据序列，由此序列计算系统的李雅普诺夫指数。现在常用的计算李雅普诺夫指数的方法有BBA方法和Wolf方法等，BBA方法见(ReggieBrown，Paul Bryant and Henry D.I Abarbanel，Computing the Lyapunov spectrum of adynamical system from an observed time series，Phys.Rev.A，Vol.43，No.6，pp.2787-2805，1991)；Wolf方法见(WolfA，et al.Determining Lyapunov exponents from a time series[J].PhysD.1985，16285-317)。其中Wolf方法仅适用于计算系统的最大李雅普诺夫指数，BBA方法可求出系统的全部李雅普诺夫指数，但运算量大，需要的数据点很多，其应用受到很大限制。
发明内容本发明要解决的技术问题是提供一种控制系统李雅普诺夫指数谱的计算方法，利用该方法，可以逼近任意非线性函数的能力进行李雅普诺夫指数指数谱的计算，不需要很多的数据点就可以得到系统的全部李雅普诺夫指数，且运算量较小。
本发明解决其技术问题所采用的技术方案是控制系统的李雅普诺夫指数谱的计算方法，包括以下步骤(1)利用相空间重构理论重现系统；(2)利用RBF神经网络的逼近能力计算所述系统的Oseledec矩阵；(3)计算所述Oseledec矩阵的特征值，求得控制系统的李雅普诺夫指数谱。
所述利用相空间重构理论重现系统的方法为对多维系统测量得到数据序列x(n)(n＝1，2，ΛN)，其中，x(n)为t0+nτ时刻记录的数据，τ为采样时间，t0为采样起始时间；根据相空间重构理论，设y(t)＝[x(t)，x(t+T)，Λx(t+(d-1)T)]，其中d为嵌入维数，T为时滞，且T是τ的整数倍，所述y(t)就是d维重构相空间中的一点，它随时间的变化形成d维欧氏空间的新的动力学系统，即y(t+T)＝F(y(t))。
计算所述系统的Oseledec矩阵的方法为所述公式y(t+T)＝F(y(t))出现微小扰动时，得到Δy(t+T)＝J(y(t))*Δy(t)，其中J(y(t))是映射F的雅可比矩阵；设映射F的第i个分量是fi，y(t)的第j个分量是xj，则J(y(t))=∂f1∂x1∂f1∂x2Λ∂f1∂xd∂f2∂x1∂f2∂x2Λ∂f2∂x2MMOM∂fd∂x1∂fd∂x2Λ∂fd∂xd,]]>根据y(t)的得到x(t+d×T)＝fd(y(t))，J(y(t))的形式为J(y(t))=010Λ0001Λ0MMMO1∂fd∂x1∂fd∂x2∂fd∂x3Λ∂fd∂xd,]]>按照该规律计算N次JN＝J(t+NT)J(t+(N-1)T)ΛJ(t)，根据Oseledec乘积遍历性定理，构造Oseledec矩阵
limN→∞[{JNT}{JN}]1/2N.]]>所述步骤(3)中计算所述Oseledec矩阵的特征值时包括如下步骤先设求解的长乘积矩阵T＝TN·TN-1ΛT1，对此长乘积矩阵，计算Ti·Q(i-1)的QR分解T1·Q(0)＝Q(1)·R(1)，T2·Q(1)＝Q(2)·R(2)，Λ，TN·Q(N-1)＝Q(N)·R(N)其中Q(i)为正交矩阵，R(i)为上三角矩阵，Q(0)为d×d阶单位矩阵，按公式limN→∞[{JNT}{JN}]1/2N]]>分解后，得到矩阵T的特征值λk=1NΣj=1N1nR(j)kk]]>，k＝1，2，Λd。
本发明的控制系统的李雅普诺夫指数谱的计算方法在计算高维动力学系统的李雅普诺夫指数时，运算量较小，需要的计算时间很少，且使用较少的数据点能很快的计算出系统的指数谱，解决了现有技术中Wolf和BBA方法计算速度慢、所需样本点多等缺点，其结果更准确，可以得到系统的全部李雅普诺夫指数谱。
具体实施方式下面结合具体实施例对本发明进一步说明。
(1)利用相空间重构理论重现系统假设对多维系统测量得到数据序列x(n)(n＝1，2，AN)，其中，x(n)表示t0+nτ时刻记录的数据，τ为采样时间，t0为采样起始时间。根据相空间重构理论，设y(t)＝[x(t)，x(t+T)，Λx(t+(d-1)T]，其中d为嵌入维数；T为时滞，是τ的整数倍。y(t)就是d维重构相空间中的一点，它随时间的变化就形成了d维欧氏空间的一个新的动力学系统，即y(t+T)＝F(y(t)) (1)通过分析(1)式的动力学系统就可以了解原系统x(n)的动力学特性。
(2)Oseledec矩阵的确定若(1)式出现微小扰动，则可以得到Δy(t+T)＝J(y(t))*Δy(t) (2)其中J(y(t))是映射F的雅可比矩阵。设映射F的第i个分量是fi，y(t)的第j个分量是xj，则J(y(t))=∂f1∂x1∂f1∂x2Λ∂f1∂xd∂f2∂x1∂f2∂x2Λ∂f2∂xdMMOM∂fd∂x1∂fd∂x2Λ∂fd∂xd---(3)]]>根据上述y(t)的定义可知x(t+d×T)＝fd(y(t)) (4)J(y(t))的形式为J(y(t))=010Λ0001Λ0MMMO1∂fd∂x1∂fd∂x2∂fd∂x3Λ∂fd∂xd---(5)]]>按照此种规律计算N次有JN＝J(t+NT)J(t+(N-1)T)ΛJ(t) (6)根据Oseledec乘积遍历性定理，可构造Oseledec矩阵limN→∞[{JNT}{JN}]1/2N---(7)]]>计算此矩阵的特征值，即可求出原数据序列的李雅普诺夫指数谱。
在计算过程中，由于(7)式定义的矩阵存在着指数和分数幂，此矩阵往往是病态的，难以直接精确计算它的全部特征值。采用长乘积矩阵分解技术，可以解决此问题。先设求解的矩阵T＝TN·TN-1ΛT1(8)对此长乘积矩阵，计算Ti·Q(i-1)的QR分解T1·Q(0)＝Q(1)·R(1)，T2·Q(1)＝Q(2)·R(2)，(9)Λ，TN·Q(N-1)＝Q(N)·R(N)式中Q(i)为正交矩阵，R(i)为上三角矩阵，Q(0)是d×d阶单位矩阵，按(7)式分解后，得到矩阵T的特征值λk=1NΣj=1N1nR(j)kk]]>，k＝1，2，Λd (10)因此，要求解数据序列所反映系统的李雅普诺夫指数，就需求出(6)式中的每个雅可比矩阵。本发明采用RBF神经网络逼近(4)式，从而得到雅可比矩阵的最后一行矩阵元。
(3)RBF神经网络RBF网络属于局部逼近式多层前向神经网络。此网络结构简单、训练简洁而且学习收敛速度快，能够逼近任意非线性函数。它采用径向基函数作为隐含层结点的基函数，隐含层对输入向量进行非线性变换，再经输出层的线性叠加输出。
本发明采用广义RBF神经网络结构，其学习算法采用自组织选取中心法。这种方法由两个阶段构成一是自组织学习阶段，即学习隐层基函数中心与方差的阶段；二是学习输出层权值的阶段。其中自组织学习阶段采用K-均值聚类算法，网络输出层权值可通过求解线性方程组来确定，即直接用伪逆的方法求解。
下面对本发明的算法进行验证以Henon映射为例，按上述方法计算它的李雅普诺夫指数谱。Henon映射表示为x(t+1)＝1-ax(t)2+y(t) (11)y(t+1)＝bx(t)其中当a＝1.4，b＝0.3时，系统成为一个典型的二维混沌动力学系统，原系统的李雅普诺夫指数分别为λ1＝0.408，λ2＝-1.62。根据“Reggie Brown，Paul Bryant，HenryD.I.Abarbanel，Computing the Lyapunov spectrum of a dynamical system from an observed timeseries.Physical ReviewA.Vol.43.No.6.pp.2787-2805.1991”的计算结果，当嵌入维d＝2时，已经可以用重构的相空间计算出原系统的李雅普诺夫指数，这说明此时的重构相空间已经可以反映原系统的动力学特性。因此可采用嵌入维d＝2，时滞T＝1的重构相空间验证本发明算法的正确性。
以初值(0.25，0.25)迭代产生的x(t)混沌序列为例进行计算，去掉前面200个从初始位置到吸引子的过渡点，每100个点拟合一次，计算系统的李雅普诺夫指数谱。
本发明采用广义RBF神经网络，隐层结点数为51个，其中有一个结点的输出恒为1，其他结点采用Gauss函数作为基函数。K-均值聚类算法前后两次中心的变化小于ε＜0.1时结束调整。
下表(表1)是在d＝2时不同样本点的比较，
在表1中，样本点数目增加1.5倍，李雅普诺夫指数的精度变为1.7％，没有很显著得提高。这说明较少的样本点数目就能算出较准确的李雅普诺夫指数，解决了BBA算法中运算量大，所需数据点多的缺点。
下表(表2)是样本点为200个时不同嵌入维数的比较，
该表的结果表明，本发明的方法在计算高维动力学系统的李雅普诺夫指数时，运算量较小，需要的计算时间很少，在一般的PC机上就可以实现。
本发明采用的算法，使用较少的数据点能很快的计算出系统的指数谱，解决了以往Wolf和BBA方法计算速度慢、所需样本点多的缺点，得到了较为满意的结果。
权利要求
1.控制系统的李雅普诺夫指数谱的计算方法，包括以下步骤(1)利用相空间重构理论重现系统；(2)利用RBF神经网络的逼近能力计算所述系统的Oseledec矩阵；(3)计算所述Oseledec矩阵的特征值，求得控制系统的李雅普诺夫指数谱。
2.根据权利要求
1所述的控制系统的李雅普诺夫指数谱的计算方法，其特征在于所述利用相空间重构理论重现系统的方法为对多维系统测量得到数据序列x(n)(n＝1，2，ΛN)，其中，x(n)为t0+nτ时刻记录的数据，τ为采样时间，t0为采样起始时间；根据相空间重构理论，设y(t)＝[x(t)，x(t+T)，Λx(t+(d-1)T)]，其中d为嵌入维数，T为时滞，且T是τ的整数倍，所述y(t)就是d维重构相空间中的一点，它随时间的变化形成d维欧氏空间的新的动力学系统，即y(t+T)＝F(y(t))。
3.根据权利要求
2所述的控制系统的李雅普诺夫指数谱的计算方法，其特征在于计算所述系统的Oseledec矩阵的方法为所述公式y(t+T)＝F(y(t))出现微小扰动时，得到Δy(t+T)＝J(y(t))*Δy(t)，其中J(y(t))是映射F的雅可比矩阵；设映射F的第i个分量是fi，y(t)的第j个分量是xj，则J(y(t))=∂f1∂x1∂f1∂x2Λ∂f1∂xd∂f2∂x1∂f2∂x2Λ∂f2∂xdMMOM∂fd∂x1∂fd∂x2Λ∂fd∂xd,]]>根据y(t)的得到x(t+d×T)＝fd(y(t))，J(y(t))的形式为J(y(t))=010Λ0001Λ0MMMO1∂fd∂x1∂fd∂x2∂fd∂x3Λ∂fd∂xd,]]>按照该规律计算N次JN＝J(t+NT)J(t+(N-1)T)ΛJ(t)，根据Oseledee乘积遍历性定理，构造Oseledec矩阵limN→∞[{JNT}{JN}]1/2N.]]>
4.根据权利要求
3所述的控制系统的李雅普诺夫指数谱的计算方法，其特征在于所述步骤(3)中计算所述Oseledec矩阵的特征值时包括如下步骤先设求解的长乘积矩阵T＝TN·TN-1ΛT1，对此长乘积矩阵，计算Ti·Q(i-1)的QR分解T1·Q(0)＝Q(1)·R(1)，T2·Q(1)＝Q(2)·R(2)，Λ，TN·Q(N-1)＝Q(N)·R(N)其中Q(i)为正交矩阵，R(i)为上三角矩阵，Q(0)为d×d阶单位矩阵，按公式limN→∞[{JNT}{JN}]1/2N]]>分解后，得到矩阵T的特征值λk=1NΣj=1NlnR(j)kk,k=1,2,Λd.]]>
专利摘要
本发明公开了一种控制系统李雅普诺夫指数谱的计算方法，属于非线性动力学系统理论应用技术。本发明的计算方法包括(1)利用相空间重构理论重现系统；(2)利用RBF神经网络的逼近能力计算所述系统的Oseledec矩阵；(3)计算所述Oseledec矩阵的特征值，求得控制系统的李雅普诺夫指数谱。本发明的控制系统的李雅普诺夫指数谱的计算方法在计算高维动力学系统的李雅普诺夫指数时，运算量较小，需要的计算时间很少，且使用较少的数据点能很快的计算出系统的指数谱，解决了现有技术中Wolf和BBA方法计算速度慢、所需样本点多等缺点，其结果更准确，可以得到系统的全部李雅普诺夫指数谱。
文档编号G05B13/02GK1996176SQ200610169885
公开日2007年7月11日 申请日期2006年12月30日
发明者吴云洁, 王卫红, 刘正华, 赵媛媛 申请人:北京航空航天大学导出引文BiBTeX, EndNote, RefMan