于 淼 劉建昌 郭 戈

(1.東北大學秦皇島分?？刂乒こ虒W院,河北秦皇島 066004;2.東北大學信息科學與工程學院,遼寧沈陽 110819;3.東北大學流程工業(yè)綜合自動化國家重點實驗室,遼寧 沈陽 110819)

1 引言

近年來,子空間辨識作為系統(tǒng)辨識的重要分支在控制領(lǐng)域受到了廣泛的關(guān)注[1–3].針對離線模型的辨識問題,子空間辨識方法已具有較好的性能.與離散時間系統(tǒng)的辨識相比較,連續(xù)時間系統(tǒng)的辨識在工業(yè)工程、航空航天、生物學、醫(yī)學及社會經(jīng)濟學等眾多領(lǐng)域應用廣泛[4].通過連續(xù)時間系統(tǒng)模型的辨識,可以較準確地判斷復雜的連續(xù)時間系統(tǒng)的現(xiàn)象和過程的內(nèi)在規(guī)律.此外,多數(shù)物理現(xiàn)象都具有連續(xù)屬性,描述它們的數(shù)學模型是微分方程[5].因此,針對連續(xù)時間系統(tǒng)的辨識問題進行深入研究,具有重大的理論意義和實際價值.

采用子空間辨識方法建立的離線模型不能有效跟蹤系統(tǒng)的動態(tài)變化.利用輸入輸出數(shù)據(jù)的更新不斷地更新所建立的模型,這樣才能及時地遞推更新模型.遞推辨識可以不斷地從系統(tǒng)獲得更新數(shù)據(jù),從而提高辨識的模型精度直至達到辨識要求.遞推辨識過程中主要使用正交三角分解(QR decomposition,QRD)和奇異值分解(singular value decomposition,SVD)等線性代數(shù)工具,雖然提高了算法的數(shù)值魯棒性,但也相應增加了子空間辨識的在線遞推困難[6–7].所以如何降低遞推子空間辨識方法的計算量是首要問題.文獻[4]提出了基于傳播算子方法的遞推預測子空間辨識方法.對于存在負載擾動的系統(tǒng),文獻[5]采用麥克勞林時間序列逼近帶有負載擾動的輸出響應,從而提出了帶有偏差消除估計性能的子空間辨識方法.在文獻[8]中,作者通過旋轉(zhuǎn)當前的信號子空間預測將來的子空間,從而遞推更新系統(tǒng)矩陣.文獻[9]提出了一種基于可變遺忘因子的穩(wěn)定遞推規(guī)范變量狀態(tài)空間方法.結(jié)合規(guī)范相關(guān)分析方法,文獻[10]提出了帶有遺忘因子的數(shù)據(jù)驅(qū)動遞推子空間辨識方法.文獻[11]研究了兩種基于遞推最小二乘的子空間辨識方法,并且分別獲得帶有負載擾動系統(tǒng)的可觀測馬爾可夫參數(shù)矩陣和負載擾動響應.針對線性參數(shù)變化系統(tǒng)(linear parameter varying,LPV),文獻[12]提出了基于預測的張量回歸子空間辨識方法,此方法保持了LPV馬爾可夫參數(shù)矩陣的固有結(jié)構(gòu)并且避免了維數(shù)災難.

上述子空間辨識方法解決了在線遞推辨識問題,它們大多研究的是系統(tǒng)階數(shù)已知的離散時間系統(tǒng)遞推更新.然而,連續(xù)時間系統(tǒng)的辨識與離散時間系統(tǒng)的辨識相比較,在參數(shù)估計、H∞控制和自校正控制等方面有更加廣泛的應用.在文獻[13]中,作者研究了基于拉蓋爾濾波器的遞推連續(xù)時間子空間辨識方法,其中采用拉蓋爾濾波器解決了連續(xù)時間系統(tǒng)中時間導數(shù)的問題.針對連續(xù)時間系統(tǒng)的辨識問題,文獻[14]提出了基于拉蓋爾濾波器的核范數(shù)子空間辨識新方法.文獻[15]采用一組泊松矩函數(shù)的線性濾波器解決了連續(xù)時間系統(tǒng)的導數(shù)和微分問題,結(jié)合泊松矩函數(shù)和子空間辨識方法提出了辨識新方法.文獻[16]提出了基于復雜頻域分析的子空間辨識方法解決了連續(xù)時間系統(tǒng)的辨識問題.采用隨機分布理論方法可以解決系統(tǒng)時間導數(shù)問題,此方法具有更少的參數(shù),并且避免了辨識過程中矩陣對數(shù)的復雜運算[6].在實際生產(chǎn)過程中,采用子空間辨識法建立的離線模型并不能有效準確地跟蹤系統(tǒng)的動態(tài)變化.采用奇異值分解等線性代數(shù)工具,增加算法數(shù)值魯棒性的同時,也相應增加了子空間辨識的在線遞推困難.并且過程噪聲和測量噪聲的存在給連續(xù)時間系統(tǒng)的在線辨識帶來了一定的困難.因此,當系統(tǒng)階數(shù)未知時,采用遞推子空間辨識方法解決連續(xù)時間系統(tǒng)的辨識問題是一個挑戰(zhàn).本文針對以上連續(xù)時間系統(tǒng)辨識過程中存在的具體問題,提出一種基于隨機分布理論的遞推子空間辨識方法.首先,采用隨機分布理論重新構(gòu)建連續(xù)時間系統(tǒng)的輸入輸出Hankel矩陣,解決了噪聲存在帶來的辨識困難以及連續(xù)時間系統(tǒng)的時間導數(shù)問題;然后固定數(shù)據(jù)矩陣“R”規(guī)模,遞推更新數(shù)據(jù)矩陣獲得系統(tǒng)階數(shù),從而降低了辨識算法的計算負擔,解決了傳統(tǒng)的奇異值分解方法帶來的在線遞推困難.

