彈性高超聲速滑翔飛行器的狀態(tài)/參數(shù)聯(lián)合估計(jì)

陳爾康，荊武興，高長(zhǎng)生

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

高超聲速滑翔飛行器一般指在大氣層內(nèi)以馬赫數(shù)大于5的高速[1]作無(wú)動(dòng)力滑翔的飛行器。為在大氣層內(nèi)高速飛行，高超聲速滑翔飛行器一般采用細(xì)長(zhǎng)體外形[2]并使用碳纖維等輕質(zhì)復(fù)合材料，結(jié)構(gòu)質(zhì)量較小，因此彈性振動(dòng)模態(tài)固有頻率較低，彈性振動(dòng)與剛體運(yùn)動(dòng)間存在不可忽略的耦合效應(yīng)[3-4]。此外，在現(xiàn)代飛行器上廣泛應(yīng)用的高增益數(shù)字控制系統(tǒng)具有較大帶寬，與高超聲速滑翔飛行器的彈性振動(dòng)模態(tài)固有頻率接近，又會(huì)引起彈性振動(dòng)和控制系統(tǒng)間相互耦合的氣動(dòng)伺服彈性問(wèn)題[5]：導(dǎo)彈的結(jié)構(gòu)彈性振動(dòng)信號(hào)與剛體運(yùn)動(dòng)信號(hào)通過(guò)傳感器送至控制系統(tǒng)，經(jīng)處理后驅(qū)動(dòng)舵面偏轉(zhuǎn)，而舵面偏轉(zhuǎn)所產(chǎn)生的力又會(huì)激勵(lì)結(jié)構(gòu)的彈性振動(dòng)。文獻(xiàn)[6]研究了彈性導(dǎo)彈的氣動(dòng)伺服彈性問(wèn)題，發(fā)現(xiàn)彈性振動(dòng)對(duì)慣導(dǎo)的測(cè)量信號(hào)有很大影響，慣導(dǎo)位置選取不當(dāng)可能導(dǎo)致失穩(wěn)。綜上，在高超聲速滑翔飛行器的總體設(shè)計(jì)和控制系統(tǒng)設(shè)計(jì)中必須考慮其結(jié)構(gòu)彈性[7]。

但是傳統(tǒng)的控制律都是基于剛體模型設(shè)計(jì)的。為抑制彈性振動(dòng)模態(tài)對(duì)控制系統(tǒng)的不利影響，可在控制回路中加入陷波器，或通過(guò)反饋撓度等測(cè)量信號(hào)進(jìn)行主動(dòng)振動(dòng)抑制控制[8]。這些方法都需要精確已知彈性振動(dòng)的模態(tài)信息，但這些信息很難準(zhǔn)確計(jì)算或直接測(cè)量。在過(guò)去的幾十年中，自適應(yīng)控制[9]等控制方法被用于彈性高超聲速飛行器的控制器設(shè)計(jì)，有效提高了控制性能。這些方法將彈性模態(tài)作為擾動(dòng)處理，雖然能夠保證穩(wěn)定性，但對(duì)控制性能有不利影響。因此，為實(shí)現(xiàn)對(duì)彈性高超聲速飛行器的高精度控制，需要對(duì)彈性狀態(tài)進(jìn)行估計(jì)，并將其用于模型預(yù)測(cè)控制[10]、滑?？刂芠11]等先進(jìn)控制方法。除狀態(tài)以外，彈性振動(dòng)模態(tài)的固有頻率和振型等參數(shù)對(duì)飛行器的動(dòng)態(tài)特性有很大影響，在狀態(tài)估計(jì)、控制器設(shè)計(jì)和故障診斷等應(yīng)用中十分重要。但這些參數(shù)很難在地面準(zhǔn)確測(cè)得，且受飛行器狀態(tài)和飛行環(huán)境影響而并非常值[12]，因此對(duì)其進(jìn)行實(shí)時(shí)估計(jì)同樣是十分必要的。針對(duì)這一問(wèn)題，文獻(xiàn)[13-14]對(duì)參數(shù)進(jìn)行估計(jì)，并將估計(jì)結(jié)果用于高超聲速飛行器的自適應(yīng)控制。但這類(lèi)方法只對(duì)參數(shù)進(jìn)行估計(jì)，無(wú)法實(shí)現(xiàn)狀態(tài)和參數(shù)的聯(lián)合估計(jì)。因此，開(kāi)展對(duì)彈性高超聲速滑翔飛行器狀態(tài)/參數(shù)聯(lián)合估計(jì)方法[15]的研究是十分必要的。

