閆循良,湯一華,徐 敏,陳士櫓

(西北工業(yè)大學(xué)航天學(xué)院,西安 710072)

0 引言

空間交會(huì)是實(shí)現(xiàn)空間操作和應(yīng)用的關(guān)鍵。航天器往往通過(guò)一系列軌道機(jī)動(dòng)實(shí)現(xiàn)交會(huì),以完成對(duì)接、空間捕獲等任務(wù)。要使航天器在各種限制條件下完成不同的空間任務(wù),進(jìn)行軌道機(jī)動(dòng)任務(wù)規(guī)劃和軌道優(yōu)化研究是至關(guān)重要的。傳統(tǒng)的軌道機(jī)動(dòng)方式為沖量機(jī)動(dòng)[1-5],它是工程問(wèn)題的簡(jiǎn)單近似。沖量變軌對(duì)研究軌道機(jī)動(dòng)仍具有重要意義。在許多情況下,沖量軌道可作為問(wèn)題的初步模型,而且所得的數(shù)值結(jié)果可為問(wèn)題的解決提供重要的參考數(shù)據(jù)。軌道機(jī)動(dòng)最優(yōu)需綜合考慮燃料和機(jī)動(dòng)時(shí)間。因此,對(duì)于沖量變軌必須進(jìn)行任務(wù)規(guī)劃,以確定最佳變軌程序。要確定能使總特征速度極小的最優(yōu)速度變化程序,一般求取解析解是不可行的,必須用計(jì)算機(jī)來(lái)進(jìn)行分析計(jì)算,以確定某個(gè)特定飛行任務(wù)的適當(dāng)發(fā)射窗口。

實(shí)際變軌過(guò)程中,推力大小是有限的,推進(jìn)也不是瞬間的。因此,必須研究有限推力下的變軌問(wèn)題。目前,研究有限推力交會(huì)軌道優(yōu)化問(wèn)題的方法主要有間接法和直接法[6-11]。間接法主要應(yīng)用龐特里亞金極值原理對(duì)優(yōu)化數(shù)學(xué)模型進(jìn)行處理,并采用數(shù)值方法求解兩點(diǎn)邊值問(wèn)題,得到最優(yōu)控制量和相應(yīng)的軌道。此方法具有較好的求解精度,但需較為準(zhǔn)確的、對(duì)沒(méi)有物理意義的協(xié)狀態(tài)初值進(jìn)行猜測(cè)。直接法是將軌道優(yōu)化轉(zhuǎn)化為參數(shù)優(yōu)化問(wèn)題,并采用非線性規(guī)劃進(jìn)行數(shù)值求解。此方法不需求解兩點(diǎn)邊值問(wèn)題,收斂半徑大,但求解精度較低,計(jì)算量大,控制曲線不光滑。

文中首先以沖量機(jī)動(dòng)作為問(wèn)題的初步模型,采用Battin-Vaughan算法對(duì)雙沖量機(jī)動(dòng)初始位置和飛行時(shí)間的組合進(jìn)行遍歷計(jì)算,通過(guò)分析和解釋特征速度等值線圖,進(jìn)行空間交會(huì)的任務(wù)規(guī)劃,為有限推力燃料最優(yōu)交會(huì)提供重要的初值條件。基于任務(wù)規(guī)劃分析,建立了有限推力燃料最優(yōu)交會(huì)的最優(yōu)控制模型;采用龐特里亞金極值原理,將最優(yōu)控制問(wèn)題轉(zhuǎn)化為兩點(diǎn)邊值問(wèn)題,將對(duì)伴隨方程的初值猜測(cè)轉(zhuǎn)化為對(duì)物理意義明顯的控制量的猜測(cè);采用共軛梯度算法,對(duì)兩點(diǎn)邊值問(wèn)題進(jìn)行數(shù)值求解。在變軌時(shí)間固定、連續(xù)變推力的情況下,以總沖最小、滿足終端位置和速度約束為指標(biāo),對(duì)推力大小和方向進(jìn)行優(yōu)化。通過(guò)仿真對(duì)以上算法進(jìn)行了數(shù)值驗(yàn)證。

