基于hp自適應(yīng)偽譜法的組合動(dòng)力可重復(fù)使用運(yùn)載器軌跡優(yōu)化

周宏宇，王小剛，崔乃剛，郎 博

（1. 哈爾濱工業(yè)大學(xué) 航天學(xué)院，哈爾濱 150001；2. 國營第六二四廠，哈爾濱 150030）

基于hp自適應(yīng)偽譜法的組合動(dòng)力可重復(fù)使用運(yùn)載器軌跡優(yōu)化

周宏宇1，王小剛1，崔乃剛1，郎 博2

（1. 哈爾濱工業(yè)大學(xué) 航天學(xué)院，哈爾濱 150001；2. 國營第六二四廠，哈爾濱 150030）

針對(duì)傳統(tǒng)化學(xué)火箭難以重復(fù)使用、發(fā)射成本高的問題，提出了一種水平起飛/降落的新一代可重復(fù)使用運(yùn)載器飛行方案，并對(duì)其動(dòng)力模式設(shè)計(jì)和上升段軌跡優(yōu)化方法等關(guān)鍵技術(shù)進(jìn)行深入研究。設(shè)計(jì)了一種渦輪沖壓發(fā)動(dòng)機(jī)結(jié)合火箭沖壓發(fā)動(dòng)機(jī)的組合動(dòng)力模式，建立了發(fā)動(dòng)機(jī)推力與高度、馬赫數(shù)等變量間的耦合模型，根據(jù)動(dòng)力形式將上升段軌跡分為兩段并采用全局搜索法確定動(dòng)力切換的最佳時(shí)機(jī)。根據(jù)分段結(jié)果，分別以燃料最省和終端速度最大為指標(biāo)，利用 hp自適應(yīng)偽譜法對(duì)兩段進(jìn)行軌跡優(yōu)化設(shè)計(jì)。該算法基于雙層策略求解最優(yōu)控制問題，兼?zhèn)鋫巫V法和有限元法的優(yōu)點(diǎn)，與打靶法、偽譜法和間接法相比，初值更易選取，收斂速度更快。

hp自適應(yīng)偽譜法；組合動(dòng)力；可重復(fù)使用運(yùn)載器；軌跡優(yōu)化

載人航天的發(fā)展歷程可分為載人飛船和航天飛機(jī)兩個(gè)階段。載人飛船發(fā)射成本高，發(fā)射準(zhǔn)備周期長，航天員受到的過載很大，只能一次性使用；航天飛機(jī)的助推火箭是一次性的，維護(hù)成本高，危險(xiǎn)性高，缺少逃生系統(tǒng)。針對(duì)上述兩種運(yùn)載方式的不足，各航天強(qiáng)國紛紛提出新的載人航天策略[1]，著手設(shè)計(jì)可重復(fù)使用的新一代載人航天運(yùn)載器，其中，水平起飛/降落的方式是目前的研究熱點(diǎn)之一。

動(dòng)力系統(tǒng)是航天運(yùn)載器的關(guān)鍵部分，當(dāng)前飛行器的主要?jiǎng)恿ρb置有渦輪、沖壓和火箭等三種形式，它們適用于不同的飛行條件，單一的某種動(dòng)力形式不能同時(shí)實(shí)現(xiàn)可重復(fù)使用和載人航天任務(wù)這兩個(gè)目標(biāo)。文獻(xiàn)[2]研究了TBCC（渦輪基組合循環(huán)發(fā)動(dòng)機(jī)）動(dòng)力裝置并聯(lián)方案，完成了渦扇/沖壓發(fā)動(dòng)機(jī)的總體性能方案設(shè)計(jì)，但只適用于0～5 Ma和H=0～30 km的工作范圍；文獻(xiàn)[3]指出RBCC（火箭基組合循環(huán)發(fā)動(dòng)機(jī)）包含多個(gè)工作模態(tài)，具有工作速域?qū)挕⒖沼驈V的優(yōu)點(diǎn)，兼具加速和巡航能力，但它必須達(dá)到一定速度后才能啟用。TBCC和 RBCC發(fā)動(dòng)機(jī)有著各自的優(yōu)缺點(diǎn)和適用環(huán)境，若將兩者結(jié)合形成組合動(dòng)力發(fā)動(dòng)機(jī)，取長補(bǔ)短，便能夠在完成航天任務(wù)的同時(shí)實(shí)現(xiàn)運(yùn)載器的可重復(fù)使用。為了充分發(fā)揮組合動(dòng)力的優(yōu)勢，必須確定兩種動(dòng)力模式的最佳切換時(shí)機(jī)，并據(jù)此對(duì)不同飛行階段進(jìn)行軌跡優(yōu)化設(shè)計(jì)，從而獲得最大的運(yùn)載能力。

