固體助推飛行器大氣層內閉環(huán)制導方法

郭 杰，宋滿金，柳 青，唐勝景

(1. 北京理工大學宇航學院，北京 100081；2. 北京空天技術研究所，北京 100074)

0 引 言

隨著助推-吸氣式高超聲速飛行器朝著多任務、多約束的方向發(fā)展[1-2],面對復雜終端約束的大氣層內上升段精確制導是目前研究的重點。為了保障助推器與載荷分離時的交班點約束,飛行器需要在上升過程中下壓彈道飛行,并在彈道終端實現(xiàn)轉平,該低彈道轉平特性會導致較大的軸向過載。此外,為了滿足載荷后續(xù)飛行條件,上升段終端攻角和動壓需保持在一定約束范圍內。而大氣層內上升段存在高度較低、氣動擾動較大、動壓變化劇烈等問題,極易導致軌跡偏移、抖動甚至發(fā)散,給制導帶來巨大挑戰(zhàn)。

針對運載器上升段制導問題,學者們已經進行了相關研究。攝動制導[3]是一類被廣泛應用的標稱軌跡制導方法,通過將實際彈道在標準彈道附近泰勒展開,控制目標量與理論值的偏差為零,實現(xiàn)約束滿足,但該方法依賴標稱軌跡,適應能力和魯棒性較差[4]。為了提高自主性和制導精度,迭代制導、凸優(yōu)化制導和能量管理等自適應制導方法陸續(xù)被提出。呂新廣等[5]、韓雪穎等[6]提出了能夠約束終端入軌姿態(tài)的迭代制導方法;王洪波等[7]基于改進的序列凸優(yōu)化設計了固體火箭入軌制導方法,有較高的求解效率;Patha等[8]提出一種交變姿態(tài)控制的能量管理方法,陳喆等[9]、陳思遠等[10]對能量管理方法進行完善,實現(xiàn)了終端姿態(tài)約束;張遷等[11-12]在傳統(tǒng)速度管控模型的基礎上引入定點制導算法,有效提高了入軌精度。上述自適應算法能夠根據飛行器當前狀態(tài)在線計算制導指令,提高運載器的自主飛行能力。然而,上述方法在進行建模時或不考慮氣動力影響,或在推導過程中對氣動力進行簡化,大多適用于高海拔大氣稀薄的飛行階段,在低海拔氣動不確定情況下會有較大偏差。

對于大氣層內上升段制導方法,除已經實際應用的標稱軌跡跟蹤制導外,國內外學者主要集中在基于最優(yōu)控制理論的運載器大氣層內軌跡在線規(guī)劃與閉環(huán)制導的研究上,其求解思路分為間接法與直接法。間接法將最優(yōu)軌跡上升問題轉化為兩點邊值問題進行求解[13-15],但存在推導繁瑣、收斂半徑小和對初值敏感等不足,且應用在運載器指數(shù)攻角轉彎后的飛行階段。直接法將最優(yōu)軌跡規(guī)劃問題轉化為有限維參數(shù)優(yōu)化問題,離散狀態(tài)量或控制變量,利用數(shù)值計算方法進行軌跡規(guī)劃,避免了大量的過程推導。崔乃剛等[16]采用兩步梯度投影算法(MGPA)非線性規(guī)劃方法進行在線軌跡規(guī)劃,并應用基于線性二次型調節(jié)器(LQR)的閉環(huán)制導方法實現(xiàn)軌跡跟蹤;賀前偉等[17]提出一種基于模型預測靜態(tài)規(guī)劃算法的自適應制導方法,實現(xiàn)不確定情況下的末段狀態(tài)高精度制導。隨著凸優(yōu)化理論日漸完善,其也被應用于大氣內上升段制導,郝澤明等[18]針對固體火箭上升段實時軌跡優(yōu)化問題,提出一種基于鄰近-牛頓-康托維奇凸規(guī)劃的軌跡優(yōu)化方法,用內點法求解二階錐規(guī)劃問題,保障軌跡優(yōu)化實時性。在智能方法方面,何潤林[19]用粒子群優(yōu)化(PSO)算法對射角、全程攻角進行優(yōu)化,并設計了自適應動態(tài)逆制導策略;周宏宇等[20]提出一種全新的飛行剖面,使用改進PSO求解剖面參數(shù),并引入強化學習機制提高求解效率。上述基于數(shù)值求解的優(yōu)化算法雖能夠求得最優(yōu)或次優(yōu)軌跡,但一般求解參數(shù)較多,且存在求解精度與計算效率的矛盾,在線應用困難。