1 問(wèn)題描述

文中研究的問(wèn)題可描述為追蹤器和目標(biāo)(或虛擬目標(biāo)點(diǎn))分別在兩異面橢圓軌道上運(yùn)動(dòng),在多種限制條件下,追蹤器進(jìn)行有限推力變軌,并實(shí)現(xiàn)與目標(biāo)點(diǎn)的燃料最優(yōu)交會(huì),以完成對(duì)接或發(fā)射有效載荷等任務(wù)。

求解有限推力燃料最優(yōu)交會(huì)的關(guān)鍵,即尋求最佳飛行軌跡和推力大小、方向隨時(shí)間變化的規(guī)律,使之滿足交會(huì)的各種約束條件,并使燃耗達(dá)到最小。為了高效進(jìn)行有限推力軌道優(yōu)化,可將沖量機(jī)動(dòng)作為問(wèn)題的初步模型,通過(guò)任務(wù)規(guī)劃確定變軌初始位置r0、速度矢量v0,交會(huì)點(diǎn)位置rf、速度矢量vf,以及變軌時(shí)間Δt,進(jìn)而利用這些重要的初始數(shù)據(jù)完成有限推力軌道優(yōu)化。

2 空間交會(huì)快速任務(wù)規(guī)劃

2.1 雙沖量交會(huì)描述

假設(shè)目標(biāo)不做機(jī)動(dòng),若已知轉(zhuǎn)移時(shí)間Δt,根據(jù)飛行軌道可計(jì)算出目標(biāo)的末位置及速度矢量,即可得到交會(huì)目標(biāo)點(diǎn)的位置和速度信息。該問(wèn)題可描述為已知初始時(shí)刻追蹤器所在初始點(diǎn)P1的信息r1、v1,目標(biāo)點(diǎn)P2的信息轉(zhuǎn)移時(shí)間Δt,求解2次沖量Δ使得在Δt時(shí)間后,追蹤器與目標(biāo)點(diǎn)實(shí)現(xiàn)交會(huì)?？梢?jiàn),該問(wèn)題實(shí)質(zhì)是高斯問(wèn)題,可通過(guò)解高斯問(wèn)題得到施加第1次沖量后的速度、施加第2次沖量前的速度以及2次沖量Δ、Δ。

2.2 B-V算法解高斯問(wèn)題

目前,求解高斯問(wèn)題的多種方法已建立,除了文獻(xiàn)[3]中的解法以外,還有原始的高斯解法、普適變量法、p迭代法等。以上方法大多需對(duì)圓錐曲線類型進(jìn)行討論,且均無(wú)法避免轉(zhuǎn)移角為π時(shí)的奇異性;由Battin等人建立的B-V算法[4]對(duì)Gauss算法作出重大改進(jìn),計(jì)算簡(jiǎn)潔,移去了轉(zhuǎn)移角為的奇異性,并大幅度改善了算法對(duì)整個(gè)轉(zhuǎn)移角范圍0～2π的收斂性。

2.2.1 Battin-Vaughan算法

設(shè)地心為O,定義弦P1P2長(zhǎng)度c=|r2-r1|,三角形OP1P2周長(zhǎng)的一半s=(r1+r2+c)/2。經(jīng)過(guò)中值點(diǎn)變換[4]后,設(shè)P2點(diǎn)的真近點(diǎn)角為f,偏近點(diǎn)角為E,軌道偏心率為e0,則開(kāi)普勒方程為

參考高斯代數(shù)方程,定義變量:

其中,D為連接P1、P2的共軛拋物軌道的中值點(diǎn)半徑:

定義y,滿足

將以上變量帶入式(1),得

式(3)、式(4)與高斯第一、第二代數(shù)方程相似。

引入無(wú)量綱參數(shù)λ,滿足

得

綜上所述,該方法消除了Δf=π處的奇異性。

為了改進(jìn)求解非線性方程(3)、(4)的逐次代入算法的收斂性,引入h1、h2:

其中

以上方程均適用于各種圓錐曲線。

最終可導(dǎo)出Battin第一和第二代數(shù)方程:

利用B-V算法求解高斯問(wèn)題的具體步驟參見(jiàn)文獻(xiàn)[4]。

2.3 軌道機(jī)動(dòng)任務(wù)規(guī)劃

2.3.1 任務(wù)規(guī)劃分析

追蹤器和目標(biāo)分別在兩異面橢圓軌道上運(yùn)動(dòng),對(duì)于某一雙沖量空間交會(huì),需進(jìn)行軌道機(jī)動(dòng)任務(wù)規(guī)劃,以確定追蹤器的最佳初始位置、機(jī)動(dòng)時(shí)機(jī)、轉(zhuǎn)移時(shí)間及為完成任務(wù)所需的特征速度Δv。其中,追蹤器轉(zhuǎn)移時(shí)間和燃料受限。追蹤器在原飛行軌道的初始位置可用當(dāng)前時(shí)刻的真近點(diǎn)角f表示。由于追蹤器初始位置和轉(zhuǎn)移時(shí)間不定,需對(duì)初始位置對(duì)應(yīng)的真近點(diǎn)角f和轉(zhuǎn)移時(shí)間Δt的不同組合遍歷求解,進(jìn)而對(duì)結(jié)果進(jìn)行分析。

仿真中,轉(zhuǎn)移時(shí)間限制Δt∈[tmin,tmax],考慮燃料限制,單次沖量大小滿足Δv≤Δvmax,計(jì)算步長(zhǎng)為t-st;追蹤器初始真近點(diǎn)角f∈[0°,360°],計(jì)算步長(zhǎng)為f-st。因此,2次沖量大小均可表示為Δ=(f,Δt),Δ=F2(f,Δt),總沖量為

對(duì)不同的數(shù)據(jù)組合(Δt,f),按照2.2節(jié)所述方法進(jìn)行數(shù)值求解,得到Δ、ΔΔv,并可作出Δv等值曲線圖,將繪制和分析Δv等值曲線圖作為任務(wù)規(guī)劃的輔助工具[5]。根據(jù)Δv等值線圖,可得到滿足要求的追蹤器初始位置點(diǎn)集(即追蹤區(qū))以及可實(shí)現(xiàn)交會(huì)的轉(zhuǎn)移時(shí)間范圍。

2.3.2 仿真算例

追蹤航天器初始軌道根數(shù):長(zhǎng)半軸a=6 871 km,偏心率e=0.001,軌道傾角i=97.38°,升交點(diǎn)赤經(jīng)Ψ=60°,近地點(diǎn)幅角ω=20°,真近點(diǎn)角滿足f∈[0°,360°]。目標(biāo)航天器軌道根數(shù):長(zhǎng)半軸=7 171 km,偏心率et=0.01,軌道傾角=100°,升交點(diǎn)赤經(jīng)Ψt=55°,近地點(diǎn)幅角ωt=30°,真近點(diǎn)角=140°。Δt∈[600 s,5 400 s],t-st=25 s,f-st=1°,追蹤器提供的單次沖量上限Δvmax=3 km/s。

通過(guò)數(shù)值計(jì)算,作出雙沖量交會(huì)的Δv等值線圖如圖1所示。

圖1 Δv等值線圖Fig.1 Contour ofΔv

由圖1可直觀地得到滿足燃料限制的(Δt,f)組合范圍,以及滿足燃料、轉(zhuǎn)移時(shí)間限制的f范圍。該f的范圍即對(duì)應(yīng)滿足要求的追蹤器初始位置范圍,或稱之為追蹤區(qū)。通過(guò)追蹤區(qū)即可得到與之對(duì)應(yīng)的發(fā)射窗口。對(duì)于燃料限制、轉(zhuǎn)移時(shí)間限制和追蹤器初始位置3個(gè)量,已知任意2個(gè)量的取值,由圖1均可迅速確定第3個(gè)量的取值范圍。例如,Δt=1 800 s時(shí),滿足燃料限制的追蹤區(qū)為f∈[0°,76°]∪[358°,360°]。同理,若追蹤器初始位置一定,即f給定,亦可得到滿足燃料限制的轉(zhuǎn)移時(shí)間范圍。例如,f=60°時(shí),滿足需求的Δt∈[1 250 s,1 800 s]。