作為一種基于模型的遞歸濾波器，擴(kuò)展卡爾曼濾波(Extended Kalman Filter,EKF)利用線性化技術(shù)對(duì)卡爾曼濾波進(jìn)行擴(kuò)展以用于非線性系統(tǒng)，現(xiàn)已廣泛用于狀態(tài)/參數(shù)估計(jì)、故障診斷和組合導(dǎo)航[16]等領(lǐng)域。但是EKF難以處理約束和參數(shù)不確定性，而這正是滾動(dòng)時(shí)域估計(jì)(Moving Horizon Estimation,MHE)的優(yōu)勢(shì)。

MHE利用一段時(shí)間內(nèi)的測(cè)量值，將估計(jì)問(wèn)題轉(zhuǎn)化為優(yōu)化問(wèn)題進(jìn)行求解，因此能夠處理帶有約束的狀態(tài)與參數(shù)聯(lián)合估計(jì)問(wèn)題，且具有更好的魯棒性和估計(jì)精度[17-18]。MHE的優(yōu)化目標(biāo)函數(shù)由過(guò)程噪聲、量測(cè)噪聲和到達(dá)代價(jià)組成[19]。到達(dá)代價(jià)描述了歷史數(shù)據(jù)對(duì)估計(jì)的影響，是保證MHE穩(wěn)定性的關(guān)鍵[20]。傳統(tǒng)的到達(dá)代價(jià)計(jì)算方法從概率和統(tǒng)計(jì)的角度出發(fā)，采用估計(jì)協(xié)方差矩陣計(jì)算到達(dá)代價(jià)[21]，但這類(lèi)方法并未充分利用滾動(dòng)時(shí)域估計(jì)的結(jié)果，在非線性較強(qiáng)的情況下存在估計(jì)精度較差和結(jié)果發(fā)散的問(wèn)題。

事實(shí)上，估計(jì)精度不僅與估計(jì)方法有關(guān)，還與傳感器布置有關(guān)。對(duì)于高超聲速飛行器常用的速率陀螺，其量測(cè)信號(hào)中同時(shí)包含剛體分量和彈性分量，可以在不增加額外傳感器的情況下用來(lái)估計(jì)彈性狀態(tài)和參數(shù)。但陀螺位置會(huì)影響量測(cè)信號(hào)中彈性分量的大小，進(jìn)而影響系統(tǒng)可觀性和估計(jì)精度[22]。因此還需要對(duì)速率陀螺的布置方案展開(kāi)研究。

為了實(shí)現(xiàn)彈性高超聲速滑翔飛行器的狀態(tài)和參數(shù)聯(lián)合估計(jì)，本文提出了一種傳感器布置策略和一種利用正交三角(QR)分解更新到達(dá)代價(jià)的滾動(dòng)時(shí)域估計(jì)算法(Moving Horizon Estimation with arrival cost updated by QR decomposition,MHE-QR)。通過(guò)對(duì)陀螺位置的優(yōu)化提高了系統(tǒng)可觀性，并在最優(yōu)布置方案下利用MHE-QR算法對(duì)狀態(tài)和參數(shù)進(jìn)行聯(lián)合估計(jì)。MHE-QR算法將到達(dá)代價(jià)的計(jì)算問(wèn)題轉(zhuǎn)化為最小二乘問(wèn)題并利用QR分解進(jìn)行求解，能夠保證算法穩(wěn)定性，有效提高了估計(jì)精度和計(jì)算速度。

1 彈性高超聲速滑翔飛行器的數(shù)學(xué)模型

本文研究的彈性高超聲速滑翔飛行器具有如圖1所示的雙錐體外形，其幾何參數(shù)如表1所示。

由于高超聲速滑翔飛行器結(jié)構(gòu)彈性的主要來(lái)源是縱向彎曲，因此本文使用的縱向動(dòng)力學(xué)模型為[5]

(1)

式中：m為飛行器質(zhì)量；V為飛行器速度大??；θ為彈道傾角；L為升力；g為重力加速度；?為俯仰角；q為俯仰角速度；M為俯仰力矩；I為轉(zhuǎn)動(dòng)慣量；α為迎角；ηi、ζi、ωi和Ni分別為第i階彈性模態(tài)的廣義坐標(biāo)、阻尼比、固有頻率和廣義力。

