李羅剛，崔祜濤，馬春飛，毛春曉，荊武興

(哈爾濱工業(yè)大學(xué)深空探測基礎(chǔ)研究中心，哈爾濱150080，llg0315@sina.com)

隨著我國嫦娥探月計劃三步走“繞、落、回”的順利實施.月面返回技術(shù)成為我國月球探測計劃必須先期解決的問題.在我國，目前對月球探測器軟著陸最優(yōu)軌道設(shè)計的工作比較多［1－3］，但對月面返回最優(yōu)軌道設(shè)計還是一片空白.而在國外，早在上個世紀60到70年代的阿波羅登月計劃中，就已經(jīng)成功將航天員送上月球并返回，在這個領(lǐng)域取得一定的成果［4－5］.但是受限于當(dāng)時控制理論限制，月面返回的控制方案還有待于進一步完善.當(dāng)?shù)诙卧虑蛱綔y高潮來臨時，部分國外學(xué)者對這方面做了更深入的研究［6－9］.大量的研究表明，對于共軛變量初值求解，一般除了特殊系統(tǒng)外，不可能求出最優(yōu)控制的明確解析表達式，需要借助數(shù)值計算方法進行大量的迭代計算.求解的方法有很多，其中，打靶法由于其各種優(yōu)點應(yīng)用較多.但是，應(yīng)用打靶法求解時需要首先猜測未知狀態(tài)變量的初值，一旦猜測偏差過大，計算過程會陷入局部極值點或發(fā)散，在月球軟著陸最優(yōu)軌道的研究中已經(jīng)有了一些可參考的成果.本文首先對有物理意義的量的初值進行猜測，通過解方程對共軛變量初值進行估計，輔助共軛變量初值的計算，用打靶法經(jīng)迭代計算得到初始共軛變量的真值，然后將狀態(tài)方程進行積分求得最優(yōu)開環(huán)控制率.

1 動力學(xué)模型

探測器在進行月面返回時，直接進入一條100 km高度的圓軌道.這里嘗試在三維空間中綜合考慮月球自傳、科里奧利力的影響，建立較為精確的動力學(xué)模型.

目前已發(fā)表的文獻中探測器的動力學(xué)模型大多都是采用二維模型，即假設(shè)月球探測器在一個固定的鉛錘面內(nèi)運動，沒有考慮側(cè)向運動，而且所采用的模型都是在忽略月球自轉(zhuǎn)的基礎(chǔ)上得到的［1－2］.由于在狀態(tài)方程中必須利用位移矢量分量、速度矢量分量共6個狀態(tài)變量實時表示從軌道坐標系到慣性坐標系轉(zhuǎn)換所使用的歐拉角，如果采用普通的三次歐拉角旋轉(zhuǎn)，歐拉角表示將會十分繁瑣.文獻［3］通過對坐標系巧妙的定義，使歐拉角化簡到兩個，實現(xiàn)了實時計算，但這樣就把其探測器軟著陸的初始軌道和終端位置都限定到白道面上，有一定局限性.

本文采用以初始軌道平面確定月心慣性坐標系的思想，即通過初始軌道平面在月球慣性空間中定義慣性坐標系的坐標軸，軌道坐標系的豎軸和慣性坐標系的兩個坐標軸共面，使歐拉角轉(zhuǎn)化為兩個，即方便了計算，也使動力學(xué)模型沒有局限性.

Ox1y1z1為軌道坐標系，原點在航天器，Oy1指向從月心到航天器的延伸線方向，Oz1延垂直于軌道平面方向，Ox1按右手坐標系規(guī)則定義.由于燃料的最優(yōu)要求初始發(fā)射點在目標軌道平面內(nèi)，即整個過程中，Oz1軸方向在月球慣性空間內(nèi)保持不變.所以定義Oxyz為月心慣性坐標系，原點在月心，Oxy平面為月球赤道面，Oz軸指向月球北極，Oy軸延赤道面和Ozz1平面交線，Ox軸按右手系確定.Ox2y2z2為月球固連坐標系，原點在月心，Ox2y2平面為月球赤道面，Ox2軸指向初始時刻Ox軸與月面的交點，Oz2軸指向月球北極與Oz軸重合，Oy2軸按右手系確定.這樣就使Oz軸、Oz1軸、Oy軸3個坐標軸處于同一平面上，如圖1所示.

