邵 楠,閆曉東
(西北工業大學航天學院,西安 710072)
火箭從高空返回并垂直定點著陸回收是火箭重復使用的一種重要方式。對于火箭垂直回收任務而言,回收過程不僅要進行減速并精確著陸,還要滿足返回過程中各種過程約束,此外還要使得燃料消耗最小。由于火箭回收的初始高度比較高,一般而言,整個返回彈道可以分為三段:動力減速段、氣動減速段和動力著陸段。動力減速段為使火箭動壓降低至滿足柵格舵工作條件的狀態;氣動減速段對射程起到調制作用,并盡可能地利用氣動力降低終端位置誤差;動力著陸段需要滿足位置、速度、姿態等終端約束實現定點垂直著陸。
由于火箭回收制導任務的復雜性,滿足多約束條件并具有快速收斂特性的制導算法一直是眾多學者研究的方向。文獻[1-3]提出了一種凸規劃算法,用于求解火星精確著陸相關的最小燃料動力下降制導問題。他們提出“無損凸化”的概念,使得非凸控制約束的軌跡優化問題轉化為一個有限維二階錐規劃問題,并在該問題的基礎上進一步引入推力指向約束,使改進的動力降落制導算法對推力約束和推力指向約束都產生了無損凸化。該方法忽略了氣動力的作用,通過線性搜索步驟確定終端時刻,無需迭代即可算出最優解[4-6]。然而其固定的終端時刻,無法保證開機-終端時刻組合的最優性。文獻[7-8]進一步提出了一種以燃料最優為指標的動力著陸問題的逐次凸化算法并給出了逐次凸化的證明。在該算法中,引入了氣動阻力和包括自由終端時間在內的各種非凸約束,通過逐次凸化、逐次線性化雖然增大了計算量,但是能夠解決更復雜的約束情況[9-11]。王勁博等[10]針對火箭動力定點垂直著陸提出一種高精度快速軌跡優化算法,算法將凸化技術與偽譜離散方法有機結合,將非凸、非線性優化問題轉化為凸優化問題,進而充分利用凸優化的求解快速性、收斂確定性以及偽譜法離散精度高的特點,實現了考慮阻力的兩階段軌跡優化。從優化結果看,雖然兩階段最優規劃方法比單獨的動力下降段最優規劃可以節省更多燃料,但是當前多階段最優軌跡優化方法僅考慮了阻力的影響。事實上,由于火箭倒飛時底阻較大,在氣動減速段存在一定的配平攻角,火箭所受的升力也十分顯著,特別是當攻角可由柵格舵調節時,升力也可控。這種情況下,過程約束和優化結果產生很大的不同。
基于此,本文考慮完整的氣動力,將氣動力、推力以及發動機點火和關機時刻均作為控制變量,基于無損凸優化方法研究了多階段垂直返回著陸軌跡優化方法,以最大程度利用現有火箭氣動減速能力,減小燃料消耗,并滿足各種過程約束和終端約束??紤]到氣動減速段飛行時間久、飛行距離長,且氣動減速段過程中升阻力控制量變化平緩,使用Legendre-Gauss-Radau偽譜法[12-15],對氣動減速段進行離散,減少計算量同時提高運算精度。同時,動力下降段由于推力突變的非光滑性,在無損凸化的基礎上采用等距離散的方法構建了凸優化問題。最后通過對多階段離散最優控制問題進行無損凸化,將原問題轉換為二階錐規劃問題,從而實現了火箭多階段定點垂直著陸最優軌跡的快速規劃。
建立著陸表面固定參考系。由于著陸段飛行距離較短,假設在整個著陸過程中地表為一個平面。
建立地面參考坐標系:坐標系原點O位于著陸點。OX軸位于水平面,指向火箭初始位置方向為正,OZ軸垂直于水平面,向上為正,OY軸與其他兩軸構成右手直角坐標系。
在地面參考坐標系下,火箭的三自由度動力學方程為:
(1)
(2)
(3)
其中,r(t),v(t)和m(t)分別為火箭的位置矢量、速度矢量和質量?;鸺l動機推力矢量T(t)方向與火箭體軸一致,Y(t)為火箭升力,D(t)為火箭阻力,定義火箭速度矢量與縱軸之間的夾角為總攻角η,升力和阻力為:
(4)
(5)
總升力方向與速度方向垂直:
YT(t)·v(t)=0
(6)
其中,SD為參考面積,大氣密度使用標準大氣模型,通過特征點選取
(7)
(8)
由于再入返回過程中箭體法向過載不能太大,一般總攻角要約束在較小的范圍內,因此有:
|η|≤ηmax
(9)
(10)
考慮到火箭返回段射程較小,因此重力加速度g可表示為:
(11)
g為重力加速度大?。?/p>
(12)
式中:Re為地球平均半徑。
偽譜離散方法布置N個離散點,可以獲得最高2N+1次代數精度。相對其他離散化方法具有離散精度高及收斂速度較快的優勢。氣動減速段飛行時間久、飛行距離長,且氣動減速段升阻力控制量連續變化,因此氣動減速段更適宜采用偽譜離散方法。使用Legendre-Gauss-Radau偽譜法,為氣動減速段布置N1個離散點,對氣動減速段進行離散:
(13)
L為Radau微分矩陣,k=1,…,N1-1;f(x,u)為式(1)~(3)中微分方程的右端函數;τk為[-1,1)區間內的配點;x=[rT,vT,m]T為滑行段的狀態變量。由于總攻角可調節,升力也可調節,可將升力作為控制變量,因此,控制為u=[TT,YT]T。
由于Radau配點不包括最后一個端點τf=1,因此還需滿足如下約束:
(14)
采用偽譜方法離散后的動力學方程為:
(15)
氣動減速段火箭發動機不開機,因此推力T[k]=0,火箭總升力方向滿足與速度垂直約束,升力阻力滿足于總攻角的約束,離散的控制量約束為:
(16)
邊界條件滿足初始時刻狀態約束:
(17)
動力著陸段推力在最大和最小推力間變化,可能出現不連續情況,使用Legendre-Gauss-Radau離散會導致計算精度下降,因此采用均勻離散方法,設離散布點數為N2,則該段終端時間為:
(18)
(19)
其中,k=N1,N1+1,…,Nf,動力學方程狀態的初始和終端約束為:
(20)
其中,mdry為火箭干重,rf和vf分別為要求的終端位置和終端速度。除狀態約束外,還需要對動力下降段的推力大小和推力指向進行約束:
(21)
(22)