綜上所述,對于上升段制導問題已有較多研究,但存在以下不足:部分方法可以對終端姿態(tài)進行約束或能夠在線應用,但是不適用于稠密大氣層內制導;或者能夠通過優(yōu)化的方法求解大氣層內上升段最優(yōu)軌跡或制導指令,但在線應用困難。本文針對固體助推器大氣層內上升段制導問題,考慮推力不可調和耗盡關機特性,重點研究了制導指令設計與終端姿態(tài)約束問題,提出一種閉環(huán)制導方法。設計了解析形式的俯仰角剖面,能夠更新剖面約束、精確控制終端攻角,且只需通過數(shù)值方法校正兩個剖面參數(shù)即可更新整個剖面并直接獲得控制量,有效提高在線計算效率;采用在線更新策略,根據飛行器當前狀態(tài)與目標位置實時更新下一制導時刻俯仰角指令,能夠在氣動參數(shù)不確定情況下有較好的魯棒性。與傳統(tǒng)方法相比,該閉環(huán)制導方法在大氣層內具有較高的制導精度,且具有工程應用潛力。

本文結構如下:首先,設計了一種多項式形式的俯仰角剖面,通過合理地選擇剖面終端值保證了終端攻角約束。然后,以數(shù)值算法校正俯仰角剖面參數(shù)保證高度和彈道傾角約束,同時設計剖面在線更新流程。最后,進行仿真分析驗證該方法的制導精度與魯棒性。

1 上升段數(shù)學模型

1.1 運動模型

本文研究固體助推器大氣層內上升段制導問題,整個上升過程高度較低、耗時較短,可忽略地球自轉,基于“瞬時平衡”假設,只考慮縱向平面運動,得到飛行器在速度坐標系下的質心運動方程為:

(1)

式中:r為位置;G為重力;T為推力;L為升力;D為阻力;m為質量;mc為秒耗量。

運動方程的控制量為攻角α。重力、推力、升力和阻力可以寫為飛行時間、位置、速度和攻角α的函數(shù)。其中,重力表達式為

(2)

式中:μ為地球引力常數(shù);R為地心到初始發(fā)射位置的矢量。

推力表達式為

T=T(t)Γ(α)Iv

(3)

式中:T為推力幅值;Γ(x)為旋轉矩陣;Iv為速度的單位矢量,表達式為

(4)

式中:v為速度矢量。Γ(x)表達式如下:

(5)

升力和阻力表達式為

(6)

式中:L,D分別為升力幅值和阻力幅值,表達式如下:

(7)

式中:ρ為大氣密度;V為速度幅值;S為飛行器參考面積;Ma為馬赫數(shù);CL,CD分別為升力系數(shù)和阻力系數(shù),是馬赫數(shù)和攻角的函數(shù),表達式如下:

(8)

1.2 約束條件

大氣層內上升段氣動環(huán)境復雜,大氣密度隨高度變化范圍大,對流層和同溫層以內大氣密度變化超過40倍;且助推段馬赫數(shù)變化劇烈,經歷亞聲速、跨聲速、超聲速甚至高超聲速等階段,為了保證飛行器飛行安全與結構完整性,飛行器大氣層內上升段需要滿足一定的過程約束。此外,飛行器在助推器關機時需要滿足交班點條件等任務要求,因此還需滿足一定的終端約束。

(1)過程約束

考慮飛行器結構強度和安全等因素,飛行過程中需要考慮過載約束:

(9)

式中:nx,ny分別為彈體坐標系下的軸向過載和法向過載;G為飛行器重量。

(2)終端約束

從飛行器的飛行任務考慮,助推段終端需要滿足高度約束和彈道傾角約束。同時為了滿足載荷后續(xù)工作條件,還需要終端攻角、馬赫數(shù)-動壓保持在一定窗口約束范圍內,綜上可得終端約束如下:

(10)

式中:h為飛行器高度;θ為彈道傾角;tf為終端時間;q為動壓,qfmin和qfmax分別為動壓窗口的下界與上界。

