利用高斯偽譜法求解小推力伴星最優(yōu)釋放軌跡

段傳輝 董云峰

（北京航空航天大學(xué)宇航學(xué)院，北京100191）

1 引言

航天器編隊飛行可以完成單個航天器所不能完成的任務(wù)，近年來成為國內(nèi)外航天領(lǐng)域的研究重點［1］。

航天器相對運動軌跡控制是編隊飛行的重要內(nèi)容，傳統(tǒng)方法一般以脈沖變軌完成。近年來由于小推力發(fā)動機(jī)技術(shù)逐漸成熟，以效率高、壽命長的特點得到充分的重視［2-4］。

求解航天器最優(yōu)軌跡問題的方法通常可分為直接法和間接法［5-7］，間接法利用Pontryagin原理，將最優(yōu)控制問題轉(zhuǎn)化為常微分方程的兩點邊值問題；間接法可以獲得很高的精度，但是因初值猜測很困難，收斂半徑小使得兩點邊值問題很難求解［5］。直接法的思想是將狀態(tài)量和控制量離散化，將最優(yōu)控制問題轉(zhuǎn)化為非線性規(guī)劃（NLP）問題，然后用常規(guī)的NLP求解方法，如系列二次規(guī)劃方法（SQP）求得近似解。直接法有許多種，各種方法在計算量、收斂性等方面各有差異，近年來由David Benson和Huntington提出的高斯偽譜法對平滑型最優(yōu)控制問題有很好的收斂性，以較少節(jié)點即可獲得較高的精度［7-8］。國內(nèi)的學(xué)者利用高斯偽譜法求解了深空探測、高超聲速飛行器的再入、月球軟著陸等最優(yōu)軌跡問題［6，9-10］。

本文針對雙星編隊的伴飛星的釋放任務(wù)，提出基于小推力變軌完成的從主星投放至形成穩(wěn)定繞飛構(gòu)型的燃料最省釋放軌跡，利用高斯偽譜法將C-W方程描述的連續(xù)型最優(yōu)控制問題轉(zhuǎn)化為NLP問題，再利用SQP方法求解此NLP問題并給出了仿真算例。

2 相對運動描述及穩(wěn)定繞飛條件

考慮目標(biāo)航天器在圓軌道或者近圓軌道上運行，追蹤航天器與目標(biāo)航天器的位置誤差遠(yuǎn)小于圓的半徑，速度差遠(yuǎn)小于目標(biāo)航天器的在軌運行速度。定義相對運動坐標(biāo)系，其原點位于目標(biāo)航天器質(zhì)心，x軸指向速度方向，y軸指向目標(biāo)飛行器矢徑方向，z軸與x和y軸構(gòu)成右手正交系。

式中A和B的定義為

考慮航天器不受軌控的繞飛情形，即令式（1）中T為零，可解出相對位置隨時間的變化關(guān)系［11］為

從式（2）可以看出，z方向和x，y是相互獨立的，可單獨分析。此外，三軸中僅x（t）中含有隨時間線性增長的項為形成長期穩(wěn)定的繞飛，則必須滿足條件

此時，x-y平面內(nèi)的方程可轉(zhuǎn)化為

由此可見，飛行器的相對運動軌跡在x-y平面投影為一橢圓，其半長軸為半短軸為半長軸的一半，進(jìn)一步分析可知，相對運動軌跡在y-z，x-z平面的投影也為橢圓，但其長軸和短軸不一定與坐標(biāo)軸重合［12］。

上述橢圓的圓心在x軸還有常值偏移量，為達(dá)到伴星對目標(biāo)星的較好觀測效果，可進(jìn)一步令

只要航天器的相對位置與相對速度初始條件滿足式（3）和式（5），在近距離且不考慮攝動情況下，伴星相對目標(biāo)星的軌跡為一個空間橢圓，且橢圓中心在目標(biāo)星質(zhì)心。

小型的伴星一般由主星搭載，入軌后通過航天員或者主星上的分離機(jī)構(gòu)釋放，初始相對位置Δr≈0，相對速度Δv≈0，此狀態(tài)一般不滿足式（3）和式（5）所確定的穩(wěn)定的橢圓軌跡繞飛的初始條件，因此伴星投放后必須實施變軌，才能進(jìn)入穩(wěn)定的繞飛狀態(tài)。

