基于高斯偽譜法的天基再入飛行器滑翔軌跡優(yōu)化設(shè)計研究

周 亮,王昊宇,尚海濱,曹 莉,蘭國峰

(1. 上海機(jī)電工程研究所，上海 201109； 2. 上海航天技術(shù)研究院北京研發(fā)中心，北京 100081； 3. 北京理工大學(xué) 宇航學(xué)院，北京 100081)

0 引 言

天基再入飛行器離開天基平臺后以極高的速度再入地球大氣層，利用高升阻比氣動特性可以在臨近空間進(jìn)行無動力跳躍滑翔飛行。天基再入飛行器具有響應(yīng)速度快、機(jī)動性能好、突防能力強(qiáng)、區(qū)域覆蓋廣等諸多優(yōu)點，可實現(xiàn)遠(yuǎn)程快速精確打擊和力量投送，被認(rèn)為是一種具有廣闊應(yīng)用前景的再入飛行器[1-2]。

天基平臺運(yùn)行于近地軌道，由于地球自轉(zhuǎn)的影響，目標(biāo)點大部分時間不在軌道平面內(nèi)，此時飛行器可到達(dá)的最大橫向距離成為評估覆蓋范圍的重要指標(biāo)[3]。天基再入飛行器的橫向距離可從兩個方面進(jìn)行控制：一是在再入飛行器離開天基平臺到進(jìn)入大氣層之前的離軌段，通過施加速度切向方向的推力來改變相對軌道面的橫向距離；二是在進(jìn)入大氣層之后的無動力滑翔階段，依靠氣動力控制飛行速度方向，從而改變橫向距離。本文主要關(guān)注無動力滑翔段的再入飛行器最遠(yuǎn)橫向滑翔軌跡優(yōu)化問題。

飛行器再入滑翔階段由于超高馬赫數(shù)導(dǎo)致氣動力變化較為劇烈，同時飛行器自身動力學(xué)模型是一組時變非線性微分方程組，無法得到解析解，給再入軌跡的設(shè)計和優(yōu)化帶來很大的難度。目前軌跡優(yōu)化問題根據(jù)求解方式不同分為間接法和直接法。間接法根據(jù)龐特里亞金極大值原理，將軌跡優(yōu)化問題轉(zhuǎn)化為兩點邊值問題進(jìn)行求解，LuPing[4]、李惠峰[5]等使用間接法對高馬赫數(shù)軌跡優(yōu)化問題進(jìn)行了分析，但間接法受限于初值猜測等問題，不能很好地應(yīng)用于工程之中。直接法是將軌跡優(yōu)化問題進(jìn)行離散化，轉(zhuǎn)換為一系列的參數(shù)優(yōu)化問題進(jìn)行求解，李瑜[6]通過直接打靶法進(jìn)行了彈道優(yōu)化，涂良輝[7]使用直接配點法研究了再入軌跡優(yōu)化問題，Rao[8]使用Legendre偽譜法研究了再入飛行器軌跡控制問題，雍恩米[9]討論了高斯偽譜法在再入軌跡優(yōu)化中的應(yīng)用。這些研究表明直接法在高馬赫數(shù)再入軌跡優(yōu)化中有很大優(yōu)勢，其中偽譜法作為配點法的一種，相比于其他方法具有參數(shù)規(guī)模小、收斂性能好等優(yōu)點，被廣泛用于高馬赫數(shù)軌跡優(yōu)化問題。尚海濱[10]提出求解解析雅可比矩陣可有效提高高斯偽譜法收斂速度和精度。

目前再入飛行器多采用升力體形式進(jìn)行軌跡優(yōu)化，本文針對軸對稱再入飛行器給出了氣動力分解方式，并采用高斯偽譜法對再入軌跡進(jìn)行優(yōu)化，同時給出了計算過程中的解析雅可比矩陣的求解方式，提高了算法的收斂精度和收斂性，最終得到了滿足約束的橫向距離最大再入軌跡。

1 軌跡優(yōu)化問題的描述

高馬赫數(shù)飛行器再入大氣層后以無動力滑翔方式飛行，主要受重力和空氣動力的影響。為方便描述飛行狀態(tài)，在半速度坐標(biāo)系中進(jìn)行再入軌跡分析。該坐標(biāo)系原點在再入飛行器的質(zhì)心o，ox軸沿再入飛行器的速度方向，oy軸在速度矢量v與地心構(gòu)成的平面內(nèi)，垂直于ox軸，向上為正，oz軸垂直于xoy平面。

