陳爾康， 荊武興， 高長生

(哈爾濱工業(yè)大學(xué) 航天學(xué)院, 哈爾濱 150001)

高速飛行器指高速飛行的有翼或無翼飛行器[1]，具有強(qiáng)非線性、不確定性、強(qiáng)耦合等特點(diǎn)[2]。此外由于大量使用輕質(zhì)材料，其結(jié)構(gòu)的固有頻率較低，機(jī)體彈性形變對(duì)其控制有較大影響[3]，控制系統(tǒng)需要對(duì)彈性進(jìn)行抑制以保證穩(wěn)定性[1]?，F(xiàn)有研究多將彈性振動(dòng)看作干擾并在控制回路中加以抑制[4]，這樣雖然能夠保證飛行器的穩(wěn)定，但難以提高響應(yīng)速度和控制精度等性能。要實(shí)現(xiàn)彈性高速飛行器的精細(xì)姿態(tài)控制[3]，需要使用動(dòng)態(tài)面控制、模型預(yù)測控制等先進(jìn)控制方法，這勢必需要引入難以直接測量的彈性狀態(tài)的反饋，因此對(duì)彈性振動(dòng)狀態(tài)進(jìn)行估計(jì)十分重要。除彈性振動(dòng)狀態(tài)外，彈性振動(dòng)的相關(guān)參數(shù)如彈性模態(tài)固有頻率等也十分重要，能夠用于飛行器的故障診斷等。固有頻率與飛行器狀態(tài)相關(guān)，在飛行過程中并非常數(shù)，因此對(duì)其進(jìn)行估計(jì)是十分必要的。綜上，需要研究彈性高速飛行器的狀態(tài)/參數(shù)聯(lián)合估計(jì)方法[5]。

作為一種基于模型的遞歸濾波器，擴(kuò)展卡爾曼濾波(EKF)利用一階線性化的方法對(duì)卡爾曼濾波(KF)進(jìn)行擴(kuò)展，是當(dāng)前較為成熟和使用廣泛的狀態(tài)/參數(shù)估計(jì)方法[5]。雖然EKF可以通過將參數(shù)擴(kuò)展為狀態(tài)量的方式對(duì)參數(shù)進(jìn)行估計(jì)，但受限于噪聲模型難以估計(jì)時(shí)變參數(shù)；此外，EKF也無法處理狀態(tài)量的約束[6]。而滾動(dòng)時(shí)域估計(jì)(MHE) 方法則從最優(yōu)控制問題的角度出發(fā)[7]，并引入滾動(dòng)時(shí)域策略[8]，通過求解優(yōu)化問題實(shí)現(xiàn)對(duì)狀態(tài)/參數(shù)的估計(jì)；因此滾動(dòng)時(shí)域估計(jì)能夠較好地處理帶約束的估計(jì)問題和狀態(tài)/參數(shù)聯(lián)合估計(jì)問題[5]，近年來得到越來越廣泛的應(yīng)用。文獻(xiàn)[9]建立了滾動(dòng)時(shí)域估計(jì)漸近穩(wěn)定和有界穩(wěn)定的充分條件，并給出了一種算法實(shí)現(xiàn)。文獻(xiàn)[10]在化學(xué)工程問題中利用滾動(dòng)時(shí)域估計(jì)進(jìn)行狀態(tài)/參數(shù)聯(lián)合估計(jì)。文獻(xiàn)[11]將滾動(dòng)時(shí)域估計(jì)用于航天器的姿態(tài)估計(jì)和傳感器參數(shù)校正。文獻(xiàn)[12]則研究了機(jī)器人的多速率采樣滾動(dòng)時(shí)域估計(jì)問題。

對(duì)滾動(dòng)時(shí)域估計(jì)的研究可分為兩大類：到達(dá)代價(jià)的計(jì)算和優(yōu)化問題的求解。到達(dá)代價(jià)的計(jì)算直接關(guān)系到滾動(dòng)時(shí)域估計(jì)的準(zhǔn)確性和穩(wěn)定性。文獻(xiàn)[7]推導(dǎo)了滾動(dòng)時(shí)域估計(jì)穩(wěn)定的充分條件，提出利用估計(jì)的協(xié)方差矩陣計(jì)算到達(dá)代價(jià)，這也是目前計(jì)算到達(dá)代價(jià)的主流思路。基于文獻(xiàn)[7]的研究成果，文獻(xiàn)[13]采用卡爾曼濾波的誤差協(xié)方差矩陣計(jì)算到達(dá)代價(jià)，文獻(xiàn)[14-15]利用EKF計(jì)算到達(dá)代價(jià)，文獻(xiàn)[16]則利用無跡卡爾曼濾波(UKF)計(jì)算到達(dá)代價(jià)。

