劉延杰，朱圣英，崔平遠

(1. 北京理工大學宇航學院，北京 100081；2.深空自主導航與控制工信部重點實驗室，北京 100081；3.飛行器動力學與控制教育部重點實驗室，北京 100081)

0 引 言

小天體蘊含著豐富的科學信息，相關(guān)研究的開展有助于人類認識生命的起源與進化，抵御外來天體的撞擊威脅[1]。近年來，隨著人類科學技術(shù)的發(fā)展和經(jīng)濟實力的提升，小天體探測方式逐漸向附著探測和采樣返回探測發(fā)展。附著探測的開展是人類從目標小天體獲得科學信息的有效途徑，也是未來人類進行小天體資源利用的前提條件。

目前，人類已經(jīng)實施了多次小天體附著任務(wù)，包括美國宇航局(NASA)的近地小行星交會任務(wù)(Near Earth asteroid rendezvous，NEAR)，日本宇航局(JAXA)的隼鳥號(Hayabusa)任務(wù)和歐空局(ESA)的羅塞塔(Rosetta)任務(wù)。此外，日本的“隼鳥2號”(Hayabusa 2)任務(wù)在2014年12月發(fā)射，預計對C類小行星1999 JU3進行為期一年半的探測，并完成采樣返回。美國的OSIRIS-REx探測器在2016年9月發(fā)射，對C類小行星101955 Bennu進行采樣返回探測，預計2023年返回地球[2]。

小天體附著軌跡的設(shè)計對探測器能否安全、準確到達預設(shè)的具有科學探測價值的目標區(qū)域起著決定性的作用，所設(shè)計的軌跡需要能夠安全、準確地到達指定著陸點，滿足初末狀態(tài)約束、路徑約束、控制約束等多重約束條件，同時使某項重要的性能指標最優(yōu)化，如燃耗、時間等。

小天體質(zhì)量體積比較小，形狀不規(guī)則，這導致其附近呈現(xiàn)不規(guī)則弱引力場。引力場建模是研究和分析小天體附近運動行為的基礎(chǔ)。經(jīng)典的球諧系數(shù)模型在鄰近小天體的空間范圍內(nèi)存在著不收斂的問題[3]。應(yīng)用多面體模型可以建立小天體附近全局有效的引力場模型，但是該方法計算效率較低[4]。通過改變球諧系數(shù)模型引力勢函數(shù)的提取因子，使勒讓德多項式的收斂條件得到改變，可以建立內(nèi)球諧引力場模型，該模型可確保引力勢函數(shù)表達式在與中心小天體相切的球體范圍內(nèi)收斂[5]。

軌跡優(yōu)化方法主要包括基于龐特里亞金極小值原理的間接法和基于參數(shù)化方法的直接法。在間接法應(yīng)用方面，文獻[6]和文獻[7]首先求解光滑程度較高的能量最優(yōu)問題，通過引入同倫參數(shù)，將能量最優(yōu)問題逐步轉(zhuǎn)化為燃耗最優(yōu)問題，最終得到燃耗最優(yōu)軌跡。同倫方法并不能從根本上解決間接法的協(xié)態(tài)初值敏感問題，求解軌跡優(yōu)化問題仍然存在一定困難。

直接法無需推導優(yōu)化問題的橫截條件，避免了求解兩點邊值問題的協(xié)態(tài)初值敏感問題，因而得到了廣泛的應(yīng)用。文獻[8]和文獻[9]在軌跡優(yōu)化中考慮了動力學參數(shù)和狀態(tài)不確定性的影響，分別將閉環(huán)敏感度矩陣方程和誤差傳播協(xié)方差矩陣引入優(yōu)化指標，進而采用高斯偽譜法對軌跡優(yōu)化問題進行解算。在小天體附著軌跡優(yōu)化方面，直接法仍然以高斯偽譜法為主，解算過程耗時較長，并且解的精度依賴于對狀態(tài)變量和控制變量的初始猜測。

