劉柯邑，岳曉奎，張 瑩

（西北工業(yè)大學(xué) 航天學(xué)院，陜西 西安 710072）

小推力轉(zhuǎn)移軌道設(shè)計(jì)與優(yōu)化問題逐步成為研究的熱點(diǎn)[1-2]，與傳統(tǒng)的火箭發(fā)動(dòng)機(jī)相比，小推力發(fā)動(dòng)機(jī)具有比沖高、質(zhì)量輕的特點(diǎn)，可有效減少燃料消耗，提高任務(wù)的有效載荷，降低發(fā)射成本。但由于其推力?。ㄍǔ閙N量級(jí)）、作用時(shí)間長(zhǎng)導(dǎo)致小推力軌道的動(dòng)力學(xué)非線性嚴(yán)重，這給小推力軌道的設(shè)計(jì)帶來(lái)了難題。小推力轉(zhuǎn)移軌道設(shè)計(jì)問題可歸結(jié)為求解有約束的連續(xù)最優(yōu)控制問題。學(xué)者們針對(duì)該問題所采用的方法包括間接法和直接法[3-4]。

間接法基于變分法將最優(yōu)控制問題最終轉(zhuǎn)換為一個(gè)兩點(diǎn)邊值問題，但這種方法收斂半徑小，對(duì)協(xié)狀態(tài)的初始猜測(cè)值很敏感，對(duì)于復(fù)雜約束難以處理。相對(duì)間接法，直接法在收斂的魯棒性和解決實(shí)際問題的適應(yīng)性上具有優(yōu)勢(shì)。直接法是將最優(yōu)控制問題離散成參數(shù)優(yōu)化問題，然后利用非線性規(guī)劃方法求解。高斯偽譜方法[5-6]是最近研究較多的直接配點(diǎn)優(yōu)化方法，相對(duì)于其他直接配點(diǎn)法[7]的優(yōu)勢(shì)是可以用較少的節(jié)點(diǎn)獲得較高的精度。麻省理工學(xué)院的Benson從理論上證明了高斯偽譜法的KKT條件與最優(yōu)控制理論中的一階必要條件是一致的；并且，由于其無(wú)需猜測(cè)協(xié)狀態(tài)的初始值，這就大大降低了求解最優(yōu)控制問題的難度。

文中研究了基于高斯偽譜法的燃料最省小推力轉(zhuǎn)移軌道設(shè)計(jì)問題。首先建立了星際小推力軌道優(yōu)化模型，對(duì)轉(zhuǎn)移軌道參數(shù)進(jìn)行了無(wú)量綱化處理以提高數(shù)值計(jì)算精度。然后利用高斯偽譜方法將轉(zhuǎn)移軌道優(yōu)化問題轉(zhuǎn)化為多約束的非線性規(guī)劃問題，并應(yīng)用SNOPT優(yōu)化軟件包求解了最優(yōu)轉(zhuǎn)移軌道的推力控制角變化規(guī)律。

1 小推力軌道優(yōu)化問題

只考慮探測(cè)器在星際巡航段的轉(zhuǎn)移軌道設(shè)計(jì)問題，采用笛卡爾坐標(biāo)系下的小推力探測(cè)器動(dòng)力學(xué)模型：

其中，r=[x，y，z]T為探測(cè)器在日心黃道慣性系中的位置矢量，v=[vx，vy，vz]T為發(fā)動(dòng)機(jī)推力在各個(gè)方向上的加速度矢量，m為飛行器和推進(jìn)劑的總質(zhì)量； 為太陽(yáng)引力系數(shù)； 為推力幅值；α=[αx，αy，αz]T為發(fā)動(dòng)機(jī)推力在各個(gè)方向的加速度矢量；g0為海平面重力加速度；Isp為發(fā)動(dòng)機(jī)比沖。

假定小推力發(fā)動(dòng)機(jī)的輸入功率P隨著探測(cè)器日心距離的增大而減小，它與探測(cè)器日心距離的平方成反比，即

其中，P0為探測(cè)器日心距離為1AU時(shí)的輸入功率，r采用天文單位。發(fā)動(dòng)機(jī)推力T為：

其中，η為推進(jìn)器的工作效率。

在推進(jìn)器連續(xù)工作條件下，探測(cè)器從地球發(fā)射與目標(biāo)星體的燃料最省轉(zhuǎn)移軌道的性能指標(biāo)可以寫成：

