利用虛擬衛(wèi)星法求解火星探測器近火點(diǎn)制動策略

劉 玥，荊武興

(哈爾濱工業(yè)大學(xué) 航天工程系，150001 哈爾濱)

從上個(gè)世紀(jì)60年代開始，火星逐漸成為了太陽系內(nèi)除月球之外最吸引各國學(xué)者研究的天體之一.最近的幾年，歐洲和美國先后發(fā)射了多顆火星探測器，在全球又掀起了新一輪的火星探測競賽，中國、印度等國家也相繼提出了自己的火星探測計(jì)劃.火星是太陽系內(nèi)與地球最相似的行星之一，但其質(zhì)量較小，引力系數(shù)只有地球的10.7%，因此火星表面逃逸速度也只有不到5 km/s.傳統(tǒng)的火星探測方式使用Hohmann 過渡的方法將探測器發(fā)射至火星附近，其雙曲進(jìn)入軌道剩余速度達(dá)到了2.8 km/s 以上，探測器必須克服這一剩余速度使自身成為火星衛(wèi)星［1－3］.若如此，則需要發(fā)動機(jī)長時(shí)間開機(jī)進(jìn)行制動，很大一部分燃料將消耗在制動過程中［4－5］.另外，火星探測任務(wù)需要預(yù)定的工作軌道，制動過程應(yīng)盡量將飛行器精確送入預(yù)定工作軌道，否則還將進(jìn)行復(fù)雜的變軌過程，進(jìn)一步消耗燃料.為此，需要優(yōu)化近火點(diǎn)制動推力策略，在節(jié)約探測器燃料的同時(shí)，使探測器精確進(jìn)入目標(biāo)軌道，提高了探測器在軌壽命，使后續(xù)任務(wù)得以順利完成［6－7］.

火星探測器近火點(diǎn)制動問題可以歸結(jié)為1 個(gè)兩點(diǎn)邊值最優(yōu)控制問題，對于此類問題的解決，一般采用極大值原理.極大值原理的實(shí)質(zhì)是求取1 個(gè)控制策略，使系統(tǒng)的Hamilton 函數(shù)達(dá)到極大值，這種方法的優(yōu)點(diǎn)在于可以精確求取連續(xù)時(shí)間內(nèi)的最優(yōu)控制策略［8－9］.在軌道控制方面，極大值原理被廣泛應(yīng)用于燃料、能量最優(yōu)軌道轉(zhuǎn)移、軌道修正與保持，交會對接等軌道控制任務(wù)中［10－12］.但是其極大值原理存在一些問題，就是需要計(jì)算共軛變量，并且共軛變量的初值需要猜測和迭代計(jì)算.由于其沒有具體的物理意義，共軛變量初值的猜測往往很困難，容易造成迭代不收斂，降低了極大值原理求解問題的效率.

本文利用虛擬衛(wèi)星方法，以最省燃料為優(yōu)化指標(biāo)，求解火星探測器近火點(diǎn)制動的推力策略，并且在使用極大值原理求解兩點(diǎn)邊值問題時(shí)，將共軛變量轉(zhuǎn)化為有實(shí)際物理意義的變量進(jìn)行猜測，使得迭代初值的選取變得簡單［13］.計(jì)算結(jié)果顯示，該策略可以在最省燃料的條件下將探測器精確送入預(yù)定工作軌道，并且邊界條件非常簡單，迭代易于收斂.

1 虛擬衛(wèi)星法原理與模型

本文使用虛擬衛(wèi)星法求解最優(yōu)制動策略.如圖1 所示，軌道I 為探測器初始軌道，S 為探測器，軌道II 為目標(biāo)軌道.假設(shè)軌道II 是1 顆虛擬衛(wèi)星S1運(yùn)動形成的1 條軌道，為了使探測器轉(zhuǎn)移到目標(biāo)軌道，只需要求解控制策略，使S 與S1的相對位置和速度均為0 即可.此方法稱為虛擬衛(wèi)星法［14-17］.

圖1 虛擬衛(wèi)星法示意

圖1 中，建立交會坐標(biāo)系S1-xy，以S1為坐標(biāo)原點(diǎn)，S1-y 軸指向火心，S1-x 軸垂直于y 軸并指向前進(jìn)方向，由于變軌時(shí)間只有數(shù)十分鐘，所以可忽略攝動力的影響，S 和S1的動力學(xué)方程分別為

其中:R、r 分別為S1和S 相對于火心的位置;μm為火星的萬有引力常數(shù);F 為發(fā)動機(jī)推力;m 為探測器總質(zhì)量.