凸優(yōu)化算法是直接法的一種，具有形式簡單，計算效率高的特點[10]。它的子類——二階錐規(guī)劃問題(SOCP)，在航天器軌跡設(shè)計與優(yōu)化方面得到了廣泛應(yīng)用[11-12]。二階錐規(guī)劃問題要求性能指標和動力學方程為線性，且約束條件滿足二階錐約束。對動力學和約束條件進行離散化以后，采取內(nèi)點法[13]對參數(shù)優(yōu)化問題進行求解，可以確保尋優(yōu)結(jié)果能夠在有限次迭代內(nèi)收斂到全局最優(yōu)解。

本文首先采用內(nèi)球諧模型對目標小天體周圍引力場進行精確建模；通過約束條件松弛、動力學方程線性化、動力學和約束條件離散化，將小天體附著多約束軌跡優(yōu)化問題轉(zhuǎn)化為一個可以迭代求解的二階錐規(guī)劃問題，并通過內(nèi)點法進行解算；最后以216 Kleopatra[14]小行星為目標小天體，對上述過程進行仿真求解，以校驗算法的可行性和有效性。

1 小天體附著動力學模型

1.1 內(nèi)球諧引力場模型

引力場建模是研究和分析小天體附近運動行為的基礎(chǔ)，本文采用內(nèi)球諧引力場模型對小天體周圍引力場進行精確建模，內(nèi)球諧坐標系下引力加速度的表達式為[5]：

(1)

(2)

(3)

(4)

式中：xi,yi,zi分別為內(nèi)球諧坐標系下探測器的位置信息。內(nèi)球諧引力系數(shù)可以通過最小二乘法進行采樣估計，計算式為：

(5)

(6)

g可以通過任意已知的引力場模型計算得到，例如多面體模型[4]、質(zhì)子群模型[15]等。為了確保式(5)為一個滿秩系統(tǒng)，數(shù)據(jù)點個數(shù)Ndata應(yīng)該至少為n2+2n個。內(nèi)球諧系數(shù)確定以后，由式(1)～(3)可以快速求出內(nèi)布里淵球內(nèi)任意點的引力加速度。所構(gòu)造內(nèi)布里淵球應(yīng)該與小天體表面的目標著陸點相切，球心選取應(yīng)保證探測器下降軌跡包含在內(nèi)布里淵球內(nèi)部。

1.2 動力學方程及坐標轉(zhuǎn)換

定義小天體固連坐標系如下：坐標原點o位于小天體的質(zhì)心，z軸沿小天體自旋軸方向，x軸沿小天體最小慣量軸方向，y軸滿足右手系法則。在小天體固連坐標系下，探測器的附著動力學方程可以寫為：

(7)

探測器附著動力學方程是在小天體固連系下推導得到，而內(nèi)球諧引力場模型的引力加速度解算是在內(nèi)球諧坐標系進行的，兩者之間需要進行坐標轉(zhuǎn)換。為了簡化計算，可以令內(nèi)球諧坐標系的三個坐標軸xi,yi,zi分別與小天體固連坐標系的x軸、y軸、z軸平行，如圖1所示。此時，兩個坐標系間的轉(zhuǎn)換關(guān)系為：

ri=r-Δr

(8)

式中：Δr表示內(nèi)球諧坐標系原點與小天體質(zhì)心之間的位置矢量差。

2 小天體附著軌跡優(yōu)化方法

本節(jié)首先建立小天體附著多約束軌跡優(yōu)化方程，進而通過約束條件松弛、動力學方程線性化、動力學和約束條件離散化，將原始軌跡優(yōu)化問題轉(zhuǎn)化為一個可以迭代求解的二階錐規(guī)劃問題，并通過內(nèi)點法進行解算。

2.1 小天體附著軌跡優(yōu)化問題描述

探測器應(yīng)該滿足如下的邊界條件：

(9)

式中：r0和v0為探測器在初始時刻的位置和速度，m0表示探測器初始時刻的質(zhì)量，rt為目標著陸點，vt為終端速度約束，對于軟著陸問題，vt=0。此外，探測器推力應(yīng)該滿足如下控制約束：

