基于Proximal-Newton-Kantorovich凸規劃的空天飛行器實時軌跡優化

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

北京航空航天大學 宇航學院，北京 100083

空天飛行器是一種可重復使用的運載系統，提供了一種經濟、便捷地進入空間的途徑。與運載火箭相比，空天飛行器在上升段中所受的氣動力更為顯著，主動利用氣動力可提升空天飛行器的運載效率。上升段中，為最大化利用氣動力的作用，需要對上升段的軌跡進行優化設計。空天飛行器上升段的軌跡優化問題已被廣泛研究[1]，軌跡優化方法包括間接法[2]、偽譜法[3]、混合優化[4]、軌道設計反方法[5]等。這些方法將軌跡優化問題轉化為非線性規劃問題，采用牛頓法、序列二次規劃方法等非線性規劃方法求解。然而，非線性規劃方法具有指數時間復雜度，計算時間難以預估，且需要設置合理的初始猜想。這些特點限制了基于非線性規劃方法的軌跡優化方法在快速發射任務和應急發射任務中的使用。

近年來，基于凸規劃方法的軌跡優化方法是一個研究熱點。凸規劃方法具有多項式時間復雜度，計算時間上界可確定的特點[6-7]。此外，對一類具有特殊結構的凸規劃問題，即線性規劃、二階錐規劃和半正定規劃問題，應用齊次內嵌自對偶內點法可高效求解，且不需要人為設置初始猜想，當問題不可行或無界時能夠在有限迭代內給出明確判據。上述特點使得凸規劃方法適合用于求解航空航天背景中的問題，若問題能夠轉換為凸規劃問題，則求解時間有明確的上界，可以滿足航空航天問題對可靠性、實時性的要求。然而，與非線性規劃方法不同，凸規劃方法對問題結構的要求較為嚴格。具體到軌跡規劃問題，凸規劃方法其要求軌跡規劃問題具有線性的性能指標、線性的運動方程、線性的等式約束和凸的不等式約束。只有很少的軌跡優化問題滿足上述要求，因此限制了凸優化方法的應用范圍。在動力下降制導[8-10]中，提出了無損凸化技術：通過變量替換和等式約束松弛，構造了一個凸規劃問題，并應用極大值原理，證明了此凸規劃問題與原軌跡規劃問題的等價性。無損凸化技術解決的是火星著陸器動力下降段的軌跡規劃問題，但此技術可直接應用至各種火箭動力飛行器的軌跡規劃問題，目前已有很多相關成果[11-13]。

然而，帶有氣動力的軌跡規劃問題如何應用凸規劃方法求解仍是難點。與火箭推力不同，氣動力呈高度非線性：三維空間中，火箭推力的可行域是一個球殼，無損凸化技術將可行域轉化為球，從而獲得凸性；氣動力的可行域則是一個復雜的曲面，且曲面形狀與狀態量(高度、馬赫數)直接相關，目前沒有成熟的無損凸化技術。一類處理氣動力的方法是以推力方向作為控制量，將氣動力列入運動方程[12,14-17]。針對運載火箭上升段的軌跡優化問題，提出了一種組合應用無損凸化技術與Newton-Kantorovich迭代的軌跡優化方法[12,15]。此類方法無法處理升力，并且假設阻力系數為常數，因此只能應用于沒有升力外形的飛行器，如火星著陸器或運載火箭，無法應用于空天飛行器這類氣動力更為顯著的飛行器。另一類方法是將推力方向和氣動系數作為控制量，并同時對兩者的控制約束進行凸化[18-19]。此類方法嘗試將無損凸化技術擴展至氣動力：通過假設升力系數和阻力系數間存在的特定函數關系，將氣動系數建模為控制量，并將對其的等式約束松弛為凸不等式約束。然而，受限于氣動力的非線性和過程約束的存在，無法理論證明凸化是無損的，即此方法的收斂結果可能不滿足飛行器的運動方程。

