閆循良,王舒眉,王培臣,劉海禮

(1.西北工業(yè)大學(xué) 航天學(xué)院,陜西 西安 710072; 2.中國(guó)商飛民用飛機(jī)試飛中心,上海 201323)

近年來(lái),高超聲速飛行器的多樣化發(fā)展和復(fù)雜任務(wù)需求對(duì)助推上升段動(dòng)力系統(tǒng)的設(shè)計(jì)提出了許多新的要求,如較強(qiáng)的短時(shí)加速、全程多次開(kāi)關(guān)機(jī)、寬域工作能力等。與傳統(tǒng)火箭動(dòng)力相比,火箭基組合循環(huán)(rocket based combined cycle,RBCC)動(dòng)力系統(tǒng)將高推重比、低比沖的火箭發(fā)動(dòng)機(jī)和低推重比、高比沖的沖壓發(fā)動(dòng)機(jī)有機(jī)地組合在一起,具有低成本、技術(shù)先進(jìn)、使用靈活等特征,可用于執(zhí)行廣空域、寬速域運(yùn)載及新質(zhì)作戰(zhàn)等飛行任務(wù)[1-2]。

然而,RBCC動(dòng)力系統(tǒng)的引入給高超聲速飛行器的上升段軌跡設(shè)計(jì)帶來(lái)了新的挑戰(zhàn)[3-4]。首先,RBCC動(dòng)力系統(tǒng)工作模態(tài)多,各模態(tài)性能迥異、動(dòng)態(tài)切換;其次,上升段狀態(tài)參數(shù)變化范圍大,動(dòng)力性能與飛行狀態(tài)耦合程度高,且面臨多種復(fù)雜的過(guò)程及終端約束限制。上述因素使得上升段軌跡設(shè)計(jì)的可行域極小、設(shè)計(jì)難度極大,傳統(tǒng)火箭助推運(yùn)載器的軌跡設(shè)計(jì)方法和經(jīng)驗(yàn)規(guī)律已難以勝任,有必要針對(duì)其動(dòng)力系統(tǒng)特點(diǎn)開(kāi)展相應(yīng)的數(shù)值優(yōu)化求解算法及策略研究。

目前,典型的上升段軌跡優(yōu)化方法可分為間接法[5]和直接法[6-9]。其中,間接法由于模型精度低、通用性不高且收斂性較差,難以應(yīng)用于具有復(fù)雜、高保真模型的RBCC上升軌跡優(yōu)化問(wèn)題。隨著計(jì)算機(jī)技術(shù)的發(fā)展,直接法中的配點(diǎn)法、偽譜法以及粒子群等數(shù)值優(yōu)化方法在上升段軌跡優(yōu)化設(shè)計(jì)中得到廣泛應(yīng)用[6-9],然而,上述方法均未能滿足上升段軌跡快速優(yōu)化的需求。因此,為提升軌跡優(yōu)化設(shè)計(jì)的計(jì)算效率,部分學(xué)者將凸優(yōu)化方法[10-14]引入上升段軌跡設(shè)計(jì)中。Szmuk等[12]以推力方向作為控制量,利用無(wú)損凸化技術(shù)求解運(yùn)載火箭上升段軌跡優(yōu)化問(wèn)題,但其忽略了運(yùn)載器所受升力。Liu等[13]將推力方向和氣動(dòng)系數(shù)作為控制量,但考慮的氣動(dòng)力模型較為簡(jiǎn)單。王嘉煒等[14]在對(duì)固體火箭助推飛行器軌跡優(yōu)化問(wèn)題進(jìn)行建模時(shí),選擇攻角作為唯一控制量,有效處理了飛行器氣動(dòng)力的計(jì)算。雖然基于凸優(yōu)化的方法已逐步用于解決上升段軌跡設(shè)計(jì)問(wèn)題,然而現(xiàn)有公開(kāi)文獻(xiàn)較少涉及RBCC動(dòng)力上升段軌跡優(yōu)化問(wèn)題。