圖1 坐標系示意圖

發(fā)動機推力F的方向與探測器縱軸重合，θ是推力方向和Oy1的夾角.α是Ox1和Ox夾角，β是Oz1和Oz的夾角.γ為月固坐標系相對于慣性坐標系的轉(zhuǎn)角.

因此，軌道坐標系到慣性坐標系的轉(zhuǎn)換矩陣可表示為

慣性坐標系到月固坐標系的轉(zhuǎn)換矩陣可表示為

根據(jù)牛頓第二定律，得到探測器在慣性坐標系中的運動方程為

其中Vx，Vy，Vz為慣性坐標系中探測器的速度矢量分量;P為探測器發(fā)動機的秒耗量;E為它的比沖;m為探測器質(zhì)量，m=mo－Pt;g為月球引力加速度.又由科里奧利加速度定理有

因此，航天器上升段在月球固連坐標系中的運動方程為

其中ω為月球自轉(zhuǎn)速度;Vxl，Vyl，Vzl是航天器在月球固連坐標系中的速度.

2017年9月—2018年5月，選取哈爾濱醫(yī)科大學(xué)第二臨床醫(yī)學(xué)院參加住院醫(yī)師規(guī)范化培訓(xùn)的醫(yī)學(xué)專碩研究生743名（2015級282人、2016級225人和2017級236人）。同時調(diào)查導(dǎo)師、輪轉(zhuǎn)科室護士長和帶教教師共210人。

2 燃料最優(yōu)控制率設(shè)計

當(dāng)達到預(yù)定軌道時，所剩質(zhì)量最大，也就是消耗質(zhì)量最小，即燃料最優(yōu)，表達式如下:

在動力上升段，當(dāng)發(fā)動機推力恒定時，發(fā)動機的燃料消耗率就已經(jīng)確定，因此性能指標就變?yōu)轱w行時間最小，飛行時間越小，燃料消耗也就越少，即

取系統(tǒng)狀態(tài)變量為

其中xl，yl，zl為月球固連坐標系中的坐標.系統(tǒng)狀態(tài)方程表示如下:

其中:θ是控制變量，μ是月球萬有引力常數(shù);aij(i，j=1，2，3)定義如下:

取共軛變量為

構(gòu)造漢密爾頓函數(shù)如下:

燃料最優(yōu)就是要找到一組容許的控制，調(diào)整推力的方向，使探測器飛行時間最小.根據(jù)極大值原理，設(shè)控制量取值范圍不受限，得到極值條件

因此，所得到的最優(yōu)控制率為

共軛方程為

3 終端約束

月面返回動力上升段初始條件靜止，各個狀態(tài)變量都為零.終端條件是要使之進入指定環(huán)月軌道，終端時刻tf自由.需要滿足一定的終端條件:1)速度向量達到一定約束條件.2)位移向量達到一定約束條件.3)速度向量與位移向量方向垂直.4)速度向量和位移向量都在預(yù)定軌道平面內(nèi).因此得到的終端約束集為

橫截條件為

其中ξ是拉格朗日乘子，將最優(yōu)控制率帶入狀態(tài)方程和共軛方程，利用初始條件進行積分，即可得到月面返回最優(yōu)軌跡，此時最優(yōu)軌道的求解就轉(zhuǎn)化成兩點邊值問題的求解.注意:終端約束集中，速度、位移各個分量都是指在慣性坐標系中，進行計算時要通過慣性坐標系和月球固連坐標系轉(zhuǎn)換矩陣轉(zhuǎn)換.

4 最優(yōu)控制問題求解

選取動力返回時終端指標函數(shù)