對(duì)于涉及復(fù)雜動(dòng)力形式和復(fù)雜高超聲速動(dòng)力學(xué)特性的優(yōu)化問題而言，打靶法、常規(guī)偽譜法和間接法等并不適用。文獻(xiàn)[4]通過簡化動(dòng)力學(xué)模型，利用間接法實(shí)現(xiàn)了高超聲速飛行器軌跡快速優(yōu)化，但發(fā)動(dòng)機(jī)模型過于簡化；文獻(xiàn)[5]針對(duì)含有復(fù)雜約束條件的非線性最優(yōu)控制問題，提出了改進(jìn)的高斯偽譜法，但存在初始不易選擇，為提高精度必須增加時(shí)間節(jié)點(diǎn)的缺點(diǎn)。

本文提出了一種TBCC配合RBCC使用的動(dòng)力模式，并建立了相應(yīng)的推力、比沖模型，使動(dòng)力系統(tǒng)在各種飛行狀態(tài)下都能獲得良好的動(dòng)力性能；然后依據(jù)動(dòng)力形式將飛行軌跡分段，針對(duì)不同階段的動(dòng)力特性、任務(wù)和環(huán)境因素選用不同的控制量、約束條件和性能指標(biāo)，采用hp自適應(yīng)偽譜法進(jìn)行分段軌跡優(yōu)化，并利用全局搜索法確定最佳動(dòng)力切換時(shí)機(jī)，提高優(yōu)化效率；最后進(jìn)行了仿真驗(yàn)證。

1 動(dòng)力系統(tǒng)建模

1.1 動(dòng)力模式

組合推進(jìn)技術(shù)不僅在理論上成立，而且已被工程實(shí)踐證明是合理可行的[3]。本文采用TBCC與RBCC相結(jié)合的動(dòng)力模式。TBCC實(shí)現(xiàn)渦輪/沖壓變循環(huán)工作過程，能在低速、低空環(huán)境下得到良好的推進(jìn)性能；RBCC將火箭和沖壓發(fā)動(dòng)機(jī)結(jié)合，在空氣稀薄處高速飛行時(shí)仍能獲得良好的推進(jìn)性能[6]。飛行器在 TBCC的作用下飛行，達(dá)到一定高度和速度后切換到RBCC，并最終完成上升段飛行，切換的高度和速度（或馬赫數(shù)）即為切換時(shí)機(jī)，分別記為Hqh和Vqh（或Mqh）。

1.2 TBCC發(fā)動(dòng)機(jī)模型

對(duì)于TBCC發(fā)動(dòng)機(jī)，認(rèn)為推力和比沖是飛行高度和飛行馬赫數(shù)的函數(shù)：

式中：P為發(fā)動(dòng)機(jī)推力，Is為發(fā)動(dòng)機(jī)比沖，H為飛行高度，Ma為飛行馬赫數(shù)。

為便于分析，參考常見的TBCC動(dòng)力特性和參數(shù)，通過多項(xiàng)式計(jì)算單位捕獲面積上的推力和比沖：

1.3 RBCC發(fā)動(dòng)機(jī)模型

對(duì)于RBCC發(fā)動(dòng)機(jī)，推力和比沖是飛行高度、飛行馬赫數(shù)和火箭發(fā)動(dòng)機(jī)秒耗量的函數(shù)。

火箭發(fā)動(dòng)機(jī)的秒耗量越大，沖壓發(fā)動(dòng)機(jī)所占比重就越小；若發(fā)動(dòng)機(jī)推力越大，則比沖就越小。為了便于分析，將RBCC發(fā)動(dòng)機(jī)單位捕獲面積上的推力和比沖看作是火箭發(fā)動(dòng)機(jī)所占比重的函數(shù)：