這些方法都是從概率與統(tǒng)計(jì)的角度出發(fā)，且并未利用滾動(dòng)時(shí)域估計(jì)的結(jié)果，因此是一類開環(huán)的到達(dá)代價(jià)計(jì)算方法，不利于進(jìn)一步提高估計(jì)精度。針對(duì)這一問題，本文從動(dòng)態(tài)規(guī)劃原理出發(fā)，將到達(dá)代價(jià)的計(jì)算轉(zhuǎn)化為最小二乘問題，并利用QR分解給出計(jì)算方法。該方法在到達(dá)代價(jià)的計(jì)算過程中使用了最新估計(jì)值，因而形成了反饋機(jī)制，有利于提高估計(jì)精度與速度。

1 彈性高速飛行器數(shù)學(xué)模型

本文主要研究彈性高速飛行器的縱向狀態(tài)/參數(shù)估計(jì)問題，因此使用如下縱向動(dòng)力學(xué)模型[17]：

(1)

本文考慮2種傳感器布置方案。方案1采用常規(guī)的傳感器輸出：角速度信號(hào)ωm和過載信號(hào)nm。

(2)

式中:φi(x)為第i階彈性模態(tài)的振型；xs為傳感器的軸向坐標(biāo)；nz為剛體模態(tài)產(chǎn)生的過載。

一般地，彈性高速飛行器僅一階彈性模態(tài)與剛體運(yùn)動(dòng)及控制系統(tǒng)頻帶耦合[17-18]，因此本文僅考慮一階模態(tài)。忽略長周期模態(tài)，且考慮彈性模態(tài)固有頻率的變化，令

(3)

式中:δ為飛行器的舵偏角。

方案2采用分別安裝在彈體前半部(距頭部xs1)和后半部(距頭部xs2)的2套陀螺[4]：

(4)

分別將上述信號(hào)作差，可得到與剛體運(yùn)動(dòng)無關(guān)的角速度信號(hào)分量ωe和受到彈性振動(dòng)信號(hào)擾動(dòng)的角速度信號(hào)分量ωr：

(5)

同樣可令

(6)

綜上，彈性高速飛行器的數(shù)學(xué)模型如下：

(7)

2 基于QR分解的滾動(dòng)時(shí)域估計(jì)

2.1 問題描述

對(duì)于式(7)所示的連續(xù)系統(tǒng)，可利用差分、多重打靶等方法轉(zhuǎn)化為式(8)所示的離散系統(tǒng)：

(8)

式中：wk和vk分別為系統(tǒng)噪聲序列和量測噪聲序列，通常情況下假設(shè)二者服從零均值高斯分布：

(9)

其中:Q和R分別為系統(tǒng)噪聲和量測噪聲的協(xié)方差矩陣。

系統(tǒng)式(8)的狀態(tài)估計(jì)問題可轉(zhuǎn)化為如下優(yōu)化問題：

(10)

式中：

(11)

優(yōu)化問題式(10)利用了全部的測量數(shù)據(jù)，因此稱全信息估計(jì)。全信息估計(jì)的計(jì)算量會(huì)隨T的增長而迅速增大到無法接受的地步，因此引入滾動(dòng)時(shí)域策略以限制全信息估計(jì)的維數(shù)，形成滾動(dòng)時(shí)域估計(jì)方法。滾動(dòng)時(shí)域估計(jì)方法的目標(biāo)函數(shù)ΦT變?yōu)?/p>

ΦT(xT-N,{wk})=CT-N(xT-N)+

(12)

式中：CT-N(xT-N)為到達(dá)代價(jià)，根據(jù)前向動(dòng)態(tài)規(guī)劃原理，其定義如下：

s.t.x(T,x0,{wk})=z

(13)

其中:x(τ,x0,{wk})表示在初值為x0且受到噪聲序列{wk}的情況下τ時(shí)刻狀態(tài)量的取值。

由式(13)可知，到達(dá)代價(jià)是一個(gè)表達(dá)式非常復(fù)雜的函數(shù)[19]。因此，一般用如下二次函數(shù)近似表示到達(dá)代價(jià)：

(14)

2.2 基于QR分解的到達(dá)代價(jià)計(jì)算方法

在T+1時(shí)刻，有

ΦT+1(xT-N+1,{wk})=CT-N+1(xT-N+1)+

(15)

因此，利用式(14)所示函數(shù)估計(jì)到達(dá)代價(jià)，根據(jù)前向動(dòng)態(tài)規(guī)劃原理可得

(16)

