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

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

0 引 言

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

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

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

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

1 模型建立

1.1 飛行器運動學模型

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

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

1.2 高斯偽譜法

高斯偽譜法[1～6]本質(zhì)上是在一系列配點處將狀態(tài)變量和控制變量進行離散，并以離散點構(gòu)造全局拉格朗日插值多項式來近似系統(tǒng)動力學方程，從而將連續(xù)的最優(yōu)控制問題轉(zhuǎn)換為非線性規(guī)劃問題。

a）時域變換。

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

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

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

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

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

c）狀態(tài)方程轉(zhuǎn)換。

狀態(tài)變量x()τ對時間τ進行求導，將動力學方程轉(zhuǎn)換為代數(shù)方程，即：

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

d）目標函數(shù)。

對于燃料最優(yōu)問題，存在目標函數(shù)：

式中ωk為高斯權(quán)重系數(shù)，k= 1 ,2,…,N。

e）轉(zhuǎn)換為NLP問題標準模式。

通過上述過程，將連續(xù)泛函最優(yōu)控制問題轉(zhuǎn)換為離散形式的NLP問題：

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

1.3 燃料最優(yōu)間接法

燃料最優(yōu)性能指標：

根據(jù)最優(yōu)控制理論，引入哈密頓函數(shù)：

式中rλ，vλ，mλ均為協(xié)態(tài)變量；μ為引力系數(shù)。

由極小值原理，推出：

帶入哈密頓函數(shù)，推出推力最優(yōu)控制函數(shù)：

2 軌道優(yōu)化方案設(shè)計

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

2.1 飛行任務(wù)規(guī)劃

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

2.2 模 型

以軌道轉(zhuǎn)移為例，目標需要在T內(nèi)完成軌道轉(zhuǎn)移以及姿態(tài)調(diào)整，假設(shè)t1∈[a,b]，飛行器完成軌道轉(zhuǎn)移，到達指定的目標。則終端約束條件為

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

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

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

性能指標取燃料最優(yōu)，有：

式中mf為飛行器燃料消耗質(zhì)量。

2.3 優(yōu)化模型

2.3.1 初值猜測

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

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

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

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

圖1 優(yōu)化過程 Fig.1 Optimization Procedure

2.3.2 有限差分法

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

考慮本文非線性方程組：

初始條件為

終端條件為

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

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

2.4 優(yōu)化方案

優(yōu)化方案如圖2所示。

圖2 優(yōu)化方案 Fig.2 Optimize Process

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

3 仿真校驗

3.1 場景設(shè)定及參數(shù)設(shè)置

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

飛行器參數(shù)見表1。

表1 飛行器參數(shù)Tab.1 Aircraft Parameters

3.2 仿真結(jié)果

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

續(xù)圖3

3.3 仿真驗證

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

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

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

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

續(xù)圖4

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

4 結(jié)束語

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