傳統(tǒng)的變軌一般由脈沖控制來完成，只需進(jìn)行有限次的脈沖即可進(jìn)入目標(biāo)軌道，此種方式計算簡便，但是誤差較大，不能精確調(diào)整。利用有限小推力發(fā)動機(jī)可以得到比脈沖控制更為精確的變軌，且可以充分利用已有的最優(yōu)控制理論，對變軌過程進(jìn)行優(yōu)化。本文選用推力恒定，方向可變的常值小推力發(fā)動機(jī)完成伴星的釋放軌跡控制，則關(guān)于求解推力器噴氣方向單位矢量u＝［ux，uy，uz］T的最優(yōu)控制問題可描述如下。

邊界點約束：x（t0）＝［r0，v0］，x（tf）＝［rf，vf］，在初始狀態(tài)選取時注意剛釋放時r0≈0，v0≈0，在終端時刻，即變軌結(jié)束點（橢圓繞飛構(gòu)型的入口點處），rf一般取為102m～103m的級別，而vf是根據(jù)rf的取值及式（3）和式（5）的約束條件來確定的。

過程約束：由于控制量u是方向矢量，必須滿足

3 高斯偽譜法求解最優(yōu)控制問題

高斯偽譜法的基本思想是通過將隨時間連續(xù)變化的狀態(tài)量和控制量在有限個數(shù)的時間點進(jìn)行離散，用這些離散點上Lagrange差值多項式近似表達(dá)狀態(tài)量和控制量函數(shù)，再利用Gauss數(shù)值積分將積分約束轉(zhuǎn)化為代數(shù)求和約束，最終將最優(yōu)控制問題轉(zhuǎn)化為NLP問題來求解。

設(shè)x（t）∈Rn為狀態(tài)變量，u（t）∈Rm為控制變量，t0和tf為起始和結(jié)束時刻，則一個典型的Bolza型的最優(yōu)控制問題［8］可描述為確定u（t），使以下指標(biāo)函數(shù)為最小。

式中Φ為終端指標(biāo)；g為積分型指標(biāo)。式（6）中各變量滿足以下約束

式中ψ為等式約束；C為不等式約束。由于Gauss積分定義在［-1，1］，需將以上t∈［t0，tf］描述變換為τ∈［-1，1］，變換方法為

對于一個定義在［-1，1］上的函數(shù)f（t），取Legendre-Gauss多項式的K個零點（LG點）τ1～τK，則可利用Gauss積分公式將積分求解轉(zhuǎn)化為近似求和運算

式中ωi為第i處的高斯權(quán)重，ωi和τi只與個數(shù)K有關(guān)，求解方法參見文獻(xiàn)［7］。

利用起始點τ＝-1及K個LG零點共計K＋1個點，對x（t）進(jìn)行K＋1次多項式Lagrange差值［8］得到：

Lagrange差值可以保證在K＋1個差值點處有x（τi）＝X（τi），其他點處二者則近似相等。

同理，對u（t）在K個LG零點處進(jìn)行Lagrange差值（因差值點個數(shù)不同

式中Dki為微分矩陣，僅與LG點的個數(shù)K有關(guān)，可在算法開始前離線確定，其定義［8］為

針對伴星釋放最優(yōu)軌跡控制的具體問題，設(shè)起始時刻t0＝0，狀態(tài)變量為x＝［x，y，z，vx，vy，vz］T，控制變量為推力器方向矢量u＝［ux，uy，uz］T，動力學(xué)微分方程的右端函數(shù)為f（x，u）。將連續(xù)量的x（τ）和u（τ）在［-1，1］區(qū)間進(jìn)行離散，其中x（τ）的離散點取為起始點、終止點及K個LG點，構(gòu)成一個6（K＋2）維的未知向量［X0，X1，…，XK，Xf］，控制量u（τ）的離散點選取在K個LG點處，構(gòu)成一個3K維的未知向量［U1，U2，…，UK］。

根據(jù)高斯偽譜法，離散變換之后，連續(xù)型的最優(yōu)控制問題的微分方程約束可化作6K維的代數(shù)約束Rk，其表達(dá)式為

邊界條件為ψ（x（t0），t0，x（tf），tf）＝0，即上文所述的釋放初始及結(jié)束時刻的位置、速度應(yīng)滿足的條件，其維數(shù)為12。

過程約束為C（Xk，Uk，τk，tf）≤0，（k＝1，2，…，K），即每個LG點處Uk模為1的條件，其維數(shù)是6K。

指標(biāo)函數(shù)為min（J），J＝tf。

