故障條件下子空間預(yù)測(cè)控制的對(duì)偶分解

王建宏， 許 鶯， 熊朝華， 徐 波

1.中國(guó)電子科技集團(tuán)公司第二十八研究所，南京210007

2.南京大學(xué)天文與空間科學(xué)院，南京210093

子空間辨識(shí)是通過最基本的矩陣奇異值分解直接將過去和將來輸入-輸出觀測(cè)數(shù)據(jù)序列構(gòu)成的矩陣進(jìn)行分解，以得到狀態(tài)空間方程的各個(gè)矩陣式.子空間預(yù)測(cè)控制來源于子空間辨識(shí)，直接從輸入-輸出觀測(cè)數(shù)據(jù)序列中構(gòu)造將來時(shí)刻的輸出預(yù)測(cè)值，而不必根據(jù)估計(jì)的狀態(tài)序列重構(gòu)狀態(tài)空間矩陣式.

子空間預(yù)測(cè)控制是一種基于數(shù)據(jù)驅(qū)動(dòng)的控制方法，能結(jié)合系統(tǒng)辨識(shí)和預(yù)測(cè)控制的優(yōu)點(diǎn)，直接由觀測(cè)數(shù)據(jù)構(gòu)造預(yù)測(cè)控制中所需的輸出預(yù)測(cè)估計(jì)值.文獻(xiàn)[1]構(gòu)造最簡(jiǎn)單子空間預(yù)測(cè)控制中的輸出預(yù)測(cè)值，文獻(xiàn)[2]比較該輸出預(yù)測(cè)值與迭代相關(guān)校正控制下的輸出預(yù)測(cè)值.文獻(xiàn)[3]采用經(jīng)典的橢球優(yōu)化算法求解無約束子空間預(yù)測(cè)控制中所出現(xiàn)的目標(biāo)準(zhǔn)則函數(shù)的全局最優(yōu)值，并利用凸分析理論中的水平集構(gòu)造該優(yōu)化算法中的初始橢球.文獻(xiàn)[4]將子空間預(yù)測(cè)控制應(yīng)用于主動(dòng)噪聲的振動(dòng)抑制，有效地解決了小型直升機(jī)懸停狀態(tài)時(shí)的顫振抑l制問題.文獻(xiàn)[5]在施加故障的上下限約束條件下，采用快速梯度算法估計(jì)故障和設(shè)計(jì)控制器.

本文對(duì)于故障條件下的狀態(tài)空間方程，根據(jù)輸入-輸出觀測(cè)數(shù)據(jù)序列建立故障與殘差間的函數(shù)關(guān)系，并以此建立一個(gè)最小二乘優(yōu)化問題來求解故障估計(jì)值.在輸出預(yù)測(cè)估計(jì)值的推導(dǎo)形式下分析殘差的統(tǒng)計(jì)分布性能，并利用矩陣直積的性質(zhì)推導(dǎo)出殘差矢量在瞬時(shí)時(shí)刻處的結(jié)構(gòu).將輸出預(yù)測(cè)估計(jì)代入預(yù)測(cè)控制的二次目標(biāo)準(zhǔn)則函數(shù)中，考慮一個(gè)帶有等式和不等式約束的規(guī)范化型的二次規(guī)劃問題.對(duì)于較復(fù)雜的二次規(guī)劃問題，提出對(duì)偶分解策略，將帶有等式和不等式約束的二次規(guī)劃問題轉(zhuǎn)化為對(duì)偶形式下的無約束優(yōu)化問題，再用迭代梯度算法進(jìn)行求解.

1 問題描述

考慮如下的隨機(jī)離散時(shí)間狀態(tài)空間模型：

式中，x(k)∈Rn,y(k)∈Rl,u(k)∈Rm,f(k)∈Rnf分別為狀態(tài)、輸出、控制和故障.矩陣A、B、C、D、E、G、K為恰當(dāng)維數(shù)的實(shí)有界矩陣，其中新息e(k)定義為

將新息e(k)的定義式(2)代入狀態(tài)空間模型式(1)可得

