楊 良,鄭宗貴,徐 衡,陳萬春,周 浩

(1.北京航空航天大學(xué) 宇航學(xué)院,北京 100191;2.第二炮兵裝備研究院,北京 100085;3.北京空天技術(shù)研究所,北京 100074)

高超聲速飛行器在低層稠密大氣中的飛行環(huán)境特別惡劣,因此,需要設(shè)計滿足多約束要求的末段制導(dǎo)方法。要實(shí)現(xiàn)對終端落點(diǎn)和落角的控制,對動力學(xué)方程進(jìn)行線化后的制導(dǎo)律能很好地滿足要求[1-3]。然而對于終端速度,由于其非線性較強(qiáng),無法線性化,控制相當(dāng)困難。對此,趙漢元提出了調(diào)節(jié)攻角的速度控制制導(dǎo)律[4],但該制導(dǎo)律需要飛行器進(jìn)行快速滾轉(zhuǎn)運(yùn)動,對于多類飛行器而言并不適用。

近年來,數(shù)值方法和計算機(jī)技術(shù)的發(fā)展大大提高了求解非線性最優(yōu)控制問題的速度。偽譜法[5-6]就是其中一種。Ross I通過勒讓德偽譜法提出了一種基于時鐘序列保持的偽譜反饋制導(dǎo)策略,并進(jìn)行了衛(wèi)星控制和再入制導(dǎo)的在線反饋研究[7-8]。崔鋒通過勒讓德偽譜法也進(jìn)行了再入制導(dǎo)的在線反饋研究[9]。偽譜反饋為終端多約束問題尤其是終端速度控制提供了一種很好的求解方法。

本文進(jìn)行了在線高斯偽譜的反饋方法研究,建立了高斯偽譜法的初值猜測策略,并應(yīng)用于高超聲速飛行器的末段飛行,成功實(shí)現(xiàn)了對落角、落速的控制,在無反饋信息的最后一個制導(dǎo)周期引入比例導(dǎo)引方法[10]保證落點(diǎn)精度。蒙特卡羅打靶實(shí)驗(yàn)表明,該制導(dǎo)方法能夠很好地滿足多項(xiàng)終端約束條件,并且具有很高的魯棒性。高斯偽譜法的微分矩陣滿秩,因此其計算結(jié)果相比勒讓德偽譜具有更高的精度[11]。

1 方法介紹

1.1 方法概述

本文采用高斯偽譜法對非線性狀態(tài)方程進(jìn)行末端軌跡優(yōu)化;在一個制導(dǎo)周期內(nèi),采用時鐘保持機(jī)制將最優(yōu)控制指令代入實(shí)際飛行過程;在同一制導(dǎo)周期上,利用上一周期的全彈道信息及當(dāng)前狀態(tài)量進(jìn)行制導(dǎo)指令的在線重構(gòu);在下一制導(dǎo)周期上,重復(fù)上述步驟,并在最后一個制導(dǎo)周期應(yīng)用比例導(dǎo)引方法擊中目標(biāo),如圖1所示。

圖1 方法流程圖

1.2 高斯偽譜反饋與指令保持

高斯偽譜反饋是通過高斯偽譜法求解基本的最優(yōu)控制問題,獲得離散的時變反饋控制指令。反饋的目的是消除由模型不確定性、外部干擾以及估計誤差造成的影響。指令的保持是根據(jù)上一步求得的以時間為自變量的控制指令,采用圖2的形式進(jìn)行保持。

圖2 指令保持示意圖[7]

圖2中,xR(ti)為當(dāng)前時刻的狀態(tài)量,u為控制量,ti和ti+1分別為當(dāng)前時刻和下一個制導(dǎo)時刻。文獻(xiàn)[5]指出成功運(yùn)用反饋控制的關(guān)鍵是能夠快速地得到一條開環(huán)控制,高斯偽譜法對光滑問題具有的指數(shù)收斂特性可以保證這種快速性。

1.3 初值估計策略

對于非線性規(guī)劃問題,為了加快尋優(yōu)算法的速度,可以充分利用上一制導(dǎo)周期的最優(yōu)結(jié)果:

(1)

式中:t0是初始時刻;tf是終端時刻,τi∈[-1,1]是高斯拉格朗日(LG)節(jié)點(diǎn);ti,tj∈[t0,tf],是LG節(jié)點(diǎn)在時域內(nèi)的分布;Li(t)為Lagrange插值多項(xiàng)式;X(ti)為上一制導(dǎo)周期優(yōu)化的LG節(jié)點(diǎn)上的狀態(tài)量和控制量。

