LPV模型的動態壓縮測量辨識算法

邱棚， 李鳴謙， 姚旭日， 翟光杰,*， 王雪艷

(1. 中國科學院國家空間科學中心， 北京 100190； 2. 中國科學院大學， 北京 100049；3. 北京信息科技大學機電工程學院， 北京 100192)

目前，被控對象的系統結構越來越復雜，而且系統的運行條件也不再是穩定的、單一的，而是一個動態的過程，這些過程都表現出很強的非線性特征。比如，多熱源的溫度控制系統、渦扇引擎控制、車輛橫向控制、高速飛行器控制等[1]。考慮到當非線性系統被描述為線性參變 (Linear Parametric Variation, LPV)模型時，很多針對線性系統的控制理論都可以直接推廣使用[2]，LPV模型在分析這類非線性系統中發揮了重要的作用。因而，LPV模型被大量應用于基于模型的控制算法中，如魯棒控制、自適應控制和最優控制等。

一般來說，一個系統的LPV模型可以通過雅可比線性化、狀態變換法以及函數替代法等理論推導的方法得到[3]。不過，在實際系統中就需要通過系統辨識的方法獲得模型參數。現有LPV模型辨識方法分為兩類：基于狀態空間模型[4]和基于輸入輸出模型[5]。針對輸入輸出模型的辨識方法有一般有兩類：全局非線性擬合的辨識算法以及局部滑動窗口辨識算法。第1類算法要求選定一組恰當的函數基底，再計算各基底的權值，如非線性最小二乘算法[6]。這類算法對參變函數的近似程度高，但是需要對系統的先驗。第2類算法屬于非參數化的方法，典型算法如基于支持向量機辨識算法[7]和基于貝葉斯的辨識算法[8]以及遞推最小二乘算法等。此類算法在參變函數近似為分段常數函數的情況下，可以準確地辨識得到系統模型。但是，如果當參變函數是線性分段或者類似正弦函數等非線性的情況下，再利用常數近似就會引入較大的函數近似誤差。另外，最小二乘算法求解的基礎是要求辨識數據個數大于未知數個數。因而，當模型規模較大時，由于模型參數大量增加，這就需要辨識數據也大量增加。一方面，對于采樣成本高或者系統響應時間長的情況，很難獲得大量的辨識數據。另一方面，LPV模型的參數是一直在變化的，因而增加測量數據就意味著需要擬合更復雜的變化情況。如果近似模型維持原狀，就會降低模型的近似程度，也就是模型的近似程度也限制了辨識數據的多少。在這種情況下，最小二乘算法很難同時提高計算精度，又減少因為數據量增多而帶來的近似誤差。但其實高階模型中的參數是存在大量零值的，而非零值是十分有限的。因此，在解決高階模型辨識問題上，利用NNG[9]和LASSO[10]等稀疏估計器是一個有效的辦法。本文基于壓縮測量辨識(Compress Measurement Identification，CMI)算法提出的動態壓縮測量辨識(Dynamic CMI，DCMI)算法是從兩個方面解決上述問題。其一，利用“勻速變化”和“非勻速變化”模型，提高對參變函數函數的近似精度，同時避免函數基底的選擇。其二，利用壓縮感知理論，通過減少所需的辨識數據個數，從而提高算法在辨識數據有限的情況下的辨識精度及參數規模。

1 壓縮感知

壓縮感知理論是在2006年由Donoho[11]、Candès和陶哲軒[12]提出的一個高效的采樣理論。近年來，該理論在信號處理領域中，已經取得了大量令人矚目的成績[13-14]。信號的稀疏性是壓縮感知的理論基礎。對于一個離散信號x，其稀疏度可以定義為信號非零項的個數。所以，信號的稀疏度越小，該信號的稀疏程度就越高，那么理論上所需的采樣數也就越少。而采樣值yi是由測量向量ri分別與信號x進行內積運算得到的，即yi=〈ri,x〉。考慮多次測量，則測量值y可以寫為矩陣形式，即y=Rx。其中，矩陣R被稱為測量矩陣。如果R是一個單位矩陣，那么壓縮感知采樣就退化為了傳統采樣。在壓縮感知理論中，測量值y的長度要遠小于原信號x。為了保證這個采樣過程是無損的，Candès和Romberg給出了有限等距性質(Restrict Isometric Property，RIP)[15]，即當測量矩陣R滿足RIP時，采樣過程就是無損的。S階有限等距常數δS定義為使得測量矩陣R任選幾列組成的子矩陣RT,T?{1,2,…,N}滿足式(1)的最小值：

