于 淼 劉建昌 王洪海 張文樂

(1.東北大學(xué)秦皇島分校控制工程學(xué)院，河北秦皇島 066004;2.東北大學(xué)信息科學(xué)與工程學(xué)院，遼寧沈陽 110819;3.東北大學(xué)流程工業(yè)綜合自動化國家重點實驗室，遼寧沈陽 110819;4.華東理工大學(xué)信息科學(xué)與工程學(xué)院，上海 200237)

1 引言

近年來，子空間辨識作為系統(tǒng)辨識的重要分支已引起控制領(lǐng)域?qū)＜覍W(xué)者的濃厚興趣[1-5].其中大多數(shù)方法都是解決離散時間系統(tǒng)的辨識問題，并依據(jù)所得模型進行數(shù)字控制器的設(shè)計.然而，與離散時間系統(tǒng)的辨識相比較，連續(xù)時間系統(tǒng)的辨識在參數(shù)估計，H∞控制和自校正控制等方面應(yīng)用更加廣泛[6-7].此外，多數(shù)物理現(xiàn)象都具有連續(xù)屬性，描述它們的數(shù)學(xué)模型是微分方程.因此，針對連續(xù)時間系統(tǒng)辨識問題提出相應(yīng)的解決方案具有重要的理論意義和實際價值.

辨識連續(xù)時間系統(tǒng)的主要方法可歸類為直接法和間接法.間接法是依據(jù)輸入輸出數(shù)據(jù)估計離散時間系統(tǒng)，然后將其轉(zhuǎn)化為連續(xù)時間系統(tǒng).然而，這種系統(tǒng)轉(zhuǎn)換過程可能存在采樣時間難以選擇，矩陣對數(shù)的運算復(fù)雜，離散時間系統(tǒng)的零點到連續(xù)時間系統(tǒng)的極點難以轉(zhuǎn)換等問題[8-9].直接法是直接依據(jù)輸入輸出數(shù)據(jù)估計連續(xù)時間系統(tǒng).然而，直接法辨識系統(tǒng)的一個難點問題是如何正確處理輸入輸出變量的時間導(dǎo)數(shù).目前有很多方法克服了重構(gòu)時間導(dǎo)數(shù)的困難.文獻[9]通過Laguerre濾波器處理輸入輸出信號，研究了基于Laguerre濾波器的子空間辨識方法.文獻[7]采用泊松矩函數(shù)(Poisson moment functionals，PMF)方法解決了輸入輸出時間導(dǎo)數(shù)和積分的估計問題，并與子空間辨識方法相結(jié)合辨識了連續(xù)時間模型.在文獻[10]中，作者利用隨機分布理論得到了分布意義下的輸入輸出信號的時間導(dǎo)數(shù)，從而推導(dǎo)出相應(yīng)的輸入輸出代數(shù)方程，并利用提出的新方法得到了連續(xù)時間隨機系統(tǒng).在文獻[11]中，作者通過Laguerre濾波和Laguerre投影處理輸入輸出信號，提出了基于Laguerre濾波的子空間辨識(subspace identification via Laguerre filters，LFSID)方法.

