霍明英，彭福軍，趙 鈞，謝少彪，3，齊乃明

(1．哈爾濱工業(yè)大學(xué)航天學(xué)院，哈爾濱150006;2．上海宇航系統(tǒng)工程研究所，上海201108;3．上海航天技術(shù)研究院，上海201109)

0 引言

谷神星(Ceres)是太陽(yáng)系內(nèi)的一顆矮行星，同時(shí)也是位于小行星帶內(nèi)最大的天體。目前，谷神星被很多學(xué)者認(rèn)為是尚存的原行星(萌芽期的行星)之一，于45．7億年前在小行星帶中形成［1］。雖然太陽(yáng)系內(nèi)大多數(shù)的原行星不是和其他的原行星合并成為類(lèi)地行星，就是被木星彈射到太陽(yáng)系外，但是谷神星被認(rèn)為是留存下來(lái)較為完整的一顆原行星［2］。因此，對(duì)谷神星的探測(cè)將非常有助于了解太陽(yáng)系的起源。不僅如此，歐洲航天局已于2014年1月確認(rèn)谷神星上有水蒸氣冒出。谷神星的紅外線光譜也顯示水合礦物非常廣泛的存在于谷神星上，這都證明在其內(nèi)部可能存在著大量的水。甚至有天體生物學(xué)方面的科學(xué)家認(rèn)為谷神星有可能是潛在的宜居性星球，因?yàn)槠浔砻娲嬖谥枰娜齻€(gè)基本條件:液態(tài)水、能量來(lái)源和某些化學(xué)成分［3］。因此，對(duì)谷神星的探測(cè)不僅能夠?qū)μ?yáng)系形成理論提供佐證，也可以對(duì)其地表是否存在水甚至是原始生命一探究竟。鑒于谷神星探測(cè)具有很大的科學(xué)意義和現(xiàn)實(shí)意義，美國(guó)航空航天局(NASA)于2007年9月發(fā)射了“黎明”號(hào)太空飛船，于2011年7月抵達(dá)灶神星(Vesta)，2015年3月抵達(dá)谷神星［4］。為了將來(lái)能更快更廉價(jià)地對(duì)谷神星進(jìn)行探測(cè)甚至取樣返回，本文提出將電動(dòng)太陽(yáng)風(fēng)帆(簡(jiǎn)稱(chēng)電動(dòng)帆)推進(jìn)技術(shù)應(yīng)用于谷神星的探測(cè)任務(wù)中。

電動(dòng)帆是由芬蘭航天專(zhuān)家Janhunen于2004年在磁帆［5］的基礎(chǔ)上提出的一種新興的無(wú)推進(jìn)劑消耗的推進(jìn)方式［6］。與更為人熟知的太陽(yáng)帆相同，電動(dòng)帆可以在不消耗任何推進(jìn)劑的情況下產(chǎn)生連續(xù)的推力，因此非常適用于長(zhǎng)期的空間探測(cè)任務(wù);與太陽(yáng)帆不同的是，電動(dòng)帆的動(dòng)力來(lái)源為太陽(yáng)風(fēng)中帶電粒子的動(dòng)能沖力，而不是太陽(yáng)光壓。電動(dòng)帆原理示意圖如圖1所示，電動(dòng)帆由許多根長(zhǎng)而細(xì)的金屬鏈所組成，這些金屬鏈通過(guò)航天器自旋進(jìn)行初始展開(kāi)和形狀保持。航天器上的太陽(yáng)能電子槍通過(guò)向外噴射電子使金屬鏈?zhǔn)冀K保持在高度的正電位。這些帶電金屬鏈會(huì)排斥太陽(yáng)風(fēng)中的正電粒子，從而最終利用太陽(yáng)風(fēng)的動(dòng)能沖力為航天器產(chǎn)生連續(xù)的推力，并使航天器在飛行過(guò)程中不消耗任何推進(jìn)劑。Janhunen對(duì)電動(dòng)帆的推進(jìn)效率與目前常規(guī)推進(jìn)技術(shù)進(jìn)行了對(duì)比分析。分析結(jié)果表明由于電動(dòng)帆在工作過(guò)程中不損耗任何推進(jìn)劑，因此能夠?yàn)轱w行器提供持續(xù)不斷的推力，工作壽命很長(zhǎng)。通過(guò)理論計(jì)算表明，在為期10年的任務(wù)中，同樣質(zhì)量的電動(dòng)帆所能產(chǎn)生的速度增量是同樣質(zhì)量化學(xué)火箭發(fā)動(dòng)機(jī)所能產(chǎn)生增量的1000倍左右，同樣質(zhì)量電火箭發(fā)動(dòng)機(jī)(SEP)所能產(chǎn)生增量的100倍左右［7］。同太陽(yáng)帆對(duì)比方面，由于太陽(yáng)帆產(chǎn)生的推力與相對(duì)太陽(yáng)距離平方成反比關(guān)系，而電動(dòng)帆產(chǎn)生的推力與相對(duì)太陽(yáng)距離一次方成反比［8］。因此，在星際遠(yuǎn)航任務(wù)中，電動(dòng)帆的推力衰減速度相對(duì)于太陽(yáng)帆的推力衰減速度較慢。而且由于電動(dòng)帆利用正電場(chǎng)對(duì)太陽(yáng)風(fēng)質(zhì)子流進(jìn)行反射，減少了反射材料質(zhì)量。在同等特征加速度要求下，電動(dòng)帆比太陽(yáng)帆質(zhì)量更輕;在同等載荷比例的情況下，電動(dòng)帆產(chǎn)生的加速度更大。