(1)

式中：x為任意稀疏信號。當有效等距常數(Restrict Isometry Constants,RIC)滿足0<δS<1時，該測量矩陣認為是滿足RIP的[16]。根據壓縮感知理論可知，欠定方程y=Rx解空間中最稀疏的解就是原信號x。由于這是一個組合問題(NP-Hard)，因此該問題需要利用貪婪的思想進行求解。而文獻[17]中提出的正交匹配追蹤(Orthogonal Matching Pursuit,OMP)算法是其中的代表，該算法在每次迭代中尋找最匹配的向量，作為信號的支撐基，從而得到原信號的最佳近似。另外，Candès和陶哲軒[12]已經證明了，該問題可以轉換為在解空間中尋找信號的最小1范數的解，如下：

s.t.y=Rx

(2)

上述問題可以進一步轉化為典型的線性規劃問題，那么也就有很多成熟的算法可以使用，比如內點法、單純形法等。除此以外， Li和Zhang還提出了可以利用TV范數重建原信號[18]。

2 LPV模型

本文中采用LPV模型的輸入輸出的表示形式為

y(k)+a1(p)y(k-1)+…+ana(p)y(k-na)=

b1(p)u(k-d-1)+…+bnb(p)u(k-d-

nb)+e(k)

(3)

式中：y(k)為系統的輸出信號；u(k)為系統的輸入信號；參數p為一個可測量的參量；ai(p)和bi(p)為關于參數p的函數；na和nb分別為模型的輸入和輸出階數；當p=k時，LPV模型也就變化為常見的線性時變(Linear Time-Variant, LTV)模型。在實際系統中，還要考慮到系統的輸入延遲以及噪聲，在此用e(k)為一個零均值的隨機噪聲；d為系統時延。上述模型可以改寫為

y(k)=φT(k)θ(p)+e(k)

(4)

式中：φ(k)∈Rna+p由輸入輸出數據組成；θ(p)∈Rna+p為系統狀態函數，由多個參變函數組成；l為輸入信號的的長度，而且有d+nb

3 動態壓縮測量辨識算法

CMI算法是利用壓縮感知理論對線性時不變系統進行辨識的一種算法。該算法在應用于LPV模型時，則需要加以改進。本文受到壓縮感知運動成像理論的啟發[19]，將參數變化函數的連續變化視作多個離散點之前的轉換。假設系統狀態函數θ(p)是滿足唯一性的關于參數p的復雜函數，即對于每個參數p都有唯一的狀態θ(p)與之對應。這個假設在大部分非線性系統中都是滿足的，比如溫度控制系統中的傳熱系數是與溫度一一對應的。但是，當系統中包含繼電器特性環節，該系統就不滿足這個假設。對函數θ(p)進行采樣就可得到一系列離散的狀態點{θ(p1),θ(p2),…,θ(pn)}。為了保證采樣的準確性，采樣過程要滿足Nyquist采樣定理，即采樣頻率要大于θ(p)帶寬的二倍以上。此時，連續函數θ(p)被n個離散點分割為n-1段函數，原辨識問題也就變為了對每段函數的辨識問題。如果可以準確地辨識得到其中一段的系統狀態，那么也就可以通過滾動迭代的方法得到局部參變函數。以下將以辨識系統其中一段，即從離散狀態θ(p1)到θ(p2)的過程為基礎。

由于參變函數可能是非線性函數，所以在離散后的每一段函數中還是非線性的，也就很難直接對其進行辨識。為了簡化這個問題，可以引入函數θ(p)是連續n階可導的假設。此時，該函數可以利用泰勒級數在p1和p2的中點處展開，此時參變函數可以表示為

(5)