式中：Tmax為發(fā)動(dòng)機(jī)最大推力，Imax為發(fā)動(dòng)機(jī)最大比沖，K∈[0,1]表示節(jié)流閥（火箭發(fā)動(dòng)機(jī)所占比重）。這樣就可以通過控制火箭發(fā)動(dòng)機(jī)秒耗量（節(jié)流閥）來調(diào)節(jié)推力。本文取Tmax=360 kN，Imax=3000s。

計(jì)算出推力和比沖后，分別乘以發(fā)動(dòng)機(jī)捕獲面積SC即得到實(shí)際推力和比沖。

2 數(shù)學(xué)模型建立

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

為便于研究軌跡設(shè)計(jì)問題，假設(shè)飛行器無側(cè)向運(yùn)動(dòng)，不考慮地球扁率與旋轉(zhuǎn)，可得飛行器質(zhì)心運(yùn)動(dòng)方程如下：

式中：V為飛行速度大小，m為飛行器質(zhì)量，α為攻角，g為重力加速度大小，γ為飛行路徑角，R為飛行器到地心的距離，X和Y分別為阻力和升力。

式中：q為飛行動(dòng)壓，S為特征面積；cx和cy分別為阻力系數(shù)和升力系數(shù)。參考常見的高超聲速飛行器的氣動(dòng)特征，本文采用多項(xiàng)式計(jì)算氣動(dòng)系數(shù)：

2.2 優(yōu)化模型

一般的非線性系統(tǒng)狀態(tài)方程為如下所示[7]：

式中：x∈Rn為狀態(tài)變量，u∈Rm為控制變量。邊界條件為：

路徑約束或控制約束為：

考慮上面的一般非線性系統(tǒng)，優(yōu)化問題可表述為：在滿足狀態(tài)方程、邊界條件和路徑約束的條件下，求解控制變量u，使以下性能指標(biāo)達(dá)到極小：

式中：Φ為終端型指標(biāo)，g為積分型指標(biāo)。

3 上升段軌跡優(yōu)化

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

hp自適應(yīng)偽譜法將Radau偽譜法與hp型有限元法結(jié)合，在一系列離散點(diǎn)上利用全局多項(xiàng)式近似微分方程，將最優(yōu)控制問題的求解轉(zhuǎn)化為對(duì)非線性規(guī)劃（Nonlinear Programing, NLP）問題[8]的求解。

3.1.1 最優(yōu)控制問題轉(zhuǎn)化

將最優(yōu)控制問題的時(shí)間區(qū)間t∈[t0,tf]劃分為K個(gè)網(wǎng)格（子區(qū)間），并將每個(gè)時(shí)間區(qū)間t∈[tk-1,tk]轉(zhuǎn)換至τ∈[-1,1]，即

式中：k=1,2,…,K ；tk=tf。

狀態(tài)量在第k個(gè)網(wǎng)格上可近似表示為：

式中：X(k)(τ)為第k個(gè)網(wǎng)格上狀態(tài)量的近似值，表示為τ的函數(shù)；Nk為第k個(gè)網(wǎng)格上的 Lagrange-Gauss-Radau配點(diǎn)數(shù)；τj為第k個(gè)網(wǎng)格的節(jié)點(diǎn)，當(dāng)j∈[1,Nk]時(shí)，τj為配點(diǎn)，當(dāng)j=Nk+1時(shí)，τj為端點(diǎn)，表示終端時(shí)刻；Xj為τj處的狀態(tài)量的值，Lj為τj處的Lagrange多項(xiàng)式。

結(jié)合式(10)和式(15)可得：

只考慮終端型指標(biāo)，則目標(biāo)函數(shù)可近似表示為：

將不等式約束在第k個(gè)網(wǎng)格上用Nk個(gè)配點(diǎn)離散化處理：

邊界約束條件可近似表示為：

為了保證網(wǎng)格的連續(xù)性，需滿足：