圖1 電動(dòng)帆原理示意圖Fig．1 Conceptual sketch of an electric sail

鑒于電動(dòng)帆在深空探測(cè)中的巨大潛力，歐盟于2010年12月組織召開(kāi)了電動(dòng)帆項(xiàng)目(項(xiàng)目編號(hào):ESAIL EU FP7)啟動(dòng)會(huì)議，資助“電動(dòng)帆推進(jìn)技術(shù)研究項(xiàng)目”170萬(wàn)歐元資金。在“歐洲第七框架計(jì)劃”的資助下，歐洲各國(guó)相關(guān)學(xué)者對(duì)電動(dòng)帆開(kāi)展了一系列理論及實(shí)驗(yàn)研究。Mengali和Quarta基于火星探測(cè)任務(wù)為背景采用間接優(yōu)化算法對(duì)電動(dòng)帆和太陽(yáng)帆的運(yùn)動(dòng)軌跡進(jìn)行了優(yōu)化，并對(duì)比了電動(dòng)帆與太陽(yáng)帆的性能。通過(guò)對(duì)比可以看出，太陽(yáng)帆產(chǎn)生的推力與相對(duì)太陽(yáng)距離平方成反比關(guān)系，而電動(dòng)帆產(chǎn)生的推力與相對(duì)太陽(yáng)距離7/6次方成反比(Janhunen已通過(guò)實(shí)驗(yàn)更正為電動(dòng)帆產(chǎn)生的推力與相對(duì)太陽(yáng)距離一次方成反比［8］)。因此，在星際遠(yuǎn)航任務(wù)中，電動(dòng)帆的推力衰減速度相對(duì)于太陽(yáng)帆的推力衰減速度較慢［9］。他們還將電動(dòng)帆應(yīng)用于太陽(yáng)系外探測(cè)任務(wù)，并通過(guò)間接優(yōu)化算法對(duì)不同特征加速度電動(dòng)帆的太陽(yáng)系外探測(cè)軌跡進(jìn)行了優(yōu)化設(shè)計(jì)。軌跡優(yōu)化結(jié)果表明，一個(gè)中等性能的電動(dòng)帆完成太陽(yáng)系外探測(cè)所用的時(shí)間為15年［10］，是旅行者1號(hào)的1/3左右。他們基于間接優(yōu)化方法還對(duì)電動(dòng)帆的各種空間飛行任務(wù)進(jìn)行了軌跡優(yōu)化，結(jié)果均表明電動(dòng)帆在深空探測(cè)領(lǐng)域方面有很大的潛力［11－13］。

本文作者在前期的研究工作中，采用高斯偽譜法對(duì)電動(dòng)帆自地球至火星的二維軌跡優(yōu)化問(wèn)題進(jìn)行了研究［14］。在研究中忽略了火星軌道面與地球軌道面之間的傾角以及星歷約束，將實(shí)際的三維軌跡優(yōu)化問(wèn)題簡(jiǎn)化為二維軌跡優(yōu)化問(wèn)題。由于谷神星公轉(zhuǎn)軌道的軌道傾角較大(10．585°)，所以這種簡(jiǎn)化假設(shè)將不適用于谷神星探測(cè)任務(wù)分析中。本文在考慮星歷約束的情況下，對(duì)電動(dòng)帆自地球至谷神星的三維軌跡優(yōu)化問(wèn)題開(kāi)展研究。并提出一種基于高斯偽譜法和遺傳算法的混合優(yōu)化算法，這種優(yōu)化算法采用遺傳算法生成高斯偽譜法中所需的狀態(tài)變量及控制變量初值，既避免了單純高斯偽譜法初值猜測(cè)繁瑣的問(wèn)題，也克服了Mengali和Quarta所采用的間接優(yōu)化算法對(duì)協(xié)態(tài)變量初值敏感的問(wèn)題?；诖藘?yōu)化算法，在不同特征推力加速度和發(fā)射時(shí)間的情況下，對(duì)電動(dòng)帆航天器自地球至谷神星的過(guò)渡軌跡進(jìn)行了優(yōu)化，并分析了上述因素對(duì)飛行任務(wù)的影響。

1 轉(zhuǎn)移軌跡優(yōu)化問(wèn)題描述

由于太陽(yáng)風(fēng)粒子流在行星引力球內(nèi)受到行星磁場(chǎng)的影響而變得較為復(fù)雜，所以本文中電動(dòng)帆軌跡優(yōu)化問(wèn)題主要集中于以太陽(yáng)為引力中心的二體問(wèn)題。為了提高數(shù)值計(jì)算效率，引入?yún)⒖季嚯x和參考時(shí)間對(duì)整個(gè)軌跡優(yōu)化問(wèn)題進(jìn)行無(wú)量綱化處理。參考距離為一個(gè)天文單位(1au)，參考時(shí)間為地球相對(duì)太陽(yáng)的公轉(zhuǎn)周期。

圖2 參考坐標(biāo)系及推進(jìn)加速度特征角Fig．2 Reference frames and propulsive acceleration's characteristic angles

1．1優(yōu)化性能指標(biāo)

傳統(tǒng)航天器軌跡優(yōu)化問(wèn)題的性能指標(biāo)大致可以分成三種，即燃料最優(yōu)、時(shí)間最優(yōu)和時(shí)間燃料混合最優(yōu)。由于電動(dòng)帆航天器在飛行過(guò)程中利用太陽(yáng)風(fēng)中帶電粒子的動(dòng)能產(chǎn)生推力而不消耗任何推進(jìn)劑，因此電動(dòng)帆航天器軌跡優(yōu)化問(wèn)題中的優(yōu)化性能指標(biāo)通常選為時(shí)間最優(yōu)［9－10］，優(yōu)化性能指標(biāo)如下:

