大氣層內(nèi)固體火箭實時軌跡優(yōu)化方法

郝澤明,張 冉,王嘉煒,李惠峰

(北京航空航天大學宇航學院,北京 100083)

0 引 言

固體火箭能量管理是指通過調(diào)整固體火箭的飛行程序角,使固體火箭對剩余燃料進行消耗,并在燃料耗盡關機時滿足指定的終端狀態(tài)。文獻[1]早期提出了一種能量管理方案,通過在主動段進行姿態(tài)調(diào)制消耗多余的能量。但該方法適用的飛行時間短,且需要離線設計,難以控制終端速度大小。文獻[2]提出了交變姿態(tài)控制能量管理方法(Alternate attitude control energy management,AEM),文獻[3]提出了一般能量管理方法(General energy management,GEM)。文獻[4]在傳統(tǒng)能量管理方法基礎上,提出了一種閉環(huán)動態(tài)逆能量管理方法,得到更有利于姿控系統(tǒng)執(zhí)行的控制指令。這幾種方法均在真空段對固體火箭進行能量管理,且難以嚴格滿足終端高度、速度與傾角的約束。文獻[5]提出了一種多約束制導方法,其可以對固體火箭速度進行管控,但仍只能在真空段進行管理。文獻[6]提出了樣條能量管理方法(Spline energy management,SEM),該方法綜合了AEM與GEM方法,能夠處理一些大氣層內(nèi)過程約束。文獻[7]對傳統(tǒng)能量管理方法進行了改進,提出了一種新的視速度模量計算方法,能夠適應低彈道多約束的能量管理問題。

大氣層內(nèi)固體火箭需要在滿足過程約束情況下機動飛行以達到指定的終端狀態(tài),因此如何嚴格滿足過程約束與終端約束是大氣層內(nèi)能量管理問題的關鍵。上述方法多數(shù)在真空段進行能量管理,僅有少數(shù)方法可用于大氣層內(nèi)的能量管理,且難以嚴格滿足全部約束。同傳統(tǒng)能量管理方法相比,軌跡優(yōu)化方法可以嚴格滿足過程與終端約束,但能量管理問題只需得到一條可行軌跡,不存在明確的性能指標。因此本文構造了性能指標,將能量管理問題轉化為軌跡優(yōu)化問題,并采用軌跡優(yōu)化方法進行在線求解,既對固體火箭進行了能量管理,同時也滿足了各項約束。針對軌跡優(yōu)化問題的求解,常見的方法主要是將問題建模為最優(yōu)控制問題進行求解,包括兩點邊值求解算法[8-9]、非線性規(guī)劃求解算法[10-13]、凸規(guī)劃[14-15]等方法。大氣層內(nèi)固體火箭能量管理問題需要有較高的實時性,而凸規(guī)劃方法有不依賴初始猜想、收斂性穩(wěn)定、求解速度快等特點[16],尤其是對于特定形式的凸規(guī)劃問題,如線性規(guī)劃、二階錐規(guī)劃和半正定規(guī)劃問題,其求解算法成熟可靠并已具備在嵌入式系統(tǒng)上運行的能力。上述特點使凸規(guī)劃方法能夠滿足航天任務中對軌跡優(yōu)化算法在可靠性、實時性與適應嵌入式系統(tǒng)等方面的要求。然而,在解決實際問題的過程中,并非所有的問題都能夠直接建模成凸形式,往往需要對原問題模型進行凸化處理。在基于凸規(guī)劃的火星著陸軌跡優(yōu)化中[17],提出了無損凸化方法,其利用變量代換與凸松弛技術對非凸動力學進行凸化,并應用極大值原理證明了凸規(guī)劃問題與原問題解的等價性,該無損凸化方法已廣泛應用于多種場景的軌跡優(yōu)化問題中[18-19]。

