線性時變系統的狀態空間模型遞推辨識研究

倪智宇， 吳志剛,2

(1. 大連理工大學 工業裝備結構分析國家重點實驗室，遼寧　大連　116024；2.大連理工大學 航空航天學院，遼寧　大連　116024)

?

線性時變系統的狀態空間模型遞推辨識研究

倪智宇1， 吳志剛1,2

(1. 大連理工大學 工業裝備結構分析國家重點實驗室，遼寧大連116024；2.大連理工大學 航空航天學院，遼寧大連116024)

摘要：針對線性時變系統中狀態空間模型的辨識問題，提出了一種新的模型參數矩陣的遞推辨識格式。不同于常用的利用奇異值分解(SVD)或者最小二乘原理計算時變狀態空間模型參數的方法，這種新的遞推方法基于信號子空間投影原理，通過重新建立輸入輸出數據之間的關系，構建新的信號子空間矩陣，從而遞推得到系統的時變狀態空間模型參數。與現有的計算時變狀態空間模型的方法相比，這種新的遞推方法由于不需要進行SVD的計算，從而大幅的減少了計算時間。特別是當系統的階次較高時，計算效率優勢更為明顯。在算例中將這種方法與經典的使用SVD的時變ERA(TV-ERA)方法從辨識結果和計算效率上進行了比較。仿真結果表明這種新的遞推算法能有效辨識狀態空間方程形式的線性時變系統的模型參數，和TV-ERA方法相比具有更高的計算效率。

關鍵詞：線性時變系統；遞推子空間方法；狀態空間模型；參數辨識

在現代工程應用中，越來越多的結構本身具有時變的特性，如大型空間結構的搬運裝配、太陽能帆板的展開與旋轉、機械臂的運動以及高速列車的行駛等，這些過程中往往都需要考慮到系統的時變特性。而在系統模型的描述中，狀態空間模型在系統極點配置、觀測器設計、最優控制等方面有著明顯的優勢。在這種情況下，如果能夠有效地辨識得到時變系統的狀態空間模型參數，這在系統控制等方面的應用有著重要的意義，所以時變結構的模型辨識問題值得關注。

針對狀態空間方程形式的系統的模型參數辨識，前人已經進行了許多研究。從20世紀80年代開始，隨著子空間系列辨識方法的提出，對時變狀態空間模型的分析有了進一步的發展[1]。子空間系列方法通過構建輸入輸出數據的結構子空間，利用奇異值分解(Singular Value Decomposition, SVD)或者QR分解的方法，辨識狀態空間模型的參數矩陣。在這基礎上，Ljung[2]將子空間方法進行了完善，并且根據不同的條件詳細討論了狀態模型參數的辨識問題。Verhaegen等[3]提出了多變量系統輸出誤差狀態空間模型辨識方法(MIMO Output-Error State Space Model Identification, MOESP)，其中概述總結了子空間辨識方法的主要流程，通過對輸入輸出數據進行子空間分解以確定系統的增廣能觀測矩陣，從而估計得到系統狀態空間模型的參數矩陣。Verhaegen等[4-5]還在此基礎上，進一步分析了直接和間接MOESP方法的性能。Juang等[6-7]在一般子空間方法的基礎上通過引入Markov參數，提出了特征系統實現算法(Eigensystem Realization Algorithm, ERA)，它可以認為是一種改進后的子空間方法，如果辨識的對象是一般的定常系統，這種方法已經被廣泛應用于工程中。如果系統是時變的，這一類子空間方法往往采用重復試驗的方式，利用多組數據進行系統的辨識。這類方法的缺點是如果系統的輸入輸出數據量很大，那么構建的Hankel矩陣的維數將會比較高，進行QR或SVD分解的計算量會非常大，從而導致這種方法對于復雜系統的計算效率較低。而遞推形式的子空間方法通過信號子空間的迭代，從而避免了較高維數的Hankel矩陣求解問題，具有計算時間短、計算效率高的優點，因此得到了越來越多的關注和廣泛應用。

