时变系统的FIR模型辩识及设计变量

文章编号:1000-5641(2001)02-0001-08
时变系统的FIR 模型辩识及设计变量
丁肇红1
, 袁震东
2
(1.上海应用技术学院自动化系,上海 200233; 2.华东师范大学数学系,上海 200062)摘要:文章考虑时变系统辨识中辨识设计变量的最优化问题。

时变系统的参数可用平稳的AR(1)随机过程描述,时变系统的FI R 模型传递函数估计的均方误差,记作M SE,它可以用一个较简单表达式来逼近它。

作者在M SE 极小的意义下,给出辨识设计变量的解析表达式。

关键词:时变系统; 系统辨识; 设计变量; 均方误差中图分类号:O231.3 文献标识码:A
在系统辨识的算法中有些参数可以看作设计变量,例如随机递度算法的步长因子,如拓广Kalman 滤波中误差协方差的比。

这些变量选得好,可以提高系统辩识的精度。

本文讨论慢时变系统辩识中设计变量的最优选取问题。

通常假设时系统的参数在随机游动情行下,利用最小均方(LMS)算法、RLS 、Kalman 滤波辨识时变系统的参数[1-3],本文讨论时变系统的参数是AR(1)随机过程时,利用递推算法辨识时变参数,求出FIR 模型的时变系统的传递函数估计的M SE 近似表达式,然后利用这个近似表达式求解辨识设计变量的最优解。

1 问题的归结
1.1 系统的描述
时变离散时间系统
y (t)=G t-1(z )u(t)+e(t )
(1)其中,u(t),y (t )分别是输入、输出,e(t)是零均值白噪声,G t -1(#)为时变系统传递函数
Ee 2(t )=r 2
(2)
设传递函数对固定的t 是稳定的,可以用FIR 模型表述
G t (z )=
6
n
k =1
g 0t (k )z
-k
(3) 系统(1)可以表示成FIR 结构,即
y (t)=U T
(t )H (t -1)+e(t)
(4)
其中
U (t)=(u (t -1),,,u (t -n))T
H (t)=(g 0t (1),,,,g 0t (n))
T
收稿日期:1999-05
作者简介:丁肇红(1966-),女,硕士,讲师.
第2期2001年6月
华东师范大学学报(自然科学版)
Journal of East China Normal University (Natural Science)
No.2 Jun.2001
定义
W n(X)=(e i X,e i2X,,,,e i n X)T(5)在时刻t的传递函数为
G t(e i X)=W*n(X)H(t)(6)
一般来说,时变系统是没有传递函数可言的,但如果系统参数变化的阶次大大低于系统变化的阶次,那么G t(e i X)可以看作时变系统在t时刻的快照。

假设时变系统的参数按AR(1)变化
H(t+1)=A H(t)+X(t+1)(7)其中,{X(t),t=1,2,}是独立的随机向量序列,均值为零。

E X(t)X T(t)=R1(t)
E X(t)e(s)=0P t,s
|A|[1
1989年S.Gunnarsson和L.Liung讨论了A=1时变系统的辨识设计,本文讨论|A|[1一般情形下的辨识设计问题。

由(7)式传递函数可写为
G t+1(e i X)=A G t(e i X)+$G t+1(e i X)(8)其中,$G t+1(e i X)=W*n(X)X(t+1)
令E[$G t(e i X)]2=r1(X,t)(9)那么r1(X,t)=W*n(X)R1(t)W n(X)(10)其中,*表示复共轭转置。

1.2输入信号的性质
假设输入信号{u(t)}是确定性的,有界且拟平稳信号,即
lim N y]1
N
6n
t=1
u(t)=0(11)
lim N y]
1
N
6n
t=1
u(t)u(t-S)=R u(S)(12)
5u(X)=6-]
S=]
R u(S)e-i XS(13)
Q>lim
N y]
1
N
6n
i=1
U(t)U T(t)(14)
2传递函数估计的平均二次误差(M SE)的表达式
给定数据{y(k),u(k),k<=t},使用递推算法可将H(t)的估计记为^H(t),从而有
G^t(e i X)=W*n(X)^H(t)(15) G^t(e i X)的好坏依赖于算法的设计变量。