至此，伴星釋放最優(yōu)控制問題轉(zhuǎn)化為針對9K＋13維的未知量Y＝［X0，…，XK，Xf，U1，U2，…，UK，tf］，在滿足Rk＝0，Rf＝0，ψ＝0，C≤0條件下的一個非線性規(guī)劃問題，可用現(xiàn)有的非線性規(guī)劃問題求解方法進(jìn)行求解，本文選用效率比較高的系列二次規(guī)劃法求解。

由于高斯偽譜法求解得到的控制量是有限點處的近似解，而完整的控制問題需要知道任意時刻的控制量，本文選用三次樣條差值法在所求解得到的離散點處的控制量［U1，U2，…，UK］上進(jìn)行差值，以此獲得控制量在整個時間歷程上隨時間連續(xù)變化的函數(shù)u（t）。

4 仿真算例

考慮目標(biāo)飛行器軌道高度為h＝500km的圓軌道，伴星質(zhì)量m＝200kg，推力器推力T＝1N，高斯偽譜法計算節(jié)點K＝100，起始點的相對位置x，y，z為2m，相對速度Vx，Vy，Vz為0.01m／s。

終端時刻的伴星在主星的軌道坐標(biāo)系下的位置即橢圓繞飛入口的坐標(biāo)為［1 000m，600m，-500m］，終端時刻的vx和vy可由式（3）和式（5）的繞飛條件來確定，再設(shè)vz＝vy，則終端的位置和速度完全確定。

在Matlab環(huán)境下，高斯偽譜法建模的開源軟件包GPOPS將最優(yōu)控制問題轉(zhuǎn)化為NLP問題，最終轉(zhuǎn)化出的NLP再用基于SQP的SNOPT［13］求解器來完成求解。在SNOPT中迭代575次，最終求解得到伴星釋放過程持續(xù)時間tf＝1 179.1s。圖1～圖3分別為釋放過程中的控制量、相對位置、相對速度隨時間的變化曲線。

從圖1可以看出，整個變軌過程中除y方向控制量有一個較小尖峰之外，基本上是一個平滑的變化過程，這種平滑的變化易于推力器的實現(xiàn)。從圖2和圖3可以看出兩個航天器的相對距離、速度從開始釋放時基本為零，平穩(wěn)地過渡到穩(wěn)定橢圓繞飛的要求值，沒有出現(xiàn)較大的波動。

圖4是伴星從釋放到穩(wěn)定的圍繞主星飛行20天的三維軌跡，可見衛(wèi)星經(jīng)過最優(yōu)變軌控制，順利進(jìn)入預(yù)定的橢圓繞飛軌道。由于在變軌結(jié)束時刻的相對位置和相對速度與橢圓繞飛所需的值之間仍有誤差，造成x軸方向仍有常值漂移項，按照所求解出的離散控制量進(jìn)行三次樣條差值，積分相對運動方程所得到tf時刻的相對位置和速度為x＝1 000m，vy＝0.553 41m／s，則由式（3）得到的x方向常值飄移速率為3vx＋6ωy＝-6.275×10-5m／s，隨著時間推進(jìn)，繞飛曲線將逐步偏離預(yù)設(shè)軌道。實際上，即使變軌過程完全沒有誤差，在攝動力作用下，伴星相對運動軌道也將逐步飄移，要完成穩(wěn)定的繞飛，需要不斷地進(jìn)行軌道修正。本文僅考慮釋放軌跡的最優(yōu)控制，對繞飛保持控制問題不做討論。

圖1 控制量的時間歷程Fig.1 Result of control

圖2 相對位置時間歷程Fig.2 Result of relative position

圖3 相對速度的時間歷程Fig.3 Result of relative velocity

圖4 釋放與繞飛三維軌跡Fig.4 3Dfigure of releasing and around-flight orbit

5 結(jié)束語

本文利用高斯偽譜法，以無控的橢圓繞飛構(gòu)型為最終目標(biāo)，求解了小推力伴星的最優(yōu)釋放軌跡，仿真算例表明高斯偽譜法無需對狀態(tài)量和控制量給出較好的猜測即可收斂，根據(jù)其求解出的離散控制量進(jìn)行擬合后得到的連續(xù)控制量能夠讓狀態(tài)變量控制到指定值。