式中，u(k)和f(k)可看作外部的確定性輸入信號(hào)，當(dāng)給定期望的輸出軌跡時(shí)，用期望輸出軌跡來衡量式(3)中將來時(shí)刻的輸出預(yù)測(cè)值，以最小化量化誤差確定控制輸入，從而形成預(yù)測(cè)控制器設(shè)計(jì)問題[6].設(shè)狀態(tài)、輸入、故障輸入矩陣分別為

當(dāng)閉環(huán)系統(tǒng)是最小實(shí)現(xiàn)時(shí)，φ=A-K C是穩(wěn)定的，且x(k),u(k)和y(k)在任何時(shí)刻k都是有界的.

2 故障估計(jì)

將狀態(tài)、輸入、故障輸入矩陣代入式(3)得

以瞬時(shí)刻k-p為初始時(shí)刻，遞推迭代地求解p個(gè)采樣瞬時(shí)直至?xí)r刻k，可得關(guān)系式

將式(5)代入(4)式可得輸出方程的迭代式

輸出方程式(5)被稱為外部輸入作用下的矢量自回歸模型，輸出方程式(5)中需要辨識(shí)的參數(shù)為

所辨識(shí)的馬爾科夫參數(shù)可由原狀態(tài)空間模型式(1)通過矩陣運(yùn)算及所定義的各個(gè)矩陣得到.由于轉(zhuǎn)換過程僅僅采用初等矩陣運(yùn)算，在閉環(huán)系統(tǒng)是最小實(shí)現(xiàn)條件下的轉(zhuǎn)換過程是唯一的[7].將式(6)中的時(shí)間指標(biāo)k分別用時(shí)間序列t,t+1,···,t+N-1來替換，將y(t),y(t+1),···,y(t+N-1)整理成一個(gè)塊行矢量，記為Yid，即

根據(jù)式(5)整理矢量Yid中的各個(gè)分量，通過計(jì)算得到

記式(6)中的馬爾科夫參數(shù)序列分別為

于是可將式(6)簡(jiǎn)寫為

選擇充分大的過去辨識(shí)時(shí)域p使得‖φp‖≈0，于是忽略式(8)的首項(xiàng)未知項(xiàng)Cφpxid，僅剩下近似式

重新整理式(9)，定義殘差矢量為

式(10)為故障fid到殘差矢量rid間的關(guān)系式.對(duì)于故障fid的估計(jì)可在最小二乘意義下求解，即構(gòu)成關(guān)于故障fid的適定最小二乘優(yōu)化問題

在式(11)中可增加關(guān)于故障fid的約束條件，故障fid估計(jì)的約束優(yōu)化問題可采用梯度算法.

3 子空間預(yù)測(cè)控制下的輸出估計(jì)

子空間預(yù)測(cè)控制的首要步驟是給出將來時(shí)刻的輸出預(yù)測(cè)值，利用式(5)可計(jì)算輸出估計(jì).因殘差在多于1個(gè)采樣瞬時(shí)的滑窗時(shí)域內(nèi)生成，即滑窗時(shí)域?yàn)閇k-L+1,k]，其中L為輸出時(shí)域水平.類似地將式(8)的時(shí)間指標(biāo)k分別用L個(gè)時(shí)間指標(biāo)來代替，采集y(k-L+1)···y(k)為一列矢量，記為

類似地定義塊矢量uk,L和ek,L，記[k-L-p+1 k-L]為過去的時(shí)域滑動(dòng)窗，p為過去的時(shí)域水平，則輸出方程可整理為

記各個(gè)塊Hankel矩陣分別為

記過去滑動(dòng)窗下的輸入-輸出為

采集過去窗下的塊輸入-輸出觀測(cè)數(shù)據(jù)序列為

將式(12)簡(jiǎn)寫為緊湊形式為

式(14)可改寫為

式(15)為子空間預(yù)測(cè)控制下的輸出預(yù)測(cè)值，但等式右邊zk,L+p包含矢量yk,L，故應(yīng)將該式進(jìn)行適當(dāng)?shù)木仃囎儞Q[8]，再代入預(yù)測(cè)控制的目標(biāo)準(zhǔn)則函數(shù)中.

4 殘差分析

根據(jù)式(15)構(gòu)造殘差矢量為

因新息e(k)是白高斯信號(hào)，殘差rk,L的統(tǒng)計(jì)分布為

在輸出方程式(5)中依次取i=1,···,L,可得

