洪 蓓 辛萬青

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

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

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

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

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

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

1.1 動(dòng)力學(xué)模型

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

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

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

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

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

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

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

1.2 優(yōu)化模型

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

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

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

(2)控制變量

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

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

終端入軌約束:

高度、動(dòng)壓和過載等過程約束分別為:

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

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

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

2.1 概述

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

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

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

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

2.2 hp－LG 離散方法

(1)時(shí)域變換

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

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

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

邊界條件:

不等式約束:

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

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

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

其中

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

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

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

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

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

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

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

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

這樣就將一個(gè)最優(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ì)化。對(duì)單元k 上的約束滿足程度進(jìn)行考核，在單元k 上任意選取L個(gè)節(jié)點(diǎn)，且滿足:

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

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

定義

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

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

3.1 仿真條件

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

表1 終端入軌約束

3.2 仿真結(jié)果分析

(1)仿真精度分析

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

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

(2)配點(diǎn)情況分析

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

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

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

圖1 高度隨時(shí)間變化曲線

圖2 速度隨時(shí)間變化曲線

圖3 控制變量隨時(shí)間變化曲線

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

4 結(jié)論

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

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

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

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

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

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

4)采用hp 自適應(yīng)偽譜法可降低對(duì)變量初始值的敏感性，具有快速、全局的優(yōu)點(diǎn)，具有良好的通用性和可擴(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è)計(jì)(上)［M］. 北京:宇航出版社，1989.

［4］楊炳尉. 標(biāo)準(zhǔn)大氣參數(shù)的公式表示［J］. 宇航學(xué)報(bào)，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.