對(duì)權(quán)重矩陣PT-N、R-1和Q-1作Cholesky分解，有

(17)

利用式(17)，式(16)可寫為

CT-N+1(xT-N+1)=

(18)

(19)

式中：

(20)

同理可得

(21)

式中：

(22)

將式(19)和式(21)代入式(18)可得

(23)

令

(24)

式(23)所示的到達(dá)代價(jià)計(jì)算問題轉(zhuǎn)換為如下所示的最小二乘問題：

(25)

為求解式(25)所示的最小二乘問題，對(duì)A作QR分解。

(26)

式中:L為正交矩陣；R1和R2為上三角矩陣；R12為矩陣。

此外為表示方便，設(shè)

(27)

將式(26)和式(27)代入式(25)可得

(28)

式(28)中的變量僅有xT-N，因此該最小二乘問題的解為

(29)

將式(29)代入式(28)可得

(30)

與式(14)比較即可得到到達(dá)代價(jià)的更新方程:

(31)

綜上，到達(dá)代價(jià)的更新計(jì)算方法如下：

步驟2線性化函數(shù)F和H。

步驟3構(gòu)造矩陣A和b。

步驟4對(duì)A進(jìn)行QR分解。

并計(jì)算

即完成對(duì)到達(dá)代價(jià)的更新。

2.3 滾動(dòng)時(shí)域估計(jì)問題的求解

由2.1節(jié)和2.2節(jié)可知，彈性高速飛行器的狀態(tài)估計(jì)問題轉(zhuǎn)化為固定維數(shù)的優(yōu)化問題。但這需要對(duì)系統(tǒng)狀態(tài)和輸入進(jìn)行采樣，因此對(duì)于彈性高速飛行器這一連續(xù)系統(tǒng)，需要選用合適方法完成系統(tǒng)方程的離散。

考慮到系統(tǒng)非線性較強(qiáng)且采樣速率恒定，因此選用離散節(jié)點(diǎn)間距恒定且精度較高的多重打靶法。多重打靶法在等間距的離散節(jié)點(diǎn)上對(duì)狀態(tài)量和控制量進(jìn)行離散，認(rèn)為相鄰節(jié)點(diǎn)間控制量恒定，并加入節(jié)點(diǎn)處狀態(tài)量相等的約束，從而實(shí)現(xiàn)連續(xù)系統(tǒng)的離散化，最終將狀態(tài)估計(jì)問題轉(zhuǎn)化為式(32)所示的非線性規(guī)劃問題。

(32)

本文采用較為成熟的序列二次規(guī)劃方法求解非線性規(guī)劃問題式(32)。

綜上，滾動(dòng)時(shí)域估計(jì)方法的步驟如下：

步驟3當(dāng)k≥N時(shí)，利用序列二次規(guī)劃方法求解非線性規(guī)劃問題式(32)。

3 仿真分析

設(shè)窗口長度N=15，采樣周期Δt=15 s，量測噪聲協(xié)方差矩陣為R=diag([2.5×10-52.5×10-5])，系統(tǒng)噪聲協(xié)方差矩陣為Q=diag([1×10-61×10-21×10-62.5×10-30.36])。系統(tǒng)輸入信號(hào)如圖1所示，為了驗(yàn)證不同輸入下方法的性能，在前25 s為正弦信號(hào)，后25 s無輸入。此外，為了驗(yàn)證滾動(dòng)時(shí)域估計(jì)方法對(duì)參數(shù)的估計(jì)性能，對(duì)一階彈性模態(tài)固有頻率加入了正弦擾動(dòng)信號(hào)。

在量測數(shù)據(jù)和估計(jì)初值相同的情況下，分別采用EKF、EKF更新到達(dá)代價(jià)的滾動(dòng)時(shí)域估計(jì)(MHE-EKF)和QR分解更新到達(dá)代價(jià)的滾動(dòng)時(shí)域估計(jì)(MHE-QR)對(duì)彈性高速飛行器的狀態(tài)和一階彈性模態(tài)的固有頻率進(jìn)行估計(jì)。

在傳感器按照方案1布置的情況下，分別對(duì)上述3種方法進(jìn)行100次Monte Carlo仿真。估計(jì)結(jié)果的均方根誤差(Root Mean Square Error, RMSE)如圖2所示。

圖1 輸入信號(hào)Fig.1 Input signal

圖2 EKF、MHE-EKF和MHE-QR方法估計(jì)結(jié)果的均方根誤差Fig.2 RMSE of estimation results of EKF, MHE-EKF and MHE-QR methods