其中kn(n=1，2，…，5)為加權(quán)系數(shù).根據(jù)初始時刻推力方向垂直向下，即θ=0，可得到λ1λ2λ3初始值對應(yīng)比例關(guān)系，大大方便了計算.利用matlab編程，以共軛變量為參數(shù)，通過優(yōu)化λ極小值終端目標，最終得到一組滿足終端約束條件的共軛變量初值.用求得的共軛變量初值和系統(tǒng)狀態(tài)變量初值進行積分，即可得到一條月面返回最優(yōu)軌跡.

5 算例及仿真結(jié)果

設(shè)探測器初始質(zhì)量為6 000 kg，定常推力發(fā)動機推力15 000 N，比沖3 000 m/s，初始速度為0，航天器初始位置北緯6°，東經(jīng)26°.并設(shè)目標軌道是一條100 km高的圓形環(huán)月停泊軌道.另外，初始時刻，γ=0，α，β通過狀態(tài)變量的值進行實時計算.可求得共軛變量初值如下:

經(jīng)過仿真計算，求得推力方向角度、終端約束指標、質(zhì)量隨時間的變化曲線，如圖2所示.

由于月球是逆時針旋轉(zhuǎn)的，可以看到推力方向角由零開始逐漸負向減小.也就是說，探測器沿月球自轉(zhuǎn)方向順行發(fā)射.終端約束指標在850 s左右時達到極小值，幾乎為零，也就是進入了預(yù)定軌道.最終消耗掉探測器一多半質(zhì)量燃料，得到最優(yōu)軌跡如圖3所示.

圖2 最優(yōu)控制率及終端約束指標

圖3 月面返回動力上升段最優(yōu)軌跡

最終得到的最優(yōu)軌道終端參數(shù)分別為

6 結(jié)論

在月球探測器月面返回問題研究中，綜合考慮月球自轉(zhuǎn)，科里奧利力的影響，建立探測器在三維空間中的較精確動力學(xué)模型.利用Pontryagin極大值原理，基于燃耗最優(yōu)的原則，設(shè)計月面返回最優(yōu)控制律.在數(shù)值計算中，采用打靶法，以共軛變量初值為參數(shù)，以月面返回動力上升段終端狀態(tài)為優(yōu)化指標，綜合考慮位置和速度約束，通過參數(shù)優(yōu)化，得到了滿足終端約束的一組共軛變量初值，同時得到了探測器月面返回的最優(yōu)軌跡，可以較高精度的實現(xiàn)探測器月面返回，對于我國探月三期工程具有實際工程應(yīng)用參考價值.

［1］王大軼，李鐵壽，馬興瑞.月球最優(yōu)軟著陸兩點邊值問題的數(shù)值解法［J］.航天控制，2000，(3):44－49.

［2］徐敏，李俊峰.月球探測器軟著陸的最優(yōu)控制［J］.清華大學(xué)學(xué)報，2001，41(8):81－89.

［3］周凈揚，周荻.月球探測器軟著陸精確建模及最優(yōu)軌道設(shè)計［J］.宇航學(xué)報，2007，28(6):1462－1471.

［4］ARMSTRONG E S，SUDDATH J H.Application of pontryagin maximum principle to the lunar orbit rendezvous problem［R］.Langley Station，Hampton，Virginia，USA:NASA Langley Research Center，1963.

［5］ARMSTRONG E S，CHILDS A G，MARKOS A T.A method for computing extremal maximum-range thrustlimited rocket trajectories with application to lunar transport［R］.［S.l.］:NASA，1970.

［6］MIELE A，WANG T，MANCUSO S.Optimal free-return trajectories for moon missions and mars missions［J］.The Journal of the Astronautical Sciences，2000，48 (2):183－206.

［7］THORNE J D，HALL C D.Minimum-time continuousthrust orbit transfer［J］.The Journal of the Astronautical Sciences，1997，45(4):411－432.

［8］WILSON J W，SIMONSEN L C，SHINN J L，et al.Radiation analysis for the human lunar return mission［R］.Langley Station，Hampton，Virginia，USA: NASA Langley Research Center，1997.

［9］RYAN P，RUSSELL，CESAR A.Geometric analysis of free-Return trajectories following a gravity-assisted flyby［J］.Journal of Spacecraft and Rocket，2005，42(1): 138－151.