王開強 張柏楠
(中國空間技術(shù)研究院載人航天總體部,北京 100094)
隨著超燃沖壓發(fā)動機的發(fā)展,高超聲速巡航飛行器逐漸進入科研人員的研究視線。其推進效率相比于傳統(tǒng)火箭得到大幅提高,可在20km 以上的高度進行馬赫數(shù)大于5 的高超聲速巡航飛行,隨后返回地面的預(yù)定機場水平著陸,在未來的空基發(fā)射入軌領(lǐng)域具有很好的發(fā)展前景。目前,僅有美國X-43 高超聲速試驗飛行器實現(xiàn)了10s時間量級、馬赫數(shù)大于7 的穩(wěn)態(tài)巡航飛行[1]。對于可返回高超聲速巡航飛行器,其飛行階段較多,且各階段均包含了眾多不同類型的約束條件。因此其飛行軌跡設(shè)計較為復(fù)雜,需要借助優(yōu)化的方法手段。
目前,國內(nèi)外已開展了一些關(guān)于巡航軌跡優(yōu)化的研究。文獻[2-3]采用差分進化算法和MATLAB 優(yōu)化工具箱等,對高超聲速乘波體飛行器在周期性巡航和穩(wěn)態(tài)巡航2 種模式下的軌跡優(yōu)化問題進行了研究。文獻[4-5]基于航空交通管制規(guī)則,采用分支定界法對民航飛機的巡航程序和軌跡進行了優(yōu)化。文獻[6]基于EUROCONTROL’S 軟件的模型和優(yōu)化算法,考慮風(fēng)的影響,得到了不同巡航高度下飛行時間和耗油量最優(yōu)的軌跡;文獻[7]采用模式搜索的確定性算法對巡航軌跡進行了優(yōu)化,以改善民航飛行管理;文獻[8]同樣應(yīng)用確定性優(yōu)化算法,對定高巡航飛行軌跡優(yōu)化算法進行了研究;文獻[9-13]僅針對高超聲速周期性巡航軌跡進行了優(yōu)化研究;文獻[14]采用遺傳算法對平流層飛艇的巡航機動航跡進行了優(yōu)化。目前已有的這些研究大多僅僅針對巡航段,采用不同的直接法、間接法,配合不同的優(yōu)化算法,進行軌跡設(shè)計優(yōu)化。
關(guān)于高超聲速巡航飛行器全軌跡優(yōu)化方面,目前的研究較少。僅有文獻[15]將其飛行軌跡分為助推段、上升段、巡航段、俯沖段4個階段進行優(yōu)化,其中巡航段事先人為確定了一個固定的巡航高度、馬赫數(shù)等,以便于進行其它飛行階段的優(yōu)化。軌跡多目標(biāo)優(yōu)化方面,國內(nèi)已有一些關(guān)于高超聲速再入飛行器的相關(guān)研究[16-18],尚無針對高超聲速巡航飛行器的軌跡多目標(biāo)優(yōu)化研究。實際上,可返回的高超聲速巡航飛行器相比于再入飛行器,其飛行階段劃分更多,各階段可能存在獨特的優(yōu)化目標(biāo),并且各階段之間可能存在一定的先后順序和耦合關(guān)聯(lián),因此有必要對該類飛行器的全飛行軌跡開展分步與多目標(biāo)優(yōu)化研究。
本文面向空基發(fā)射入軌任務(wù),對配置超燃沖壓發(fā)動機的可返回地面預(yù)定機場的高超聲速巡航飛行器全軌跡分步與多目標(biāo)優(yōu)化問題進行研究。首先,對本文使用的飛行動力學(xué)模型、飛行器模型和地球物理模型進行介紹。其次,對整個全軌跡分步多目標(biāo)優(yōu)化問題進行描述,對軌跡分步多目標(biāo)優(yōu)化的策略步驟、優(yōu)化算法、優(yōu)化建模等內(nèi)容進行說明。然后,對剩余段軌跡多目標(biāo)優(yōu)化中采用的第三代直接搜索域(Directed Search Domain-III,DSD-III)的經(jīng)典多目標(biāo)優(yōu)化算法進行簡要描述。最后給出優(yōu)化算例和仿真結(jié)果,并基于得到的多目標(biāo)優(yōu)化解集在設(shè)計空間中的分布特性,進一步分析設(shè)計解選擇決策的建議,為實際工程設(shè)計提供參考。
本文采用經(jīng)度—緯度—高度3個變量描述飛行器質(zhì)心在三維空間中運動的位置信息,采用的動力學(xué)方程如下(其中考慮了地球扁率,以及地球自轉(zhuǎn)引起的牽連加速度和科氏慣性加速度)
式中t為飛行時間;v、γ、Az、m、h、λ、δ分別為飛行速度、航跡傾斜角、速度方位角(也稱航向角)、飛行器質(zhì)量、高度、經(jīng)度和緯度,這7個變量構(gòu)成了本文飛行器飛行軌跡的狀態(tài)變量;F為發(fā)動機推力;m˙為發(fā)動機質(zhì)量流量;α為飛行攻角;g為地球重力加速度;ωe=7.292×10–5rad/s 為地球自轉(zhuǎn)角速度;re為地球半徑;D、L分別為飛行過程中受到的氣動阻力和升力,可由公式計算得到
式中Ar是飛行器的參考面積;CD、CL分別為阻力系數(shù)、升力系數(shù),一般由馬赫數(shù)、飛行攻角計算得到,或在各系數(shù)隨馬赫數(shù)和攻角的二維插值表中進行線性插值得到;q為動壓為大氣密度)。
本文飛行器在助推爬升段和巡航飛行段分別使用固體火箭發(fā)動機、超燃沖壓發(fā)動機進行飛行。對于固體火箭發(fā)動機,為便于分析計算,假設(shè)其推力和比沖均為一恒定值,因此有
式中ve為發(fā)動機比沖,單位為m/s。
對于超燃沖壓發(fā)動機,本文根據(jù)桑格爾的發(fā)動機數(shù)據(jù)表[19],通過線性插值得到推力和質(zhì)量流量。其中,插值變量為馬赫數(shù)Ma、高度h、攻角α。
飛行器包含固體火箭助推器(包含過渡段)和巡航飛行器本體2 部分,整個組合體在發(fā)射時刻的初始質(zhì)量m0如下
式中m1d、m1p分別為助推器干質(zhì)量和推進劑質(zhì)量,m2d、m2p分別為飛行器本體干質(zhì)量和燃料質(zhì)量。
當(dāng)固體助推器結(jié)束工作轉(zhuǎn)入無動力滑行段時,固體助推器被拋離,此時無動力滑行段初始時刻飛行器的質(zhì)量m02為
本文采用J2引力模型描述重力加速度和地球半徑
式中μe=398600km3/s2為地球引力常數(shù);J2=1.08264×10–3為J2引力常數(shù);ra=6378.135km 和rb=6356.912km分別為地球赤道半徑和極半徑。
大氣模型采用美國SA76 模型,其根據(jù)飛行高度計算得到大氣密度ρ和大氣溫度T,此時一定高度上的聲速vs可由公式(7)[15]計算
本文中飛行任務(wù)剖面與美國X43 的飛行過程類似,不同的是本文中的高超聲速巡航飛行器需安全返回至地面預(yù)定的機場著陸。整個飛行剖面主要分為以下5個階段,如圖1 所示。
圖1 飛行剖面Fig.1 Flight profile
1)爬升段:整個固體火箭助推器與飛行器的組合體由載機運送至高空進行空中發(fā)射,然后在固體助推器的作用下進行爬升,直至推進劑燃盡后關(guān)機,并與飛行器本體分離。
2)無動力滑行段:飛行器本體與助推器分離后,進行無動力滑行至規(guī)定的巡航高度和馬赫數(shù)。
3)巡航段:飛行器的超燃沖壓發(fā)動機點火,進行定高定速的穩(wěn)態(tài)巡航飛行,直至燃料燃盡后關(guān)機,巡航段結(jié)束。
4)返回段:飛行器向預(yù)定的著陸機場進行無動力的返回飛行,并達到規(guī)定的進場著陸入口飛行狀態(tài)。
5)進場著陸段:飛行器瞄準著陸機場,通過變換下降飛行的軌跡策略,逐步降低高度和速度,并調(diào)整航跡傾斜角,直至在機場著陸。
本文軌跡優(yōu)化中,將爬升段和無動力滑行段合并為爬升滑行段進行優(yōu)化。此外,進場著陸段的飛行軌跡設(shè)計方法較為特殊,與前面各飛行階段完全不同,不屬于本文研究的范疇,因此本文的軌跡優(yōu)化僅包含爬升滑行段、巡航段和返回段。
本文將整個軌跡優(yōu)化分為第一步巡航段優(yōu)化和第二步剩余段優(yōu)化兩步依次進行,各段有其對應(yīng)的優(yōu)化目標(biāo)。主要有以下兩方面考慮:
1)從飛行任務(wù)層面,空基發(fā)射入軌任務(wù)中,巡航段是最主要的飛行階段,入軌飛行器的空基投放就是在該段完成。因此,巡航軌跡與飛行主任務(wù)直接相關(guān),應(yīng)首先予以優(yōu)化確定。
2)從軌跡優(yōu)化步驟層面,剩余飛行段的軌跡優(yōu)化須基于巡航飛行狀態(tài)進行。其中爬升滑行段需要初始巡航飛行狀態(tài)作為其軌跡優(yōu)化的終端狀態(tài)約束,從而在該段結(jié)束時順利轉(zhuǎn)入高超聲速巡航飛行;而返回段的軌跡優(yōu)化則需要以巡航段結(jié)束時刻的飛行狀態(tài)作為狀態(tài)變量的初始值。因此,巡航飛行狀態(tài)也是其它飛行段軌跡優(yōu)化的前提基礎(chǔ)。
整個軌跡分步多目標(biāo)優(yōu)化的流程如圖2 所示。首先對巡航段飛行軌跡進行單目標(biāo)優(yōu)化,得到最優(yōu)巡航初始飛行狀態(tài)參數(shù)。然后基于此進行剩余段的軌跡多目標(biāo)優(yōu)化,其首先根據(jù)最優(yōu)初始巡航狀態(tài)進行爬升滑行段的軌跡優(yōu)化;然后進行巡航段飛行仿真,得到巡航段結(jié)束時刻的飛行狀態(tài)參數(shù);最后基于此進行返回段的軌跡優(yōu)化,最終得到整個飛行軌跡。
圖2 分步軌跡優(yōu)化流程Fig.2 Procedure of multistep trajectory optimization
優(yōu)化算法方面,第一步巡航段軌跡單目標(biāo)優(yōu)化采用序列二次規(guī)劃算法(Sequential Quadratic Programming,SQP);第二步剩余段則采用經(jīng)典多目標(biāo)優(yōu)化方法中的第三代直接搜索域法DSD-III[20-21]。
(1)巡航段優(yōu)化目標(biāo)
巡航段的優(yōu)化目標(biāo)對應(yīng)于飛行任務(wù)級目標(biāo)。在本文針對的空基發(fā)射入軌任務(wù)中,為入軌飛行器提供更多的初始機械能是其追求的目標(biāo)。
(2)剩余段優(yōu)化目標(biāo)
剩余段軌跡優(yōu)化的目標(biāo)是減小爬升滑行段和返回段的飛行難度,在保證達到規(guī)定巡航飛行狀態(tài)要求的前提下,減小整個飛行任務(wù)執(zhí)行的難度。本文主要從駐點峰值熱流密度和軌跡振蕩程度兩方面評價飛行難度,并作為優(yōu)化目標(biāo)考慮:
1)駐點峰值熱流密度最?。簽榱藴p小高超聲速巡航時的阻力,巡航飛行器的外形非常扁平,因而其前緣駐點處的半徑很小,此時飛行過程中的駐點熱流密度較大。根據(jù)X43 飛行器的相關(guān)數(shù)據(jù)[22],其在30km以馬赫數(shù)7 和2°攻角進行巡航時的駐點熱流約為6.2MW/m2。因此,本文研究的這一類巡航飛行器中,駐點峰值熱流問題較為突出,應(yīng)予以優(yōu)化減小,從而減輕對結(jié)構(gòu)熱防護的設(shè)計壓力。
2)軌跡振蕩最?。号郎卸魏头祷囟蔚娘w行可能出現(xiàn)一定的軌跡振蕩,頻繁、大幅的軌跡振蕩對飛行控制不利。尤其在返回段,軌跡振蕩可能影響轉(zhuǎn)進場著陸段時的飛行狀態(tài),進而對最終的著陸過程產(chǎn)生不利影響。因此,以軌跡振蕩最小為優(yōu)化目標(biāo),有助于降低整個任務(wù)的飛行難度。
考慮到上述兩個優(yōu)化目標(biāo)存在此消彼長的制約關(guān)系,本文采用多目標(biāo)優(yōu)化方法求解剩余段軌跡優(yōu)化問題,以期得到多目標(biāo)優(yōu)化解集(Pareto 非劣解集)。
巡航段優(yōu)化目標(biāo)為
式中vcr和hcr分別為巡航速度和高度。該目標(biāo)可以在空基發(fā)射入軌任務(wù)中為入軌飛行器提供更多的初始機械能。
巡航段設(shè)計變量為:1)巡航初始時刻的高度h0cr;2)巡航初始時刻的馬赫數(shù)Ma0cr;3)巡航期間的飛行攻角αcr。在定高定速的穩(wěn)態(tài)巡航模式中,巡航高度和馬赫數(shù)基本不變。因此由上述3個設(shè)計變量即可確定巡航段的飛行狀態(tài)。三者的取值范圍須處于超燃沖壓發(fā)動機的工作窗口內(nèi)。本文參考“桑格爾”沖壓發(fā)動機[19]的窗口參數(shù)給出取值范圍如表1 所示。
表1 設(shè)計變量取值范圍Tab.1 Ranges of the design variables
巡航段約束條件為:
1)根據(jù)穩(wěn)態(tài)巡航的推阻平衡和升質(zhì)量平衡條件,得巡航期間飛行馬赫數(shù)、航跡傾斜角的理論變化量應(yīng)為零,但在實際飛行中很難保證上述變化量絕對為零,需要給定一個很小的容許變化范圍,本文對最大、最小巡航馬赫數(shù)和以及巡航期間最大航跡傾斜角最小航跡傾斜角作如下限定:
2)動壓約束,實際飛行動壓q應(yīng)小于動壓約束上邊界
3)駐點峰值熱流約束,駐點熱流Qs的最大值即駐點峰值熱流Qsm應(yīng)小于上邊界
剩余段軌跡優(yōu)化為二維多目標(biāo)優(yōu)化
Δγsum為軌跡振蕩總量,用航跡傾斜角的全程累積變化量表示
式中tf為全軌跡的總飛行時間。
剩余段軌跡優(yōu)化的設(shè)計變量為:
1)爬升滑行段飛行攻角αcg,取值范圍為–10°≤αcg≤25°;
2)返回段飛行攻角αrt,取值范圍為0°≤αrt≤5°。
剩余段軌跡優(yōu)化的約束條件為:
1)動壓約束:q≤qub;
2)過載約束:實際過載nload應(yīng)小于許用過載上限
4)爬升滑行段終端飛行狀態(tài)約束。爬升滑行段結(jié)束時,其高度hfcg和馬赫數(shù)Mafcg必須與第一步優(yōu)化得到的巡航段入口高度h0cr和馬赫數(shù)Ma0cr相等,且此時航跡傾斜角γfcg應(yīng)基本為0,這里同樣設(shè)置一個很小的容許范圍0.1°,則有以下約束:
5)進場著陸初始飛行狀態(tài)約束。本文飛行器最終在地面航程已定的指定機場著陸,因此返回段結(jié)束時須滿足一定的關(guān)于高度、速度、航跡傾斜角和航程的進場著陸條件
式中hfr為飛行高度fh的要求值;和分別為飛行速度vf的上、下邊界;和分別為航跡傾斜角γf的上、下邊界;和則分別限定了航程xf的上、下邊界。
一個標(biāo)準的多目標(biāo)優(yōu)化問題可由以下形式表達:
式中x是設(shè)計變量;f(x)為多目標(biāo)函數(shù)向量;fi(i=1, 2, …,n)為第i個目標(biāo)函數(shù);n為目標(biāo)函數(shù)數(shù)量;D*是可行設(shè)計空間。各目標(biāo)函數(shù)構(gòu)成了目標(biāo)空間。滿足約束條件的目標(biāo)空間稱為可行目標(biāo)空間,記為Y*。
多目標(biāo)優(yōu)化目標(biāo)函數(shù)中,任意兩個目標(biāo)之間都是此消彼長的制約關(guān)系,因此其解不唯一。Pareto 優(yōu)化解,也稱多目標(biāo)優(yōu)化非劣解,其對應(yīng)的所有n個目標(biāo)函數(shù)取值不可能同時小于任何其它優(yōu)化解對應(yīng)的相應(yīng)目標(biāo)函數(shù)值。否則,該解為劣解,不屬于Pareto 優(yōu)化解。
本文采用第三代直接搜索域DSD-III 的經(jīng)典多目標(biāo)優(yōu)化方法。其基于錨點進行多目標(biāo)尋優(yōu)。每個錨點對應(yīng)于一個優(yōu)化目標(biāo)的單目標(biāo)最優(yōu)解,因此n個目標(biāo)函數(shù)對應(yīng)n個錨點。例如,在圖3 所示的二維多目標(biāo)優(yōu)化問題中,f1和f2為兩個目標(biāo)函數(shù),因而總共存在兩個錨點。其中,左上方的錨點對應(yīng)于目標(biāo)f1最優(yōu)(最小)時的單目標(biāo)優(yōu)化解,而右下方的錨點則對應(yīng)于目標(biāo)f2最優(yōu)時的單目標(biāo)優(yōu)化解。過所有錨點可構(gòu)建一個目標(biāo)空間,稱為理想超平面。在圖3 中,該超平面實際為經(jīng)過兩個錨點的直線(圖3 中的虛線)。圖中左下方黑色曲線段為Pareto 前緣,其中的點構(gòu)成了Pareto 前緣點集。
圖3 多目標(biāo)優(yōu)化的基本定義Fig.3 Basic definitions of MOO
DSD-III 可以較為高效的求得多目標(biāo)優(yōu)化Pareto 前緣點集,且該點集能較為全面和均勻的覆蓋Pareto前緣。Pareto 點集的分布均勻性,能夠以有限的點表征優(yōu)化目標(biāo)之間的量化制約變化信息,為進一步對優(yōu)化解的選擇和決策提供支持。有關(guān)該方法的詳細描述可參閱原始文獻[21]。
算例中的助推器及巡航飛行器各部分質(zhì)量如表2 所示。巡航飛行器由載機進行空基投放發(fā)射,投放高度為12km,投放初速度為236m/s,初始航跡傾斜角為5°。其中,固體火箭助推器總推力為240kN,比沖為2893m/s。
表2 助推器及巡航飛行器質(zhì)量Tab.2 Mass of booster and the cruise vehicle單位:kg
約束條件方面,動壓約束設(shè)置為不超過qub=130kPa,過載不超過5,駐點峰值熱流小于5MW/m2,進場著陸初始狀態(tài)約束參考航天飛機[23]設(shè)置如下:
本文采用MATLAB R2016a 軟件進行整個優(yōu)化過程的仿真計算。其中采用了fmincon 優(yōu)化器,選用SQP 算法作為底層的單目標(biāo)優(yōu)化算法。
(1)第一步優(yōu)化結(jié)果
第一步巡航段單目標(biāo)最優(yōu)軌跡如表3 所示,飛行器巡航期間高度和馬赫數(shù)基本不變,滿足穩(wěn)態(tài)巡航的狀態(tài)要求。
表3 巡航段優(yōu)化解Tab.3 The cruise optimal solution
(2)第二步優(yōu)化結(jié)果
基于上述巡航飛行狀態(tài),得到的剩余段軌跡二維多目標(biāo)優(yōu)化Pareto 前緣點集如圖4 所示。其中,所有點均為有效點,未得到冗余的優(yōu)化點。同時,該前緣點集較為均勻的勾勒出了該多目標(biāo)優(yōu)化問題的Pareto 前緣,其變化趨勢具有較強的非線性。分析該Pareto 點集的分布可知,在圖中左側(cè)駐點峰值熱流Qsm較小的區(qū)域中,隨著航跡傾斜角全程累積變化量Δγsum的大幅降低,Qsm的增加量很小。而隨著Qsm的增加,航跡傾斜角累積變化量的降低幅度逐漸減小。在圖4 中的右側(cè),當(dāng)Qsm>4.35MW/m2時,降低Δγsum需要付出的熱流代價較高。因此在設(shè)計解的選擇和決策時,建議在Pareto 前緣點集分布圖中相對靠左的區(qū)域(Qsm<4.35MW/m2)選擇設(shè)計解,避免熱流代價很高的設(shè)計解區(qū)域。最終的設(shè)計解,還需進一步結(jié)合防熱技術(shù)發(fā)展現(xiàn)狀和實際工程具體設(shè)計需求進行分析確定。
圖4 剩余段軌跡優(yōu)化Pareto 前緣點集Fig.4 Pareto frontier set of rest trajectory optimization
圖4 中位于最左上方的Pareto 點對應(yīng)于駐點峰值熱流Qsm最小的單目標(biāo)優(yōu)化解,最右下方的Pareto 點則對應(yīng)于軌跡振蕩程度Δγsum最小的單目標(biāo)優(yōu)化解。兩者的對比如圖5 所示,主要軌跡參數(shù)如表4 所示。
圖5 兩個軌跡優(yōu)化解對比Fig.5 Comparison of the two optimal trajectories
表4 軌跡參數(shù)Tab.4 Trajectory parameters
由表4 所示數(shù)據(jù)可得,所有約束條件均滿足,且剩余段軌跡多目標(biāo)優(yōu)化結(jié)果體現(xiàn)了兩個優(yōu)化目標(biāo)的最優(yōu)性以及兩者之間量化制約變化關(guān)系。其中,將Qsm最優(yōu)解和Δγsum最優(yōu)解進行相互對比可知,兩者分別將各自的優(yōu)化目標(biāo)Qsm和Δγsum相對削減了8.6%和35.7%。
此外,結(jié)合圖5 的飛行攻角曲線可以看出,在兩條最優(yōu)軌跡中,作為設(shè)計變量的爬升滑行段(100s以前)飛行攻角存在明顯的不同;而返回段(120s 以后)的飛行攻角在初期較為近似,350s 后逐漸顯現(xiàn)出差異。因此,兩條軌跡的設(shè)計變量對于優(yōu)化目標(biāo)尋優(yōu)的期望值差異主要集中在爬升滑行段和返回段中期以后的飛行過程。
本文面向空基發(fā)射入軌任務(wù),對配置超燃沖壓發(fā)動機的可返回高超聲速巡航飛行器的軌跡多目標(biāo)優(yōu)化問題進行了構(gòu)建、分析和求解,該飛行器須在巡航結(jié)束之后減速返回地面預(yù)定機場水平著陸。首先將高超聲速巡航飛行器的軌跡設(shè)計劃分為巡航段軌跡單目標(biāo)優(yōu)化和剩余段軌跡多目標(biāo)優(yōu)化兩步先后進行,其次對各段飛行軌跡優(yōu)化模型進行了定義,最后應(yīng)用SQP 優(yōu)化算法和DSD-III 多目標(biāo)方法對整個軌跡優(yōu)化問題進行了求解。
研究結(jié)果表明,針對高超聲速巡航飛行器,采用分步進行多目標(biāo)優(yōu)化的策略,可以協(xié)調(diào)處理其不同階段飛行軌跡之間的關(guān)聯(lián)關(guān)系和各飛行階段的不同優(yōu)化需求,并將整個軌跡設(shè)計問題轉(zhuǎn)換為具體模型化表達的優(yōu)化問題。通過采用相應(yīng)的優(yōu)化算法,求解各段飛行軌跡優(yōu)化問題,最終綜合形成完整的飛行軌跡優(yōu)化解。優(yōu)化方法層面,在剩余段軌跡多目標(biāo)優(yōu)化中,采用DSD-III 方法優(yōu)化得到的Pareto 前緣點集無冗余優(yōu)化點,且較為均勻的描繪出了多目標(biāo)優(yōu)化問題的非線性Pareto 前緣;兩個優(yōu)化目標(biāo)均取得了不同幅度的優(yōu)化效果,所有約束條件均滿足要求,可有效求解該軌跡多目標(biāo)優(yōu)化問題。本文的研究成果可為類似的實際工程優(yōu)化設(shè)計提供一種從問題構(gòu)建、到采用相關(guān)優(yōu)化算法進行求解、再到根據(jù)多目標(biāo)優(yōu)化結(jié)果進行制約權(quán)衡設(shè)計的設(shè)計思路和參考范式。