在本文中，應用遞推子空間辨識方法來實現系統的狀態空間模型參數的辨識。第一類遞推子空間方法是通過遞推方法刷新系統的增廣能觀測矩陣[8]。Longman等[9]構建了多變量定常結構系統的狀態空間模型。Bosse等[10]提出了一種利用URV分解的遞推子空間方法。這一類的方法通過遞推形式得到各個時刻的Hankel矩陣，然后對其進行QR或SVD分解，由于仍需要逐一時刻進行子空間分解操作，所以計算量依然比較大。另一類遞推方法是基于信號子空間追蹤的概念，通過使用陣列信號處理算法來遞推刷新系統的能觀性矩陣[11-12]。龐世偉等[13]利用該方法辨識得到了移動質量-簡支梁的時變模態參數。Real等[14]為了提高辨識計算的速度，提出了快速估計子空間方法(Fast Approximate Subspace Tracking, FAST)。針對原有的遞推子空間方法中投信號投影子空間不嚴格正交的問題，Badeau等[15]對信號子空間的投影乘以修正后的正交矩陣，在此基礎上得到了逼近冪迭代(Approximated Power Iteration, API)方法，并且在此基礎上提出了快速API方法(Fast API, FAPI)[16]。丁峰等[17]提出了狀態空間模型遞階辨識方法，這種方法首先假設系統狀態是已知的，根據該狀態估計和系統輸入輸出數據遞歸計算系統的模型參數，然后基于系統輸入輸出數據和獲得的參數估計，遞歸計算系統的狀態估計。楊華等[18]針對時變系統的狀態模型參數辨識，將梯度型算法引入子空間跟蹤方法中，實現對廣義能觀測矩陣的估計，避免了子空間近似帶來的估計有偏性。與第一類遞推子空間方法相比，第二類方法計算效率更高，但是目前已知的工作主要集中在通過狀態空間方程中信號子空間矩陣的遞推來得到系統的模態參數值[13,19]，對于如何辨識得到完整的一組時變狀態空間模型{A(k),B(k),C(k)}并沒有提及。

本文對遞推子空間方法的應用作進一步的擴展，在前人利用遞推子空間方法辨識時變模態參數的基礎上，關注于時變的狀態空間模型的辨識。通過重構輸入輸出關系，實現了對系統的時變狀態空間模型的矩陣參數{A(k),B(k),C(k)}的識別。數值仿真選用二階彈簧-集中質量和旋轉機械臂連桿系統這兩個經典的參數時變模型，辨識系統的時變狀態空間模型。最后通過計算對比辨識得到的狀態模型和原系統模型的響應值，證明了利用該方法辨識得到的狀態模型能夠反映出原系統的響應特性，本文提出的遞推辨識狀態空間模型的方法是有效的。文中還將該方法與已有的基于SVD分解的時變ERA方法[20](Time-Varying ERA, TV-ERA)進行了計算效率的比較，結果證明與采用多組輸入輸出數據的重復試驗方法[20]相比，遞推方法在滿足精度要求的前提下，具有更高的計算效率。

1遞推子空間方法概述

1.1數據的預處理

子空間跟蹤算法在進行辨識之前，需要對狀態空間方程中的輸入輸出數據進行適當的前處理。利用系統的輸入輸出數據矩陣，構建對應時刻的狀態量變化，將其作為子空間跟蹤中的輸入信號參數，然后在子空間追蹤方法中不斷更新該狀態量，實現對系統的信號子空間(能觀性矩陣)的追蹤求解，進而得到時變系統的狀態空間模型。

這里首先考慮如下的離散時變系統：

x(k+1)=A(k)x(k)+B(k)u(k)

(1)

y(k)=C(k)x(k)

(2)

式中,系統的階數為n，輸入的維數為r，輸出維數為m。則輸入輸出矩陣的Hankel矩陣形式分別可以寫為

(3)

其中

下標k為對應的時刻，M為選取的合適的Hankel矩陣維數。利用文獻[10]中的輸入輸出數據的遞推更新方法，構建得到的狀態量z(k)的更新格式為

(4)

(5)

(6)

[U(k)UT(k)]-1=

[U(k-1)UT(k-1)]-1-

(7)

[Y(k)U?(k)]=

[Y(k-1)U?(k-1)]-z(k)wT(k)

(8)

1.2遞推API方法