递推算法如下:
^H(t)=A^H(t-1)+r(t)R-1(t)h(t,^H(t-1),E(t
E(t)=y(t)-U T(t)^H(t-1)
R(t)=R(t-1)+r(t)H(t,R(t,R(t-1),U(t)(16)
2华东师范大学学报(自然科学版)2001年
其中
|A |[1
h(t ,^H (t -1)),E (t),U (t)=U (t)E (t).(A )随机梯度算法
^H (t)=A ^H (t -1)+r (t)U (t)E (t)R (t)S I
H (t ,R (t -1),U (t))S 0
(17)
(B)Kalman 滤波算法
^H (t)=A ^H (t -1)+
1A 2
^r 2
[P(t)-r 2
(t)R ^1(t )]U (t)E (t)(18)P (t)=A 2
P(t -1)-P (t -1)U (t)U T
(t)P (t -1)U T (t )P (t -1)U (t)+^r 2
+r 2
(t)R ^1(t)(19)
下面考虑传递函数估计的均方误差(MSE):
P n (X ,t )>E |G
^t (e i X )-G t (e i X )|2(20)时变参数估计误差
H (t)= H (t)-H (t)
(21)利用(7)、(16)式得
H (t)=(A I -r (t)R -1(t)U (t )) H (t -1)+r (t)R -1(t)U (t )e(t)-X (t)(22)令
Z(t )=r (t)R
-1
(t)
(23)其中对不同算法,Z(t )可用不同的表达式,如Z(t)=rI
(24)或
Z(t )=
1
A 2
^r 2
[P(t)-r 2R ^1(t)](25)
则
H (t )=(A I -r (t )Z (t)U (t)U T (t)) H (t -1)+Z (t)U (t)e(t)-X (t )
(26) 令协方差矩阵 n (t)=E [ H (t) H T
(t)]
(27)
由(26)(27)式得
n (t)=(A I -Z (t)U (t )U T (t)) n (t -1)(A
I -Z(t )U (t)U T
(t ))T +r 2Z(t)U (t )U T
(t)Z T (t)+R 1(t )
(28)由(6),(15),(21)式得
G t (e i X
)=W *
n (X )(^H (t)-H (t))
(29)因而传递函数辨识的M SE 为:
P n (X ,t)=W *n (X ) n (t )(t)W n (X )(30)
由(28)(30)式及z (t)的表达式,可以得到MSE 的精确表达式,但非常复杂,下面用个
近似表达来刻划M SE 。

3 MSE 的近似表达式
3.1 随机梯度算法
3
第2期丁肇红,等:时变系统的FI R 模型辩识及设计变量
讨论系统变化缓慢和参数自校正缓慢时, n(t)的渐近必当Z(t)=rI时, n(t)的表达式可以重新表述。

当系统变化缓慢,即R1(t)小且当自适应缓慢,即r小。

定义
n(t)=A2 n(t-1)-A r Q n(t)-A r n(t-1)Q+r2r2Q+R1(t)(31)注意 n(t)和 n(t)都依赖于r和R1(t)。

由文献[2,4]启发,有下面的引理。

引理1若 n(t)、 n(t)由(28)式及(31)式给出,当r和+R1(#)+趋向零时,有(r++R1(#)+)-D+ n(t,r,R1(#)- n(t,r,R1(#)+y00([D[1)(32)证明见附录。

从引理1知,当r y0且+R1(#)+y0时, n(t)可以用 n(t)近似,因而进一步有
P n(X,t)U W*n(X) n(t)W n(X)(33)为了进一步简化P n(X,t)表达式,首先需引入矩阵类N的定义[4]。

引入方程
^P n(X,t)=(A2-2A r5u(X))^P n(X,t-1)+nr2r25u(X)+r1(X,t)(34)定理1假设^P n(X,t), n(t)分别满足(34)和(31)式,又假设R1(t), n(0)是T oeplitz结构,当n趋向无穷大时,它们是扩充的Toeplitz结构且R1(t), n(0)属于N,那么1
n
|^P n(X,t)-W*n(X) n(t)W n(X)|y0.n y](35)证明见附录。

由定理1及(33)式知
P n(X,t)U^P n(X,t)(36)现在就得到了在随机梯度算法下的MSE的近似表达式,这个近似表达式是在缓慢自适应和高介模型下得到的,在定理1中n可以看作形式参数而非系统模型的阶次,所以定理1可视为数学陈述的代数形式,对于较大的且有限的n,可以用式(34)来近似。

3.2Kalman滤波算法
令
r(t)S r,Z(t)=
P(t-1)
^r2+U T(t)P(t-1)U(t)
类似随机梯度算法,考虑系统变化缓慢和自适应缓慢时 n(t)的性质,在高阶模型下可以得到M SE的近似表达式。

P(t)能用
P(t)近似且
P(t)满足
P(t)=A2[ P(t-1)-1
^r2
P(t-1)Q
P(t-1)+r2R^1(t)](37)
Z(t)可以近似Z(t)且
Z(t)满足
Z(t)=
P(t)
^r2
(38)
进一步 n(t)可以用 n(t)近似, n(t)满足差分方程
n(t)=A2 n(t-1)-A Z(t)Q n(t-1)
-A n(t-1)Q Z(t)+r2 Z(t)Q Z(t)+R1(t)(39)定义^r1(X,t)=W*n(X)R^1(t)W n(X)(40) 4华东师范大学学报(自然科学版)2001年
引入方程:
^P n (X ,t )=(A 2
-2A 5u (X ) Z (X ,t))^P n (X ,t -1)
+nr 25u (X ) Z 2
(X ,t)+r 1(X ,t))(41)其中 Z (X ,t )=A 2
Z (X ,t -1)-A 2
5u (X ) Z 2
(X ,t -1)+r 2
n^r 2^r 1
(X ,t)
(42)
定理2 假设^P n (X ,t), n (t)分别满足(41),(39)时,又假设R 1(t),R ^1(t)是T oeplitz 矩阵, n (0), Z (0I N ),则
1n
|^P n (X ,t)-W *n (X ) n (t)W n (X )|y 0,n y ](43)
证明见附录。

若^r 1(X ,t )不依赖于时间t,即^r 1(X ,t)=^r 1(X )时,(42)式的平稳解为
Z (X )=1-1A 2125u (X )+1-1A 22125u (X )2
+1A 25u (X )
r 2
^r 1(X )n^r 2(44)利用Kalman 滤波法辨识时变系统的参数,类同随机梯度法,在缓慢自适应和高阶模型
下,得到了M SE 的近拟表达式,它满足了一个简单差分方程,下面将讨论^P n (X ,t )如何依赖于算法的设计变量。

4 最优设计变量
以上给出了M SE 的近似表达式,^P n (X ,t)满足一个简单的差分方程。

M SE 如何依赖于设计变量就转化为^P n (X ,t)如何依赖于算法的设计变量。

仍然讨论时变系统变化是缓慢的,假设真实系统的随机变化具有平稳性,即r 1(X ,t)与时间t 无关,记为r 1(X ,t)=r 1(X )。

4.1 随机梯度算法
^P n (X ,t )可以表述为:
^P n (X ,t )=(A 2-2A r 5u (X ))^P n (X ,t -1)+nr 2r 25u (X )+r 1(X )
差分方程的平稳解为
^P n (X )=nr 2r 2
5u (X )+r 1(X )
1+2A r 5u (X )-A
2
(45)要判断传递函数估计的均方误差这个准则的好坏情况,只需使得MSE 的近似表达式的平稳解^P n (X )尽可能的小,此处,算法的设计变量是步长因子r ,考虑min r
^P n (X )对^P n (X )求导,即5^P n (X )
5r =0,得 r
op t
=125u (X )A -1
A +A -1
A
2
125u (X )
2
+
r 1(X )nr 25u (X )
(46)
且当r =r o p t
时,52^P n (X )
5r
2
>0,因而r o p t 是最优的。

注意:r op t
对所有频率X 并非最优的,它仅对固定X 才是最优的。

4.2 Klaman 滤波算法
^P n (X ,t )可表述为:
5
第2期丁肇红,等:时变系统的FI R 模型辩识及设计变量
^P n(X,t)=(A2-2A5u(X)
Z(X,t))^P n(X,t-1)+nr25u(X)
Z2(X,t)+r1(X)(47)将
Z(X,t)的平稳解(44)代入(47)式可得^P n(X,t)的平稳解如下:
^P n(X)=nr25u(X) Z2(X)+r1(X)
1+2A r5u(X)
Z(X)-A2
(48)
其中,
Z(X)满足(44)式
要使得^P n(X)极小,可得
Z o p t(X)=
1
25u(X)
A-1A+A-1A
21
25u(X)
2
+
r1(X)
nr25u(X)
(49)
这个算法中需设计变量是r2^r1(X)和^r2,而在Kalm an滤波算法中已转化为依赖于比值r2^r1(X)
^r2,令
^F=r2^r1(X)
^r2
(50)
比较(44)式和(49)式可得,当^F op t(X)=r1(X)
r2
时,即Z^op t(X)=A Z^(X)时,可以使得
^P n(X)极小,也就是^F op t(X)=r1(X)
r2
是所求得的最优设计变量。

Klam an滤波算法所依赖的设计变量是^F(X),对所有频率,X,^F op t(X)都是最优的,这和随机梯度法不一样。

5结束语
本文对FIR结构时变系统,在缓慢时变,缓慢自适应和系统阶次充分高的情形下,得到了M SE的近似表达式,这个近似表达展示了输入谱、噪声方差、模型阶、系统变化以及算法的步长因子等设计变量对MSE的影响。

同时,它的平稳解^P n(X),也给出了干扰灵敏度和跟踪性能的权衡,最后,对两种不同的递推算法,分别给出了固定X的最优步长因子和对所有频率X的最优参数^F(X)解。

[参考文献]
[1]袁震东著.自适应控制理论及应用[M].上海:华东师范大学出版社,1988.
[2]S Gunnars son,L Ljung.Frequency domain tracking characteristics of adaptive algori thms[J].IEEE T rans On Acous-
tic,S peech and Sign PROC.1989,37(7):1072~1089.
[3]A Lsaksson.Identification of time varying sys tem through adaptive Kalman Filtering[A].IPC Com mi ttee of IFAC.In:
Preprints10th IFAC w orld congress[C].M unich:Pergm on Pres s,Vol.X,1987,305~311.
[4]S Gunnarsson.Frequency domai n aspects of modeling and control in adaptive system,Linkopi ng Studies in Sci ence and
T echnoi ogy T hesis[D].No194,Linkoping University,S-581,83Sw eden,1988.
6华东师范大学学报(自然科学版)2001年
Identification of FIR Model and Design Variable for Time -varying System
DING Zhao -hong 1
, YU AN Zhen -dong
2
(1.S hanghai College o f App lie d T ec h nology ,Sha nghai 200233,China ;
2.Department of M athematics,E ast China Nor mal Univer sity ,Shanghai 200062,China)
Abstract: In this art icle t he time-varying system is described as a FIR model where the parameters modelled as AR(1)process.T he mean square error of t ransfer function estimat ions for time-varying F IR model is denoted as M SE.T he approx imat ion for MSE is obtained and used t o solve the optimiza -t ion problem for design variables.
Key wor ds: time-varying syst em; syst e m ident ific ation; design variable; M SE(mean square error)
附 录
引理1证明:令
n (t )=r n (t )
由(28)得
n (t )=A 2
n (t -1)-A
r U (t)U T
(t) n (t -1)-A r n (t -1)U (t )U T
(t )+r 2U (t )U T (t ) n (t -1)U (t )U T (t )+r r 2U (t)U T (t)+
1
r R 1
(t )(51) n (t )= n (t -1)+r [H 1( n (t -1),U (t))+H 2(U (t ),R 1)](52)
其中
H 1=A 2-1
r
n (t -1)-A U (t)U T (t) n (t -1)-A n (t -1)U (t )U T (t )+r U (t )U T (t ) n (t -1)U (t )U T (t )
(53)H 2=r 2U (t )U T (t )+1
r 2
R 1(t )(54)
则
EH 1=
A 2-1
r
n -A Q n -A n Q (55)
E H 2=r 2Q +1
r
2R 1(56)
因而(52)中的 n (t)渐近趋于 D n ,且 D
n 满足连带微分方程:
d d S D n (S )=A 2-1r D n (S )-A Q D n (S )-A D
n (S )Q +r 2Q +1r 2R 1(57)
令 n (t)=r U n (t )则
7
第2期丁肇红,等:时变系统的FI R 模型辩识及设计变量
U
n(t)=U
n(t-1)+r
A2-1
r
U
n(t-1)-A Q
U
n(t-1)-A
U
n(t-1)Q+r2Q+
1
r2
R1(t)(58)
因而(58)的解U
n
(t)渐近趋于 D n,且 D n满足连带微分方程(57)。

所以
+ n(t)- n(t)+=r+ n(t)-U
n (t)+=r+( n(t)- D n)+( D n-U
n
(t))+
[(r++R1(Ó)+(+ n(t)- D n+++U
n
(t)- D n+)即
(r++R1(Ó)+-D+ n(t)- n(t)+
[(r++R1(Ó)+1-D(+ n(t)- D n+++U
n
(t)- D n+)
y0,当r y0,+R1(Ó)+y0时(0[D[1)命题得证。

定理1证明:由文献[2]附录的引理A.1及A.2有
n(t)I N P t
lim n y]1
n
W*n(X) n(t)W n(X)= P(X,t)
定义
^P n(X,t)=n P(X,t)因此在(31)式两边作用极限算子
lim n y]1
n
W*n(X)(Ó)W n(X)
根据[2]的引理A.3即可证得命题且
^P n(X,t)
满足(34)式。

定理2证明:由文献[2]附录的引理A.1及A.2有
n(t)I N, Z(t)I N P t
lim n y]1
n
W*n(X) n(t)W n(X)= P(X,t)
lim n y]1
n
W*n(X)Z n(t)W n(X)= Z(X,t)
由(37),(38),(39)两边作用极限算子
lim n y]1
n
W*n(X)(Ó)W n(X)
根据[2]的引理A.3即可证得命题且
^P n(X,t)和 Z n(X,t)
满足(41)和(42)式。

8华东师范大学学报(自然科学版)2001年。