對于大氣層內(nèi)上升段問題,固體火箭飛行受氣動力影響,且飛行過程中存在攻角、動壓、過載等過程約束,使原軌跡規(guī)劃問題具有較強的非線性。氣動力與火箭推力不同,其通常與攻角、馬赫數(shù)等直接有關,而這兩者又與狀態(tài)量和控制量直接耦合,因而尚無成熟完備的無損凸化技術可以處理強非線性動力學與過程約束。針對這種情況,通常在給定參考軌跡的情況下對非線性動力學與約束進行線性化處理,并根據(jù)求解結果更新參考軌跡進行迭代求解,這種方法被稱為序列凸規(guī)劃方法[20-22]。對于固體火箭,其推力與秒耗量通常為時間的函數(shù),即飛行過程中任意時刻的推力幅值與質(zhì)量均為已知。由于這種推力不可調(diào)的特性,在氣動力與過程約束存在的情況下,傳統(tǒng)的無損凸松弛技術[17]難以滿足,可能會導致問題的求解結果不滿足真實動力學。

本文的主要創(chuàng)新點在于采用軌跡優(yōu)化的思路求解固體火箭能量管理問題,其能夠系統(tǒng)的處理約束并增強算法的適應性。在將能量管理問題轉化為軌跡優(yōu)化問題時,通過規(guī)則化方法主動構造了性能指標,并采用鄰近-牛頓-康托維奇凸規(guī)劃、控制平滑以及虛擬控制等技術,改善了算法性能：首先,在性能指標中引入推力方向的鄰近規(guī)則化項,顯著提高了算法的收斂性；然后,在優(yōu)化問題中加入了推力方向變化角速率約束與控制平滑指標,并在離散系統(tǒng)中通過數(shù)值微分表示推力方向的變化,同傳統(tǒng)方法相比,該方法避免出現(xiàn)奇異弧并提高了求解速度；最后,在推力方向幅值約束與動壓約束中引入虛擬控制變量并在性能指標中對其進行懲罰,保證了凸規(guī)劃子問題在迭代過程中的可行性,降低算法對初始猜想的依賴。仿真實例表明,采用所提出的軌跡優(yōu)化方法求解能量管理問題是有效的,能夠嚴格滿足各項約束并實現(xiàn)高精度終端；此外,算法在地面平臺具有百毫秒級的運算速度,具有優(yōu)異的實時性,能夠滿足能量管理對求解速度的要求。

1 大氣層內(nèi)固體火箭能量管理問題描述

本文研究固體火箭上升段在縱向平面內(nèi)的質(zhì)心運動,假設發(fā)動機推力始終沿體軸方向,發(fā)動機推力幅值和秒耗量全程已知,忽略地球自轉并采用球形重力場。定義平面大地坐標系OXYZ：以發(fā)射點為原點O,X軸指向火箭發(fā)射方向,Y軸沿鉛錘方向指向天空,Z軸由右手定則確定。固體火箭在平面大地坐標系下無量綱化(下文中的變量均為無量綱化的)運動學方程為：

(1)

式中：r∈3表示位置；V∈3表示速度；T表示發(fā)動機推力幅值;m表示固體火箭質(zhì)量,其均為時間的函數(shù)且已知；1T∈3表示發(fā)動機推力方向并滿足3為縱向平面內(nèi)與1T垂直的單位矢量；A與N分別表示軸向氣動力大小與法向氣動力大小；g∈3表示火箭所受重力加速度；軸向氣動力大小與法向氣動力大小的表達式如下：

(2)

式中：ρ(r)表示r處的大氣密度；Sref表示參考面積；V表示速度模值;Cx與Cy分別為軸向力系數(shù)與法向力系數(shù),其均為馬赫數(shù)Ma與攻角α的二元函數(shù)。以位置與速度作為狀態(tài)量,定義：

(3)

以推力方向作為控制量,將運動學方程寫為緊湊形式：

(4)

固體火箭上升段在飛行過程中應嚴格滿足過程約束,包括攻角約束、動壓約束以及法向過載約束,其表示方法如下：

|α|≤αmax

(5)

(6)

(7)

式中：αmax表示允許的最大攻角；qmax表示允許的最大動壓；nmax表示允許的最大法向過載。值得注意的是,由于攻角可以表示為推力方向與速度的函數(shù)：cosα=VT1T/V,因此攻角范圍約束可以表示為：

VT1T≥Vcosαmax

(8)

本文考慮的終端約束包括終端高度范圍約束、終端兩軸速度范圍約束,其可以表示如下：

(9)

