等再入航程返回軌道模糊混合優化設計

陳偉躍,王國軍,王陽,陳蒙,張治國

中國空間技術研究院,北京 100094

1 引言

長期運行的航天器的軌道高度是變化的,如國際空間站的運行軌道高度在340～483km的范圍內。對于長期運行的航天器變高度返回的工程問題,通過返回軌道設計使得再入段過載、駐點熱流密度等彈道特性參數類似,對返回再入的工程實施具有重要意義。等航程再入彈道具有再入角一致,再入段航程近似不變、再入段過載和駐點熱流密度等特性類似的優點,有利于返回再入控制,因此可以利用等航程再入彈道來確保制動高度變化時返回確定著陸點的再入彈道特性一致。

在再入控制過程中不論是標稱軌道制導法還是預測校正制導算法,都需要用到標稱返回軌道的先驗信息[1-2],標稱返回軌道的優化對于返回再入飛行任務設計和制導及控制系統系統設計具有重要意義[3-5]。目前標稱返回軌道的設計方法多采用工程試算調試法[2]、有界試位迭代法[6]、微分修正法[7],此類方法依賴于設計經驗確定初始值,需要進行反復調試才能確定初值。Gauss偽譜法[8]、局部配點法[9]、凸優化[10]、序列凸規劃[11]等離散節點類優化算法則是標稱返回軌道優化使用較多的方法。采用Gauss偽譜法、局部配點法進行標稱返回軌道優化,為了達到較高的計算精度,需要增加離散點的數量,計算量較大。凸優化、序列凸規劃對求解的問題有凸性的要求,需要將彈道優化問題中原始非凸問題轉化為凸問題,轉化過程較為復雜。上述離散節點類方法節點數量多,處理過程較復雜,計算量大。

對于標稱返回軌道優化問題,序列二次規劃算法(SQP)[12-13]對初值敏感且容易陷入局部最優解,遺傳算法[14]不依賴于設計經驗確定初始值,但為了獲得收斂精度較高的結果,種群數量大,迭代次數多,計算量大。為了克服梯度類算法對初值的敏感性,并提高遺傳算法的計算精度,遺傳與序列二次規劃混合算法(GA-SQP)被應用于彈道優化設計[15-16]。文獻[15]在遺傳算法初始種群生成時采用混沌搜索,在執行選擇操作前選取符合單目標的優秀個體作為初值。文獻[16]在遺傳算法適應度計算時采用模擬退火罰函數處理約束。現有的遺傳與序列二次規劃混合算法(GA-SQP)[15-16]以遺傳算法給出的準全局最優解作為SQP的輸入,對遺傳算法的求解精度要求較高,種群大小較大,迭代次數較多。

如何減小種群大小,減少迭代次數,提高搜索速度,降低對遺傳算法搜索全局最優解的精度要求,同時保證遺傳算法得出的結果對于SQP而言是一個“較好”的初值,是改進遺傳與序列二次規劃混合算法并將其應用于工程實踐的關鍵問題。針對此問題,本文基于模糊評價[17]提出了改進的穩態遺傳算法[18]與SQP混合的求解方法,并在等航程返回軌道設計中進行了應用。此方法不需要提供初始值,有效地提高了計算效率,適用于近地及深空返回再入軌跡優化,具備良好的工程應用前景。

2 等再入航程返回軌道

從不同高度返回時,若再入段航程相等,則稱之為等再入航程。當制動點的初始條件變化時,可通過恰當的選擇制動點和制動速度增量,使得再入段的航程偏差控制在一定范圍內,如±2km,在工程實際中可以認為再入段航程相等。即使制動點的初始條件產生較大變化,再入段航程仍相等,再入環境類似,有利于再入控制。

令航程為S,速度為v,當地速度傾斜角為θ,地球平均半徑為Rp,再入過程中高度為h,開傘點的高度為hf,航天器質量為m,升力為L。采用一階近似解進行分析,可以忽略重力和離心力對再入彈道的影響[19],則有:

(1)

對于近地航天器再入彈道,h?Rp,則有