探測(cè)器的發(fā)射時(shí)間和到達(dá)時(shí)間均不固定，對(duì)于任何的發(fā)射時(shí)間t0，探測(cè)器受到地球星歷約束，探測(cè)器的位置和速度需要滿足下面的等式約束

其中，rsc，vsc是發(fā)射時(shí)刻探測(cè)器在日心黃道坐標(biāo)系中的位置和速度；rE，vE是地球的位置和速度；v∞是地球逃逸雙曲線超速，其大小已知，方向不確定。由公式（4），（5）構(gòu)成了探測(cè)器軌道的初始邊界條件：

當(dāng)探測(cè)器與目標(biāo)星體交會(huì)時(shí)，探測(cè)器的位置和速度需要滿足下面的等式約束：

其中，tf是探測(cè)器與目標(biāo)星體的交會(huì)時(shí)間；rt（tf），vt（tf）是探測(cè)器交會(huì)時(shí)目標(biāo)星體的位置和速度。由公式（7）、（8）構(gòu)成了探測(cè)器軌道的末端邊界條件：

在轉(zhuǎn)移過程中，推力方向矢量必須滿足如下條件：

則星際探測(cè)任務(wù)的燃料最省轉(zhuǎn)移軌道可以描述為：尋求最優(yōu)的發(fā)射時(shí)間t0，連續(xù)的軌道狀態(tài) xsc（t），連續(xù)的推力方向矢量α（t）和交會(huì)時(shí)間tf，使探測(cè)器軌道在滿足動(dòng)力學(xué)約束（1），邊界約束（6）和（9），推力方向矢量約束（10）的條件下，使性能指標(biāo)（3）達(dá)到最小。

尋優(yōu)參數(shù)可以表示為：

為提高計(jì)算精度，避免利用高斯偽譜方法求解過程中出現(xiàn)病態(tài)矩陣，需要對(duì)軌道參數(shù)進(jìn)行無(wú)量綱化處理，使各狀態(tài)參數(shù)的取值范圍都在一個(gè)相近的數(shù)值范圍內(nèi)。為此，文中取量綱質(zhì)量單位[M]=m0；量綱距離單位[L]=AU；量綱時(shí)間單位[T]=。下面將采用高斯偽譜法對(duì)上述小推力軌道優(yōu)化問題進(jìn)行處理。

2 高斯偽譜方法

高斯偽譜法的基本求解思路為：將未知的軌道狀態(tài)和控制時(shí)間歷程在一系列高斯點(diǎn)上離散化，而后用這些離散的狀態(tài)與控制分別構(gòu)造Lagrange插值多項(xiàng)式去逼近真實(shí)的軌道狀態(tài)與控制，再通過對(duì)狀態(tài)量求導(dǎo)來(lái)代替動(dòng)力學(xué)微分方程。

由于高斯點(diǎn)分布在（-1，1）區(qū)間，若要把動(dòng)力學(xué)方程在高斯點(diǎn)上進(jìn)行離散，首先需要把所研究問題的時(shí)間區(qū)間t∈[t0，tf]轉(zhuǎn)換到τ∈[-1，1]，這個(gè)轉(zhuǎn)化可以通過下式完成：

高斯偽譜法取 N 個(gè)高斯點(diǎn) τ1，τ2，…，τN和初始端點(diǎn) τ0=-1上的離散狀態(tài)構(gòu)造Lagrange插值多項(xiàng)式去近似狀態(tài)的時(shí)間歷程，即

式中，x（τ）為真實(shí)的狀態(tài)時(shí)間歷程，X（τ）為由 Lagrange插值多項(xiàng)式近似得到的狀態(tài)時(shí)間歷程；Li（τ）為L(zhǎng)agrange插值基函數(shù)，i=0，1，…，N。

類似地，控制量的時(shí)間歷程也可以用離散的控制構(gòu)造Lagrange插值多項(xiàng)式來(lái)近似，即

式中，u（τ）為真實(shí)的控制時(shí)間歷程，U（τ）為由 Lagrange插值多項(xiàng)式近似得到的控制時(shí)間歷程；（τ）為L(zhǎng)agrange插值基函數(shù)，i=1，…，N。

需要注意的是，構(gòu)造狀態(tài)用了N個(gè)高斯點(diǎn)與初始端點(diǎn)共N+1個(gè)點(diǎn)，而構(gòu)造控制只用了N個(gè)高斯點(diǎn)，因而方程（13）的腳標(biāo)是從0開始，而方程（15）的腳標(biāo)卻從1開始。

式（13）對(duì)時(shí)間求導(dǎo)，得

