伴隨方法在亞軌道飛行器預(yù)測(cè)制導(dǎo)中的應(yīng)用

馬廣富，佘智勇

(哈爾濱工業(yè)大學(xué)航天學(xué)院，150001哈爾濱，magf@hit.edu.cn)

亞軌道飛行器具有較高的飛行速度和高度，但不具備軌道運(yùn)行的速度和高度.相對(duì)入軌飛行器而言，亞軌道飛行器飛行具有高度低，任務(wù)目標(biāo)多樣化等特點(diǎn)，因此成為近年來(lái)工程技術(shù)人員廣泛關(guān)注的問(wèn)題，現(xiàn)有文獻(xiàn)多為亞軌道飛行器再入及上升軌跡優(yōu)化及制導(dǎo)，但多數(shù)為采用間接參數(shù)優(yōu)化方法求解最優(yōu)制導(dǎo)軌跡.本文重點(diǎn)研究以固體火箭發(fā)動(dòng)機(jī)為動(dòng)力系統(tǒng)的亞軌道飛行器軌跡優(yōu)化與制導(dǎo)指令的實(shí)現(xiàn)問(wèn)題，首先基于最優(yōu)控制理論，利用改進(jìn)的伴隨方法，間接解析的在線生成最優(yōu)制導(dǎo)指令，其次利用滾動(dòng)時(shí)域控制理論，快速實(shí)現(xiàn)對(duì)最優(yōu)指令的跟蹤.

1 改進(jìn)型伴隨方法

1.1 軌跡優(yōu)化方法回顧

直接離線方法［1］是近年來(lái)飛行器軌跡優(yōu)化方法的主流，此方法將控制變量進(jìn)行離散化，同時(shí)按照優(yōu)化準(zhǔn)則，直接調(diào)整控制變量的值，直至獲得滿足性能指標(biāo)要求的最優(yōu)解，但求解過(guò)程需要的先驗(yàn)知識(shí)多，計(jì)算量大，計(jì)算耗時(shí)，并且容易陷入局部極小值而使優(yōu)化失敗.間接法利用最優(yōu)理論求解問(wèn)題［2］，間接求解控制變量，可以避免局部極小值的產(chǎn)生，但傳統(tǒng)方法收斂速度難以保證，且對(duì)協(xié)態(tài)變量初始值比較敏感［3］.與傳統(tǒng)的利用伴隨方法求解最優(yōu)控制靈敏度問(wèn)題［4］不同，增強(qiáng)型伴隨方法對(duì)前述2種方法進(jìn)行了有機(jī)的結(jié)合，對(duì)于具有一定先驗(yàn)知識(shí)的初始解，按照間接法逐次迭代，具有收斂精度高，收斂速度快，控制變量平滑等特點(diǎn).滾動(dòng)時(shí)域控制技術(shù)［5］是近年來(lái)興起的一種基于非線性模型的時(shí)變跟蹤控制方法，具有一定的預(yù)測(cè)與校正能力，能夠保證全局的收斂性與穩(wěn)定性，對(duì)于制導(dǎo)指令的實(shí)現(xiàn)具有其他控制方法無(wú)法比擬的優(yōu)勢(shì).

1.2 改進(jìn)型伴隨方法

對(duì)于具有以下形式的多輸入多輸出系統(tǒng):

其中，Z∈Rn×1，F(xiàn)∈Rn×1，對(duì)于某一近似滿足初始邊界與終端邊界的軌跡，按照某種指標(biāo)進(jìn)行逐次的迭代，直至尋得最優(yōu)軌跡為止.

設(shè)Zk(t)，k≥0為前次迭代軌跡，那么，如果第(k+1)次迭代所產(chǎn)生的Zk+1(t)與Zk(t)滿足如下關(guān)系:

其中I∈Rn×n為單位陣，β∈Rn×n，0＜‖β‖≤1，將兩次迭代的誤差表示成變分δZk=Zk+1－Zk，那么在Zk對(duì)F(Zk+1)做一階泰勒展開(kāi)，即

其中?F/?ZT為雅克比矩陣.可以構(gòu)造輔助方程=－(?F/?ZT)Tη，其中η∈Rn×1，那么

