基于Laguerre濾波器的核范數子空間辨識

于 淼 劉建昌 王洪海 張文樂

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

1 引言

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

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

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

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

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

2 問題描述

考慮如下連續時間隨機系統[18]:

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

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

根據式(1)，得到代數關系如下:

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

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

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

3 輸入輸出代數方程的建立

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

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

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

其中符號w表示算子.

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

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

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

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

至此獲得一組Laguerre濾波器，它的結構如圖1所示.

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

根據式(5)和式(9)，有以下關系:

根據式(5)，系統(1)可以轉換為以下形式狀態空間模型:

其中:初始狀態X0=0，

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

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

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

其協方差矩陣為

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

其中:

4 核范數子空間辨識

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

其中y是最優變量.

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

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

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

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

根據式(22)，有

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

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

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

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

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

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

進一步，

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

系統矩陣可以通過以下方程求解:

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

噪聲強度為

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

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

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

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

步驟3根據式(20)，構造秩最小化問題;

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

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

5 仿真研究

本節通過數值例子的仿真實驗驗證了提出NLFSID方法的有效性、精確性及比較優勢.考慮如下兩輸入兩輸出連續時間隨機系統[8]:

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

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

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

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

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

圖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

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

圖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 結論

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