假設(shè)S1的真近點(diǎn)角為θ1，與x 軸同向的單位向量為I，與y 軸同向的單位向量為J，S 相對于S1的位置矢量為ρ，S 在S1-xy 系中的坐標(biāo)為(x，y)，則有如下關(guān)系:

由旋轉(zhuǎn)矢量的微分公式

將R=-RJ 帶入得到分量形式的動力學(xué)方程

式中Fx、Fy為推力F 的分量，eI、pI、hI分別是目標(biāo)偏心率、目標(biāo)半正焦弦和目標(biāo)軌道角動量.

此模型的優(yōu)越性在于終端條件極其簡單，因?yàn)樘綔y器與虛擬衛(wèi)星的終端相對位置和速度只需為零.

2 探測器狀態(tài)方程

假設(shè)在［t0，tf］內(nèi)發(fā)動機(jī)持續(xù)工作，探測器質(zhì)量變化即成為已知函數(shù)

探測器動力學(xué)方程可以用狀態(tài)方程的形式表達(dá)如下:

下面推導(dǎo)初始狀態(tài)與真近角的關(guān)系.設(shè)A 點(diǎn)在初始軌道上對應(yīng)的真近角為Ω，在目標(biāo)軌道上對應(yīng)的真近點(diǎn)角為ΩI，點(diǎn)火初始時(shí)刻，探測器的真近角為θ(t0)，虛擬衛(wèi)星的真近角為θI(t0)，令

稱α 為點(diǎn)火角，順軌道測量為正.不難得到

設(shè)真實(shí)衛(wèi)星和虛擬衛(wèi)星的飛行路徑角分別為γ 和γI，則在S1-xy 坐標(biāo)系上

利用如下關(guān)系:

得到相對速度的表達(dá)式

其中e，p，h 分別為軌道I 的偏心率、半正焦弦和軌道角動量;eI，pI，hI分別為軌道II 的偏心率、半正焦弦和軌道角動量;σ=θI-ΩI+Ω，ωI=hI/R2.

以上兩組式子確定了真近角與探測器的初始狀態(tài)之間的關(guān)系，任意確定一組出事的θ，θI，即可確定唯一與之對應(yīng)的一組初始狀態(tài)x，y，，.

3 最省燃料推力方向

最省燃料控制問題可表述為:求取合適的X(t0)，V(t0)，θI(t0)和推力方向策略u(t)，使得系統(tǒng)在t=tf時(shí)，X(tf)=V(tf)=0，θI無約束.且使最優(yōu)指標(biāo)

達(dá)到最大值.由Pontryagin 極大值原理，此指標(biāo)極大與探測器剩余燃料質(zhì)量極大是等價(jià)的.

當(dāng)X(t0)，V(t0)，θI(t0)均為給定值時(shí)，上述最優(yōu)控制問題即為1 個(gè)兩點(diǎn)邊值問題，但是其初值是由探測器與虛擬衛(wèi)星的真近點(diǎn)角θ 和θI共同確定的，而事實(shí)上，探測器初始點(diǎn)火位置的選擇是任意的，即初始真近點(diǎn)角θ 可以取任意值，而虛擬衛(wèi)星初始位置同樣可以任取.所以，最優(yōu)指標(biāo)可以表達(dá)為

θ(t0)，θI(t0)允許在(-π，π)上變動，于是求解上述問題將分為兩部分，最優(yōu)控制的求取與參數(shù)優(yōu)化，應(yīng)用極大值原理，可以將其統(tǒng)一為1 個(gè)參數(shù)優(yōu)化問題.

系統(tǒng)的Hamilton 函數(shù)定義為

其中λx，λv，λθI分別為X，V，θI的共軛狀態(tài).

由Pontryagin 極大值原理，求解的推力矢量方向u(t)應(yīng)使Hamilton 函數(shù)H 達(dá)到極大，可得

可見，速度共軛矢量λv的方向可以決定最優(yōu)推力方向.

系統(tǒng)協(xié)狀態(tài)方程如下:

如果已知θ(t0)，θI(t0)，λx(t0)，λv(t0)，λθI(t0)，則最優(yōu)推力矢量方向u(t)就可以唯一確定，因此只需求解參數(shù)優(yōu)化問題即可解決上述最優(yōu)控制問題.

4 橫截條件與迭代變量的選擇

前一節(jié)指出，需要確定5 個(gè)參數(shù)即可求出最優(yōu)推力方向，其中θ(t0)，θI(t0)具有實(shí)際的物理意義，但是另外3 個(gè)參量沒有實(shí)際意義，很難猜測初始值，這是Pontryagin 極大值原理固有的缺陷，下面將給出在本文背景下解決此缺陷的方法.