則成立如下關(guān)系：

利用矩陣直積的性質(zhì)[9]有

類似有

在式(8)的左右兩邊同乘以偽逆矩陣z+id，得

計(jì)算馬爾科夫參數(shù)矩陣E1的誤差矩陣為

對(duì)式(22)取矢量化運(yùn)算可得

類似有

將y(k-L+i)按照i=1,···,L的順序排列成矢量，利用以上各等式可得

式中，數(shù)據(jù)矩陣z0可表示為

在殘差矢量的推導(dǎo)中，式(17)給出了殘差矢量rk,L的統(tǒng)計(jì)分布，式(25)給出了rk,L的具體計(jì)算結(jié)構(gòu).

5 子空間預(yù)測(cè)控制的對(duì)偶分解

子空間預(yù)測(cè)控制中將來時(shí)刻控制器輸入的求解可回歸到一個(gè)優(yōu)化問題的最優(yōu)解[10].假設(shè)期望的輸出軌跡序列γ已知，為

預(yù)測(cè)控制中常用的二次代價(jià)性能準(zhǔn)則函數(shù)為

式中，兩矩陣Q1和R可選擇為對(duì)稱正定加權(quán)矩陣，其中決策變量選擇為

為展開式(28)，需將第1項(xiàng)展開整理成uk,L的顯式形式[11].根據(jù)式(15)可得將來時(shí)刻輸出觀測(cè)數(shù)據(jù)序列的預(yù)測(cè)值為

將預(yù)測(cè)輸出值改寫為

引入3個(gè)矩陣來簡(jiǎn)記式(30)為

三矩陣Λ1、Λ2、Λ3為zk-L,p、uk,L、fk,L+p前各自的矩陣乘積.將式(31)代入預(yù)測(cè)控制二次準(zhǔn)則函數(shù)中，構(gòu)成優(yōu)化問題[12]

整理式(32)可得

對(duì)于式(33)的求解，采用快速梯度算法來解決限制控制輸入幅值的特例.對(duì)于包含有控制輸入uk,L的等式和不等式約束的最優(yōu)化問題，如果從對(duì)偶角度來闡述，就可以通過對(duì)偶分解將原約束優(yōu)化問題轉(zhuǎn)化為一個(gè)無約束優(yōu)化問題[13].將優(yōu)化問題式(33)改寫成

在式(33)中，增加的規(guī)范項(xiàng)在于保證決策變量不會(huì)發(fā)生較大的突變[14]，即構(gòu)成約束優(yōu)化問題

引入拉格朗日乘子矢量μ1、μ2、μ3，構(gòu)成對(duì)偶問題中的拉格朗日函數(shù)

對(duì)拉格朗日函數(shù)整理可得

先求解關(guān)于α的下確界運(yùn)算得

式中，[]i表示矢量的第i個(gè)元素，為求解式(36)關(guān)于uk,L的最小化運(yùn)算，設(shè)

式(36)關(guān)于uk,L的最小化運(yùn)算為

對(duì)偶問題為

定義負(fù)對(duì)偶函數(shù)為

因?yàn)閒(μ)包含有正半定海賽矩陣的二次項(xiàng)和線性項(xiàng)，則f(μ)是可微的，其對(duì)應(yīng)的梯度為

由凸優(yōu)化理論中的極小-極大定理可知[15]：▽f(μ)的最小Lipschitz常數(shù)ζ為

對(duì)于偶問題式(39)，用最近鄰梯度法可得迭代值

式中，t為迭代次數(shù)，Pμ表示μ所在集合上的投影，新迭代值μt+1是在上次迭代值的基礎(chǔ)上加上一負(fù)梯度投影方向.原優(yōu)化問題的控制輸入迭代值為

將式(37)中的vt值代入▽f(μ)中得

定義

將式(41)并行化處理為

對(duì)于等式和不等式約束的最優(yōu)化問題，雖然對(duì)偶迭代求解增加了對(duì)偶步驟，但對(duì)偶過程將約束優(yōu)化轉(zhuǎn)化為無約束優(yōu)化，從而簡(jiǎn)化了復(fù)雜度.

6 仿真算例