定義微分近似矩陣D，其中元素為：

式中，D 為 N×（N+1）二維矩陣，k=1，…，N，i=0，1，…，N。

由式（17）、（18）將軌道動(dòng)力學(xué)約束轉(zhuǎn)化為一系列代數(shù)約束：

式中，Xk=X（τk），Uk=U（τk）；k=1，2，…，N， f為軌道動(dòng)力學(xué)方程右端函數(shù)。

終端狀態(tài)由高斯求積公式得到：

式中，X0=X（-1），ωk為高斯積分權(quán)重，有

式（19）將軌道動(dòng)力學(xué)約束轉(zhuǎn)化為只與節(jié)點(diǎn)處離散狀態(tài)和控制變量相關(guān)的代數(shù)約束，同樣，性能指標(biāo)和約束條件也可以處理為只與節(jié)點(diǎn)處狀態(tài)和控制變量相關(guān)的形式。

如此，星際小推力軌道優(yōu)化問題可以轉(zhuǎn)化為一個(gè)多約束的參數(shù)優(yōu)化問題，尋優(yōu)參數(shù)為發(fā)射時(shí)間、交會(huì)時(shí)間以及各節(jié)點(diǎn)上的狀態(tài)和控制變量：

性能指標(biāo)可以表示為：

始末端邊界條件可以表示為：

推力矢量約束可以表示為：

離散后的小推力軌道最優(yōu)控制問題為典型的非線性規(guī)劃問題，本文采用TOMLAB/PROPT優(yōu)化軟件包進(jìn)行求解。

3 數(shù)值仿真與分析

為了驗(yàn)證本文所提設(shè)計(jì)方法的正確性和有效性，對(duì)小推力探測(cè)器的星際轉(zhuǎn)移軌道進(jìn)行設(shè)計(jì)與優(yōu)化。探測(cè)器主要參數(shù)如表1所示。

表1 探測(cè)器主要參數(shù)表Tab.1 Parameters of space detector

算例一：地球-火星交會(huì)任務(wù)

采用直接打靶法和高斯偽譜法對(duì)地球-火星交會(huì)任務(wù)的燃料最省小推力轉(zhuǎn)移軌道分別進(jìn)行設(shè)計(jì)。

探測(cè)任務(wù)要求探測(cè)器在2009年1月1日到2009年12月31日之間從地球出發(fā)，出發(fā)時(shí)刻認(rèn)為探測(cè)器的日心位置和速度與地球相同。地球-火星交會(huì)的初始轉(zhuǎn)移時(shí)間設(shè)為365天。采用直接打靶法和高斯偽譜法優(yōu)化的軌道參數(shù)如表2所示。

由表2可以看出，采用直接打靶法和高斯偽譜法的優(yōu)化結(jié)果基本一致，這表明采用高斯偽譜法求解行星際小推力轉(zhuǎn)移軌道問題是有效地。采用直接打靶法求解該問題時(shí)，所用計(jì)算時(shí)間為1 712.343 8 s；采用高斯偽譜法時(shí)計(jì)算時(shí)間為56.140 6 s，這主要是由于直接打靶法在迭代過程中需要對(duì)狀態(tài)動(dòng)力學(xué)方程進(jìn)行積分造成的，同時(shí)也驗(yàn)證了高斯偽譜法在求解小推力轉(zhuǎn)移軌道問題上的優(yōu)勢(shì)。

表2 地球-火星軌道優(yōu)化結(jié)果對(duì)比Tab.2 Orbit optimization result of earth-mars

對(duì)于小推力軌道設(shè)計(jì)問題，通常采用俯仰和偏航兩個(gè)控制角來(lái)定義推力的指向。定義俯仰控制角α為推力矢量在軌道平面內(nèi)的投影與當(dāng)?shù)厮街g的夾角；偏航控制角β為推力矢量與密切軌道平面之間的夾角，則推力的單位方向矢量就可以表示為：a=[sinαcos β cosαcos β sin β]。 直接打靶法和高斯偽譜法的控制角隨時(shí)間的變化規(guī)律如圖1所示。

由于火星軌道面與地球軌道面之間的夾角很小，偏航控制角在整個(gè)飛行過程中均較小，最大時(shí)約為20°。探測(cè)器的推力方向隨飛行軌跡的變化如圖2所示。

圖1 直接打靶法和高斯偽譜法控制軌線對(duì)比Fig.1 Compare of orbital control trajectory

圖2 探測(cè)器飛行軌跡與推力方向Fig.2 Trajectory and thrust orientation

