基于hp自適應偽譜法的組合動運載器發(fā)射窗口拓展能力分析

周宏宇，王小剛，崔乃剛，許河川

（1. 哈爾濱工業(yè)大學 航天學院，哈爾濱 150001；2. 哈爾濱建成機械廠，哈爾濱 150030）

發(fā)射窗口（Launch Window, LW）是指完成航天任務所允許的發(fā)射時間間隔[1]。傳統(tǒng)運載火箭在發(fā)射后基本只能在射面內(nèi)飛行，因而有效載荷在主動段結束后所進入的軌道平面主要取決于發(fā)射位置、發(fā)射時刻和發(fā)射方向，其中，發(fā)射位置和發(fā)射方向決定了軌道傾角，發(fā)射時刻決定了軌道升交點赤經(jīng)。通常發(fā)射位置和發(fā)射方向很容易保證，但發(fā)射時刻會由于各種問題而造成延遲。一旦錯過了預定發(fā)射時刻，載荷進入的實際軌道和預定軌道間將產(chǎn)生一定的夾角，即軌道偏差，此時若要消除入軌偏差，就需要消耗上面級載荷的燃料進行變軌，導致載荷運行壽命減低，或是重新制定發(fā)射計劃，等待下一個發(fā)射窗口的到來。因此，為順利完成任務，采用運載火箭的發(fā)射方式必須在規(guī)定的較小的時間區(qū)間內(nèi)完成點火發(fā)射，因而發(fā)射窗口很窄，在靈活性上具有很大局限性。

基于組合動力的水平起飛飛行器（Combined Cycle Launch Vehicle, CCLV）具有巡航飛行的能力，其高機動性使其能夠在巡航段進行橫向機動[2]，改變飛行方向，進而改變?nèi)胲壍能壍绤?shù)，即使錯過了發(fā)射時刻也能通過橫向機動進入預定軌道面，由此實現(xiàn)了對發(fā)射窗口的拓寬。

為分析CCLV對發(fā)射窗口的拓寬能力，首先要通過軌跡優(yōu)化獲知其最大橫向機動能力。由于飛行過程中氣動環(huán)境復雜、動力性能受飛行環(huán)境影響較大、運動耦合性強、約束條件嚴格、控制量可能存在尖點[3]，傳統(tǒng)的打靶法、Gauss偽譜法和間接法等難以解決這類優(yōu)化問題[4]。文獻[5]通過簡化運動模型，利用間接法實現(xiàn)了高超聲速飛行器軌跡快速優(yōu)化；文獻[6]針對含有復雜約束條件的非線性最優(yōu)控制問題，提出了改進的高斯偽譜法，但為提高精度必須增加節(jié)點，由此增大了計算量，甚至導致問題病態(tài)；文獻[7]結合 hp自適應偽譜法和模式搜索法，解決了水平起飛組合動力飛行器的上升段軌跡優(yōu)化問題。

針對傳統(tǒng)發(fā)射方式發(fā)射窗口窄、發(fā)射時刻不靈活的問題，提出了采用組合動力飛行器作為運載器的新思路。建立了組合動力運載器的數(shù)學模型，包括動力學模型和組合動力發(fā)動機模型；給出了發(fā)射窗口的計算方法，推導了橫向機動范圍與發(fā)射窗口的關系；給出了hp自適應偽譜法的算法流程，并針對爬升段和巡航段分別選取控制參數(shù)和性能指標，建立了優(yōu)化模型。最后通過仿真驗證了優(yōu)化算法的有效性，并根據(jù)優(yōu)化結果分析了發(fā)射窗口能力。

1 組合動力運載器數(shù)學模型

組合動力運載器的動力系統(tǒng)由多種動力模式組成：首先采用渦輪基組合動力發(fā)動機（Turbo based Combined Cycle, TBCC）水平起飛、爬升和加速，當高度和速度達到一定值時，切換到火箭基組合動力發(fā)動機（Rocket Based Combined Cycle, RBCC），進一步加速和爬升。由于TBCC和RBCC發(fā)動機有著各自的優(yōu)缺點和適用環(huán)境，組合動力可以將二者結合使用，從而實現(xiàn)單級入軌和重復使用。