利用MHE-QR方法分別在傳感器布置方案1(MHE-QR1)和傳感器布置方案2(MHE-QR2)下進(jìn)行100次Monte Carlo仿真，估計(jì)結(jié)果的均方根誤差如圖3所示。

表1為上述方法的均方根誤差均值。為驗(yàn)證估計(jì)彈性模態(tài)頻率的必要性，表1還顯示了只估計(jì)狀態(tài)量的滾動(dòng)時(shí)域估計(jì)方法(MHE-S)和EKF(EKF-S)的均方根誤差均值。

由表1可知，MHE-QR和MHE-EKF的精度明顯高于EKF。且MHE-QR對(duì)ω1的估計(jì)精度最高，MHE-EKF的估計(jì)誤差大于MHE-QR，而EKF對(duì)ω1的估計(jì)出現(xiàn)了較大的誤差。這是由于在滾動(dòng)時(shí)域估計(jì)中，ω1只是一個(gè)優(yōu)化變量，而噪聲項(xiàng)只是提供了改變?chǔ)?取值的手段，對(duì)ω1的估計(jì)并不依賴具體的模型；而EKF將ω1看作狀態(tài)量，對(duì)其的估計(jì)精度依賴于模型的準(zhǔn)確度，但這類參數(shù)的變化并不存在具體的模型，因此EKF很難準(zhǔn)確估計(jì)時(shí)變的參數(shù)。此外，只估計(jì)狀態(tài)量的MHE-S和EKF-S的誤差均明顯大于其他3種同時(shí)估計(jì)參數(shù)和狀態(tài)量的方法，說明對(duì)變化的參數(shù)進(jìn)行估計(jì)是十分必要的。

表2顯示了3種方法的計(jì)算耗時(shí)，仿真在Windows 10系統(tǒng)(CPU為i5-7400，主頻為3.00 GHz)中MATLAB R2017a環(huán)境下完成。EKF的耗時(shí)明顯短于MHE-QR和MHE-EKF。MHE-QR的平均計(jì)算耗時(shí)短于MHE-EKF，且MHE-QR的最大計(jì)算耗時(shí)低于采樣周期，具備應(yīng)用潛力；而MHE-EKF雖然平均計(jì)算耗時(shí)低于采樣周期，但最大計(jì)算耗時(shí)高于采樣周期，表明EKF更新到達(dá)代價(jià)在計(jì)算效率上不如QR分解。

圖3 不同方案時(shí)MHE-QR方法估計(jì)結(jié)果的均方根誤差Fig.3 RMSE of estimation results of different schemes using MHE-QR method

表1 不同方法估計(jì)結(jié)果的均方根誤差均值Table 1 Average RMSE mean values of estimation results of different methods

表2 不同方法的計(jì)算耗時(shí)Table 2 Run time of different methods

4 結(jié) 論

本文提出了一種基于QR分解的到達(dá)代價(jià)更新方法，并將其用于彈性高速飛行器的滾動(dòng)時(shí)域估計(jì)中，實(shí)現(xiàn)了狀態(tài)/參數(shù)聯(lián)合估計(jì)。

1) 狀態(tài)/參數(shù)聯(lián)合估計(jì)方法的精度遠(yuǎn)高于只估計(jì)狀態(tài)的方法。由于彈性高速飛行器彈性模態(tài)的固有頻率并非常數(shù)，會(huì)隨飛行器狀態(tài)變化而變化，因此對(duì)其進(jìn)行在線估計(jì)是非常必要的，能夠有效提高狀態(tài)估計(jì)的精度。

2) 滾動(dòng)時(shí)域估計(jì)的精度明顯高于EKF。相較于傳統(tǒng)的EKF更新方法，QR分解更新到達(dá)代價(jià)在精度類似的同時(shí)，具有更快的計(jì)算速度(最大計(jì)算耗時(shí)優(yōu)于EKF)。這得益于QR分解更新到達(dá)代價(jià)的策略利用了滾動(dòng)時(shí)域估計(jì)的結(jié)果，形成了反饋機(jī)制，并通過直接求解優(yōu)化問題更新到達(dá)代價(jià)。

3) 傳感器采用布置方案2時(shí)的滾動(dòng)時(shí)域估計(jì)結(jié)果好于布置方案1。這是由于方案2通過引入有效信息預(yù)估而進(jìn)一步提升了估計(jì)效果。

4) QR分解更新到達(dá)代價(jià)的滾動(dòng)時(shí)域估計(jì)方法的最長計(jì)算耗時(shí)低于采樣速率，具有實(shí)際應(yīng)用的潛力，后續(xù)應(yīng)繼續(xù)研究更快的優(yōu)化算法，提高計(jì)算速度。