式中：θ(n)(p)為θ(p)的n階導數；Rn(p)為泰勒公式的余項。如果只保留一階導數，得到原函數的一個近似線性函數為

(6)

該線性函數和原函數的近似誤差主要取決于高階項的大小。如果采樣頻率遠大于原信號頻率或者原信號是緩慢變化的，那么利用式(6)來近似就可以得到比較好的結果。然而，式(6)等號右側的各項都是未知的，所以可以將式(6)改寫為

(7)

相應的，系統的輸出也可以視為由各離散狀態以及輸入信號而產生的。此時，如果使用傳統的辨識算法，辨識結果將無法收斂。為了保證近似模型的精度就要提高系統采樣頻率，但系統的參數一直在改變，所以對于每個狀態下所獲得的辨識數據M的個數可能小于待辨識的參數個數N。所以，該辨識問題相當于求解一個欠定的方程組。顯然，用傳統的辨識算法都很難處理這種情況。這就需要用到壓縮感知理論，結合該理論的CMI算法正好適用于這種采樣數不足的情況。不過，利用CMI算法對該模型進行辨識，仍然會存在很大的誤差。這是因為CMI算法適用于模型不變的系統，而此處系統狀態是一直在改變的。因此，還需要引入“運動”的思想，即系統狀態一直在變化的，是從一個離散狀態“運動”到另一個離散狀態的過程。而系統每一時刻的輸入是對這個“運動”狀態的一次觀測，系統的輸出就是狀態在“運動”過程中所產生的。也就是說，系統在任意時刻的輸出，都可以表示為由相鄰的兩個狀態與輸入信號的卷積所產生的。從之前的敘述可知，系統某一時刻的采樣值y(k)可以認為是由臨近的兩個離散狀態所共同影響，并且兩個狀態對系統輸出影響的比例是一直在變動的。此時，式(4)可以改寫為

y(k)=(1-λ)φT(k)θ(p1)+

λφT(k)θ(p2)+e(k)

(8)

考慮M次采樣，系統輸出可以表示為

y=(P*Φ)θ(p1)+(Q*Φ)θ(p2)+e=

[P*Φ,Q*Φ][θ(p1),θ(p2)]+e

(9)

式中：pi和qi分別為在k時刻狀態θ(p1)和θ(p2)在系統輸出中所占的權重；Φ代表測量矩陣；P、Q為權重矩陣；e為噪聲向量。其中，權重矩陣P、Q與Φ的“*”運算為Hadamard內積，即對應項相乘。而且，從本節描述可知矩陣P、Q與變化速度v相關。對于“勻速”變化，pi和qi的大小由變化速度v所決定；對于非勻速變化，只需要按照 “非勻速”速度計算得出各項值即可。在此，本文將[P*Φ,Q*Φ]稱為比例測量矩陣，用B表示。而[θ(p1),θ(p2)]整體變為了一個新的待辨識參量，其長度為原信號的兩倍，用x表示。那么，式(9)就變為了常見的壓縮感知的形式，y=Bx+e。通過上述方法，就可以在采集到所有辨識數據后獲得系統模型，即離線的對系統進行辨識。而此方法是可以很容易推廣到在線辨識的情況的。因為，本文對假設的各離散狀態的位置是任意的。所以，只要在采集到了最新的M個辨識數據后，就可以辨識得到兩個離散狀態。而且，通過滾動辨識的方法就可以實時的獲取最新的系統狀態，利用該狀態得到的系統輸出預測值也是最準確的。在此給出，系統的一步預測器：

y*(k+1)=φ(k+1)θ*(k)

(10)

式中：y*(k+1)為一步預測器在k時刻的預測值；φ(k+1)為在k+1時刻的系統輸入輸出組成的觀測向量；θ*(k)為辨識算法得到第k時刻的模型參數。