圖1 彈性高超聲速滑翔飛行器幾何外形Fig.1 Flexible hypersonic glide vehicle geometry

表1 幾何參數(shù)Table 1 Geometry parameters

為估計(jì)飛行器狀態(tài)、彈性模態(tài)固有頻率和振型信息，彈性高速飛行器的傳感器輸出為4個(gè)速率陀螺的角速度信號(hào)qm,j[12]：

(2)

(3)

式中：w(t)為過(guò)程噪聲；v(t)為量測(cè)噪聲；函數(shù)f描述了系統(tǒng)動(dòng)力學(xué)；函數(shù)h描述了系統(tǒng)輸出。

由于MHE處理的是離散系統(tǒng)，因此利用多重打靶法將系統(tǒng)(3)離散:

(4)

2 傳感器布置策略

傳感器布置可視為一個(gè)優(yōu)化問(wèn)題，其中優(yōu)化指標(biāo)的選擇較為關(guān)鍵。由陀螺的量測(cè)模型(2)可知，陀螺測(cè)得的角速度由剛體分量和彈性分量組成。當(dāng)陀螺位于振型斜率為零處時(shí)，其測(cè)得的角速度中不含彈性分量，這對(duì)彈性狀態(tài)和參數(shù)的估計(jì)是非常不利的。較好的傳感器布置方案應(yīng)使得測(cè)量信號(hào)中包含較多的彈性分量以減小估計(jì)誤差。因此應(yīng)選擇能夠反映系統(tǒng)可觀測(cè)度的優(yōu)化指標(biāo)。傳統(tǒng)的可觀性判據(jù)只能判斷狀態(tài)是否可觀，不適合作為優(yōu)化指標(biāo)。在這方面，Popov-Belevitch-Hautus(PBH)判據(jù)的改進(jìn)形式能夠反映系統(tǒng)可觀測(cè)度，可選為性能指標(biāo)[22-23]。該性能指標(biāo)由PBH判據(jù)導(dǎo)出，PBH判據(jù)的具體內(nèi)容如下所述。

考慮式(5)所示的系統(tǒng):

(5)

若對(duì)于Ae的每一個(gè)特征值λi，式(6)均成立，則該系統(tǒng)是可觀的。

rank(PBH(λi,Ae,Ce))=n

(6)

式中：n為Ae的維數(shù)，PBH(λi,Ae,Ce)為式(7)所示的PBH矩陣：

(7)

式中：In為n維單位矩陣。

PBH矩陣的每個(gè)條件數(shù)能夠描述其對(duì)應(yīng)模態(tài)的可觀測(cè)度。每個(gè)條件數(shù)越大，其對(duì)應(yīng)的模態(tài)就越接近不可觀[21]。因此，優(yōu)化指標(biāo)J選為所有PBH矩陣條件數(shù)中的最大值，即

J=max(cond(PBH(λi,Ae,Ce)))i=1,2,…,n

(8)

式中：Ae和Ce可由式(1)～式(3)經(jīng)線性化得到。

在僅考慮一階彈性振動(dòng)模態(tài)的情況下，傳感器位置與性能指標(biāo)間的關(guān)系如圖2所示。當(dāng)陀螺位于飛行器中部時(shí)，性能指標(biāo)最大，可觀測(cè)度最差。此時(shí)陀螺測(cè)量信號(hào)中幾乎不含彈性分量，對(duì)狀態(tài)/參數(shù)的估計(jì)非常不利。而布置在飛行器頭部或尾部時(shí)，陀螺測(cè)量信號(hào)中彈性分量較大，對(duì)狀態(tài)/參數(shù)估計(jì)較為有利。因此上述優(yōu)化指標(biāo)的選擇是合理的。

在考慮三階彈性模態(tài)的情況下，單傳感器位置與性能指標(biāo)間的關(guān)系如圖3所示。受二階和三階彈性模態(tài)的影響，一些位置的可觀測(cè)度變差，傳感器位置與性能指標(biāo)間的關(guān)系變的復(fù)雜。

圖2 單傳感器位置與性能指標(biāo)間的關(guān)系(一階)Fig.2 Relationship between single sensor position and performance index (one order)

圖3 單傳感器位置與性能指標(biāo)間的關(guān)系(三階)Fig.3 Relationship between single sensor position and performance index(three orders)