由仿真數(shù)據(jù)亦可得到任一Δt對(duì)應(yīng)交會(huì)問(wèn)題的燃料最優(yōu)數(shù)值解Δvmin,相應(yīng)的Δt-Δvmin曲線,以及滿足近似最優(yōu)解的Δt-f曲線,見(jiàn)圖2。

利用以上數(shù)值解及相應(yīng)的近似解曲線,不需進(jìn)行繁瑣低效的智能算法尋優(yōu),便可確定滿足任意給定初始狀態(tài)、燃耗及時(shí)間限制的最佳發(fā)射窗口,近似最小燃耗,以及燃料最優(yōu)交會(huì)軌道,即可迅速方便進(jìn)行軌道機(jī)動(dòng)任務(wù)規(guī)劃。整個(gè)計(jì)算過(guò)程耗時(shí)小于10 s。由本節(jié)仿真結(jié)果可知,在給定的限制條件下,最優(yōu)的Δvmin≈1.021 km/s,對(duì)應(yīng)的轉(zhuǎn)移時(shí)間Δt≈1 h,追蹤器初始真近點(diǎn)角f≈280°,該組結(jié)果為有限推力交會(huì)提供了必需的初始條件。

圖2 燃料最優(yōu)交會(huì)f、Δvmin與Δt的關(guān)系曲線Fig.2 f vs t andΔv vs tdeterm ined by numerical solutions

3 有限推力燃料最優(yōu)交會(huì)

3.1 動(dòng)力學(xué)模型

考慮地球扁率攝動(dòng),在服從平方反比律的中心引力場(chǎng)和有限推力假設(shè)下,以相對(duì)于地心赤道慣性系的位置、速度為狀態(tài)量的變軌運(yùn)動(dòng)動(dòng)力學(xué)方程為

式中 r=[x,y,z]T為航天器地心距矢量;v=[x·,y·,z·]為航天器速度矢量;f為軌道攝動(dòng)加速度矢量;μ為地球引力常數(shù);F為發(fā)動(dòng)機(jī)推力大小;m0為航天器初始質(zhì)量;|m·|為發(fā)動(dòng)機(jī)燃料秒耗量;t為飛行時(shí)間;u=[ux,uy,uz]T為推力方向在地心赤道慣性坐標(biāo)系中的單位矢量;Isp為發(fā)動(dòng)機(jī)比沖;gn為標(biāo)準(zhǔn)重力加速度值。

為計(jì)算需要,推力方向矢量u可用推力的2個(gè)控制方位角α和β表示,即

式中 α為俯仰控制角,定義為推力矢量與赤道面的夾角;β為偏航控制角,定義為推力矢量在赤道面的投影與ECI坐標(biāo)系X軸的夾角。

3.2 最優(yōu)控制

3.2.1 最優(yōu)控制數(shù)學(xué)模型

基于上述的任務(wù)規(guī)劃,最優(yōu)控制的目標(biāo)為尋求最優(yōu)的控制量α、β、F,使得經(jīng)過(guò)固定時(shí)間Δt實(shí)現(xiàn)平臺(tái)與目標(biāo)軌道上預(yù)定目標(biāo)點(diǎn)軟交會(huì),且燃料消耗達(dá)到最小。

變軌過(guò)程發(fā)動(dòng)機(jī)推力連續(xù)可變,t0和tf為工作起始和終止時(shí)間,均可由上述任務(wù)規(guī)劃獲得。不失一般性,設(shè)初始變軌時(shí)刻t0=0。交會(huì)過(guò)程中燃料消耗最少等價(jià)于總沖量最小,性能指標(biāo)為

同時(shí),滿足狀態(tài)變量邊界約束:

滿足控制變量約束:

3.2.2 最優(yōu)控制必要條件

引入哈密爾頓函數(shù):

根據(jù)極小值原理,最優(yōu)推力方向應(yīng)使H達(dá)到最小,同時(shí)滿足uTu=1,得到

伴隨方程為

