張 凱 楊小龍 鐘 震

北京宇航系統(tǒng)工程研究所，北京，100076

0 引言

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

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

1 飛行器再入模型

設(shè)飛行器再入過(guò)程中地球?yàn)閳A球，飛行器為剛體，主要受氣動(dòng)升力、氣動(dòng)阻力和引力作用，則飛行器三維的運(yùn)動(dòng)學(xué)、動(dòng)力學(xué)方程描述如下[12]

(1)

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

升力和阻力加速度的表達(dá)式為

(2)

整理可得

(4)

其中

(5)

設(shè)R為再入段航程，則有

(6)

根據(jù)上式可以通過(guò)V-D曲線、速度-航跡傾角曲線來(lái)預(yù)測(cè)航程。

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

(7)

其中

(8)

飛行器不同的再入傾角所引起的阻力加速度變化不同，同時(shí)，有些飛行器再入段結(jié)束時(shí)對(duì)航跡傾角有要求，如果采用單獨(dú)的再入段航跡傾角調(diào)節(jié)控制，會(huì)增加軌跡設(shè)計(jì)的復(fù)雜性，因此需要找出航跡傾角與規(guī)劃的V-D曲線之間的關(guān)系，保證所規(guī)劃的軌跡滿足航跡傾角約束。根據(jù)前文推導(dǎo)，將軌跡傾角表達(dá)式移項(xiàng)整理可進(jìn)一步得到

(9)

其中等式右端所有的量都可知，因此當(dāng)?shù)玫匠跏蓟蚪K端航跡傾角時(shí)，就可以計(jì)算出初始、終端的V-D曲線的斜率。記

(10)

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

3 基于標(biāo)稱軌跡的線性化模型

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

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

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

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

(11)

將V-D曲線分成等距的N個(gè)點(diǎn)，將全N個(gè)點(diǎn)都考慮進(jìn)去，可得

(12)

將上式簡(jiǎn)記為A1x1=x2。其中

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

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

(13)

分別簡(jiǎn)記為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為原規(guī)劃曲線的阻力加速度、控制量，Dmax、Dmin為再入走廊的阻力加速上下限，umax、umin為控制能力的上下限，其含義為修正軌跡后保證新軌跡不超過(guò)再入走廊的邊界，同時(shí)所需控制力不超過(guò)飛行器的控制能力極限。其中Dmax、Dmin、umax和umin也可根據(jù)實(shí)際情況進(jìn)行調(diào)整。

4 基于橢球近似的解析優(yōu)化算法

基于飛行器再入的標(biāo)稱軌跡得到控制能力和阻力加速度的線性關(guān)系即其應(yīng)滿足的約束，約束條件為

Ax=b，xmini≤xi≤xmaxi

(16)

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

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

(17)

λi為設(shè)置的權(quán)系數(shù)，其中

Ui=[xmini,xmaxi]

(18)

(19)

該優(yōu)化指標(biāo)的意義是讓x的每個(gè)分量都盡可能的遠(yuǎn)離其取值范圍的邊界。

記滿足約束條件中等式約束的x的解集為Ω1，記滿足第二個(gè)約束的x的集合為Ω2，易知Ω2為凸集。該指標(biāo)的意義是在集合Ω2中尋找一點(diǎn)x，使該點(diǎn)到Ω2邊界的距離最遠(yuǎn)，若結(jié)合等式約束，則完整意義是在集合Ω1和集合Ω2的交集中尋找一點(diǎn)x，使該點(diǎn)到Ω2邊界的距離最遠(yuǎn)。原指標(biāo)為無(wú)法求導(dǎo)的非線性函數(shù)，因此考慮將原指標(biāo)改寫為可導(dǎo)函數(shù)。在二維空間中考慮原指標(biāo)，如圖2所示，Ω2的邊界即原指標(biāo)意義下的等值線，對(duì)等值線進(jìn)行縮放，當(dāng)與Ω1只有一個(gè)交點(diǎn)時(shí)，該點(diǎn)即原指標(biāo)下的最優(yōu)解。

圖2 原指標(biāo)二維空間示意圖

將原指標(biāo)函數(shù)改寫為可導(dǎo)函數(shù)，理論上改寫后指標(biāo)等值線越接近原指標(biāo)等值線則優(yōu)化結(jié)果越接近原指標(biāo)，因此考慮將指標(biāo)寫為如下形式

(20)

其中a為自然數(shù)，a越大則其等值線越接近原指標(biāo)的等值線，但函數(shù)形式也越復(fù)雜，求解時(shí)運(yùn)算量也更大，因此采用a=1時(shí)的指標(biāo)，如下

(21)

圖3 二維空間示意圖

圖4 新指標(biāo)與原指標(biāo)等值線

如前所述，該優(yōu)化問(wèn)題的最優(yōu)解即Ω1與等值線相切的點(diǎn)，Ω1為一仿射集，在Ω1與等值線相切的點(diǎn)處，等值線梯度方向垂直于Ω1，即垂直于Ω1對(duì)應(yīng)子空間的基向量，由矩陣論相關(guān)知識(shí)可知，矩陣(I-A+A)中線性無(wú)關(guān)的列向量即Ω1對(duì)應(yīng)子空間的基向量。最優(yōu)指標(biāo)等值線梯度為

(22)

其中

(23)

由上式可知梯度是關(guān)于x的一次函數(shù)，因此滿足下列方程的解即為最優(yōu)解

(24)

整理可得

(25)

其中

整理得

(26)

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

本文中制導(dǎo)控制方法采用PID控制，方法如下

(27)

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

5 仿真

為了驗(yàn)證本文方法的有效性，根據(jù)某參考對(duì)象數(shù)據(jù)進(jìn)行仿真驗(yàn)證，速度-攻角曲線如圖5。

圖5 再入攻角曲線

5.1 再入傾角范圍

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

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

由圖中可知，再入初始段V-D曲線斜率隨再入傾角變化比較劇烈，當(dāng)再入傾角比較大時(shí)，V-D曲線的斜率為負(fù)，隨著再入傾角的減小，V-D曲線斜率也不斷減小，但是當(dāng)再入傾角進(jìn)一步減小到-20°以下時(shí)，在某個(gè)再入傾角下V-D曲線斜率突變?yōu)檎@是由于高空再入傾角進(jìn)一步減小的時(shí)候，由于阻力加速度比較小，沿速度方向重力加速度的分量大于阻力加速度，這會(huì)導(dǎo)致初始再入段飛行器的速度增加、阻力加速度增加，反映在V-D曲線斜率上即斜率為正。

設(shè)再入段結(jié)束時(shí)刻阻力加速度為50m/s2，速度為1800m/s，則不同的軌跡傾角對(duì)應(yīng)的末端V-D曲線的斜率如圖7。

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

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

5.2 優(yōu)化算法的有效性驗(yàn)證

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

圖8 優(yōu)化前后阻力加速度曲線

圖9 優(yōu)化前后需用控制力

圖10 優(yōu)化前后航跡傾角

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

圖11 阻力加速度曲線跟蹤

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

圖13 航跡傾角

圖14 傾側(cè)角變化

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

表1 仿真案例初始條件

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

表2 優(yōu)化前后仿真案例對(duì)比

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

6 結(jié) 論

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