假設(shè)地球為勻速自轉(zhuǎn)的標(biāo)準(zhǔn)球體，則軸對稱再入飛行器在半速度坐標(biāo)系中的動力學(xué)模型可以表示為[11]

(1)

式中：r為再入飛行器質(zhì)心和地心的距離；γ和φ分別為再入飛行器的地心經(jīng)度和緯度；v為再入飛行器速度；θ為彈道傾角(再入飛行器速度方向和當(dāng)?shù)厮矫鎶A角)；σ為彈道偏角(再入飛行器速度方向與正北夾角，順時針為正)；we為地球自轉(zhuǎn)速度；Fx、Fy和Fz分別為氣動阻力、升力和側(cè)向力，控制量攻角α、側(cè)滑角β隱含在其中，計算方式為

(2)

式中：S為軸對稱飛行器特征面積；ρ為飛行器飛行高度處的大氣密度；Cx、Cy和Cz分別為氣動阻力系數(shù)、升力系數(shù)和側(cè)向力系數(shù)。

氣動力系數(shù)是馬赫數(shù)M、攻角α、側(cè)滑角β和高度h的函數(shù)，下面給出其推導(dǎo)方式。

定義飛行器本體縱軸指向頭部方向與飛行器速度方向之間的夾角η為總攻角，則阻力系數(shù)和總升力系數(shù)可以分別表示為

(3)

式中：CN和CA分別是飛行器本體坐標(biāo)系下的法向力氣動系數(shù)和軸向力氣動系數(shù)。

根據(jù)球面三角形公式可推出總攻角η與攻角α和側(cè)滑角β的關(guān)系為

cosη=cosαcosβ

(4)

由于氣動力在半速度坐標(biāo)系中的分解為

(5)

式中：D為阻力；L為升力。

因此可以得到半速度坐標(biāo)系中三軸的氣動系數(shù)分別為

(6)

CL和CD是馬赫數(shù)M和攻角α的函數(shù)，根據(jù)氣動參數(shù)表進(jìn)行二元擬合得到，具體擬合式為

(7)

實際任務(wù)中，再入飛行器的目標(biāo)點一般不在初始軌道面內(nèi)，此時需要飛行器利用氣動力做出橫向機(jī)動，橫向機(jī)動的能力代表飛行器的有效覆蓋范圍，性能指標(biāo)為[11]

minJ=-|sinδ|

(8)

式中，δ表示飛行器再入軌跡的橫程對應(yīng)的地心角，δ和橫程Rcr可分別由式(9)、式(10)求出。

sinδ=-sinσ1sinφ2cosφ1cosΔλ-cosσ1cosφ1sinΔλ+

sinσ1cosφ2sinφ1

(9)

Rcr=RE|δ|

(10)

式中：σ1、φ1是再入點的飛行器彈道偏角和緯度；φ2是終點的緯度；Δλ是終點經(jīng)度與起始點經(jīng)度之差；RE為地球平均半徑。

本文討論的是橫程最大問題，因此|δ|越大，則對應(yīng)的橫程越大。由于再入飛行器的機(jī)動能力是有限的，有|δ|<90°。在該條件下，|δ|越大，則|sinδ|越大，因此可以取公式(8)作為性能指標(biāo)。

(11)

(12)

(13)

(14)

式中：R為球頭半徑；ρs為自由流密度；vs為自由流速度；ρ0為海平面大氣密度1.225 kg/m3；v0為第一宇宙速度7 900 m/s。

對于終端點不確定的橫向飛行距離最大問題，終端約束為終端高度、速度和彈道傾角等，可表示為

(15)

根據(jù)上述分析，軸對稱再入飛行器橫向飛行距離最大的優(yōu)化問題P0可表示為

minJ=-|sinδ| s.t. (1),(11)～(15)成立

2 高斯偽譜法原理分析

