基于粒子群-偽譜凸優化的RBCC中段組合軌跡優化方法

摘 要：""""" 針對RBCC中段組合軌跡優化設計問題， 提出了一種基于粒子群-偽譜凸優化的嵌套優化方法。 首先， 根據飛行任務需求給出了中段飛行方案， 并對組合軌跡優化問題進行了描述； 其次， 通過分析組合軌跡各段耦合機理， 將組合軌跡優化問題轉化為段間銜接靜態參數尋優與子段軌跡優化問題， 并設計基于粒子群-偽譜凸優化的雙層嵌套優化策略對該問題進行求解； 上層通過粒子群算法確定靜態參數， 在此基礎上， 下層采用偽譜凸優化方法分段進行軌跡優化設計， 通過偽譜離散和凸化技術的有機結合， 將非凸、 非線性優化問題轉化為離散凸優化問題， 并設計了基于信賴域收縮的序列凸優化求解策略， 在保證各段軌跡最優性的同時， 實現了中段組合軌跡優化問題的快速求解； 最終， 以某RBCC動力概念飛行器為例， 完成了中段組合軌跡優化設計仿真， 驗證了所提方法的可行性和有效性。

關鍵詞："""" RBCC； 中段組合軌跡； 雙層嵌套優化； 粒子群； 偽譜凸優化

中圖分類號：""""" TJ760； V412.4

文獻標識碼：""" A

