邵 楠，閆曉東

(西北工業(yè)大學(xué)航天學(xué)院，西安 710072)

0 引 言

火箭從高空返回并垂直定點(diǎn)著陸回收是火箭重復(fù)使用的一種重要方式。對于火箭垂直回收任務(wù)而言，回收過程不僅要進(jìn)行減速并精確著陸，還要滿足返回過程中各種過程約束，此外還要使得燃料消耗最小。由于火箭回收的初始高度比較高，一般而言，整個(gè)返回彈道可以分為三段：動力減速段、氣動減速段和動力著陸段。動力減速段為使火箭動壓降低至滿足柵格舵工作條件的狀態(tài)；氣動減速段對射程起到調(diào)制作用，并盡可能地利用氣動力降低終端位置誤差；動力著陸段需要滿足位置、速度、姿態(tài)等終端約束實(shí)現(xiàn)定點(diǎn)垂直著陸。

由于火箭回收制導(dǎo)任務(wù)的復(fù)雜性，滿足多約束條件并具有快速收斂特性的制導(dǎo)算法一直是眾多學(xué)者研究的方向。文獻(xiàn)[1-3]提出了一種凸規(guī)劃算法，用于求解火星精確著陸相關(guān)的最小燃料動力下降制導(dǎo)問題。他們提出“無損凸化”的概念，使得非凸控制約束的軌跡優(yōu)化問題轉(zhuǎn)化為一個(gè)有限維二階錐規(guī)劃問題，并在該問題的基礎(chǔ)上進(jìn)一步引入推力指向約束，使改進(jìn)的動力降落制導(dǎo)算法對推力約束和推力指向約束都產(chǎn)生了無損凸化。該方法忽略了氣動力的作用，通過線性搜索步驟確定終端時(shí)刻，無需迭代即可算出最優(yōu)解[4-6]。然而其固定的終端時(shí)刻，無法保證開機(jī)-終端時(shí)刻組合的最優(yōu)性。文獻(xiàn)[7-8]進(jìn)一步提出了一種以燃料最優(yōu)為指標(biāo)的動力著陸問題的逐次凸化算法并給出了逐次凸化的證明。在該算法中，引入了氣動阻力和包括自由終端時(shí)間在內(nèi)的各種非凸約束，通過逐次凸化、逐次線性化雖然增大了計(jì)算量，但是能夠解決更復(fù)雜的約束情況[9-11]。王勁博等[10]針對火箭動力定點(diǎn)垂直著陸提出一種高精度快速軌跡優(yōu)化算法，算法將凸化技術(shù)與偽譜離散方法有機(jī)結(jié)合，將非凸、非線性優(yōu)化問題轉(zhuǎn)化為凸優(yōu)化問題，進(jìn)而充分利用凸優(yōu)化的求解快速性、收斂確定性以及偽譜法離散精度高的特點(diǎn)，實(shí)現(xiàn)了考慮阻力的兩階段軌跡優(yōu)化。從優(yōu)化結(jié)果看，雖然兩階段最優(yōu)規(guī)劃方法比單獨(dú)的動力下降段最優(yōu)規(guī)劃可以節(jié)省更多燃料，但是當(dāng)前多階段最優(yōu)軌跡優(yōu)化方法僅考慮了阻力的影響。事實(shí)上，由于火箭倒飛時(shí)底阻較大，在氣動減速段存在一定的配平攻角，火箭所受的升力也十分顯著，特別是當(dāng)攻角可由柵格舵調(diào)節(jié)時(shí)，升力也可控。這種情況下，過程約束和優(yōu)化結(jié)果產(chǎn)生很大的不同。

