周 洋，閆 野，黃 煦

（國防科學(xué)技術(shù)大學(xué) 航天科學(xué)與工程學(xué)院，湖南 長沙 410073）

0 引言

解決有限推力航天器最優(yōu)交會問題主要有間接法和直接法［1］。直接法是將連續(xù)的最優(yōu)控制問題離散并參數(shù)化，直接用數(shù)值方法尋優(yōu)，由于計算機(jī)技術(shù)進(jìn)展獲得了很大的發(fā)展［2］。間接法是基于Pontryagin極大值原理將最優(yōu)控制問題轉(zhuǎn)化為一個邊值問題，但協(xié)態(tài)變量的初值難以選取。直接配點(diǎn)法將控制變量和狀態(tài)變量同時離散，也稱作直接配點(diǎn)非線性規(guī)劃（DCNLP）［3］。文獻(xiàn)［4］用非線性規(guī)劃法研究了共面航天器的最優(yōu)交會問題，但其設(shè)計變量多達(dá)159個。文獻(xiàn)［5］先求取變軌問題的雙脈沖最優(yōu)解，再用Gauss偽譜法求解有限推力解。文獻(xiàn)［6］用高階多項式以獲得更高的精度，但增加了非線性規(guī)劃（NLP）問題求解的難度。文獻(xiàn)［7］用五階多項式擬合狀態(tài)量解決了最小時間攔截等簡單問題，具較高的精度。

本文基于改進(jìn)配點(diǎn)法，對有限推力航天器最省燃料交會進(jìn)行了研究。

1 問題描述

航天器交會過程中若目標(biāo)航天器運(yùn)行于圓軌道，且相對距離遠(yuǎn)小于目標(biāo)軌道半徑時，則可由CW方程近似描述兩個航天器的相對運(yùn)動。C-W方程建立在目標(biāo)軌道坐標(biāo)系中，坐標(biāo)原點(diǎn)O位于目標(biāo)質(zhì)心，Ox軸沿目標(biāo)的地心矢徑方向，Oz軸與目標(biāo)軌道角動量方向重合，Oy軸由右手法則確定。C-W方程具體形式為

式中：x，y，z別為追蹤航天器的位置狀態(tài)；aTx，aTy，aTz分別為3個方向的推力加速度；n為目標(biāo)軌道角速度。

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

由于采用有限推力，推力加速度aT滿足約束

DCNLP中最常見是三階Simpson法。具體為：先將系統(tǒng)整個時間過程分為N段，每段的兩個端點(diǎn)為節(jié)點(diǎn)，在兩節(jié)點(diǎn)間用三次Hermite多項式代表狀態(tài)變量隨時間的變化關(guān)系，并假定控制量變化為線性［3］。

每一段區(qū)間為［ti，ti＋1］，記

節(jié)點(diǎn)間狀態(tài)量用三次Hermite多項式表示

則多項式的系數(shù)矩陣

形成由

缺陷向量構(gòu)成等式約束。此處：

2 改進(jìn)配點(diǎn)法

文獻(xiàn)［7］指出，用更高階Gauss-Lobatto多項式表示節(jié)點(diǎn)間的狀態(tài)量時有更高的精度，但會加大求解NLP問題的難度。但如系統(tǒng)方程有

形式，就可用五階多項式擬合狀態(tài)量的變化，同樣假設(shè)控制量的變化為線性，能在不增加約束條件的情況下提高精度。對C-W方程，有

式中：X1，X2分別為位置和速度狀態(tài)，且X1＝［xyz］T，X2＝ ［vxvyvz］T；A，B為系數(shù)矩陣，且

狀態(tài)量X1用五階多項式表示為

則五階多項式的系數(shù)為

與三階直接配點(diǎn)法類似，形成缺陷向量

其中，中間狀態(tài)量X1（tc1），X1（tc2），X2（tc1），X2（tc2）滿足

用SQP算法優(yōu)化時直接將缺陷向量作為約束條件，效率非常低。將約束改寫為

的線性形式后，優(yōu)化效率可明顯提高。此處：

為需優(yōu)化的狀態(tài)量和控制量；