因此,本文以RBCC動(dòng)力高超聲速飛行器為研究對(duì)象,提出了一種考慮復(fù)雜約束、高非線性及強(qiáng)耦合模型限制的動(dòng)力上升段軌跡快速優(yōu)化方法。針對(duì)攻角控制系統(tǒng)是否存在二階滯后情況,構(gòu)建了完整的RBCC動(dòng)力上升段軌跡優(yōu)化模型;設(shè)計(jì)了基于序列凸優(yōu)化的上升段軌跡快速優(yōu)化策略及算法。仿真結(jié)果表明,該方法在保證可行性的同時(shí),突破了傳統(tǒng)軌跡優(yōu)化方法求解該類復(fù)雜約束問(wèn)題時(shí)效率低或無(wú)法找到可行解的局限,能夠提高RBCC動(dòng)力飛行器的軌跡優(yōu)化設(shè)計(jì)效率。

1 上升段軌跡優(yōu)化問(wèn)題描述

1.1 上升段運(yùn)動(dòng)模型

考慮飛行器的上升段運(yùn)動(dòng)主要保持在縱向平面內(nèi),故不考慮側(cè)向運(yùn)動(dòng)。假設(shè)地球?yàn)椴恍D(zhuǎn)圓球,建立縱向平面內(nèi)上升段質(zhì)心運(yùn)動(dòng)方程

(1)

式中:r,V,θ,m,α,g,ms,P,D,L分別為地心距、速度、當(dāng)?shù)貜椀纼A角、飛行器質(zhì)量、攻角、當(dāng)?shù)刂亓铀俣?、發(fā)動(dòng)機(jī)燃料秒流量、推力、氣動(dòng)阻力及升力。對(duì)于上述模型而言,其狀態(tài)量x=[r,V,θ,m]T,控制量u=[α,ms]T。

若考慮攻角控制存在二階滯后,則攻角實(shí)際值與指令值并不完全一致。此時(shí),補(bǔ)充如下方程

(2)

氣動(dòng)力D,L可由(3)式計(jì)算得到

(3)

式中:q=0.5ρV2為飛行動(dòng)壓;ρ為大氣密度;Sref為飛行器氣動(dòng)參考面積;升、阻力系數(shù)CL,CD均為高度、攻角和馬赫數(shù)的函數(shù),可通過(guò)插值得到。

本文研究的高超聲速飛行器采用RBCC動(dòng)力系統(tǒng)助推,該系統(tǒng)由火箭發(fā)動(dòng)機(jī)與沖壓發(fā)動(dòng)機(jī)分系統(tǒng)組合而成,因此,其推力與燃料秒流量可表示為

(4)

式中:msr和msa分別為火箭發(fā)動(dòng)機(jī)和沖壓發(fā)動(dòng)機(jī)的燃料秒流量;PR和PA分別代表火箭發(fā)動(dòng)機(jī)與沖壓發(fā)動(dòng)機(jī)推力。

火箭發(fā)動(dòng)機(jī)推力PR計(jì)算如(5)式所示

(5)

式中,Ispr為火箭發(fā)動(dòng)機(jī)比沖,可表征為關(guān)于高度的函數(shù)fIspr(h)。

沖壓發(fā)動(dòng)機(jī)推力計(jì)算為

(6)

式中,S,Er,φ,Ispa分別為沖壓發(fā)動(dòng)機(jī)進(jìn)氣道橫截面積、燃油流量比、流量系數(shù)及比沖。可見(jiàn),沖壓推力與高度、馬赫數(shù)、攻角及發(fā)動(dòng)機(jī)參數(shù)均相關(guān),即動(dòng)力性能與飛行狀態(tài)存在高度的非線性和強(qiáng)耦合關(guān)系。

1.2 約束條件描述

(7)

從實(shí)際可行性以及安全飛行角度考慮,狀態(tài)變量亦需要滿足一定的約束,即有

xmin≤x≤xmax

(8)

綜合考慮飛行性能、控制系統(tǒng)性能和發(fā)動(dòng)機(jī)性能,可建立控制變量約束模型,有

umin≤u≤umax

(9)

此外,對(duì)于本文研究的問(wèn)題而言,端點(diǎn)約束可分為初始約束與終端約束,即有

(10)

1.3 優(yōu)化問(wèn)題描述

