吳愛國,鞏志浩
哈爾濱工業(yè)大學(xué)(深圳),機(jī)電工程與自動(dòng)化學(xué)院,深圳 518055
燃料最優(yōu)是目前世界各國發(fā)射火星探測器的一項(xiàng)重要指標(biāo)。自Howard于1962年首次提出利用氣動(dòng)力輔助軌道轉(zhuǎn)移以來[1],人們開始研究利用大氣阻力對(duì)火星探測器進(jìn)行減速,使其成為環(huán)火星探測器,這種方式能夠大大減少火星探測器在捕獲過程中的燃料消耗。火星探測中的氣動(dòng)輔助技術(shù)主要有氣動(dòng)捕獲和氣動(dòng)剎車2種,其中氣動(dòng)捕獲技術(shù)是指探測器僅在一次穿越火星大氣后即可在大氣阻力的作用下到達(dá)目標(biāo)軌道,而氣動(dòng)剎車則需要多次穿越大氣。相比于氣動(dòng)剎車來說,氣動(dòng)捕獲可以大大縮短軌道轉(zhuǎn)移的時(shí)間。氣動(dòng)捕獲軌道的優(yōu)化問題是關(guān)系到氣動(dòng)捕獲能否成功以及最終效果的關(guān)鍵問題,其實(shí)際上是探測器在穿越火星大氣時(shí),以燃料最優(yōu)或能量最優(yōu)為性能指標(biāo),滿足一系列狀態(tài)約束和邊界約束的非線性優(yōu)化問題。將火星探測器由捕獲軌道以接近最小的燃料消耗和相應(yīng)的進(jìn)入條件轉(zhuǎn)移到指定的目標(biāo)軌道,根據(jù)軌道力學(xué)和可控集理論可以求得入口條件和出口條件[2]。氣動(dòng)捕獲軌道的優(yōu)化問題一般來說沒有解析解,現(xiàn)有的軌跡優(yōu)化方法主要包括間接法、直接法等[3]。
在利用間接法解決軌跡優(yōu)化問題時(shí),首先需要依據(jù)Pontryagin極大值原理將其轉(zhuǎn)化為一個(gè)Hamiltonian兩點(diǎn)邊值問題,然后進(jìn)行求解[4-5]。采用間接法求解軌跡優(yōu)化問題能夠保證最優(yōu)解滿足一階最優(yōu)性必要條件,但是間接法存在求解過程復(fù)雜、收斂域小、對(duì)初值依賴性高且難以估計(jì)等缺點(diǎn)[6]。
相對(duì)于間接法來說,直接法在軌跡優(yōu)化中的應(yīng)用更為廣泛。其利用離散變換將問題轉(zhuǎn)化為非線性規(guī)劃問題,然后根據(jù)優(yōu)化目標(biāo)對(duì)其進(jìn)行尋優(yōu)。其中配點(diǎn)法是廣泛應(yīng)用的間接法,它同時(shí)離散控制變量和狀態(tài)變量,并通過多項(xiàng)式插值對(duì)控制變量和狀態(tài)變量進(jìn)行擬合,通過優(yōu)化配點(diǎn)處的控制變量搜尋滿足約束的最優(yōu)軌跡[7-8]。正交配點(diǎn)法是最常用的一種配點(diǎn)法,又稱作偽譜法。配點(diǎn)一般選擇正交多項(xiàng)式的根,利用全局插值多項(xiàng)式得到狀態(tài)的表達(dá)式,通過對(duì)插值多項(xiàng)式求導(dǎo)將微分方程轉(zhuǎn)化為代數(shù)約束,通過積分得到終端狀態(tài)[9-11]。根據(jù)所選用的插值多項(xiàng)式,可以將偽譜法分為Legendre偽譜法、Gauss偽譜法、Chebyschev偽譜法、Radau偽譜法[6]。
近年來,啟發(fā)式智能優(yōu)化算法憑借其魯棒性、適應(yīng)性以及全局優(yōu)化性能受到越來越多的關(guān)注,已經(jīng)被應(yīng)用到再入軌跡優(yōu)化和高超聲速飛行器上升軌跡的優(yōu)化等問題中,其無需使用哈密頓函數(shù)、也不用計(jì)算其導(dǎo)數(shù)即可獲得接近最佳的結(jié)果[12-13]。另外,也有學(xué)者將啟發(fā)式智能優(yōu)化算法和配點(diǎn)法相結(jié)合,使用較高階多項(xiàng)式來近似控制量,并將此控制變量在配點(diǎn)處進(jìn)行離散,使用啟發(fā)式智能算法優(yōu)化正交配點(diǎn)的值[14-15]。
基于鴿子獨(dú)特的導(dǎo)航能力,文獻(xiàn)[16]提出了一種新型的群體智能優(yōu)化方法:鴿群算法。文獻(xiàn)[16]利用所提出的鴿群算法給出了空間機(jī)器人路徑規(guī)劃的優(yōu)化方法。文獻(xiàn)[17]引入捕食逃逸機(jī)制對(duì)鴿群算法進(jìn)行改進(jìn),并將其用于無人機(jī)緊密編隊(duì)協(xié)同控制的內(nèi)環(huán)控制器的優(yōu)化設(shè)計(jì)。文獻(xiàn)[18]提出了基于自適應(yīng)種群變異的鴿群算法,并將其應(yīng)用于航天器集群的軌道優(yōu)化,提高了種群演化深度和收斂速度。文獻(xiàn)[19]則將鴿群算法與正交配點(diǎn)法相結(jié)合,用于再入航天器的軌道設(shè)計(jì)。文獻(xiàn)[20] 針對(duì)基本鴿群算法的尋優(yōu)缺陷利用量子行為規(guī)則對(duì)其進(jìn)行了改進(jìn),并利用改進(jìn)后的鴿群算法用于無人機(jī)緊密編隊(duì)控制。
本文考慮火星探測器的氣動(dòng)捕獲軌道的優(yōu)化問題。以火星探測器穿越火星大氣時(shí)的動(dòng)力學(xué)方程為基礎(chǔ),分析并確定了能量最優(yōu)氣動(dòng)捕獲軌道,給出了優(yōu)化的性能指標(biāo)。以傾側(cè)角為控制序列,將氣動(dòng)捕獲軌道優(yōu)化問題轉(zhuǎn)化為多參數(shù)優(yōu)化問題。為了解決該多參數(shù)優(yōu)化問題,本文通過引入?yún)?shù)對(duì)基本鴿群算法進(jìn)行了改進(jìn),并將其用于火星探測器的氣動(dòng)捕獲軌道優(yōu)化。
假設(shè)火星大氣相對(duì)于火星靜止,并且認(rèn)為飛行過程中探測器的質(zhì)量變化可忽略,則探測器在火星大氣內(nèi)飛行的動(dòng)力學(xué)方程為[21]
(1)
式中:探測器的狀態(tài)由6個(gè)變量r、λ、θ、v、γ、ψ(分別為火心距、經(jīng)度、緯度、飛行速度、飛行路徑角及偏航角)確定;gM為火星重力加速度;σ為滾轉(zhuǎn)角;D為阻力加速度;L為升力加速度。D、L、gM的表達(dá)式分別為
式中:Sr為探測器的有效截面積;m為探測器的質(zhì)量;CD為阻力系數(shù);CL為升力系數(shù);μ為火星引力常數(shù);ρ為火星大氣密度,不同高度的火星大氣密度不同,火星大氣密度與高度之間的關(guān)系為
ρ=ρ0exp(-h/hs)
(2)
式中:hs=8 805.7 m為火星大氣標(biāo)高;ρ0=0.0147 4 kg/m3為火星大氣參考密度。又由式(1) 可知,火心距、飛行速度和飛行路徑角不改變軌道平面,并且其變化與其他3個(gè)狀態(tài)無關(guān),所以當(dāng)不考慮軌道平面的變化時(shí),上述大氣內(nèi)飛行動(dòng)力學(xué)方程可以簡化為
(3)
利用火星大氣對(duì)探測器進(jìn)行減速,通過調(diào)整姿態(tài)進(jìn)而改變所受升力的大小來控制飛行軌跡,可以認(rèn)為此過程是不消耗能量的無動(dòng)力飛行過程,這里的能量最優(yōu)指的是由捕獲軌道轉(zhuǎn)移到目標(biāo)軌道時(shí)所消耗的能量最少。由捕獲軌道轉(zhuǎn)移到目標(biāo)軌道所需要的速度增量為[2]
(4)
(5)
總速度增量為
ΔV=ΔV1+ΔV2
(6)
則
(7)
又
(8)
(9)
式中:rca為捕獲軌道遠(yuǎn)火點(diǎn)半徑;rcp為捕獲軌道近火點(diǎn)半徑;rta為目標(biāo)軌道遠(yuǎn)火點(diǎn)半徑;rtp為目標(biāo)軌道近火點(diǎn)半徑。由式(9)可知,對(duì)于給定的rcp,ΔV在區(qū)間rca∈[0,rta]內(nèi)單調(diào)遞減,在區(qū)間rca∈[rta,+∞]內(nèi)單調(diào)遞增。因此,當(dāng)rca=rta時(shí),ΔV取得最小值,即當(dāng)探測器穿越大氣后的后捕獲軌道的遠(yuǎn)地點(diǎn)高度與目標(biāo)軌道高度相等,也就是二者相切時(shí),軌道轉(zhuǎn)移所需的能量消耗最少。
由動(dòng)量矩守恒原理和能量守恒原理可知
rtavca=ravfcosγf
(10)
式中:vf為探測器飛出大氣時(shí)的速度;γf為探測器飛出大氣時(shí)的航跡角;ra為火星大氣的高度;vca為捕獲軌道遠(yuǎn)火點(diǎn)處速度。根據(jù)能量守恒有
(11)
聯(lián)立式(10)和式(11)可得
(12)
由式(12)可得
(13)
因此,當(dāng)捕獲軌道遠(yuǎn)火點(diǎn)與目標(biāo)軌道遠(yuǎn)火點(diǎn)相切時(shí)有
(14)
式中:vta為目標(biāo)軌道遠(yuǎn)火點(diǎn)處的速度,當(dāng)目標(biāo)軌道確定時(shí)其為一常值。由式(14)可知,探測器到達(dá)目標(biāo)軌道高度時(shí)所需脈沖變軌的大小與飛出角度γf密切相關(guān)??梢?,飛出角度γf越小,在到達(dá)目標(biāo)軌道高度時(shí)所需的脈沖ΔV也就越小。
氣動(dòng)捕獲軌道優(yōu)化問題實(shí)質(zhì)上就是對(duì)探測器穿越火星大氣時(shí)的軌跡進(jìn)行控制,使得探測器在飛出大氣時(shí)恰好到達(dá)或接近能量最優(yōu)氣動(dòng)捕獲軌道。由前述分析可知,要使氣動(dòng)捕獲軌道為最優(yōu)氣動(dòng)捕獲軌道,則氣動(dòng)捕獲軌道的遠(yuǎn)地點(diǎn)應(yīng)與目標(biāo)軌道相切,飛出時(shí)的航跡角盡量小,且飛出速度與飛出時(shí)的航跡角滿足式(13)。取最優(yōu)控制指標(biāo)為
(15)
式中:
為對(duì)應(yīng)飛出航跡角滿足式(13)的飛出速度。進(jìn)入條件為
(16)
(17)
式中:γfmax為保證捕獲軌道與目標(biāo)軌道相切的最大飛出航跡角。探測器大氣內(nèi)飛行時(shí)的總飛行約束為
(18)
為將氣動(dòng)捕獲軌道優(yōu)化問題轉(zhuǎn)化為多參數(shù)優(yōu)化問題,選取傾側(cè)角σ作為控制變量,將控制變量在配點(diǎn)處進(jìn)行離散,并以這些離散點(diǎn)作為節(jié)點(diǎn)構(gòu)造Lagrange插值多項(xiàng)式來近似擬合控制變量,選取Chebyschev多項(xiàng)式的根作為配點(diǎn),有
τl=cos(πl(wèi)/n)l=0,1,…,n
(19)
式中:n為插值多項(xiàng)式的階數(shù),并且插值點(diǎn)落在[-1,1]之間,對(duì)t∈(t0,tf)進(jìn)行時(shí)域轉(zhuǎn)換
(20)
在這里由于大氣內(nèi)飛行的時(shí)間未知,所以將飛出時(shí)間tf也作為待優(yōu)化的變量。由于控制量σ是關(guān)于時(shí)間t的連續(xù)函數(shù),很難得到其精確表達(dá)式。為此,按照式(19)和式(20)的時(shí)間轉(zhuǎn)換方式對(duì)控制量σ進(jìn)行離散化,即
σl=σ(t(τl))l=0,1,…,n
因此優(yōu)化參數(shù)為
(21)
于是,對(duì)于氣動(dòng)捕獲軌道的優(yōu)化問題就轉(zhuǎn)化為對(duì)大氣內(nèi)飛行時(shí)間tf以及控制序列σ在配點(diǎn)處的離散值的參數(shù)優(yōu)化問題。針對(duì)這個(gè)多參數(shù)優(yōu)化問題,第2節(jié)將提出一種改進(jìn)鴿群算法對(duì)上述待優(yōu)化的參數(shù)向量進(jìn)行優(yōu)化。
鴿群算法是受鴿子覓食的啟發(fā)提出的一種仿生學(xué)算法[16]。當(dāng)鴿子距離目標(biāo)較遠(yuǎn)時(shí),通過感知磁場在腦海中形成地圖,來不斷接近目標(biāo);隨著鴿子不斷接近目的地,其導(dǎo)航工具由磁場變?yōu)槟康牡馗浇牡貥?biāo)。熟悉地標(biāo)的鴿子直接飛往目的地,其他鴿子則跟隨那些熟悉地標(biāo)的鴿子向目的地飛行。
設(shè)某優(yōu)化問題的目標(biāo)函數(shù)為J(ξ),其中優(yōu)化變量為ξ。鴿群算法在處理優(yōu)化問題時(shí)模擬鴿子導(dǎo)航的過程。假設(shè)在優(yōu)化時(shí)需要N只鴿子,每只鴿子的位置Xi,i=1,2,…,N表示優(yōu)化變量ξ的值,在優(yōu)化時(shí)化還需要對(duì)應(yīng)鴿子的速度Vi,i=1,2,…,N。假設(shè)優(yōu)化變量是m維的,則每只鴿子的位置和速度分別記作在采用鴿群算法優(yōu)化時(shí),設(shè)鴿群中有N只鴿子,它們的位置和速度分別為Xi和Vi。該算法求解優(yōu)化問題時(shí)包括2部分迭代過程,分別稱作地圖迭代和地標(biāo)迭代。在地圖迭代中,鴿群中的每一個(gè)體通過種群中的最優(yōu)信息來更新自身的位置和速度,更新公式為
(22)
式中:R為迭代因子;rand為[0,1]區(qū)間上的隨機(jī)值;t=1,2,…,tmax為當(dāng)前迭代次數(shù),tmax為迭代總次數(shù);i=1,2,…,N代表第i只鴿子;Vi(t)為第i個(gè)鴿子在第t次迭代的速度;Xi(t)為第i個(gè)鴿子在第t次迭代的位置;G(t-1)為Xi(t-1) (i=1,2,…,N)中的某一個(gè),其滿足:
當(dāng)?shù)螖?shù)到達(dá)最大迭代次數(shù)的75%以后,進(jìn)入地標(biāo)迭代的過程,選取Xc(t)為第t次迭代鴿群的中心點(diǎn),作為鴿子飛行的參考方向?qū)γ恳恢圾澴舆M(jìn)行更新,更新公式為
(23)
式中:
其中:fitness[Xi(t)]為當(dāng)前迭代下鴿子的適應(yīng)度;N(t)為第t次迭代下種群的數(shù)量;ε為防止分母為零而附加的充分小的數(shù)。
可見,在地標(biāo)迭代過程中,每次迭代舍去部分個(gè)體,既保證了算法的較優(yōu)信息又提高了迭代收斂的速率。原始鴿群算法的全局搜索能力主要體現(xiàn)在地圖迭代的過程中,動(dòng)態(tài)權(quán)值e-RN體現(xiàn)了算法的探索能力,在鴿群向當(dāng)前最優(yōu)值收斂的過程中保持鴿群的多樣性。然而這種方式具有局限性,在鴿群數(shù)目一定的情況下,R值的選取在很大程度上決定了全局搜索的結(jié)果。
圖1為R在不同取值下動(dòng)態(tài)因子的變化情況。可以看出,當(dāng)R取較大值0.3時(shí),迭代到20次e-Rt就已衰減接近為零,之后的迭代過程主要體現(xiàn)在向最優(yōu)值靠近,失去了探索能力,算法易早熟;當(dāng)R取較小值0.002 6時(shí),在迭代完成時(shí)e-Rt仍然具有較大的取值,算法不易收斂;為了平衡二者的關(guān)系,可以通過指定迭代完成時(shí)期望的動(dòng)態(tài)因子的取值自適應(yīng)地選取R值。令迭代完成時(shí)期望的動(dòng)態(tài)因子的值為α,有e-RT=α,則R=-(1/T)lnα, 其中T為迭代次數(shù)。
圖1 R的取值對(duì)動(dòng)態(tài)權(quán)值的影響Fig.1 Influence of R on dynamic weight
通過自適應(yīng)的方式選取R值,雖然能夠在一定程度上緩解易早熟和不易收斂的問題,但是卻不能從根本上很好地平衡全局搜索與局部搜索之間的關(guān)系。在實(shí)際應(yīng)用中,在迭代初期,往往更注重于全局搜索能力,而在迭代后期則希望局部搜索能力更強(qiáng)一些。由圖1可以看出,當(dāng)取α=0.01 時(shí),在迭代后期動(dòng)態(tài)因子較小,算法主要在最優(yōu)值附近進(jìn)行局部搜索,然而在迭代初期,動(dòng)態(tài)因子下降很快,并不能很好地進(jìn)行全局搜索;當(dāng)取α=0.4時(shí),雖然使得動(dòng)態(tài)因子下降速度減緩,全局搜索能力增強(qiáng),但是卻削弱的算法的局部搜索能力。
本節(jié)在原有鴿群算法的負(fù)指數(shù)動(dòng)態(tài)權(quán)值的基礎(chǔ)上,為了平衡全局搜索與局部搜索提出一種廣義的動(dòng)態(tài)權(quán)值項(xiàng)。即將式(23)修改為
(24)
式中:
(25)
且有
0≤f(t)≤1
(26)
為了更加方便直觀地確定參數(shù)ω、b的取值, 引入另外2組參數(shù)(t1,α1)、(t2,α2),令
(27)
表示當(dāng)?shù)螖?shù)為t1時(shí),動(dòng)態(tài)權(quán)值取值為α1;當(dāng)?shù)螖?shù)為t2時(shí),動(dòng)態(tài)權(quán)值取值為α2。即
(28)
求解可得
(29)
因此,可以通過指定k、(t1,α1)、(t2,α2)來唯一確定動(dòng)態(tài)權(quán)值。
下面對(duì)各參數(shù)的取值情況進(jìn)行討論。對(duì)式(25) 求導(dǎo)可得
(30)
當(dāng)k≥1時(shí),對(duì)式(25)取極限可得
(31)
此時(shí)k、(t1,α1)、(t2,α2)可任意取值都可以滿足式(26)的約束。
當(dāng)k<1時(shí),此時(shí)k的取值與(t1,α1)、(t2,α2)相互制約,由式(26)可得
(32)
即
1-e-ω b≤k
(33)
又由式(29)可得
(34)
(35)
不妨取t1=T/2、t2=T,則有
(36)
(37)
即
(38)
又
(39)
(40)
(41)
(42)
由式(40)可知此時(shí)k≥0。當(dāng)取k=0時(shí),由式(29) 有b=0,又由式(25)可得
(43)
此時(shí)改進(jìn)算法退化為原始的鴿群算法,也就是說原始鴿群算法是改進(jìn)后算法的一種特殊情況。
改進(jìn)后的鴿群算法能夠通過調(diào)整動(dòng)態(tài)權(quán)值很好地平衡全局搜索與局部搜索之間的關(guān)系,既保證了算法的收斂性,又大大加強(qiáng)了算法的探索能力。圖2為改進(jìn)后的鴿群算法在k、(t1=N/2,α1)、(t2=N,α2=0.1)的不同取值下動(dòng)態(tài)權(quán)值的變化情況,其中N=200。通過調(diào)節(jié)α1和k的值可得到不同的動(dòng)態(tài)權(quán)值曲線。
可以看出,當(dāng)k=0且b=0時(shí),算法即為原始的鴿群算法,即原算法是改進(jìn)后的一種特殊情況。當(dāng)k≠0時(shí),可以通過調(diào)節(jié)(t1,α1)、(t2,α2)選擇合適的動(dòng)態(tài)權(quán)值曲線,在搜索前期使動(dòng)態(tài)權(quán)值緩慢下降,保持搜索活力,提高全局搜索的能力;在搜索后期,使得動(dòng)態(tài)權(quán)值取得較小的值,保證算法的收斂。相比于原算法,改進(jìn)后的算法不僅兼容了原算法中指數(shù)形式下降的動(dòng)態(tài)權(quán)值,而且對(duì)其進(jìn)行了擴(kuò)展,在形式選擇上提供了更大的自由度,提升了全局搜索與局部搜索之間的平衡。
圖2 動(dòng)態(tài)權(quán)值隨迭代次數(shù)的變換Fig.2 Evolution of dynamic weights
采用改進(jìn)的鴿群算法對(duì)氣動(dòng)捕獲軌道進(jìn)行優(yōu)化。鴿群中的每一只鴿子所在的位置代表所求問題的一個(gè)可能解,其初值在搜索空間內(nèi)隨機(jī)生成。選擇鴿群中鴿子的數(shù)目為m,則鴿群的位置X可表示為
(44)
式中:Xk(k=1,2,…,m)為每一只鴿子當(dāng)前所得到的待優(yōu)化變量的值,包括探測器飛出大氣的時(shí)間tf,以及以傾側(cè)角為控制序列在配點(diǎn)處的離散值。鴿群的更新速度V與X同維,首先對(duì)每一只鴿子所優(yōu)化得到的控制量進(jìn)行Lagrange插值得到控制序列,然后根據(jù)所得到的控制序列和優(yōu)化得到的飛出時(shí)間tf對(duì)動(dòng)力學(xué)方程式(3)進(jìn)行積分得到各個(gè)狀態(tài)的值,計(jì)算此時(shí)的目標(biāo)函數(shù)值,然后進(jìn)行迭代尋優(yōu)。
鴿群算法本身并不具備處理約束的能力,在進(jìn)行迭代尋優(yōu)的過程中為了保證滿足過程約束和終端約束,將這些約束條件以罰函數(shù)的方式加到目標(biāo)函數(shù)當(dāng)中,經(jīng)過修正的目標(biāo)函數(shù)為
(45)
表1 火星探測器參數(shù)Table 1 Parameters of Mars explorers
表2 火星參數(shù)Table 2 Parameters of Mars
采用改進(jìn)鴿群算法,取傾側(cè)角為五階插值多項(xiàng)式時(shí),改進(jìn)鴿群算法的迭代收斂曲線以及優(yōu)化結(jié)果如圖3所示,圖3(a)為探測器飛行高度隨時(shí)間變化的情況,可見進(jìn)入高度和飛出高度相等,都等于火星大氣高度,并且最小飛行高度大于所要求的35 000 m。圖3(b)為飛行速度隨時(shí)間變化的曲線,飛出速度與飛出航跡角基本滿足式(13)。圖3(c)為航跡角隨時(shí)間變化的曲線,圖3(d)和圖3(e) 為路徑約束,分別表示過載和熱流率隨時(shí)間的變化情況,可見大氣內(nèi)飛行過程可以滿足過程約束,圖3(f)為控制變量隨時(shí)間變化的曲線。綜上,通過改進(jìn)鴿群算法,以傾側(cè)角為控制變量可以得到一條接近能量最優(yōu)的氣動(dòng)捕獲軌道。
圖3 氣動(dòng)捕獲軌道優(yōu)化結(jié)果Fig.3 Performance of optimized aerocapture orbit
1) 在對(duì)火星探測器的氣動(dòng)捕獲過程分析的基礎(chǔ)上,從脈沖變軌的角度給出了氣動(dòng)捕獲軌道的優(yōu)化指標(biāo)。
2) 對(duì)原始鴿群算法的局限性進(jìn)行分析,提出一種能夠平衡收斂速度與全局搜索關(guān)系的改進(jìn)算法。
3) 將氣動(dòng)捕獲軌道優(yōu)化問題轉(zhuǎn)化為多參數(shù)優(yōu)化問題,并利用提出的改進(jìn)個(gè)鴿群算法對(duì)其進(jìn)行優(yōu)化。