以無人直升機(jī)在空中的懸停狀態(tài)為例，采用本文所提出的方法進(jìn)行預(yù)測(cè)控制器的設(shè)計(jì).仿真時(shí)無人直升機(jī)帶有舵回路和清洗網(wǎng)絡(luò)的全狀態(tài)變量模型的各個(gè)系統(tǒng)矩陣式分別取為

控制輸入是方向舵和副翼的伺服輸入信號(hào).輸出矩陣C中的值57.2958是角度與弧度的轉(zhuǎn)換值.根據(jù)系統(tǒng)的狀態(tài)矩陣A可解出系統(tǒng)的開環(huán)荷蘭滾模態(tài)有極點(diǎn)-0.442 5±j3.063，可見開環(huán)荷蘭滾的阻尼是不足的，而開環(huán)螺旋模態(tài)有極點(diǎn)-0.016 31，可采用輸出反饋調(diào)節(jié)器來設(shè)計(jì)無人直升機(jī)的懸停控制.原狀態(tài)空間系統(tǒng)中過程噪聲e(k)為零均值白噪聲，其方差矩陣為Ce=10-4I5.預(yù)測(cè)和控制時(shí)域水平分別取為p=L=100.對(duì)于加性的截?cái)嗥鞴收希ＴO(shè)E=B,G=D=0.故障f(k)的估計(jì)可采用最小二乘法來辨識(shí).無人直升機(jī)中的加性故障常見于傳感器故障和截?cái)嗥鞴收?，其中的傳感器故障并不直接影響無人機(jī)的飛行特性，只是傳感器在測(cè)量過程中所出現(xiàn)的各種由于大氣干擾所導(dǎo)致的誤差.截?cái)嗥魇轻槍?duì)飛行試驗(yàn)中不可避免地受到大氣紊流激勵(lì)的影響.大氣紊流激勵(lì)作為一種不可測(cè)激勵(lì)，它所產(chǎn)生的隨機(jī)響應(yīng)將作為過程噪聲包含在測(cè)量輸出方程中.仿真中取故障估計(jì)為

選取3 000個(gè)輸入-輸出采樣數(shù)據(jù)用于閉環(huán)試驗(yàn)的預(yù)測(cè)控制器設(shè)計(jì).采樣周期取為0.2 s，初始輸出反饋控制器選擇為

圖1 二維故障估計(jì)Figure 1 Two dimensions fault estimations

參考信號(hào)γ取為零均值的白噪聲，其對(duì)應(yīng)的方差為4×10-4I.過去和將來的水平線取為p=60，從而有‖φp‖2≈0，Q1=R=0.01I，等式和不等式約束條件分別取為

增加的規(guī)范項(xiàng)為

對(duì)于子空間預(yù)測(cè)控制的優(yōu)化式(34)，采用對(duì)偶分解算法求解最優(yōu)預(yù)測(cè)控制器.為驗(yàn)證閉環(huán)系統(tǒng)性能，系統(tǒng)初始狀態(tài)選擇為x0=10I，利用操縱桿控制無人直升機(jī)俯仰姿態(tài)的具體過程如下：首先推前飛桿，使直升機(jī)向下低頭約3°，推桿動(dòng)作持續(xù)3 s后回零位；松桿停5 s，再向后拉前飛桿，使直升機(jī)向上抬頭約3°，拉桿持續(xù)3 s后再回零位，在上述操縱過程下閉環(huán)系統(tǒng)的仿真響應(yīng)曲線如圖2所示.無人直升機(jī)的姿態(tài)角響應(yīng)和角速率響應(yīng)均能達(dá)到預(yù)先期望的軌跡.對(duì)操縱桿的控制輸入曲線如圖3所示，推桿的作用力大小為時(shí)變的持續(xù)變化力.

懸停狀態(tài)下進(jìn)行仿真時(shí)，給定懸停點(diǎn)的高空離地面坐標(biāo)為(-19 m,-83 m),此時(shí)無人直升機(jī)的航向?yàn)?97.8°.無人直升機(jī)進(jìn)入懸停狀態(tài)后，利用指令盤依次發(fā)送如下各個(gè)指令：前移5 m、右移5 m、后移5 m、左移5 m，各指令發(fā)送的間隔均為15 s.利用子空間預(yù)測(cè)控制器對(duì)閉環(huán)系統(tǒng)進(jìn)行輸出反饋控制，得到無人直升機(jī)位移變化如圖4所示.