基于此，本文考慮完整的氣動力，將氣動力、推力以及發(fā)動機(jī)點(diǎn)火和關(guān)機(jī)時(shí)刻均作為控制變量，基于無損凸優(yōu)化方法研究了多階段垂直返回著陸軌跡優(yōu)化方法，以最大程度利用現(xiàn)有火箭氣動減速能力，減小燃料消耗，并滿足各種過程約束和終端約束?？紤]到氣動減速段飛行時(shí)間久、飛行距離長，且氣動減速段過程中升阻力控制量變化平緩，使用Legendre-Gauss-Radau偽譜法[12-15]，對氣動減速段進(jìn)行離散，減少計(jì)算量同時(shí)提高運(yùn)算精度。同時(shí)，動力下降段由于推力突變的非光滑性，在無損凸化的基礎(chǔ)上采用等距離散的方法構(gòu)建了凸優(yōu)化問題。最后通過對多階段離散最優(yōu)控制問題進(jìn)行無損凸化，將原問題轉(zhuǎn)換為二階錐規(guī)劃問題，從而實(shí)現(xiàn)了火箭多階段定點(diǎn)垂直著陸最優(yōu)軌跡的快速規(guī)劃。

1 模型建立

建立著陸表面固定參考系。由于著陸段飛行距離較短,假設(shè)在整個(gè)著陸過程中地表為一個(gè)平面。

建立地面參考坐標(biāo)系：坐標(biāo)系原點(diǎn)O位于著陸點(diǎn)。OX軸位于水平面，指向火箭初始位置方向?yàn)檎?，OZ軸垂直于水平面，向上為正，OY軸與其他兩軸構(gòu)成右手直角坐標(biāo)系。

在地面參考坐標(biāo)系下，火箭的三自由度動力學(xué)方程為：

(1)

(2)

(3)

其中，r(t),v(t)和m(t)分別為火箭的位置矢量、速度矢量和質(zhì)量。火箭發(fā)動機(jī)推力矢量T(t)方向與火箭體軸一致，Y(t)為火箭升力，D(t)為火箭阻力，定義火箭速度矢量與縱軸之間的夾角為總攻角η，升力和阻力為：

(4)

(5)

總升力方向與速度方向垂直：

YT(t)·v(t)=0

(6)

其中，SD為參考面積，大氣密度使用標(biāo)準(zhǔn)大氣模型，通過特征點(diǎn)選取

(7)

(8)

由于再入返回過程中箭體法向過載不能太大，一般總攻角要約束在較小的范圍內(nèi)，因此有：

|η|≤ηmax

(9)

(10)

考慮到火箭返回段射程較小，因此重力加速度g可表示為：

(11)

g為重力加速度大?。?/p>

(12)

式中：Re為地球平均半徑。

2 多階段離散優(yōu)化模型

2.1 氣動減速段

偽譜離散方法布置N個(gè)離散點(diǎn)，可以獲得最高2N+1次代數(shù)精度。相對其他離散化方法具有離散精度高及收斂速度較快的優(yōu)勢。氣動減速段飛行時(shí)間久、飛行距離長，且氣動減速段升阻力控制量連續(xù)變化，因此氣動減速段更適宜采用偽譜離散方法。使用Legendre-Gauss-Radau偽譜法，為氣動減速段布置N1個(gè)離散點(diǎn)，對氣動減速段進(jìn)行離散：

(13)

L為Radau微分矩陣，k=1,…,N1-1；f(x,u)為式(1)～(3)中微分方程的右端函數(shù);τk為[-1,1)區(qū)間內(nèi)的配點(diǎn);x=[rT,vT,m]T為滑行段的狀態(tài)變量。由于總攻角可調(diào)節(jié)，升力也可調(diào)節(jié)，可將升力作為控制變量，因此，控制為u=[TT,YT]T。

由于Radau配點(diǎn)不包括最后一個(gè)端點(diǎn)τf=1，因此還需滿足如下約束：

(14)

采用偽譜方法離散后的動力學(xué)方程為：

(15)

氣動減速段火箭發(fā)動機(jī)不開機(jī)，因此推力T[k]=0，火箭總升力方向滿足與速度垂直約束，升力阻力滿足于總攻角的約束，離散的控制量約束為：

(16)

邊界條件滿足初始時(shí)刻狀態(tài)約束：

(17)

2.2 動力著陸段

動力著陸段推力在最大和最小推力間變化，可能出現(xiàn)不連續(xù)情況，使用Legendre-Gauss-Radau離散會導(dǎo)致計(jì)算精度下降，因此采用均勻離散方法，設(shè)離散布點(diǎn)數(shù)為N2，則該段終端時(shí)間為：

(18)

(19)

其中，k=N1,N1+1,…,Nf，動力學(xué)方程狀態(tài)的初始和終端約束為：