因此，多傳感器的布置方案需要通過(guò)優(yōu)化確定。最終，多傳感器布置問(wèn)題被轉(zhuǎn)化為如下優(yōu)化問(wèn)題：

s.t.xmin

xm,i+1-xm,i>lmin

(9)

式中：xmin和xmax分別為陀螺到飛行器頭部的最小和最大距離；lmin為相鄰兩個(gè)陀螺間的最小距離。本文中xmin、xmax和lmin的取值如表2所示。

利用序列二次規(guī)劃算法求解優(yōu)化問(wèn)題(9)即可得到最優(yōu)傳感器布置方案，如表3所示。

表2 傳感器布置參數(shù)Table 2 Sensor placement parameters

表3 傳感器布置方案Table 3 Sensor placement scheme

3 滾動(dòng)時(shí)域估計(jì)

3.1 問(wèn)題描述

在本文中，為表示方便，采用如下記法：

(10)

(11)

式中：a為m維向量；W為m×m維矩陣。

記滾動(dòng)窗口長(zhǎng)度為N。MHE利用當(dāng)前時(shí)刻之前N個(gè)時(shí)間點(diǎn)的數(shù)據(jù)，將狀態(tài)/參數(shù)估計(jì)問(wèn)題轉(zhuǎn)化為優(yōu)化問(wèn)題：

(12)

s.t.

(13)

式中：wk為過(guò)程噪聲；Q為過(guò)程噪聲協(xié)方差矩陣，vk為量測(cè)噪聲；R為量測(cè)噪聲協(xié)方差矩陣。ΘT-N(xT-N)為到達(dá)代價(jià)，到達(dá)代價(jià)概括了滾動(dòng)窗口之前的數(shù)據(jù)，對(duì)MHE的性能和穩(wěn)定性有很大的影響。

3.2 到達(dá)代價(jià)更新算法

根據(jù)MHE問(wèn)題的定義和前向動(dòng)態(tài)規(guī)劃原理，到達(dá)代價(jià)的定義為[20]：

(14)

s.t.

xk+1=F(xk,uk)+wkk=0,1,…,T-N-1

(15)

由式(14)和式(15)可知，到達(dá)代價(jià)的計(jì)算較為復(fù)雜，對(duì)于約束非線性系統(tǒng)更是無(wú)法得到解析表達(dá)式[20]。因此一般采用二次函數(shù)估算到達(dá)代價(jià)，即

(16)

PT-N=Q+GPT-N-1GT-

GPT-N-1C(R+CPT-N-1CT)CPT-N-1GT

(17)

(18)

在本文中，利用式(17)更新到達(dá)代價(jià)的MHE算法稱(chēng)為利用EKF更新到達(dá)代價(jià)的滾動(dòng)時(shí)域估計(jì)方法(Moving Horizon Estimation with arrival cost updated by EKF,MHE-EKF)。MHE-EKF算法未能充分利用MHE的估計(jì)結(jié)果，難以進(jìn)一步提高M(jìn)HE的估計(jì)精度和計(jì)算速度。對(duì)此，本文在確定性框架下給出一種利用QR分解更新到達(dá)代價(jià)的算法。

由式(14)可知，T+1時(shí)刻的到達(dá)代價(jià)為

(19)

根據(jù)前向動(dòng)態(tài)規(guī)劃原理和式(14)，式(19)可改寫(xiě)為

(20)

s.t.

(21)

利用近似到達(dá)代價(jià)代替到達(dá)代價(jià)，將式(16)代入式(19)可得

(22)

為表示方便，對(duì)權(quán)重矩陣作Cholesky分解可得

(23)

(24)

(25)

式中：

(26)

將式(25)代入式(24)可得

(27)

至此，到達(dá)代價(jià)的更新問(wèn)題被轉(zhuǎn)化為如下最小二乘問(wèn)題：

(28)

式中：

(29)

為保證數(shù)值穩(wěn)定性和提高計(jì)算速度，使用QR分解求解這一最小二乘問(wèn)題。對(duì)A作QR分解可得

(30)

式中：E為正交矩陣；F1和F2為上三角矩陣。設(shè)

(31)

由式(29)～式(31)可得

(32)

顯然可得此最小二乘問(wèn)題的解為

(33)

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

(34)

由式(34)可得到達(dá)代價(jià)更新的遞推方程為

(35)

