基于高斯偽譜法和有限差分結合的飛行器 燃料最優軌道轉移優化

周珺儀，劉佳琪，焦勝海，段紅亮，陳曉光

（1.北京航天長征飛行器研究所，北京，100076；2.試驗物理與計算數學國家級重點實驗室，北京，100076）

0 引 言

隨著戰場態勢的日趨復雜，空間任務逐漸多樣化，要求飛行器能夠在短時間內快速完成多項任務。因此需要運載器在攜帶有限燃料的條件下，具備快速、多次機動至不同軌道完成任務的能力，且每次任務消耗的燃料越少，運載器能夠執行的任務越多，其作戰能力就越強。為了提高燃料的有效利用率，飛行器必須根據任務目標，自主規劃燃料最優變軌方案。

燃料最優變軌[1～3]作為一個連續過程的優化問題，其本質是泛函求解極值的問題。常用的方法有直接法、間接法和混合法。直接法是對參數進行離散化，將問題轉化為參數優化求解；間接法是基于極小值原理，將問題轉化為兩點邊值求解；混合法則是將問題轉化為具有約束的參數求解，用非線性規劃方法進行求解。

偽譜法作為經典的直接優化法，已經廣泛應用于深空探測衛星軌道優化[3]、臨近空間高超聲速飛行器制導優化[4]、日-火Halo轉移軌道快速優化設計[5]、二級助推火箭多階段軌跡優化[6]等領域。間接法[7～10]計算精度高，包括用于最優月球軟著陸軌道的隱式打靶法、軌道轉移的多重打靶法、全電推進衛星軌道優化的同倫解法等。

偽譜法收斂速度快，但控制曲線存在突變，不能直接應用于工程；有限差分法作為間接法的一種，結果精度高，但對初值猜測敏感，初值猜測誤差過大會導致發散。為了能夠快速得到任意空間任務需求下燃料消耗最少的精確軌跡，本文將高斯偽譜法和有限差分法結合，進行小推力變軌優化設計，提出考慮任務規劃下常值推力短時間快速變軌最優燃料消耗的精確解，并給出相應的仿真結果。

1 模型建立

1.1 飛行器運動學模型

假設飛行器在中段飛行過程處于瞬時平衡狀態；忽略除發動機推力之外的攝動，飛行器控制系統處于理想工作情況，不存在延時，則小推力飛行器動力學模型[3]為

式中m為飛行器總質量；T為發動機推力大小；g0為重力加速度；Isp為比沖；r為飛行器的位置矢量；v為飛行器的速度矢量；α為推力方向的單位矢量。

1.2 高斯偽譜法

高斯偽譜法[1～6]本質上是在一系列配點處將狀態變量和控制變量進行離散，并以離散點構造全局拉格朗日插值多項式來近似系統動力學方程，從而將連續的最優控制問題轉換為非線性規劃問題。

a）時域變換。

偽譜法求解時域為[-1，1]，因此要將實際求解時間[t0,tf]進行變換，引入轉換變量τ∈[-1，1]，即：

式中t0為求解初始時刻；tf為求解末時刻。

b）狀態變量x()τ、控制變量u()τ離散化。

通過對選取的高斯點構造插值多項式對連續狀態變量x()τ、控制變量u()τ進行逼近，即：

式中pi(τ)為插值基函數，插值基函數及其求導結果為：

c）狀態方程轉換。

狀態變量x()τ對時間τ進行求導，將動力學方程轉換為代數方程，即：

式中Dni為離散矩陣元素；f(Xn,Un,τn;t0,tf)為運動學方程。

d）目標函數。

對于燃料最優問題，存在目標函數：

式中ωk為高斯權重系數，k= 1 ,2,…,N。

e）轉換為NLP問題標準模式。

通過上述過程，將連續泛函最優控制問題轉換為離散形式的NLP問題：

式中J為性能指標函數；Φ為邊界條件；C為不等式約束條件。

1.3 燃料最優間接法

燃料最優性能指標：

根據最優控制理論，引入哈密頓函數：

式中rλ，vλ，mλ均為協態變量；μ為引力系數。

由極小值原理，推出：

帶入哈密頓函數，推出推力最優控制函數：

2 軌道優化方案設計

本文以飛行器多次自主機動進行軌道位置轉移為研究背景，仿真選取 2個機動目標位置、一次實施機動變軌為例進行軌道位置轉移分析優化。

2.1 飛行任務規劃

已知飛行器初始位置R0、速度V0及初始姿態，在飛行過程中期望于t1內由目標A位置機動變軌到目標B位置附近伴飛，并保持相對靜止時長t2。

2.2 模 型

以軌道轉移為例，目標需要在T內完成軌道轉移以及姿態調整，假設t1∈[a,b]，飛行器完成軌道轉移，到達指定的目標。則終端約束條件為

式中v(b)，x(b)分別為飛行器轉移完成后在b時刻的速度和位置；VB(b)，XB(b)分別為b時刻慣性飛行器的速度和位置；ε為誤差允許量。