1.1 組合動力運載器動力學模型

忽略地球扁率與旋轉，飛行器質(zhì)心運動方程如下：

式中：V為飛行速度大小，γ為飛行路徑角，P為推力大小，α為攻角，g為重力加速度大小，σ為傾側角，ψ為飛行航向角，r為飛行器到地心的距離，m為飛行器質(zhì)量，m˙為發(fā)動機秒耗量，X和Y分別為氣動阻力和升力。

式中：q為飛行動壓，S為特征面積；xc和yc分別為阻力系數(shù)和升力系數(shù)，是攻角和馬赫數(shù)的函數(shù)，可通過多項式擬合獲得。

1.2 發(fā)動機模型

對于TBCC，推力和比沖是飛行高度和飛行馬赫數(shù)的函數(shù)，為便于分析，通過多項式計算推力和比沖：

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

對于RBCC，火箭發(fā)動機秒耗量越大，則沖壓發(fā)動機所占比例越小，對應推力越大、比沖越小[8]。為便于分析，將推力和比沖視作火箭發(fā)動機所占比重（記為K）的函數(shù)：

式中： Tmax為最大推力， Imax為最大比沖， K ∈[0,1]稱為相對流量。本文取 Tmax=360 kN， Imax=3000 s。

2 發(fā)射窗口建模

2.1 發(fā)射窗口計算方法

在此定義飛行器在主動段結束后所進入的軌道為初始軌道，期望進入的軌道為目標軌道。下面推導無橫向機動情況下飛行器的發(fā)射窗口寬度計算方法及橫向機動拓寬發(fā)射窗口的原理。

若飛行器始終在初始射面內(nèi)飛行，則當發(fā)射方位角和發(fā)射點緯度一定時，可根據(jù)式(5)計算初始軌道的軌道傾角：

式中：i為軌道傾角，A為發(fā)射方位角即發(fā)射方向，φ為發(fā)射點緯度。

如果同時確定了發(fā)射時刻，則由式(6)可求得初始軌道的升交點赤經(jīng)：

式中：sΩ為初始軌道升交點赤經(jīng)，Lθ為發(fā)射時刻對應的恒星時。

當初始軌道和目標軌道的軌道傾角不同時，兩軌道間的夾角為

式中： i1為初始軌道的傾角， i2為目標軌道的傾角，ΔΩs為初始軌道和目標軌道間的升交點赤經(jīng)差。

因此，若要使最終的軌道面偏差不超過Δα，則發(fā)射窗口為

式中：eω地球自轉角速度。

圖1 機動拓寬發(fā)射窗口示意圖Fig.1 Sketch for extending the launch window

假設預定的發(fā)射時刻為 Tc，則在無橫向機動的情況下，飛行器必須在時間區(qū)間[Tc-Δ t / 2 ,Tc+Δt /2]內(nèi)點火發(fā)射才能保證入軌精度。

由上文敘述可知，發(fā)射時刻上的偏差最終將導致升交點赤經(jīng)上的偏差，即軌道面偏差，而通過橫向機動可以調(diào)整軌道面，進而彌補發(fā)射時刻上偏差，實現(xiàn)拓寬發(fā)射窗口。如圖1所示，初始軌道和目標軌道的軌道傾角相同，但升交點赤經(jīng)由于發(fā)射時刻不準確而有所差異，兩個軌道面間存在夾角ΔΩs= Ωs1- Ωs2；通過橫向機動，飛行器從初始軌道轉移到目標軌道，在不改變軌道傾角的情況下改變了升交點赤經(jīng)。因此，由發(fā)射時刻上的提前或推遲所造成的軌道面偏差可以通過橫向機動修正回來，進而獲得更大的發(fā)射窗口和發(fā)射靈活性。

2.2 發(fā)射窗口與橫向機動的關系建立