在進行輸入輸出數據的預處理后，時變系統的辨識轉變成為子空間跟蹤問題。在多種遞推子空間方法中，本文在這里采用遞推API方法進行系統的遞推辨識。和經典的投影估計子空間跟蹤方法(Projection Approximation Subspace Tracking, PAST)相比，API方法的計算量略高于PAST方法，但是API方法在子空間的遞推方面增加了正交項，保證了投影矩陣的正交性，從而確保了遞推時的收斂性。

由于遞推子空間方法的介紹并非本文工作的重點內容，所以在這里只給出API方法的基本步驟，詳細的API方法的推導過程請參考文獻[15-16]：

s(k)=WH(k-1)z(k)

(9)

h(k)=Z(k-1)s(k)

(10)

g(k)=h(k)(β+sH(k)h(k))-1

(11)

e(k)=x(k)-W(k-1)s(k)

(12)

Θ(k)=(In+‖e(k)‖2g(k)g(k)H)-1/2

(13)

g(k)sH(k))Z(k-1)Θ-H(k)

(14)

W(k)=(W(k-1)+e(k)gH(k))Θ(k)

(15)

2時變狀態空間模型參數的辨識

2.1系統矩陣A(k)和輸出矩陣C(k)的辨識

式中上標“^” 為系統的辨識值。

(17)

(18)

2.2輸入矩陣B(k)的辨識

(19)

(20)

(21)

即

?(k)y(k)

(22)

那么每個時刻的遞推形式可以寫為

?(k)y(k)

(23)

將式(22)和(23)寫成廣義矩陣的形式：

(24)

式中

(25)

那么

(26)

(27)

(28)

式中：

圖1　遞推子空間方法辨識時變狀態空間模型的流程圖Fig.1 The flow chart of identify the time-varying state space model using the recursive subspace method

3數值仿真

為了驗證這種遞推方法對于狀態空間模型參數辨識的效果，數值算例首先選用一個二階彈簧-質量系統，如圖2所示，其中滑塊的質量m1=m2=1 kg，阻尼系數E1=1 Ns/m，E2=0.8 Ns/m。采樣時間Δt=0.000 2 s。輸入選為隨機激勵信號，質量塊的水平位移yi(i=1,2)作為輸出。定義輸出信號的信噪比(SNR)為

SNR=σr/σnr

(29)

式中：σr表示原輸出信號的標準差，σnr表示含有噪聲的輸出信號標準差。在仿真中，定義SNR=20，式(10)的Hankel矩陣中M=10。考慮剛度矩陣突變情況，即

(30)

圖2　兩自由度集中質量模型Fig.2 2DOF lumped mass model

圖3　原系統和辨識得到的系統的位移響應值比較(SNR=20)Fig.3 Comparison of displacement response value of the original and identified systems (SNR=20)

圖4則給出了在不同信噪比下(SNR分別為50、30和15時)、時間0～1 s內的辨識得到的狀態空間模型參數所對應的響應曲線，從圖中的結果可以看出，當信噪比較低(例如SNR=15時)，仍能夠有效計算得到系統的狀態空間模型；而當SNR逐漸增加時，辨識精度得到進一步的提高，響應曲線也更加接近理論值。圖4的結論說明，在不同的信噪比下能夠辨識得到時變系統的狀態空間模型參數，證明這種方法是具有一定的抗噪聲能力的。

圖4　在不同信噪比下的原系統和辨識得到的系統的位移響應值比較Fig.4 Comparison of displacement response value of the original and identified systems for different SNRs

為了進一步驗證這種方法對其它系統模型的辨識效果，再選用一個文獻[21]中的空間機械臂的例子，結構如圖5所示。原文中對系統的模態參數進行了辨識，這里我們利用本文提出的方法，辨識它的模型參數。u1和 u2為在關節上施加的力矩。關節1和關節2的阻尼系數為d1和d2，轉動剛度為k1和k2，均質連桿的質量均為m，長度均為l。考慮施加一個時變的作用力f(t)在連桿的自由端，令f(t)始終垂直于水平方向x軸，即與水平方向的夾角φ3=90°。連桿與水平方向的初始夾角為φ10和φ20。當連桿受到擾動時，連桿將會在它們的平衡位置附近振動，則實際角度為φ1=φ10+φ11和φ2=φ20+φ21。假定角度的振動是很小的，系統的一般形式的動力學方程可以寫為

(31)

式中：

Δφ0=φ10-φ20， Δφ1=φ10-φ3， Δφ2=φ20-φ3。