其中，t0為任務(wù)初始時(shí)刻，tf為任務(wù)終端時(shí)刻。在電動(dòng)帆軌跡優(yōu)化問(wèn)題中，終端時(shí)刻tf通常為一個(gè)未知變量，而初始時(shí)刻t0可以是一個(gè)已知變量，也可以當(dāng)做一個(gè)未知變量來(lái)處理。當(dāng)初始時(shí)刻t0為未知變量時(shí)，軌跡優(yōu)化問(wèn)題不僅要得出最優(yōu)的過(guò)渡軌跡，還要得出使性能指標(biāo)最優(yōu)的任務(wù)開(kāi)始時(shí)間。

1．2動(dòng)力學(xué)微分約束

在描述電動(dòng)帆動(dòng)力學(xué)方程時(shí)主要涉及到兩個(gè)參考系，即日心黃道參考系T⊙(x，y，z)和軌道參考系TO(xo，yo，zo)。如圖2所示，日心黃道參考系的原點(diǎn)為太陽(yáng)中心，正x軸指向歷元J2000．0時(shí)刻平春分點(diǎn)方向，正z軸垂直于J2000．0時(shí)刻黃道面并指向黃道北極方向，y軸與x軸和z軸構(gòu)成右手系。軌道參考系的原點(diǎn)位于電動(dòng)帆航天器質(zhì)心，正zo軸為太陽(yáng)－電動(dòng)帆航天器的矢量方向，yo軸與zo軸和日心黃道參考系中的z軸垂直，方向指向飛行運(yùn)動(dòng)方向，xo軸與yo軸和zo軸構(gòu)成右手系。電動(dòng)帆航天器在日心黃道參考系下的動(dòng)力學(xué)方程為:

式中:μ⊙為太陽(yáng)的引力常數(shù);r=［rx，ry，rz］T為航天器在日心黃道參考系下描述的位置矢量;r=r為航天器與太陽(yáng)的相對(duì)距離;v=［vx，vy，vz］T為航天器在日心黃道參考系下描述的速度矢量;a為電動(dòng)帆的推進(jìn)加速度矢量。

電動(dòng)帆的推力是由太陽(yáng)風(fēng)粒子的速度、密度以及帶電金屬鏈的電壓、長(zhǎng)度和根數(shù)所決定。當(dāng)電動(dòng)帆遠(yuǎn)離(接近)太陽(yáng)時(shí)，所能利用的太陽(yáng)風(fēng)粒子速度和密度將減小(增大)，Janhunen通過(guò)理論推導(dǎo)及實(shí)驗(yàn)發(fā)現(xiàn)電動(dòng)帆產(chǎn)生的推力與相對(duì)太陽(yáng)距離一次方成反比［8］。另外，當(dāng)電動(dòng)帆金屬鏈長(zhǎng)度和根數(shù)一定時(shí)，電動(dòng)帆可以通過(guò)太陽(yáng)能電子槍向外噴射電子調(diào)整金屬鏈的電壓，進(jìn)而實(shí)現(xiàn)對(duì)推力的調(diào)整。參考Janhunen通過(guò)實(shí)驗(yàn)獲得的電動(dòng)帆推力模型［8，15］以及圖2，電動(dòng)帆推力加速度矢量在日心黃道參考系下可以寫(xiě)作:

式中:a⊕為電動(dòng)帆的特征加速度，即電動(dòng)帆距離太陽(yáng)r⊕=1au處所能產(chǎn)生的最大加速度值;κ∈［0，1］為電動(dòng)帆推力開(kāi)關(guān)系數(shù)，可以通過(guò)電子槍調(diào)整金屬鏈的電壓來(lái)調(diào)整電動(dòng)帆整體的推力，這一特性與太陽(yáng)帆有很大的不同;當(dāng)κ=0時(shí)，電子槍處于關(guān)閉狀態(tài)，電動(dòng)帆無(wú)推力輸出;當(dāng)κ=1時(shí)，電動(dòng)帆以最大特征加速度進(jìn)行工作;α∈［0，αmax］為推進(jìn)錐角，即電動(dòng)帆推進(jìn)加速度矢量與zo軸之間的夾角，Janhunen通過(guò)實(shí)驗(yàn)表明，電動(dòng)帆航天器的推進(jìn)錐角不能超過(guò)一個(gè)安全穩(wěn)定閾值αmax，否則會(huì)造成電動(dòng)帆航天器構(gòu)型的不穩(wěn)定;β∈［－π，π］為推進(jìn)鐘角，即電動(dòng)帆推進(jìn)加速度矢量在xoyo平面投影分量與xo軸之間的夾角，逆時(shí)針為正;θ∈［0，π］和φ∈［－π，π］為電動(dòng)帆的極角和方位角，與電動(dòng)帆位置矢量在日心黃道系下坐標(biāo)［r］T⊙?［rx，ry，rz］T的關(guān)系為

將式(3)和式(4)代入式(1)和式(2)，便可獲得電動(dòng)帆航天器的動(dòng)力學(xué)方程:

式中:狀態(tài)變量為電動(dòng)帆的位置及速度，即x?［rx，ry，rz，vx，vy，vz］T;控制變量為影響電動(dòng)帆推進(jìn)加速度矢量的兩個(gè)推進(jìn)角α，β和推力開(kāi)關(guān)系數(shù)κ，即u?［α，β，κ］T;狀態(tài)微分函數(shù)f(x(t)，u(t)，t)簡(jiǎn)記為［f1，f2，f3，f4，f5，f6］T。

1．3邊界約束

電動(dòng)帆航天器軌跡優(yōu)化問(wèn)題中的邊界條件約束主要包括初始狀態(tài)約束和終端狀態(tài)約束。假設(shè)電動(dòng)帆航天器在初始時(shí)刻位于地球逃逸拋物線軌跡上，且逃逸剩余能量為零，初始狀態(tài)約束可寫(xiě)作:

其中rx⊕(t0)，ry⊕(t0)，rz⊕(t0)和vx⊕(t0)，vy⊕(t0)，vz⊕(t0)為地球t0時(shí)刻在日心黃道參考系下的位置和速度，可通過(guò)美國(guó)噴氣推進(jìn)實(shí)驗(yàn)室(JPL)發(fā)布的DE405星歷計(jì)算得出［16］。

為了完成電動(dòng)帆自地球至谷神星的過(guò)渡，須在終點(diǎn)時(shí)刻tf時(shí)的狀態(tài)變量與谷神星狀態(tài)變量一致，所以終端約束條件可寫(xiě)作:

式 中:rxΔ(tf)，ryΔ(tf)，rzΔ(tf)和vxΔ(tf)，vyΔ(tf)，vzΔ(tf)為谷神星tf時(shí)刻在日心黃道參考系下的位置和速度，可通過(guò)谷神星的軌道根數(shù)計(jì)算得出。

1．4路徑約束

由于航天器距離太陽(yáng)越近，電動(dòng)帆帶電金屬鏈所接觸到太陽(yáng)風(fēng)高能粒子的溫度、速度和密度都會(huì)變得越大，電動(dòng)帆的使用壽命會(huì)越短，所以為了保證電動(dòng)帆中的金屬鏈能夠長(zhǎng)期穩(wěn)定地工作，要求電動(dòng)帆在飛行過(guò)程中相對(duì)太陽(yáng)的距離應(yīng)大于電動(dòng)帆航天器允許的最小距離rmin。若電動(dòng)帆金屬鏈材質(zhì)為鋁合金，一般假設(shè)rmin=0．5au。綜上所述，電動(dòng)帆航天器在飛行過(guò)程中的路徑約束可寫(xiě)作:

2 基于高斯偽譜法和遺傳算法的軌跡優(yōu)化策略

針對(duì)單純高斯偽譜法［17］初值賦值繁瑣的問(wèn)題，本文提出一種結(jié)合高斯偽譜法和遺傳算法的混合優(yōu)化方法?；旌蟽?yōu)化算法的優(yōu)化策略如下:首先以較少離散點(diǎn)進(jìn)行高斯偽譜法離散化，并通過(guò)罰函數(shù)將有約束優(yōu)化問(wèn)題轉(zhuǎn)化成無(wú)約束優(yōu)化問(wèn)題，應(yīng)用遺傳算法對(duì)此優(yōu)化問(wèn)題進(jìn)行全局尋優(yōu)，并對(duì)遺傳算法獲得的結(jié)果進(jìn)行三次樣條插值;然后以較多離散點(diǎn)進(jìn)行高斯偽譜法離散化，并以可行解插值結(jié)果作為最優(yōu)解計(jì)算的初值，最后采用序列二次規(guī)劃算法在此初值基礎(chǔ)上計(jì)算得出最優(yōu)解，基于混合優(yōu)化算法的電動(dòng)帆軌跡優(yōu)化流程圖如圖3所示。

圖3 混合優(yōu)化算法計(jì)算流程Fig．3 The flow chart of the hybrid optimization method

2．1基于高斯偽譜法離散化

最優(yōu)控制問(wèn)題的時(shí)間區(qū)間為［t0，tf］，而高斯偽譜法的離散點(diǎn)分布在區(qū)間［－1，1］之間［18］。所以需要通過(guò)引入一個(gè)新的時(shí)間變量τ，將上述定義在［t0，tf］區(qū)間上的最優(yōu)控制問(wèn)題轉(zhuǎn)化成定義在［－1，1］區(qū)間內(nèi)，時(shí)間變換關(guān)系為:

由于電動(dòng)帆航天器軌跡優(yōu)化問(wèn)題中的性能指標(biāo)函數(shù)不包含積分項(xiàng)，所以在離散化后的非線性規(guī)劃問(wèn)題中仍可以寫(xiě)成非常簡(jiǎn)單的形式:

為了將原無(wú)限維的連續(xù)非線性最優(yōu)控制問(wèn)題轉(zhuǎn)化成有限維的非線性規(guī)劃問(wèn)題，需要用N階Lagrange插值多項(xiàng)式對(duì)原問(wèn)題中連續(xù)的狀態(tài)變量和控制變量進(jìn)行插值:

式中:X(τ)和U(τ)分別為狀態(tài)變量和控制變量的N階近似多項(xiàng)式;Lk(τ)，k=1，2，…，N，為N階Lagrange插值多項(xiàng)式:

式中:τk，k=1，2，…，N為N個(gè)Legendre-Gauss(LG)點(diǎn)，狀態(tài)變量和控制變量是在這些非均勻分布的LG點(diǎn)上進(jìn)行離散化的。這N個(gè)LG點(diǎn)是N次Legendre多項(xiàng)式的零點(diǎn)，N次Legendre多項(xiàng)式的表達(dá)式如下:

但是這些離散點(diǎn)未包括初始時(shí)刻及終端時(shí)刻所對(duì)應(yīng)的點(diǎn)，所以定義τ0=－1和τN+1=1分別對(duì)應(yīng)初始時(shí)刻和終端時(shí)刻。為了表達(dá)的方便，簡(jiǎn)記狀態(tài)變量在離散點(diǎn) τi，i=0，．．．，N+1上的值X(τi)為Xi，控制變量在離散點(diǎn) τi，i=0，…，N+1上的值U(τi)為Ui，時(shí)間在離散點(diǎn) τi，i=0，…，N+1上的值t(τi)為ti。

為了將連續(xù)的非線性最優(yōu)控制問(wèn)題轉(zhuǎn)化成非線性規(guī)劃問(wèn)題，需要把狀態(tài)變量的微分約束轉(zhuǎn)化成等式約束。對(duì)式(11)求導(dǎo)可得:

式中:Dk，i，i=0，1，…，N，k=1，2，…，N是狀態(tài)微分矩陣D中的一個(gè)元素，可通過(guò)下式進(jìn)行計(jì)算:

至此，便可將式(5)中描述的6個(gè)狀態(tài)變量微分約束轉(zhuǎn)化成下面所示的6N代數(shù)等式約束:

通過(guò)高斯積分可獲得電動(dòng)帆航天器在終點(diǎn)時(shí)刻tf的狀態(tài)，進(jìn)而將式(7)描述的邊界條件約束轉(zhuǎn)化為6個(gè)等式約束:

式中:wk，k=1，2，…，N，為高斯積分中的積分權(quán)重。

將路徑約束式(8)在LG點(diǎn)上進(jìn)行離散化，便可得到離散化的路徑約束:

綜合式(10)和式(17)～(19)，便可將式(5)～(8)描述的連續(xù)最優(yōu)控制問(wèn)題轉(zhuǎn)化成下述的非線性規(guī)劃問(wèn)題:

對(duì)于初始時(shí)間t0不固定的電動(dòng)帆航天器軌跡優(yōu)化問(wèn)題，通過(guò)高斯偽譜法離散化所得到的非線性規(guī)劃問(wèn)題中共有9N+2個(gè)設(shè)計(jì)變量，分別是6N個(gè)狀態(tài)變量，3N個(gè)控制變量，以及初始時(shí)間t0和終端時(shí)間tf。等式約束個(gè)數(shù)為6N+6個(gè)，不等式約束個(gè)數(shù)為N個(gè)。

2．2基于遺傳算法的初值生成

對(duì)于式(20)所描述的非線性規(guī)劃問(wèn)題，序列二次規(guī)劃算法是最常用的參數(shù)化尋優(yōu)方法，但是在實(shí)際應(yīng)用中卻存在一定困難。當(dāng)選取LG點(diǎn)個(gè)數(shù)N較多時(shí)，設(shè)計(jì)變量數(shù)目會(huì)比較龐大，為序列二次規(guī)劃算法給定設(shè)計(jì)變量初值會(huì)十分繁瑣，且在無(wú)先驗(yàn)知識(shí)情況下有一定難度。假設(shè)LG點(diǎn)個(gè)數(shù)N=60，根據(jù)上一節(jié)討論的內(nèi)容可知設(shè)計(jì)變量個(gè)數(shù)為542。對(duì)于如此多的設(shè)計(jì)變量，在序列二次規(guī)劃算法應(yīng)用中給定設(shè)計(jì)變量初值的工作會(huì)比較繁瑣，且不恰當(dāng)?shù)某踔禃?huì)使問(wèn)題收斂到不可行解。不僅如此，序列二次規(guī)劃算法不具備全局尋優(yōu)的能力，所得到的解主要決定于所猜測(cè)的設(shè)計(jì)變量初值，所以多是靠近初值最近的局部最優(yōu)解。

因此，本文提出采用遺傳算法全局尋優(yōu)獲得高斯偽譜法中非線性規(guī)劃問(wèn)題狀態(tài)變量及控制變量的初值。通常來(lái)說(shuō)遺傳算法只能用于處理無(wú)約束的參數(shù)優(yōu)化問(wèn)題，然而本節(jié)中的非線性規(guī)劃問(wèn)題是有等式約束及不等式約束的。因此，優(yōu)化問(wèn)題的目標(biāo)函數(shù)將不只包括時(shí)間最優(yōu)問(wèn)題，也應(yīng)該通過(guò)增加罰函數(shù)將等式約束考慮在內(nèi)?；谶z傳算法優(yōu)化問(wèn)題的性能指標(biāo)函數(shù)(適應(yīng)度函數(shù))可寫(xiě)成:

式中:M為懲罰系數(shù)。

2．3基于序列二次規(guī)劃的最優(yōu)解計(jì)算