文章編號："""" 1673-5048（2025）01-0115-11

DOI： 10.12132/ISSN.1673-5048.2024.0067

0 引" 言

火箭基組合循環（Rocket Based Combined Cycle， RBCC）動力系統^［1］將火箭發動機大推力和沖壓發動機高比沖的性能優勢相結合， 具有巡航效率高、 可重復開關機和經濟性好等特點， 可用于執行大空域、 寬速域以及遠程精確投送的飛行任務， 尤其適合作為中段寬空域機動飛行的動力裝置^［2］。 中段飛行是整個飛行過程中飛行距離最長、 軌跡形式（涉及滑翔、 巡航以及縱向機動等）最為靈活多變的飛行階段， 可綜合體現RBCC動力飛行器高效性、 機動性、 靈活性及復雜任務適應性等特點^［3］。 因此， 中段組合軌跡的優化設計是發揮其組合動力綜合性能優勢， 提升飛行效能的重要手段。 然而， RBCC動力飛行器的組合軌跡優化面臨著多重挑戰。 首先， 復雜多樣的飛行任務和臨近空間飛行環境使得中段飛行剖面存在多個典型階段， 各階段面臨多種復雜的過程、 終端及發動機工作條件約束限制； 其次， RBCC動力系統工作模態多、 各模態性能差異大， 切換時機未知， 且沖壓模態與飛行狀態耦合性強， 工作條件極為苛刻^［4］。 此外， 組合軌跡優化涉及多個子系統和模塊， 需進行全面的系統整合和驗證， 以確保各部分協同工作正常。 上述因素使得RBCC動力飛行器中段組合軌跡優化的可行域極小， 設計難度極大， 有必要針對其多階段復雜約束和動力系統特點開展相應的優化求解算法及策略研究。

中段組合軌跡優化設計的本質是求解多階段復雜約束下的最優控制問題。 目前， 間接法和直接法是最常用的兩類方法^［5］。 間接法由于推導復雜、 對初值猜測敏感， 難以應用于這類復雜的多階段、 多約束軌跡優化問題。 而直接法中的偽譜法、 啟發式算法等方法依托其較高的收斂性和可靠性等優勢， 在軌跡優化設計中得到了廣泛應用。 易文雙^［6］以RBCC動力飛行器為研究對象， 利用對初值敏感較低的Gauss偽譜法進行了軌跡設計以及分段和全軌跡射程優化研究。 鄭雄等^［7］研究了RBCC動力飛行器爬升-巡航組合軌跡優化問題，" 采用粒子群-偽譜法的嵌套優化策略得到燃料最省的飛行軌跡。 Chai等^［8］為提升吸氣式飛行器航程以及突防能力， 提出了一種助推跳躍式軌跡， 并采用hp高斯自適應偽譜法進行了優化求解。 近年來， 凸優化方法因其具有多項式復雜度、 理論全局最優性等優勢^［9］， 在火箭回收^［10］、 再入制導^［11］、 交會對接^［12］以及火星著陸^［13］等復雜約束下的軌跡優化問題中取得了廣泛應用。 為了進一步提高凸優化方法的求解精度和效率， 王勁博^［14］針對火箭子級精確軟著陸軌跡優化問題， 通過凸化技術與偽譜離散方法的有機結合， 將非凸、 非線性優化問題轉化為凸優化問題， 實現了考慮阻力的兩階段軌跡優化。 Sagliano^［15］將偽譜法和凸優化相結合形成偽譜凸優化方法， 并將其應用于火星著陸場景， 提高了問題的求解精度。 然而， 針對RBCC動力飛行器組合軌跡優化問題， 目前鮮有基于序列凸優化方法的研究成果。

此外， 目前RBCC動力飛行器軌跡優化的研究一般針對特殊任務進行， 通用性較差， 相關研究主要集中在爬升段和巡航段等場景相對單一的軌跡優化問題中。 值得注意的是， 對于組合軌跡優化問題的求解， 一般將中段組合軌跡優化問題分段進行優化設計， 并采用上述數值優化算法求解， 雖然經典優化方法局部最優性和收斂性較好， 但由于分段進行軌跡優化求解， 且算法多依賴初始猜測， 難以保證軌跡全局最優性。 啟發式算法雖然在全局搜索方面具有優勢， 但求解效率低下， 且由于高速動力學非線性強， 各段之間往往是強耦合的， 計算復雜度高， 難以實際應用。 與經典優化方法對軌跡分段求解后再拼接的全段軌跡設計方法不同， 本文根據實際工程需要， 對軌跡分段進行約束設計， 利用啟發式算法搜索最優段間銜接參數， 在此基礎上， 將原中段組合軌跡優化問題轉化為求解各子段最優軌跡的優化問題， 從而得到全飛行過程最優軌跡， 可解決經典優化方法中分段最優但全過程非最優的問題， 從而提升軌跡的全局性能。

基于上述分析， 本文針對RBCC動力飛行器中段組合軌跡優化問題， 提出了一種基于粒子群-偽譜凸優化的嵌套優化方法。 該方法上層通過粒子群算法確定各段間最優銜接參數， 下層采用偽譜凸優化方法對各段分別進行軌跡優化設計， 將偽譜離散和凸化技術有機結合， 完成非線性動力學和約束的近似凸化和離散處理， 并設計了基于信賴域收縮的序列凸優化求解策略， 實現了多階段、 多控制以及多約束的中段組合軌跡優化問題的快速求解。 最終， 以某RBCC動力概念飛行器中段滑翔-巡航拉起軌跡為例進行仿真， 驗證了所提方法的有效性和計算效率優勢。

1 問題描述

1.1 飛行方案

在軍事偵察、 情報獲取以及遠程投送等飛行任務中， 飛行器的燃料經濟性是首要考慮的因素。 因此， 本文將無動力滑翔與巡航的優勢相結合， 設計了滑翔-巡航拉起的中段組合軌跡方案。 如圖1所示， 飛行器首先由軌跡頂點無動力滑翔至巡航拉起段交班點， 飛行過程中規避相應的禁飛區約束； 當高度、 速度以及彈道傾角滿足交班點條件時， RBCC發動機點火工作使飛行器保持一定高度巡航飛行； 最終， 在目標點附近一定范圍內， 飛行器進行縱向拉起機動， 為俯沖段提供更好的初始條件， 降低被攔截的概率。

采用這種飛行方案的優勢在于： （1）可充分利用無動力滑翔增程， 節省燃料； （2）利用RBCC發動機進行高速巡航， 增加續航和機動能力； （3）巡航末段進行跳躍拉起， 增大軌跡預測難度， 進一步提高末段突防能力； （4）飛行軌跡具有靈活性和對復雜任務的適應性， 從而能夠更好地滿足多樣化的任務需求。

1.2 軌跡優化數學模型

1.2.1 質心運動模型

根據飛行剖面特性， 本文將RBCC動力飛行器中段軌跡分為無動力滑翔段和巡航拉起段。 各階段描述飛行器質心運動的微分方程組均一致， 但推力和燃料秒流量等參數有所不同。 本文假設地球為無旋轉的均質圓球， 飛行器為可控質點， 其姿態控制系統理想工作。 基于上述假設， RBCC動力飛行器質心運動的微分方程組為

r·=vsinθ

λ·=vcosθsinψrcos

·=vcosθcosψr

v·=-D+Pcosαm-gsinθ

θ·=（L+Psinα）cosυmv-gvcosθ+vrcosθ

ψ·=（L+Psinα）sinυmvcosθ+vrcosθsinψtan

m·=-ms （1）

式中： r為地心距， r=h+Re， h為飛行高度， Re為地球半徑； λ為經度； 為緯度； v為速度； θ為彈道傾角； ψ為彈道偏角； m為飛行器質量； α為攻角； υ為傾側角； ms為發動機燃料秒流量； P為發動機推力； g為重力加速度， g=R2eg0/r2， g0為海平面處重力加速度； D和L分別為氣動阻力和升力， 表達式為

D=qSrefCD

L=qSrefCL （2）

式中： q=0.5ρv2為飛行動壓， ρ為大氣密度， ρ=ρ0e^-h/hs， ρ0為海平面處的大氣密度， hs為標準大氣模型的參考高度； Sref為飛行器參考面積； 升、 阻力系數CL、 CD均為攻角和馬赫數的函數， 可由氣動數據庫插值得到。

RBCC動力系統由火箭發動機和沖壓發動機組合而成， 總推力和秒流量由兩部分構成， 具體可表示為

P=PR+PA

ms=msr+msa （3）

式中： PR和PA分別為火箭發動機與沖壓發動機推力； msr和msa分別為火箭發動機與沖壓發動機燃料秒流量。

火箭發動機推力PR計算如下：

PR=msrIspr（4）

式中： Ispr為火箭發動機比沖， 可表征為關于高度的函數fIspr（h）。

沖壓發動機推力PA計算為

PA=msaIspamsa=ρvSaΦEr （5）

式中： Sa， Φ， Er分別為沖壓發動機進氣道橫截面積、 流量系數和燃油當量比； 且Φ可表征為馬赫數和攻角的非線性函數fΦ（Ma， α）； 比沖Ispa可表征為馬赫數的非線性函數fIspa（Ma）。 可見， 沖壓推力與高度、 馬赫數、 攻角及發動機參數Sa、 Er均相關， 即飛行器動力性能與飛行狀態存在強非線性和強耦合關系。

1.2.2 控制變量選擇

以往軌跡優化設計研究中， 均假設RBCC飛行器能瞬時響應指令， 導致設計的軌跡過于理想， 與實際應用存在較大差異， 為了提升控制量變化的連續性和可行性， 本文將火箭發動機燃料秒耗量msr、 攻角變化率uα·、 傾側角變化率uυ·和燃油當量比變化率uEr作為控制量， 攻角α、 傾側角υ和燃料當量比Er為新增狀態量， 并引入補充微分方程， 具體如下：

α·=uα·

υ·=uυ·

E·r=uEr（6）

1.2.3 約束條件

（1） 端點約束

RBCC動力飛行器中段組合軌跡初始和終端約束可表示為

x（t0）=x0

x（tf）=xf（7）

式中： t0， tf分別為中段飛行初始和終端時刻； x0， xf為對應時刻的狀態量約束。

此外， 考慮到軌跡終端為巡航拉起頂點， 為保證飛行器在該點的氣動控制能力， 需滿足最小動壓約束， 即

q（tf）=0.5ρv2（tf）≥qmin（8）

式中： q為動壓， qmin為最小動壓約束值， 在縱向拉起頂點高度一定的情況下， 上述約束可轉化為頂點最小速度約束進行處理。

（2） 段間銜接約束

在軌跡的各段之間， 狀態變量需保持連續， 故引入段間銜接約束， 具體如下：

x（t2， 0）=x（t1， f）（9）

式中： t1， f為無動力滑翔段終端時刻； t2， 0為巡航拉起段初始時刻。

（3） 控制約束

綜合考慮飛行性能、 控制系統性能和發動機性能， 可建立控制變量約束模型， 定義u為控制變量， 則相應控制約束為

umin≤u≤umax（10）

式中： umin和umax分別為控制量u的最小值和最大值。

（4） 狀態約束

選取高度、 速度和彈道傾角等描述飛行器運動的物理量以及攻角和傾側角等描述控制作用的物理量作為狀態量， 建立相應的狀態約束模型， 若用x表示狀態變量， 則相應狀態約束為

xmin≤x≤xmax（11）

式中： xmin和xmax分別為狀態量x的最小值和最大值。

（5） 沖壓發動機工作走廊約束

作為RBCC動力系統重要的組成部分， 沖壓發動機須在一定的高度、 馬赫數約束條件下進行工作， 具體約束如下：

Hlow（Ma）≤h≤Hup（Ma）（12）

式中： Hup（Ma）和Hlow（Ma）分別為沖壓發動機工作走廊的上下邊界。

（6） 過程約束

為保證飛行安全， 在飛行過程中需滿足動壓q、 法向過載n和駐點熱流Q·等約束條件， 具體約束模型如下：

q=0.5ρv2≤qmax

n=Lcosα+Dsinαmg0≤nmax

Q·=KQρv^3.15≤Q·max（13）

式中： qmax， nmax， Q·max為對應的動壓、 法向過載和駐點熱流約束允許最大值； KQ為常值， 其取值與RBCC動力飛行器外形相關。

（7） 禁飛區約束

在飛行過程中， 對于敵對勢力反導區域、 天基武器殺傷區域、 雷達識別區域以及由于地緣政治原因無法通行區域組成的禁飛區， 飛行器需進行規避機動飛行。 本文研究假設禁飛區為無限高“圓柱”區域， 則禁飛區約束可以表示為

fN（λ， ）=R2N-［（λ-λN）2+（-N）2］≤0（14）

式中： RN為禁飛區半徑； λN， N為禁飛區中心的經度和緯度。

圖2給出了RBCC動力飛行器中段組合軌跡優化設計的約束邊界示意圖， 結合上述約束建模和分析可知， 本文所研究的組合軌跡優化問題面臨狀態、 過程、 端點和沖壓工作走廊等復雜約束限制， 這使得軌跡優化的可行域極小， 求解難度極大， 有必要針對該復雜問題開展軌跡快速優化設計方法研究。

1.2.4 優化問題描述

本文考慮飛行器燃料的經濟性， 將中段滑翔-巡航拉起組合軌跡性能指標選取為

J=m（tf）→max（15）

綜上所述， 本文研究的RBCC動力飛行器中段軌跡優化問題P0描述為： 尋找最優的狀態量x和控制量u， 使得在滿足運動方程約束式（1）和（6）、 端點約束式（7）、 段間銜接條件式（9）、 控制約束式（10）、 狀態約束式（11）、 沖壓發動機工作條件約束式（12）、 過程約束式（13）以及禁飛區約束式（14）的條件下， 性能指標式（15）達到最大值。 為了提升優化求解的收斂性， 往往需要對模型進行無量綱化處理， 限于篇幅， 此處不再贅述。

2 算法設計

2.1 優化策略

對于中段整體飛行過程， 考慮到攜帶燃料有限以及存在最佳燃料利用率等情況， 將中段軌跡劃分為無動力滑翔段和巡航拉起段， 分段進行軌跡優化設計。 其中， 無動力滑翔段需為RBCC發動機開機提供合適的初始高度、 速度彈道傾角及航向角條件， 而巡航拉起段則需快速抵近指定目標點， 并在巡航末段實現縱向拉起機動。 基于上述分析， 引入不同子段終端約束如下：

無動力滑翔段終端約束為

h（t1， f）=h

v（t1， f）=v

θ（t1， f）=θ1f

ψ（t1， f）=ψLOS （16）

式中： t1， f為滑翔段的終端時刻； h， v分別為滑翔段的期望終端高度、 速度； θ1f為期望的彈道傾角； ψLOS為期望的目標方位角。

巡航拉起段終端約束為

h（t2， f）=hpeak

v（t2， f）≥vpeak

θ（t2， f）=θpeak

（λ（t2， f）， （t2， f））=（λT， T） （17）

式中： t2， f為巡航拉起段終端時刻； hpeak為巡航拉起頂點高度約束； vpeak為頂點速度約束； θpeak為頂點彈道傾角約束； λT， T分別為目標點經度和緯度。

為使中段整體軌跡最優， 即飛行器終端剩余質量最大， 需合理選取各段的性能指標。 為降低巡航拉起段燃料消耗， 滑翔段的性能指標設置為滑翔終端與目標點距離最短， 具體為

J1=λ（t1， f）-λT+（t1， f）-T→min（18）

巡航拉起段性能指標同式（15）一致， 此處不再贅述。

本文無動力滑翔段的狀態量和控制量為

x=［r， λ， ， v， θ， ψ， m， α， υ］T

u=［uα·，" uυ·］T（19）

巡航拉起段的狀態量和控制量分別為

x=［r， λ， ， v， θ， ψ， m， α， υ， Er］T

u=［msr， uα·， uυ·， uEr］T（20）

根據飛行方案可知， 無動力滑翔段終端需為巡航拉起段提供合適的交班條件（包括終端高度h和速度v等）， 其余狀態量作為下一階段的初始條件輸入。 在巡航拉起段， 在目標點指定且飛行器燃料質量一定的情況下， 存在合適的巡航高度和速度使得燃料利用率最優。 可見， 段間銜接高度h和速度v的取值對RBCC動力飛行器中段軌跡最優性影響較大。 因此， 選擇二者作為優化變量， 其所受約束如下：

hmin≤h≤hmax

vmin≤v≤vmax（21）

式中： hmin和hmax分別為銜接點高度的下限和上限； vmin和vmax分別為銜接點速度的下限和上限， 可根據沖壓工作條件予以確定。

綜上， 考慮到粒子群算法具有概念簡單、 智能背景深刻以及全局收斂性好的特點， 本文采用粒子群算法確定段間銜接高度、 速度等靜態參數。 在確定段間銜接靜態參數的基礎上， 采用偽譜凸優化方法分段求解中段滑翔-巡航拉起多階段、 多約束的軌跡優化問題。

基于上述分析， 本文針對中段組合軌跡優化設計中的參數與最優控制混合優化問題， 設計了基于粒子群靜態參數尋優和偽譜凸優化動態軌跡優化嵌套的參數/軌跡一體化優化求解策略。 其中， 上層通過粒子群算法確定靜態參數， 下層在參數優化結果的基礎上采用偽譜凸優化方法進行軌跡優化設計， 進而通過嵌套框架的內外環迭代， 實現參數/軌跡一體化優化設計。 基于本文提出的“粒子群-偽譜凸優化”嵌套優化求解策略， 可同時優化段間銜接高度、 速度以及相應的最優軌跡控制指令， 確保中段組合軌跡的最優性。

2.2 粒子群算法

設粒子群算法中有K個粒子， 其中， 每個粒子k（k=1， 2， …， K）表示優化問題解空間中的一個備選解。 所有粒子都有自己的位置x（k）=［x1（k）， …， xn（k）］T和速度v（k）=［v1（k）， …， vn（k）］T以及由目標函數J=f（x（k））計算得到的適應度， 其中n為待優化變量的個數， 通過比較各粒子的適應度來判斷粒子的優劣。 根據每個優化變量的取值范圍xi∈［ai， bi］， i=1， 2， …， n， 對粒子位置和速度的限制如下：

ai≤xi（k）≤bi

|vi（k）|≤|ai-bi| （22）

式中： i=1， 2， …， n； k=1， 2， …， K。

基本粒子群算法的粒子和位置更新公式參考文獻［16］， 在此不再贅述。 此外， 本文粒子群算法的適應度同巡航拉起段的性能指標式（15）一致， 粒子群算法結束條件則根據具體問題選擇達到最大迭代次數或搜索到的最優位置滿足預期精度。

2.3 基于偽譜凸優化的組合軌跡優化算法

基于前文的分析， 通過2.2節粒子群算法確定段間銜接高度h和速度v等靜態參數后， 中段組合軌跡優化問題P0可采用分段優化策略進行軌跡優化設計。 然而， P0是一個高度非線性的優化問題， 需要對該問題進行凸化、 離散化處理， 進而采用序列凸優化算法求解。 為方便后文描述， 將動力學方程式（1）、 （6）改寫為：

x·=fp（x）+Bpu（23）

式中： p=1， 2； 函數fp（x）為各段的微分方程右端函數； 矩陣Bp為相應的常值矩陣。

因此， 本節首先通過偽譜離散方法對問題P0中的微分方程式（23）、 過程約束式（12）～（13）、 禁飛區約束式（14）以及目標函數式（5）、 （18）進行離散處理， 再通過逐次線性化將非凸約束轉化為凸約束， 進而得到離散凸優化問題， 最終， 通過設計基于信賴域收縮策略的序列凸優化方法實現對該問題的分段求解。

2.3.1 偽譜離散化處理

目前， 偽譜離散方法是相對成熟的技術， 本文研究并非圍繞偽譜法本身展開， 而是著眼于偽譜法與凸優化方法有機結合。 因此， 有關偽譜法的技術細節在此不再討論， 與本文研究相關的Radau偽譜法， 詳細內容見文獻［17］。

首先， 對于微分動力學約束， 直接給出其Radau偽譜離散化形式如下：

∑Ni=0Djix（τi）=tf-t02f［x（τj）， u（τj）， τ;t0， tf］（24）

式中： j=1， 2， …， N； Dji為偽譜微分矩陣D-∈R^N×（N+1）中的元素； τi（0≤i≤N）是用來離散狀態變量的節點， 其中包括N個配點和τ0=-1。

為了將問題P0轉化為終端時間自由的軌跡優化問題， 將各段的終端飛行時間tf作為擴展控制變量， 結合式（24）在此直接給出針對某一飛行階段的系統動力學方程離散形式如下：

2∑Ni=0Djixi+faug， p（yi）=0（25）

式中： j=1， 2， …， N； p=1， 2； faug， p（yi）=（t0-tf）fp（x（τi）， u（τi））， t0和tf分別為當前離散階段內的初始和終端時刻； yi={x（τi）， u（τi）， tf}， x（τi）和u（τi）分別為當前離散階段內的離散狀態和控制。

對于原始軌跡優化問題中的性能指標以及其他過程、 禁飛區或端點約束， 其離散化過程較為簡單， 因此， 有關偽譜法的技術細節在此不再討論， 與

本文研究相關的Radau偽譜法可參考文獻［17］。

2.3.2 問題的凸化

（1） 動力學方程凸化

可以看出， faug， p中存在非線性項， 使用一階線性化對式（25）進行凸化處理， 可得

2∑Ni=0Djixi+faug， p（xk， uk， tkf）+Ap（xk， tkf）·（x-xk）+

Bp（tkf）·（u-uk）+Tp（xk， uk）·（tf-tkf）=0（26）

化簡得到

2∑Ni=0Djixi+Ap（xk， tkf）·x+Bp（tkf）·u+

"Tp（xk， uk）·tf+Wp（xk， uk， tkf）=0

Wp（xk， uk， tkf）=faug， p（xk， uk， tkf）-Ap（xk， tkf）·

xk-Bp（tkf）·uk-Tp（xk， uk）·tkf（27）

式中： xk， uk和tkf為第k次迭代時線性化展開的參考點， faug， p（xk， uk， tkf）為相應的參考軌跡； 矩陣Ap=faug， p/x， Bp=faug， p/u和Tp=faug， p/tf。 以巡航拉起段為例， 具體為

A2=000a14a1500000

a210a23a24a25a260000

a3100a34a35a360000

a4100a44a450a47a4800

a5100a54a550a57a58a590

a610a63a64a65a66a67a68a690

a7100a74000a780a710

0000000000

00000000000000000000（28）

B2=0000

0000

0000

b41000

b51000

b61000

-（t0-tf）000

0（t0-tf）00

00（t0-tf）0

000（t0-tf）（29）

T2=t11t21t31t41t51t61t71t81t91t101T（30）

上述矩陣Ap、 Bp、 Tp（p=1， 2）中的各元素為faug， p（x， u， tf）對狀態量、 控制量以及終端時間tf求一階偏導數的結果， 對于難以直接求取解析偏導數的項， 可在參考軌跡（xk， uk， tkf）附近利用有限差分法近似求取偏導數， 限于篇幅， 各系數具體形式在此不再贅述。

此外， 為保證泰勒展開一階線性化的準確性， 引入信賴域約束， 具體如下：

x-xk≤δx

u-uk≤δu

tf-tkf≤δt（31）

式中： δx， δu， δt分別為針對狀態量以及控制變量設計的信賴域。

（2）" 過程約束凸化

原問題模型P0中的過程約束均為非線性約束， 同樣利用上述方法對離散的過程約束進行線性化， 將非凸約束近似凸化。 考慮到過程約束均為地心距r、 速度v和攻角α的函數， 故在基準狀態（rk， vk， αk）附近進行線性化展開。

駐點熱流約束、 動壓約束與地心距r、 速度v相關， 則有

fi（r， v）≈fi（rk， vk）+f′i（rk， vk）·［r-rk;v-vk］（32）

過載約束與地心距r、 速度v、 攻角α相關， 則有

fn（r， v， α）≈fn（rk， vk， αk）+f′n（rk， vk， αk）·［r-rk; v-vk; α-αk］（33）

綜上可得

fQ·（rk， vk）+fQ·r（r-rk）+fQ·v（v-vk）≤Q·max

fq（rk， vk）+fqr（r-rk）+fqv（v-vk）≤qmax

fn（rk， vk， αk）+fnr（r-rk）+fnv（v-vk）+fnα（α-αk）≤nmax（34）

式中： fQ·（r， v）， fq（r， v）， fn（r， v， α）分別代表駐點熱流、 動壓與過載約束。

（3） 禁飛區約束凸化

假設禁飛區為無限高“圓柱”區域， 其表達式如式（14）所示。 將禁飛區約束在基準狀態（λk， k）附近進行線性化， 則有

fN（λk， k）+fNλ（λ-λk）+fN（-k）≤0（35）

式中： 偏導數fNλ=2（λN-λk）， fN=2（N-k）。

（4） 目標函數凸化

對于式（18）所示的目標函數， 引入兩個松弛變量η1和η2， 將目標函數等價轉化為

J1=η1+η2（36）

松弛變量η1、 η2服從于附加約束：

|λ（t1， f）-λT|≤η1

|（t1， f）-T|≤η2（37）

目標函數式（36）關于松弛變量η1和η2是線性的， 且服從于線性約束式（37）。

至此， 已完成問題的離散化、 凸化處理。

綜上所述， 以巡航拉起段為例， 可得偽譜凸優化問題P1， 即

minx， u J2=-m（t2， f）

s.t.

2∑Ni=0Djixi+A2（xk， tkf）·x+B2（tkf）·u+T2（xk， uk）·

tf+W2（xk， uk， tkf）=0

x（t0）=x0， x（tf）=xf

x∈［xmin， xmax］

α∈［αmin， αmax］， |uα·|≤uα·max

υ∈［υmin， υmax］， |uυ·|≤uυ·max

msr∈［msrmin， msrmax］， |uEr|≤uErmax

fQ·（rk， vk）+fQ·r（r-rk）+fQ·v（v-vk）≤Q·max

fq（rk， vk）+fqr（r-rk）+fqv（v-vk）≤qmax

fn（rk， vk， αk）+fnr（r-rk）+fnv（v-vk）+

fnα（α-αk）≤nmax

fN（λk， k）+fNλ（λN-λk）+fN（N-k）≤0

|x-xk|≤δx， |u-uk|≤δu， |tf-tkf|≤δt

（38）

2.3.3 偽譜序列凸優化求解策略

凸優化問題P1與原問題P0存在線性化誤差， 需通過設計序列凸優化算法對P1進行迭代求解， 以逼近原始問題P0的解。

本文算法迭代收斂條件為

max|x^k+1-xk|≤ε（39）

式中： ε為常值收斂誤差向量。

（1） 信賴域收縮策略

在確定初始猜測（x0， u0， t0f）的基礎上， 為了確保本文凸優化迭代求解的收斂性， 對信賴域約束進行調整， 提出一種信賴域收縮策略^［18］： 記初始迭代時選取的信賴域半徑為δ0， 在后續迭代求解的過程中， 對比相鄰兩代求解所得的性能指標， 分別記為J^k-1和Jk， 則

δk=βδ^k-1Jk≤J^k-1

δ^k-1Jkgt;J^k-1 （40）

式中： β為信賴域收縮系數， 其范圍為0lt;βlt;1。

由于初始信賴域半徑δ0是在保持線性化精度的前提下選取的， 因此上述信賴域收縮策略能在保證凸化可靠性的前提下改善問題的迭代收斂性。

（2） 信賴域收縮的偽譜序列凸優化求解流程

綜上所述， 本文所提方法求解流程如圖3所示， 具體步驟為：

Step1： k=0， 設置算法參數δ0、 β以及ε， 由線性初值給定參考軌跡（x0， u0， t0f）。

Step2： k=k+1， 利用參考軌跡（xk， uk， tkf）建立凸優化問題P1， 并利用凸優化算法進行求解， 獲得當前最優解（x^k+1， u^k+1， t^k+1f）。

Step3： 判斷最優解是否滿足式（39）。 若滿足， 轉入Step5， 否則轉入Step4。

Step4： 根據式判斷是否更新信賴域約束， 同時， 令k=k+1， 生成下一代的基準軌跡（xk， uk， tkf）， 返回Step2。

Step5： 迭代終止， 輸出軌跡優化的收斂解（xk， uk， tkf）。

2.4 基于粒子群-偽譜凸優化的嵌套優化求解流程

綜上， 本文通過上層粒子群算法確定靜態參數， 優化變量為段間銜接高度和速度， 下層在參數優化結果的基礎上采用偽譜凸優化方法分段進行軌跡優化設計， 進而通過嵌套框架進行內外環迭代， 當迭代步數達到最大值， 或者相鄰兩次種群的最優適應度之差小于某一小量時， 迭代停止， 輸出最優參數以及最優軌跡， 否則繼續迭代尋優， 從而實現參數/軌跡一體化優化設計。 該求解流程如圖4所示。

嵌套優化策略的具體步驟如下：

Step1： 設定粒子群算法、 偽譜凸優化方法相應參數。

Step2： 確定段間銜接高度和速度的上下限， 隨機初始化粒子的位置和速度。 每個粒子的位置包含一組{h， v}， 以此作為滑翔段的終端高度、 速度約束以及巡航拉起段的初始條件。

Step3： 在每一組{h， v}下， 采用偽譜凸優化方法對滑翔-巡航拉起軌跡進行分段優化設計， 得到相應的最優性能指標， 并將其作為粒子的適應度。

Step4： 評估每個粒子的適應值， 比較粒子當前適應值與個體最優值的大小， 若當前適應函數值較小， 則替換為個體最優值。

Step5： 評估粒子當前適應值與群體最優值的大小關系， 若當前適應值較小， 則替換為群體最優值。

Step6： 判斷是否滿足結束條件， 若否， 更新種群中所有粒子的速度與位置， 并進入下一輪尋優， 返回Step3； 反之， 則結束， 輸出最優設計參數及最優軌跡。

3 算例仿真及分析

3.1 仿真對象

本文以某RBCC動力概念飛行器為研究對象開展仿真分析， 該飛行器采用面對稱氣動布局。 圖5給出了該飛行器在不同攻角和馬赫數條件下的升阻比曲線， 由圖可知， 升阻比與馬赫數、 攻角的關系并不單調， 在Malt;3的情況下， 升阻比隨攻角的增大而增大， 最大升阻比攻角為10°； 在Ma≥3的情況下， 升阻比隨攻角的增大先增大后減小， 最大升阻比攻角為6°； 全局最大升阻比出現在馬赫數為4， 攻角為6°的情況下， 此時飛行器的氣動增程能力達到最大。 圖6則給出了不同馬赫數下沖壓模態的比沖曲線， 由圖可知， 隨馬赫數的增大， 沖壓模態比沖先增大后減小， 并在Ma=4.5附近取得最大值， 此時飛行器的動力續航能力達到最大。

基于上述氣動、 動力性能參數影響因素分析， 綜合考慮飛行器氣動增程和動力續航能力， 該飛行器的最佳巡航馬赫數應處于4和4.5之間， 最佳滑翔和巡航攻角應保持在6°附近。 此外， 圖7則進一步給出了沖壓發動機工作走廊邊界曲線， 由圖可知， 隨著馬赫數的增大， 沖壓發動機工作高度邊界整體增大， 可行高度范圍小幅縮小。

3.2 仿真條件

RBCC動力飛行器軌跡中段各段軌跡的初始條件和終端條件如表1所示， 飛行過程中各段軌跡的主要狀態約束如表2所示。 飛行器的初始攻角α0=6.5°， 初始傾側角υ0=0°； 控制變量邊界約束為： α∈［-6°， 10°］， υ∈［-60°， 60°］， |uα·|≤0.15（°）/s， |uυ·|≤2（°）/s， msr∈［0， 40］ kg/s， 0≤Er≤1以及|uEr|≤0.1； 過程約束中： 最大熱流密度Q·max=1 000 kW/m2， 最大動壓qmax=200

kPa、" 最大過載nmax=5g；" 禁飛區約束中心位置分別為（20°， 5°）， （35°， 10°）， 半徑分別為400 km， 200 km。

基本粒子群算法參數設定如下： 種群數K=15， 最大迭代次數iter=30， 種群維度D=2， 慣性權重w0=0.729， 個體學習因子c1=2， 群體學習因子c2=2； 優化變量范圍： 段間銜接高度h∈［30， 40］ km， 速度v∈［1 400， 1 500］ m/s。

考慮到組合軌跡的飛行時間較長， 本文將偽譜序列凸優化算法中離散節點數設為101， 初始猜測為線性初值， 信賴域收縮系數β=0.7， 初始信賴域約束和收斂條件設置為

δ=［10 km， 30°， 30°， 1 000 m/s， 30°， 30°， 30°， 30°， 200 s， 0.05］

ε=［0.2 km， 0.05°， 0.05°， 10 m/s， 0.05°， 0.05°， 0.05°， 0.05°， 1 s， 0.05］

3.3 結果分析

基于3.2節仿真條件開展RBCC動力飛行器中段組合軌跡優化設計仿真， 所得仿真結果如圖8和表3所示。 優化所得銜接點高度h為35.2 km， 速度v為1 432.5 m/s， 對應馬赫數為4.21， 終端剩余質量為1 124.34 kg。 優化所得巡航馬赫數在4～4.5范圍內， 進一步驗證了3.1節氣動、 動力性能參數分析的合理性。 由表3可知， 飛行器拉起頂點終端高度、 速度及位置偏差均較小， 滿足各項終端約束要求。

由圖8（a）～（c）可知， 無動力滑翔、 巡航拉起段均滿足各項狀態、 端點約束要求， 完成了禁飛區的有效規避， 實現了預設的飛行方案。 其中， 巡航拉起段的飛行高度、 馬赫數長期穩定在35 km、 4.2附近， 以實現高效巡航， 節省燃料消耗， 同時滿足沖壓發動機工作條件約束。 由圖8（d）可知， 飛行過程中， 動壓、 過載以及熱流密度均滿足給定的過程約束限制， 符合設計要求。

由圖8（e）～（f）可知， 飛行過程中攻角和傾側角滿足相應控制約束要求， 在滑翔段， 攻角大部分時間保持在6°最大升阻比攻角附近， 以實現遠距離滑翔， 再次驗證了3.1節氣動參數分析的合理性； 而在巡航拉起段， 隨著飛行器燃料消耗， 飛行器質量減小， 維持巡航所需的升力也相應減小， 因此攻角逐漸減?。?最后， 為實現跳躍拉起機動， 飛行器將攻角提升至最大攻角10°， 以充分發揮大攻角進氣和增升作用， 為拉起機動提供充足的法向過載。 由圖8（f）可知， 飛行器在滑翔段通過2次傾側反轉以實現禁飛區的有效規避， 在巡航段則通過3次傾側反轉使飛行軌跡始終處于射面附近， 盡可能降低燃料消耗， 提升最優性。

由圖8（g）～（i）可知， 火箭和沖壓發動機秒流量均滿足相應約束要求。 在巡航拉起段起點， 火箭發動機短暫開機提供推力， 快速調整飛行狀態以滿足高效巡航要求； 隨后， 火箭發動機關機， 由沖壓發動機提供推力， 使高度、 馬赫數分別維持在35 km、 4.2附近， 充分發揮沖壓發動機比沖高， 巡航效率高的優勢， 進一步印證了3.1節的動力性能分析結論； 最后， 為實現拉起跳躍機動， 火箭發動機再次開機， 并配合大攻角沖壓進氣， 產生組合大推力， 迅速提高飛行器高度和速度， 完成末段縱向拉起機動。

以上結果表明， 所提組合軌跡優化算法可根據不同飛行階段特點， 充分發揮RBCC動力綜合性能和氣動性能優勢， 具有較好的可行性和最優性。

為進一步體現本文方法在計算效率方面的優勢， 將本文方法與基于hp偽譜法^［8］的組合軌跡優化方法進行對比。 2種方法的仿真條件與3.2節一致， hp偽譜法采用GPOPS優化工具包實現， 設置求解器為snopt， 網格收斂精度為1×10^-5。 所得仿真結果如圖9和表4所示。

由圖9可知， 兩種算法的優化結果曲線基本一致； 由表4可知， 兩種方法剩余質量基本相當，" 而本文方法的計算耗時僅為hp偽譜法的4.72%， 進一步說明了本文方法在保證可行性和最優性的同時， 具有更高的計算效率優勢。

4 結" 論

針對RBCC動力飛行器中段組合軌跡優化設計問

題， 本文研究了一種基于粒子群-偽譜凸優化嵌套優化的參數/軌跡一體化優化求解方法， 仿真結果表明：

（1） 所提雙層嵌套優化策略可用于處理含靜態參數、 動態軌跡優化的中段組合軌跡優化設計問題， 具有較好的可行性和最優性。

（2） 所提下層偽譜序列凸優化方法可用于快速求解多約束、 多階段、 強非線性的RBCC動力飛行器軌跡優化問題， 對復雜問題適應性較好。

（3） 與現有方法相比， 所提方法在優化設計結果與優化目標相當的情況下， 優化耗時降低了95.28%， 進一步體現了所提算法在計算效率方面的優勢。

參考文獻：

［1］ 王亞軍，" 何國強，" 秦飛，" 等. 火箭沖壓組合動力研究進展［J］. 宇航學報，" 2019，" 40（10）： 1125-1133.

Wang Yajun，" He Guoqiang，" Qin Fei，" et al. Research Progress of Rocket Based Combined Cycle Engines［J］. Journal of Astronautics，" 2019，" 40（10）： 1125-1133.（in Chinese）

［2］ 曾家，" 黃輝，" 朱平平，" 等. 火箭基組合動力研究進展與關鍵技術［J］. 宇航總體技術，" 2022，" 6（3）： 49-57.

Zeng Jia，" Huang Hui，" Zhu Pingping，" et al. Research Progress and Key Technology Analysis of Rocket Based Combined Cycle Engines［J］. Astronautical Systems Engineering Technology，" 2022，" 6（3）： 49-57.（in Chinese）

［3］ 張遠龍，" 謝愈. 滑翔飛行器彈道規劃與制導方法綜述［J］. 航空學報，" 2020，" 41（1）： 023377.

Zhang Yuanlong，" Xie Yu. Review of Trajectory Planning and Gui-dance Methods for Gliding Vehicles［J］. Acta Aeronautica et Astronautica Sinica，" 2020，" 41（1）： 023377.（in Chinese）

［4］ 閆循良， "王培臣，" 王舒眉，" 等. 基于混沌多項式的RBCC飛行器上升段魯棒軌跡快速優化［J］. 航空學報，" 2023，" 44（21）： 306-323.

Yan Xunliang，" Wang Peichen，" Wang Shumei，" et al. Rapid Robust Trajectory Optimization for RBCC Vehicle Ascent Based on Polynomial Chaos［J］. Acta Aeronautica et Astronautica Sinica，" 2023，" 44（21）： 306-323.（in Chinese）

［5］ 崔乃剛，" 郭冬子，" 李坤原，" 等. 飛行器軌跡優化數值解法綜述［J］. 戰術導彈技術，" 2020（5）： 37-51.

Cui Naigang，" Guo Dongzi，" Li Kunyuan，" et al. A Survey of Numerical Methods for Aircraft Trajectory Optimization［J］. Tactical Missile Technology，" 2020（5）： 37-51.（in Chinese）

［6］ 易文雙. 動力滑翔高超飛行器軌跡設計與射程優化［D］. 哈爾濱： 哈爾濱工業大學，" 2014.

Yi Wenshuang. Trajectory Design and Range Optimization of Gliding Hypersonic Vehicle with Powered［D］.Harbin： Harbin Institute of Technology，" 2014. （in Chinese）

［7］ 鄭雄，" 劉竹生，" 楊勇，" 等. 火箭基組合循環高超聲速飛行器爬升-巡航全局軌跡優化研究［J］. 導彈與航天運載技術，" 2018（2）： 1-8.

Zheng Xiong，" Liu Zhusheng，" Yang Yong，" et al. Research on Climb-Cruise Global Trajectory Optimization for RBCC Hypersonic Vehicle［J］. Missiles and Space Vehicles，" 2018（2）： 1-8.（in Chinese）

［8］ Chai D，" Fang Y W，" Wu Y L，" et al. Boost-Skipping Trajectory Optimization for Air-Breathing Hypersonic Missile［J］. Aerospace Science and Technology，" 2015，" 46： 506-513.

［9］ Boyd S，" Vandenberghe L. Convex Optimization［M］.Beijing： World Publishing Corporation，" 2004.

［10］ 邵楠，" 閆曉東. 火箭垂直回收多階段最優軌跡規劃方法［J］. 宇航學報，" 2019，" 40（10）： 1187-1196.

Shao Nan，" Yan Xiaodong. Multi-Stage Trajectory Optimization for Vertical Pin-Point Landing of a Reusable Launch Vehicle［J］. Journal of Astronautics，" 2019，" 40（10）： 1187-1196.（in Chinese）

［11］ Liu X F，" Shen Z J，" Lu P. Entry Trajectory Optimization by Se-cond-Order Cone Programming［J］. Journal of Guidance， Control， and Dynamics，" 2016，" 39（2）： 227-241.

［12］ Lu P，" Liu X F. Autonomous Trajectory Planning for Rendezvous and Proximity Operations by Conic Optimization［J］. Journal of Guidance， Control， and Dynamics，" 2013，" 36（2）： 375-389.

［13］ Acikmese B，" Ploen S R. Convex Programming Approach to Po-wered Descent Guidance for Mars Landing［J］. Journal of Guidance， Control， and Dynamics，" 2007，" 30（5）： 1353-1366.

［14］ 王勁博. 可重復使用運載火箭在線軌跡優化與制導方法研究［D］. 哈爾濱： 哈爾濱工業大學，" 2019.

Wang Jinbo. Research on Online Trajectory Optimization and Guidance Methods for Reusable Rocket［D］.Harbin： Harbin Institute of Technology，" 2019. （in Chinese）

［15］ Sagliano M. Pseudospectral Convex Optimization for Powered Descent and Landing［J］. Journal of Guidance， Control， and Dynamics，" 2018，" 41（2）： 320-334.

［16］ Eberhart R，" Kennedy J. A New Optimizer Using Particle Swarm Theory［C］∥ the Sixth International Symposium on Micro Machine and Human Science， IEEE，" 1995： 39-43.

［17］ Kameswaran S，" Biegler L T. Convergence Rates for Direct Transcription of Optimal Control Problems Using Collocation at Radau Points［J］. Computational Optimization and Applications，" 2008，" 41（1）： 81-126.

［18］ 宋瑞，" 朱勇，" 徐廣通，" 等. 基于序列凸優化的高超聲速飛行器協同再入軌跡規劃［J］. 戰術導彈技術，" 2020（6）： 7-16.

Song Rui，" Zhu Yong，" Xu Guangtong，" et al. Cooperative Reentry Trajectory Planning of Hypersonic Vehicle Based on Sequential Convex Programming［J］. Tactical Missile Technology，" 2020（6）： 7-16.（in Chinese）

RBCC Mid-section Combined Trajectory Optimization Method

Based on

Particle Swarm-Pseudospectral Convex Optimization

Yang Yuxuan1， Fei Wanghua2， Liu Haili3， Wang Peichen1， Yan Xunliang^1*

（1. School of Astronautics， Northwestern Polytechnical University，" Xi’an 710072，" China;

2. Research Development Center of China Academy of Launch Vehicle Technology，" Beijing 100076，" China;

3. Jiangnan Electromechanical Design Institute，" Guiyang 550009，" China）

Abstract： In order to solve the problem of combined trajectory optimization of RBCC mid-section，" a nested optimization method based on particle swarm-pseudospectral" convex optimization is proposed. Firstly，" according to the requirements of the mission，" the mid-course flight scheme is given，" and the combined trajectory optimization problem is described. Secondly，" by analyzing the coupling mechanism of each segment of the combined trajectory，" the combined trajectory optimization problem is transformed into the inter-segment convergence static parameter optimization and sub-segment trajectory optimization problem，" and a double-layer nested optimization strategy based on particle swarm-pseudospectral convex optimization is designed to solve the problem. The static parameters are determined by the particle swarm optimization algorithm in the upper layer，" and on this basis，" the pseudospectral convex optimization method is used to optimize the trajectory design in segments，" and the non-convex and nonlinear optimization problems are transformed into convex optimization problems by the organic combination of pseudo-spectral discrete method and convex technology，" and the sequential convex optimization solution strategy based on trust domain shrinkage is designed，" which not only ensures the optimality of each segment of trajectories，" but also realizes the rapid solution of the middle segment combined trajectory optimization problem. Finally，" the simulation of mid-section combined trajectory optimisation design is completed with an RBCC-powered concept vehicle as an example，" which verifies the feasibility and effectiveness of the proposed method.

Key words： RBCC; mid-section combined trajectory; double-layer nesting optimization; particle swarm; pseudospectral convex optimization