上述離散過程將最優(yōu)控制問題轉(zhuǎn)化為非線性規(guī)劃問題，可以利用SNOPT[9]等軟件包進(jìn)行求解。

3.1.2 hp自適應(yīng)更新判定準(zhǔn)則

在每個(gè)網(wǎng)格內(nèi)設(shè)定最大容許誤差 εc，如果第k個(gè)網(wǎng)格的最大誤差，則需對(duì)該網(wǎng)格進(jìn)行重新細(xì)化。細(xì)化方法有兩種：1）增加網(wǎng)格內(nèi)插值多項(xiàng)式的階次p（即Nk）；2）在網(wǎng)格內(nèi)插入h個(gè)網(wǎng)點(diǎn)，將該網(wǎng)格細(xì)分為（h+1）個(gè)網(wǎng)格。

當(dāng)網(wǎng)格需重新細(xì)化時(shí)，首先要判定應(yīng)該增加p還是h。設(shè)網(wǎng)格內(nèi)的第i個(gè)狀態(tài)分量的曲率函數(shù)為：

3.1.3 更新變量計(jì)算公式

若增加p，則新的插值多項(xiàng)式的階次p1為

式中：p0為更新前的插值多項(xiàng)式階次，ceil()為正向舍入函數(shù)，A為自定義參數(shù)。

若增加h，則新增網(wǎng)點(diǎn)數(shù)h1為

式中：B為自定義參數(shù)。

3.1.4 hp自適應(yīng)偽譜法計(jì)算流程

hp自適應(yīng)偽譜法計(jì)算步驟如下所示，算法流程圖如圖1所示。具體流程為：

1）按照劃分好的網(wǎng)格，將最優(yōu)控制問題轉(zhuǎn)化為NLP問題，并對(duì)NLP問題進(jìn)行求解；

3）如果rk＜rkmax，將網(wǎng)格內(nèi)插值多項(xiàng)式的階次增加為p1，否則將網(wǎng)數(shù)增加為h1；

4）如果狀態(tài)量和路徑約束的相對(duì)誤差滿足要求，則優(yōu)化結(jié)束，否則返回步驟1）。

圖1 hp自適應(yīng)偽譜法流程圖Fig.1 Flow chart of hp-adaptive pseudo-spectral method

從以上過程可以看出，該算法較偽譜法而言，減少了多項(xiàng)式階次（配點(diǎn)）的增加，進(jìn)而減小了NLP問題的規(guī)模，有利于提高收斂速度；另外，多項(xiàng)式階次過高會(huì)使得輸出對(duì)輸入過于敏感，導(dǎo)致問題病態(tài)，該算法恰好可以避免這一點(diǎn)。

3.2 軌跡分段方法

TBCC發(fā)動(dòng)機(jī)特性如圖2所示，可以看出：在高度20 km以上時(shí)，推力隨馬赫數(shù)變化不大，且數(shù)值較小，因此Hqh不宜超過25 km；但Hqh過低，大氣密度會(huì)較大，導(dǎo)致動(dòng)壓很大，不適合RBCC工作，不利于后續(xù)爬升，故設(shè)定Hqh不低于10 km。

此外，TBCC的工作馬赫數(shù)一般不超過5。因此，本文將切換參數(shù)的范圍限定為：10≤Hqh≤25，2≤Mqh≤4。通過分析發(fā)動(dòng)機(jī)特性進(jìn)而縮小切換時(shí)機(jī)的選擇范圍，有利于提高計(jì)算速度和效率。

圖2 TBCC發(fā)動(dòng)機(jī)特性Fig.2 Feature of the engine with turbo-based combined cycle

4 仿真分析

4.1 仿真初始條件

采用 hp自適應(yīng)偽譜法對(duì)飛行器上升軌跡進(jìn)行優(yōu)化，并采用全局搜索法[10]在之前劃定的范圍內(nèi)確定發(fā)動(dòng)機(jī)切換狀態(tài)Hqh和Mqh。

飛行器特征面積63 m2，發(fā)動(dòng)機(jī)捕獲面積17 m2，TBCC段性能指標(biāo)為燃料消耗最少，控制變量為攻角α；RBCC段性能指標(biāo)為終端速度最大，控制變量為攻角α和發(fā)動(dòng)機(jī)相對(duì)流量K。飛行器初始狀態(tài)、終端狀態(tài)，以及搜索出的切換狀態(tài)如表1所示。