(10)

式中：Tmax為發(fā)動機推力最大值。以燃耗作為優(yōu)化指標，可以表述如下：

(11)

式中：tf為附著所需時間。

優(yōu)化指標(11)，動力學方程(7)，邊界約束(9)及控制約束(10)構(gòu)成小天體附著燃耗最優(yōu)軌跡優(yōu)化方程。

2.2 序列凸優(yōu)化算法

(12)

(13)

松弛后的優(yōu)化問題可以寫為：

minJT，Γ=∫tf0Γdτ s．t． r=v=-2ω×v-ω×(ω×r)+g+Tm=-ΓIspger(0)=r0v(0)=v0m(0)=m0r(tf)=rfv(tf)=03×1T2x+T2y+T2z≤Γ0≤?！躎maxì?í??????????????????(14)

二階錐規(guī)劃問題要求系統(tǒng)的動力學方程必須是線性的，所以需要對式(14)中的微分方程組進行線性化處理。定義一組新變量如下：

(15)

將以上變量代入式(14)。此時，動力學方程為：

(16)

邊界約束條件為：

(17)

控制約束條件為：

(18)

式中：p0(t)=ln(m0-Tmaxt/Ispge)，優(yōu)化指標可以寫為：

(19)

通過以上線性化過程，優(yōu)化變量從[T,Γ]轉(zhuǎn)換為[u,σ]，從而消除了原動力學方程中分母變量m帶來的非線性問題。

其中,

(20)

式(20)中各矩陣定義為：

其中，uk,σk, ▽Vk分別表示u,σ, ▽V在時間節(jié)點tk時刻的值。離散化后的軌跡優(yōu)化問題可以概括如下：

minJuk，σk=∑Nk=0σk s．t． X0=[rT0，vT0，p0]TMXN=[rTf，01×3]Tuk≤σk0≤σk≤Tmaxe-p0(tk)[1-(pk-p0(tk))]ì?í??????(21)

式中：M=[I6,06×1]，k從0取到N。經(jīng)過離散化處理以后，微分方程約束被轉(zhuǎn)化為代數(shù)方程約束。原來的軌跡優(yōu)化問題轉(zhuǎn)換為一個離散的參數(shù)優(yōu)化問題。

式(21)中的終端狀態(tài)XN需要由式(20)遞推得到。遞推過程中，需要求解每個時間節(jié)點處的引力加速度信息，通常情況下，引力加速度是關(guān)于位置r的非線性函數(shù)，這使得式(21)不是一個標準的二階錐規(guī)劃問題，不能直接采用內(nèi)點法進行求解。為此，本文采用一種迭代算法來進行優(yōu)化問題解算，即首先假設(shè)引力加速度為一個常量，例如▽V0，此時式(21)成為一個二階錐規(guī)劃問題，采用內(nèi)點法可以得到此時的最優(yōu)軌跡；以得到的軌跡作為參考軌跡，由式(1)～(3)計算該軌跡上每個時間節(jié)點的引力加速度值，再代入式(21)，采用內(nèi)點法進行求解，得到一條新的軌跡，將得到的新軌跡作為下一次迭代的參考軌跡；經(jīng)過多次迭代，直到軌跡收斂，即前后兩次迭代得到的優(yōu)化軌跡之差小于某一個固定值ε，從而得到小天體附著燃耗最優(yōu)軌跡。以上每次迭代過程，都需要求解一次二階錐規(guī)劃問題，這種迭代算法也被稱作序列凸優(yōu)化算法。

3 仿真結(jié)果及分析

本節(jié)以216 Kleopatra小行星為目標小天體進行仿真分析，以校驗算法的可行性和有效性。216 Kleopatra小行星尺寸約為217×94×81 km，光譜類型為M類。216 Kleopatra小行星多面體模型如圖2所示，仿真參數(shù)設(shè)置如表1所示。

時間x/my/mz/mvx/(m·s-1)vy/(m·s-1)vz/(m·s-1)m/kgt=0-1501086010-1034100500t=tf-1126006340-12520000-