為了保證采樣的準確性，就需要測量矩陣滿足RIP。在壓縮感知理論中，Candès和陶哲軒[12]證明了獨立同分布的隨機測量矩陣是滿足RIP的。不過，由于線性系統的輸出是用卷積的形式表示的，所以測量矩陣是一個有結構的Toeplitz矩陣。Bajwa等證明了當Toeplitz矩陣各項服從高斯分布或者伯努利分布時，Toeplitz矩陣有很大可能滿足RIP，只要測量數M>O(S2lb(N/S))[20]，N為信號總長度。隨后Rauhut等利用Dudley不等式證明了，當測量數M>O(S1.5·lb(N1.5))[21]，測量矩陣也是滿足RIP的。不過，隨機測量矩陣對測量數的要求只要M>O(Slb(N/S))[22]，這是因為Toeplitz矩陣各列之間的相關度要大于隨機測量矩陣。對于本文中的比例測量矩陣，該矩陣是由兩個Toeplitz矩陣分別乘以不同的系數所組成，因而該矩陣相比于測量矩陣，其各列之間的相關性被進一步放大。因此，要想獲得同等準確的采樣值，所需的測量數肯定要更多。

4 仿真試驗

假設待辨識系統模型參數是稀疏的，即在{ai(p)}和{bi(p)}中一共只有S個位置是非零值，其他位置均為0。鑒于實際系統中S個非零項的大小和位置都是未知的，因此在每一次仿真中非零項的位置和大小在都會隨機改變。由于DCMI方法并不要求提前設定準確的模型階數，因而在仿真中選擇的模型階數足夠大即可。仿真中設置稀疏度S=6，輸入信號階數na與輸出信號階數nb均為20，即信號總長度N=40。

系統的輸入信號u(k)各項是從標準的正態分布中隨機選取。為了接近真實情況，觀測值是包含噪聲的。該噪聲的大小由信噪比(SNR)所衡量。本文比較了DCMI算法、CMI算法以及最小二乘算法對θ(p1)和θ(p2)的變化過程的辨識結果。此外，由于最小二乘算法在欠定條件下無法獲得唯一解，因而本文在目標函數中增加了正則項來解決。其中，DCMI和CMI的恢復算法均選用OMP算法來重建模型參數信號。為了保證仿真結果的普適性，本文進行了200次蒙特卡羅仿真，并且在每次仿真中都會重新隨機生成一個新的LPV系統模型。考慮到真實系統中的模型真值是未知的，因而本文除了比較模型本身外，還比較了利用模型預測的系統輸出以及實測系統輸出的差距來評價模型的準確性。具體指標為均方誤差(Mean Square Error，MSE)和平均接近程度(Average Fit Rate，AFR)，該指標的具體形式如下：

(11)

4.1 勻速變化模型

假設兩個系統狀態θ(p1)和θ(p2)之間的變化是線性的，LPV模型真值和不同算法辨識結果的對比如圖1所示。雖然已知狀態變化是線性，但是狀態中的各參數的變化速度是不同的，而且其大小也是是未知的。

此時，權重矩陣P和Q由系統測量值個數M給出：

(12)

圖1中：Idx表示θ(p)向量中的位置，f(p)表示θ(p)中下標為Idx處的參數值大小。該圖中展現了200次試驗中的兩次結果，其中圖1(a1)、(b1)為LPV模型θ(p)的真值，其他各分圖分別表示各辨識算法計算得到的模型參數值。可以看出，最小二乘算法和CMI算法的結果是用常數來近似表示線性函數，只有DCMI算法辨識結果與模型真值是一致的。另外，最小二乘算法的結果在很多原本是零的位置，得到了大量的非零值，即導致了模型的過擬合，引入了大量的誤差。

鑒于在高階LPV模型的辨識問題中，系統輸入輸出的測量值個數是主要的影響因素，因而隨后重點分析測量數對各算法的影響。由于在實際系統中，模型真值一般是未知的，因而也就無法直接對模型進行比較。另一個評價的辦法就是給定一個輸入信號，對比系統輸出的實測值和利用模型得到的預測值，該預測值由第3節中提到的一步預測器給出。圖2顯示了3種算法在不同測量數條件下的預測結果。其中，最小二乘算法效果最差，而且隨著測量數的增加，預測精度反而下降。此外，該算法在測量數M=40時，由于過擬合導致極大的預測誤差。相比較而言，CMI算法要強于最小二乘算法，但同樣隨著測量數的增加，算法的預測精度反而逐漸下降。這同樣是由于算法沒有考慮到參數是動態變化的所導致的。從圖中可以看到，DCMI算法的預測精度顯然是最高的，但該算法在測量數不足時表現較差。因為，該算法需要同時重建兩個時刻的模型參數，其長度是其他兩個算法的2倍。由于觀測本身存在誤差，而且一步預測器利用的是上一時刻的模型參數，所以DCMI的預測值與觀測值還是存在差距的。同時，由于算法計算時間同測量數多少直接相關。因而，在仿真中記錄了各算法所花費的時間，并發現算法計算花費時間均在在1 ms以內，但會隨著模型規模的增加而增加。