采用高斯偽譜法對上述優(yōu)化模型進(jìn)行分析求解。高斯偽譜法作為直接法中的一種，是基于配點法發(fā)展而來的[11]，其核心思想是采用插值多項式基函數(shù)的實根作為配點，同時對控制變量和狀態(tài)變量進(jìn)行離散，利用離散點上的軌跡狀態(tài)值構(gòu)造拉格朗日多項式，對狀態(tài)和控制進(jìn)行全局的擬合，將連續(xù)空間的最優(yōu)控制問題轉(zhuǎn)化為非線性規(guī)劃問題，從而方便使用數(shù)值方法求解。

2.1 時域變換

采用勒讓德-高斯-洛巴托(Legendre-Gauss-Lobatto)點，簡稱LGL點，對連續(xù)的優(yōu)化問題進(jìn)行離散化，離散區(qū)間為[-1,1]。原動力學(xué)中狀態(tài)變量和控制變量的時域為[t0,tf]，需要進(jìn)行時域變換，將原優(yōu)化問題的時間區(qū)間從t∈[t0,tf]轉(zhuǎn)換到τ∈[-1,1]，建立仿射關(guān)系

(17)

可將動力學(xué)方程(1)表示為簡化形式

則原動力學(xué)方程通過時域變換后的形式為

(18)

2.2 插值計算

定義集合Γ為N+1個LGL點，且集合中的 LGL點嚴(yán)格單調(diào)遞增，以LGL點為插值節(jié)點，分別構(gòu)建N+1個拉格朗日插值多項式對狀態(tài)變量和控制變量進(jìn)行插值擬合，如式(19)～(20)所示。

(19)

(20)

式中：x(τi)和u(τi)分別為τi點對應(yīng)的真實狀態(tài)變量和控制變量；Li(τ)為N次拉格朗日插值基函數(shù)，表示為

(21)

式中：PN為N次勒讓德多項式。

顯然Li(τk)具有如下性質(zhì)

(22)

代入式(19)～(20)中，可得到在LGL插值點上估計值是精確等于真實值的，即

X(τi)=x(τi)

(23)

U(τi)=u(τi)

(24)

拉格朗日插值擬合的狀態(tài)方程(19)的微分形式為

(25)

(26)

式中，k=0,1,2,…，N。

根據(jù)式(18)和式(25)構(gòu)建微分偏差方程

(27)

原優(yōu)化問題中的動力學(xué)約束可轉(zhuǎn)為ξk=0的約束，即只與LGL節(jié)點處狀態(tài)變量和控制變量有關(guān)的約束。同理原優(yōu)化問題中的過程約束(11)～(14)也可按照上述方法進(jìn)行插值并轉(zhuǎn)化為與LGL節(jié)點處狀態(tài)變量和控制變量有關(guān)的約束[12]，定義為

gk(x(τk),u(τk),τk)≤0

(28)

根據(jù)式(6)～(7)可看出，性能指標(biāo)只與初始狀態(tài)變量和末端狀態(tài)變量有關(guān)，已為離散狀態(tài)。

至此原優(yōu)化問題P0轉(zhuǎn)換為非線性規(guī)劃問題P1，優(yōu)化變量為每個離散點的狀態(tài)變量x(τi)、控制變量u(τi)，可以表示為

(29)

2.3 解析雅可比矩陣

求解上述非線性規(guī)劃問題時，需要多次對性能指標(biāo)和約束雅可比矩陣進(jìn)行求解，因此提高雅可比計算精度可以有效地提升非線性規(guī)劃問題的收斂性[13]，在此給出解析雅可比矩陣的推導(dǎo)過程。

由式(9)可以看出，性能指標(biāo)只與始末狀態(tài)的經(jīng)緯度和彈道偏角有關(guān)，因此性能指標(biāo)對優(yōu)化變量的雅可比矩陣可表示為

(30)

根據(jù)式(27)，ξk與每個LGL點的優(yōu)化變量都有關(guān)，以第k點的約束為例，雅可比矩陣可表示為

(31)

式中：Fk為第k點的狀態(tài)和控制變量。

對于式(28)中的過程約束，由原約束式(11)～(14)可以看出，各個約束量離散后都只與當(dāng)前離散點的優(yōu)化變量有關(guān)，約束式(15)只與最后一個離散點有關(guān)，同理可以求出解析雅可比矩陣，此處不再列舉具體形式。

3 仿真與分析

表1 初始狀態(tài)Tab. 1 Initial state