本文提出一種基于Proximal-Newton-Kantorovich凸規劃的軌跡優化方法。在對空天飛行器軌跡優化問題進行建模時，選擇攻角作為唯一控制量，推力和氣動力都由攻角確定。Proximal-Newton-Kantorovich凸規劃方法首先應用Newton-Kantorovich迭代方法將軌跡優化問題轉換為一系列的子問題，子問題中的性能指標、運動方程和約束都是線性的。隨后，針對子問題求解時控制量始終處在邊界，導致迭代無法收斂的問題，提出鄰近規則化方法予以解決，應用最優控制理論證明了Proximal-Newton-Kantorovich迭代方法的收斂結果一定是軌跡優化問題的局部最優解；最后，將規則化子問題離散為二階錐規劃問題，并應用內點法進行求解。此方法外環為Proximal-Newton-Kantorovich迭代方法，內環為凸規劃方法，稱作Proximal-Newton-Kantorovich凸規劃方法。

本文的主要貢獻是提出了Proximal-Newton-Kantorovich迭代方法，并驗證了此方法是一種實時求解非線性軌跡規劃問題的可行途徑。說明除無損凸化方法以外，Newton-Kantorovich迭代方法也可將非線性軌跡規劃問題轉化為凸優化問題。與無損凸化方法相比，此方法對軌跡規劃問題結構要求較低，不需要做連續可微以外的假設。實例表明Proximal-Newton-Kantorovich迭代方法可求解空天飛行器大氣層內上升段軌跡優化這一非線性的軌跡優化問題，收斂結果與偽譜法提供的最優參考相近，且計算時間在毫秒級。

1 空天飛行器上升段軌跡優化問題建模

本文研究固體火箭助推、垂直起飛的空天飛行器在縱向平面內的質心運動，忽略地球的旋轉和扁平率，對飛行器做如下假設：

假設1推力幅值和秒耗量是飛行時間的函數。

假設2推力方向沿飛行器的結構縱軸指向前。

假設3升力系數和阻力系數是馬赫數和攻角的函數，并且連續可微。

空天飛行器在縱向平面內的質心運動方程為

(1)

將運動方程的控制量選為攻角。重力、推力、升力和阻力可寫作飛行時間、位置、速度和攻角的函數。重力G可寫為

圖1 空天飛行器的運動和受力狀態Fig.1 Motion and forces of a hypersonic vehicle

(2)

式中：μ為地球引力常數。

推力T可寫為

(3)

式中：T為推力幅值；Φ(·)為旋轉矩陣，表達式為

(4)

(5)

(6)

(7)

本文研究的空天飛行器上升段軌跡優化問題為：給定空天飛行器的初始位置和速度，在滿足某些終止約束和攻角約束的條件下，最小化某性能指標，即：

問題1空天飛行器上升段軌跡優化問題

2 Newton-Kantorovich迭代方法

空天飛行器上升段軌跡優化問題是一個非線性最優控制問題，可應用Newton-Kantorovich迭代方法求解。首先簡要介紹Newton-Kantorovich迭代方法[20-21]。Newton-Kantorovich迭代方法是Newton迭代方法向泛函空間的擴展：Newton迭代方法求解的是代數方程，而Newton-Kantorovich迭代方法求解的是微分方程。兩者的核心依據都是Taylor展開，下面應用類比的方式說明。Newton 迭代方法求解非線性代數方程

w(x)=0

(8)

核心依據是非線性函數w(x)的Taylor展開：

(9)

Newton迭代方法忽略Taylor展開中的高階項，求解關于δ的線性代數方程:

(10)

獲得δ。求解完成后，令

x+=x+δ

(11)

直到迭代收斂，式中：x+為下一次迭代中的Taylor展開點。相應的用Newton-Kantorovich迭代方法求解非線性微分方程

(12)

時，核心依據是非線性算子

(13)

的廣義Taylor展開：

(14)

(15)

Newton-Kantorovich迭代方法忽略廣義Taylor展開中的高階項，求解關于δ的線性算子方程：

(16)

(17)

求解此線性微分方程獲得δ后，令

x+=x+δ

(18)

直到迭代收斂。可見，Newton-Kantorovich迭代方法將非線性微分方程轉化為一系列的線性微分方程求解。與求解非線性微分方程的其他方法相比，區別在于線性化和離散化的順序不同：Newton-Kantorovich迭代方法先將非線性微分方程轉化為線性微分方程(線性化)，再離散為線性代數方程進行求解(離散化)；其他方法先將非線性微分方程離散為非線性代數方程(離散化)，再應用Newton迭代方法或其變種進行求解(線性化)。Newton-Kantorovich迭代方法先進行線性化的優勢在于其顯式地計算了Jacobian矩陣。如果先進行離散化，非線性代數方程由Newton迭代方法或其變種求解時會隱式地應用數值差分方法計算Jacobian矩陣，而數值差分引起的誤差會降低優化方法的收斂速度[22]。