(23)


綜上,以燃料消耗最優為性能指標的火箭自主返回垂直著陸軌跡優化問題(問題1)為:
maxJ=m[Nf]
s.t. 式(13)~(17),式(19)~式(23)
(24)
問題1為燃料消耗最優火箭自主返回垂直著陸軌跡優化問題的原始問題,由于分母的質量項、氣動力非線性以及自由時間變量的引入導致問題1是非凸的。
由于問題1的約束是非凸的,不能應用凸優化技術來求解,因此需要對原問題的約束進行凸化處理。針對推力矢量約束的非凸性,本節通過無損凸化方法進行凸化,非線性動力學方程和其它非線性約束通過逐次線性化的方法轉化為線性約束,最終原問題被轉化為二階錐規劃問題。
針對原非凸問題提出對應的松弛凸優化問題,通過增加問題的維度,將原非凸可行域適當松弛放大,構成高維的凸可行域,并通過極大值原理保證松弛問題的最優解仍在原問題的可行域內,凸化原理如圖2所示。
引入松弛變量Γ(t)對原控制約束進行松弛,將原控制變量約束做如下變換:
(25)
Tmin≤Γ(t)≤Tmax
(26)
在動力學方程和其它非線性約束中,非線性主要來源于氣動升力/阻力、分母項中包含的質量項以及狀態與發動機開機時間和發動機關機時間的乘積。這些非線性特性通過逐次線性化的方法予以線性化。
1)動力學方程約束線性化
在動力學方程約束(13)中,矩陣L是常數矩陣,非線性主要來源于包含終端時間的動力學方程:
(27)
對式(27)在第i次迭代結果處進行一階泰勒級數展開,取其線性項得:
(28)
(29)
式(28)中,在離散點處的常數矩陣為:
將式(28)代入式(13),得到線性化的動力學方程約束為:
(30)
對于等距離散的動力學方程(19),采用相似的離散化方法,此處不再贅述。
由于引入升力控制量,式(16)增加了一個非線性約束:
YT[k]·v[k]=0
(31)
采用一階泰勒級數展開線性化為:
YT[k]·v[k]=Yi[k]·vi[k]+YT[k]·
vi[k]+YiT[k]·v[k]
(32)
2)信賴域約束
由于線性化是在上次迭代結果上進行的,因此若兩次迭代結果相差比較大時,可能使問題無解。為了減少這種風險,希望確保線性化所涉及的變量不會與前一個迭代中獲得的值發生顯著偏離,因此需要引入信賴域約束來保證線性化的可行性。
將第i+1次迭代變量與第i次迭代變量結果之間的差約束到一個范圍內有:
|xi+1-xi|≤ex
(33)
|ui+1-ui|≤eu
(34)
其中,ex和eu分別表示容許的偏差值,可提前設定為固定常數,也可以根據線性化偏差不斷更新。
3)松弛
在逐次迭代的初期,特別是當選取的初始彈道與最優解相差較大時,可能由于線性化造成“偽不可行”現象,即原問題存在可行解,線性化后退化為不可行問題。在這種情況下,不可行性阻礙了迭代過程,阻礙了收斂。一種解決方法是對式(30)添加松弛變量,并把松弛變量加入到性能指標中最小化,松弛約束為:
(35)
式中:Ω(τk)為松弛變量。
由于引入升力控制量,式(16)經過一階泰勒級數展開線性化為(32),該式過于嚴格,在線性化迭代求解中需要進行松弛以獲得更好的迭代結果:
YT[k]·v[k]<=y(k)
(36)
(37)
式中:kΩ,k是系數向量,調節該系數向量使得松弛變量趨于0。
問題1經過無損凸化、線性化和松弛,得到由不等式約束、等式約束和二階錐約束構成的問題2:
s.t.
(38)
原始問題1經過對推力凸化處理后,得到了新的問題2。可以證明,該凸化過程是一種無損凸化,因此問題2與問題1是等價的[1]。
文獻[7]針對非線性動力學、非線性過程約束以及二階錐約束的序列凸化算法進行了一階條件下的收斂性證明。本文中經凸化處理后問題2是一個典型的序列二階錐規劃問題,根據凸優化理論[16],問題2的解序列至少存在一個聚點,該聚點為原問題的駐點。因此,在一階條件下,證明序列凸優化過程是收斂的,并且當信賴域約束動態調整時,具有超線性收斂特性。
為了比較升力對射程的影響,以射程最大和最小為性能指標,構建問題3,其邊界約束與過程約束與問題2有所不同,修改如下。
修改性能指標(38),分別以最大射程和最小射程為性能指標:
最小射程:
最大射程:
終端位置約束式(20)中,由于航向位置作為性能指標,因此去除航向位置約束:
(39)
由于終端位置不再是原點,修改下滑斜率約束式(23),使下滑斜率約束隨終端位置變化:
(40)
雖然性能指標和約束有所不同,但問題3依舊是一個典型的序列二階錐規劃問題。
問題2和問題3是一個序列二階錐規劃問題,需要通過序列凸化算法進行逐次迭代求解。在迭代求解時,首先需要提供一組狀態初始剖面,一組良好的初始剖面有助于迭代過程快速收斂到最優解。