由輔助微分方程，利用已知的η(tf)，將其從tf到t0進(jìn)行反向積分，可求取輔助變量的變化歷程，那么可以求得原方程狀態(tài)的變分初始值為δZ(t0)，并以此為初始值積分，如下所示:

則由已知的當(dāng)前狀態(tài)可以求取更為滿足約束的優(yōu)化軌跡狀態(tài)為

反復(fù)進(jìn)行該過(guò)程，可以得到任意精度的最優(yōu)軌跡.

1.3 改進(jìn)型伴隨方法的收斂性

由于存在關(guān)系式(2)，那么下式則必然滿足:

2 數(shù)學(xué)模型

2.1 飛行器動(dòng)力學(xué)模型

采用非線性剛體無(wú)因次動(dòng)力學(xué)模型以便于進(jìn)行數(shù)值計(jì)算［2］，即

其中:V為速度;h為高度;θ為彈道傾角;α為攻角;P為單位重量的推力;N、A分別為單位重量的法向與軸向氣動(dòng)力.高度用R0無(wú)因次化，R0= 6 378 137 m;速度用無(wú)因次化，其中g(shù)0= 9.806 65 m/s2;時(shí)間用無(wú)因次化;重力加速度g=.

2.2 邊界條件與性能指標(biāo)

初始邊界為

終端邊界為

選擇BTT轉(zhuǎn)彎方式，不產(chǎn)生相對(duì)于風(fēng)的側(cè)滑角，只對(duì)攻角幅值以及變化率進(jìn)行限制，為

以飛行時(shí)間最小(即燃料最省)作為最優(yōu)性能指標(biāo)

2.3 最優(yōu)控制問(wèn)題

以攻角作為最優(yōu)控制量，求取其隨時(shí)間的變化歷程，在式(8)～(10)的約束條件下，最小化性能指標(biāo)式(11)，這是典型的兩點(diǎn)邊值問(wèn)題，可以利用改進(jìn)的伴隨方法來(lái)求解.

2.4 哈密頓函數(shù)及協(xié)狀態(tài)

其中

為協(xié)態(tài)變量方程.

2.5 考慮約束哈密頓函數(shù)及協(xié)狀態(tài)

不失一般性，令

為約束方程，且令S=［S1S2］T，

則利用最優(yōu)控制理論，可以得到哈密頓函數(shù)為

式中λ為協(xié)狀態(tài)向量，滿足

且λα為約束拉格朗日乘子

并且滿足

其中?i=1，2.

注1:對(duì)于滿足約束的情況，λα被屏蔽;對(duì)于在約束邊界的軌跡，λα被激活，同時(shí)最優(yōu)控制量由約束方程S獲得，λα滿足駐點(diǎn)最優(yōu)控制條件.

3 滾動(dòng)時(shí)域控制

滾動(dòng)時(shí)域控制(Receding horizon control)［6－7］是基于在線計(jì)算，不斷地根據(jù)當(dāng)前測(cè)得的系統(tǒng)狀態(tài)求解系統(tǒng)的最優(yōu)控制問(wèn)題的一種控制技術(shù).設(shè)當(dāng)前時(shí)刻為t，狀態(tài)x，在時(shí)域［t，t+T］內(nèi)，關(guān)于系統(tǒng)的最優(yōu)控制問(wèn)題為

其中T為有限時(shí)域.將當(dāng)前狀態(tài)測(cè)量值作為初始條件，并視為拉格朗日型軌跡優(yōu)化問(wèn)題進(jìn)行求解，在線計(jì)算出最優(yōu)控制解，在時(shí)域T內(nèi)執(zhí)行控制，直到系統(tǒng)獲得新的狀態(tài)測(cè)量值，并將其作為新的初始條件.以相同的方法計(jì)算出下一有限時(shí)域的最優(yōu)控制解，該過(guò)程不斷反復(fù)進(jìn)行直到滿足要求，便得到一組狀態(tài)反饋控制律.顯然，該控制方案最重要的是使性能指標(biāo)不斷減少.滾動(dòng)時(shí)域控制只要求對(duì)系統(tǒng)當(dāng)前軌跡中的狀態(tài)求解最優(yōu)控制，避免了哈密頓雅可比方法的全局性和難于計(jì)算等問(wèn)題.滾動(dòng)時(shí)域控制的系統(tǒng)漸進(jìn)穩(wěn)定性可以參考文獻(xiàn)［6－7］.滾動(dòng)時(shí)域控制實(shí)現(xiàn)的一般過(guò)程如下:

1)在［t，t+T］時(shí)間段，按二次型性能指標(biāo)求取控制u=［utut+T/N… ut+T］;

2)只取ut作為當(dāng)前時(shí)刻的系統(tǒng)輸入，驅(qū)動(dòng)系統(tǒng)沿時(shí)間軸線前進(jìn);

3)在［t+T/N，t+T/N+T］時(shí)間段求取控制u=［ut+T/Nut+2T/N… ut+T/N+T］，只取ut+T/N作為當(dāng)前時(shí)刻的系統(tǒng)輸入，驅(qū)動(dòng)系統(tǒng)沿時(shí)間軸線前進(jìn);

4)在不斷推進(jìn)的時(shí)間窗口中反復(fù)應(yīng)用滾動(dòng)時(shí)域控制技術(shù)，直至飛行任務(wù)結(jié)束.

3.1 RHC在線求解流程［8］

取權(quán)矩陣為

式中Q取為最大狀態(tài)誤差平方的倒數(shù)，R取為最大控制變量誤差平法的倒數(shù).

將時(shí)間域T做N等分，即h=T/N.將狀態(tài)變量在每個(gè)時(shí)間步長(zhǎng)上做一階泰勒展開(kāi)，即

式中:

如果令

那么性能指標(biāo)可以寫(xiě)為

如果將控制序列寫(xiě)為

那么，性能指標(biāo)可以整理為

相應(yīng)各變量均為向量，且易于從式(12)得到.如果有下式成立

那么可以參考文獻(xiàn)［10］，控制變量具有如下式所示的解形式.

那么，取v的第1個(gè)信號(hào)為當(dāng)前控制，同時(shí)在下1個(gè)仿真周期，繼續(xù)上述過(guò)程.

由求解流程可以看到，RHC控制將微分方程求解轉(zhuǎn)化為代數(shù)方程，且是確定代數(shù)方程，不存在解方程問(wèn)題，這大大加快了計(jì)算效率，因此可以作為在線求解的有效方式.

4 最優(yōu)制導(dǎo)指令的生成及實(shí)現(xiàn)

將攻角α作為最優(yōu)制導(dǎo)指令，對(duì)飛行器做在線的軌跡優(yōu)化及跟蹤，利用一階泰勒展開(kāi)［9－10］，將系統(tǒng)分解為前饋標(biāo)稱輸入與反饋控制輸入2部分，其系統(tǒng)框圖如下.

圖1 制導(dǎo)系統(tǒng)框圖

4.1 最優(yōu)制導(dǎo)指令的生成

將原系統(tǒng)狀態(tài)進(jìn)行增廣，與協(xié)狀態(tài)組成新的狀態(tài)量為

做增廣方程

對(duì)其實(shí)施改進(jìn)的伴隨方法，得到系統(tǒng)前饋標(biāo)稱輸入.

取ηV(tf)=［1 0 0 0 0 0］T，那么

同理Δh，Δθ可得，那么由關(guān)系式(4)～(6)可以得到最優(yōu)系統(tǒng)軌線.同時(shí)由最優(yōu)控制駐點(diǎn)條件?H/?α=0及上述狀態(tài)變量，可以求得最優(yōu)攻角指令.

4.2 最優(yōu)制導(dǎo)流程

最優(yōu)制導(dǎo)流程如下所示:

1)選取初值Z0(t);

2)積分方程˙Z=F(Z，t);

3)求取雅克比矩陣，構(gòu)造輔助伴隨方程;

4)取ηV(tf)=［1 0 0 0 0 0］T，對(duì)輔助伴隨方程進(jìn)行反向積分可得到ηV(t);

5)ηh(tf)=［0 1 0 0 0 0］T及ηθ(tf)=［0 01 0 00］T可以確定另外2組條件;