(3)約束轉化

助推器上升過程中,法向過載大小決定著細長軸體所受彎矩大小,過大會導致軸體斷裂,造成嚴重后果,因此法向過載是重點關注的約束量。由式(7)和式(9)可得:

(11)

定義CLcosα+CDsinα為“過載系數(shù)”,式(11)可寫為

(12)

上升段為一壓低彈道的過程,在大動壓段攻角基本為負值,且過載系數(shù)是速度和攻角的函數(shù),其變化規(guī)律如圖1所示,表明在一定范圍內,過載系數(shù)與速度和攻角呈現(xiàn)單調函數(shù)關系。在速度一定時,過載系數(shù)絕對值與攻角絕對值正相關,因此可通過限制攻角大小來滿足法向過載約束。

圖1 過載系數(shù)與速度和攻角的關系

在飛行器實際飛行過程中,已知當前高度H和速度V,可通過求解式(13)得到滿足約束的當前最小攻角值。

(13)

若當前α不滿足幅值約束,則通過式(14)進行限幅處理

(14)

式中:αt為實際攻角值;αmin為當前滿足約束的攻角最小值。

2 上升段閉環(huán)制導方法

2.1 控制量選擇

由固體助推飛行器的動力學方程可知,攻角和俯仰角都可作為飛行控制量,因此需要從數(shù)值計算的角度分析攻角和俯仰角對飛行軌跡的影響能力,確定剖面選擇[21]。

圖2為攻角改變時飛行器體軸以及軌跡改變情況示意圖,攻角從α1改變到α2,彈道傾角從θ1改變到θ2。除飛行器體軸改變造成的正增量Δα2外,還有速度方向改變產生的攻角負增量Δα1;圖3為俯仰角改變時飛行器體軸以及軌跡改變情況示意圖,俯仰角從φ1改變到φ2,彈道傾角同樣從θ1改變到θ2,但是過程中只有體軸變化產生的俯仰角正增量Δφ。

圖2 攻角改變時軌跡變化

圖3 俯仰角改變時軌跡變化

2.2 飛行剖面設計

基于上節(jié)分析,對于大氣層內上升段制導,考慮大氣密度分布特點和飛行器速度規(guī)律,本文設計了如下分段多項式形式的俯仰角剖面:

(15)

式中:ki(i=1,2,3,4)分別為3次多項式系數(shù);t1為垂直上升段結束時間;tf為飛行總時間,剖面形狀如圖4所示。確定式(15)中的3次多項式系數(shù)需要4組(t,φ)形式的約束條件,其求解方法如下。

圖4 俯仰角剖面

初始上升段中,t1可通過人為選擇給定,考慮到飛行器發(fā)射高度低,大氣密度大,初始飛行階段處于大氣稠密區(qū),氣動力和推力之間存在較強耦合,選擇合適的t1可以使飛行器在跨聲速之后進行轉彎,避免受到較大的氣動載荷和氣動擾動,降低載荷、控制系統(tǒng)設計難度。根據上升段姿態(tài)角關系φ=α+θ,可通過終端彈道傾角約束與終端攻角約束得到終端俯仰角大小:

(16)

設θ′f為上一制導周期預測得到的終端彈道傾角值,則θ′f會隨著制導進程逐漸接近實際彈道傾角終端值θf,因此可在第一個制導周期內采用式(16)計算剖面終端值,后續(xù)采用下式更新剖面終端值,實現(xiàn)Δθf的誤差消除:

φ(tf)=θ′f+(αfmin+αfmax)/2

(17)

綜上,可得到兩組約束條件:(t1,π/2)、(tf,φ(tf))。

另外兩組約束條件通過設計剖面參數(shù)獲得。在上述分析的基礎上,增加(t2,φ2)、(t3,φ3)兩個控制點。其中t2,t3為人為給定的時間節(jié)點,且滿足t1≤t2≤t3≤tf;φ2,φ3為2個待設計剖面參數(shù)。一旦確定φ2與φ3,則可得到如下方程組:

(18)

求解式(18)即可得到式(15)中的3次多項式系數(shù),進而確定整個俯仰角剖面。該剖面分為兩段,其中第一段常值俯仰角用于將飛行器渡過跨聲速階段;第二段多項式形式俯仰角充分利用飛行器性能,進行轉彎以及調整飛行狀態(tài)來滿足終端約束。