綜上所述,本節(jié)定義了大氣層內(nèi)上升段固體火箭能量管理問題。解決此類問題一般采用開環(huán)姿態(tài)調(diào)制或閉環(huán)能量管理的方法,此類方法通常難以處理過程約束,因而一般選擇在真空段進行能量管理。軌跡優(yōu)化方法能夠滿足嚴格的過程約束與終端約束,其可以用于解決大氣層內(nèi)上升段問題。因此,本文在下一節(jié)中將固體火箭能量管理問題轉化為軌跡優(yōu)化問題進行求解。

2 大氣層內(nèi)固體火箭能量管理軌跡優(yōu)化問題建模

針對傳統(tǒng)能量管理方法難以嚴格滿足大氣層內(nèi)過程約束的不足,本節(jié)提出一種對控制量模值積分進行懲罰的規(guī)則化方法,該方法主動構造了一個性能指標,將大氣層內(nèi)固體火箭上升段能量管理問題轉化為了軌跡優(yōu)化問題。

(10)

并通過規(guī)則化[23-24]方法構造性能指標：

(11)

式中：t0為飛行初始時刻；tf為飛行終端時刻；γ1>0為懲罰因子。因此,大氣層內(nèi)固體火箭上升段軌跡優(yōu)化問題為：

問題1:大氣層內(nèi)固體火箭上升段軌跡優(yōu)化問題

(12)

式中：t0為飛行初始時刻；tf為飛行終端時刻；fa(α,Ma)為氣動合力,其只與攻角α和馬赫數(shù)Ma相關；r0為初始位置；V0為初始速度。

該規(guī)則化方法對控制量的等式約束進行了松弛,將非凸控制約束轉化為了凸約束,易于軌跡優(yōu)化問題的求解。同時,當懲罰因子γ1足夠大時,可以證明這種松弛方式是無損的。

證.大氣層內(nèi)固體火箭上升段軌跡優(yōu)化問題是一個最優(yōu)控制問題,其Hamiltonian函數(shù)為：

(13)

根據(jù)極大值原理,有:

(14)

不妨設協(xié)態(tài)變量λV與控制量1T的夾角為Θ1,與氣動合力fa的夾角為Θ2,則其示意圖如圖1所示。

圖1 控制變量與協(xié)態(tài)變量示意圖Fig.1 Illustration of control variables and costate variables

則問題1的Hamiltonian函數(shù)可以表示為：

(15)

當懲罰因子γ1足夠大時,有:

(16)

3 鄰近-牛頓-康托維奇凸規(guī)劃方法

根據(jù)Newton-Kantorovich迭代方法[25],對問題1的運動學方程進行廣義泰勒展開,得到：

(17)

考慮推力方向變化角速率約束,其可以在后續(xù)離散化過程中直接處理；而對于攻角約束,由于迭代過程中固體火箭的狀態(tài)量變化相對劇烈,而控制量變化相對平緩,因此為保證問題的收斂性,僅對控制量進行廣義泰勒展開。定義x(i)為當前迭代對應的廣義泰勒展開點,則攻角約束為：

(18)

考慮除上述約束外的過程約束,將其表示為緊湊形式：

ψp(t,x,1T)≤0

(19)

終端不等式約束的緊湊形式為：

ψf(tf,xf)≤0

(20)

將上述約束進行廣義泰勒展開,得到：

(21)

對性能指標進行廣義泰勒展開,得到：

(22)

在性能指標中加入鄰近項,形成鄰近-牛頓-康托維奇子問題：

問題2：鄰近-牛頓-康托維奇子問題

(23)

式中：γ2>0為鄰近項對應的懲罰因子。

鄰近-牛頓-康托維奇凸規(guī)劃方法的廣義泰勒展開點更新方式如下：

x(i+1)=x(i)+δx

(24)

(25)

若鄰近-牛頓-康托維奇子問題的迭代求解收斂,則其可收斂至原軌跡優(yōu)化問題的局部最優(yōu)解。將大氣層內(nèi)固體火箭上升段軌跡優(yōu)化的最優(yōu)控制問題轉化為鄰近-牛頓-康托維奇子問題迭代求解,可以將強非線性的運動學方程與約束條件轉化為線性形式,便于采用凸規(guī)劃方法實時求解。

4 實時軌跡優(yōu)化算法實現(xiàn)