由前文可知，只要已知x，θI，另外的初始狀態(tài)就可以確定，因此這些初始狀態(tài)并不是相互獨(dú)立的，還應(yīng)存在3 個(gè)約束條件.由式(2)，Δα 可表示為θI和x 的函數(shù)，從式(1)和(2)可以得出3 個(gè)初始狀態(tài)約束條件:

終端條件為

由上述條件，根據(jù)變分法原理，可以得到共軛協(xié)狀態(tài)的初始和終端條件為

其中k1，k2，k3為待定拉格朗日乘子.則共軛變量初值的猜測可轉(zhuǎn)化為對拉格朗日乘子初值的猜測，但這些變量初值仍然難以猜測，所以要進(jìn)一步分析.

由式(4)、(7)可以看出，k2，k3與初始推力方向有如下關(guān)系:

如用推力仰角φ 表示U(規(guī)定推力方向與oy軸成銳角時(shí)為負(fù))，則

如果不考慮質(zhì)量變化，由極大值原理，在初始時(shí)刻應(yīng)有‖λv(t0)‖=1，實(shí)驗(yàn)證實(shí)，考慮質(zhì)量變化時(shí)將有較小的變化，一般不超過± 1，所以‖λv(t0)‖和具有物理意義的推力仰角初值φ 可以代替k2，k3作為迭代初值.

由式(5)可得

由于ωI的數(shù)值比較小，代入?yún)f(xié)狀態(tài)的表達(dá)式可以發(fā)現(xiàn)，推力方向變化率的初始值主要由k1決定，一般來說，希望推力方向(與λv(t0)有關(guān))變化較平穩(wěn)，所以可取k1≈0 作為初始猜測.

綜上所述，迭代變量的初始猜測值可選擇為θ(t0)，θI(t0)，φ，‖λv(t0)‖和k1.

在最優(yōu)軌跡末端，哈密頓函數(shù)應(yīng)滿足:H(tf)=0，由式(3)、(6)、(7)可得到‖λv(tf)‖=1.至此可以看出，5 個(gè)迭代變量有3 個(gè)是具有實(shí)際意義的量，另外兩個(gè)可以有較為確定的猜測值，為迭代算法的收斂性打下了良好的基礎(chǔ).

5 仿真算例分析

考慮火星探測器初始質(zhì)量為500 kg，攜帶490 N 發(fā)動機(jī)一臺，比沖為290 s.假設(shè)探測器預(yù)定工作軌道為一共面橢圓軌道，近火點(diǎn)高度為300 km，偏心率為0.5.迭代初值選擇如表1 所示.

表1 迭代初值選擇

圖2 迭代殘差變化

迭代最終結(jié)果顯示，發(fā)動機(jī)工作1 253 s，探測器質(zhì)量m，真近角θ，虛擬衛(wèi)星真近角θI，探測器軌道偏心率e 和近火點(diǎn)高度Hp的計(jì)算結(jié)果列于表2.

表2 計(jì)算結(jié)果

真實(shí)衛(wèi)星與虛擬衛(wèi)星的終端位置誤差為0.37 m，速度誤差約為10－3m/s.可見，使用本文中的算法可以達(dá)到將探測器一次精確送入工作軌道的目的.推力仰角隨時(shí)間的變化規(guī)律如圖3 所示，探測器與虛擬衛(wèi)星在慣性系下的軌跡如圖4所示.

圖3 推力仰角變化規(guī)律

圖4 探測器飛行軌跡

由圖3 可以看到，發(fā)動機(jī)推力仰角隨時(shí)間的變化近似為線性的，變軌過程中姿態(tài)控制比較容易.由圖4 可以看到，探測器在最優(yōu)推力的控制下，其軌跡最終與虛擬衛(wèi)星軌跡重合，兩者實(shí)現(xiàn)了軟交會，即完成了精確入軌過程.

綜上所述，利用虛擬衛(wèi)星法求解火星探測器近火點(diǎn)制動推力策略是可行的.在實(shí)際計(jì)算過程中，如果考慮各種攝動因素，如大型天體引力攝動，火星非球形攝動等，結(jié)果的收斂性也很好，迭代均能在10 次以內(nèi)收斂，說明該方法對于模型的誤差有一定的魯棒性.但是，此方法仍然有局限性.例如，由于假設(shè)發(fā)動機(jī)一直開機(jī)，所以一般情況下，該方法只能處理1 次入軌的情況，若要處理多次入軌的情況，則需要設(shè)定一些過渡軌道，多次變軌，其間分別使用虛擬衛(wèi)星方法;另外，虛擬衛(wèi)星方法的收斂性與發(fā)動機(jī)推進(jìn)效果相關(guān)，若換裝更小推力的發(fā)動機(jī)，或者增加飛行器的初始質(zhì)量，則會導(dǎo)致制動時(shí)間增加、收斂范圍變小.經(jīng)大量仿真驗(yàn)證，對于本文中所使用的490 N 發(fā)動機(jī)，虛擬衛(wèi)星方法可以滿足初始質(zhì)量在1 200 kg 以內(nèi)的飛行器進(jìn)行近火點(diǎn)制動，且迭代收斂性較好.