圖2給出了飛行器橫向機動的俯視圖，其中，ab為初始發(fā)射縱平面，飛行器經(jīng)橫向機動后到達c點。近似認為飛行器的轉彎軌跡為一個圓弧，則圓弧上 a點和c點處的切線相交于d點，且有ad = c d。β為機動后的軌道面和初始軌道面間的夾角，且有

式中：cL和dL分別為機動過程中飛行器飛過的橫程和縱程，cL即為橫向機動距離。

圖2 β計算方法示意圖Fig.2 The calculation method of β

若 Lc是在考慮實際運動狀態(tài)和約束條件下所獲得的最大橫程，則β便是對初始軌道面的最大調(diào)整角度。因此，當給定所允許的最大軌道面偏差Δα時，具有橫向機動能力的飛行器可允許的初始軌道面偏差為 Δα ′ = Δα + β 。若認為實際軌道和目標軌道間不存在軌道傾角偏差，則由式(8)可得對應的發(fā)射窗口：

3 軌跡優(yōu)化方法

3.1 hp自適應偽譜法

偽譜法在配點上利用P階多項式近似狀態(tài)量，為了提高近似精度，常采用增加節(jié)點（h方法[9]）或增加多項式階次（p方法[10]）的方式。h方法通過細化網(wǎng)格降低差值多項式的階次，但網(wǎng)格過密會增大非線性規(guī)劃（NLP）問題的規(guī)模，影響求解速度；p方法即增加多項式階次，具有設計變量少、初值易選取等優(yōu)點，但同樣會影響求解速度。

hp自適應偽譜法通過估計各時間區(qū)間上的狀態(tài)量曲率來判斷是采用h方法還是p方法來提高求解精度，其中采用h法還能夠自行確定新增節(jié)點數(shù)量及位置。hp自適應偽譜法結合了h方法和p方法的優(yōu)點，避免了不必要的網(wǎng)格細化因而減小了NLP的規(guī)模，利于提高求解速度[11]，同時具備求解精度高、收斂快、初值易選擇、魯棒性強等特點。

3.2 算法流程

1） 時間區(qū)間轉換

將時間區(qū)間 t ∈[t0,tf]分為K個區(qū)間，并將每個區(qū)間轉換至 τ ∈[-1 ,1]：

式中： t ∈ [ tk-1,tk]，k = 1 ,2,… ,K 。

2）計算離散誤差和約束誤差

將最優(yōu)控制問題離散化，計算各區(qū)間內(nèi)各狀態(tài)量的微分方程離散誤差 e(k)和過程約束誤差 b(k)。

3）精度提升方法選取

如果第k個區(qū)間內(nèi)的 e(k)和b(k)都小于設定閾值εd，則認為該區(qū)間能夠精確近似原最優(yōu)控制問題，優(yōu)化結束；否則進一步優(yōu)化網(wǎng)格。

如果 rk＜rkmax（rkmax為設定的正常數(shù)），則認為該區(qū)間內(nèi)各狀態(tài)量的曲率差別較大，一致性不足，需要增加網(wǎng)格內(nèi)插值多項式的階次來提高精度；否則增加網(wǎng)格內(nèi)的配點數(shù)目。

4）多項式階次的確定

設 D1、D0分別為更新前后插值多項式階次，為 ε(k)=[e(k), b(k)]中的最大值，ceil()為正向取整函數(shù)，則：

式中：A為任意常整數(shù)。

5）新增節(jié)點數(shù)與節(jié)點位置的確定

更新后配點數(shù)的計算方法如下：

式中：Y為任意正常數(shù)。

新節(jié)點位置通過曲率密度函數(shù)來決定。曲率密度函數(shù)定義為

式中：c為設定常數(shù)，滿足 F （ 1） = 1 。

3.3 優(yōu)化模型建立

飛行器先后分別在 TBCC和 RBCC的作用下飛行，故將飛行分為TBCC段和RBCC段?？紤]到不同飛行階段的動力模式特點、任務需求和飛行環(huán)境，需針對各飛行階段選取不同的優(yōu)化變量、約束條件和性能指標。