本節(jié)介紹鄰近-牛頓-康托維奇凸規(guī)劃的具體實現(xiàn)方法。首先,采用等時間間隔離散方法：

(26)

運動學方程的離散采用一階保持器,其假設控制量在相鄰兩離散點間是線性變化的,即:

(27)

因此,離散后的動力學方程為：

(28)

其中,狀態(tài)轉移矩陣Φ(tk+1,tk)滿足：

(29)

(30)

式中：

(31)

針對過程約束,其離散形式為：

(32)

針對性能指標,采用歐拉法進行數(shù)值積分得到離散化的性能指標：

(33)

在實際工程應用中,軌跡優(yōu)化所得結果應易于被姿態(tài)控制系統(tǒng)跟蹤,因此考慮加入推力方向變化角速率約束,使得優(yōu)化所得推力方向不會出現(xiàn)較大的跳變。針對該約束,通常在連續(xù)系統(tǒng)中以推力方向的導數(shù)(或變化角速度)為控制量,以推力方向為狀態(tài)量建立最優(yōu)控制問題。然而這種處理方式增加了問題的維度,并且引入了額外的狀態(tài)約束,使問題的求解變得更加困難。本文的處理方式是在離散系統(tǒng)中以數(shù)值微分對推力方向進行約束,在引入規(guī)則化方法的情況下,其可以避免奇異弧且具有更快的求解速度,類似的處理方式也出現(xiàn)在文獻[26] 中。離散系統(tǒng)中的推力方向幾何示意圖如下：

圖2 推力方向幾何示意圖Fig.2 Geometry schematic of thrust direction

根據(jù)圖中幾何關系,有：

(34)

式中：ωφ為推力方向變化角速率。因此,推力加速度變化加速度約束可以表示為：

(35)

其為二階錐約束,可以采用凸優(yōu)化方法直接求解。

此外,為使優(yōu)化所得推力方向變化盡可能平滑,設計如下性能指標：

(36)

因此,離散化后的子問題為：

問題3:離散Proximal-Newton-Kantorovich子問題

(38)

根據(jù)鄰近-牛頓-康托維奇凸規(guī)劃迭代求解方法,需要給定初始軌跡以開始迭代。而由于原問題具有較為嚴苛的過程約束,因而任意給定的初始軌跡易使問題出現(xiàn)不可行,從而導致迭代終止。因此,在控制約束和過程約束中引入虛擬控制變量κ1、κ2：

(39)

并在性能指標中加入相應的懲罰項：

maxJ2=-γ4κ1-γ5κ2

(40)

式中：γ4,5為正實數(shù)。二階錐規(guī)劃問題為：

問題4:二階錐規(guī)劃問題

maxJ=J1+J2

(41)

從數(shù)學角度考慮,引入的虛擬控制變量實際上是一種精確罰函數(shù),其具有如下性質(zhì)：

性質(zhì)1.若原問題存在可行解,且γ4與γ5足夠大時,虛擬控制變量κ1與κ2的收斂解為零,引入的虛擬控制變量對原問題的求解結果無影響。該性質(zhì)的證明如下。

(42)

(43)

式中：

(44)

為使Lagrange對偶函數(shù)有界,需要滿足：

(45)

因此,當γ4與γ5滿足：

(46)

則虛擬控制變量κ1與κ2的為零。

文中描述的二階錐規(guī)劃問題,其可以采用ECOS軟件[27]求解。針對上述算法設計迭代終止條件為：控制量增量與虛擬控制變量小于容許誤差,即

σ2≤ε1,κ1,2≤ε2

(47)

式中：ε1為控制量增量容許誤差；ε2為虛擬控制變量容許誤差。

5 仿真校驗

本節(jié)采用大氣層內(nèi)固體火箭上升段軌跡優(yōu)化問題的具體實例對所提出的基于鄰近-牛頓-康托維奇凸規(guī)劃的軌跡優(yōu)化方法進行仿真校驗。本文所采用的仿真數(shù)據(jù)為：地球半徑為6371004 m,地球引力常數(shù)為3.986005×1014m3/s2,大氣密度模型和聲速模型擬合自美國標準大氣1976；固體火箭的初始質(zhì)量為5400 kg,發(fā)動機推力為130 kN,秒耗量為50 kg/s,參考面積為0.5 m2,飛行時間為75 s,軸向力系數(shù)和法向力系數(shù)為：