(2)

式中:CD為阻力系數;s為參考面積;D為阻力;ρ為大氣密度,ρ=ρ0e-βh,β=1.3889×10-5,ρ0為高度為0處的大氣密度。

忽略重力和離心力,根據二階近似解有[19]

(3)

式中:θe為再入角;ρe為再入點的大氣密度。將式(3)代入式(2)可得

(4)

其中

對于近地軌道的半彈道式和彈道式返回,|θe|<5°,因此有cosθe≈1。再入點的大氣稀薄,可以認為ρe≈0,因此b≈-1,則積分(4)式可得

S=(θf-θe)/β+[cot(θf/2)-cot(θe/2)]/β

(5)

由(3)式可知,當開傘點的高度hf、升阻比、阻力系數、參考面積、質量和大氣模型確定的情況下,開傘點的速度傾角θf由再入角θe確定。再結合(5)式可知,再入段航程S由再入角θe確定。等再入航程返回軌道的基本特征是再入角一致,再入航程一致,再入段的彈道特性保持近似不變。若設計的返回軌道對某項特性具有最優性,在初始條件變化時,采用等再入航程返回軌道該項特性仍可保持最優性,有利于再入控制。

一般出現最大加速度時,Ma較大,可認為升力系數CL和阻力系數CD為常數。忽略重力和離心力對再入彈道的影響[19],最大過載Nmax可表述為

(6)

采用經驗公式[20]計算駐點熱流密度:

(7)

式中:Rn為駐點曲率半徑;vc為航天器近地軌道環繞速度;c為經驗系數。

在升阻比、阻力系數、參考面積、質量、大氣模型確定的情況下,再入角決定最大熱流密度出現時的速度傾角、高度以及速度與再入速度之比。等航程再入彈道最大熱流密度出現的高度一致。

3 等再入航程返回軌道設計方法

3.1 等再入航程返回軌道設計問題的數學描述

將航天器星下點軌跡經過指定著陸區域的圈次稱為返回圈。航天器返回圈升交點的軌道六根數采用σ表示,航天器的軌道運動表示如下:

(8)

航天器的軌道運動方程可根據需要考慮非球形引力攝動、大氣阻力攝動等軌道攝動項[21-22]。

以返回圈升交點時刻為時間起點,設制動開機時刻相對于返回圈升交點時刻的相對時間為Tb。Tb>0表示在返回圈升交點之后制動,Tb<0表示在返回圈升交點之前制動。

在制動開機點建立返回坐標系o0x0y0z0。原點是航天器質心在制動時刻的地心矢量與標準地球橢球體表面的交點o0,o0y0軸沿o0點與制動時刻航天器質心連線的方向,o0x0在返回制動時刻航天器運行軌道平面內且與o0y0相互垂直,并指向航天器的運動方向,o0z0與o0x0、o0y0軸組成右手直角坐標系;在再入點(高度為100km)建立再入坐標系oRxRyRzR,原點在再入時刻的地心矢量與標準地球橢球體表面的交點oR,oRyR軸沿oR點與再入時刻航天器質心連線的方向,oRxR在返回制動時刻航天器運行軌道平面內且與oRyR相互垂直,并指向航天器的再入運動方向,oRzR與oRxR、oRyR軸組成右手直角坐標系。再入運動方向由再入點地心矢量與返回瞄準點地心矢量構成的大圓在再入點的地心方位角描述;半速度坐標系為o1xhyhzh,該坐標系的原點o1為航天器質心,o1xh軸沿速度方向,與速度坐標系o1xv方向重合,o1yh在返回坐標系x0o0y0平面內垂直于o1xh軸,o1xhyhzh構成右手直角坐標系。體坐標系為o1xByBzB,該坐標系的原點o1為航天器質心,o1xB為航天器的縱軸,指向頭部,o1yB在航天器主對稱面內,該平面在返回制動時刻與返回坐標系x0o0y0平面重合,o1yB垂直于o1xB,o1zB垂直于主對稱面,沿運動方向看o1zB指向右方。坐標系關系如圖1所示。