基于QR分解的到達(dá)代價(jià)更新算法根據(jù)上一時(shí)刻MHE的估計(jì)結(jié)果計(jì)算到達(dá)代價(jià)，且能夠保證數(shù)值穩(wěn)定性，能夠提高M(jìn)HE的估計(jì)精度和計(jì)算速度。此外，該到達(dá)代價(jià)更新算法還具有如下性質(zhì)。

(36)

(37)

注2該到達(dá)代價(jià)更新算法給出的到達(dá)代價(jià)是有界的，說(shuō)明歷史數(shù)據(jù)對(duì)估計(jì)結(jié)果的影響是有限度的，不會(huì)無(wú)限增大。即對(duì)于任意向量v，有

(38)

式中：λmax為矩陣Q-1最大的特征值。

由式(29)可得

(39)

由式(30)可得

(40)

比較式(39)和式(40)等號(hào)右側(cè)矩陣第2行第2列元素可得

(41)

由式(35)和式(41)可知，對(duì)于任意向量v，有

(42)

注3該到達(dá)代價(jià)更新算法能夠保證MHE的穩(wěn)定性。

根據(jù)式(20)，該到達(dá)代價(jià)算法具有如下性質(zhì)：

(43)

反復(fù)應(yīng)用式(43)可得

(44)

式(44)表明該達(dá)到代價(jià)更新算法滿足文獻(xiàn)[20] 中的條件C2：

(45)

因此根據(jù)文獻(xiàn)[20]中的命題3.4，該到達(dá)代價(jià)算法能夠保證MHE的穩(wěn)定性。

3.3 MHE的數(shù)值解法

在3.2節(jié)提出的到達(dá)代價(jià)更新算法的基礎(chǔ)上，利用序列二次規(guī)劃算法求解式(46)所示的優(yōu)化問(wèn)題即可實(shí)現(xiàn)對(duì)狀態(tài)和參數(shù)的估計(jì)，形成MHE-QR算法。

(46)

s.t.

(47)

MHE算法的流程如圖4所示。

圖4 滾動(dòng)時(shí)域估計(jì)算法流程圖Fig.4 Flow chart of MHE

4 仿真分析

在仿真中，窗口長(zhǎng)度設(shè)為N=15，采樣周期為Δt=0.05 s，量測(cè)噪聲方差矩陣為R=diag(2.5×10-5,2.5×10-5,2.5×10-5,2.5×10-5)，系統(tǒng)噪聲方差矩陣為Q=diag(1×10-6,1× 10-2,1×10-6,2.5×10-3)。對(duì)一階彈性模態(tài)固有頻率和振型斜率加入了正弦擾動(dòng)信號(hào)。為了驗(yàn)證不同輸入下算法的性能，系統(tǒng)輸入信號(hào)如圖5 所示。

為驗(yàn)證傳感器布置策略，并分析MHE-QR、MHE-EKF和EKF 3種算法的性能，分別在表4所示的仿真場(chǎng)景下進(jìn)行仿真。

在相同條件下分別使用MHE-QR(場(chǎng)景1)、MHE-EKF(場(chǎng)景2)和EKF(場(chǎng)景3)3種算法進(jìn)行蒙特卡羅仿真，均方根誤差(Root Mean Square Error,RMSE)如圖6和圖7所示。

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

表4 仿真場(chǎng)景Table 4 Simulation cases

圖6 不同算法下?、q和η1的均方根誤差Fig.6 RMSE of ?、qandη1 with different algorithms

圖7 不同算法下和的均方根誤差 different algorithms

同樣使用MHE-QR算法，分別在傳感器最優(yōu)布置(場(chǎng)景1)、傳感器非最優(yōu)布置(場(chǎng)景4)和不估計(jì)振型斜率(場(chǎng)景5)的情況下進(jìn)行蒙特卡羅仿真，均方根誤差如圖8和圖9所示。