(48)

式中：

(49)

固體火箭初始時刻以10 m/s的速度垂直起飛；在飛行過程中,要求飛行攻角絕對值全程不大于10°,終端攻角絕對值不大于1°；飛行過程中的動壓不超過70 kPa,法向過載不超過4。在用凸規(guī)劃求解該問題時,取離散點數(shù)量N為40,控制量增量容許誤差ε1取1×10-1；虛擬控制變量容許誤差ε2取1×10-4。仿真選擇如表1所示的三種工況,其分別表示能量管理能達到的最小、中間及最大終端速度,這三種工況的仿真結果如下。

表1 仿真工況Table 1 Chemical composition of experimental alloy

圖3給出了固體火箭三種工況下的縱向平面內(nèi)的位置與速度曲線。從圖中可知,固體火箭能夠嚴格滿足終端高度約束與終端兩軸速度約束。圖4分別給出了固體火箭飛行過程中的俯仰角與攻角曲線。從圖中可知,固體火箭在工況一情況下需要達到較小的終端速度,因而軌跡優(yōu)化方法優(yōu)化得到了擺動形式的姿態(tài)角。軌跡優(yōu)化所得結果與傳統(tǒng)能量管理方法的姿態(tài)調(diào)制結果相似,但能夠嚴格滿足終端約束,并且無需人為設計調(diào)制方式,在大氣層內(nèi)上升段其終端速度具有約10%的調(diào)節(jié)能力。

圖3 縱向平面內(nèi)位置與速度軌跡Fig.3 Position and velocity trajectories in longitudinal plane

圖4 俯仰角與攻角曲線Fig.4 Pitch angel and angel-of-attack

圖5表示三種工況下的過載與動壓曲線。從圖中可知,固體火箭飛行全程均滿足過載與動壓約束,表示軌跡優(yōu)化方法所得軌跡滿足過程約束。同時,為了驗證軌跡優(yōu)化方法的優(yōu)化精度,根據(jù)三種工況下軌跡優(yōu)化所得俯仰角曲線,對其進行開環(huán)數(shù)值積分,得到如表2所示的優(yōu)化誤差。從表中可知,軌跡優(yōu)化方法所得的終端高度與速度均具有米級精度,表示軌跡優(yōu)化方法所得結果滿足動力學約束。

圖5 過載與動壓曲線過載系數(shù)Fig.5 Overload and dynamic pressure

表2 軌跡優(yōu)化誤差Table 2 Trajectory optimization errors

為進一步驗證本文所提出軌跡優(yōu)化方法的正確性與有效性,以工況1為例進行仿真與對比。圖6所示為攻角與動壓對比結果。首先考慮本文提出的控制平滑方法,從圖中可知,加入了控制平滑指標的飛行攻角相比未加入的更為平滑,這是由于控制平滑指標是在離散系統(tǒng)中對推力方向二階差分模值積分進行了懲罰,這種處理方式在連續(xù)系統(tǒng)中相當于最小化推力方向二階導數(shù)的變化率。因此在這種處理方式下,推力方向的變化更為平緩,這在飛行攻角曲線上的體現(xiàn)為：加入了控制平滑指標的攻角曲線其拐點較少,而不加入該指標的攻角曲線拐點較多,且加入了控制平滑指標的曲線拐點兩側攻角變化角速度突變更小。從數(shù)學角度考慮,光滑函數(shù)要求函數(shù)無窮階連續(xù)可微,而離散系統(tǒng)中的拐點表示該點連續(xù)但不可微,因此可以認為拐點較少的攻角曲線其光滑程度更好。從實際工程角度考慮,拐點前后的攻角角速度若具有較大的突變,姿態(tài)控制系統(tǒng)在跟蹤該指令時,其跟蹤效果會受到較大影響。若軌跡優(yōu)化所得攻角曲線具有較多拐點且拐點前后攻角角速度具有較大突變,則該軌跡在實際工程中不利于姿態(tài)控制系統(tǒng)跟蹤。因此,本文考慮的控制平滑指標可以提高軌跡優(yōu)化所得控制量的光滑程度,使其在實際應用中更易于被姿態(tài)控制系統(tǒng)跟蹤。

