RBCC高超聲速飛行器上升段軌跡快速優化

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

(1.西北工業大學 航天學院,陜西 西安 710072; 2.中國商飛民用飛機試飛中心,上海 201323)

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

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

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

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

1 上升段軌跡優化問題描述

1.1 上升段運動模型

考慮飛行器的上升段運動主要保持在縱向平面內,故不考慮側向運動。假設地球為不旋轉圓球,建立縱向平面內上升段質心運動方程

(1)

式中:r,V,θ,m,α,g,ms,P,D,L分別為地心距、速度、當地彈道傾角、飛行器質量、攻角、當地重力加速度、發動機燃料秒流量、推力、氣動阻力及升力。對于上述模型而言,其狀態量x=[r,V,θ,m]T,控制量u=[α,ms]T。

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

(2)

氣動力D,L可由(3)式計算得到

(3)

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

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

(4)

式中:msr和msa分別為火箭發動機和沖壓發動機的燃料秒流量;PR和PA分別代表火箭發動機與沖壓發動機推力。

火箭發動機推力PR計算如(5)式所示

(5)

式中,Ispr為火箭發動機比沖,可表征為關于高度的函數fIspr(h)。

沖壓發動機推力計算為

(6)

式中,S,Er,φ,Ispa分別為沖壓發動機進氣道橫截面積、燃油流量比、流量系數及比沖?？梢?沖壓推力與高度、馬赫數、攻角及發動機參數均相關,即動力性能與飛行狀態存在高度的非線性和強耦合關系。

1.2 約束條件描述

(7)

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

xmin≤x≤xmax

(8)

綜合考慮飛行性能、控制系統性能和發動機性能,可建立控制變量約束模型,有

umin≤u≤umax

(9)

此外,對于本文研究的問題而言,端點約束可分為初始約束與終端約束,即有

(10)

1.3 優化問題描述

為挖掘飛行器運載潛力,將優化目標選取為上升段末端機械能最大,優化目標可以表示為

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

(11)

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

為了提升優化求解的收斂性,往往需要對模型進行無量綱化處理,限于篇幅,此處不再贅述。

2 基于凸優化的求解框架及策略設計

典型凸優化問題即為尋找最優控制量,使得

(12)

2.1 時間自由問題轉化處理

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

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

(13)

同時,將原自變量時間作為新的狀態變量,則有

(14)

(15)

(16)

2.2 凸化處理

考慮到前述上升段軌跡優化模型是非凸的,將其進行凸化處理是應用凸優化技術數值求解的關鍵。本節以姿態控制理想情況下上升段軌跡優化問題為例采用逐次線性化等[11]技術進行模型凸化。

1) 動力學方程凸化

(17)

取其作為初始參數序列,將動力學方程(16)在該參考狀態序列下進行一階泰勒展開,得到線性化的動態方程約束,即

(18)

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

(19)

(20)

其動態方程約束可類比得到,此處不做贅述。

2) 過程約束處理

(21)

對其進行線性化可得

(22)

3) 指標函數凸化處理

記機械能公式

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

(23)

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

(24)

此外,注意到狀態量和控制量表達(8)～(10)式均為凸形式,故無需做凸化處理。至此,該上升段軌跡優化即轉化為一典型的凸問題。

2.3 離散化處理

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

(25)

上標k表示第k次迭代,m表示離散點編號,m=0,1,…,M;Δτ=(τf-τ0)/M。

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

(26)

控制變量約束、狀態變量約束的離散形式為

(27)

狀態變量的端點條件可以進一步表述為

(28)

綜上所述,原軌跡優化問題轉化為參數化凸問題,可表示為如下形式

(29)

2.4 凸優化求解策略改進及流程

1) 終端約束處理

考慮終端高度等端點等式約束,若在優化求解過程中直接限制為固定值,初期迭代過程中很難滿足該約束,甚至會出現不可行解進而導致優化失敗。因此,可將這類約束轉化為性能指標中的懲罰項[11]

(30)

式中:常系數ci>0,i為終端等式約束序號;I為約束總數量?？紤]該懲罰項為非凸形式,故通過引入松弛矢量R,將其轉化為(31)式的形式引入性能指標中

p(R)=c·R

(31)

對應的松弛約束為

(32)

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

于是,上述參數化凸問題可松弛為

(33)

2) 優化求解流程

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

為保證線性化的有效性和解的收斂性及精度,迭代過程以當前優化軌跡不斷更新參考軌跡,并以相鄰2次迭代的凸優化解對應的狀態量最大偏差作為收斂準則,即

(34)

式中,ε為收斂誤差限。

圖1 基于序列凸優化的上升段軌跡優化求解流程圖

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

3 算例仿真及分析

3.1 仿真條件

綜合考慮計算效率與解的精度,本文仿真過程中共取61個離散點,序列凸優化算法的初始信賴域和收斂誤差限設置為

此外,回溯搜索算法中的2個參數均取0.8[15]。

所有仿真均在搭載Intel Core i7-8700 3.20 GHz Intel處理器的臺式機完成,仿真環境為MATLAB 2016b平臺,基于CVX工具包進行軌跡優化算法開發,并調用SDPT3求解器求解凸優化子問題。

3.2 攻角控制系統理想情況下的上升段軌跡優化

將控制系統理想情況下軌跡優化模型簡記為模型一,針對該模型進行終端機械能最大的上升段軌跡優化。凸優化算法經過12輪迭代后收斂,每輪迭代約耗時1.8 s,優化過程共耗時19.88 s,優化結果如圖2～3所示。

圖2 模型一對應的部分狀態量優化結果

圖3 模型一對應的控制量優化結果

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

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

此外,為探究離散點數選取對凸優化算法精度及計算效率的影響,取31和91個離散點分別進行軌跡優化,并利用所得優化控制量進行積分,所得優化及積分結果對比如表1所示。

表1 不同離散點數量優化及積分結果對比

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

3.3 考慮攻角控制系統二階滯后上升段軌跡優化

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

圖4 考慮控制滯后的狀態量優化結果

由結果可知,上升段共飛行136.20 s,剩余質量為800 kg,終端高度為62.18 km、飛行馬赫數為15.59,滿足終端約束,動壓、過載、駐點熱流密度均滿足給定的過程約束限制。

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

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

圖5 考慮控制滯后的控制量優化結果

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

表2 不同ωn取值優化結果

圖6 不同ωn取值優化結果攻角對比

4 結 論

針對RBCC高超聲速上升段軌跡快速優化問題,本文設計了一種基于序列凸優化的軌跡優化算法和策略,相關結論如下:

1) 所設計優化算法可以快速、有效地求解復雜工作模態下RBCC高超聲速上升這一復雜、非線性、強耦合的非凸最優控制問題;

2) RBCC發動機火箭模態峰-谷-峰的工作過程能夠最大限度地結合并發揮火箭發動機工作范圍廣和沖壓發動機經濟性高的優勢,從而使飛行器在燃料更省的情況下獲得更大的終端高度、馬赫數。

3) 相較于理想情況,考慮控制系統二階滯后特性軌跡優化結果的最優性有所降低,且降低程度與自然頻率ωn有關,設計姿態控制系統時應依據實際需要進行控制參數選擇。