TBCC段的主要任務是爬升和加速，在該階段飛行器不進行橫向機動。該階段優(yōu)化變量為攻角，性能指標為燃料最省，約束條件為α＜30°，終端約束為動力切換時的高度、速度和飛行路徑角。

RBCC段的任務是繼續(xù)加速爬升，并通過橫向機動進行軌道調(diào)整。該階段的優(yōu)化變量為攻角、傾側角和相對流量，性能指標為橫程最大，約束條件為α＜8°、σ＜ 40°、0 ≤ K ≤ 1。該階段對攻角的約束更加嚴格，這是由于RBCC發(fā)動機必須保證足夠的進氣量，過大的攻角會導致發(fā)動機熄火。終端約束為主動段結束時的飛行路徑角、高度和飛行器質(zhì)量，即在給定燃料質(zhì)量下完成飛行。

因此，共有4個優(yōu)化變量和2個性能指標，且各階段的約束條件也有所不同。此外，在 TBCC段和RBCC段，還需對法向過載和動壓進行約束。本文統(tǒng)一規(guī)定動壓不超過90 kPa，法向過載不超過5.0。

4 仿真分析

4.1 仿真初始條件

首先采用優(yōu)化方法進行軌跡優(yōu)化。由式(3)可知，當高度大于25 km時，TBCC推力趨于平穩(wěn)，推力值較小，且TBCC的工作馬赫數(shù)一般不超過 5.0。故本文將動力切換時機設定為高度20 km，速度1100 m/s。飛行器起飛質(zhì)量90 t，燃料質(zhì)量80 t，入軌點高度為110 km。其它仿真條件見表1。

表1 各階段初始及終端運動狀態(tài)Tab.1 Initial and final states in different phases

此外，由式(10)可知，發(fā)射窗口與目標軌道傾角有關，因此需要分析同樣機動能力在不同目標軌道傾角下對發(fā)射窗口的拓展能力。本文所選擇的目標軌道傾角范圍為 20°~80°，優(yōu)化出最大橫程后，利用 2.2節(jié)中的方法求出發(fā)射窗口。αΔ取1°，即允許1°的軌道面偏差。

4.2 仿真結果

TBCC段燃料消耗5 506.5 kg，飛行73.4 s；RBCC段最大橫程達到441.1 km，對應縱程652.4 km，飛行200.1 s。仿真結果表明，本文所采用的優(yōu)化算法能夠滿足終端約束及過程約束，并且僅在需要的地方增加了節(jié)點，由此減小了NLP問題的規(guī)模。

圖3中攻角在動力切換時發(fā)生了突變，這是由于攻角作為控制量并未對其連續(xù)性進行約束，而法向過載是與攻角有關的量，故圖7中法向過載也在同一時刻產(chǎn)生了突變。實際應用中，可以在TBCC段結束時先通過姿控將攻角調(diào)整至RBCC段的起始攻角，然后開啟RBCC。

對于橫程441.1 km、縱程652.4 km，各目標軌道傾角對應的發(fā)射窗口如表2所示?？梢钥闯觯繕塑壍纼A角越低，同樣的橫向機動距離獲得的發(fā)射窗口更寬，這一點也可由式(10)推導得出；同時，組合動力飛行器的橫向機動能力大大的拓寬了發(fā)射窗口。

此外，由式(5)可知，傳統(tǒng)發(fā)射方式的初始軌道傾角不能小于發(fā)射點緯度，要進入低傾角軌道必須消耗上面的燃料級進行變軌。如果采用本文的方案，則能夠通過橫向機動在上升段進入低傾角軌道，避免消耗上面級燃料。

圖3 攻角變化情況Fig.3 Time history of attack angle

圖5 RBCC段傾側角變化情況Fig.5 Time history of leaning angle in RBCC

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

圖6 RBCC段相對流量變化情況Fig.6 Time history of relative flow in RBCC