圖1 坐標系關系

圖1中右邊子圖的角度定義為:θ=arctg(vy/vx),σ=arcsin(-vz/v)。vx,vy,vz為返回坐標系下航天器相對于地球的飛行速度分量,v為航天器相對于地球的飛行速度大小,θ為速度傾角,σ為航跡偏航角。

航天器返回再入質心動力學方程描述如下[1]:

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

等再入航程返回軌道設計的數學描述為,求解制動開機時刻Tb、制動開機時長Tp、再入段傾側角序列γi(i=1,2,3,4),優化目標函數為開傘高度處的航程偏差[1]ΔRL與橫程偏差[1]ΔRZ的平方和最小:

(19)

并使得再入角θe、再入過載N、再入段航程S和過載超過過載界限的累計時間tg滿足如下約束[23]:

Smin≤S≤Smax,θe=θc,N≤Nmax,tg≤tmax

(20)

式中:Smin為允許的再入段航程最小值;Smax為允許的再入段航程最大值;θc為滿足再入航天器氣動熱設計要求的目標再入角;tmax為過載超過過載界限G0的累計時間的最大允許值。

3.2 等再入航程返回軌道設計問題的求解方法

本文給出一種基于模糊評價的改進穩態遺傳算法與序列二次規劃算法的混合求解方法。首先采用穩態遺傳算法搜索等再入航程問題的粗略解,采用模糊評價方法確定搜索結果的好壞。將改進的穩態遺傳算法搜索的粗略解作為序列二次規劃算法的初始值,采用序列二次規劃算法獲得最終解。此方法不需要設計者給出初值,可自行搜索到較好的粗略解,同時避免陷入局部最優解,并可以獲得較高的求解精度。

3.3 等再入航程返回軌道設計問題的求解步驟

(1)采用改進的穩態遺傳算法求解初值

采用穩態(Steady State)遺傳算法[18]求解等再入航程返回軌道設計問題時,根據優化變量Tb、制動開機時長Tp、再入段傾側角序列γi(i=1,2,3,4)的數值范圍采用二進制編碼隨機產生初始種群。此后克隆上一代的原始種群中一定比例或一定數量的個體以產生臨時種群,臨時種群按照一定的概率交叉、突變后與原始種群合并,合并后按照評價分數依次移除評分最低的個體,種群中個體數量恢復至原始種群的數量時,完成新一代種群的產生。以每一代及之前qc代的最好分數與上一代及其之前qc代的最好分數的比值作為收斂度,當收斂度達到指定的值時,遺傳算法搜索停止。

使用穩態(steady state)遺傳算法求解等再入航程返回軌道設計問題的粗略解時有3個關鍵的步驟:①隨機產生初始種群,②根據最差評分移除劣質個體,③根據評分計算收斂度,依據收斂度停止迭代。由于初始種群是隨機產生的,在第一代種群中產生的劣質基因需要在步驟②中根據評分移除,穩態遺傳算法新一代種群與上一代種群存在一定比例的重疊,完全剔除劣質基因的影響需要一定的遺傳代數,第一代種群的劣質基因會減緩收斂速度。步驟②和③依賴于遺傳算法適應度函數的形式,對于多目標優化問題,鮮見采用模糊思想建立評分函數。以下針對等再入航程返回軌道設計問題,給出初始種群的基因檢測方法、基于模糊隸屬度的評分函數建立方法和搜索過程中劣質基因的評分方法。

為適應較大的運行軌道高度變化范圍,等再入航程返回軌道設計的優化變量數值存在較大的變化范圍。隨機生成的初值中,制動開機時刻初值取值不合適可能導致再入點在瞄準點之后;制動開機時長初值取值不合適可能導致制動完成后滑行軌道的近地點高度較高,無法形成再入彈道。初始基因檢測方法用于對隨機產生的初始種群的健康狀態進行檢查,存在上述兩種情況中的任何一種均標記為不健康個體,對于不健康的個體需要重新隨機生成。基因檢查方法如下:對于隨機產生的初始種群個體X0l(l=1,2,…,l,l為種群大小),計算式(8)～(18),檢查式(21)(22)。其中Hp為制動后滑行軌道的近地點高度,H0為制動后滑行軌道的近地點高度的允許上限,Ae為再入點地心方位角。若如下2式之一成立