beq＝ ［（X0）TQ1Q1…Q1（Xtf）T］T.此處：Z1＝（X11）T；Z2＝（X21）T；Z3＝（U1）T；Z4＝（X12）T；ZN＋1＝ （X1N＋1）T；Z2N＋1＝ （X2N＋1）T；Z2N＋2＝ （U2N＋1）T；C＝ ［I6×606×3］；CL＝［C1C2C3］3×9；CU＝［C4C5C6］3×9；Q1＝03×4。其中：C1～C6分別為X1a，X2a，Ua，X1b，X2b，Ub的系數(shù)矩陣。

3 仿真與分析

本文算例是初始狀態(tài)和終端狀態(tài)固定的非共面時間自由交會，驗證改進(jìn)配點(diǎn)法的有效性。設(shè)追蹤航天器的最大推力加速度0.8m／s2，目標(biāo)航天器處于軌道高度500km的圓軌道。初始相對狀態(tài)和終端相對狀態(tài)分別為

式中：x0＝1 000m；y0＝2 000m；z0＝1 500m；v0x＝－10m／s；v0y＝－20m／s；v0z＝－15m／s；xtf＝y(tǒng)tf＝ztf＝0m；vtfx＝vtfy＝vtfz＝0m／s。將整個過程分為N＝10段，用改進(jìn)配點(diǎn)法將問題轉(zhuǎn)為NLP問題，再用SQP算法進(jìn)行參數(shù)優(yōu)化。優(yōu)化結(jié)果為：整個交會時間229.82s，消耗速度增量27.26m／s。交會過程中推力加速度如圖1所示，狀態(tài)量如圖2所示。

圖1 推力加速度Fig.1 Thrust acceleration

圖2 狀態(tài)量Fig.2 State

為分析改進(jìn)配點(diǎn)法的優(yōu)化精度，用優(yōu)化的推力加速度對C-W方程進(jìn)行數(shù)值積分，結(jié)果如圖3所示。由圖可知：最終交會的精度很高，位置精度10－3m，速度精度10－7m／s，均較直接配點(diǎn)法高出1個量級。

4 結(jié)束語

本文基于C-W方程，用改進(jìn)配點(diǎn)法求解最省燃料交會問題。給出了最優(yōu)控制問題的數(shù)學(xué)模型，并基于C-W方程推導(dǎo)了五階多項式擬合時的缺陷向量及內(nèi)點(diǎn)的表達(dá)式。然后將缺陷向量的等式約束改寫成線性約束的形式，利于用SQP進(jìn)行參數(shù)優(yōu)化。五階多項式擬合狀態(tài)變量能在不增加約束的前提下提高精度。用SQP算法時將約束寫成線性約束的形式可顯著提高優(yōu)化效率。需指出的是，改進(jìn)配點(diǎn)法適于一類具有特殊形式的方程，而并不局限于CW方程。改進(jìn)配點(diǎn)法能用于絕對軌道轉(zhuǎn)移、橢圓軌道交會等問題。

圖3 狀態(tài)量數(shù)值積分變化曲線Fig.3 State components vs.time via numerical integration

［1］ CONWAY B A.A Survey of methods available for the numerical optimization of continuous dynamic systems［J］.J Optim Theory Appl，2012（152）：271-306.

［2］ 唐國金，羅亞中，雍恩米.航天器軌跡優(yōu)化理論、方法及應(yīng)用［M］.北京：科學(xué)出版社，2012.

［3］ 雍恩米，陳 磊，唐國金.飛行器軌跡優(yōu)化數(shù)值方法綜述［J］.宇航學(xué)報，2009，29（2）：397-406.

［4］ 王 華，唐國金.用非線性規(guī)劃求解有限推力最優(yōu)交會［J］.國防科技大學(xué)學(xué)報，2003，25（5）：9-13.

［5］ 桑 艷，周 進(jìn).基于偽譜方法的有限推力變軌策略研究［J］.飛行力學(xué)，2010，28（2）：71-74.

［6］ HERMAN A L，CONWAY B A.Direct optimization using collocation based on high-order Gauss-Lobatto quadrature rules［J］.J Guid Control Dyna，1996，19：592-599.

［7］ HU G S，ONG C J，TEO C L.An enhanced transcribing scheme for the numerical solution of a class of optimal control problems［J］.Engineering Optimization，2002，34（2）：155-173.