基于橢球近似解析法的飛行器再入軌跡優化

張 凱 楊小龍 鐘 震

北京宇航系統工程研究所，北京，100076

0 引言

為使飛行器跨越復雜的大氣環境安全返回，滿足多種約束的飛行器再入軌跡規劃一直是研究熱點[1]。按照是否跟蹤標準軌跡，飛行器再入過程主要分為標準軌跡制導和預測校正制導兩大類[2]。標準軌跡制導形式簡潔跟蹤穩定有較高應用價值[3]，如何設計標準軌跡是該方法的重要研究方向[4]。經典的阻力加速度設計方法在阻力加速度剖面內得到再入走廊，引入航程與阻力加速度的關系，將航程轉化成再入走廊內分段阻力加速度曲線，大大簡化了軌跡的設計流程，最早在航天飛機再入中獲得成功應用[5]。文獻[6]討論了阻力加速度方法在低升阻比火星再入飛行器的應用，文獻[7]將過程約束轉化到再入走廊中，通過規劃滿足約束的阻力加速度速度剖面生成航天飛機參考軌跡，文獻[8]將航天飛機的二維再入軌跡規劃擴展到三維情況，采用基于降階的運動方程和最優控制理論實現縱橫向參考加速度的在線規劃，文獻[9]改進了傳統阻力加速度剖面，采用兩段折線與走廊邊界快速設計標準軌跡，提高了標準軌跡制導方法的實時性，文獻[10]將阻力加速度剖面與預測校正方法結合，由走廊上下邊界插值獲得標準軌跡，插值參數由預測校正方法確定，并在三維空間對標準軌跡進行校正，提升了再入飛行能力，文獻[11]考慮了速度傾角約束，通過在阻力加速度-能量剖面終端加入速度傾角控制滿足了速度傾角約束。

不同于現有文獻直接采用阻力加速度的軌跡規劃方法，本文推導了根據阻力加速度曲線快速計算再入段軌跡所需控制能力的公式，分析了再入段航跡傾角與阻力加速度曲線的關系，可以根據再入航跡傾角約束規劃相應的曲線，保證飛行器在再入初始段和結束段滿足相應約束，提出了將原軌跡阻力加速度剖面線性化的方法，并在此基礎上提出了一種基于橢球近似解析的優化算法，根據該算法可以將標稱軌跡修正為與其控制能力匹配的軌跡，通過仿真驗證了算法的有效性。

1 飛行器再入模型

設飛行器再入過程中地球為圓球，飛行器為剛體，主要受氣動升力、氣動阻力和引力作用，則飛行器三維的運動學、動力學方程描述如下[12]

(1)

其中r是飛行器地心距，r=R0+h，h是飛行器當前高度，R0為地球平均半徑，ωe為地球自轉角速度，θ、φ為經緯度，V為速度，γ為航跡傾角，ψ為航跡偏角，L、D為升力、阻力加速度，σ為傾側角。

升力和阻力加速度的表達式為

(2)

整理可得

(4)

其中

(5)

設R為再入段航程，則有

(6)

根據上式可以通過V-D曲線、速度-航跡傾角曲線來預測航程。

2 飛行器控制力與再入傾角

(7)

其中

(8)

飛行器不同的再入傾角所引起的阻力加速度變化不同，同時，有些飛行器再入段結束時對航跡傾角有要求，如果采用單獨的再入段航跡傾角調節控制，會增加軌跡設計的復雜性，因此需要找出航跡傾角與規劃的V-D曲線之間的關系，保證所規劃的軌跡滿足航跡傾角約束。根據前文推導，將軌跡傾角表達式移項整理可進一步得到

(9)

其中等式右端所有的量都可知，因此當得到初始或終端航跡傾角時，就可以計算出初始、終端的V-D曲線的斜率。記

(10)

使規劃曲線的初始段和結束段滿足上述公式得到的斜率即可保證航跡傾角約束。

3 基于標稱軌跡的線性化模型

再入段標稱軌跡描述采用曲線-直線-曲線的三段式曲線，形式如下：

圖1 阻力加速度曲線示意圖

其具體規劃方法本文不再詳述。通過該方法規劃的初始軌跡所需控制能力可能大大超過了飛行器的容許控制能力，傳統方法在這種情況下需要不斷修正V-D曲線，且不能保證所修正的曲線滿足控制能力約束，理想的方法是修正該曲線使其所需控制能力在飛行器的控制能力極限之內。