圖7 動壓和法向過載變化情況Fig.7 Time history of dynamic pressure and load

圖8 飛行剖面Fig.8 The flight profile

表2 發(fā)射窗口計算結果Tab.2 The result of launch window

5 結 論

本文針對載人航天任務未來的發(fā)展趨勢，設計了一種組合動力模式，并采用優(yōu)化方法進行了軌跡優(yōu)化設計，分析了橫向機動能力和發(fā)射窗口拓展能力間的關系，具體如下：

1）采用基于組合動力的水平起飛運載器，可以克服傳統(tǒng)發(fā)射方式窗口窄的問題，發(fā)射方式和時機更加靈活，還能在上升段進入軌道傾角小于發(fā)射點緯度的軌道。

2）hp自適應偽譜法能根據(jù)優(yōu)化進程自適應調(diào)整離散網(wǎng)格，自行選擇采用h方法或是p方法，在提高精度的同時降低了不必要的網(wǎng)格加密，以免NLP問題規(guī)模過大，因此提高了收斂速度，能夠有效求解軌跡優(yōu)化問題。

3）目標軌道傾角越低，發(fā)射窗口越易拓寬。

（

）：

[1] Deerwester J M, Mclaughlin J F, Wolfe J F. Earthdeparture plane change and launch window considerations for interplanetary missions[J]. Journal of Spacecraft &Rockets, 2015, 3(2): 169-174.

[2] Pu Z Q, Yuan R Y, Tan X M, et al. Active robust control of uncertainty and flexibility suppression for air-breathing hypersonic vehicles[J]. Aerospace Science & Technology,2015, 42: 429-441.

[3] Bowcutt K G.. Multidisciplinary optimization of air-breathing hypersonic vehicles[J]. Journal of Propulsion &Power, 2015, 17(6): 1184-1190.

[4] Takashima N, Lewis M J. Optimization of waverider based hypersonic cruise vehicles with off-design considerations[J]. Journal of Aircraft, 2015, 36(1): 235-245.

[5] Murillo O, Lu P. Fast ascent trajectory optimization for hypersonic air-breathing vehicles[C]//AIAA Guidance,Navigation, and Control Conference. 2013: 1-24.

[6] 孫勇. 基于改進 Gauss偽譜法的高超聲速飛行器軌跡優(yōu)化與制導[D]. 哈爾濱工業(yè)大學, 2012: 41-59.Sun Y. Trajectory optimization and guidance of hypersonic vehicle based on improved gauss pseudospectral method[M]. Harbin Institute of Technology, 2012: 41-59 .

[7] 周宏宇, 王小剛, 崔乃剛, 等. 基于hp自適應偽譜法的組合動力可重復使用運載器軌跡優(yōu)化[J]. 中國慣性技術學報, 2016, 24(6): 832-837.Zhou H Y, Wang X G, Cui N G, et al. Trajectory optimazation of reusable vehicle with combined power based on hp adaptive pseudospectral algorithm[J]. Journal of Chinese Inertial Technology, 2016, 24(6): 832-837.

[8] Wang Y J, Jiang L I, Qin F, et al. Research on thermal adjustment method of RBCC in ramjet mode[J]. Journal of Solid Rocket Technology, 2016, 39(2): 151-159.

[9] Gandhi M. Trajectory optimization algorithm studies[J].Mathematics, 2015, 15(2): 267-281.

[10] Ross I M, Fahroo F. Pseudospectral knotting methods for solving non-smooth optimal control problems[J]. Journal of Guidance Control & Dynamics, 2004, 27(3): 397-405.

[11] 夏紅偉, 李秋實, 李莉, 等. 基于hp自適應偽譜法的飛行器再入軌跡優(yōu)化與制導[J]. 中國慣性技術學報, 2015,23(6): 818-823.Xia H W, Li Q S, Li L, et al. Trajectory optimization and guidance for reentry craft based on hp-adaptive pseudospectral method[J]. Journal of Chinese Inertial Technology, 2015, 23(6): 818-823.