本文其他部分構(gòu)成如下:第2節(jié)介紹問題描述;第3節(jié)介紹連續(xù)時間系統(tǒng)的遞推子空間辨識;第4節(jié)通過仿真驗證所提方法的有效性和精確性;第5節(jié)給出結(jié)論.

2 問題描述

考慮如下連續(xù)時間系統(tǒng):

其中:X(t)∈Rn為狀態(tài)向量,U(t)∈Rm為連續(xù)輸入向量,Y(t)∈Rl為連續(xù)輸出向量,A ∈Rn×n,B ∈Rn×m,C ∈Rl×n,D ∈Rl×m為系統(tǒng)矩陣,W(t)∈Rn,V(t)∈Rl分別為系統(tǒng)和觀測方程的連續(xù)高斯白噪聲,均值為零,協(xié)方差矩陣為

其中:E(·)是期望值算子,Q(t),S(t),R(t)是噪聲強度,δ(t ?τ)是δ函數(shù).U(t)與隨機噪聲過程W(t),V(t)都獨立.本文的目的是根據(jù)輸入輸出數(shù)據(jù)遞推更新系統(tǒng)的階數(shù)和系統(tǒng)矩陣A,B,C,D.在辨識連續(xù)時間系統(tǒng)問題之前,首先考慮以下確定性的情況:

其中輸入Ud(t)是足夠次可微的.考慮Yd(t)在每次不同時刻t1,t2,···,tN(不必是等距的)的直到i ?1階高階導數(shù),有以下輸入輸出矩陣關(guān)系:

通過對輸入輸出投影矩陣的奇異值分解可以得到矩陣Γi的列空間以及系統(tǒng)階數(shù),進而得到系統(tǒng)矩陣和噪聲強度[7].辨識過程中,在增加算法數(shù)值魯棒性的同時,也增加了子空間辨識的在線遞推困難.辨識算法要求構(gòu)建輸入輸出數(shù)據(jù)關(guān)系的充分條件是輸入U(t)和輸出Y(t)是足夠次可微的,但是在討論連續(xù)時間系統(tǒng)的情況下,連續(xù)高斯白噪聲W(t),V(t)處處不可微,而且狀態(tài)向量X(t)、輸出向量Y(t)處處不可導.所以在隨機噪聲存在的情況下,如何設(shè)計辨識算法遞推更新連續(xù)時間系統(tǒng)是一個難點問題.本文針對以上問題深入研究基于隨機分布理論的遞推子空間辨識方法(recursive subspace identification based on random distribution theory,RDT–RSID).

3 連續(xù)時間系統(tǒng)的遞推子空間辨識

本節(jié)利用測試函數(shù)和廣義分布函數(shù)理論推導出連續(xù)時間系統(tǒng)的輸入輸出矩陣方程,并采用將輸入輸出數(shù)據(jù)矩陣“R”規(guī)模固定的方法更新連續(xù)時間系統(tǒng).具體地,首先,對連續(xù)隨機分布函數(shù)進行求導,然后,得到連續(xù)時間系統(tǒng)的輸入輸出矩陣方程,最后遞推更新連續(xù)時間系統(tǒng).

3.1 連續(xù)隨機分布函數(shù)的微分計算