(19)

非線性算子的廣義泰勒展開為

(20)

(21)

(22)

忽略廣義Taylor展開中的高階項，線性算子方程

(23)

可寫作線性微分方程：

(24)

為簡明起見，函數f及其Jacobian矩陣的自變量被隱去。

在空天飛行器上升段軌跡優化問題中，將關于r和α的非線性微分方程替換為關于κ和τ的線性微分方程，并將性能指標φ和終止約束ψ替換為其一階Taylor展開，即得到Newton-Kantorovich子問題：

問題2Newton-Kantorovich子問題

r+=r+κ

(25)

(26)

α+=α+τ

(27)

Newton-Kantorovich子問題在離散化后形成線性規劃問題。線性規劃問題是凸規劃問題的一種，可采用內點法在多項式時間內求解具有實時性保障[6-7]。然而，Newton-Kantorovich子問題是一個線性最優控制問題，且不等式約束僅施加在控制量τ上，所以其最優解的控制量一定處于邊界。攻角增量τ的飽和會導致下次迭代中攻角α處于邊界，進一步阻礙了迭代的收斂。這一問題出現的原因是Newton-Kantorovich迭代方法忽略了微分方程中的高階信息。本文在第3節中提出Proximal-Newton-Kantorovich迭代方法解決這一問題。

3 Proximal-Newton-Kantorovich迭代方法

本文提出Proximal-Newton-Kantorovich迭代方法，在Newton-Kantorovich子問題的性能指標中加入鄰近規則化項[23]構成Proximal-Newton-Kantorovich子問題：

問題3Proximal-Newton-Kantorovich子問題

鄰近規則化項的加入使得子問題成為非線性最優控制問題，避免了線性最優控制問題控制量始終處于邊界的問題。Proximal-Newton-Kantorovich子問題離散化后構成二次規劃問題，可應用凸規劃方法求解，此時稱作Proximal-Newton-Kantorovich凸規劃方法。

(28)

1) 微分方程

(29)

(30)

(31)

2) 邊界條件

(32)

ψ=0

(33)

(34)

(35)

3) 控制條件

λα1≥0,α1-α≤0,λα1(α1-α)=0

(36)

λα2≥0,α-α2≤0,λα2(α-α2)=0

(37)

(38)

Newton-Kantorovich迭代方法收斂后，子問題的必要條件和原問題的必要條件存在緊密聯系，可通過子問題的Lagrange乘子構造原問題的Lagrange乘子。子問題的Hamiltonian函數為

(39)

1) 微分方程

(40)

(41)

(42)

2) 邊界條件

(43)

(44)

(45)

(46)

3) 控制條件

(47)

(48)

(49)

4 數值實現

本節介紹Proximal-Newton-Kantorovich凸規劃方法的數值實現方法。Proximal-Newton-Kantorovich凸規劃方法的步驟可總結為

步驟2求解 Proximal-Newton-Kantorovich 子問題。

下面分別介紹初始猜想選擇、子問題求解方法和終止條件。

應用Euler法將Proximal-Newton-Kantorovich子問題離散為二階錐規劃問題。離散點數為N，離散時間均勻分布：

(50)

優化變量選擇為位置增量、速度增量和攻角增量：

(51)

并引入一個松弛變量，用于將規則化方法中的二次項轉化為二階錐約束：

σ1,σ2,…,σN-1

(52)

離散化后的問題為

問題4二階錐規劃問題

式中：k為規則化系數。通過求解此問題，可獲得 Proximal-Newton-Kantorovich 子問題的近似解。二階錐規劃問題應用 ECOS 求解[24]。

迭代終止條件設計為：位置增量、速度增量和攻角增量均小于設置的容差，即

(53)

(54)

(55)

5 數值實驗和結果分析

本節中應用 Proximal-Newton-Kantorovich 凸規劃方法求解空天飛行器上升段軌跡優化問題的一個實例，以說明方法的有效性。