圖5　二連桿機械臂模型Fig.5 Model of two-link robotic manipulator

在仿真中，首先將式(31)改寫為狀態空間方程形式。算例中的各參數給定如下：l=1 m，m=1 kg，d1=d2=0.8 Nms/rad，k1=k2=80 Nm/rad，初始角φ10=0°，φ20=90°，φ3=90°。參數矩陣中的作用力f(t)為：

(32)

式中：f0=20 N，Δf=10 N。而輸入u(t)采用高斯分布的隨機激勵信號。輸出矩陣直接給定為C=I，即y=x，輸出信號的信噪比SNR=50，采樣時間Δt=0.01 s。利用本文提出的方法辨識系統的狀態空間模型，其中構建的Hankel矩陣中M=20，遺忘因子β=0.98。

圖6　原系統和辨識得到的系統的響應值比較Fig.6 Comparison of response value of the original and identified systems

Hankel行塊矩陣數目M計算時間/s本文遞推方法時變ERA方法M=107.7523.01M=209.8448.05M=3014.0183.62M=4016.12129.31

4結論

針對線性時變系統的狀態空間模型的辨識，本文從遞推子空間理論出發，提出了一種遞推辨識系統矩陣的格式。這種新的遞推格式通過重建狀態空間方程中的輸入輸出關系，主要思想是構建新的信號子空間，從而得到狀態空間模型參數矩陣。文中通過二階質量塊和旋轉機械臂的例子，驗證了本文提出的這一方法。最后在數值仿真中通過比較原系統和辨識得到的狀態空間模型的響應值，驗證了所辨識模型的精度。仿真算例結果表明本文提出的這種遞推辨識方法能夠有效地辨識線性時變系統的狀態空間模型參數。與采用多組數據的重復實驗方法相比，遞推方法具有更高的計算效率，數值仿真的結果證明了這一點。

參 考 文 獻

[ 1 ] 楊麗芳, 于開平, 龐世偉, 等. 用于線性時變結構系統辨識的子空間方法比較研究[J]. 振動與沖擊,2007,26(3): 8-12.

YANG Li-fang, YU Kai-ping, PANG Shi-wei, et al.Comparison study on identification methods applied to linear time varying structures [J]. Journal of Vibration and Shock, 2007, 26(3): 8-12.

[ 2 ] Ljung L. System identification: theory for the user [M].Upper Saddle River,NJ:Prentice Hall, 1999.

[ 3 ] Verhaegen M, Dewilde P. Subspace model identification part 1. The output-error state-space model identification class of algorithms [J]. International Journal of Control,1992,56(5): 1187-1210.

[ 4 ] Verhaegen M, Dewilde P. Subspace model identification part 2. Analysis of the elementary output-error state-space model identification algorithm [J]. International Journal of Control, 1992, 56(5): 1211-1241.

[ 5 ] Verhaegen M. Subspace model identification part 3. Analysis of the ordinary output-error state-space model identification algorithm [J]. International Journal of control,1993,58(3): 555-586.

[ 6 ] Juang J N, Pappa R S. An eigensystem realization algorithm for modal parameter identification and model reduction [J]. Journal of Guidance, Control, and Dynamics, 1985, 8(5): 620-627.

[ 7 ] 祝志文, 顧明. 基于自由振動響應識別顫振導數的特征系統實現算法[J]. 振動與沖擊, 2006, 25(5): 28-31.

ZHU Zhi-wen, GU Ming. Identification of flutter derivatives by using eigensystem realization algorithm [J]. Journal of Vibration and Shock, 2006, 25(5): 28-31.

[ 8 ] Oku H, Kimura H. Recursive 4SID algorithms using gradient type subspace tracking [J]. Automatica, 2002, 38(6): 1035-1043.

[ 9 ] Longman R W, Juang J N. Recursive form of the eigensystem realization algorithm for system identification [J]. Journal of Guidance, Control, and Dynamics, 1989, 12(5): 647-652.

[10] Bosse A, Tasker F, Fisher S. Real-time modal parameter estimation using subspace methods: applications [J]. Mechanical Systems and Signal Processing, 1998, 12(6): 809-823.

[11] Yang B. Projection Approximation subspace tracking [J]. IEEE Transactions on Signal Processing,1995, 43(1):95-107.

