黃金峰,張合新,張 植
(1.第二炮兵工程學(xué)院自動(dòng)控制系,陜西西安 710025;2.第二炮兵駐二00廠軍代室,北京 100854)
近年來,子空間算法[1,2]因?qū)ο闰?yàn)知識(shí)要求較少以及數(shù)值計(jì)算上的優(yōu)勢,對多變量系統(tǒng)辨識(shí)廣泛適用,受到控制和辨識(shí)領(lǐng)域的廣泛關(guān)注.隨著對其研究的不斷拓展,遞推子空間辨識(shí)算法[3,4]也成為一個(gè)熱點(diǎn)研究問題,并取得了一些成果[5-7].
Mercère G等[5]基于投影近似子空間跟蹤(PAST)實(shí)現(xiàn)對SVD的更新,以遞推方式獲取擴(kuò)展能觀測矩陣并成功應(yīng)用于Wiener系統(tǒng)的辨識(shí),但由于投影的近似性,估計(jì)值不一定能收斂到真值.Oku H等[4]利用梯度型子空間跟蹤算法實(shí)現(xiàn)了廣義能觀測矩陣的遞推形式,但沒有考慮引入遺忘因子機(jī)制.楊華等[6]提出一種新的梯度型遞推子空間辨識(shí)算法,實(shí)現(xiàn)基于RLS-like遺忘因子數(shù)值Hankel矩陣的構(gòu)造,提高了遞推估計(jì)的辨識(shí)速度,但采用固定遺忘因子的方法,當(dāng)A的特征值發(fā)生變化時(shí),算法不能同時(shí)保證快速跟蹤和跟蹤效果.
本文以MOESP[5-7]中的OM算法為研究對象,針對一類時(shí)變系統(tǒng),提出一種變因子遞推子空間算法.該算法利用系統(tǒng)矩陣A的特征值空間歐氏距離[8,9]信息實(shí)現(xiàn)變因子的步驟,在保證跟蹤效果的條件下提高了算法對參數(shù)變化的跟蹤能力.同時(shí)針對變動(dòng)檢測中子空間跟蹤問題中W的列向量可能出現(xiàn)不相交的情況,運(yùn)用OPAST算法遞推估計(jì)廣義能觀測矩陣,保證了W的列向量收斂于主子空間的正交基上.仿真結(jié)果驗(yàn)證了算法的有效性.
考慮一類子空間模型:
式中:u(t)∈Rm,y(t)∈Rl,x(t)∈Rn分別為輸入、輸出和狀態(tài)向量;w(t)∈Rn,v(t)∈Rl分別為過程和測量噪聲.為保證系統(tǒng)的可辨識(shí)性和收斂性,假設(shè)滿足(A,B)能控,(A,C)能觀測,系統(tǒng)為最小實(shí)現(xiàn);且外部輸入和噪聲不相關(guān),輸入滿足充分激勵(lì)條件.同時(shí)為研究簡便,假設(shè)w(t)=0,即狀態(tài)方程中不含過程噪聲.
子空間辨識(shí)的本質(zhì)就是直接從數(shù)據(jù)對中估計(jì)廣義能觀測矩陣,為此Verhaegen.M等設(shè)計(jì)了一系列的MOESP算法.本文在其中的OM算法(Ordinary MOESP)的基礎(chǔ)上設(shè)計(jì)了基于歐氏空間變因子的遞推子空間辨識(shí)算法,其步驟如下:
Step1:構(gòu)造帶有遺忘因子的Hankel矩陣:
給定初始輸入輸出數(shù)據(jù)對{ui,yi}(i=1,2,…,N),構(gòu)造初始Hankel矩陣(以輸入Hankel矩陣為例)[10]
據(jù)式(3)形式,系統(tǒng)輸入輸出有如下關(guān)系:
Step2:OPAST[11]法更新廣義能觀測矩陣及估計(jì)系統(tǒng)參數(shù)矩陣:
在0(N序列)步,對數(shù)據(jù)矩陣進(jìn)行QR分解
遞推第1步,獲取新的數(shù)據(jù)對后,對更新的增廣數(shù)據(jù)矩陣做QR分解
經(jīng)一系列Givens旋轉(zhuǎn),則式(8)R部分變?yōu)?/p>
R22(1)可以利用數(shù)據(jù)更新QR分解獲得_
為實(shí)現(xiàn)廣義能觀測矩陣子空間的重構(gòu),首先對R22(0)進(jìn)行SVD分解,則有
為遞推更新(j),運(yùn)用子空間跟蹤問題,形成如下無約束最優(yōu)化問題求解[12]
基于Yang的準(zhǔn)則,形成了不同方法的遞歸計(jì)算W(j),如輔助變量PAST算法[13]和 IV-PAST算法[5].理論上,最小化準(zhǔn)則J(W)中W的列向量是正交的,PAST算法指數(shù)級(jí)收斂于主子空間的正交基上;但在某些情況下,PAST算法將會(huì)在兩個(gè)矩陣間來回?cái)[動(dòng),這樣就不會(huì)收斂于主子空間的正交基.文獻(xiàn)[11]給出了這個(gè)問題的一個(gè)證明.如果W列向量不是正交的,這就表示在估計(jì)的過程中,不能夠保證在同一狀態(tài)空間坐標(biāo)下的估計(jì)一致性,尤其是運(yùn)用于變動(dòng)檢測上[3].為此,本文采用文獻(xiàn)[11]提出的正交化PAST算法(簡稱OPAST算法),運(yùn)用矩陣求逆引理得到W的最小二乘方法計(jì)算步驟
把上述的向量zj用Givens變換后的觀測向量(j)代替,獲得的W矩陣就是子空間辨識(shí)所需的廣義能觀測矩陣估計(jì)(j).變遺忘因子機(jī)制下輸入輸出Hankel矩陣不能直接用因子矩陣和數(shù)據(jù)矩陣乘積描述,應(yīng)按式(14)近似處理.
因子.則有該機(jī)制下Hankel矩陣存在如下關(guān)系:
則矩陣估計(jì)^C(j)可直接從(j)中提取
由式(15)可知,估計(jì)A^(j)有如下關(guān)系:
當(dāng)遺忘因子β為固定值β0時(shí),M(j)方陣就退化為.而矩陣B,D的估計(jì)值包含在矩陣(j)中,為得到其估計(jì)值,需從式(15a)中消去(j)(j)項(xiàng).定義矩陣(j)為(j)的正交補(bǔ),式(15a)兩邊左乘PG(j)右乘)可得
求解系統(tǒng)矩陣^B(j)和 ^D(j)估計(jì)的余下步驟見文獻(xiàn)[10].
Step3:固定尺度歐氏距離準(zhǔn)則更新遺忘因子:
采用遺忘因子方法進(jìn)行遞推計(jì)算的準(zhǔn)則是:對快時(shí)變參數(shù)選較小的遺忘因子,而對慢時(shí)變參數(shù)選較大的遺忘因子.基于該準(zhǔn)則,本文利用系統(tǒng)矩陣A的特征值空間歐氏距離大小來判斷參數(shù)變化的快慢,但運(yùn)用相鄰估計(jì)距離差會(huì)因變因子更新頻繁導(dǎo)致辨識(shí)精度變差.因此本文采用固定區(qū)間[(k-1)p+1,kp]內(nèi)最大值和最小值的距離差來判斷參數(shù)變化的快慢從而避免了變因子變化過于頻繁,提高了算法的穩(wěn)定性.令其中p表示固定步長,k=1,2,3,….即
這里λi為矩陣A的第i個(gè)特征值.利用上述的距離信息來更新遺忘因子 βj
式中:n為矩陣A的維數(shù);c為弱化因子.βmin為遺忘因子的下限值.
仿真1:建立時(shí)不變仿真模型,選擇系統(tǒng)矩陣參數(shù)為
仿真中采用方差為1的零均值白噪聲作為激勵(lì)信號(hào),由四階Butterworth濾波器產(chǎn)生,截止頻率為0.8倍奈奎斯特頻率,輸出被不相關(guān)的零均值白噪聲污染,信噪比為 25 dB.仿真參數(shù)為:N=20,i=3,βmin=β0=0.9.分別選擇固定尺度p=2,5,20進(jìn)行仿真.
圖1表明:增加固定尺度p,變因子β越晚接近于1,跟蹤越快,但跟蹤效果越易受噪聲影響.反之則 β越早接近于1,跟蹤能力越慢,但跟蹤效果越不易受噪聲影響.仿真結(jié)果顯示慢變系統(tǒng)辨識(shí)選擇p=5較為合適.
仿真2:建立慢變模型,k<1 000時(shí)模型參數(shù)同仿真1.k≥1 000時(shí),A的參數(shù)為
分別運(yùn)用變因子法、不加遺忘因子和固定因子0.9進(jìn)行了仿真.仿真結(jié)果如圖2所示.
圖1 時(shí)不變系統(tǒng)不同尺度跟蹤曲線Fig.1 Different scales’parameters of time invariant system
圖2 時(shí)變模型不同方法辨識(shí)參數(shù)曲線Fig.2 Different identification methods parameters of Time-varying model
仿真2結(jié)果表明:對慢變系統(tǒng)定因子0.9法雖能夠較好地跟蹤曲線變化,但跟蹤精度差;無遺忘因子的穩(wěn)定性較好,但在參數(shù)發(fā)生變化時(shí),跟蹤速度慢;而本文研究的固定尺度變因子法不但保持了較好的跟蹤速度而且保持了較好的跟蹤效果.
本文針對傳統(tǒng)的OM子空間辨識(shí)算法采用固定遺忘因子方法進(jìn)行變動(dòng)檢測時(shí)不能同時(shí)保證快速跟蹤和跟蹤效果等問題,提出了一種基于歐氏距離的變因子遞推子空間辨識(shí)算法.該算法分為 3個(gè)步驟:首先引入變因子構(gòu)造與更新Hankel矩陣和觀測向量;其次為保證廣義能觀測矩陣的列向量收斂于主子空間的正交基上,采用OPAST算法遞推估計(jì)廣義能觀測矩陣,并由廣義能觀測矩陣估計(jì)系統(tǒng)參數(shù)矩陣;最后用的特征值空間距離信息實(shí)現(xiàn)變因子,此算法具有的自適應(yīng)能力.應(yīng)用于一類時(shí)變系統(tǒng),仿真結(jié)果驗(yàn)證算法的有效性.
[1]Trnka P.Subspace like identification incorporating prior information[J].Automatica,2009,46(4):1046-1091.
[2]楊華,李少遠(yuǎn).一種完全數(shù)據(jù)驅(qū)動(dòng)的子空間辨識(shí)與魯棒預(yù)測控制器設(shè)計(jì)[J].控制理論與應(yīng)用,2007,24(5):732-736.
Yang Hua,LiShaoyuan.A novel robust predictive controller design based on data-driven subspace identification[J].Control Theory&Applications,2007,24(5):732-736.(in Chinese)
[3]Mercère G,Bako L,et al.Propagator-based methods for recursive subspace model identification[J].Signal Processing,2008,88(3):468-491.
[4]Oku H,Kimura H.Recursive 4SID algorithms using gradient type subspace tracking[J].Automatica,2002,38(6):1035-1043.
[5]Lovera M,Gustafsson T,Verhaegen M.Recursive subspace identificationof linear and non-linear Wiener state-space models[J].Automatica,2000,36(11):1639-1650.
[6]楊華,李少遠(yuǎn).一種新的基于遺忘因子的遞推子空間算法[J].控制理論與應(yīng)用,2009,26(1):69-72.
Yang Hua,Li Shaoyuan.A novel recursive MOESP subspace identification algorithm based on forgetting factor[J].Control Theory&Applications,2009,26(1):69-72.(in Chinese)
[7]Bontempi G,Birattari M,Bersini H.Lazy learning for local modeling and control design[J].International Journal of Control,1999,72(7/8):643-658.
[8]鄭建柏,朱永利.基于歐氏聚類和支持向量機(jī)的變壓器故障診斷[J].電力科學(xué)與工程,2008,24(4):13-15.
Zheng Jianbai,Zhu Yongli.Transformer fault diagnosis based on Euclidean clustering and support vector machines[J].Electric Power Science and Engineering,2008,24(4):13-15.(in Chinese)
[9]Zheng W X.Parameter Estimation of Stochastic Linear Systems with Noisy Input[J].International Journal of Systems Science,2004,35(3):185-190.
[10]Zhang C,Bitmead R.Subspace system identification for training-based MIMO channel estimation[J].Automatica,2005,41(9):1623-1632.
[11]Abed-Meraim K,Chkeif A,Hua Y.Fast orthonormal PAST algorithm[J].IEEE Signal Processing letters,2000,7(3):60-61.
[12]Yan F B.Asymptotic convergence analysis of the projection approximation subspace tracking algorithms[J].Signal Processing,1996,50(1/2):123-136.
[13]Gustafsson T.Instrumental variable subspace tracking using projection approximation[J].IEEE Transactions on signal processing,1998,46(33):669-681.