6)為計(jì)算考慮，為終端誤差增加調(diào)整因子，求解方程組為

其中0＜Ci≤1，i=1，2，3;

8)重新積分增廣狀態(tài)方程，如滿足精度要求則結(jié)束，否則返回2).

4.3 最優(yōu)制導(dǎo)指令的實(shí)現(xiàn)

在干擾存在條件下，將系統(tǒng)進(jìn)行如下建模

5 算例

5.1 真空環(huán)境優(yōu)化仿真設(shè)置

為驗(yàn)證改進(jìn)型伴隨方法的快速有效性，忽略飛行器所受的氣動(dòng)力，以前蘇聯(lián)RD-214火箭發(fā)動(dòng)機(jī)為動(dòng)力模型，研究亞軌道飛行器上升段的軌跡優(yōu)化問(wèn)題.各仿真參數(shù)設(shè)置見(jiàn)表1.

表1 真空環(huán)境仿真參數(shù)

表1中，狀態(tài)量的初始與終端邊界，協(xié)態(tài)變量的初始值和初始預(yù)估時(shí)間用于改進(jìn)型伴隨方法對(duì)飛行器軌跡進(jìn)行優(yōu)化，時(shí)變系統(tǒng)狀態(tài)矩陣和控制輸入矩陣的擾動(dòng)為驗(yàn)證滾動(dòng)時(shí)域控制的控制效果.

5.2 真空環(huán)境優(yōu)化仿真結(jié)果

利用伴隨方法的軌跡優(yōu)化迭代誤差分布見(jiàn)表2，其迭代收斂過(guò)程如圖2～4所示，表中No.為迭代次數(shù)，Ve為終端速度誤差，He為終端高度誤差，Se為終端傾角誤差，T為終端時(shí)間.由表2可見(jiàn)，經(jīng)過(guò)有限的6次迭代，系統(tǒng)軌線最終收斂到最優(yōu)軌線，在微型計(jì)算機(jī)上進(jìn)行運(yùn)算，周期不超過(guò)1 s，具有完全在線優(yōu)化的潛力.

表2 迭代誤差及飛行時(shí)間

圖2～4可見(jiàn)整個(gè)軌跡的收斂趨勢(shì)，在不是十分理想的初始估計(jì)前提下，改進(jìn)型伴隨方法達(dá)到了令人滿意的收斂效果，同時(shí)具有自動(dòng)尋找最優(yōu)工作時(shí)間的能力.并且與直接優(yōu)化方法相比，其控制過(guò)程更平滑，且不會(huì)陷入局部極小值，這對(duì)于制導(dǎo)控制系統(tǒng)實(shí)現(xiàn)對(duì)制導(dǎo)指令的跟蹤具有極大的優(yōu)勢(shì).

圖2 高度收斂過(guò)程

圖3 速度收斂過(guò)程

5.3 大氣層內(nèi)仿真結(jié)果

各空氣動(dòng)力表達(dá)式如下:

式中:ρ=ρ0e－QhHN為大氣密度;CN為法向力系數(shù); CA為由攻角引起的軸向力系數(shù).各參數(shù)值如表3所示.

圖4 彈道傾角收斂過(guò)程

表3 大氣層內(nèi)仿真參數(shù)

由圖5～6可見(jiàn)，伴隨方法能夠在有空氣存在情況下得到收斂的上升段最優(yōu)軌跡.

圖5 大氣內(nèi)帶有控制約束的俯仰角收斂過(guò)程

圖6 大氣內(nèi)帶有控制約束的高度收斂過(guò)程

5.4 滾動(dòng)時(shí)域控制仿真結(jié)果

在有干擾及不確定的情況下，對(duì)優(yōu)化軌線進(jìn)行實(shí)時(shí)跟蹤，擾動(dòng)參數(shù) ‖Δ~A‖、‖Δ~B‖、ΔCx、 ΔCy、ΔP的值分別為5%、5%、10%、10%、5%.

分別對(duì)小擾動(dòng)狀態(tài)方程的系統(tǒng)矩陣和輸入矩陣、阻力系數(shù)、升力系數(shù)和推力進(jìn)行擾動(dòng)，分別利用PID算法和RHC算法對(duì)最優(yōu)軌線進(jìn)行跟蹤，結(jié)果見(jiàn)圖7.

