线性时变系统的状态空间模型递推辨识研究

倪智宇， 吴志刚,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