為挖掘飛行器運(yùn)載潛力,將優(yōu)化目標(biāo)選取為上升段末端機(jī)械能最大,優(yōu)化目標(biāo)可以表示為

J=[0.5m·V2+m·g·h]|tf

(11)

綜上,RBCC上升段軌跡優(yōu)化問(wèn)題可以描述為:尋找最佳狀態(tài)量x(t)及控制量u(t),使得目標(biāo)函數(shù)(11)式最大,同時(shí)滿足狀態(tài)方程約束(1)或(2)式,以及各類約束(7)～(10)式。

為了提升優(yōu)化求解的收斂性,往往需要對(duì)模型進(jìn)行無(wú)量綱化處理,限于篇幅,此處不再贅述。

2 基于凸優(yōu)化的求解框架及策略設(shè)計(jì)

典型凸優(yōu)化問(wèn)題即為尋找最優(yōu)控制量,使得

(12)

2.1 時(shí)間自由問(wèn)題轉(zhuǎn)化處理

對(duì)于上升段總飛行時(shí)間未知的這類時(shí)間自由問(wèn)題,首先定義新的自變量τ∈[0,1]和控制量ut=tf-t0,并且將原問(wèn)題的時(shí)間區(qū)間映射到[0,1]上,得到

t=t0+(tf-t0)τ,τ∈[0,1]

(13)

同時(shí),將原自變量時(shí)間作為新的狀態(tài)變量,則有

(14)

(15)

(16)

2.2 凸化處理

考慮到前述上升段軌跡優(yōu)化模型是非凸的,將其進(jìn)行凸化處理是應(yīng)用凸優(yōu)化技術(shù)數(shù)值求解的關(guān)鍵。本節(jié)以姿態(tài)控制理想情況下上升段軌跡優(yōu)化問(wèn)題為例采用逐次線性化等[11]技術(shù)進(jìn)行模型凸化。

1) 動(dòng)力學(xué)方程凸化

(17)

取其作為初始參數(shù)序列,將動(dòng)力學(xué)方程(16)在該參考狀態(tài)序列下進(jìn)行一階泰勒展開(kāi),得到線性化的動(dòng)態(tài)方程約束,即

(18)

同時(shí),為保證線性化的有效性和精度,需引入信賴域約束:

(19)

(20)

其動(dòng)態(tài)方程約束可類比得到,此處不做贅述。

2) 過(guò)程約束處理

(21)

對(duì)其進(jìn)行線性化可得

(22)

3) 指標(biāo)函數(shù)凸化處理

記機(jī)械能公式

J=[0.5m·V2+m·g·h]|tf=f10

(23)

為非凸形式。同理,采用逐次線性化方法,可將其轉(zhuǎn)化為

(24)

此外,注意到狀態(tài)量和控制量表達(dá)(8)～(10)式均為凸形式,故無(wú)需做凸化處理。至此,該上升段軌跡優(yōu)化即轉(zhuǎn)化為一典型的凸問(wèn)題。

2.3 離散化處理

采用梯形法則[11]對(duì)(18)式進(jìn)行離散,可得

(25)

上標(biāo)k表示第k次迭代,m表示離散點(diǎn)編號(hào),m=0,1,…,M;Δτ=(τf-τ0)/M。

類似地,(22)式的離散化形式為

(26)

控制變量約束、狀態(tài)變量約束的離散形式為

(27)

狀態(tài)變量的端點(diǎn)條件可以進(jìn)一步表述為

(28)

綜上所述,原軌跡優(yōu)化問(wèn)題轉(zhuǎn)化為參數(shù)化凸問(wèn)題,可表示為如下形式

(29)

2.4 凸優(yōu)化求解策略改進(jìn)及流程

1) 終端約束處理

考慮終端高度等端點(diǎn)等式約束,若在優(yōu)化求解過(guò)程中直接限制為固定值,初期迭代過(guò)程中很難滿足該約束,甚至?xí)霈F(xiàn)不可行解進(jìn)而導(dǎo)致優(yōu)化失敗。因此,可將這類約束轉(zhuǎn)化為性能指標(biāo)中的懲罰項(xiàng)[11]

(30)

