基于hp自適應(yīng)偽譜法的固體運載火箭軌跡優(yōu)化

洪 蓓 辛萬青

北京宇航系統(tǒng)工程研究所，北京100076

多級固體運載火箭(Multistage Solid Launch Vehicle，MSLV)具有響應(yīng)快速、機(jī)動性強(qiáng)、成本低、可靠性高的特點，目前已成為各國研制的熱點。運載能力是運載火箭設(shè)計的重要指標(biāo)，而主動段軌跡優(yōu)化設(shè)計對提高M(jìn)SLV 運載能力具有重要的意義。

MSLV 主動段軌跡優(yōu)化實質(zhì)上是一類終端狀態(tài)固定、終端時刻自由，帶有狀態(tài)約束、控制約束的多階段多約束非線性最優(yōu)控制問題。最優(yōu)控制問題的求解方法一般可分為間接法和直接法兩種［1］。間接法根據(jù)極小值原理，將最優(yōu)控制問題轉(zhuǎn)化為兩點邊值問題，其解的精度高，但收斂半徑小，且協(xié)態(tài)變量由于沒有物理意義而很難準(zhǔn)確估計。直接法不需對協(xié)態(tài)變量初值進(jìn)行猜測，采用離散化方法將最優(yōu)控制問題轉(zhuǎn)化為非線性規(guī)劃問題(NLP)，其解的收斂半徑大，但解的最優(yōu)性不能保證。

作為直接法的一種，全局偽譜法近年來得到了廣泛的應(yīng)用，其通過較少的節(jié)點便可獲得較高精度的解。但在求解大型復(fù)雜問題及非光滑問題時，其運用效果并不理想。hp 自適應(yīng)偽譜法融合了全局偽譜法與有限元法的優(yōu)點，采用雙層策略對最優(yōu)控制問題進(jìn)行求解，相比于全局偽譜法其占用更少的計算時間卻可以獲得更高精度的解［2］。

MSLV 采用三級固體助推+液體末修的推進(jìn)方式，整個飛行過程分為5個階段:一級、二級、三級動力飛行段、無動力滑行段和末級飛行段。對hp 自適應(yīng)偽譜法在MSLV 主動段軌跡優(yōu)化中的應(yīng)用進(jìn)行了研究，仿真結(jié)果表明各級之間過渡平滑，終端入軌條件和過程約束得到很好滿足，hp 自適應(yīng)偽譜法具有很高的計算精度和較快的收斂速度，是一種行之有效的新型優(yōu)化算法。

1 主動段軌跡優(yōu)化模型

1.1 動力學(xué)模型

在地心慣性坐標(biāo)系OE－ XIYIZI下建立MSLV動力學(xué)模型，簡單且利于計算。地心慣性坐標(biāo)系定義［3］為:原點位于地心OE，OEXI軸于赤道平面內(nèi)指向春分點;OEZI軸垂直于赤道平面，與地球自轉(zhuǎn)角速度方向一致;OEYI軸由右手定則確定。MSLV 三自由度質(zhì)點動力學(xué)方程如式(1)所示，該方程基于下列前提或假設(shè):

1)地球是一個繞自身軸旋轉(zhuǎn)的均勻圓球;

2)推力方向總是與飛行器體軸方向重合;

3)MSLV 在飛行過程中以零側(cè)滑角飛行。

由上式可知，該動力學(xué)方程共包含7個狀態(tài)變量，其中r=［rxryrz］T，v =［vxvyvz］T分別為位置和速度矢量，m 為運載火箭質(zhì)量;此外，μ =GM為引力常數(shù)，T 為發(fā)動機(jī)推力，u=［uxuyuz］T為發(fā)動機(jī)推力方向，η 為發(fā)動機(jī)秒流量，D 為氣動阻力，其計算公式如下:

式(2)中，Cd為阻力系數(shù)，Sref為運載火箭參考面積，ρ 為大氣密度，其計算采用楊炳尉在“標(biāo)準(zhǔn)大氣參數(shù)的公式表示”一文中給出的擬合公式［4］，vref為相對運動速度，即:

ωe為地球相對慣性空間的旋轉(zhuǎn)角速度。

1.2 優(yōu)化模型

(1)性能指標(biāo)

在發(fā)動機(jī)性能參數(shù)一定的情況下，MSLV 主動段軌跡優(yōu)化的目的是為了提高運載能力，由于MSLV 前三級均為耗盡關(guān)機(jī)，因此選取末級助推消耗推進(jìn)劑質(zhì)量最小為優(yōu)化指標(biāo)，其可以等價描述為末級助推結(jié)束后MSLV 剩余質(zhì)量最大，即狀態(tài)變量中質(zhì)量的終端時刻值為最大，故指標(biāo)泛函選擇如下:

式中，tf為MSLV 主動段飛行結(jié)束的終端時刻。

(2)控制變量

選取發(fā)動機(jī)推力方向u =［uxuyuz］T及無動力滑行時間thx為待優(yōu)化的控制變量。

(3)約束條件控制變量約束:

終端入軌約束:

高度、動壓和過載等過程約束分別為:

狀態(tài)變量連接點約束:

式(5)～(9)中，tmin，tmax分別為無動力滑行時間的上下限，Hf，Vf，Θf分別為入軌點高度、速度和當(dāng)?shù)貜椀纼A角，af，ef，if分別為軌道的半長軸、偏心率和軌道傾角，可根據(jù)終端位置rf和速度vf求得，Re為地球半徑，qmax為動壓上限值，nmax為法向過載上限值，ms為分離時拋掉的結(jié)構(gòu)質(zhì)量。

2 hp 自適應(yīng)偽譜法

2.1 概述

近年來，在求解最優(yōu)控制問題時，直接法得到了越來越廣泛的應(yīng)用，文獻(xiàn)［5］從有限元法的思想出發(fā)，將直接法分為h 方法、p 方法和hp 方法3 種。

h 方法先將優(yōu)化問題進(jìn)行單元剖分，每個單元上使用低階插值基函數(shù)進(jìn)行離散，在優(yōu)化過程中保證基函數(shù)階次不變而通過逐步加密網(wǎng)格剖分對單元長度h 進(jìn)行細(xì)分(即h－細(xì)化)。常用的有限元法一般采用h 方法。采用h 方法離散后的NLP 一般都是稀疏的，計算效率高，但其為多項式級的收斂速度，在求解大型復(fù)雜問題時將帶來維數(shù)災(zāi)難問題。

p 方法正好相反，其將優(yōu)化問題作為一個單元(或少量單元)處理，每個單元上采用高階插值基函數(shù)進(jìn)行離散，優(yōu)化過程中保持單元剖分不變，通過增大基函數(shù)的階次p 以提高精度(即p－細(xì)化)。全局偽譜法就是一種典型的p 方法，其具有簡單的結(jié)構(gòu)、較高的精度及指數(shù)性的收斂速度，但對于非光滑問題，收斂速度非常慢，即使采用高階的插值基函數(shù)也可能無法得到滿足精度的解。

hp 方法允許單元長度h 和基函數(shù)階次p 在優(yōu)化過程中同時進(jìn)行自適應(yīng)調(diào)節(jié)。作為hp 方法的一種，hp 自適應(yīng)偽譜法將有限元法與全局偽譜法的優(yōu)點進(jìn)行融合，運用雙層策略決定單元數(shù)和基函數(shù)的階次以滿足精度要求。該方法包含兩部分內(nèi)容:hp－LG 離散方法和hp－網(wǎng)格細(xì)化方法。

2.2 hp－LG 離散方法

(1)時域變換

假設(shè)整個最優(yōu)控制問題在t∈［t0，tf］上被分成K個單元，即t0＜t1＜… ＜tK=tf。對于任意單元k，通過下式將時間區(qū)間由t∈［tk－1，tk］轉(zhuǎn)換到τ∈［－1，1］。

轉(zhuǎn)換后的性能指標(biāo)為:

系統(tǒng)微分方程變?yōu)?

邊界條件:

不等式約束:

此外，還存在內(nèi)點約束為:

(2)狀態(tài)變量與控制變量的近似

設(shè)x(k)(τ)和u(k)(τ)分別表示單元k 上的狀態(tài)變量和控制變量。在單元k 上以Nk個LG 點及τ0= －1為配點，構(gòu)造Nk+1 階Lagrange 插值多項式并以此為基函數(shù)近似狀態(tài)變量為:

其中

類似的，控制變量也可以近似如下:

(3)動力學(xué)微分方程約束轉(zhuǎn)換為代數(shù)約束

定義微分近似矩陣D(k)為:

系統(tǒng)微分方程動力學(xué)約束可進(jìn)一步轉(zhuǎn)換成一系列代數(shù)形式:

(4)離散條件下的終端狀態(tài)約束

由于在狀態(tài)變量的近似中沒有考慮終端時刻τf=1，因此將終端狀態(tài)約束離散并近似為:

(5)離散條件下的性能指標(biāo)

將Bolza 性能指標(biāo)中的Lagrange 積分項用Gauss 積分近似，即

這樣就將一個最優(yōu)控制問題轉(zhuǎn)化為參數(shù)優(yōu)化問題，選擇合適的非線性規(guī)劃算法即可進(jìn)行求解。

2.3 hp－網(wǎng)格細(xì)化方法

hp－網(wǎng)格細(xì)化方法的核心在于如何判定單元k上到底是進(jìn)行h－細(xì)化還是p －細(xì)化。對單元k 上的約束滿足程度進(jìn)行考核，在單元k 上任意選取L個節(jié)點，且滿足:

其中，i=1，…，n;j=1，…，s;l=1，…，L。

取e(k)的算術(shù)平均值為:

定義

定義誤差指標(biāo)ρ，下面分3 種情況進(jìn)行討論:

3 仿真計算及結(jié)果分析

3.1 仿真條件

目標(biāo)軌道選擇為700km 高度的太陽同步圓軌道，終端入軌約束條件如表1 所示。

表1 終端入軌約束

3.2 仿真結(jié)果分析

(1)仿真精度分析