Hp≥H0

(21)

Ae<0

(22)

則標記為不健康個體,剔除不健康個體并重新隨機生成。

對于式(19)～(20)描述的優化問題,基于模糊隸屬度[17]建立遺傳算法的評分函數。

對于式(19)采用正態分布描述開傘點偏差的好壞程度:

μ(d)=e-(d/a)2

(23)

對于再入段航程約束采用正態分布描述約束滿足的好壞程度:

(24)

對于再入角等式約束采用“Λ”型模糊隸屬度函數[17]描述再入角約束滿足的好壞程度:

(25)

式中:θ0根據需要取值,如θ0=0.1°;θc為再入角目標值,認為再入角偏離目標值±0.1°,再入角的優化結果不滿足,隸屬度為0。

對于最大過載約束N≤Nmax和過載超限時長約束tg≤tmax,采用戒上型隸屬度函數[17]描述。對于最大過載約束滿足程度的好壞描述為:

μ(N)=

(26)

式中:N0根據需要取值,如N0=0.1g,即認為最大過載超過允許值的數值達到0.1g時,最大過載約束不滿足,隸屬度為0。

對于過載超限時間約束滿足的好壞程度如下:

(27)

式中:t0根據需要取值,如t0=10s,即認為最大過載超過過載界限G0累計時間達到10s,過載超限時長約束不滿足,隸屬度為0。

遺傳算法的評分函數建立如下:

V=100[λ1μ(d)+λ2μ(S)+λ3μ(θe)+

λ4μ(N)+λ5μ(tg)]

(28)

式中:λj(j=1,2,3,4,5)為模糊綜合評分的權重[17],須滿足如下約束:

(29)

式(28)將遺傳算法解的隸屬度轉化成了百分制評分函數,利用式(28)采用穩態遺傳算法求解等再入航程返回軌道設計問題的粗略解。

在初始種群生成過程中進行了基因檢測,對于不健康個體重新隨機生成。由于在搜索過程中每一代種群個體的基因存在變異,仍然可能產生劣質基因。以下給出搜索過程中劣質基因的評分方法。

制動開機點基因變異可能導致再入點在瞄準點之后,令d、S、θe、N、tg的取值分別為:

d=1000km,S=-Smax,θe=2θc,

N=2Nmax,tg=2tmax

(30)

將式(30)代入式(28),帶有劣質基因的個體評分結果為0,將在搜索過程中作為最差的個體被移除。

制動開機時長基因變異可能導致制動完成后滑行軌道的近地點高度較高無法形成再入彈道,令d、S、θe、N、tg的取值分別為:

d=1000km,S=0,θe=0,

N=2Nmax,tg=2tmax

(31)

將式(31)代入式(28),帶有劣質基因的個體評分結果為0,將在搜索過程中作為最差的個體被移除。

(2)采用序列二次規劃算法求解

將遺傳算法計算得到的粗略解作為SQP的初始值,采用文獻[24]給出的序列二次規劃算法計算等再入航程返回軌道設計問題的最終解。對于(20)式中再入航程約束、最大過載約束和過載超限時間約束,需要轉化為g(x)≥0的形式,g1(S)與g2(S)為航程約束,g3(θe)為再入角約束,g4(N)為最大過載約束,g5(tg)為過載超限時間約束,約束函數的形式如下:

g1(S)=S-Smin,g2(S)=-(S-Smax),

g3(θe)=θe-θc,g4(N)=-(N-Nmax),

g5(tg)=-(tg-tmax)

(32)

對于制動點不合適導致再入點在瞄準點之后的特殊情況,優化目標函數(19)和約束函數(32)中的S、θe、N、tg按照式(30)取值,d=2000km。