式中:常系數(shù)ci>0,i為終端等式約束序號(hào);I為約束總數(shù)量。考慮該懲罰項(xiàng)為非凸形式,故通過(guò)引入松弛矢量R,將其轉(zhuǎn)化為(31)式的形式引入性能指標(biāo)中

p(R)=c·R

(31)

對(duì)應(yīng)的松弛約束為

(32)

其他終端等式約束即可采用該方法進(jìn)行處理。

于是,上述參數(shù)化凸問(wèn)題可松弛為

(33)

2) 優(yōu)化求解流程

考慮到上述凸化后的問(wèn)題與原軌跡優(yōu)化問(wèn)題之間存在偏差,通過(guò)序列凸優(yōu)化算法[15]迭代求解上述凸問(wèn)題可逐步逼近原問(wèn)題的解。

為保證線性化的有效性和解的收斂性及精度,迭代過(guò)程以當(dāng)前優(yōu)化軌跡不斷更新參考軌跡,并以相鄰2次迭代的凸優(yōu)化解對(duì)應(yīng)的狀態(tài)量最大偏差作為收斂準(zhǔn)則,即

(34)

式中,ε為收斂誤差限。

圖1 基于序列凸優(yōu)化的上升段軌跡優(yōu)化求解流程圖

3) 判斷是否滿足(34)式和回溯搜索條件,若滿足,則優(yōu)化結(jié)束;否則,令k=k+1,轉(zhuǎn)入步驟2);

3 算例仿真及分析

3.1 仿真條件

綜合考慮計(jì)算效率與解的精度,本文仿真過(guò)程中共取61個(gè)離散點(diǎn),序列凸優(yōu)化算法的初始信賴域和收斂誤差限設(shè)置為

此外,回溯搜索算法中的2個(gè)參數(shù)均取0.8[15]。

所有仿真均在搭載Intel Core i7-8700 3.20 GHz Intel處理器的臺(tái)式機(jī)完成,仿真環(huán)境為MATLAB 2016b平臺(tái),基于CVX工具包進(jìn)行軌跡優(yōu)化算法開(kāi)發(fā),并調(diào)用SDPT3求解器求解凸優(yōu)化子問(wèn)題。

3.2 攻角控制系統(tǒng)理想情況下的上升段軌跡優(yōu)化

將控制系統(tǒng)理想情況下軌跡優(yōu)化模型簡(jiǎn)記為模型一,針對(duì)該模型進(jìn)行終端機(jī)械能最大的上升段軌跡優(yōu)化。凸優(yōu)化算法經(jīng)過(guò)12輪迭代后收斂,每輪迭代約耗時(shí)1.8 s,優(yōu)化過(guò)程共耗時(shí)19.88 s,優(yōu)化結(jié)果如圖2～3所示。

圖2 模型一對(duì)應(yīng)的部分狀態(tài)量?jī)?yōu)化結(jié)果

圖3 模型一對(duì)應(yīng)的控制量?jī)?yōu)化結(jié)果

分析仿真結(jié)果可知,上升段共飛行135.66 s,終端質(zhì)量為800 kg,終端高度為62.08 km,速度為4 866.99 m/s,滿足終端約束;動(dòng)壓、過(guò)載、駐點(diǎn)熱流密度均滿足給定的過(guò)程約束限制。將所得控制量進(jìn)行插值后代入運(yùn)動(dòng)方程進(jìn)行積分,所得積分結(jié)果與優(yōu)化結(jié)果對(duì)應(yīng)的狀態(tài)相對(duì)偏差小于0.1%,驗(yàn)證了優(yōu)化結(jié)果的可行性,表明本文所提序列凸優(yōu)化方法可逼近原軌跡優(yōu)化問(wèn)題的解。圖2～3中虛線為凸優(yōu)化所使用的初始猜測(cè)軌跡,本文僅將初始猜測(cè)軌跡取為初、末端狀態(tài)猜測(cè)值的連線,這種情況下,凸優(yōu)化方法仍能有效、快速地收斂至優(yōu)化結(jié)果?？梢?jiàn),本文所設(shè)計(jì)的優(yōu)化算法和求解策略可以快速、有效地求解RBCC高超聲速上升這一復(fù)雜、非線性、強(qiáng)耦合的非凸最優(yōu)控制問(wèn)題。