3.3 數(shù)值求解

3.3.1 約束處理

對(duì)終端條件約束,懲罰函數(shù)法(SUMT方法)是較有效的方法之一。由于靜態(tài)罰函數(shù)適應(yīng)性較差,懲罰效率不高,利用動(dòng)態(tài)SUMT方法進(jìn)行約束處理。令X=[x,y,z,x·,y·,z·]T,引入罰函數(shù):

則性能指標(biāo)變?yōu)?/p>

式中 ki為罰因子,其取值原則應(yīng)使每一個(gè)罰函數(shù)都與積分項(xiàng)相當(dāng)。

借鑒模擬退火思想,設(shè)計(jì)了具有自適應(yīng)性的退火懲罰函數(shù)[10-11]。令ki=1/Ti,其中:

3.3.2 橫截條件

3.3.3 共軛梯度法

共軛梯度法是一種求解最優(yōu)控制問(wèn)題的行之有效的方法,通常采用以下公式:

其中,βk為共軛系數(shù),其選取對(duì)應(yīng)了不同的共軛梯度法,PRP方法是目前認(rèn)為數(shù)值表現(xiàn)最好的共軛梯度方法之一,故選取該方法,有=(gk-gk-1,gk)/(gk-1,gk-1),(＊,＊)表示矢量?jī)?nèi)積。

為了避免對(duì)無(wú)物理意義的協(xié)態(tài)量初值進(jìn)行猜測(cè)而引起的初值敏感問(wèn)題,文中選擇對(duì)控制量初值進(jìn)行猜測(cè),采用共軛梯度法進(jìn)行求解。具體步驟如下:

(1)在容許控制集中適當(dāng)猜測(cè)初始控制uk。置k=0;設(shè)定迭代計(jì)算精度ε1、ε2。

(2)用uk和初始條件x(t0)=x0從t0到tf正向積分狀態(tài)方程,得到xk。利用橫截條件λ()=λf,從tf到t0反向積分共軛方程,得到λk。

(3)用uk及計(jì)算得到的xk和λk計(jì)算梯度gk=

如果‖gk≤ε1‖,停止;否則,進(jìn)行下一步。

(5)用一維精確線性搜索確定最優(yōu)搜索步長(zhǎng)αk＞0,使得

(6)引入控制約束算子

3.4 仿真算例

利用前述任務(wù)規(guī)劃的數(shù)值結(jié)果知,追蹤器初始真近點(diǎn)角f=280°,其余軌道根數(shù)與2.3.2節(jié)相同。追蹤器初始質(zhì)量為3 000 kg,發(fā)動(dòng)機(jī)比沖300 s。交會(huì)時(shí)間3 600 s,交會(huì)精度要求:相對(duì)位置±1 km,相對(duì)速度±1 m/s?？刂屏砍踔颠x為α=0°,β=0°,F=1 200 N。仿真曲線見(jiàn)圖3～圖6。

由仿真結(jié)果可知,滿足交會(huì)精度要求時(shí),消耗特征速度1.976 km/s,追蹤器與目標(biāo)的相對(duì)距離為916.4 m,相對(duì)速度0.999 8 m/s,在滿足各項(xiàng)限制條件的同時(shí),達(dá)到了預(yù)設(shè)的遠(yuǎn)程交會(huì)技術(shù)指標(biāo)。有限推力交會(huì)消耗的特征速度與雙沖量交會(huì)相比,增加了92.9%。主要原因有:

(1)沖量交會(huì)假設(shè)推力無(wú)限大,推力作用是瞬間完成的,而實(shí)際的有限推力持續(xù)時(shí)間較長(zhǎng),存在較大的重力損失。

(2)由于該仿真中,發(fā)動(dòng)機(jī)推力較小,將沖量交會(huì)時(shí)間作為有限推力交會(huì)時(shí)間,存在一定誤差,進(jìn)而產(chǎn)生特征速度損失。但將沖量變軌作為初步模型進(jìn)行空間交會(huì)任務(wù)規(guī)劃,可迅速高效地為有限推力最優(yōu)變軌提供重要的參考數(shù)據(jù),這對(duì)注重任務(wù)規(guī)劃時(shí)間的情況(如空間打擊)而言,具有積極而重要的意義。