對于當(dāng)前制導(dǎo)周期,去除上一制導(dǎo)周期所用時間tc,可以得到新的時間節(jié)點(diǎn):

進(jìn)而可以使用上一制導(dǎo)周期的優(yōu)化結(jié)果來估計當(dāng)前制導(dǎo)周期的初值:

通過后面的仿真結(jié)果可以看出,將X(tj)作為初值代入非線性規(guī)劃問題求解器進(jìn)行計算具有較高的計算效率。

2 飛行器動力學(xué)模型

2.1 動力學(xué)方程

考慮大氣相對于地球相對靜止的無動力三自由度質(zhì)點(diǎn)彈道運(yùn)動學(xué)方程[12]:

式中:x,y,z為位置坐標(biāo);v為速度;m為飛行器質(zhì)量;γ為彈道傾角;φ為航向角;α為攻角;ε為側(cè)傾角;g=μ/r2是重力加速度,r為地心距,μ為地球重力加速度常數(shù);u1,u2為控制量;升力FL和阻力FD為

式中:q=ρv2/2為動壓,CL為升力系數(shù),CD為阻力系數(shù),Sref為飛行器參考面積;ρ為大氣密度,其表達(dá)式為

ρ=ρ0e-βH

(6)

式中:ρ0為海平面大氣密度,H為高度,β為密度常數(shù)。

2.2 約束條件

從結(jié)構(gòu)和熱防護(hù)的角度出發(fā),要求高超聲速飛行器在飛行過程中滿足動壓和過載約束:

式中:qmax和nz,max根據(jù)具體飛行器來選取。

考慮到高超聲速飛行器在末段飛行過程中的惡劣環(huán)境,除了要求攻角α和傾側(cè)角ε應(yīng)該保持在一定的范圍之內(nèi),其變化率也必須保持在一定的范圍之內(nèi),即

2.3 優(yōu)化性能指標(biāo)

為使優(yōu)化得到的攻角和傾側(cè)角足夠光滑,選取各控制量的平方和最小作為性能指標(biāo)。同時,為了增加尋優(yōu)算法的收斂區(qū)域,增強(qiáng)制導(dǎo)算法的魯棒性,把末端狀態(tài)量的約束加入到性能指標(biāo)中,具體形式為

式中:K為加權(quán)系數(shù)。

對于某些特殊的擦邊界最優(yōu)控制問題,由于原問題在擦邊界的過程中不再是光滑的最優(yōu)控制問題,偽譜法不再具有指數(shù)收斂特性,即失去解算的快速性。這時通過改變性能指標(biāo)可以將該問題轉(zhuǎn)化為光滑問題,進(jìn)而獲得指數(shù)收斂特性,可以進(jìn)行在線優(yōu)化。

3 仿真結(jié)果及分析

3.1 飛行器模型

本文采用的飛行器模型為美國在“獵鷹”計劃支持下發(fā)展的新一代高超聲速遠(yuǎn)程快速精確打擊武器系統(tǒng)——通用航空飛行器(CAV)。該飛行器具有乘波體構(gòu)型,依靠氣動力在大氣層內(nèi)滑行,具有升阻比大、飛行馬赫數(shù)高、飛行空間環(huán)境惡劣、對控制系統(tǒng)要求較高等特點(diǎn)。NASA公布的CAV的總體參數(shù)如表1所示。

表1 飛行器總體參數(shù)

CAV在低層稠密大氣的氣動系數(shù)公式為

其中各氣動參數(shù)可參考文獻(xiàn)[13]。

3.2 任務(wù)要求

飛行器的任務(wù)要求是在初始拉偏和過程拉偏的情況下滿足以下終端約束條件:

攻角、側(cè)傾角以及控制量的變化范圍為

任務(wù)要求的主要難點(diǎn)在于速度控制,飛行器的速度對升力、阻力非常敏感,也是長周期的控制過程,調(diào)節(jié)能力非常有限,因此要求制導(dǎo)律在末段飛行段進(jìn)行精確的控制。

3.3 不同初始條件下的制導(dǎo)結(jié)果

考慮飛行器不同初始條件下的無干擾彈道,可以驗(yàn)證高斯偽譜反饋制導(dǎo)方法對不同初始條件的適應(yīng)性。表2列出了5種初始條件,初速均為1 718 m/s。表3為終端狀態(tài)約束。