序列二次規(guī)劃算法是目前應(yīng)用最廣泛且適用性最好的一種基于梯度的非線性規(guī)劃問(wèn)題求解方法。序列二次規(guī)劃法的優(yōu)點(diǎn)是計(jì)算效率高、收斂速度快且所得解精度較高，缺點(diǎn)是需要提供初值猜測(cè)，且所得的解對(duì)設(shè)計(jì)變量初值猜測(cè)具有一定依賴(lài)性，非常容易得到局部最優(yōu)解。因此，本文通過(guò)遺傳算法在無(wú)任何初值猜測(cè)的情況下通過(guò)全局搜索生成較少LG離散點(diǎn)的初解，然后通過(guò)插值得到序列二次規(guī)劃算法所需的初值猜測(cè)。這樣就實(shí)現(xiàn)了在無(wú)任何初始輸入的情況下，完成對(duì)電動(dòng)帆航天器軌跡優(yōu)化問(wèn)題最優(yōu)解的求解。避免了在無(wú)先驗(yàn)知識(shí)的情況下，對(duì)初值進(jìn)行猜測(cè)十分繁瑣且困難的情況。而且由于所使用的設(shè)計(jì)變量初值是通過(guò)全局搜索而生成的，所以所獲得的解更接近于全局最優(yōu)解。本文中序列二次規(guī)劃算法所要處理的非線性規(guī)劃問(wèn)題如式(15)所示，這類(lèi)非線性規(guī)劃問(wèn)題可通過(guò)非線性規(guī)劃求解器SNOPT［19－20］進(jìn)行求解。

3 計(jì)算結(jié)果分析

3．1仿真算例

本節(jié)在電動(dòng)帆航天器自地球至谷神星轉(zhuǎn)移軌跡優(yōu)化中，假設(shè)電動(dòng)帆航天器在初始時(shí)刻位于地球逃逸拋物線軌跡上，且逃逸剩余能量C3=0 km2/s2。谷神星的位置和速度通過(guò)軌道根數(shù)(如表1所示)計(jì)算得出，計(jì)算中忽略了谷神星軌道的章動(dòng)影響。電動(dòng)帆航天器的特征加速度為a⊕=1 mm/s2，最大推進(jìn)錐角αmax=35°，最小允許相對(duì)太陽(yáng)距離rmin=0．5au。飛行初始時(shí)刻的選擇范圍為2020年1月1日(JD2458849．5)到2030年12月31日(JD2462866．5)。在基于遺傳算法的初值計(jì)算中，LG離散點(diǎn)個(gè)數(shù)為10，種群大小為200，迭代次數(shù)為50，懲罰系數(shù)為100。在基于序列二次規(guī)劃算法的最優(yōu)解計(jì)算中，LG離散點(diǎn)個(gè)數(shù)為60，約束允許誤差為10－7。

基于上述參數(shù)對(duì)電動(dòng)帆航天器自地球－谷神星飛行軌跡進(jìn)行了優(yōu)化，數(shù)學(xué)仿真運(yùn)行的環(huán)境為MATLAB(R2010a)，運(yùn)行平臺(tái)為雙核2．7GHz主頻的個(gè)人計(jì)算機(jī)，初值生成所用時(shí)間為22分45秒，得出最優(yōu)解所用時(shí)間為14分22秒。由于所采用的軌跡優(yōu)化方法利用遺傳算法進(jìn)行初值生成，其優(yōu)化計(jì)算效率比單純的高斯偽譜法略低，但是對(duì)于離線軌跡優(yōu)化來(lái)說(shuō)仍在可接收的范圍之內(nèi)。電動(dòng)帆航天器自地球－火星飛行軌跡如圖4所示，若電動(dòng)帆航天器于2030年8月3日離開(kāi)地球影響球，將歷時(shí)879天于2032年12月29日抵達(dá)谷神星引力影響范圍。由電動(dòng)帆航天器自地球至谷神星的控制變量－時(shí)間曲線(圖5)可以看出，混合優(yōu)化算法得出的控制變量能夠滿(mǎn)足控制變量約束，且變化相對(duì)平緩，利于工程實(shí)現(xiàn)。

表1 谷神星軌道參數(shù)Table 1 Orbital parameters of Ceres

圖4 地球－谷神星過(guò)渡軌跡(a⊕=1 mm/s2)Fig．4 Earth－Ceres optimal transfer trajectory(a⊕=1 mm/s2)

為了驗(yàn)證軌跡優(yōu)化算法的精度，將所得到的控制變量插值后代入到控制模型中進(jìn)行積分。積分公式為四階－五階龍哥庫(kù)塔積分公式，積分函數(shù)為ode45，積分相對(duì)精度為10－9，積分得出的狀態(tài)變量結(jié)果繪制于圖6和圖7。由這兩個(gè)圖可以看出，混合優(yōu)化算法得到的軌跡與相同控制變量積分得到的軌跡差別較小(在圖6和圖7中兩條軌跡都已畫(huà)出，積分軌跡為實(shí)線，優(yōu)化得出的原軌跡為虛線，但是兩者在圖中過(guò)于接近已經(jīng)基本重合)。這是由于高斯偽譜法所采用的高斯積分公式具有較高的精度，這也正是高斯偽譜法近年來(lái)得到相關(guān)學(xué)者重視的一個(gè)原因。通過(guò)積分結(jié)果對(duì)交會(huì)處的相對(duì)位置和相對(duì)速度進(jìn)行分析，位置終端誤差為35616 km，小于谷神星的引力球半徑(約為47000 km);速度終端誤差為23．526 m/s，相對(duì)谷神星的逃逸能量小于零，能夠?qū)崿F(xiàn)與谷神星的交會(huì)。在真實(shí)的任務(wù)軌跡優(yōu)化中，可通過(guò)增加高斯偽譜法的LG點(diǎn)個(gè)數(shù)來(lái)進(jìn)一步提高交會(huì)精度。

圖5 控制變量－時(shí)間曲線(a⊕=1 mm/s2)Fig．5 Time histories of the control variables(a⊕=1 mm/s2)

圖6 位置矢量－時(shí)間曲線(a⊕=1 mm/s2)Fig．6 Time histories of the position vector(a⊕=1 mm/s2)