[12] 張家濱, 陳國平. 基于隨機子空間的遞推在線模態識別算法[J]. 振動與沖擊, 2009, 28(8): 42-45.

ZHANG Jia-bin, CHEN Guo-ping.Stochastic subspace based on-line recursive modal identification method [J]. Journal of Vibration and Shock, 2009, 28(8): 42-45.

[13] 龐世偉, 于開平, 鄒經湘. 用于時變結構模態參數識別的投影估計遞推子空間方法[J]. 工程力學, 2005, 22(5): 115-119.

PANG Shi-wei, YU Kai-ping, ZOU Jing-xiang. A projection approximation recursive subspace method for identification of modal parameters of time-varying structures[J]. Engineering Mechanics, 2005, 22(5): 115-119.

[14] Real E C, Tufts D W, Cooley J W. Two algorithms for fast approximate subspace tracking [J]. IEEE Transactions on Signal Processing, 1999, 47(7): 1936-1945.

[15] Badeau R, Richard G, David B, et al. Approximated power iterations for fast subspace tracking [C]//Signal Processing and Its Applications, 2003. Proceedings. Seventh International Symposium on. IEEE, 2003, 2: 583-586.

[16] Badeau R, David B, Richard G. Fast approximated power iteration subspace tracking [J]. IEEE Transactions on Signal Processing, 2005, 53(8): 2931-2941.

[17] 丁鋒, 蕭德云. 多變量系統狀態空間模型的遞階辨識[J]. 控制與決策, 2005, 20(8): 848-853.

DING Feng, XIAO De-yun.Hierarchical identification of state space models for multivariable systems [J]. Control and Decision, 2005, 20(8): 848-853.

[18] 楊華, 李少遠. 一種新的基于遺忘因子的遞推子空間辨識算法[J]. 控制理論與應用, 2009, 26(1): 69-72.

YANG Hua, LI Shao-yuan. A novel recursive MOESP subspace identification algorithm based on forgetting factor[J]. Control Theory and Applications, 2009, 26(1): 69-72.

[19] 楊凱, 于開平, 劉榮賀, 等. 兩種新的基于子空間跟蹤的時變模態參數快速辨識算法[J]. 工程力學,2012,29(10):294-300.

YANG Kai, YU Kai-ping, LIU Rong-he, et al.Two new fast identification algorithms of time-varying modal parameters based on subspace tracking [J]. Engineering Mechanics, 2012, 29(10): 294-300.

[20] Majji M, Juang J N, Junkins J L. Time-varying eigensystem realization algorithm [J]. Journal of Guidance, Control, and Dynamics, 2010, 33(1): 13-28.

[21] Liu K. Identification of linear time-varying systems [J]. Journal of Sound and Vibration, 1997, 206(4): 487-505.

Recursive identification for state space model of a linear time-varying system

NIZhi-yu1,WUZhi-gang1,2

(1. State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, China;2. School of Aeronautics and Astronautics, Dalian University of Technology, Dalian 116024, China)

Abstract:A novel recursive identification form for identifying state space model of a linear time-varying system was presented here. Being different from identification methods based on the singular value decomposition (SVD) and the least square estimation, the proposed recursive method was derived based on the signal subspace projection theory. The time-varying state space model of a system was obtained with the new signal subspace matrix by reconstructing the relation of input and output data. Comparing with the existing identification methods, the computation time of the proposed approach decreased because the recursive method did not require the SVD calculation. Particularly, when the system’s order was high, the advantage of computational efficiency for the recursive method was significant. In numerical simulation examples, the identified results and computational efficiency were compared with those using the classical time-varying eigensystem realization algorithm (TV-ERA) based on SVD. The simulation results showed that the proposed approach can be applied to identify the state space model of a linear time-varying system and it has a higher computational efficiency than TV-ERA does.

Key words:linear time-varying system; recursive subspace method; state space model; parametric identification

中圖分類號:O321；O324；TB123

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.04.002

通信作者吳志剛 男，教授，博士生導師，1971年生

收稿日期：2014-09-18修改稿收到日期：2015-03-11

基金項目：國家自然科學基金資助項目(11072044;11372056)；高等學校博士點基金資助項目(20110041130001)

第一作者 倪智宇 男，博士生，1985年生

E-mail：wuzhg@dlut.edu.cn