根據(jù)探測器的初始位置和目標著陸點，構(gòu)造包含附著軌跡的內(nèi)球諧引力場模型。采用多面體模型作為參考，階次n選為30，所構(gòu)造內(nèi)球諧引力場模型誤差分布如圖3所示。從圖3可以看出，誤差由內(nèi)布里淵球球心沿半徑方向逐漸增大，最高不超過6%，可以滿足軌跡優(yōu)化精度需求。

小行星216 Kleopatra的自旋角速度為3.2411×10-4rad/s,發(fā)動機比沖Isp=300 s,推力幅值Tmax=100 N,仿真時間tf=1000 s,離散時間節(jié)點N=500。采用序列凸優(yōu)化算法，進行軌跡優(yōu)化解算，內(nèi)點法由Matlab環(huán)境下的MOSEK[16]軟件實現(xiàn)，令ε=0.3 m，則經(jīng)過四次迭代以后，軌跡最終收斂，優(yōu)化軌跡如圖4所示。優(yōu)化得到的末端時刻探測器的質(zhì)量為479.9 kg，即消耗燃料20.1 kg。

探測器推力矢量大小隨時間變化曲線如圖5所示?？梢园l(fā)現(xiàn)，推力大小滿足最優(yōu)控制的bang-bang控制形式，推力變化大概可以分為三個階段，在前期，發(fā)動機處于飽和狀態(tài)實現(xiàn)對探測器動力加速；在中間階段，將發(fā)動機關(guān)閉，令探測器自由運動，達到節(jié)省燃料消耗的目的；在最終階段，發(fā)動機全開提供反向推力，實現(xiàn)探測器的軟著陸。在發(fā)動機工作的時段，總推力的大小基本保持在100 N，即發(fā)動機在點火工作時處于最大推力狀態(tài)，只是推力方向不斷變化。

迭代次數(shù)引力加速度解算耗時/s內(nèi)點法解算耗時/s100．5622．840．5932．790．6343．240．59

每次迭代引力加速度估計和優(yōu)化耗時如表2所示。由于內(nèi)球諧引力場模型為冪級數(shù)形式，避免了求解線積分和面積分的過程，使得計算效率較多面體模型有了極大地提升。此外，將原始軌跡優(yōu)化問題轉(zhuǎn)換為二階錐規(guī)劃問題以后，形式簡單，求解快速，盡管增加了迭代過程，但總耗時僅為11.24 s。

為了進一步驗證凸優(yōu)化算法在軌跡優(yōu)化計算效率方面的優(yōu)勢，采用高斯偽譜法進行對比分析，離散節(jié)點個數(shù)N同樣設(shè)為500。高斯偽譜法采用Matlab環(huán)境下的GPOPS-II軟件實現(xiàn)[17]，兩種方法沿x方向的優(yōu)化軌跡如圖6所示。由圖6可知，高斯偽譜法優(yōu)化軌跡與序列凸優(yōu)化軌跡幾乎吻合，優(yōu)化得到的末端時刻探測器的質(zhì)量同為479.9 kg，高斯偽譜法的優(yōu)化用時達到了32830.49 s。仿真結(jié)果表明，本文所采用的序列凸優(yōu)化算法可以在不降低優(yōu)化性能的前提下，極大地提升小天體附著軌跡優(yōu)化問題的求解效率。

4 結(jié)束語

本文針對小天體附著多約束軌跡優(yōu)化問題，提出了一種基于序列凸優(yōu)化的軌跡優(yōu)化算法。首先采用內(nèi)球諧引力場模型，對形狀不規(guī)則的小天體進行引力場精確建模；通過約束松弛、線性化和離散化處理，將小天體附著多約束軌跡優(yōu)化方程轉(zhuǎn)換為一個可以迭代求解的二階錐規(guī)劃問題，并由內(nèi)點法得到燃耗最優(yōu)軌跡。數(shù)學仿真結(jié)果顯示，采用序列凸優(yōu)化算法得到的附著軌跡以零速度到達預定著陸點，符合各項約束條件，并實現(xiàn)燃耗最優(yōu)的目標。