圖7 速度矢量－時(shí)間曲線(a⊕=1 mm/s2)Fig．7 Time histories of the velocity vector(a⊕=1 mm/s2)

3．2起始時(shí)間對(duì)探測(cè)任務(wù)的影響

為了進(jìn)一步驗(yàn)證混合優(yōu)化算法的有效性以及考察起始時(shí)間t0對(duì)探測(cè)任務(wù)的影響，本節(jié)對(duì)不同起始時(shí)間情況下電動(dòng)帆航天器自地球至谷神星的飛行軌跡進(jìn)行優(yōu)化。本節(jié)中起始時(shí)間t0不作為設(shè)計(jì)變量，而是自2020年1月1日(JD2458849．5)到2030年12月31日(JD2462866．5)每30天為兩個(gè)仿真的時(shí)間間隔進(jìn)行指定起始時(shí)間的軌跡優(yōu)化，共計(jì)134個(gè)軌跡優(yōu)化算例，其余仿真參數(shù)均與上一節(jié)一致，飛行時(shí)間(tf－t0)總結(jié)在圖8中。

由圖8可以看出，電動(dòng)帆航天器自地球至谷神星的飛行時(shí)間隨著起始時(shí)間的變化呈周期性波動(dòng)，在3720天內(nèi)共波動(dòng)8次，平均變化周期為465天，這一數(shù)值與谷神星和地球的會(huì)合周期466．7天十分接近。出現(xiàn)這種現(xiàn)象的原因是當(dāng)電動(dòng)帆航天器從地球軌道到達(dá)谷神星軌道時(shí)，谷神星也必須同時(shí)達(dá)到，才能使航天器與谷神星實(shí)現(xiàn)交會(huì)。所以航天器發(fā)射時(shí)地球與谷神星的相對(duì)位置將直接影響飛行所需時(shí)間，而兩者相對(duì)位置近似成周期性變化，變化周期既是谷神星和地球的會(huì)合周期。這就是電動(dòng)帆航天器自地球至谷神星飛行時(shí)間變化周期與交會(huì)周期相近的原因。另外，由于谷神星軌道傾角及偏心率等因素的影響，相鄰波動(dòng)段內(nèi)的飛行時(shí)間并不完全相等，而是存在一定差異。

圖8 不同起始時(shí)間情況下的最小飛行時(shí)間(a⊕=1 mm/s2)Fig．8 Minimum transfer time as a function of the starting date(a⊕=1 mm/s2)

3．3特征加速度對(duì)探測(cè)任務(wù)的影響

為了考察電動(dòng)帆航天器特征加速度a⊕對(duì)所需飛行時(shí)間的影響，以不同特征加速度下電動(dòng)帆自地球至谷神星的過(guò)渡軌跡進(jìn)行了優(yōu)化計(jì)算?？紤]的特征加速度范圍為0．4 mm/s2～2 mm/s2，每個(gè)仿真之間的特征加速度間隔為0．1 mm/s2，其余仿真參數(shù)與第3．1節(jié)中一致。不同電動(dòng)帆特征加速度下地球－谷神星過(guò)渡時(shí)間如圖9所示，電動(dòng)帆航天器的特征加速度越小，完成過(guò)渡所需要的飛行時(shí)間越長(zhǎng)。當(dāng)特征加速度小于0．6 mm/s2時(shí)，飛行時(shí)間有顯著的增加。通過(guò)對(duì)軌跡的分析可知，航天器在自地球至谷神星的過(guò)渡過(guò)程中可粗略的劃分成兩個(gè)調(diào)整任務(wù)，即軌道半長(zhǎng)軸的增大和軌道傾角的調(diào)整。電動(dòng)帆特征加速度的大小對(duì)軌道半長(zhǎng)軸調(diào)整所用時(shí)間影響較大，而對(duì)軌道傾角調(diào)整時(shí)間影響相對(duì)較小。當(dāng)電動(dòng)帆加速度大于0．6 mm/s2時(shí)，軌道半長(zhǎng)軸調(diào)整任務(wù)先完成，影響任務(wù)完成時(shí)間的主要限制為軌道傾角的調(diào)整，所以過(guò)渡所用時(shí)間變化相對(duì)平緩。當(dāng)電動(dòng)帆加速度小于0．6 mm/s2時(shí)，軌道傾角調(diào)整任務(wù)先完成，影響任務(wù)完成時(shí)間的主要限制為軌道半長(zhǎng)軸的增長(zhǎng)，所以過(guò)渡所用時(shí)間變化相對(duì)較大。

考慮一個(gè)依據(jù)現(xiàn)有技術(shù)比較合理的電動(dòng)帆特征加速度a⊕=0．6 mm/s2，電動(dòng)帆航天器將歷時(shí)1014天完成自地球至谷神星的過(guò)渡，比采用電推進(jìn)技術(shù)的黎明號(hào)空間探測(cè)器(歷時(shí)1388天完成自地球至灶神星的過(guò)渡，之后又歷時(shí)約920天完成自灶神星至谷神星的過(guò)渡)所用時(shí)間較短。不僅如此，由于電動(dòng)帆在飛行過(guò)程中不消耗任何推進(jìn)劑，所以在完成上述任務(wù)后仍可繼續(xù)進(jìn)行其他飛行任務(wù)，如取樣返回和其他天體探測(cè)，甚至太陽(yáng)系邊界探測(cè)等，減少飛行任務(wù)成本。