［1］張育林，曾國強(qiáng)，王兆魁，等.分布式衛(wèi)星系統(tǒng)理論及應(yīng)用［M］.北京：科學(xué)出版社，2008：3-11.ZHANG YULIN，ZENG GUOQIANG，WANG ZHAOKUI，et al.Theory and application of distribute satellite system ［M］.Beijing：Science Press，2008：3-11.

［2］顧大可，段廣仁，張卯瑞.有限推力交會的最省燃料軌跡［J］.宇航學(xué)報，2010，31（1）：75-81.GU DAKE，DUAN GUANGREN，ZHANG MAORUI.Fuel-optimal trajectories for finite-thrust rendezvous［J］.Journal of Astronautics，2010，31 （1）：75-81.

［3］師鵬，李寶軍，趙育善.有限推力下的航天器繞飛軌道保持與控制［J］.北京航空航天大學(xué)學(xué)報.2007，33（7）：757-760.SHI PENG，LI BAOJUN，ZHAO YUSHAN.Orbital maintenance and control of spacecraft fly-around with finite-thrust［J］.Journal of Beijing University of Aeronautics and Astronautics，2007，33 （7）：757-760.

［4］吳文華，曹喜濱，張世杰，等.編隊衛(wèi)星相對軌道與姿態(tài)一體化耦合控制［J］.南京航空航天大學(xué)學(xué)報，2010，42 （1）：13-20.WU WENHUA，CAO XIBIN，ZHANG SHIJIE，et al.Relative orbit and attitude integrated coupled control for formation satellite［J］.Journal of Nanjing University of Aeronautics ＆ Astronautics.2010，42 （1）：13-20.

［5］周文雅，楊滌，李順利.利用高斯偽譜法求解具有最大橫程的再入軌跡［J］.系統(tǒng)工程與電子技術(shù)，2010，32 （5）：1039-1042.ZHOU WENYA，YANG DI，LI SHUNLI.Solution of reentry trajectory with maximum cross range by using Gauss pseudospectral method［J］.Systems Engineering and Electronics，2010，32 （5）：1039-1042.

［6］尚海濱，崔平遠(yuǎn)，徐瑞，等.基于高斯偽光譜法的星際小推力轉(zhuǎn)移軌道快速優(yōu)化［J］.宇航學(xué)報，2010，31 （4）：1005-1011.SHANG HAIBIN，CUI PINGYUAN，XU RUI，et al.Fast optimization of interplanetary low-thrust transfer trajectory based on Gauss pseudospectral algorithm ［J］.Journal of Astronautics，2010，31 （4）：1005-1011.

［7］DAVID BENSON.A Gauss pseudospectral transcription for optimal control［D］.Cambridge：Department of Aeronautics and Astronautics，MIT，2004.

［8］HUNTINGTON GEOFFREY TODD.Advancement and analysis of a Gauss pseudospectral transcription for optimal control problems［D］.Cambridge：Department of Aeronautics and Astronautics，MIT，2007

［9］雍恩米，唐國金，陳磊.基于Gauss偽譜法的高超聲速飛行器再入軌跡快速優(yōu)化［J］.宇航學(xué)報，2008，29 （6）：1766-1772.YONG ENMI，TANG GUOJIN，CHEN LEI.Rapid trajectory optimization for hypersonic reentry vehicle via Gauss pseudospectral method ［J］.Journal of Astronautics，2008，29 （6）：1766-1772.

［10］彭祺擘，李海陽，沈紅新.基于高斯偽譜法的月球定點著陸軌道快速優(yōu)化設(shè)計［J］.宇航學(xué)報，2010，31 （4）：1012-1016.PENG QIBO，LI HAIYANG，SHEN HONGXIN.Rapid lunar exact-landing trajectory optimization via Gauss pseudospectral method［J］.Journal of Astronautics，2010，31 （4）：1012-1016.

［11］劉暾，趙駿.空間飛行器動力學(xué)［M］.哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2003：83-89.LIU TUN，ZHAO JUN.Spacecraft dynamics［M］.Harbin：Press of Harbin Institute of Technology，2003：83-89.

［12］肖業(yè)倫.航空航天器運動的建模［M］.北京：北京航空航天大學(xué)出版社，2003：118-122.XIAO YELUN.Model of aircraft and spacecraft dynamics［M］.Beijing：Beihang University Press，2003：118-122.

［13］GILL PHILIP E，MURRAY WALTER，SAUNDERS MICHAEL A.SNOPT：an SQP algorithm for large-scale constrained optimization ［J］.SIAM Journal on Optimization，2002，12 （4）：979-1006.