表1 初始及終端運(yùn)動(dòng)狀態(tài)設(shè)定Tab.1 Setting of the craft states

對(duì)于吸氣式高超聲速飛行器，為防止發(fā)動(dòng)機(jī)熄火，需要將攻角約束在一定范圍內(nèi)，以保證足夠的進(jìn)氣量；考慮氣動(dòng)熱約束及RBCC工作環(huán)境，需要對(duì)動(dòng)壓進(jìn)行約束；出于結(jié)構(gòu)強(qiáng)度的考慮，還需對(duì)過載進(jìn)行約束，具體約束條件見表2。

表2 過程約束設(shè)定Tab.2 Setting of the course constrains

4.2 仿真結(jié)果

TBCC段燃料消耗14 037 kg，飛行108.57 s；RBCC段結(jié)束時(shí)速度為4 951.53 m/s，飛行300 s。相關(guān)仿真結(jié)果圖像如圖3～9所示。

圖3 飛行攻角Fig.3 Time history of attack angle

圖4 速度變化情況Fig.4 Time history of velocity

圖5 高度變化情況Fig.5 Time history of height

圖6 飛行路徑角變化情況Fig.6 Time history of path angle

圖7 飛行剖面Fig.7 Flight profile

圖8 飛行過程變量Fig.8 Time history of course variety

圖9 相對(duì)流量變化情況Fig.9 Time history of relative flow

從仿真結(jié)果可以看出，優(yōu)化出的軌跡能夠滿足過程約束、控制約束和終端約束，且結(jié)果合理可行。

圖3中的攻角在動(dòng)力切換時(shí)發(fā)生了突變，這是由于攻角作為控制量在仿真中并未對(duì)其連續(xù)性進(jìn)行約束，而法向過載是與攻角有關(guān)的量，因此圖8中的法向過載也在同一時(shí)刻產(chǎn)生了突變。在實(shí)際應(yīng)用中，可以在TBCC段結(jié)束時(shí)，先通過姿控將飛行器攻角調(diào)整至RBCC段的起始攻角，然后再開機(jī)。

5 結(jié) 論

本文針對(duì)載人航天任務(wù)未來的發(fā)展趨勢，設(shè)計(jì)了一種組合動(dòng)力模式，并采用優(yōu)化方法進(jìn)行了軌跡優(yōu)化設(shè)計(jì)。具體如下：

1）對(duì)飛行器的動(dòng)力模式進(jìn)行了分析、設(shè)計(jì)，建立了與飛行環(huán)境相耦合的發(fā)動(dòng)機(jī)模型，更利于軌跡優(yōu)化設(shè)計(jì)。

2）為處理復(fù)雜優(yōu)化問題，采用了hp自適應(yīng)偽譜法，根據(jù)不同動(dòng)力形式將上升段飛行分段并結(jié)合各段的特點(diǎn)，分別選取不同的控制量、過程約束和性能指標(biāo)進(jìn)行分段優(yōu)化。

3）為提高優(yōu)化效率，盡可能提高性能指標(biāo)，關(guān)鍵在于確定各段之間的切換時(shí)機(jī)。通過分析飛行器動(dòng)力特點(diǎn)縮小了切換時(shí)機(jī)選取范圍，然后采用全局搜索法結(jié)合hp自適應(yīng)偽譜法的方法確定切換時(shí)機(jī)，并據(jù)此進(jìn)行分段優(yōu)化。

（References）：

[1] 孫龍, 劉映國. 2011年世界載人航天發(fā)展綜合分析[J].載人航天, 2012, 18(1): 92-96. Sun Long, Liu Ying-guo. A comprehensive analysis of the world manned spaceflight development in 2011[J]. Manned Spaceflight, 2012, 18(1): 92-96.

[2] Zhang H J, Liu X G, Guo R W. Design of turbo diffuser for TBCC inlet based on characteristics of turbo mode[J]. Journal of Aerospace Power, 2014(1): 181-191.

