序列凸優(yōu)化的小天體附著軌跡優(yōu)化

劉延杰，朱圣英，崔平遠(yuǎn)

(1. 北京理工大學(xué)宇航學(xué)院，北京 100081；2.深空自主導(dǎo)航與控制工信部重點(diǎn)實(shí)驗(yàn)室，北京 100081；3.飛行器動(dòng)力學(xué)與控制教育部重點(diǎn)實(shí)驗(yàn)室，北京 100081)

0 引 言

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

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

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

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

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

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

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

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

1 小天體附著動(dòng)力學(xué)模型

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

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

(1)

(2)

(3)

(4)

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

(5)

(6)

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

1.2 動(dòng)力學(xué)方程及坐標(biāo)轉(zhuǎn)換

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

(7)

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

ri=r-Δr

(8)

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

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

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

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

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

(9)

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

(10)

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

(11)

式中：tf為附著所需時(shí)間。

優(yōu)化指標(biāo)(11)，動(dòng)力學(xué)方程(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)的動(dòng)力學(xué)方程必須是線性的，所以需要對(duì)式(14)中的微分方程組進(jìn)行線性化處理。定義一組新變量如下：

(15)

將以上變量代入式(14)。此時(shí)，動(dòng)力學(xué)方程為：

(16)

邊界約束條件為：

(17)

控制約束條件為：

(18)

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

(19)

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

其中,

(20)

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

其中，uk,σk, ▽Vk分別表示u,σ, ▽V在時(shí)間節(jié)點(diǎn)tk時(shí)刻的值。離散化后的軌跡優(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)換為一個(gè)離散的參數(shù)優(yōu)化問題。

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語

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

[1] 崔平遠(yuǎn), 喬棟. 小天體附近軌道動(dòng)力學(xué)與控制研究現(xiàn)狀與展望[J]. 力學(xué)進(jìn)展, 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] 崔平遠(yuǎn), 袁旭, 朱圣英, 等. 小天體自主附著技術(shù)研究進(jìn)展[J]. 宇航學(xué)報(bào), 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] 胡海靜，高艾，朱圣英，等. 考慮跟蹤制導(dǎo)的小天體著陸軌跡閉環(huán)優(yōu)化方法[J]. 宇航學(xué)報(bào), 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.