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

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

0 引 言

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

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

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

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

1 模型建立

1.1 飛行器運(yùn)動(dòng)學(xué)模型

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

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

1.2 高斯偽譜法

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

a）時(shí)域變換。

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

式中t0為求解初始時(shí)刻；tf為求解末時(shí)刻。

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

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

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

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

狀態(tài)變量x()τ對(duì)時(shí)間τ進(jìn)行求導(dǎo)，將動(dòng)力學(xué)方程轉(zhuǎn)換為代數(shù)方程，即：

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

d）目標(biāo)函數(shù)。

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

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

e）轉(zhuǎn)換為NLP問題標(biāo)準(zhǔn)模式。

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

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

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

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

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

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

由極小值原理，推出：

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

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

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

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

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

2.2 模 型

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

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

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

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

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

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

2.3 優(yōu)化模型

2.3.1 初值猜測

簡化飛行器運(yùn)動(dòng)模型[11]：

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

對(duì)于常值推力發(fā)動(dòng)機(jī)，最短時(shí)間飛行問題可以等價(jià)于燃料最優(yōu)，因此，最優(yōu)控制性能指標(biāo)簡化為

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

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

2.3.2 有限差分法

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

考慮本文非線性方程組：

初始條件為

終端條件為

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

將區(qū)間離散化，在節(jié)點(diǎn)應(yīng)用差分公式，得到新的代數(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)成的邊值問題進(jìn)行求解，通過有限差分法繼續(xù)優(yōu)化。當(dāng)滿足精度要求時(shí)，結(jié)束優(yōu)化，輸出最優(yōu)控制推力、姿態(tài)指令。

3 仿真校驗(yàn)

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

假設(shè)初始時(shí)刻有2個(gè)飛行器A、B，A和B相距5 km（方向隨機(jī)設(shè)定），A、B逐漸拉開距離。C飛行器在初始時(shí)刻于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ā)動(dòng)機(jī)開關(guān)控制量隨時(shí)間變化的關(guān)系如圖3a所示；飛行器質(zhì)量變化同發(fā)動(dòng)機(jī)開關(guān)控制量隨時(shí)間變化的關(guān)系如圖3b所示，發(fā)動(dòng)機(jī)開機(jī)時(shí)間燃料持續(xù)消耗，飛行器質(zhì)量呈線性下降，發(fā)動(dòng)機(jī)關(guān)機(jī)時(shí)間，無燃料消耗，飛行器質(zhì)量不變；飛行器的優(yōu)化軌跡如圖3c所示，飛行器在起始階段發(fā)動(dòng)機(jī)工作向預(yù)定目標(biāo)軌道轉(zhuǎn)移，中間發(fā)動(dòng)機(jī)關(guān)機(jī)自由飛行，到預(yù)定目標(biāo)附近發(fā)動(dòng)機(jī)再次工作實(shí)現(xiàn)伴飛（圖3中實(shí)線為控制變量，虛線為發(fā)動(dòng)機(jī)控制量）。

續(xù)圖3

3.3 仿真驗(yàn)證

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

表2 方案對(duì)比Tab.2 Optimization Scheme Comparison

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

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

續(xù)圖4

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

4 結(jié)束語

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