專利名稱:控制系統的李雅普諾夫指數譜的計算方法
技術領域：
本發明屬于非線性動力學系統理論應用技術，涉別涉及由測量所得的實驗數據序列計算系統的李雅普諾夫指數譜，從而確定系統的混沌特性的方法。
背景技術：
由故障誘發而成的故障混沌控制系統廣泛地存在于航空、航天、化工、發電等領域的控制系統中，它的存在嚴重影響控制系統的正常工作。如果對這種現象有深入研究，知道產生的內在規律，并能進一步尋求將處于無序運動狀態的故障控制系統改變成在平衡狀態的較大臨域內穩定的系統，或者使平衡狀態成為相空間中對初始條件不敏感的吸引子，將大大提高控制系統的可靠性。
李雅普諾夫(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