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

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

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

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

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

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

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

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

圖1 坐標(biāo)系示意圖

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

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

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

根據(jù)牛頓第二定律，得到探測(cè)器在慣性坐標(biāo)系中的運(yùn)動(dòng)方程為

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

因此，航天器上升段在月球固連坐標(biāo)系中的運(yùn)動(dòng)方程為

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

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

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

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

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

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

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

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

取共軛變量為

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

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

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

共軛方程為

3 終端約束

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

橫截條件為

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

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

選取動(dòng)力返回時(shí)終端指標(biāo)函數(shù)

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

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

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

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

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

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

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

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

6 結(jié)論

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

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

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

［3］周凈揚(yáng)，周荻.月球探測(cè)器軟著陸精確建模及最優(yōu)軌道設(shè)計(jì)［J］.宇航學(xué)報(bào)，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.