此外,結(jié)合仿真結(jié)果可知,對(duì)于該問(wèn)題而言,上升段飛行主要可分為3段。第一段為0～22 s,此時(shí)飛行馬赫數(shù)較低,因此沖壓發(fā)動(dòng)機(jī)推力小、效率低,RBCC動(dòng)力系統(tǒng)工作在沖壓+火箭混合模態(tài),且以火箭加速為主,此時(shí)馬赫數(shù)逐漸增加,但飛行高度變化較小。第二段為22～91 s,該段飛行高度、馬赫數(shù)都較適宜沖壓發(fā)動(dòng)機(jī)工作,故此時(shí)火箭發(fā)動(dòng)機(jī)關(guān)機(jī),RBCC動(dòng)力系統(tǒng)工作在沖壓模態(tài),飛行馬赫數(shù)和高度均緩慢變化,直至爬升至約30 km,飛行器在該段將達(dá)到動(dòng)壓約束邊界并保持約17 s。第三段對(duì)應(yīng)為91 s至上升段結(jié)束,此時(shí)由于高度較高,沖壓發(fā)動(dòng)機(jī)效率逐步降低,RBCC動(dòng)力系統(tǒng)先工作在沖壓+火箭混合模態(tài),后逐步過(guò)渡到純火箭模態(tài)直至燃料耗盡,此時(shí)飛行高度和馬赫數(shù)快速增加。事實(shí)上,這種火箭動(dòng)力峰-谷-峰的工作過(guò)程能夠最大限度地發(fā)揮火箭發(fā)動(dòng)機(jī)工作范圍廣和沖壓發(fā)動(dòng)機(jī)推進(jìn)效率高的優(yōu)勢(shì),從而使飛行器能夠在燃料一定的情況下獲得更大終端高度、馬赫數(shù),即實(shí)現(xiàn)終端機(jī)械能最大。

此外,為探究離散點(diǎn)數(shù)選取對(duì)凸優(yōu)化算法精度及計(jì)算效率的影響,取31和91個(gè)離散點(diǎn)分別進(jìn)行軌跡優(yōu)化,并利用所得優(yōu)化控制量進(jìn)行積分,所得優(yōu)化及積分結(jié)果對(duì)比如表1所示。

表1 不同離散點(diǎn)數(shù)量?jī)?yōu)化及積分結(jié)果對(duì)比

由表1可知,3組優(yōu)化結(jié)果的迭代次數(shù)基本一致,而隨著離散點(diǎn)數(shù)的增加,計(jì)算耗時(shí)有所增加,優(yōu)化解與積分解之間的誤差顯著降低,表明優(yōu)化解的可行性逐漸提升。相較于離散點(diǎn)取61與91的情況,離散點(diǎn)取31時(shí),終端狀態(tài)的優(yōu)化解與積分解誤差較大,表明該組優(yōu)化結(jié)果不可行;而相較于離散點(diǎn)91的情況,離散點(diǎn)取61的優(yōu)化結(jié)果在精度基本相當(dāng)情況下,計(jì)算耗時(shí)減少了30%。因此,綜合考慮計(jì)算效率和精度,取離散點(diǎn)為61較為合適。

3.3 考慮攻角控制系統(tǒng)二階滯后上升段軌跡優(yōu)化

將考慮控制系統(tǒng)二階滯后的軌跡優(yōu)化模型簡(jiǎn)記為模型二,針對(duì)該模型進(jìn)行上升段軌跡優(yōu)化。仿真中,取阻尼比ξ=0.7,固有頻率ωn=0.5。凸優(yōu)化算法經(jīng)過(guò)13輪迭代后收斂,每輪迭代約耗時(shí)1.8 s,優(yōu)化過(guò)程共耗時(shí)23.08 s,優(yōu)化結(jié)果如圖4～5所示。

圖4 考慮控制滯后的狀態(tài)量?jī)?yōu)化結(jié)果

由結(jié)果可知,上升段共飛行136.20 s,剩余質(zhì)量為800 kg,終端高度為62.18 km、飛行馬赫數(shù)為15.59,滿足終端約束,動(dòng)壓、過(guò)載、駐點(diǎn)熱流密度均滿足給定的過(guò)程約束限制。