(20)

其中，mdry為火箭干重，rf和vf分別為要求的終端位置和終端速度。除狀態(tài)約束外，還需要對動力下降段的推力大小和推力指向進(jìn)行約束：

(21)

(22)

(23)

2.3 最優(yōu)控制問題模型

綜上，以燃料消耗最優(yōu)為性能指標(biāo)的火箭自主返回垂直著陸軌跡優(yōu)化問題(問題1)為：

maxJ=m[Nf]

s.t. 式(13)～(17)，式(19)～式(23)

(24)

問題1為燃料消耗最優(yōu)火箭自主返回垂直著陸軌跡優(yōu)化問題的原始問題，由于分母的質(zhì)量項(xiàng)、氣動力非線性以及自由時(shí)間變量的引入導(dǎo)致問題1是非凸的。

3 問題的凸化

由于問題1的約束是非凸的，不能應(yīng)用凸優(yōu)化技術(shù)來求解，因此需要對原問題的約束進(jìn)行凸化處理。針對推力矢量約束的非凸性，本節(jié)通過無損凸化方法進(jìn)行凸化，非線性動力學(xué)方程和其它非線性約束通過逐次線性化的方法轉(zhuǎn)化為線性約束，最終原問題被轉(zhuǎn)化為二階錐規(guī)劃問題。

3.1 無損凸化

針對原非凸問題提出對應(yīng)的松弛凸優(yōu)化問題，通過增加問題的維度,將原非凸可行域適當(dāng)松弛放大，構(gòu)成高維的凸可行域，并通過極大值原理保證松弛問題的最優(yōu)解仍在原問題的可行域內(nèi)，凸化原理如圖2所示。

引入松弛變量Γ(t)對原控制約束進(jìn)行松弛，將原控制變量約束做如下變換：

(25)

Tmin≤Γ(t)≤Tmax

(26)

3.2 逐次凸化

在動力學(xué)方程和其它非線性約束中，非線性主要來源于氣動升力/阻力、分母項(xiàng)中包含的質(zhì)量項(xiàng)以及狀態(tài)與發(fā)動機(jī)開機(jī)時(shí)間和發(fā)動機(jī)關(guān)機(jī)時(shí)間的乘積。這些非線性特性通過逐次線性化的方法予以線性化。

1)動力學(xué)方程約束線性化

在動力學(xué)方程約束(13)中，矩陣L是常數(shù)矩陣，非線性主要來源于包含終端時(shí)間的動力學(xué)方程：

(27)

對式(27)在第i次迭代結(jié)果處進(jìn)行一階泰勒級數(shù)展開，取其線性項(xiàng)得：

(28)

(29)

式(28)中，在離散點(diǎn)處的常數(shù)矩陣為：

將式(28)代入式(13)，得到線性化的動力學(xué)方程約束為：

(30)

對于等距離散的動力學(xué)方程(19)，采用相似的離散化方法，此處不再贅述。

由于引入升力控制量，式(16)增加了一個(gè)非線性約束：

YT[k]·v[k]=0

(31)

采用一階泰勒級數(shù)展開線性化為：

YT[k]·v[k]=Yi[k]·vi[k]+YT[k]·

vi[k]+YiT[k]·v[k]

(32)

2)信賴域約束

由于線性化是在上次迭代結(jié)果上進(jìn)行的，因此若兩次迭代結(jié)果相差比較大時(shí)，可能使問題無解。為了減少這種風(fēng)險(xiǎn)，希望確保線性化所涉及的變量不會與前一個(gè)迭代中獲得的值發(fā)生顯著偏離，因此需要引入信賴域約束來保證線性化的可行性。

將第i+1次迭代變量與第i次迭代變量結(jié)果之間的差約束到一個(gè)范圍內(nèi)有：

|xi+1-xi|≤ex

(33)

|ui+1-ui|≤eu

(34)

其中，ex和eu分別表示容許的偏差值，可提前設(shè)定為固定常數(shù)，也可以根據(jù)線性化偏差不斷更新。

3)松弛