控制量U=[αx,αy,αz,u]約束條件為

式中xα，yα，zα分別為推力沿x，y，z3個方向的分量；u為發動機控制量。

性能指標取燃料最優，有：

式中mf為飛行器燃料消耗質量。

2.3 優化模型

2.3.1 初值猜測

簡化飛行器運動模型[11]：

引力加速度采用平方反比引力場模型時，對于飛行器位置變換，引力加速度變化不大，因此采用平均引力假設來近似引力計算。

對于常值推力發動機，最短時間飛行問題可以等價于燃料最優，因此，最優控制性能指標簡化為

因此，對于偽譜法求解最優燃料消耗問題時，控制量可以通過基于脈沖推力的軌道求解獲得近似解，加快求解速度，具體過程如圖1所示。

圖1 優化過程 Fig.1 Optimization Procedure

2.3.2 有限差分法

最優控制模型可以簡化為含未知參數的兩點邊值問題求解，針對本文求解14個方程組的邊值問題，選擇有限差分法進行優化。

考慮本文非線性方程組：

初始條件為

終端條件為

將區間t∈[t0,tf]等分為N個子區間，y(t)在ti處Taylor展開取t=ti+1=t+ih，忽略二階以上部分，得一階導數的前向差分近似：

將區間離散化，在節點應用差分公式，得到新的代數方程組，以及方程組的近似Jocabi矩陣。

2.4 優化方案

優化方案如圖2所示。

圖2 優化方案 Fig.2 Optimize Process

如圖2所示，首先，根據初始條件按Lambert問題求解近似控制變量初值，并將求得的變量初值帶入Gauss偽譜法求解方程，取誤差小于10-4快速求解協態變量近似值。然后，將求解的協態變量近似初值和初始狀態變量帶入間接法構成的邊值問題進行求解，通過有限差分法繼續優化。當滿足精度要求時，結束優化，輸出最優控制推力、姿態指令。

3 仿真校驗

3.1 場景設定及參數設置

假設初始時刻有2個飛行器A、B，A和B相距5 km（方向隨機設定），A、B逐漸拉開距離。C飛行器在初始時刻于A附近伴飛，要求C飛行器6 s后向B飛行器轉移，40 s內完成軌道轉移在B飛行器附近伴飛。

飛行器參數見表1。

表1 飛行器參數Tab.1 Aircraft Parameters

3.2 仿真結果

C飛行器40 s完成軌道轉移后，飛行器剩余質量830.3 kg。仿真得到推力沿x、y、z三軸的控制變量同發動機開關控制量隨時間變化的關系如圖3a所示；飛行器質量變化同發動機開關控制量隨時間變化的關系如圖3b所示，發動機開機時間燃料持續消耗，飛行器質量呈線性下降，發動機關機時間，無燃料消耗，飛行器質量不變；飛行器的優化軌跡如圖3c所示，飛行器在起始階段發動機工作向預定目標軌道轉移，中間發動機關機自由飛行，到預定目標附近發動機再次工作實現伴飛（圖3中實線為控制變量，虛線為發動機控制量）。

續圖3

3.3 仿真驗證

根據3.2節的仿真實例，采用高斯偽譜法同本文優化方案進行對比，驗證本文方案能夠快速收斂，能夠應用于在線規劃燃料最優轉移軌道。表2為不同優化方案的結果、計算速度對比。通過表2可以發現，采用本文提出的優化方法，改進初值后的計算時長比隨機初值的高斯偽譜法大大縮短，減少了計算量。同時，3種方法的飛行器剩余質量誤差處于允許范圍，驗證了本文方法的準確性。

表2 方案對比Tab.2 Optimization Scheme Comparison

選取3種方法得到沿x軸的控制變量以及發動機開關曲線隨時間的變化關系作為對比，仿真結果如圖4所示（圖4中實線為控制變量，虛線為發動機控制量）。

圖4 x軸控制變量及發動機開關曲線變化關系Fig.4 Relationship between x-axis Control Variables and Engine Switch Curve

續圖4

由圖4可知，本文采用的結合法，控制曲線連續光滑；采用精度較低的高斯偽譜法，雖然能夠快速計算出結果，但是推力控制曲線為變推力曲線，同實際情況不符，誤差較大；采用高斯偽譜法同本文方法結果相差不大，但是曲線存在突變，計算時間過長，不能很好地適應快速變化的戰場態勢。

4 結束語

本文研究了飛行器在空間自主軌道轉移優化的方法，提出了改進初值猜測的高斯偽譜法和有限差分法相結合的方法，分析了常值推力機動變軌燃料最優的軌跡優化問題，基于設定的場景和參數開展了軌道轉移仿真分析。研究結果表明：相比高斯偽譜法，本文所用方法控制過程穩定無誤差，收斂速度提高，計算結果精度更高，可滿足戰場空間任務快速、多樣化的需求。