表2 初始條件

表3 終端狀態(tài)約束

圖3為空間彈道曲線,可以看出不同初始條件下的彈道均能滿足終端狀態(tài)約束,并且以相似形式擊中目標(biāo)。

圖4和圖5分別為不同初始條件下的速度曲線和彈道傾角曲線。由于終端狀態(tài)相同,不同初始條件的彈道在末段飛行后期差異不大,因此對于彈道的調(diào)整主要集中在末段飛行的前期。雖然末段飛行前期的速度和彈道傾角差異較大,但都能滿足終端馬赫數(shù)為3和終端彈道傾角為-80°的約束條件。

圖3 空間彈道

圖4 速度曲線

圖5 彈道傾角曲線

圖6為不同初始條件下的攻角變化曲線。攻角在末段飛行后期差異不大,在末段飛行前期差異很大,彈道的調(diào)節(jié)也主要是集中在末段飛行的前期。不同初始條件下的終端攻角都能達(dá)到0°的約束條件。

圖7為每一制導(dǎo)周期內(nèi)高斯偽譜法單次優(yōu)化的計算時間。5種初始條件情況均在主頻3.3 GHz計算機(jī)上使用Matlab 7.0進(jìn)行仿真,計算時間tcal均集中在0.1 s左右。由于不存在過程干擾,優(yōu)化時間主要用于克服LG節(jié)點(diǎn)之間的龍格現(xiàn)象。結(jié)果表明,選取0.8 s的制導(dǎo)周期可以滿足指令解算的需要,同時相比50 s以上的末段飛行時間也可以保證合理的制導(dǎo)精度。

圖6 攻角曲線

圖7 計算時間

表4為不同初始條件下的終端狀態(tài)偏差。5種初始條件下的脫靶量M和速度偏差都能保證在百分之一的數(shù)量級,終端彈道傾角和攻角偏差能保證在萬分之一的數(shù)量級。

表4 終端狀態(tài)偏差

3.4 擾動模型

建立風(fēng)速、大氣密度和溫度隨高度變化的拉偏模型:

式中:v(H),ρ(H),T(H)分別為風(fēng)速、大氣密度和溫度干擾隨高度變化的函數(shù),v(H)最大能達(dá)到50 m/s,ρ(H)最大能達(dá)到0.09ρst,T(H)最大能達(dá)到30 K。Δv,Δρ,ΔT為[-200%,200%]的拉偏系數(shù)。

飛行器的氣動拉偏模型是在標(biāo)準(zhǔn)氣動升阻力系數(shù)的基礎(chǔ)上加上偏差值:

考慮擾動模型后,飛行器速度控制的難度增加,需要在末段飛行段之前的巡航段對飛行器進(jìn)行關(guān)機(jī)控制。假設(shè)飛行器在巡航段勻速飛行,可以利用高斯偽譜法尋找最優(yōu)關(guān)機(jī)時刻,并作為初值代入末段飛行段中。

3.5 過程擾動下的蒙特卡洛打靶結(jié)果

為了測試高斯偽譜反饋制導(dǎo)方法的魯棒性,在帶關(guān)機(jī)控制的情況下進(jìn)行1 000次蒙特卡羅打靶實(shí)驗(yàn),其中考慮初始狀態(tài)拉偏項(xiàng)和3.4中所述拉偏項(xiàng),如表5所示,表中,μ為均值,3σ為均方差。由于速度對升力、阻力較為敏感,且速度控制能力有限,因此升力、阻力拉偏相對較小。

表5 蒙特卡洛打靶拉偏分布

圖8和圖9為終端位置在x和y方向上的誤差PΔx,PΔy概率分布,各個方向的誤差都在0.1 m以內(nèi)。圖10和圖11為終端速度和終端彈道傾角的誤差概率PΔv,PΔγ分布,終端速度和終端彈道傾角大部分分布在約束范圍之內(nèi)。圖12為終端攻角的誤差概率PΔα分布,終端攻角主要集中在-0.7°附近,滿足約束條件。

圖8 終端位置x方向誤差概率分布

圖9 終端位置y方向誤差概率分布

圖10 終端速度誤差概率分布

圖11 終端彈道傾角誤差概率分布

圖12 終端攻角誤差概率分布