在逐次迭代的初期，特別是當(dāng)選取的初始彈道與最優(yōu)解相差較大時(shí)，可能由于線性化造成“偽不可行”現(xiàn)象，即原問題存在可行解，線性化后退化為不可行問題。在這種情況下，不可行性阻礙了迭代過程，阻礙了收斂。一種解決方法是對式(30)添加松弛變量，并把松弛變量加入到性能指標(biāo)中最小化，松弛約束為：

(35)

式中:Ω(τk)為松弛變量。

由于引入升力控制量，式(16)經(jīng)過一階泰勒級數(shù)展開線性化為(32)，該式過于嚴(yán)格，在線性化迭代求解中需要進(jìn)行松弛以獲得更好的迭代結(jié)果：

YT[k]·v[k]<=y(k)

(36)

(37)

式中：kΩ,k是系數(shù)向量，調(diào)節(jié)該系數(shù)向量使得松弛變量趨于0。

4 無損凸化的最優(yōu)控制問題及其求解

4.1 無損凸化以后的優(yōu)化問題

問題1經(jīng)過無損凸化、線性化和松弛，得到由不等式約束、等式約束和二階錐約束構(gòu)成的問題2：

s.t.

(38)

原始問題1經(jīng)過對推力凸化處理后，得到了新的問題2?？梢宰C明，該凸化過程是一種無損凸化，因此問題2與問題1是等價(jià)的[1]。

文獻(xiàn)[7]針對非線性動力學(xué)、非線性過程約束以及二階錐約束的序列凸化算法進(jìn)行了一階條件下的收斂性證明。本文中經(jīng)凸化處理后問題2是一個(gè)典型的序列二階錐規(guī)劃問題，根據(jù)凸優(yōu)化理論[16]，問題2的解序列至少存在一個(gè)聚點(diǎn)，該聚點(diǎn)為原問題的駐點(diǎn)。因此，在一階條件下，證明序列凸優(yōu)化過程是收斂的，并且當(dāng)信賴域約束動態(tài)調(diào)整時(shí)，具有超線性收斂特性。

4.2 以射程為性能指標(biāo)的凸優(yōu)化問題

為了比較升力對射程的影響，以射程最大和最小為性能指標(biāo)，構(gòu)建問題3，其邊界約束與過程約束與問題2有所不同，修改如下。

修改性能指標(biāo)(38)，分別以最大射程和最小射程為性能指標(biāo)：

最小射程：

最大射程：

終端位置約束式(20)中，由于航向位置作為性能指標(biāo)，因此去除航向位置約束：

(39)

由于終端位置不再是原點(diǎn)，修改下滑斜率約束式(23)，使下滑斜率約束隨終端位置變化：

(40)

雖然性能指標(biāo)和約束有所不同，但問題3依舊是一個(gè)典型的序列二階錐規(guī)劃問題。

4.3 問題求解過程

問題2和問題3是一個(gè)序列二階錐規(guī)劃問題，需要通過序列凸化算法進(jìn)行逐次迭代求解。在迭代求解時(shí)，首先需要提供一組狀態(tài)初始剖面，一組良好的初始剖面有助于迭代過程快速收斂到最優(yōu)解。

基于該初值，針對問題2進(jìn)行迭代求解，直到兩次優(yōu)化結(jié)果狀態(tài)量的差滿足收斂條件為止。圖3為迭代求解流程圖。

5 仿真校驗(yàn)

仿真算例分別計(jì)算了無升力軌跡優(yōu)化和升力可控軌跡優(yōu)化，對比說明了加入升力后對最優(yōu)燃料性能指標(biāo)的影響。選擇某型火箭作為仿真模型，任務(wù)初值與終端約束如表1所示?；鸺w行參數(shù)與氣動參數(shù)如表2所示。

參數(shù)值v0/(m·s-1)[-500 -10 -70]Tvf/(m·s-1)[0 0 0]Tr0/km[41.5 0.2 47.2]Trf/km[0 0 0]T

表2 火箭參數(shù)Table 2 Parameters of rocket

本文所有數(shù)值計(jì)算在MATLAB 2017b環(huán)境下編制及運(yùn)行，采用CVX優(yōu)化工具箱和MOSEK 8.0.0.60作為算例的凸優(yōu)化求解工具。

5.1 無升力算例