由圖4可知,模型一、二的優(yōu)化結(jié)果較為接近,而由圖4d)的優(yōu)化曲線可知,模型二的彈道傾角變化存在一定滯后,尤其是在飛行前中期的幾個(gè)拐點(diǎn)位置,考慮是控制系統(tǒng)滯后導(dǎo)致實(shí)際攻角變化滯后于模型一,這一點(diǎn)亦可從圖4e)的攻角指令曲線得以印證。

由圖5a)可知,模型二的攻角變化率總是提前模型一變化,考慮是模型二控制系統(tǒng)存在二階滯后,需要提前改變控制指令,以補(bǔ)償滯后效應(yīng),使實(shí)際攻角曲線盡可能逼近模型一的最優(yōu)結(jié)果。此外,模型一終端機(jī)械能為952 380 kJ,模型二終端機(jī)械能為952 310 kJ,略低于模型一?？梢?jiàn),模型二的控制滯后犧牲了一定的最優(yōu)性。

圖5 考慮控制滯后的控制量?jī)?yōu)化結(jié)果

為進(jìn)一步探究自然頻率ωn取值對(duì)軌跡的影響,分別取ωn=0.1,ωn=0.9進(jìn)行仿真,結(jié)果如圖6及表2所示。由表2可知,相較于理想情況,考慮控制系統(tǒng)二階滯后特性軌跡優(yōu)化結(jié)果的最優(yōu)性有所降低,且最優(yōu)性隨著ωn的增大而提升。由圖6可知,當(dāng)ωn較小時(shí),攻角曲線較為光滑,此時(shí)實(shí)際攻角對(duì)攻角指令的響應(yīng)速度較慢,攻角變化較指令值滯后情況明顯,因此需要提前給出指令信號(hào)以更好地操縱飛行。當(dāng)ωn較大時(shí),實(shí)際攻角對(duì)攻角指令的跟蹤效果更好,此時(shí)飛行器具有較好的操縱性。通常來(lái)說(shuō),飛行器自然頻率由其靜穩(wěn)定性決定,但設(shè)計(jì)控制系統(tǒng)時(shí),一般情況下要求自然頻率低于自動(dòng)駕駛儀頻率以免出現(xiàn)共振。同時(shí),由于控制系統(tǒng)的超調(diào)量及振蕩次數(shù)均只與阻尼比ξ相關(guān),而實(shí)際攻角的最大峰值并未超過(guò)指令攻角,因此設(shè)計(jì)姿態(tài)控制系統(tǒng)時(shí)應(yīng)依據(jù)實(shí)際需要進(jìn)行控制參數(shù)選擇。

表2 不同ωn取值優(yōu)化結(jié)果

圖6 不同ωn取值優(yōu)化結(jié)果攻角對(duì)比

4 結(jié) 論

針對(duì)RBCC高超聲速上升段軌跡快速優(yōu)化問(wèn)題,本文設(shè)計(jì)了一種基于序列凸優(yōu)化的軌跡優(yōu)化算法和策略,相關(guān)結(jié)論如下:

1) 所設(shè)計(jì)優(yōu)化算法可以快速、有效地求解復(fù)雜工作模態(tài)下RBCC高超聲速上升這一復(fù)雜、非線性、強(qiáng)耦合的非凸最優(yōu)控制問(wèn)題;

2) RBCC發(fā)動(dòng)機(jī)火箭模態(tài)峰-谷-峰的工作過(guò)程能夠最大限度地結(jié)合并發(fā)揮火箭發(fā)動(dòng)機(jī)工作范圍廣和沖壓發(fā)動(dòng)機(jī)經(jīng)濟(jì)性高的優(yōu)勢(shì),從而使飛行器在燃料更省的情況下獲得更大的終端高度、馬赫數(shù)。

3) 相較于理想情況,考慮控制系統(tǒng)二階滯后特性軌跡優(yōu)化結(jié)果的最優(yōu)性有所降低,且降低程度與自然頻率ωn有關(guān),設(shè)計(jì)姿態(tài)控制系統(tǒng)時(shí)應(yīng)依據(jù)實(shí)際需要進(jìn)行控制參數(shù)選擇。