圖7 滾動(dòng)時(shí)域控制結(jié)果

標(biāo)稱攻角曲線利用改進(jìn)型伴隨方法獲取，實(shí)際攻角曲線為利用滾動(dòng)時(shí)域控制技術(shù)和PID控制所得到的校正控制.仿真結(jié)果表明，在存在系統(tǒng)狀態(tài)矩陣不確定和控制輸入擾動(dòng)的情況下，利用滾動(dòng)時(shí)域控制技術(shù)，可以實(shí)現(xiàn)對(duì)標(biāo)稱最優(yōu)軌線的在線跟蹤與預(yù)測(cè)校正.并且控制增益比傳統(tǒng)PID控制要小.

6 結(jié)論

本文提出了一種改進(jìn)型的伴隨方法，基于該方法，對(duì)亞軌道飛行器進(jìn)行了高精度、快速軌跡優(yōu)化.與間接參數(shù)化方法相比，改進(jìn)型的伴隨方法具有收斂精度高、收斂速度快、控制變量平滑等特點(diǎn)，能夠避免局部極小值的產(chǎn)生.與傳統(tǒng)伴隨方法相比，減小了對(duì)協(xié)態(tài)變量初始值的依賴，增大了協(xié)態(tài)變量的收斂域，具有直接快速在線優(yōu)化飛行器運(yùn)動(dòng)軌跡，實(shí)現(xiàn)閉環(huán)制導(dǎo)的潛力.對(duì)于利用改進(jìn)型伴隨方法得到的制導(dǎo)指令，利用滾動(dòng)時(shí)域控制技術(shù)，在外界干擾存在情況下，實(shí)現(xiàn)最優(yōu)軌線的高精度在線跟蹤與預(yù)測(cè)校正.仿真實(shí)驗(yàn)表明基于改進(jìn)的伴隨方法的滾動(dòng)時(shí)域控制技術(shù)有效、可靠，對(duì)外界干擾具有一定的魯棒性.

［1］丁洪波.亞軌道飛行器上升段軌跡優(yōu)化與快速重規(guī)劃［J］.宇航學(xué)報(bào)，2009，30(3):877－883.

［2］BURROWS R R.Variational problems and their solution by the adjoint method［R］.Huntsville，Ala:Technical Memorandum，1966，NASA-TM-X-53459.

［3］趙吉松，谷良賢，潘雷.月球最優(yōu)軟著陸兩點(diǎn)邊值問(wèn)題的數(shù)值解法［J］.中國(guó)空間技術(shù)，2009，29(4): 21－27.

［4］SILLY-CARETTE J，LAUTRU D，GATI A，et al.Optimization of the homogenization of tissues using the adjoint method and the FDTD［C］//2008 IEEE MTT-S International on Microwave Symposium Digest.Atlanta，GA:［s.n.］，2008:1361－1364.

［5］羅德林.基于滾動(dòng)時(shí)域控制的尋的末制導(dǎo)律研究［J］.南京航空航天大學(xué)學(xué)報(bào)，2005，37(1):52－56.

［6］KATSARGYRI G E，KOLMANOVSKY I，MICHELINI J，et al.Path dependent receding horizon control policies for Hybrid Electric Vehicles［C］//2009 IEEE conference on Control Applications＆Intelligent Control.St.Petersburg，:［s.n.］，2009:607－612.

［7］MAYNE D Q，MICHALSKA H.Receding horizon control of nonlinear systems［J］.IEEE transaction on automatic control，1990，35(7):814－824.

［8］LU P.Closed-form control laws for linear time-varying systems［J］.IEEE Transactions on Automatic Control，2000，45(3):537－542.

［9］WU F，PACKARD A，BALAS G.LPV control design for pitch-axis missile autopilots［C］//Proceedings of the 34th IEEE Conference on Decision and Control.New Orleans，LA:［s.n.］，1995:188－193.

［10］ZHU J J，BANKER B D.X-33 ascent flight control design by trajectory linearization-a singular perturbation approach［R］.Denver，CO:AIAA，2000:1－19.