首先，對空天飛行器上升段軌跡優化問題進行建模。地球半徑取6 371.004 km；地球引力常數μ取3.986 005×1014km；大氣密度模型和聲速模型擬合自美國標準大氣1976[25]，大氣密度模型在0～130 km范圍內擬合誤差不大于4%，聲速模型在0～86 km范圍內擬合誤差不大于2%。空天飛行器的總飛行時間為75 s；初始質量為4 000 kg；秒耗量為40 kg/s；推力幅值為100 kN；參考面積為0.8 m2；升力系數和阻力系數擬合為

CL=(4+4Ma-0.4Ma2)α

(56)

CD=0.1+Ma-0.2Ma2+0.02Ma3+

(5+0.5Ma)α2

(57)

升力系數和阻力系數也可由其他形式的擬合公式或連續可微的插值方法計算得到。攻角范圍為-10°～10°。初始時，飛行器以100 m/s的速度垂直起飛；終止時，飛行器以27 km的高度水平飛行；性能指標為最大化終止速度，即

(58)

(59)

式中：hf為27 km。

圖2 軌跡對比Fig.2 Comparison of trajectories

圖3 速度對比Fig.3 Comparison of velocities

圖4 攻角對比Fig.4 Comparison of angles-of-attack

圖5 飛行路徑角對比Fig.5 Comparison of flight path angles

圖2給出了空天飛行器的軌跡，可見飛行器垂直起飛，逐漸轉入水平飛行，終止高度為27 km，滿足了終止約束。圖3給出了速度的軌跡，軌跡優化的性能指標是最大化終止速度，收斂結果與最優參考的終止速度一致，說明了軌跡規劃方法的最優性。圖4給出了攻角的軌跡，收斂結果與最優參考的攻角變化規律一致，但存在一定差異，這是離散誤差導致的：Euler法的離散精度不及偽譜法，且本文的離散點數較少(40個)，導致了攻角軌跡的不同。圖5給出了飛行路徑角的軌跡，收斂結果與最優參考十分相近，飛行器由垂直起飛轉入水平飛行。總而言之，所有的終止約束和過程約束都得到滿足，且收斂結果與最優參考相近，這說明了Proximal-Newton-Kantorovich迭代方法的有效性和最優性。

圖6～圖8給出了初始猜想和前5次迭代的軌跡。可見初始猜想與最優參考的差距較大，但Proximal-Newton-Kantorovich迭代方法能夠收斂至最優參考。每一次迭代后軌跡都更接近最優參考，而第5次迭代后求得的解已與最優參考十分相近。工程應用中，若是采用閉環控制律跟蹤優化出的軌跡，實際上可放寬軌跡優化的容差，以進一步降低計算時間。

圖9給出了迭代中位置增量、速度增量和攻角增量。可見所有的增量都單調下降。第14次迭代后，所有增量小于設置的容差，迭代停止。

圖6 軌跡的迭代過程Fig.6 Trajectories in iterations

圖7 速度的迭代過程Fig.7 Velocities in iterations

圖8 攻角的迭代過程Fig.8 Angles-of-attack in iterations

圖9 增量的迭代過程Fig.9 Residuals in iterations

6 結 論

本文提出了一種基于Proximal-Newton-Kantorovich凸規劃的軌跡優化方法，用于實時求解空天飛行器大氣層內上升段軌跡優化問題。結論包括：

1) Proximal-Newton-Kantorovich迭代方法中加入的鄰近規則化項改善了Newton-Kantorovich迭代方法的收斂性。應用最優控制理論可證明Proximal-Newton-Kantorovich迭代方法的收斂結果一定是軌跡優化問題的局部最優解。數值實驗表明此方法的收斂結果與最優參考相近。

2) 應用Proximal-Newton-Kantorovich迭代方法將軌跡優化問題轉化為一系列的Proximal-Newton-Kantorovich子問題，子問題離散化后構成二階錐規劃問題，可應用凸規劃方法中的內點法求解，從而保障了軌跡優化方法的實時性。數值實驗表明此方法的計算時間在毫秒級。

上述結論表明，本文提出的Proximal-Newton-Kantorovich凸規劃方法是一種實時求解非線性軌跡規劃問題的可行途徑。