5 結(jié)論

本文利用虛擬衛(wèi)星法設(shè)計(jì)火星探測器近火點(diǎn)制動的最優(yōu)推力方向策略.采用Pontryagin 極大值原理求取燃料最優(yōu)指標(biāo)下的最優(yōu)推力方向，并且給出了協(xié)狀態(tài)的選擇方法以使迭代易于收斂.通過優(yōu)化的制動策略使真實(shí)衛(wèi)星與虛擬衛(wèi)星實(shí)現(xiàn)“軟交會”，達(dá)到精確制動入軌的目的.所選用的虛擬衛(wèi)星方法收斂性好，對模型具有一定的魯棒性，對迭代初值不敏感，可以適用于一般的常規(guī)推力近心點(diǎn)制動任務(wù).

［1］李爽，彭玉明，陸宇平.火星EDL 導(dǎo)航、制導(dǎo)與控制技術(shù)綜述與展望［J］.宇航學(xué)報(bào)，2010，3(3):621－627.

［2］HUBBARD G S.The exploration of Mars，historical context ＆ current results［C］//42nd AIAA Aerospace Sciences Meeting and Exhibit.Reno，NV:American Institute of Aeronautics and Astronautics，2004:1－10.

［3］OLESON S R，MCGUIRE M L，LAURA B.Mars earth return vehicle (MERV)propulsion options［C］//46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference＆ Exhibit.Nashville，TN:American Institute of Aeronautics and Astronautics，2010:1－20.

［4］OZIMEK M T，HOWELL K C.Low-thrust transfers in the earth-moon system，including applications to libration point orbits［J］.Journal of Guidance，Control，and Dynamics，2010，33(2):533－549.

［5］WILLIAMS P.Optimal control of electro dynamic tether orbit transfers using timescale deparation［J］.Journal of Guidance，Control，and Dynamics，2010，33(1):88－97.

［6］ULYBYSHEV Y.Continuous thrust orbit transfer optimization using large-scale linear programming［J］.Journal of Guidance，Control，and Dynamics，2007，30(2):427－436.

［7］JAMISON B R，COVERSTONE V.Analytical study of the primer vector and orbit transfer switching function［J］.Journal of Guidance，Control，and Dynamics，2010，33(1):235－245.

［8］CHEN Yumin，SHEU Donglong.Parametric optimization analysis for minimum-fuel low-thrust coplanar orbit transfer ［J］.Journal of Guidance，Control，and Dynamics，2006，29 (6):1446－1451.

［9］BANDO M，ICHIKAWA A.Periodic orbits of nonlinear relative dynamics along an eccentric orbit［J］.Journal of Guidance，Control，and Dynamics，2010，33(2):385－395.

［10］王小軍，吳德隆，余夢倫.最少燃料消耗的固定推力共面軌道變軌研究［J］.宇航學(xué)報(bào)，1995，16(4):9－16.

［11］岳新成，楊瑩，耿志勇.有限推力能量、燃料最優(yōu)軌道轉(zhuǎn)移控制［J］，南京航空航天大學(xué)學(xué)報(bào)，2009，41(6):715－719.

［12］梁新剛，楊滌.應(yīng)用非線性規(guī)劃求解異面最優(yōu)軌道轉(zhuǎn)移問題［J］.宇航學(xué)報(bào)，2006，27(3):365－370.

［13］荊武興，吳瑤華.衛(wèi)星軌道圓化的有限推力點(diǎn)火控制策略［J］.宇航學(xué)報(bào)，1997，4(4):1－6.

［14］荊武興，吳瑤華，楊滌.基于交會概念的最省燃料共面有限推力軌道轉(zhuǎn)移方法［J］.哈爾濱工業(yè)大學(xué)學(xué)報(bào)，1997，4(4):132－135.

［15］王明春，荊武興，楊滌，等.能量最省有限推力同平面軌道轉(zhuǎn)移［J］.宇航學(xué)報(bào)，1992，3(3):24－31.

［16］荊武興，吳瑤華.基于交會概念的最省燃料異面有限推力軌道轉(zhuǎn)移研究［J］.哈爾濱工業(yè)大學(xué)學(xué)報(bào)，1998，30(2):124－128.

［17］劉玥.火星探測器近心點(diǎn)制動與軌道保持優(yōu)化設(shè)計(jì)［D］.哈爾濱:哈爾濱工業(yè)大學(xué)，2011:18－29.