雖然上述子空間辨識方法解決了連續(xù)時間系統(tǒng)的時間導(dǎo)數(shù)求解問題，但是都使用了奇異值分解去估計系統(tǒng)階數(shù)，這種處理方式帶來了新的問題，即當存在過程噪聲和測量噪聲時，如何獲取最優(yōu)的系統(tǒng)階數(shù)[12].這種構(gòu)建低秩矩陣逼近的秩最小化問題是一個非確定性多項式困難(non-deterministic polynomial-hard，NP-hard)問題[13]，核范數(shù)最小化方法是解決此問題的有力工具.由于具有逼近矩陣的秩最小化和相應(yīng)的線性性質(zhì)，核范數(shù)最小化方法已經(jīng)成功地應(yīng)用于系統(tǒng)辨識的研究之中[14].文獻[15]通過半正定規(guī)劃描述核范數(shù)最小化問題，利用內(nèi)點法求解此優(yōu)化問題，并且將提出的方法應(yīng)用于線性系統(tǒng)的辨識中.在文獻[16]中，作者研究了基于加權(quán)核范數(shù)最小化的子空間辨識方法，并采用交替方向乘子法求解此優(yōu)化問題.文獻[17]基于文獻[18]研究了核范數(shù)系統(tǒng)辨識方法，并且將其應(yīng)用于缺失輸入輸出數(shù)據(jù)的隨機系統(tǒng)中.在文獻[12]中，作者提出了頻域的核范數(shù)最小化子空間辨識方法，并且研究了模型階數(shù)和殘差的Pareto最優(yōu)軌跡問題.文獻[3]探索了基于Pareto優(yōu)化的核范數(shù)子空間辨識方法，并建立了新息形式的多變量狀態(tài)空間模型.

可見，核范數(shù)最小化方法能夠有效處理NP-hard問題，從而獲得低秩矩陣逼近.目前核范數(shù)最小化方法對于帶有過程噪聲和測量噪聲的連續(xù)時間隨機系統(tǒng)的辨識涉及較少.在實際的生產(chǎn)過程中，過程噪聲和測量噪聲的存在給連續(xù)時間系統(tǒng)的辨識帶來了一定的困難，如難以確定系統(tǒng)的階數(shù)等.因此，當系統(tǒng)階數(shù)未知時，應(yīng)用核范數(shù)最小化方法去解決帶有過程噪聲和測量噪聲的連續(xù)時間隨機系統(tǒng)的辨識問題是一個挑戰(zhàn).

本文對基于Laguerre濾波器的核范數(shù)子空間辨識方法進行深入研究，有效地解決了帶有過程噪聲和測量噪聲的連續(xù)時間隨機系統(tǒng)的辨識問題.本文其他部分構(gòu)成如下:2)問題描述;3)采用Laguerre濾波器獲得系統(tǒng)的輸入輸出代數(shù)方程;4)通過核范數(shù)子空間辨識方法確定系統(tǒng)的最優(yōu)階數(shù)，并且分別得到模型的系統(tǒng)矩陣和噪聲強度;5)通過數(shù)值仿真驗證所提方法的有效性和精確性;6)給出結(jié)論.

2 問題描述

考慮如下連續(xù)時間隨機系統(tǒng)[18]:

其中:X(t)∈Rn是狀態(tài)向量;Y(t)∈Rl是輸出向量;A ∈Rn×m，C ∈Rl×n是系統(tǒng)矩陣;W(t)∈Rn，V(t)∈Rl是高斯白噪聲，其協(xié)方差矩陣為

其中:E(·)是期望算子，Δ(t-τ)是Dirac delta函數(shù).本文的目的是根據(jù)輸入輸出數(shù)據(jù)確定系統(tǒng)階數(shù)和系統(tǒng)矩陣A，C.

根據(jù)式(1)，得到代數(shù)關(guān)系如下:

矩陣Γi∈Rli×n(i ＞n)是廣義能觀性矩陣，矩陣Hi∈Rli×mi是低階分塊三角矩陣.它們定義為[19]

一方面，系統(tǒng)矩陣A，C可以通過矩陣Γi和矩陣Hi導(dǎo)出.對于式(4)，Yi存在至少(i-1)階導(dǎo)數(shù)是必要條件.然而對于隨機情況，噪聲W(t)和V(t)是高斯白噪聲過程，處處不可微.并且狀態(tài)X(t)和輸出Y(t)也不存在時間導(dǎo)數(shù).另一方面，通過對輸入輸出投影矩陣的奇異值分解可以得到矩陣Γi的列空間.理論上，奇異值的個數(shù)等于系統(tǒng)的階數(shù)[20].然而，由于噪聲的影響，在構(gòu)造低秩矩陣逼近的過程中如何獲取最優(yōu)的系統(tǒng)階數(shù)是一個難點問題.