由圖8和圖9可知，由于傳感器并非最優(yōu)布置，場(chǎng)景4的RMSE存在較多的突然增大的情況，估計(jì)誤差大于采用最優(yōu)傳感器最優(yōu)布置的場(chǎng)景1。而在15～20 s的時(shí)間段內(nèi)，由于系統(tǒng)狀態(tài)已經(jīng)收斂至穩(wěn)態(tài)值，雖然場(chǎng)景1的均方根誤差大于場(chǎng)景4和場(chǎng)景5，但場(chǎng)景4和場(chǎng)景5均存在均方根誤差突然增大的情況，這對(duì)實(shí)際應(yīng)用非常不利。因此，仿真結(jié)果表明傳感器非最優(yōu)布置對(duì)估計(jì)精度有不利影響。此外，受不對(duì)振型斜率進(jìn)行估計(jì)的影響，場(chǎng)景5中彈性狀態(tài)和彈性模態(tài)固有頻率的RMSE均顯著大于場(chǎng)景1和場(chǎng)景4，說(shuō)明振型斜率的偏差對(duì)估計(jì)結(jié)果影響較大。因此，傳感器布置策略和包含振型斜率的狀態(tài)/參數(shù)聯(lián)合估計(jì)能夠有效提高估計(jì)精度。

圖8 不同傳感器布置下?、q和η1的均方根誤差Fig.8 RMSE of ?、q and η1 with different sensor placements

圖9 不同傳感器布置下的均方根誤差Fig.9 RMSE of with different sensor placements

對(duì)于估計(jì)算法的實(shí)際應(yīng)用，計(jì)算耗時(shí)是非常重要的。表5給出了不同場(chǎng)景下的計(jì)算耗時(shí)，仿真在Windows 10系統(tǒng)(CPU為i5-7400，主頻為3.00 GHz)中MATLAB R2017a環(huán)境下完成。其中EKF(場(chǎng)景3)的耗時(shí)明顯短于其他算法，而MHE-QR(場(chǎng)景1、4、5)的平均耗時(shí)短于MHE-EKF。雖然同樣使用MHE-QR算法，但傳感器非最優(yōu)布置的場(chǎng)景4的最大耗時(shí)是所有場(chǎng)景中最大的，說(shuō)明傳感器最優(yōu)布置能夠縮短計(jì)算耗時(shí)。由于估計(jì)狀態(tài)數(shù)量最小，場(chǎng)景5的平均耗時(shí)和最大耗時(shí)都是使用MHE-QR算法的所有場(chǎng)景中最小的，但其估計(jì)精度較差，難以實(shí)際應(yīng)用。而場(chǎng)景1則在估計(jì)精度和計(jì)算耗時(shí)之間做到了較好的平衡。得益于傳感器最優(yōu)布置和基于QR分解的到達(dá)代價(jià)更新策略，場(chǎng)景1的估計(jì)精度是最高的，但其平均耗時(shí)和最大耗時(shí)都僅長(zhǎng)于使用EKF的場(chǎng)景2和不估計(jì)振型斜率的場(chǎng)景5，且其最大耗時(shí)并未超出采樣周期，表明其具備實(shí)時(shí)應(yīng)用潛力。因此，本文提出的傳感器布置策略和到達(dá)代價(jià)更新策略不僅能夠提高估計(jì)精度，也能夠提高計(jì)算速度。

表5 計(jì)算耗時(shí)Table 5 CPU time consumption

5 結(jié) 論

1)為提高彈性高超聲速滑翔飛行器的狀態(tài)/參數(shù)估計(jì)精度，本文提出了一種傳感器布置優(yōu)化方法。該優(yōu)化方法選擇PBH矩陣條件數(shù)的最大值為優(yōu)化指標(biāo)。仿真結(jié)果表明，傳感器布置優(yōu)化方法能夠有效提高狀態(tài)/參數(shù)的估計(jì)精度。

2)針對(duì)彈性高超聲速滑翔飛行器的狀態(tài)/參數(shù)聯(lián)合估計(jì)問(wèn)題，本文提出一種新型的滾動(dòng)時(shí)域估計(jì)算法——MHE-QR算法。該算法將到達(dá)代價(jià)的估計(jì)問(wèn)題轉(zhuǎn)化為最小二乘問(wèn)題，并給出一種新型的利用QR分解更新到達(dá)代價(jià)的算法。MHE-QR算法限制了到達(dá)代價(jià)過(guò)度增長(zhǎng)，能夠保證滾動(dòng)時(shí)域估計(jì)算法的穩(wěn)定性。仿真結(jié)果表明，MHE-QR算法的估計(jì)精度高于EKF和MHE-EKF算法。此外，MHE-QR的計(jì)算耗時(shí)小于采樣周期，具備實(shí)時(shí)應(yīng)用的潛力。

后續(xù)需要研究滾動(dòng)時(shí)域估計(jì)算法的數(shù)值求解方法，以進(jìn)一步提高計(jì)算速度，滿足實(shí)時(shí)應(yīng)用的要求。