線性偽譜法中將狀態方程在標稱軌跡附近線性化[13]，根據此思路，本文以標稱V-D曲線為基礎，考慮將控制力和V-D曲線的對應關系在小范圍內進行線性化。求u對D的偏導數，得到V-D曲線附近Δu和ΔD之間的線性關系。將V-D曲線用一系列等距的V對應D的散點表示，V-u曲線同理。其求導運算(D對V的導數)用差分計算來替代，求一階導數時，需要前后相鄰的兩個量，求二階導數需要前后相鄰的5個量。由于計算u時最多需要D對V的二階導數，因此只需要計算u對其前后5個D值的偏導數即可，這5個偏導數的計算采用Matlab軟件符號運算進行推導，不再展開結果。由于改變D導致的u的改變量可以表示如下

(11)

將V-D曲線分成等距的N個點，將全N個點都考慮進去，可得

(12)

將上式簡記為A1x1=x2。其中

x1=[ΔD1…ΔDN]T，x2=[Δu1…ΔuN]T。

考慮增加航程約束、高度約束、航跡傾角約束如下，

(13)

分別簡記為A2x1=b2、A3x1=b3、A4x1=b4。綜合可得

(14)

式(14)記為Ax=b。變量x1、x2中元素為ΔDi、Δui，需滿足如下約束

Dmin-Di≤ΔDi≤Dmax-Di
umin-ui≤Δui≤umax-ui

(15)

其中Di、ui為原規劃曲線的阻力加速度、控制量，Dmax、Dmin為再入走廊的阻力加速上下限，umax、umin為控制能力的上下限，其含義為修正軌跡后保證新軌跡不超過再入走廊的邊界，同時所需控制力不超過飛行器的控制能力極限。其中Dmax、Dmin、umax和umin也可根據實際情況進行調整。

4 基于橢球近似的解析優化算法

基于飛行器再入的標稱軌跡得到控制能力和阻力加速度的線性關系即其應滿足的約束，約束條件為

Ax=b，xmini≤xi≤xmaxi

(16)

其中矩陣A為(N+5)×2N矩陣，x為2N維列向量，b為N+5維列向量，2N>N+5。優化指標為

J=max(min(λi·dis(xi,Ui))i=1,2，…，2N)

(17)

λi為設置的權系數，其中

Ui=[xmini,xmaxi]

(18)

(19)

該優化指標的意義是讓x的每個分量都盡可能的遠離其取值范圍的邊界。

記滿足約束條件中等式約束的x的解集為Ω1，記滿足第二個約束的x的集合為Ω2，易知Ω2為凸集。該指標的意義是在集合Ω2中尋找一點x，使該點到Ω2邊界的距離最遠，若結合等式約束，則完整意義是在集合Ω1和集合Ω2的交集中尋找一點x，使該點到Ω2邊界的距離最遠。原指標為無法求導的非線性函數，因此考慮將原指標改寫為可導函數。在二維空間中考慮原指標，如圖2所示，Ω2的邊界即原指標意義下的等值線，對等值線進行縮放，當與Ω1只有一個交點時，該點即原指標下的最優解。

圖2 原指標二維空間示意圖

將原指標函數改寫為可導函數，理論上改寫后指標等值線越接近原指標等值線則優化結果越接近原指標，因此考慮將指標寫為如下形式

(20)

其中a為自然數，a越大則其等值線越接近原指標的等值線，但函數形式也越復雜，求解時運算量也更大，因此采用a=1時的指標，如下

(21)

圖3 二維空間示意圖

圖4 新指標與原指標等值線

如前所述，該優化問題的最優解即Ω1與等值線相切的點，Ω1為一仿射集，在Ω1與等值線相切的點處，等值線梯度方向垂直于Ω1，即垂直于Ω1對應子空間的基向量，由矩陣論相關知識可知，矩陣(I-A+A)中線性無關的列向量即Ω1對應子空間的基向量。最優指標等值線梯度為

(22)

其中

(23)

由上式可知梯度是關于x的一次函數，因此滿足下列方程的解即為最優解

(24)

整理可得

(25)

其中

整理得

(26)

該方法可直接求最優解。解出新軌跡后可計算優化后的軌跡航程R′，如果R′與原軌跡航程差距較大，則更新b2，令b2=b2+R′-R，再計算上式解出x，重復1-3次即可。