針對上述帶有過程和測量噪聲的連續(xù)時間隨機系統(tǒng)的辨識問題，本文提出了基于Laguerre濾波器的核范數(shù)子空間辨識(nuclear norm subspace identification based on Laguerre filters，NLFSID)方法.具體地，推導(dǎo)了系統(tǒng)的輸入輸出代數(shù)方程，詳見第3節(jié);解決了獲得系統(tǒng)最優(yōu)階數(shù)的問題，詳見第4節(jié).

3 輸入輸出代數(shù)方程的建立

本小節(jié)通過Laguerre濾波器解決輸入輸出時間導(dǎo)數(shù)難以求解的問題.為了獲得輸入輸出代數(shù)方程，給出以下定義.

定義L2(R+)表示平方可積空間和時間區(qū)間0 ＜t＜∞的Lebesgue可測函數(shù)，內(nèi)積定義為

空間H2是L2(R+)的閉子空間[9].一階全通濾波器為

其中符號w表示算子.

定義li(t)是i階全通濾波器的時域表示.[liy](t)是y(t)和li(t)的卷積，表示為

一組Laguerre濾波器可以通過全通濾波器得到.零階Laguerre濾波器為

式(7)乘式(5)，得到一階Laguerre濾波器

以此類推，第i階Laguerre濾波器為

至此獲得一組Laguerre濾波器，它的結(jié)構(gòu)如圖1所示.

并且L0(s)至L4(s)Laguerre濾波器的脈沖響應(yīng)函數(shù)如圖2所示.

根據(jù)式(5)和式(9)，有以下關(guān)系:

根據(jù)式(5)，系統(tǒng)(1)可以轉(zhuǎn)換為以下形式狀態(tài)空間模型:

其中:初始狀態(tài)X0=0，

圖1 一組Laguerre濾波器的框圖Fig.1 Block diagram of a Laguerre filter bank

圖2 L0(s)至L4(s)Laguerre濾波器的脈沖響應(yīng)函數(shù)Fig.2 Impulse reponse functions of the Laguerre filters L0(s)up to L4(s)

新的噪聲過程Ww和Vw表示為

其協(xié)方差矩陣為

根據(jù)式(10)-(11)，得到輸入輸出代數(shù)方程

其中:

4 核范數(shù)子空間辨識

通過時刻{tk，k=1，2，···，N}(-∞＜tk＜∞)的輸出數(shù)據(jù)，得到矩陣Yw:

其中y是最優(yōu)變量.

式(20)是一個NP-hard問題[21]，而核范數(shù)最小化方法可以解決該問題，基于此上式轉(zhuǎn)化為

其中:‖·‖*是核范數(shù)，Φw(y)是線性函數(shù)，β是加權(quán)參數(shù)，ym是通過Laguerre濾波器得到的測量輸出序列.

優(yōu)化問題(21)可以通過交替方向乘子法求解.定義新變量Zw=-Φw(y)，則問題(21)的增廣拉格朗日函數(shù)為

由Zw-最小化，y-最小化和Δ-更新構(gòu)成交替方向乘子法.

根據(jù)式(22)，有

根據(jù)文獻[22]的迭代算法及式(23)，可得

其中:U，V 可以通過奇異值分解得到;

其中:rp是原始殘差，rd是對偶殘差，ep是原始承受量，ed是對偶承受量，ea是絕對誤差，er是相對誤差，n是測量輸出數(shù)據(jù)長度，矩陣是Zw之前迭代數(shù)值.

采用以下更新方程代替采用固定懲罰參數(shù)ρ，改進交替方向乘子法的收斂速度:

其中μ＞1和τ ＞1是參數(shù)(一般μ=10，τ=2).

從方程(24)的最優(yōu)解Φw(y)，計算奇異值分解:

進一步，

根據(jù)式(25)及文獻[6]，卡爾曼狀態(tài)序列為