圖9 不同特征加速度情況下的最小飛行時(shí)間Fig．9 Minimum transfer time as a function of the characteristic acceleration

通過(guò)以上的軌跡優(yōu)化算例可知，所提出的混合優(yōu)化算法是有效的，能夠在無(wú)任何初值猜測(cè)的情況下完成電動(dòng)帆航天器飛行軌跡的優(yōu)化。另外，一個(gè)具有中等特征加速度的電動(dòng)帆航天器便能在可接受的時(shí)間內(nèi)完成自地球至谷神星的過(guò)渡，所以電動(dòng)帆航天器是適用于谷神星探測(cè)任務(wù)的。

4 結(jié)論

為了將來(lái)能更快更廉價(jià)地對(duì)谷神星進(jìn)行探測(cè)，提出將電動(dòng)帆推進(jìn)技術(shù)應(yīng)用于谷神星探測(cè)任務(wù)中。并提出一種結(jié)合高斯偽譜法和遺傳算法的混合優(yōu)化算法對(duì)電動(dòng)帆自地球至谷神星的軌跡進(jìn)行了優(yōu)化。數(shù)學(xué)仿真結(jié)果表明:1)電動(dòng)帆航天器自地球至谷神星的飛行時(shí)間隨著起始時(shí)間的變化成周期性波動(dòng)，波動(dòng)周期基本與地球－谷神星會(huì)合周期基本一致;2)電動(dòng)帆航天器的特征加速度越小，完成過(guò)渡所需要的飛行時(shí)間越長(zhǎng)，且一個(gè)具有中等特征加速度的電動(dòng)帆航天器便能在可接受的時(shí)間內(nèi)完成自地球至谷神星的過(guò)渡;3)所提出的混合優(yōu)化算法是有效的，能夠在無(wú)任何初值猜測(cè)的情況下完成電動(dòng)帆航天器飛行軌跡的優(yōu)化，避免了單純高斯偽譜法初值賦值繁瑣的問(wèn)題。

［1］Petit J，Morbidelli A．The primordial excitation and clearing of the asteroid belt［J］．Icarus，2001，153(2):338–347．

［2］McCord T B，Sotin C．Ceres:evolution and current state［J］．Journal of Geophysical Research，2005，110(E5):1－14．

［3］ZOL新聞中心．太陽(yáng)系第二大“蓄水池”谷神星或存在外星生命［EB/OL］．2015－01－02［2015－01－02］．http://news．zol．com．cn/article/368722．html．

［4］Rayman M D，Mase R A．Dawn's operations in cruise from Vesta to Ceres［J］．Acta Astronautica，2014，103:113－118．

［5］Zubrin M，Andrews G．Magnetic sails and interplanetary travel［J］．Journal of Spacecraft and Rockets，1991，28(2):197–203．

［6］Janhunen P．Electric sail for spacecraft propulsion［J］．Journal of Propulsion and Power，2004，20(4):763–764．

［7］Janhunen P．Status report of the electric sail in 2009［J］．Acta Astronautica，2011，68(5－6):567－571．

［8］Janhunen P．The electric solar wind sail status report［C］．European Planetary Science Congress，Rome，Italy，Sep 19－24，2010．

［9］Mengali G，Quarta A A，Janhunen P．Electric sail performance analysis［J］．Journal of Spacecraft and Rockets，2008，45(1):122－129．

［10］Quarta A A，Mengali G．Electric sail mission analysis for outer solar system exploration［J］．Journal of Guidance Control and Dynamics，2010，33(3):740－755．

［11］Quarta A A，Mengali G，Janhunen P．Optimal interplanetary rendezvous combining electric sail and high thrust propulsion system［J］．Acta Astronautica，2011，68(5－6):603－621．

［12］Mengali G，Quarta A A，Janhunen P．Considerations of electric sailcraft trajectory design［J］．Journal of British Interplanetary Society，2008，61(8):326－329．

［13］Mengali G，Quarta A A．Escape from elliptic orbit using constant radial thrust［J］．Journal of Guidance Control and Dynamics，2009，32(3):1018–1022．

［14］ 齊乃明，霍明英，袁秋帆．電動(dòng)帆軌跡優(yōu)化及其性能分析［J］．宇航學(xué)報(bào)，2013，34(5):634－641．［Qi Nai-ming，Huo Ming-ying，Yuan Qiu-fan．The electric sail trajectory optimization and performance analysis［J］．Journal of Astronautics，2013，34(5):634－641．］

［15］Janhunen P．Increased electric sail thrust through removal of trapped shielding electrons by orbit chaotisation due to spacecraft body［J］．Annales Geophysica，2009，27(8):3089－3100．

［16］Standish E M．JPL planetary and lunar ephemerides，DE405/LE405［EB/OL］．2008［2008］．http://iau-com m4．jpl．nasa．gov/de405iom/de405iom．pdf．

［17］David B．A Gauss pseudospectral transcription for optimal control［D］．Cambridge:Massachusetts Institute of Technology，2005．

［18］ 唐國(guó)金，羅亞中，雍恩米．航天器軌跡優(yōu)化理論、方法及應(yīng)用［M］．北京:科學(xué)出版社，2011:149－151．

［19］Gill P E，Murray W，Saunders M A．User’s guide for SNOPT V7:software for large-scale nonlinear programming［EB/OL］．2008［2008］．http://ccom．ucsd．edu/～peg．

［20］Gill P E，Murray W．SNOPT:an SQP algorithm for large-scale constrained optimization［J］．SIAM Journal on Optimization，2002，12(4):979–1006．