基于普適變量法的雙沖量行會(huì)優(yōu)化建模與仿真

基于普適變量法的雙沖量行會(huì)優(yōu)化建模與仿真

0最優(yōu)雙沖量廟會(huì)空間交流技術(shù)是空間操作的關(guān)鍵技術(shù)之一，它是空間系統(tǒng)、空間平臺(tái)和空間運(yùn)輸系統(tǒng)的安裝、維護(hù)、配件、維護(hù)和救援等先決條件。最初,人們主要考慮單圈Lambert雙沖量交會(huì)問題,即在固定時(shí)間、給定始末位置下,在始末時(shí)刻(亦即始末位置)分別對(duì)追蹤航天器施加脈沖,使其在單圈下與目標(biāo)航天器對(duì)接。隨著研究的逐步深入,多圈Lambert最優(yōu)雙沖量交會(huì)問題倍受關(guān)注。例如,文獻(xiàn)針對(duì)圓軌道交會(huì)問題,提出了一種多圈Lambert最優(yōu)雙沖量交會(huì)的求解方法;文獻(xiàn)分析了Lambert交會(huì)問題中轉(zhuǎn)移時(shí)間與轉(zhuǎn)移軌道參數(shù)及變軌速度增量之間的關(guān)系;文獻(xiàn)提出了一種多圈Lambert變軌的求解算法,并把它用于尋求航天器交會(huì)的燃料最優(yōu)軌道。此外,為了進(jìn)一步優(yōu)化所需能量,在給定始末位置、固定時(shí)間下,人們開始尋找兩次脈沖的最佳施加時(shí)刻使得完成交會(huì)所需的沖量最小。例如,文獻(xiàn)基于CW方程研究了始末軌道為圓軌道且第一次施加脈沖時(shí)刻不確定的最優(yōu)雙沖量交會(huì)問題;文獻(xiàn)利用遺傳算法研究了第一次施加脈沖時(shí)刻不確定的雙沖量軌道轉(zhuǎn)移的數(shù)值求解;文獻(xiàn)與針對(duì)始末軌道為圓軌道的情況,提出了在將固定交會(huì)時(shí)間劃分為三段(即最初漂移時(shí)間段、軌道轉(zhuǎn)移時(shí)間段、終端停泊時(shí)間段)的前提下分析能量最優(yōu)問題的想法。本文基于普適變量研究了兩個(gè)共面軌道的最優(yōu)雙沖量交會(huì)問題。具體地,在給定交會(huì)時(shí)間、交會(huì)始末位置且在始末位置施加脈沖的條件下,我們考慮了在將給定時(shí)間段劃分為三階段(即初始飄移階段、軌道轉(zhuǎn)移階段與終端停泊階段)的前提下,圓軌道和拱線相同的橢圓軌道的軌道轉(zhuǎn)移階段為單圈轉(zhuǎn)移的最優(yōu)雙沖量交會(huì)問題,并對(duì)之進(jìn)行數(shù)學(xué)建模和數(shù)值求解。1la東南角調(diào)整模型我們首先回顧一下Lambert雙沖量交會(huì)問題。為此,我們假設(shè)追蹤航天器C的軌道與目標(biāo)航天器A的軌道共面且拱線相同、在初始時(shí)刻兩航天器在各自的開普勒軌道上同按逆(順)時(shí)針運(yùn)行,且建立如下坐標(biāo)系(注:在全文中我們將采用統(tǒng)一的坐標(biāo)系,因此在下面幾節(jié)中我們不再重新建立坐標(biāo)系):以地心為坐標(biāo)原點(diǎn),x軸指向近拱點(diǎn),y軸與x軸構(gòu)成右手系。進(jìn)一步假設(shè):在t=0時(shí)刻,追蹤航天器C的位置為→rC0、速度為→vC0;兩航天器在t=T時(shí)刻交會(huì),此時(shí)目標(biāo)航天器A的位置為→rAΤ、速度為→vAΤ。Lambert雙脈沖交會(huì)問題就是在雙脈沖作用下,求解追蹤航天器C在t=0時(shí)刻所需的速度增量→Δv1和t=T時(shí)刻所需的速度增量→Δv2,使得兩航天器在t=T時(shí)刻、→rAΤ處交會(huì)。求解Lambert交會(huì)問題的一個(gè)經(jīng)典求法是Gauss算法。具體地,假設(shè)追蹤航天器C在轉(zhuǎn)移軌道上→rC0處的速度為→v1、r→AΤ處的速度為v→2,則基于普適變量z?v→1和v→2由下列方程決定:v→1=r→AΤ-fr→C0g?v→2=g˙r→AΤ-r→C0g(1)其中u為地球萬有引力常數(shù),Δθ為r→C0與r→AΤ間小于π弧度的夾角。由公式(1)可得:Δv→1=v→1-v→C0?Δv→2=v→AΤ-v→2提一下,上面的D的表達(dá)式有正負(fù)號(hào)。這是因?yàn)榻粫?huì)可以通過“短程”實(shí)現(xiàn)或者“長程”實(shí)現(xiàn)。它們分別對(duì)應(yīng)著追蹤航天器C從r→C0首次到達(dá)r→AΤ的角度改變量為Δθ或從r→C0首次到達(dá)r→AΤ的角度改變量為(2π-Δθ)的情況。另外,由于交會(huì)發(fā)生在T時(shí)刻,Lambert交會(huì)問題還應(yīng)滿足:uΤ=(y(z)B)3S+Dy(z)文獻(xiàn)給出了一個(gè)求解多圈Lambert最優(yōu)交會(huì)問題的直接方法。這里,我們將借助普適變量法建立一個(gè)關(guān)于普適變量z的全局優(yōu)化模型,使其覆蓋多圈Lambert交會(huì)優(yōu)化問題。為此,我們首先引入一個(gè)選擇因子d∈{-1,1},使得d=1對(duì)應(yīng)的交會(huì)任務(wù)是通過“短程”轉(zhuǎn)移軌道來實(shí)現(xiàn)、d=-1對(duì)應(yīng)的交會(huì)任務(wù)是通過“長程”轉(zhuǎn)移軌道來實(shí)現(xiàn),從而將D的計(jì)算式改進(jìn)為D=drC0rAΤsinΔθ1-cosΔθ。進(jìn)一步,令ΔV(z?d)=|v→C0-r→AΤ-fr→C0g|+|v→AΤ-g˙r→AΤ-r→C0g|?其中z∈R+,d∈R,我們把優(yōu)化問題轉(zhuǎn)化為考慮目標(biāo)函數(shù)ΔV(z,d)在約束條件為d2-1=0?z＞0?D≠0?uΤ=(y(z)B)3S+Dy(z)下的最小值問題。記這個(gè)問題為優(yōu)化模型(1)。顯然地,通過求解上述含有約束條件的優(yōu)化模型(注:關(guān)于優(yōu)化模型的數(shù)值求解將在第四節(jié)介紹),我們將得到滿足速度增量最優(yōu)的全局解。其全局解本質(zhì)上就是文獻(xiàn)中的最優(yōu)解。2各航天器各時(shí)間t時(shí)的雙沖量優(yōu)化模型在這一節(jié)中,我們將針對(duì)始末軌道為圓軌道的情況,在給定交會(huì)時(shí)間、交會(huì)地點(diǎn)的前提下,通過尋求施加雙脈沖的時(shí)刻進(jìn)一步優(yōu)化單圈下兩航天器雙沖量交會(huì)問題。首先,假設(shè)在t=0時(shí)刻追蹤航天器C和目標(biāo)航天器A分別運(yùn)行在兩個(gè)共面的圓軌道上,半徑分別為rC和rA,rC≤rA,角速度分別為wC和wA且同為逆時(shí)針或瞬時(shí)針,追蹤航天器C的位置為r→C0、速度為v→C0、真近點(diǎn)角為θC0(0≤θC0<2π),目標(biāo)航天器A的位置為r→A0、真近點(diǎn)角為θA0(0≤θA0<2π);在t=T時(shí)刻,目標(biāo)航天器A的位置為r→AΤ、速度為v→AΤ、真近點(diǎn)角為θAΤ=θA0+ΤurA3,這里T為最終交會(huì)時(shí)刻,r→AΤ為交會(huì)地點(diǎn)。接著,為了進(jìn)一步優(yōu)化,我們將整個(gè)交會(huì)過程劃分為三個(gè)階段:初始漂移階段、軌道轉(zhuǎn)移階段、終端停泊階段。為此,我們引入兩個(gè)新的時(shí)間變量t1,t2,當(dāng)t∈[0,t1)時(shí),追蹤航天器C和目標(biāo)航天器A運(yùn)行在各自初始軌道上;當(dāng)t=t1時(shí),施加脈沖,使得追蹤航天器C開始執(zhí)行軌道轉(zhuǎn)移;當(dāng)t∈(t1,t1+t2)時(shí),追蹤航天器C沿著新的軌道運(yùn)行;當(dāng)t=t1+t2時(shí),施加脈沖,使得追蹤航天器C與目標(biāo)航天器A提前實(shí)現(xiàn)交會(huì);當(dāng)t∈(t1+t2,T]時(shí),兩航天器運(yùn)行在目標(biāo)航天器A的軌道上,共同到達(dá)r→AΤ。下面我們將基于時(shí)間變量t1,t2建立優(yōu)化數(shù)學(xué)模型,從而尋找施加脈沖力的最佳時(shí)刻t1與t2+t1,使得兩航天器交會(huì)過程速度增量最小。經(jīng)過數(shù)學(xué)上的嚴(yán)格推導(dǎo),我們得出:若給定的交會(huì)時(shí)間T滿足Τ≥th=wCΤz2+θA0-θC0-ΝπwC-wA,其中Τz=2π(rA+rC)38u?Ν∈{-3,-1,1}是滿足關(guān)系式t1=wAΤz2+θA0-θC0-ΝπwC-wA≥0的最大整數(shù),則使得速度增量最小的轉(zhuǎn)移軌道為霍曼轉(zhuǎn)移軌道,其在整個(gè)交會(huì)過程的所需的沖量為:ΔVmin=|u(2rC-2rC+rA)-urC|+|u(2rC+rA-2rA)-urA|此時(shí)t2=Τz2。若T<th,則不存在霍曼轉(zhuǎn)移。此時(shí)其最佳時(shí)刻t1與t2+t1很難確定。下面,我們主要針對(duì)這種情況建立雙沖量優(yōu)化模型。令ΔV(z,t1,t2,d)=|v→Ct1-r→A(t1+t2)-fr→Ct1g|+|v→A(t1+t2)-g˙r→A(t1+t2)-r→Ct1g|其中r→A(t1+t2)為目標(biāo)航天器A在t1+t2時(shí)刻的位置,v→A(t1+t2)為目標(biāo)航天器A在t1+t2時(shí)刻的速度,r→Ct1為追蹤航天器C在t1時(shí)刻的位置,v→Ct1為追蹤航天器C在t1時(shí)刻的速度,f=1-y(z)rC0?g=Dy(z)u?g˙=1-y(z)rAΤ?y(z)=rC+rA-D(1-zS)B?B=1-coszz?D=drCrAsinΔθ1-cosΔθ,S=z-sinzz3,Δθ為r→Ct1與r→A(t1+t2)兩矢量間小于等于π弧度的夾角。根據(jù)上述標(biāo)號(hào),我們把優(yōu)化問題轉(zhuǎn)化為考慮求解目標(biāo)函數(shù)ΔV(z,t1,t2,d)在約束條件為d2-1=0、ut2=(yB)3S+Dy、t1+t2≤T、t1≥0、t2≥0、0<z<(2π)2下的最小值問題。記這個(gè)問題為優(yōu)化模型(2)。注意到,優(yōu)化模型(2)中的目標(biāo)函數(shù)存在奇點(diǎn)g=0。同時(shí),我們又注意到,當(dāng)g趨向于0時(shí),ΔV(z,t1,t2,d)的極限存在。因此,我們可以在g=0處補(bǔ)充定義,使得ΔV(z,t1,t2,d)在g=0處連續(xù)。然后,根據(jù)函數(shù)的連續(xù)性,我們可加限制|g|≥δ(其中δ為一小常數(shù),其值由數(shù)值求解所允許的誤差來適當(dāng)選擇),用ΔV(z,t1,t2,d)在|g|=δ的取值來近似代替ΔV(z,t1,t2,d)在奇點(diǎn)g=0處的取值,從而去掉奇點(diǎn)。按照以上思路,我們把優(yōu)化問題近似轉(zhuǎn)化為求解目標(biāo)函數(shù)ΔV(z,t1,t2,d)在約束條件為d2-1=0、ut2=(yB)3S+Dy、t1+t2≤T、t1≥0、t2≥0、0<z<(2π)2、|g|≥δ(其中δ為一小常數(shù),其值由數(shù)值求解所允許的誤差來適當(dāng)選擇)的情況下的最小值問題。記這個(gè)問題為優(yōu)化模型(2a)。(在模型(2)中,若將t1+t2≤T改為t1+t2=T,則模型變?yōu)橹挥谐跏计齐A段和軌道轉(zhuǎn)移階段的雙沖量交會(huì)對(duì)接的優(yōu)化模型;若令t1=0且將t1+t2≤T改為t2≤T,則模型變?yōu)橹挥熊壍擂D(zhuǎn)移階段和終端漂移階段的雙沖量交會(huì)對(duì)接的優(yōu)化模型;若令t1=0,t2=T,且g≠0,則模型(2)變?yōu)閱稳ο碌哪Ｐ?1)。)從優(yōu)化模型的建立過程我們可以看到,優(yōu)化模型(2)并不涉及到T與th的關(guān)系問題。因此優(yōu)化模型(2)也包含了T≥th的最優(yōu)能量情況。易見,優(yōu)化模型(1)是一個(gè)二維非線性規(guī)劃數(shù)學(xué)模型,變量為z、d;而優(yōu)化模型(2a)是一個(gè)四維非線性規(guī)劃數(shù)學(xué)模型,變量為z、d、t1、t2。由于優(yōu)化問題是一個(gè)NP問題,隨著變量個(gè)數(shù)的增加,其優(yōu)化求解問題的計(jì)算復(fù)雜度會(huì)隨著指數(shù)增加。同樣,我們將借助求解非線性規(guī)劃的數(shù)學(xué)軟件lingo進(jìn)行求解。其求解問題將在第4節(jié)介紹。3追蹤航天器a的目標(biāo)函數(shù)在這一節(jié)中,我們將針對(duì)始末軌道為橢圓軌道的情況,在給定交會(huì)時(shí)間、交會(huì)地點(diǎn)的前提下,通過尋求施加雙脈沖力的時(shí)刻進(jìn)一步優(yōu)化單圈下兩航天器最優(yōu)雙沖量交會(huì)問題。首先,假設(shè)在t=0時(shí)刻,追蹤航天器C和目標(biāo)航天器A都運(yùn)行在兩拱線相同的橢圓軌道上,其軌道方程分別為rC=aC(1-eC2)1+eCcosθC和rA=aA(1-eA2)1+eAcosθA,角速度分別為和θ˙A=uaA3(1-eA2)3(1+eAcosθA)2,并且追蹤航天器C的位置為r→C0=aC(1-eC2)(1+eCcosθC0)、真近點(diǎn)角為θC0、偏近點(diǎn)角為EC0,目標(biāo)航天器A的位置為r→A0=aA(1-eA2)(1+eAcosθA0)、真近點(diǎn)角為θA0、偏近點(diǎn)角為EA0,其中0≤θC0,EC0,θA0,EA0<2π;在t=T時(shí)刻,目標(biāo)航天器A的位置為r→AΤ、真近點(diǎn)角為θAT、速度為v→AΤ,最終交會(huì)地點(diǎn)為r→AΤ。接著,我們將追蹤航天器C和目標(biāo)航天器A的真近點(diǎn)角表達(dá)成關(guān)于時(shí)間的近似函數(shù):θC=ΜC+(2eC-0.25eC3)sinΜC+1.25eC2sin(2ΜC)+1312×eC3sin(3ΜC)+Ο(eC4)θA=ΜA+(2eA-0.25eA3)sinΜA+1.25eA2sin(2ΜA)+1312eA3sin(3ΜA)+Ο(eA4)其中ΜC=uaC3t+EC0-sinEC0為追蹤航天器C在t時(shí)刻的平近點(diǎn)角,ΜA=uaA3t+EA0-sinEA0為目標(biāo)航天器A在t時(shí)刻的平近點(diǎn)角。根據(jù)這些標(biāo)號(hào),我們把優(yōu)化問題轉(zhuǎn)化為求解目標(biāo)函數(shù)ΔV(z,t1,t2,d)在約束條件為d2-1=0、ut2=(yB)3S+Dy、t1+t2≤T、t1≥0、t2≥0、0<z<(2π)2的情況下的最小值問題。記這個(gè)問題為優(yōu)化模型(3)(注:類似于第二節(jié),我們可以通過改變模型(3)的約束條件得到三種不同的優(yōu)化模型)。同樣地,優(yōu)化模型(3)中的目標(biāo)函數(shù)存在奇點(diǎn)g=0。按照前一節(jié)提到的奇點(diǎn)處理思路,我們類似地把優(yōu)化問題近似轉(zhuǎn)化為考慮求解目標(biāo)函數(shù)ΔV(z,t1,t2,d)在約束條件為d2-1=0、ut2=(yB)3S+Dy、t1+t2≤T、t1≥0、t2≥0、0<z<(2π)2、|g|≥δ(其中δ為一小常數(shù),其值由數(shù)值求解所允許的誤差來適當(dāng)選擇)的最小值問題。記這個(gè)問題為優(yōu)化模型(3a)。顯然,模型(3a)也是一個(gè)四維問題,其變量分別為t1,t2,z,d。我們也同樣地借助求解非線性規(guī)劃的優(yōu)秀數(shù)學(xué)軟件Lingo進(jìn)行求解。4航天器a運(yùn)行在變軌前運(yùn)行在前面三節(jié)中,我們建立了最優(yōu)雙脈沖交會(huì)的三個(gè)非線性規(guī)劃數(shù)學(xué)模型,下面我們將利用非線性規(guī)劃軟件Lingo分別對(duì)其進(jìn)行數(shù)值求解。Lingo是一個(gè)利用線性規(guī)劃和非線性規(guī)劃來簡潔地闡述、解決和分析復(fù)雜問題的簡便工具。其特點(diǎn)是程序執(zhí)行速度很快,易于輸入、修改、求解和分析一個(gè)數(shù)學(xué)規(guī)劃問題。關(guān)于Lingo的詳細(xì)介紹請(qǐng)參見文獻(xiàn)。首先,為了便于計(jì)算,在整個(gè)計(jì)算過程中,我們約定DU與TU分別為距離單位和時(shí)間單位,即1DU=6378145米,1TU=806.8118744秒。此時(shí)萬有引力常數(shù)為u=1DU3/TU2。下面,我們將分別針對(duì)共面圓軌道和拱線相同的共面橢圓軌道上兩航天器交會(huì)的具體實(shí)例,借助數(shù)學(xué)規(guī)劃軟件Lingo對(duì)模型(2a)與(3a)進(jìn)行求解,給出雙沖量交會(huì)的全局最優(yōu)解。(1)假設(shè)追蹤航天器C在變軌前運(yùn)行在半徑為rC=1.3DU的圓軌道上,在t=0時(shí)刻的位置為r→C0=[1.3,0]DU,真近點(diǎn)角為0;目標(biāo)航天器A運(yùn)行在半徑為rA=1.5DU的圓軌道上,在t=0時(shí)刻的位置為r→A0=[0.7500,1.2990]DU,真近點(diǎn)角3.14159263。①假設(shè)交會(huì)時(shí)間為T=15TU,則最終交會(huì)位置為r→AΤ=[-1.46622,0.31653]DU,真近點(diǎn)角為9.21216,代入優(yōu)化模型(2a),得到四種情況(即同時(shí)含有初始漂移時(shí)間段和終端時(shí)間段、只含初始漂移時(shí)間段、只含終端漂移時(shí)間段、不含初始漂移時(shí)間段和終端時(shí)間段)的數(shù)值結(jié)果(見表1)。②設(shè)交會(huì)時(shí)間為T=3TU,則交會(huì)位置為r→AΤ=[-1.34314,0.667806]DU,真近點(diǎn)角為2.68019,代入優(yōu)化模型(2),得到四種情況(即同時(shí)含有初始漂移時(shí)間段和終端時(shí)間段、只含初始漂移時(shí)間段、只含終端漂移時(shí)間段、不含初始漂移時(shí)間段和終端時(shí)間段)的數(shù)值結(jié)果(見表2)。(2)假設(shè)追蹤航天器C在變軌前運(yùn)行在橢圓軌道rC=1.22×(1-0.12)(1+0.1cosθC),在t=0時(shí)刻的位置為r→C0=[1.098,0]DU,真近點(diǎn)角為0;目標(biāo)航天器A運(yùn)行在半徑為rA=1.26×(1-0.152)1+0.15cosθA的橢圓軌道上,在t=0時(shí)刻的位置為r→A0=[0.902192,0.622872]DU,偏近點(diǎn)角為3.14159266。①設(shè)交會(huì)時(shí)間為T=10TU,則交會(huì)位置為r→AΤ=[0.045245,1.224027]DU,真近點(diǎn)角為7.817,代入優(yōu)化模型(3a),得到四種情況(即同時(shí)含有初始漂移時(shí)間段和終端時(shí)間段、只含初始漂移時(shí)間段、只含終端漂移時(shí)間段、不含初始漂移時(shí)間段和終端時(shí)間段)的數(shù)值結(jié)果(見表3)。②設(shè)交會(huì)時(shí)間為T=3TU,則交會(huì)位置為r→AΤ=[-1.294923,0.5969354]DU,真近點(diǎn)