系統(tǒng)矩陣可以通過以下方程求解:

其中ρω和ρυ是Kalman濾波殘差.

噪聲強度為

根據(jù)式(28)的系統(tǒng)矩陣和式(12)，得到系統(tǒng)矩陣A，C，由式(25)估計系統(tǒng)(1)階數(shù).至此完成了連續(xù)時間隨機系統(tǒng)的辨識過程.

本文基于Laguerre濾波器的核范數(shù)子空間辨識(NLFSID)方法可歸納如下:

步驟1依據(jù)式(9)，得到通過Laguerre濾波器的輸入信號;

步驟2通過式(10)-(11)，得到系統(tǒng)的輸入輸出代數(shù)方程;

步驟3根據(jù)式(20)，構(gòu)造秩最小化問題;

步驟4由式(22)-(24)，采用交替方向乘子法對優(yōu)化問題進行求解;

步驟5由式(25)-(29)，獲得系統(tǒng)階數(shù)以及連續(xù)時間隨機系統(tǒng)的模型.

5 仿真研究

本節(jié)通過數(shù)值例子的仿真實驗驗證了提出NLFSID方法的有效性、精確性及比較優(yōu)勢.考慮如下兩輸入兩輸出連續(xù)時間隨機系統(tǒng)[8]:

e(t)為零均值單位方差的高斯白噪聲過程，輸出信號為根據(jù)系統(tǒng)(30)得到的仿真輸出，N=1000，輸入輸出信號連續(xù)采樣時間為Ts=0.01 s，執(zhí)行100次蒙特卡羅仿真.

分別采用NLFSID 與LFSID 方法對系統(tǒng)(30)進行辨識，得到的規(guī)范化奇異值指標的平均值如圖3所示，圖中給出了分別由兩種方法辨識得到的式(19)中矩陣Φw的規(guī)范化奇異值平均值.

根據(jù)最大奇異值的個數(shù)確定系統(tǒng)階數(shù).可以看出，LFSID的奇異值很接近，這使得模型的階數(shù)確定非常困難.相比之下，NLFSID的奇異值具有顯著性差異，從而確定模型階數(shù)更容易.

分別采用NLFSID與LFSID方法對系統(tǒng)(30)進行辨識，得到相應(yīng)的方法辨識模型的平均伯德圖如圖4-5所示.

可見無論頻率高低，NLFSID比LFSID得到的模型更接近真實模型的伯德圖曲線(30)、系統(tǒng)辨識的精度更高.

圖3 NLFSID與LFSID的最大奇異值指標平均值Fig.3 The averaged index of largest singular values index of the NLFSID and LFSID

圖4 NLFSID模型的平均伯德圖Fig.4 Averaged bode plot of NLFSID model

采用兩種方法獲得系統(tǒng)(30)的模型后，得到相應(yīng)的NLFSID與LFSID的零極點平均值如圖6-7所示.從圖形中可見NLFSID比LFSID得到的結(jié)果更接近真實模型的零極點、辨識模型更精確.

圖5 LFSID模型的平均伯德圖Fig.5 Averaged bode plot of LFSID model

圖6 NLFSID模型的零極點平均值Fig.6 Pole-zero average of NLFSID model

圖7 LFSID模型的零極點平均值Fig.7 Pole-zero average of LFSID model

6 結(jié)論

本文提出了一種基于Laguerre濾波器的連續(xù)時間隨機系統(tǒng)的核范數(shù)子空間辨識方法.利用Laguerre濾波器對數(shù)據(jù)進行處理，避免了輸入輸出Hankel矩陣的時間導(dǎo)數(shù)計算.通過核范數(shù)最小化方法獲得最優(yōu)的系統(tǒng)階數(shù)，并采用交替方向乘子法解決此優(yōu)化問題.最后，通過數(shù)值例子的比較仿真實驗驗證了提出方法的有效性和精確性，并說明了本方法更容易獲得模型階數(shù)、辨識精度更高的優(yōu)勢.