基于該初值,針對問題2進行迭代求解,直到兩次優化結果狀態量的差滿足收斂條件為止。圖3為迭代求解流程圖。
仿真算例分別計算了無升力軌跡優化和升力可控軌跡優化,對比說明了加入升力后對最優燃料性能指標的影響。選擇某型火箭作為仿真模型,任務初值與終端約束如表1所示。火箭飛行參數與氣動參數如表2所示。

參數值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 火箭參數Table 2 Parameters of rocket
本文所有數值計算在MATLAB 2017b環境下編制及運行,采用CVX優化工具箱和MOSEK 8.0.0.60作為算例的凸優化求解工具。
該算例不考慮升力作用,在氣動減速段僅有氣動阻力作用。氣動減速段和動力著陸段各為30個離散點。火箭飛行時間作為優化變量,優化后火箭總飛行時間為135.7 s,其中氣動減速飛行時間為123.7 s,動力著陸段飛行時間為12.0 s。關機后火箭質量剩余14.89 t,總燃料消耗為1.33 t。
圖4為逐次凸化迭代后氣動減速段的三維軌跡,表明火箭成功在預定落點降落,圖中所示圓錐為下滑斜率約束面。從圖5、圖6可以看出,位置和速度均收斂到預定終端狀態。經過空氣阻力減速,發動機開機點的坐標為[114.7, -140.8, 1089] m,初始速度為[-19.38, -0.33, -180.3] m/s。
圖7為火箭發動機推力大小和三向推力曲線。從仿真結果可以看出,氣動減速段推力為0,氣動減速段后,火箭先以最小推力工作,然后切換到最大推力狀態。
本算例考慮氣動升力,且氣動升力大小能夠由總攻角進行調節,總攻角約束為:|η|≤15°。優化后火箭總飛行時間為144.0 s,其中氣動減速段飛行時間為135.1 s,動力著陸段飛行時間為8.9 s。關機時刻火箭剩余質量為15.09 t,燃料消耗為1.13 t。相比無升力軌跡優化燃料節省0.2 t。
圖8、圖9表明速度和位置滿足終端約束。動力著陸段初始的坐標為[222.1, -105.3, 534.8] m,初始速度為[-54.03, 25.65, -125.9] m/s。
從仿真結果可以看出,有升力情況下發動機的開機高度明顯降低,速度矢量的方向更趨近于減少位置誤差的方向,從而導致發動機工作時間縮短,更加節省燃料。
圖10為火箭總攻角η隨時間變化曲線,最大攻角為15°,可以看出攻角滿足約束條件。
通過求解問題3,可以獲得回收過程的最小射程和最大射程。對于無升力情況,最小射程為41.232 km,最大射程為41.867 km,落點射程調節范圍為0.623 km。由于無法對升力進行調節,射程調制主要由推力調節來實現,因此圖11僅給出動力下降段的縱向軌跡曲線。不難發現,由于存在最大燃料、推力指向、推力大小以及下降斜率約束,射程調節范圍較小。
圖12為以升力可控情況的最小、最大射程軌跡。火箭飛行最小射程為36.97 km,最大射程為47.171 km,落點射程調節范圍為10.201 km。相較于無升力控制算例,可大幅提升可達域范圍。
本文針對火箭垂直回收提出一種將氣動減速段與動力著陸段統一規劃的多階段燃料最優軌跡規劃算法。該算法通過將原問題進行無損凸化,并經過線性化和離散化,得到序列SOCP問題,可快速迭代求解。通過仿真校驗,可得出如下結論:
1)將氣動力作為控制變量時,仍可對原問題進行無損凸化。
2)將發動機開、關機時間作為優化變量有助于獲得全局燃料最優解。
3)在氣動減速段調節升力,可以獲得更大的射程調節范圍,并可以進一步節省動力著陸段的燃料。