表2 末端狀態(tài)約束Tab. 2 Terminal state restriction

圖1 高度隨時間變化曲線Fig.1 Time history of height

圖2 速度隨時間變化曲線Fig.2 Time history of velocity

圖3 彈道傾角隨時間變化曲線Fig.3 Time history of trajectory oblique angle

圖1～3分別為飛行器再入過程中高度r、速度v以及彈道傾角θ隨時間變化的仿真曲線。由圖1可以看出軸對稱飛行器再入過程呈現(xiàn)出跳躍式軌跡；結(jié)合圖2可以發(fā)現(xiàn)再入初始階段飛行器速度較快，飛行器高度變化較為劇烈；隨著速度的降低，根據(jù)圖3可以看出飛行器進(jìn)入相對平穩(wěn)滑翔階段，彈道傾角保持在[-1°,1°]之間波動，飛行的最后階段由于要滿足末端點的狀態(tài)約束，飛行器在滑行飛行末段進(jìn)行了姿態(tài)調(diào)整，彈道傾角波動范圍增大，末端點處為-2°，滿足所給出的約束。

3種不同初始狀態(tài)下的仿真結(jié)果對比如表3所示，可以看出再入時間與再入速度呈正相關(guān)，即再入速度越大，滿足約束條件的再入時間越長，同時優(yōu)化性能指標(biāo)橫向距離也越大。

表3 仿真結(jié)果Tab.3 Results of simulation

圖4～5為控制變量隨時間變化曲線，初始階段彈道傾角為負(fù)數(shù)，因此采用大攻角再入，快速調(diào)整彈道傾角使熱流密度峰值滿足約束要求，之后滑翔階段攻角緩慢增加，符合約束要求并且利于工程實現(xiàn)。初始階段以較大側(cè)滑角再入可以快速改變彈道偏角，使飛行器速度轉(zhuǎn)向最大橫向距離方向，之后側(cè)滑角逐漸減小到0°，利用氣動力進(jìn)行速度方向調(diào)整，使得飛行器再入的橫向距離增大。圖6～8分別給出了熱流密度、過載和動壓3個主要路徑約束隨時間的變化曲線。

圖4 攻角隨時間變化曲線Fig.4 Time history of attack angle

圖5 側(cè)滑角隨時間變化曲線Fig.5 Time history of sideslip angle

圖6 熱流密度隨時間變化Fig.6 Time history of heat flux

圖7 過載隨時間變化Fig.7 Time history of overload

圖8 動壓隨時間變化Fig.8 Time history of kinetic pressure

可以看出，3種再入狀態(tài)下的路徑約束都滿足了設(shè)定的約束條件，說明高斯偽譜法均能得到滿足約束并且使橫向距離最大的最優(yōu)控制曲線。

3個路徑約束只有熱流密度達(dá)到了約束邊界，說明影響飛行過程的最主要因素為熱流密度峰值約束。由熱流密度約束式(11)可以看出熱流密度主要受大氣密度和飛行器速度影響，初始階段飛行器高度快速下降導(dǎo)致大氣密度增加，而飛行器速度下降并不明顯，導(dǎo)致熱流密度在此處達(dá)到峰值，因此飛行器初始階段飛行條件較為惡劣。此階段對飛行器控制系統(tǒng)、隔熱系統(tǒng)等有較高要求，這也符合工程應(yīng)用實際，從而進(jìn)一步說明此方法得到的最優(yōu)解接近實際情況。

綜上分析，對于給定的3種初始狀態(tài)，高斯偽譜法均能得到滿足約束并且使橫向距離最大的最優(yōu)控制曲線，說明了算法的有效性和魯棒性。

4 結(jié)束語

本文以天基再入飛行器軌跡優(yōu)化設(shè)計為背景，對軸對稱高馬赫數(shù)飛行器再入問題進(jìn)行建模分析，提出用高斯偽譜法解決此類問題并對算法進(jìn)行詳細(xì)說明。對不同的再入狀態(tài)進(jìn)行仿真計算，得到了較為平滑的控制曲線，并且滿足各種約束條件，實現(xiàn)了軸對稱高馬赫數(shù)飛行器再入軌跡優(yōu)化問題的求解。因此，本方法對于高馬赫數(shù)飛行器再入軌跡的設(shè)計和優(yōu)化有一定的工程意義。