對終端狀態(tài)偏差進(jìn)行統(tǒng)計,脫靶量的期望為0.092 m,落入概率為100%;速度誤差的期望為-6.913 m/s,落入概率為81.2%;彈道傾角偏差的期望為0.040 40°,落入概率為93.3%;攻角偏差的期望為-0.611°,落入概率為99.7%。因此,飛行器終端脫靶量、速度偏差、彈道傾角偏差以及攻角偏差的期望都很小,滿足約束條件。值得注意的是,攻角、彈道傾角和速度需要進(jìn)行全程能量管理才能保證,這反映出在線高斯偽譜制導(dǎo)的效果。落入概率為終端狀態(tài)偏差落入約束范圍內(nèi)的概率。終端速度偏差的落入概率最低,說明終端速度最難控制。蒙特卡羅打靶結(jié)果表明,高斯偽譜反饋制導(dǎo)方法在惡劣的飛行環(huán)境下具有很好的制導(dǎo)能力和魯棒性。

4 結(jié)束語

本文針對高超聲速飛行器,基于高斯偽譜法,建立了一種在末段飛行段滿足多種約束條件的在線反饋制導(dǎo)方法。使用上一制導(dǎo)周期優(yōu)化結(jié)果作為初值猜測的策略,在解算過程中可以有效地加速算法收斂。在最后一個制導(dǎo)周期加入比例導(dǎo)引方法可以保證終端命中點(diǎn)具有較高的精度。不同初始條件下的制導(dǎo)結(jié)果說明,該方法具有很好的適用性。各種拉偏情況下的蒙特卡洛打靶實(shí)驗(yàn)表明,該方法在惡劣的飛行環(huán)境下具有較好的穩(wěn)定性和魯棒性,能有效地對終端速度進(jìn)行控制。

[1] KIM M,GRIDER K V.Terminal guidance for impact attitude angel constrained flight trajecories[J].IEEE Transactions on Aerospace and Electronic Systems,1973,10(6):852-859.

[2] LEE Yong-in,RYOO C K,KIM E.Optimal guidance with constraints on impact angle and terminal acceleration[R].AIAA 2003-5795,2003.

[3] LU P,DOMAN D B,SCHIERMAN J D.Adaptive terminal guidance for hypervelocity impact in specified direction[J].Journal of Guidance,Control,and Dynamics,2006,29(2):269-278.

[4] 趙漢元.飛行器再入動力學(xué)與制導(dǎo)[M].長沙:國防科技大學(xué)出版社,1997.

ZHAO Han-yuan.Entry dynamics and guidance[M].Changsha:National University of Defence Technology Press,1997.(in Chinese)

[5] DAVID A B,GEOFFREY T H.Direct trajectory optimization and costate estimation via an orthogonal collocation method[J].Journal of Guidance,Control,and Dynamics,2006,29(6):1 435-1 440.

[6] FAHROO F,ROSS I M.Costate estimation by a Legendre pseudospectral method[R].AIAA-98-4 222,1998.

[7] ROSS I M,SEKHAVAT P.Optimal feedback control:foundations,examples,and experimental results for a new approach[J].Journal of Guidance,Control,and Dynamics,2008,31(2):307-321.

[8] BOLLONO K P,ROSS I M.A pseudospectral feedback method for real-time optimal guidance of reentry vehicles[C]//Proceedings of the 2007 American Control Conference.New York City,USA:The American Antomatic Control Council,2007:11-13.

[9] 崔鋒.采用偽譜法的再入飛行器最優(yōu)反饋制導(dǎo)方法[J].中國制造業(yè)信息化,2011,40(19):42-45.

CUI Feng.The optimal feedback guidance control method of entry vehicles based on pseudospectral method[J].Manufacture Information Engineering of China,2011,40(19):42-45.(in Chinese)

[10] PAUL Z.Tactical and strategic missile guidance[M].Cambridge,Massachusetts:Charles Stark Draper Laboratory,Inc.,1997.

[11] GARG D,PATTERSON M A,HAGER W W,et al.A unified framework for the numerical solution of optimal control problems using pseudospectral methods[J].Automatica,2010,46(11):1 843-1 851.

[12] VINH N X.Flight mechanics of high-performance aircraft[M].Cambridge,England,U.K:Cambridge Univ Press,1993.

[13] TERRY H P.A common aero vehicle model,description,and employment guidance[R].Arlington:Schafer Corporation for AFRL and AFSPC,2003.