定義(?,F,P)為概率空間,空間D是在R上的C∞-函數(shù)φ(t)所構(gòu)成的空間,函數(shù)φ(t)為測試函數(shù).一個連續(xù)隨機過程{y(t,ω),?∞

如果對于任意函數(shù)φ(t)∈D,隨機變量y(φ)是高斯的,即隨機分布y(φ)稱為高斯過程,亦為規(guī)則廣義分布函數(shù).

對于?φ ∈D,隨機過程y(t)在分布意義下的一階二階導數(shù)可以計算為

隨機過程y(t)的k階導數(shù)為

測試函數(shù)(φ ∈S)可以從更廣泛的一類函數(shù)來選擇,要求函數(shù)有足夠階的導數(shù)并且在區(qū)間(?∞,∞)上為快速遞減函數(shù)[17–18].

定義{tj,j1,2,···,N}(?∞

隨機過程y(t)在分布意義下的一階、二階導數(shù)為

隨機過程y(t)的k階導數(shù)為

3.2 輸入輸出矩陣方程的建立

定義{tj,j1,2,···,N}(?∞

并且X(φ)(tj),U(φ)(tj),V(φ)(tj)同樣定義.

對方程(1)關(guān)于t進行微分,然后對兩邊結(jié)果同時乘φ(t;tj),并且在(?∞,∞)上進行積分,所以方程(1)在分布意義下的一階導數(shù)為

代入方程(1),并結(jié)合方程(7),可以得到

結(jié)合方程(8),二階導數(shù)以相同方式給出

結(jié)合方程(9),重復求導(i ?1)次,則第(i ?1)階導數(shù)為

通過列向量來表示方程(18),可以得到輸入輸出矩陣方程

U(φ)(tj),W(φ)(tj)和V(φ)(tj)與Y(φ)(tj)相似定義,Γi和Hi與方程(4)相似定義,Σi定義如下:

因此,通過輸出列向量表示的輸入輸出矩陣方程為

其中φ ∈D,狀態(tài)矩陣為

輸出矩陣為

并且Ui(φ),Wi(φ),Vi(φ)和Yi(φ)相似定義[19].

3.3 系統(tǒng)矩陣和噪聲強度的遞推估計

根據(jù)新的輸入輸出數(shù)據(jù){ut+1,yt+1},構(gòu)造數(shù)據(jù)向量yi(t+1)和ui(t+1)作為更新向量:

所以,有以下關(guān)系:

其中:ui(t ?j+1)和yi(t ?j+1)與yi(t+1)定義類似,u2i(t+1)和u2i(t ?j+1)是輔助變量,由輸入輸出數(shù)據(jù)構(gòu)建,

Yi(·|·),Ui(·|·)是塊Hankel矩陣,U2i(·|·)是由輸入輸出數(shù)據(jù)組成的輔助變量矩陣.

對于式(24)左邊的子矩陣在時刻t處進行LQ分解,有

對于式(24)右邊的子矩陣在時刻t+1處進行類似的LQ分解,有以下關(guān)系:

根據(jù)式(32)–(36),得到

通過式(39)以及對L32(t)(t)執(zhí)行奇異值分解,并依據(jù)文獻[14,18]辨識過程中最小二乘法和殘差分析法可得系統(tǒng)階數(shù)n以及系統(tǒng)矩陣.通過不斷添加新數(shù)據(jù)來更新矩陣L32(t)(t),從而遞推更新模型(1)中系統(tǒng)矩陣A,B,C,D和系統(tǒng)階數(shù)n.至此完成了連續(xù)時間系統(tǒng)(1)的遞推辨識過程.

本文基于隨機分布理論的遞推子空間辨識(RDT–RSID)方法可歸納如下:

步驟1依據(jù)式(11)–(13),得到輸入信號的一階至k階導數(shù);

步驟2通過式(21),得到系統(tǒng)的輸入輸出矩陣方程;

步驟3根據(jù)式(24),構(gòu)造數(shù)據(jù)矩陣等式關(guān)系;

步驟4由式(26),對新得到的輸入輸出向量執(zhí)行LQ分解得到{Lj(t),j0,1,···,4};

步驟5由式(39),計算L31(t)執(zhí)行奇異值分解L32(t)遞推更新模型系統(tǒng)矩陣和系統(tǒng)階數(shù).

4 仿真研究

4.1 數(shù)值仿真

本節(jié)通過數(shù)值例子的仿真實驗驗證了提出方法RDT–RSID的有效性及比較優(yōu)勢.考慮如下連續(xù)時間系統(tǒng)[20]:

e(t)為零均值白噪聲,由MATLAB加入且噪聲水平為20 dB,輸入信號選擇偽隨機數(shù)二進制序列(pseudo-random binary sequence,PRBS)信號,采樣時間為Ts0.01 s,執(zhí)行100次蒙特卡羅仿真.

分別采用方法文獻[9]中的遞推子空間辨識(recursive subspace identification,RSID)與RDT–RSID對系統(tǒng)(40)進行辨識,分別得到RDT–RSID與真實系統(tǒng)的平均伯德圖以及RSID與真實系統(tǒng)的平均伯德圖如圖1–2所示.可見無論頻率高低,RDT–RSID比RSID得到的模型更接近真實模型(40)的伯德圖,說明所提方法RDT–RSID的辨識精度更高.

圖1 RDT–RSID與真實系統(tǒng)的平均伯德圖Fig.1 The averaged bode diagram of the RDT–RSID and true system

圖2 RSID與真實系統(tǒng)的平均伯德圖Fig.2 The averaged bode diagram of the RSID and true system

4.2 連續(xù)攪拌釜反應過程應用研究

本節(jié)采用連續(xù)攪拌釜加熱器系統(tǒng)(continuous stirred tank heater,CSTH)的辨識數(shù)據(jù)進行仿真研究,驗證提出方法的有效性和精確性.連續(xù)攪拌釜加熱器是加拿大阿爾伯塔大學設(shè)計的一個小型試驗裝置,其基本配置如圖3所示[21].攪拌釜中的冷水和熱水通過加熱管的蒸汽進行加熱,混合的液體通過一個長管子排出.假設(shè)攪拌釜里的冷熱水均勻混合,水的溫度和輸出流的水溫相同.有3個比例積分控制器來控制冷水液位、溫度和流量.冷水的液位和溫度是可以被控制的,所有信號都是統(tǒng)一量程,為(4～20)mA.

圖3 連續(xù)攪拌釜加熱器的配置圖Fig.3 The configuration sketch of CSTH

本實驗選擇液位環(huán)和溫度環(huán)的設(shè)定值作為輸入,攪拌釜的液位和溫度作為輸出.溫度和液位分別受到幅值為0.5 mA和1 mA的隨機二進制干擾,此干擾信號由MATLAB系統(tǒng)辨識工具箱中的idinput函數(shù)產(chǎn)生.采樣時間設(shè)定為Ts1 s,采樣數(shù)據(jù)長度N800.

采用方法文獻[13]中的基于Laguerre濾波器的遞歸子空間辨識(recursive subspace identification using Laguerre filters,RSILF)和所提方法RDT–RSID進行辨識,通過RDT–RSID模型和RSILF模型的預測輸出分別如圖4–5所示.從圖中可以看出,RDT–RSID預測數(shù)據(jù)和測量數(shù)據(jù)的軌跡具有高度的一致性.為了更方便直觀的觀察,給出了兩種方法所對應的實際測量值和預測輸出值之間的誤差如圖6–7所示.從圖中可以看出,RDT–RSID模型輸出值與測量值的誤差較小,表明所提出的遞推辨識方法具有較高的辨識精度.

圖4 RDT–RSID模型的輸出Fig.4 The outputs of the RDT–RSID

圖5 RSILF模型的輸出Fig.5 The outputs of the RSILF

圖6 RDT–RSID模型輸出值與測量值的誤差Fig.6 The errors between the outputs of the RSILF and measurements

圖7 RSILF模型輸出值與測量值的誤差Fig.7 The errors between the outputs of the RSILF and measurements

根據(jù)文獻[7],辨識模型的預測誤差可以通過以下方程進行計算.均方根誤差(root mean squared error,RMSE)和預測誤差(sum of squares of prediction error,SSPE)平方和分別如下:

基于方法RDT–RSID和基于RSILF的液位和溫度測量變量的均方根誤差(RMSE)和預測誤差平方和(SSPE)如表1所示.

表1 RDT–RSID 和RSILF 方法建模誤差的RMSE和SSPETable 1 The RMSE and SSPE of prediction error from RDT–RSID and RSILF

從表中結(jié)果可以看出,基于RDT–RSID的液位和溫度測量變量的RMSE和SSPE值都比基于RSILF方法的小,這表明提出方法RDT–RSID有更好的模型辨識性能.

5 結(jié)論

本文提出了基于隨機分布理論的遞推子空間辨識方法.利用隨機分布理論構(gòu)造輸入輸出數(shù)據(jù)Hankel矩陣,得到分布意義下的輸入輸出矩陣方程,從而解決了連續(xù)時間系統(tǒng)的時間導數(shù)計算問題.采用將輸入輸出數(shù)據(jù)矩陣“R”規(guī)模固定的方法,降低了辨識算法的計算量和存儲成本.通過仿真實驗驗證了提出方法的有效性和精確性.下一步計劃在本文基礎(chǔ)上進一步考慮輸入信號和噪聲的相關(guān)性對算法的影響,通過理論分析結(jié)合實驗仿真說明算法的辨識效果.