同時(shí),文中采用的數(shù)值求解技術(shù),降低了初值猜測(cè)難度,提高了求解兩點(diǎn)邊值問(wèn)題的效率,所得最優(yōu)解平滑,精度較高。但由于間接法求解最優(yōu)化問(wèn)題的固有特性,文中方法仍無(wú)法完全避免初值敏感問(wèn)題。

圖3 推力方向變化曲線Fig.3 Curves of thrust direction vs time

圖4 推力大小變化曲線Fig.4 Curves of thrustm agnitude vs time

圖5 相對(duì)距離變化曲線Fig.5 Curves of relative distance vs time

圖6 相對(duì)速度變化曲線Fig.6 Curves of relative velocity vs time

4 結(jié)論

(1)采用B-V算法求解高斯問(wèn)題,不需對(duì)圓錐曲線類型進(jìn)行討論,計(jì)算簡(jiǎn)潔,移去了轉(zhuǎn)移角為π的奇異性,并大幅度改善了算法對(duì)整個(gè)轉(zhuǎn)移角范圍的收斂性。

(2)對(duì)空間交會(huì)等任務(wù)分析而言,將繪制和分析等值線圖作為任務(wù)規(guī)劃的輔助工具,是普遍適用和有效的,能直觀快速地為工程應(yīng)用提供參考。

(3)將沖量變軌作為初步模型進(jìn)行空間交會(huì)任務(wù)規(guī)劃,可迅速高效地為有限推力最優(yōu)變軌提供重要的參考數(shù)據(jù)。

(4)將有限推力交會(huì)最優(yōu)控制伴隨變量初值猜測(cè)轉(zhuǎn)化為對(duì)物理意義明顯的控制量初值的猜測(cè),極大地降低了求解兩點(diǎn)邊值問(wèn)題的難度,所得最優(yōu)解平滑,精度較高,可作為工程應(yīng)用的相關(guān)參考。

[1] Prussing JE,Chiu JH.Optimal multiple-impu lsetime-fixed rendezvous between circu lar orbits[R].AIAA-84-2036.

[2] 張立佳,郭繼峰,崔乃剛.伴飛航天器燃料最優(yōu)交會(huì)[J].宇航學(xué)報(bào),2007,28(4):870-874.

[3] 陳統(tǒng),徐世杰.基于遺傳算法的最優(yōu)Lambert雙脈沖轉(zhuǎn)移[J].北京航空航天大學(xué)學(xué)報(bào),2007,33(3):273-277.

[4] Richard H Battin,Robin M Vaughan.An elegant lambertalgorithm[J].J.guidance,7(6):662-670.

[5] Haijun Shen,Panagiotis Tsiotras.Optimal two-impulse rendezvous using multip le-revolution lambert solutions[J].Journalof Guidance,Control and Dynamics,2003,26(1):50-61.

[6] John TBetts.Survey of numerical methods for trajectory optimization[J].Journal of Guidance,Control and Dynamics,1998,21(2):193-207.

[7] Gao Y,Kluever C A.Low-thrust interplanetary orbit transfers using hybrid trajectory optimization method with multiple shooting[R].AIAA 2004-5088.

[8] Paul J Enright,Brace A Conwayt.Optimal finite-thrust spacecraft trajectories using collocation and non linear programming[J].J.Guidance,1991,14(5):981-985.

[9] Yuri Ulybyshev.Continuous thrustorbit transfer optimization using large-scale linear programming[J].Journal of Guidance,Control and Dynamics,2007,30(2):427-436.

[10] 朱建豐,徐世杰.基于自適應(yīng)模擬退火遺傳算法的月球軟著陸軌道優(yōu)化[J].航空學(xué)報(bào),2007,28(4):806-812.

[11] 任遠(yuǎn),崔平遠(yuǎn),欒恩杰.基于退火遺傳算法的小推力軌道優(yōu)化問(wèn)題研究[J].宇航學(xué)報(bào),2007,28(1):162-166.