算例二：地球-金星交會(huì)任務(wù)

為了進(jìn)一步考察高斯偽譜法求解行小推力軌道優(yōu)化問題的有效性，采用高斯偽譜法對(duì)地球-金星交會(huì)任務(wù)的燃料最省小推力轉(zhuǎn)移軌道進(jìn)行設(shè)計(jì)。探測(cè)任務(wù)要求探測(cè)器在2021年1月1日到2021年12月31日之間從地球出發(fā)，出發(fā)時(shí)刻認(rèn)為探測(cè)器的日心位置和速度與地球相同，即探測(cè)器的地球逃逸速度為零。初始假設(shè)探測(cè)器飛行 天到達(dá)金星，此時(shí)要求探測(cè)器的日心速度等于金星的日心速度，與金星進(jìn)行交會(huì)。采用高斯偽譜法優(yōu)化的軌道參數(shù)如表3所示。

表3 地球-金星軌道優(yōu)化結(jié)果Tab.3 Orbit optimization result of earth-mars

表2和表3中轉(zhuǎn)移時(shí)間和剩余質(zhì)量的優(yōu)化結(jié)果對(duì)比，說明了發(fā)動(dòng)機(jī)開機(jī)時(shí)間長(zhǎng)短對(duì)太陽(yáng)能電推進(jìn)器的燃料消耗的影響不大；由于發(fā)動(dòng)機(jī)推力大小與探測(cè)器日心距離的平方成反比，相比于地球-火星段軌道，探測(cè)器地球-金星軌道段的推力比較大，所以剩余質(zhì)量相差不大。

探測(cè)器在地球-金星交會(huì)任務(wù)中的飛行軌跡如圖3所示。

圖3 探測(cè)器飛行軌跡Fig.3 Trajectory of space detector

圖4 地球-金星軌道控制軌線Fig.4 Control orbital trajectory of earth-venus

從飛行軌跡圖3可以看出，在轉(zhuǎn)移軌道初始階段，發(fā)動(dòng)機(jī)推力方向主要指向軌道內(nèi)以減小軌道的半長(zhǎng)軸；軌道中段，推力方向主要指向軌道面外以增加軌道傾角；在探測(cè)器接近金星的最后階段，推力方向主要指向軌道外以降低探測(cè)器的速度，避免撞擊金星。探測(cè)器在地球-金星交會(huì)任務(wù)中的發(fā)動(dòng)機(jī)推力矢量控制角隨時(shí)間的變化如圖4所示。由于金星軌道面與地球軌道面之間的夾角相對(duì)較大，在整個(gè)飛行過程中，偏航控制角最大時(shí)約為。

4 結(jié) 論

文中研究了基于高斯偽譜法的小推力轉(zhuǎn)移軌道設(shè)計(jì)與優(yōu)化問題。通過對(duì)燃料消耗最省火星和金星交會(huì)任務(wù)的軌道設(shè)計(jì)與分析表明：采用高斯偽譜法進(jìn)行求解時(shí)，仿真初值選取相對(duì)自由，并且收斂速度快，采用此優(yōu)化算法，可保證所得解最優(yōu)的基礎(chǔ)上，有效提高問題的計(jì)算效率，具有良好的應(yīng)用前景。

[1]Rayman M D，Williams S N.Design of the first interplanetary solar electric propulsion mission[J].Journal of Spacecraft and Rockets，2012，39（4）:589-595.

[2]Kugelberge J，Bodin P，Persson S，et al.Accommodating electric propulsion on SMART-1[J].Acta Astronautica，2004，55（2）:121-130.

[3]Ranieri C L，Ocampo C A.Indirect optimization of spiral trajectories[J].Journal of Guidance， Control and Dynamics，2006，29（6）:1360-1366.

[4]Betts J T，Erb S O.Optimal low thrust trajectories to the moon[J].SIAM Journal of Applied Dynamical Systems，2003，2 （2）:144-170.

[5]Benson D A.A gauss pseudospectral transcription for optimal control[D].Massachusetts:Department of Aeronautics and Astronautics， Massachusetts Institute of Technology，2005.

[6]Benson A，Thorvaldsen T，Rao V.Direct trajectory optimization and costae estimation via an orthogonal collocation method[J].Journal of Guidance， Control and Dynamics，2006，29 （6）:1435-1440.

[7]Betts J T.Very low-thrust trajectory optimization using a direct SQP method[J].Journal of Computation&Applied Mathmatics，2009，120（1）:27-40.