對于制動時長過短導致無法形成再入彈道的特殊情況,按照約束函數連續可微分的原則,d、S、θe、N、tg取值如下:

d=1000km,S=0,θe=0,N=0,tg=0

(3)等再入航程返回軌道設計問題的求解流程

等航程再入彈道設計問題的求解流程包括兩個主要步驟,首先采用改進的穩態遺傳算法獲得粗略解,然后將其作為序列二次規劃求解的初始值,采用序列二次規劃獲得等再入航程彈道設計問題的精確解。改進的穩態GA-SQP的具體計算步驟如下:

1)穩態遺傳算法參數設置。設置種群大小、突變概率、交叉概率、收斂度界限、計算收斂度使用的代數。

2)隨機生成初始種群。

3)初始種群個體的基因檢查,若個體基因不健康則重新隨機生成個體。

4)采用輪盤賭選擇算法按照一定比例或指定個體數量從上一代種群中選擇臨時種群的父母基因。

5)交叉、變異。按照概率交叉、變異,生成臨時種群。

6)將臨時種群與上一代種群合并。

7)依次剔除評分最差的個體至種群數量恢復至原始值。

8)檢查收斂度是否滿足,若不滿足繼續生成新一代種群,否則輸出等再入航程彈道設計問題的初略解。

9)序列二次規劃算法參數設置。以穩態遺傳算法得出的粗略解作為優化變量的初始值,定義優化目標函數、約束函數。

10)將帶約束的優化問題轉化為拉格朗日擴展的優化問題。

11)在當前解Xk處將約束線性化獲得拉格朗日擴展優化問題的二次逼近。

12)求解二次逼近子問題。

13)以擴展的拉格朗日函數作為優化目標函數線性搜索下一步解Xk+1。

14)判斷是否滿足收斂精度要求,如果滿足精度要求則輸出等再入航程彈道設計問題的精確解,否則采用BGFS方法更新二次逼近子問題的二階偏導數矩陣[24],轉至步驟11)。

GA-SQP算法計算流程如圖2所示。

圖2 GA-SQP算法計算流程

4 等再入航程返回軌道設計實例

以CEV從國際空間站ISS軌道以半彈道式返回至愛德華空軍基地(W 117°50′06.77",N 34°55′53.02")為例,說明等再入航程返回軌道的設計。國際空間站ISS的初始軌道根數如表1所示,CEV與返回相關的參數如表2所示[25]。

表1 ISS返回圈升交點J2000坐標系下的初始軌道根數

表2 CEV與返回相關的參數

為對本文提出的改進的穩態GA-SQP與未改進的穩態遺傳-序列二次規劃算法(UGA-SQP)進行比較,選取UGA-SQP的適應度函數為[26]:

采用UGA-SQP優化時,根據所得的優化解采用式(28)得出評分值,在改進的穩態GA-SQP與UGA-SQP評分值偏差小于0.5分的情況下進行對比。為使得UGA-SQP給出的解的評分值與改進的穩態GA-SQP給出的解的評分值偏差小于0.5分,需要調整UGA-SQP的種群大小和收斂代數,其他參數兩種方法取值一致。兩種方法的SQP部分均考慮對再入點在瞄準點之后和制動時長過短兩種特殊情況的處理。

采用VC++編程語言實現本文提出的改進的穩態GA-SQP和UGA-SQP,運行環境為Win7系統,8核CPU,主頻3.40GHz。半長軸為6715.8km時,采用改進的穩態GA-SQP和UGA-SQP分別計算5次的結果如表3和表4所示,表3和表4中p為種群大小,q為遺傳算法的遺傳代數、SQP算法的迭代次數。采用改進的穩態GA-SQP優化時,種群大小為20,收斂代數為40。在改進的穩態GA-SQP與UGA-SQP評分值偏差小于0.5分的前提下,采用UGA-SQP優化時,種群大小為80,收斂代數為80,采用UGA-SQP優化時的種群大下與文獻[15]一致。

表3 改進的穩態GA-SQP等航程再入彈道的設計結果