本文中制導控制方法采用PID控制，方法如下

(27)

其中ur為參考軌跡的標稱控制量，解算出σ即可，Dr是參考阻力加速度，δD=D-Dr。

5 仿真

為了驗證本文方法的有效性，根據某參考對象數據進行仿真驗證，速度-攻角曲線如圖5。

圖5 再入攻角曲線

5.1 再入傾角范圍

設再入高度為60km，阻力加速度為4m/s2，再入速度為5380m/s，則不同的再入傾角對應的初始V-D曲線的斜率如下圖

圖6 再入初始段軌跡傾角與V-D曲線斜率

由圖中可知，再入初始段V-D曲線斜率隨再入傾角變化比較劇烈，當再入傾角比較大時，V-D曲線的斜率為負，隨著再入傾角的減小，V-D曲線斜率也不斷減小，但是當再入傾角進一步減小到-20°以下時，在某個再入傾角下V-D曲線斜率突變為正，這是由于高空再入傾角進一步減小的時候，由于阻力加速度比較小，沿速度方向重力加速度的分量大于阻力加速度，這會導致初始再入段飛行器的速度增加、阻力加速度增加，反映在V-D曲線斜率上即斜率為正。

設再入段結束時刻阻力加速度為50m/s2，速度為1800m/s，則不同的軌跡傾角對應的末端V-D曲線的斜率如圖7。

圖7 再入末段航跡傾角與V-D曲線斜率

結束時刻阻力加速度一般比較大，重力加速度分量不會改變速度變化方向，由圖中可見此時的V-D曲線斜率隨軌跡傾角的變化在一個小范圍內線性變化。根據以上仿真結果，再入時刻的軌跡傾角取0°～-10°，再入結束時刻的軌跡傾角范圍可以大一些。

5.2 優化算法的有效性驗證

初始段軌跡規劃時Vd-Vc∶Vc-Vb∶Vb-Va=1∶2∶1，再入高度為60km，初始速度為4500m/s，再入初始段阻力加速度2.9m/s2，再入終點阻力加速度為50m/s2，終點速度為1800m/s，航程350km，再入傾角-1°，末端傾角0°，最大可用控制能力為控制能力極限的0.6倍。采樣間隔為10m/s，航程迭代3次，運行環境為windos7，i5-6500，Matlab2017a，優化前后的曲線對比如下

圖8 優化前后阻力加速度曲線

圖9 優化前后需用控制力

圖10 優化前后航跡傾角

優化前曲線的航程為350km，優化后為351.3km，可見航程誤差小，線性化假設合理。通過上圖可得知該方法保證了再入初始段和再入末段的航跡傾角滿足約束，并改善了控制能力，使所需控制能力在飛行器控制能力極限之內。下面通過仿真跟蹤優化前后的曲線驗證優化算法的有效性，采用PD控制，其中控制參數的阻尼、頻率分別選為1、0.8，可用的最大、最小傾側角為175°、5°。優化前后的曲線跟蹤仿真結果如下。

圖11 阻力加速度曲線跟蹤

圖12 阻力加速度曲線跟蹤誤差

圖13 航跡傾角

圖14 傾側角變化

由以上仿真可知，優化后的軌跡航程接近原軌跡，且與飛行器的控制能力更加匹配，優化后軌跡的跟蹤精度好于優化前，再入段末端可以更好地滿足傾角約束。進行了4組仿真，驗證該方法的有效性，表1中列出了4組仿真案例的條件，分別是再入速度、再入段初始阻力加速度、再入段末端阻力加速度、初始航跡傾角、末端航跡傾角、航程。

表1 仿真案例初始條件

表2 中列出了4組仿真案例下軌跡優化前和優化后的跟蹤效果，優化時間均在0.3s內，有3個指標，分別是最大阻力加速度偏差、末端阻力加速度偏差、末端航跡傾角偏差。

表2 優化前后仿真案例對比

根據上述結果可知，該優化方法可以在有限時間內將原軌跡修正為控制能力符合飛行器真實控制能力的軌跡，滿足再入走廊和航程、航跡傾角等約束。

6 結 論

針對飛行器再入段軌跡規劃問題，本文提出了一種基于橢球近似解析法的再入軌跡優化算法。經仿真，該方法可以較好地滿足再入傾角約束和控制能力約束等軌跡約束，運算速率快，有一定應用潛力。