圖1 LPV模型真值和不同算法辨識結果對比Fig.1 Comparison between truth value of LPV model and identification result of different algorithms

圖2 測量數對不同算法預測精度的影響Fig.2 Influence of measurement number on prediction accuracy of different algorithms

接下來，本文對3種算法的魯棒性進行了測試。為了避免測量數在結果中的影響，在此設置測量數M=50，仿真結果如圖3所示。從圖中可以看出，CMI和DCMI兩種算法對白噪聲都有一定的抗干擾能力，而最小二乘算法的對噪聲比較敏感，其結果偏差較大。當噪聲分貝降低到40 dB左右時，噪聲基本不會影響DCMI算法的預測的準確性。這是因為CMI算法考慮到了信號本身是稀疏的，所以大部分噪聲在辨識結果中被抑制了。而且重建算法OMP采用了貪婪的思想，即只考慮了測量矩陣中影響最大的幾列，從而降低了噪聲對恢復結果的干擾。

圖3 觀測噪聲對不同算法預測精度的影響Fig.3 Influence of measurement noise on prediction accuracy of different algorithms

4.2 非勻速變化模型

在4.1節中考慮了參數“勻速”變化的情況，本節將以一個參數以“勻加速”變化的系統作為例子來說明該方法在參數非線性變化的LPV系統中的使用。先假設系統兩個相鄰的離散狀態θ(p1)和θ(p2)，參數變化速度v以加速度a均勻增加。由于參數變化是非勻速的，所以權重矩陣也就需要根據測量數的不同而重新計算。已知，兩個離散狀態之間的任意狀態可以由式(7)表示。那么，對于測量初始時刻ts的狀態θs和最后時刻te的狀態θe有：

(13)

式中：λe和λs分別為初始時刻和最后時刻兩個狀態的影響比例。

將式(13)整理可得

(14)

再將式(14)代入式(7)，可得θs和θe之間任意狀態θ的表示為

(15)

根據式(15)就可以計算得出DCMI算法所需的權重矩陣P和Q。由于4.1節已經分析了噪聲的影響，所以為了減少噪聲對結果的影響，設置噪聲大小為40 dB。將圖4與圖2對比可以發現，在線性變化的系統中最小二乘算法和CMI算法的AFR與DCMI算法差距在10%左右，而在非線性變化的系統中的該差距被拉大到20%以上。這是因為最小二乘和CMI算法的本質都是靜態算法，因而算法辨識出的模型參數是各狀態的均值。而該參數均值與當前值的差距，在非線性變化過程中被拉大了，所以利用該均值得到的預測值與觀測值的差距也就被拉大了。這同樣影響到了DCMI算法的準確性，相比于線性變化情況，此處最佳的AFR要降低了2%左右。而且，可以看出此時算法所需的測量數要高于線性變化的情況。其原因可能是由于非線性的權重矩陣使得測量矩陣性質變差。

圖4 線性變化情況下測量數對不同算法預測精度的影響Fig.4 Influence of measurement number on output prediction accuracy of different algorithms in case of linear variation

5 結 論

1) 根據仿真試驗結果表明，DCMI算法中使用的勻速變化模型可以無損地對線性參變函數近似，而“非勻速變化”模型對非線性參變函數近似程度明顯好于傳統算法。

2) 本文提出的DCMI算法在測量數據量不足的情況下，仍然能夠得到準確的LPV模型。另外，DCMI算法可以很好地避免參數過擬合問題，而無需根據先驗知識選定恰當的模型階數。

3) DCMI算法可以有效地抵抗白噪聲對辨識結果的干擾。