對于t2,t3的選擇,理論上只要滿足t2≠t3的條件,就能夠求解式(18)得到完整的俯仰角剖面。但若|t2-t3|在數(shù)值上過小,求解過程中可能會出現(xiàn)矩陣奇異現(xiàn)象,導致無解;同理,|t1-t2|和|t3-tf|的值也不能過小。因此,本文選擇t2,t3的策略如下:

(19)

2.3 剖面參數(shù)更新策略

傳統(tǒng)基于在線軌跡規(guī)劃的閉環(huán)制導思路為:在飛行過程中根據當前狀態(tài)和終端約束規(guī)劃軌跡,將軌跡設計由離線過程變?yōu)樵诰€過程,并設計制導律對軌跡進行跟蹤[22]。這類方法在實際應用中需要合理選擇軌跡規(guī)劃周期與制導周期,才能兼顧計算效率與制導精度。本文所設計的俯仰角剖面為解析形式,且完全由2個剖面參數(shù)(φ2,φ3)確定,并且實際飛行中飛行器可直接根據剖面得到制導指令,即將閉環(huán)制導問題轉化為剖面參數(shù)(φ2,φ3)的搜索問題,大幅提高求解效率。

在每個制導周期內,給定初始剖面參數(shù)(φ2,φ3),求解得到完整的俯仰角剖面;以當前狀態(tài)為初始狀態(tài),積分式(1)得到終端狀態(tài),從而得到終端高度h(tf)和終端彈道傾角θ(tf)。由于(φ2,φ3)決定了整個俯仰角剖面,也即決定了終端高度與終端彈道傾角,因此可以利用二元非線性方程組求根的牛頓迭代法搜索得到滿足終端約束的剖面參數(shù)。

在當前時刻與終端時刻確定的情況下,終端高度和終端彈道傾角是兩個剖面參數(shù)(φ2,φ3)的函數(shù):

(20)

本文在每個制導周期內,以高度、彈道傾角雙約束,基于俯仰角剖面同時對兩個剖面參數(shù)進行校正:

(21)

式中:hf為期望終端高度,θf為期望終端彈道傾角,二者可由終端約束得到。上述兩個變量在每個制導周期內均為已知量,則方程組式(21)可化為關于俯仰角參數(shù)(φ2,φ3)的二元非線性方程組:

(22)

式中:F(φ2,φ3)=h(tf)-hf,G(φ2,φ3)=θ(tf)-θf。利用牛頓迭代法可快速求解上式,且能得到收斂解[23-24]：

(23)

式中:F,G分別為式(22)中F(φ2,φ3),G(φ2,φ3)函數(shù);φ2,i和φ3,i為第i次迭代時的φ2,φ3值;F′φ2,F′φ3,G′φ2,G′φ3為F,G函數(shù)分別對φ2,φ3的偏導數(shù)。

考慮到非線性函數(shù)的復雜性,無法直接求得準確的偏導數(shù)表達式,本文采用有限差分法近似得到偏導數(shù)值,且精度滿足如下要求。

(24)

式中:Δ為一大于0的小量。實際計算過程中,通過式(23)得到下一次迭代的φ2,φ3值,重復迭代直至F<ε1且G<ε2,即可得到滿足精度要求的φ2,φ3解。ε1,ε2為可接受的終端高度與終端彈道傾角誤差。若在飛行過程中因參數(shù)擾動過大,或因控制量限幅導致后續(xù)制導周期得到不滿足精度的解,此時取目標函數(shù):

(25)

利用參數(shù)優(yōu)化的牛頓迭代法可快速求解使目標函數(shù)式(25)最小的剖面參數(shù)[25]：

(26)

(27)

式中:a1,a2,a3,b1和b2為濾波器待定參數(shù)。式(27)的時域表達式如下:

(28)

式中:x(t)為輸入信號;y(t)為輸出信號。將上式轉換成對應的離散域差分方程,并將時域離散信號轉變?yōu)樾蛄行盘柨傻?

y(n)=λ1x(n)+λ2x(n-1)+λ3x(n-2)+

λ4y(n-1)+λ5y(n-2)

(29)