[3] 何國強(qiáng), 秦飛, 魏祥庚, 等. 火箭沖壓組合發(fā)動(dòng)機(jī)燃燒的若干基礎(chǔ)問題研究[J]. 實(shí)驗(yàn)流體力學(xué), 2016, 30(1): 1-14. He Guo-qiang, Qin Fei, Wei Xiang-geng, et al. Investigation of several fundamental combustion problems in rocket based combined-cycle engines[J]. Journal of Experiments in Fluid Mechanics, 2016, 30(1): 1-14.

[4] Murillo O. Fast ascent trajectory optimization for hypersonic air-breathing vehicles[C]//AIAA Guidance, Navigation, and Control Conference. 2010: 14-38.

[5] 孫勇. 基于改進(jìn) Gauss偽譜法的高超聲速飛行器軌跡優(yōu)化與制導(dǎo)[D]. 哈爾濱工業(yè)大學(xué), 2012: 42-46. Sun Yong. Trajectory optimization and guidance of hypersonic vehicle based on improved gauss pseudospectral method[D]. Harbin, China: Harbin Institute of Technology, 2012: 42-46.

[6] Gong C, Chen B, Gu L. Design and optimization of RBCC powered suborbital reusable launch vehicle[C]// AIAA International Space Planes and Hypersonic Systems and Technologies Conference. 2014: 26-47.

[7] Song C, Zhao G, Chen J. Reentry trajectory optimization for hypersonic vehicle based on GPM[C]//Proceedings of the 2012 Second International Conference on Electric Information and Control Engineering. 2012: 1183-1186.

[8] Tohidi E, Lotfi N S. An efficient legendre pseudospectral method for solving nonlinear quasi bang-bang optimal control problems[J]. Journal of Applied Mathematics, Statistics and Informatics, 2012, 8(2): 73-85.

[9] Conte J P. OpenSees-SNOPT framework for finite elementbased optimization of structural and geotechnical systems[J]. Journal of Structural Engineering, 2012, 137 (6): 822-834.

[10] Crouser R J, Schmidt M C, Kelley S, et al. Global pattern search at scale[C]//IEEE International Symposium on Technologies for Homeland Security. IEEE, 2015: 1-6.

Trajectory optimization of reusable vehicle with combined power based on hp adaptive pseudospectral algorithm

ZHOU Hong-yu1, WANG Xiao-gang1, CUI Nai-gang1, LANG Bo2
(1. Department of Astronautics, Harbin Institute of Technology, Harbin 150001, China; 2. Stale-owned Factory No.624, Harbin 150030, China)

In view that traditional chemical rocket is hard to be reused and suffers from high launch cost, a novel flight scheme based on horizontally taking off and landing is proposed, by which a vehicle can be completely reusable. The power mode and the optimization algorithm of the ascend trajectory are designed. A combined power mode of turbo-based combined cycle and rocket-based combined cycle is presented, and the engine model is built based on the relationship between thrust and height/velocity. Then the ascend trajectory is divided into two phases depending on the power mode. A universal searching method is employed to confirm the best opportunity in changing the engine mode, and an hp-adaptive pseudospectral method (HPM) is developed, which takes into account the influences of various constrains and value functions on the results of the universal searching method in different phases. The HPM adopts dual-layer strategy to solve the optimal control problem, and has both the merits of the finite element method and the pseudo-spectral method. The HPM can choose the initial value more easily, and has even faster convergence speed, compared with the shooting method, the pseudo- spectral method and the indirect method.

hp-adaptive pseudospectral; combined power; reusable vehicle; trajectory optimization

V412.4

：A

2016-09-02；

：2016-11-30

國家自然科學(xué)基金（61304236）

周宏宇（1991—），男，博士研究生，研究方向?yàn)轱w行器制導(dǎo)與優(yōu)化。E-mail: 826867857@qq.com

聯(lián) 系 人：王小剛（1980—），男，副教授，研究方向?yàn)榉蔷€性濾波理論及應(yīng)用。E-mail: wangxiaogang@hit.edu.cn

1005-6734(2016)06-0832-06

10.13695/j.cnki.12-1222/o3.2016.06.024