hp 自適應(yīng)偽譜法得到的終端入軌參數(shù)與目標(biāo)值相比，偏差很小，詳見表2 所示，說明該算法具有非常高的計算精度。從圖1 和圖2 中可以看出，各階段之間銜接過渡平滑，高度和速度的終端約束得到了很好的滿足，三級與末級之間加入無動力滑行段，起到了提升高度和降低速度的作用。在無動力滑行段，由于發(fā)動機(jī)推力為零，推力方向不再是控制量，在滿足控制約束條件下可以任意取值，為方便起見，在優(yōu)化過程中將uz置為1，ux，uy置為0，因此控制變量在無動力滑行段出現(xiàn)了跳變。

表2 優(yōu)化軌道根數(shù)與目標(biāo)軌道根數(shù)的偏差

(2)配點情況分析

Gauss 偽譜法的配點為LG 點，該配點呈現(xiàn)出兩邊密中間疏的特點［6］，而hp 自適應(yīng)偽譜法在優(yōu)化過程中能夠?qū)卧獢?shù)和基函數(shù)的階次進(jìn)行自動調(diào)節(jié)，配點分布與Gauss 偽譜法有很大的不同。在無動力滑行段，MSLV 為純慣性飛行，運動相對簡單，配點數(shù)也較少，在其余的4個飛行階段，運動相對復(fù)雜，因此配點數(shù)也較多。從配點分布情況來看，配點不再呈現(xiàn)出兩邊密中間疏的特點，而是在梯度變化小的地方配點數(shù)較少，變化大的地方配點數(shù)相對較多，說明該方法具有良好的自適應(yīng)調(diào)節(jié)配點數(shù)的能力。

(3)最優(yōu)性分析

圖4 ～圖5 為狀態(tài)變量對應(yīng)的協(xié)態(tài)變量曲線，hp 自適應(yīng)偽譜法可以對協(xié)態(tài)變量特別是協(xié)態(tài)初值進(jìn)行準(zhǔn)確的估計。將協(xié)態(tài)變量初值與狀態(tài)變量初值一起代入狀態(tài)變量與協(xié)態(tài)變量組成的微分方程組，積分求解即可得到主動段最優(yōu)軌跡，將該軌跡與優(yōu)化得到的軌跡進(jìn)行比較可知兩者基本一致。圖6 為Hamilton 函數(shù)變化曲線，在每一個階段，Hamilton 函數(shù)都保持常值，這也從反面驗證了hp自適應(yīng)偽譜法與最優(yōu)控制理論中的一階必要條件是一致的。

圖1 高度隨時間變化曲線

圖2 速度隨時間變化曲線

圖3 控制變量隨時間變化曲線

圖4 協(xié)態(tài)變量λx，λy，λz隨時間變化曲線

4 結(jié)論

本文以MSLV 運載能力最大為優(yōu)化目標(biāo)，建立了基于hp 自適應(yīng)偽譜法的固體運載火箭主動段軌跡優(yōu)化模型，并進(jìn)行了仿真計算。仿真結(jié)果表明:

1)優(yōu)化后的軌跡各級之間連接過渡平滑，嚴(yán)格遵循過程約束條件，且終端入軌約束滿足精度高;

圖5 協(xié)態(tài)變量λvx，λvy，λvz隨時間變化曲線

圖6 Hamilton函數(shù)隨時間變化曲線

2)hp 自適應(yīng)偽譜法的配點分布能根據(jù)實際情況進(jìn)行自適應(yīng)調(diào)節(jié)，相比全局偽譜法兩端密中間疏的情況，該分布更為合理;

3)hp 自適應(yīng)偽譜法對協(xié)態(tài)變量的初值能夠進(jìn)行準(zhǔn)確的估計，可將該算法用于偽譜最優(yōu)反饋控制，實現(xiàn)對軌跡的實時跟蹤;

4)采用hp 自適應(yīng)偽譜法可降低對變量初始值的敏感性，具有快速、全局的優(yōu)點，具有良好的通用性和可擴(kuò)展性。

［1］Betts J T. Survey of Numerical Methods for Trajectory Optimization［J］. Journal of Guidance，Control and Dynamics，1998，21(2):193-207.

［2］Darby C L，Hager W W，Rao A V. An hp-Adaptive Pseudospectral Method for Solving Optimal Control Problems［J］. Optimal Control Applications and Methods，2011，32(4):476-502.

［3］龍樂豪，等.總體設(shè)計(上)［M］. 北京:宇航出版社，1989.

［4］楊炳尉. 標(biāo)準(zhǔn)大氣參數(shù)的公式表示［J］. 宇航學(xué)報，1983，(1):83-86.(Yang bingwei.Formulization of Standard Atmospheric Parameters［J］. Journal of Astronautics，1983，(1):83-86.)

［5］Darby C L，Hager W W，Rao A V.Direct Trajectory Optimization Using a Variable Low-Order Adaptive Pseudospectral Method［J］. Journal of Spacecraft and Rockets，2011，48(3):433-445.

［6］David Benson. A Gauss Pseudospectral Transcription for Optimal Control［D］. Cambridge，Massachusetts:Massachusetts Institute of Technology，2005.