式中:λi(i=1,2,3,4,5)為相應的系數(shù),表達式如下:

(30)

式中:Ts為采樣周期。使用式(29)對所求得的俯仰角進行濾波,最終可以得到平滑的制導指令。整個基于俯仰角剖面更新的閉環(huán)制導算法總體流程如圖5所示。

圖5 閉環(huán)制導方法流程

在制導過程中,為了提高預測環(huán)節(jié)計算效率,減少對箭上計算能力的占用,前期可采用定步數(shù)積分,積分步長Δh取(tf-t)/N,t為當前時刻,N為積分總步數(shù),可取200;當后期Δh減小到與制導周期一致時,保持不變。另外,為減少剖面更新迭代次數(shù),可取上一制導周期剖面參數(shù)作為迭代初值,提高搜索效率;還可以在上升段前期多個制導周期更新一次俯仰角剖面。隨著飛行時間增加,剖面更新耗時減少,可減少為每個制導周期內更新一次剖面。本文采取剖面更新頻率如下:

(31)

上式表示κ個制導周期更新一次剖面參數(shù)。

上述基于俯仰角剖面參數(shù)更新的閉環(huán)制導核心思路是將牛頓迭代與剖面設計相結合,將每個制導周期內制導指令求解轉化為終端高度、終端彈道傾角雙約束的雙剖面參數(shù)搜索問題。這種思路還可以適用于離線軌跡規(guī)劃,既可基于攻角-時間剖面,也可基于俯仰角-高度剖面,甚至基于俯仰角-速度剖面進行設計。剖面參數(shù)也可以選擇任意兩個同時影響整個剖面形狀的參數(shù)。

3 仿真校驗

以固體火箭助推飛行器為研究對象,地球半徑取6 371 km;地球引力常數(shù)μ取3.986 005×1014m3/s2;g為重力加速度;聲速模型采用美國標準大氣1976,在50 km以下具有較高的擬合精度。

飛行器總飛行時間為90 s;垂直上升段時間t1取20 s;初始質量m0為4 150 kg;秒耗量mc為30 kg/s;推力幅值T為95 kN;參考面積S為0.8 m2;制導周期取0.1 s;初始剖面參數(shù)的迭代初值由線性插值得到,φ2取60°,φ3取30°;濾波器參數(shù)a1,a2,a3,b1和b2分別取0.1, 0.2, 1, 1, 1.4,采樣周期Ts取0.1 s。

以飛行器出發(fā)射筒時狀態(tài)作為初始狀態(tài),如表1所示;發(fā)動機關機時飛行器需處于轉平狀態(tài),終端狀態(tài)約束如表2所示。過載約束為:軸向過載最大值nxmax為3;法向過載最大值nymax為2。

表1 初始狀態(tài)

表2 終端狀態(tài)

設置3組仿真條件:第1組為標準條件;第2組為大氣密度、升力系數(shù)和阻力系數(shù)各拉偏-10%,推力拉偏-5%,質量秒耗量拉偏-2%;第3組為大氣密度、升力系數(shù)和阻力系數(shù)各拉偏+10%,推力拉偏+5%,質量秒耗量拉偏+2%。

本文數(shù)值計算在Intel i5 3.10 GHz臺式機電腦下進行,程序均在MATLAB 2020b環(huán)境下編譯運行。

3.1 濾波效果分析

為驗證本文所提出的閉環(huán)制導方法中的濾波效果,分別進行有濾波器和無濾波器的對比仿真。第3組條件下仿真曲線如圖6所示,所有組仿真結果見表3。

表3 濾波前后仿真結果

圖6 濾波前后對比仿真曲線

由表3可知,濾波器使用前后,飛行器終端狀態(tài)接近且均具有較高精度,而由圖6中(a)、(b)可知,濾波器使用前,俯仰角和攻角出現(xiàn)震蕩現(xiàn)象,這是因為在干擾條件下,為保證終端精度,相鄰兩次剖面更新所得的俯仰角剖面存在差異,導致時序制導指令產生波動。加入濾波器后,俯仰角和攻角較為平滑,符合實際飛行需求。

3.2 對比仿真分析