圖6 攻角與動壓對比Fig.6 Comparison of angle-of-attack and dynamic pressure

考慮固體火箭上升段需要滿足的過程約束。以動壓約束為例,從圖6中可知,在不加入動壓約束情況下,飛行過程中的最大動壓約為75 kPa,高于任務要求的70 kPa。在加入動壓約束情況下,固體火箭全程動壓均不超過動壓上限約束,表明動壓約束產(chǎn)生實際作用,本文提出的軌跡優(yōu)化方法能夠嚴格滿足過程約束。

圖7給出了軌跡、攻角與虛擬控制變量的部分迭代過程。從圖中可知,初始軌跡與最終收斂軌跡的差異較大,但本文所提出的軌跡優(yōu)化方法能夠在初始軌跡偏差較大的情況下使軌跡收斂。值得注意的是,第一次迭代中的虛擬控制變量均不為零,表示初始軌跡偏差較大,若不加入虛擬控制變量則問題無可行解,導致迭代無法進行。

圖7 軌跡、攻角與虛擬控制變量的迭代過程Fig.7 Trajectory,angle-of-attack and virtual control variables in iterations

傳統(tǒng)能量管理方法多數(shù)為制導方法,通常其難以直接處理過程約束,且多數(shù)在固體火箭三級進行能量管理。本文采用軌跡優(yōu)化則與之不同,其能夠直接對過程約束進行處理,并能夠在固體火箭一級進行能量管理。因此為驗證所提出軌跡優(yōu)化方法的性能,以工況1為例,采用偽譜法求解能量管理問題得到可行軌跡,與本文的軌跡優(yōu)化結果進行對比。

本文以縱向平面內(nèi)的為位置、速度、俯仰角以及俯仰角速率作為狀態(tài)量,以俯仰角角加速度作為控制量,分別以最小化俯仰角角加速度模值積分為性能指標,采用GPOPS軟件[28]求解軌跡優(yōu)化問題,其仿真對比結果如圖8所示。

從仿真結果可以看出,采用偽譜法能夠求解本文所提出的軌跡優(yōu)化問題,實現(xiàn)了對固體火箭的能量管理,同時滿足了動壓、過載等過程約束。這表明本文通過構造性能指標將能量管理問題轉化為軌跡優(yōu)化問題的方法是正確可行的。從仿真對比結果可以看出,偽譜法優(yōu)化所得結果與本文優(yōu)化結果存在差異,這是由于能量管理問題存在多個局部最優(yōu)解,兩種不同的非線性規(guī)劃方法可以得到不同的可行解。最終,在i7-8565 CPU 1.8 GHz環(huán)境下,利用GPOPS軟件求解能量管理問題,其離散點數(shù)量為329,求解時間為6.79 s；與此相比,在同一環(huán)境下運行本文所提出的軌跡優(yōu)化算法,其平均計算耗時為83 ms。因此,本文所提出方法在對固體火箭進行能量管理的同時,具有良好的收斂性與更快的計算速度。

6 結 論

本文提出一種鄰近規(guī)則化技術,將固體火箭能量管理問題轉化為軌跡優(yōu)化問題,并提出一種鄰近-牛頓-康托維奇凸規(guī)劃的軌跡優(yōu)化方法,實時求解了大氣層內(nèi)固體火箭上升段軌跡優(yōu)化問題。結論包括：

本文提出的規(guī)則化技術將固體火箭能量管理問題轉化為軌跡優(yōu)化問題,并應用最優(yōu)控制理論對其進行了理論分析,證明軌跡優(yōu)化問題的解是能量管理問題的可行解；本文提出一種鄰近-牛頓-康托維奇凸規(guī)劃方法,避免算法求解出現(xiàn)奇異弧、提高了求解速度并保證了凸規(guī)劃子問題在迭代過程中的可行性。數(shù)值仿真結果表明,采用該方法求解軌跡優(yōu)化問題能夠滿足動力學與各項約束,優(yōu)化所得控制量易于被姿態(tài)控制系統(tǒng)跟蹤,算法具有良好的收斂特性與實時性且對初始猜想不敏感。