圖2 俯仰姿態(tài)控制響應(yīng)曲線Figure 2 Pitch attitude control response curves

圖3 最優(yōu)控制輸入Figure 3 Optimal control input

由圖4可得：定點(diǎn)懸停閉環(huán)控制在子空間預(yù)測(cè)控制器的作用下具有良好的控制響應(yīng)特性，穩(wěn)態(tài)定位偏差小于0.5 m，均符合無人直升機(jī)的控制精度要求.前移和后移時(shí)，其側(cè)向位移最大值小于0.5 m；左移和右移時(shí)，其前向上的位移最大值小于0.5 m，縱向和側(cè)向間的耦合較小.

7 結(jié)語(yǔ)

本文研究故障條件下的子空間預(yù)測(cè)控制問題，推導(dǎo)殘差矢量的統(tǒng)計(jì)分布具體形式.針對(duì)等式和不等式約束的預(yù)測(cè)控制器最優(yōu)化問題，通過對(duì)偶運(yùn)算并采用最近鄰梯度法即可求解最優(yōu)的預(yù)測(cè)控制器.關(guān)于算法的漸近性分析是下一步的研究重點(diǎn).

[1]CHIUSO A.The role of vector autoregressive modeling in predictor based subspace identif ication[J].Automatica,2007,43(6):1034-1048.

[2]CHIUSOA.On the relation between CCA and predictor based subspace identif ication[J].IEEE Transactions of Automatic Control,2007,52(10):1795-1811.

[3]王建宏，王道波.子空間預(yù)測(cè)控制算法在主動(dòng)噪聲振動(dòng)中的應(yīng)用[J].振動(dòng)與沖擊，2011,30(10):129-135.

WANGJianhong,WANGDaobo.Subspace predictive control applied to active noise and vibration control[J].Journal of Vibration and Shock,2011,30(10):129-135.(in Chinese)

[4]王建宏，王道波.子空間預(yù)測(cè)控制中的橢球優(yōu)化及其應(yīng)用[J].應(yīng)用科學(xué)學(xué)報(bào)，2010,28(4):424-429.

WANG Jianhong,WANG Daobo.Application of the ellipsoid optimization algorithm in subspace predictive control[J].Journal of Applied Sciences,2010,28(4):424-429.(in Chinese)

[5]王建宏.故障估計(jì)下子空間預(yù)測(cè)控制的快速梯度算法[J].上海交通大學(xué)學(xué)報(bào)：自然科學(xué)版，2013,47(7):1015-1021.

WANGJianhong.Fast gradient algorithm in subspace predictive control under fault estimation[J].Journal of Shanghai Jiaotong University:Natural Science Edition,2013,47(7):1015-1021.(in Chinese)

[6]LJUNG L.System identif ication:theory for the user[M].[S.l.]:Prentice Hall,1999.

[7]BOYD S,VANDENBERGHE L.Convex optimization[M].UK:Cambridge University Press,2008.

[8]ZEILINGER M.Real time suboptimal model predictive control using a combination of explicit MPC and online optimization[J].IEEE Transactions of Automatic Control,2011,56(7):1524-1534.

[9]JAKOB K.A design algorithm using external perturbation to improve iterative feedback tuning convergence[J].Automatica,2011,47(2):2665-2670.

[10]OHLSSON H.Identif ication of switched linear regression models using sum of norms regularization[J].Automatica,2013,49(4):1045-1050.

[11]Van MULDERSA.Identif ication of systems with localised nonlinearity:from state space to block structured models[J].Automatica,2013,49(5):1392-1396.

[12]FELLERC.An improved algorithm for combinatorial multi-parameteric quadratic programming[J].Automatica,2013,49(5):1370-1376.

[13]FIOERH.Observers for interconnected nonlinear and linear systems[J].Automatica,2012,48(7):1339-1346.

[14]RAKO L.Identif ication of switched linear systems via sparse optimization[J].Automatica,2011,47(4):668-677.

[15]PAUL M J.Identif ication of dynamic models in complex networks with prediction error methods[J].Automatica,2013,49(10):2994-3006.