該算例不考慮升力作用，在氣動減速段僅有氣動阻力作用。氣動減速段和動力著陸段各為30個(gè)離散點(diǎn)?；鸺w行時(shí)間作為優(yōu)化變量，優(yōu)化后火箭總飛行時(shí)間為135.7 s，其中氣動減速飛行時(shí)間為123.7 s，動力著陸段飛行時(shí)間為12.0 s。關(guān)機(jī)后火箭質(zhì)量剩余14.89 t，總?cè)剂舷臑?.33 t。

圖4為逐次凸化迭代后氣動減速段的三維軌跡，表明火箭成功在預(yù)定落點(diǎn)降落，圖中所示圓錐為下滑斜率約束面。從圖5、圖6可以看出，位置和速度均收斂到預(yù)定終端狀態(tài)。經(jīng)過空氣阻力減速，發(fā)動機(jī)開機(jī)點(diǎn)的坐標(biāo)為[114.7， -140.8， 1089] m，初始速度為[-19.38， -0.33， -180.3] m/s。

圖7為火箭發(fā)動機(jī)推力大小和三向推力曲線。從仿真結(jié)果可以看出，氣動減速段推力為0，氣動減速段后，火箭先以最小推力工作，然后切換到最大推力狀態(tài)。

5.2 氣動升力作為控制量仿真算例

本算例考慮氣動升力，且氣動升力大小能夠由總攻角進(jìn)行調(diào)節(jié)，總攻角約束為：|η|≤15°。優(yōu)化后火箭總飛行時(shí)間為144.0 s，其中氣動減速段飛行時(shí)間為135.1 s，動力著陸段飛行時(shí)間為8.9 s。關(guān)機(jī)時(shí)刻火箭剩余質(zhì)量為15.09 t，燃料消耗為1.13 t。相比無升力軌跡優(yōu)化燃料節(jié)省0.2 t。

圖8、圖9表明速度和位置滿足終端約束。動力著陸段初始的坐標(biāo)為[222.1, -105.3, 534.8] m，初始速度為[-54.03, 25.65， -125.9] m/s。

從仿真結(jié)果可以看出，有升力情況下發(fā)動機(jī)的開機(jī)高度明顯降低，速度矢量的方向更趨近于減少位置誤差的方向，從而導(dǎo)致發(fā)動機(jī)工作時(shí)間縮短，更加節(jié)省燃料。

圖10為火箭總攻角η隨時(shí)間變化曲線，最大攻角為15°，可以看出攻角滿足約束條件。

5.3 升力對射程調(diào)制對比

通過求解問題3，可以獲得回收過程的最小射程和最大射程。對于無升力情況，最小射程為41.232 km，最大射程為41.867 km，落點(diǎn)射程調(diào)節(jié)范圍為0.623 km。由于無法對升力進(jìn)行調(diào)節(jié)，射程調(diào)制主要由推力調(diào)節(jié)來實(shí)現(xiàn)，因此圖11僅給出動力下降段的縱向軌跡曲線。不難發(fā)現(xiàn)，由于存在最大燃料、推力指向、推力大小以及下降斜率約束，射程調(diào)節(jié)范圍較小。

圖12為以升力可控情況的最小、最大射程軌跡?；鸺w行最小射程為36.97 km，最大射程為47.171 km，落點(diǎn)射程調(diào)節(jié)范圍為10.201 km。相較于無升力控制算例，可大幅提升可達(dá)域范圍。

6 結(jié) 論

本文針對火箭垂直回收提出一種將氣動減速段與動力著陸段統(tǒng)一規(guī)劃的多階段燃料最優(yōu)軌跡規(guī)劃算法。該算法通過將原問題進(jìn)行無損凸化，并經(jīng)過線性化和離散化，得到序列SOCP問題，可快速迭代求解。通過仿真校驗(yàn)，可得出如下結(jié)論：

1)將氣動力作為控制變量時(shí)，仍可對原問題進(jìn)行無損凸化。

2)將發(fā)動機(jī)開、關(guān)機(jī)時(shí)間作為優(yōu)化變量有助于獲得全局燃料最優(yōu)解。

3)在氣動減速段調(diào)節(jié)升力，可以獲得更大的射程調(diào)節(jié)范圍，并可以進(jìn)一步節(jié)省動力著陸段的燃料。