[1] 崔平遠, 喬棟. 小天體附近軌道動力學與控制研究現(xiàn)狀與展望[J]. 力學進展, 2013(5): 526-539. [Cui Ping-yuan, Qiao Dong. State-of-the-art and prospects for orbital dynamics and control near small celestial bodies[J]. Advances in Mechanics, 2013(5): 526-539.]

[2] 崔平遠, 袁旭, 朱圣英, 等. 小天體自主附著技術(shù)研究進展[J]. 宇航學報, 2016, 37(7): 759-767. [Cui Ping-yuan, Yuan Xu, Zhu Sheng-ying, et al. Research progress of small body autonomous landing techniques[J]. Journal of Astronautics, 2016, 37(7): 759-767.]

[3] Scheeres D J, Williams B G, Miller J K. Evaluation of the dynamic environment of an asteroid: applications to 433 Eros[J]. Journal of Guidance, Control, and Dynamics, 2000, 23(3): 466-475.

[4] Werner R A, Scheeres D J. Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of asteroid 4769 Castalia[J]. Celestial Mechanics and Dynamical Astronomy, 1996, 65(3): 313-344.

[5] Takahashi Y, Scheeres D J, Werner R A. Surface gravity fields for asteroids and comets[J]. Journal of Guidance, Control, and

Dynamics, 2013, 36(2): 362-374.

[6] Ren Y, Shan J J. Reliability-based soft landing trajectory optimization near asteroid with uncertain gravitational field[J]. Journal of Guidance, Control, and Dynamics, 2015, 38(9): 1810-1820.

[7] Yang H, Baoyin H. Fuel-optimal control for soft landing on an irregular asteroid[J]. IEEE Transactions on Aerospace and Electronic Systems, 2015, 51(3): 1688-1697.

[8] 胡海靜，高艾，朱圣英，等. 考慮跟蹤制導的小天體著陸軌跡閉環(huán)優(yōu)化方法[J]. 宇航學報, 2015, 36(12): 1384-1390. [Hu Hai-jing, Gao Ai, Zhu Sheng-ying, et al. Closed-loop trajectory optimization for precision landing on small bodies considering tracking guidance[J]. Journal of Astronautics, 2015, 36(12): 1384-1390.]

[9] Hu H, Zhu S, Cui P. Desensitized optimal trajectory for landing on small bodies with reduced landing error[J]. Aerospace Science and Technology, 2016, 48: 178-185.

[10] Liu X, Lu P, Pan B. Survey of convex optimization for aerospace applications[J]. Astrodynamics, 2017, 1(1): 23-40.

[11] Acikmese B, Ploen S R. Convex programming approach to powered descent guidance for Mars landing[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(5): 1353-1366.

[12] Blackmore L, Acikmese B, Scharf D P. Minimum-landing-error powered-descent guidance for Mars landing using convex optimization[J]. Journal of guidance, control, and dynamics, 2010, 33(4): 1161-1171.

[13] Nesterov Y E, Todd M J. Self-scaled barriers and interior-point methods for convex programming[J]. Mathematics of Operations Research, 1997, 22(1): 1-42.

[14] Descamps P, Marchis F, Berthier J, et al. Triplicity and physical characteristics of asteroid (216) Kleopatra[J]. Icarus, 2011, 211(2): 1022-1033.

[15] Park R S, Werner R A, Bhaskaran S. Estimating small-body gravity field from shape model and navigation data[J]. Journal of Guidance, Control, and Dynamics, 2010, 33(1): 212-221.

[16] Andersen E D, Roos C, Terlaky T. On implementing a primal-dual interior-point method for conic quadratic optimization[J]. Mathematical Programming, 2003, 95(2): 249-277.

[17] Patterson M A, Rao A V. GPOPS-II: A MATLAB software for solving multiple-phase optimal control problems using hp-adaptive gaussian quadrature collocation methods and sparse nonlinear programming[J]. ACM Transactions on Mathematical Software, 2014(41): 1.