為驗證本文所提方法的有效性,選取基于導引系數(shù)快速求取的攝動制導[3]、基于PD的最優(yōu)軌跡跟蹤[26]和本文閉環(huán)制導3種制導方法進行仿真對比。其中最優(yōu)軌跡以終端速度最大作為目標函數(shù),通過Gauss偽譜+序列二次規(guī)劃(SQP)方法在標準條件下離線規(guī)劃。第3組條件下仿真曲線如圖7所示,所有組仿真結果見表4。

表4 不同仿真條件下的仿真結果

圖7 第3組拉偏情況下仿真曲線

由表4中3組仿真結果的高度誤差可知,3種制導方法在標準狀況下均有較高精度,但在干擾較大的情況下,基于導引系數(shù)快速求取的攝動制導高度誤差接近80 m,說明簡化氣動力的該方法不適用于大氣層內上升段制導。基于PD的最優(yōu)軌跡跟蹤和閉環(huán)制導方法均能較好地約束終端高度和終端彈道傾角,其中閉環(huán)制導方法的終端高度誤差在10 m以內,終端彈道傾角誤差不超過0.1°,相比于PD跟蹤均提高了一個數(shù)量級。這是因為本文提出的閉環(huán)制導方法,能夠在每一個制導周期內,通過搜索雙剖面參數(shù)實現(xiàn)對終端高度和彈道傾角的精確控制,魯棒性較強。由馬赫數(shù)和動壓值仿真結果可知,飛行全程動壓值較小,3種制導方法均滿足終端窗口約束;3組終端攻角也均滿足窗口約束,且在標準條件下攻角波動均不超過0.1°。在有干擾的情況下,閉環(huán)制導的終端攻角波動比PD跟蹤更小,說明本文提出的更新剖面終端值方法能更好地約束終端攻角。

由圖7(a)(d)可知,標稱最優(yōu)軌跡與3種制導方法所得軌跡平滑、變化平穩(wěn);圖7(b)(e)(f)表明所有軌跡均滿足過程約束;圖7(c)表明4種方法得到的攻角曲線連續(xù)且沒有抖振出現(xiàn),其中基于PD的最優(yōu)軌跡跟蹤攻角曲線變化趨勢與標稱最優(yōu)攻角接近,且與閉環(huán)制導攻角曲線差異較大。這是因為PD跟蹤依賴標稱軌跡,在干擾存在時調整能力有限,而閉環(huán)制導自主性較強,能夠根據狀態(tài)實時調整控制量,充分發(fā)揮飛行器飛行能力;俯仰角剖面平均更新時間為28 ms,證明該方法具有在線應用價值。

3.3 制導精度與魯棒性分析

針對氣動系數(shù)偏差、大氣密度偏差、推力偏差和質量秒耗量偏差進行200次蒙特卡洛仿真。假設各項參數(shù)擾動符合正態(tài)分布,擾動參數(shù)如表5所示。仿真結果如圖8所示,終端約束平均值見表6。

表5 擾動參數(shù)

表6 終端約束蒙特卡洛仿真平均值

由圖8(a)(d)可知,高度和彈道傾角曲線平滑、軌跡無震蕩;圖8(b)(e)(f)表明在氣動不確定情況下,所有仿真均滿足過程約束,且動壓、過載變化平緩,終端動壓均在窗口范圍內;圖8(c)表明所有仿真攻角曲線無抖振出現(xiàn),終端攻角均滿足窗口約束;圖8(g)(h)表明200組仿真終端高度誤差均不超過10 m、終端彈道傾角誤差均不超過0.15°,具有較高精度;圖8(i)為終端攻角分布,各種隨機擾動下,所有終端攻角的最大波動值不超過0.3°,證明本文所提方法能夠較好地約束終端攻角,且具有較強的魯棒性。

4 結 論

本文將具有復雜終端約束的固體火箭助推飛行器大氣層內上升段制導問題,通過設計制導指令與動態(tài)更新策略,轉化為參數(shù)實時搜索問題,提出一種閉環(huán)制導方法。理論分析和仿真結果表明:

1) 閉環(huán)制導方法能夠在線生成控制量,且控制量曲線平滑,滿足飛行器上升段飛行任務需求,具有在線應用價值。

2) 方法對終端攻角有較強的約束能力,能夠為飛行器后續(xù)工作創(chuàng)造良好環(huán)境。

3) 所設計的方法在擾動條件下具有一定的精度和魯棒性。