表4 未改進的GA-SQP等航程再入彈道的設計結果

采用改進的穩態GA-SQP算法計算半長軸變化時的結果如表5所示,圖3與圖4為表5中編號為1的工況對應的再入彈道設計結果,圖4中的經度范圍為0°～360°。

圖3 橫向偏差及傾側角曲線

圖4 CEV返回過程的星下點軌跡

從表3可以看出,在參數設置不變的情況,由于遺傳算法隨機搜索的特性,每次計算遺傳算法(GA)得出的收斂程度是不一致的,模糊綜合評分結果V的百分數值可以表示遺傳算法搜索結果的好壞程度。對于改進的穩態遺傳算法的搜索結果,SQP均可找到滿足各項約束的設計結果。表3中5次計算的制動開機時刻與平均值的最大偏差為146ms,制動開機時長與平均值的最大偏差為5ms,開傘點最大偏差1.141km,再入角偏差小于0.001°,再入航程約束的上限偏差為1.234km,滿足最大過載約束,超過3.0g的時間最大為9.2s。采用改進的穩態GA-SQP計算最小耗時為169s,最大耗時為261s。從表4可以看出,表4中5次計算的制動開機時刻與平均值的最大偏差為720ms,制動開機時長與平均值的最大偏差為3ms,除第二次計算開傘點偏差略大外其他幾次計算開傘點偏差均滿足要求,再入角偏差小于0.001°,再入航程約束的下限偏差為-5.652km,過載約束存在未滿足的情況,最大過載偏差為0.017g,超過3.0g的時間最大為15.4s。采用UGA-SQP計算最小耗時1161s,最大耗時為2027s,采用UGA-SQP計算遺傳代數會明顯增加,導致計算耗時顯著增長。上述對比仿真表明,在不降低計算精度的前提下改進的穩態GA-SQP種群數量減小,迭代次數減少34.7%,計算耗時明顯縮短。改進的穩態GA-SQP與工程試算調試法[2]、有界試位迭代法[6]相比,不需要設計者給出初值,可自行搜索到較好的解,計算精度和速度均達到工程適用程度。

從表3可以看出傾側角γ1～γ4的大小整體呈現依次減小的規律,不同次數傾側角的計算結果存在一定程度的變化,為再入控制系統的設計留有一定的自由度。從表5可以看出,在固定再入角的條件下,當ISS軌道高度增大時,制動點的位置沿著返回圈的星下點軌跡后退,制動點高度增加,再入彈道的再入點位置不變,制動時長和速度增量隨軌道高度的增高而增大。再入段航程與平均值的最大偏差為2km,實現了制動高度變化時再入段航程一致的設計目標,且再入段最大過載基本一致。

在ISS軌道高度發生變化的情況下,再入彈道的特性保持不變,對再入控制與防熱是十分有利的。

5 結論

本文推導了再入航程與再入角的關系,分析了等航程再入彈道的優點,提出了基于模糊綜合評分的穩態遺傳與序列二次規劃混合求解算法。針對多目標等航程再入彈道設計問題,設計了穩態遺傳算法初始種群的基因檢測方法,給出了基于模糊隸屬度的評分函數建立方法、搜索過程中的劣質基因評分方法,針對序列二次規劃求解過程中再入點在瞄準點之后和制動時長過短的特殊情況,給出了約束的處理方法。考慮了再入過載超限時間約束、傾側角翻轉的最大角速度限制等實際工程約束。

以CEV從ISS軌道返回為例,說明了等航程再入彈道設計方法和設計結果的特性。采用改進的穩態遺傳算法的初始搜索結果,SQP算法均可找到滿足各項約束的設計結果,計算精度達到工程適用程度,為控制系統的設計留有一定的自由度,實現了制動高度變化時再入彈道特性一致的設計目標,計算速度達到工程適用程度。

針對軌道高度發生變化的情況,等航程再入彈道能保證再入彈道的基本特性不變,為返回再入控制、防熱創造良好的條件。本文的結果可為中國空間站階段的返回軌道及返回方案設計提供參考。