張蕊 孫國慶 楊澤山
(北京航空航天大學(xué)宇航學(xué)院,北京 100191)
臨近空間飛行器再入時(shí),由于需要較長一段時(shí)間在大氣層中飛行,將受到熱流密度、動(dòng)壓、過載的嚴(yán)格約束,以及需要進(jìn)行主動(dòng)姿態(tài)控制,給飛行器材料和內(nèi)部設(shè)施熱絕緣帶來巨大的挑戰(zhàn),所以安全可靠的再入成為重要問題,再入軌跡的優(yōu)化也成為一項(xiàng)關(guān)鍵技術(shù)。
現(xiàn)有的軌跡設(shè)計(jì)方法有工程規(guī)劃法,依據(jù)工程經(jīng)驗(yàn)將軌跡分為數(shù)段,然后分段求解[1]。軌跡優(yōu)化設(shè)計(jì)主要有基于極大值原理的間接法[2-3]和基于非線性規(guī)劃理論的直接法,直接法中的偽譜法[4-5]應(yīng)用最為廣泛,但無論是直接法還是間接法,均采用傳統(tǒng)的尋優(yōu)方法,根據(jù)工程經(jīng)驗(yàn)設(shè)置參數(shù)初值,然后通過試湊法對(duì)這些參數(shù)進(jìn)行調(diào)整。由于目標(biāo)函數(shù)對(duì)待優(yōu)化參數(shù)的微小變化十分敏感,這樣在全局尋優(yōu)時(shí)給軌跡優(yōu)化帶來很大問題,很難對(duì)整個(gè)的解空間進(jìn)行搜索并找到全局最優(yōu)解,即使能找到全局最優(yōu)解,也需花費(fèi)大量的時(shí)間尋優(yōu),而且優(yōu)化結(jié)果的好壞,在很大程度上取決于某些參數(shù)的初始給定值,因而不便于工程應(yīng)用。除了傳統(tǒng)的直接法和間接法之外,還有許多其他優(yōu)化方法應(yīng)用也較廣泛,如動(dòng)態(tài)規(guī)劃方法[6]、快速搜索 隨機(jī)樹 法[7]、遺傳 算法[8]等,這些算法在OPTIMUS軟件平臺(tái)中都有體現(xiàn)。
OPTIMUS軟件平臺(tái)是一種過程集成與多學(xué)科優(yōu)化的軟件平臺(tái),主要功能包括過程自動(dòng)化、設(shè)計(jì)空間開發(fā)、先進(jìn)的概率和統(tǒng)計(jì)方法、數(shù)值優(yōu)化算法等。OPTIMUS已經(jīng)被廣泛應(yīng)用于航空航天領(lǐng)域的結(jié)構(gòu)優(yōu)化、渦輪氣動(dòng)優(yōu)化設(shè)計(jì)[9]、彈道仿真控制系統(tǒng)參數(shù)優(yōu)化[10]等方面。在再入軌跡優(yōu)化問題中,基于OPTIMUS平臺(tái),可以充分應(yīng)用其集成的數(shù)值優(yōu)化算法,快速確定初始猜想值,并選取最佳優(yōu)化方案進(jìn)行求解。
本文針對(duì)再入軌跡優(yōu)化中初值選取和計(jì)算效率的問題,利用OPTIMUS 尋優(yōu)平臺(tái),找出軌跡優(yōu)化的初值點(diǎn),并確定最佳的優(yōu)化算法,再依據(jù)算法編寫優(yōu)化程序,這種計(jì)算機(jī)集成的輔助優(yōu)化設(shè)計(jì)可以取代大量的低效手動(dòng)設(shè)計(jì)。仿真結(jié)果表明,基于OPTIMUS的軌跡優(yōu)化方案計(jì)算快速、尋優(yōu)準(zhǔn)確,且計(jì)算結(jié)果精度高。
飛行器再入過程中,由于各種狀態(tài)變量的絕對(duì)值相差較大,這對(duì)于再入運(yùn)動(dòng)方程的數(shù)值計(jì)算非常不利,為提高計(jì)算精度、增強(qiáng)算法收斂性與快速性,需將實(shí)際再入運(yùn)動(dòng)方程進(jìn)行歸一化處理。因此,在考慮地球旋轉(zhuǎn)的情況下,飛行器再入無量綱運(yùn)動(dòng)方程組為
上述再入運(yùn)動(dòng)方程組是關(guān)于無量綱時(shí)間τ=t/的微分方程組,其中,R0為地球半徑,g0是海平面重力加速度。式(1)~(6)方程組中:r為飛行器質(zhì)心-地心的距離與地球半徑R0的比值;θ和φ分別為再入過程中的經(jīng)度緯度;V為飛行器相對(duì)地球的速度與歸一化變量(圓周速度)的比值;γ為航跡傾角;ψ為航向角;Ω為地球自轉(zhuǎn)角速度,σ為滾轉(zhuǎn)角;L為飛行器的總升力加速度;D為阻力加速度。L、D計(jì)算公式為
式中:ρ為大氣密度;m為飛行器質(zhì)量;Sref為飛行器的參考面積;CL和CD分別為飛行器的升力系數(shù)和阻力系數(shù),其與飛行器的馬赫數(shù)Ma和攻角α有關(guān),也即CL=CL(Ma,α),CD=CD(Ma,α)。
再入軌跡優(yōu)化問題,即尋找最優(yōu)控制律σ(t),使飛行器在再入過程中,滿足要求的初始狀態(tài)、目標(biāo)狀態(tài)和運(yùn)動(dòng)學(xué)方程,并使給定的性能泛函具有極值。本文中設(shè)計(jì)的性能泛函為在給定縱程下使再入橫程最大。
為簡(jiǎn)化問題,再入軌跡優(yōu)化時(shí)采用以下假設(shè)條件:重力加速度g為常值g0;滿足準(zhǔn)平衡滑翔條件γ≈0,γ≈0;由于高度相對(duì)地球半徑很小,令r為常值且r=1;不考慮地球自轉(zhuǎn),將哥氏加速度與牽連加速度影響因素忽略。簡(jiǎn)化后的再入運(yùn)動(dòng)方程為
過程約束考慮氣動(dòng)加熱、過載及動(dòng)壓3種不等式約束。熱流率約束為
動(dòng)壓約束q主要取決于機(jī)體結(jié)構(gòu)和運(yùn)載器表面防熱材料的強(qiáng)度,qmax為允許的最大動(dòng)壓值。
過載約束n主要取決于飛行器結(jié)構(gòu)強(qiáng)度和器載設(shè)備的過載范圍,nmax為允許的最大過載值。
給定縱程下最大橫程問題的描述如下(見圖1):設(shè)飛行器的初始再入點(diǎn)為O點(diǎn),其對(duì)應(yīng)的經(jīng)緯度為θO和φO;要求到達(dá)的目標(biāo)點(diǎn)為T,其對(duì)應(yīng)的經(jīng)緯度為θT和φT;P點(diǎn)為經(jīng)過O點(diǎn)和T點(diǎn)的大圓弧上的點(diǎn),且滿足O點(diǎn)到P點(diǎn)的距離為給定的縱程值SP,其對(duì)應(yīng)的經(jīng)緯度為θP和φP;F點(diǎn)即為給定縱程值SP所對(duì)應(yīng)的最大橫程點(diǎn),其對(duì)應(yīng)的經(jīng)緯度為θF和φF。再入軌跡優(yōu)化問題即為尋求最優(yōu)滾轉(zhuǎn)角控制量σ(t),使飛行器在滿足初始狀態(tài)(O點(diǎn))和給定的縱程SP的情況下,到達(dá)最大橫程點(diǎn)F處,且滿足再入過程約束和再入運(yùn)動(dòng)方程。
圖1 最大橫程示意圖Fig.1 Maximum crossrange diagram
本文主要是基于最優(yōu)控制理論,設(shè)計(jì)滾轉(zhuǎn)角控制量,通過選擇合適的求解方法求解最優(yōu)控制問題,使再入飛行器在給定縱程SP下,獲得最大橫程SF。優(yōu)化目標(biāo)是橫程SF最大,優(yōu)化變量為滾轉(zhuǎn)角σ(|σ|<π/2),性能泛函即為
終點(diǎn)約束條件主要有2個(gè):①達(dá)到預(yù)設(shè)的末端能量;②保證連接∠OPF為直角。
式中:rF為終端高度;VF為終端速度;eF為終端能量,哈密爾頓函數(shù)H寫為
式中:pθ,pφ,pV,pψ為共軛狀態(tài)變量,共軛狀態(tài)方程為
根據(jù)最優(yōu)條件,最優(yōu)滾轉(zhuǎn)角控制律要滿足
根據(jù)運(yùn)動(dòng)的積分以及整理,具體參見文獻(xiàn)[11],可得到滾轉(zhuǎn)角σ最優(yōu)控制律為
在滾轉(zhuǎn)角表達(dá)式(22)中,可以看出分子分母是常數(shù)c1,c2和c3的線性表達(dá)式,由此可以得出,滾轉(zhuǎn)角的值只與其中2個(gè)獨(dú)立變量有關(guān),最終優(yōu)化目標(biāo)也即通過搜索找到2個(gè)獨(dú)立變量的值。假設(shè):c2=sinλ,c3=cosλ。終點(diǎn)約束條件2可以改寫為
至此,給定縱程下最大橫程優(yōu)化問題可表述為:在滿足終點(diǎn)約束條件的前提下,利用尋優(yōu)算法求解最優(yōu)控制參數(shù)c1和λ的值。c1作用域?yàn)椋?1,1],λ作用域?yàn)椋?π,π]。對(duì)實(shí)際飛行器軌跡優(yōu)化而言,可以取
其中K1,K2為權(quán)值,可取為K1=1,K2=1。需要尋找兩個(gè)參數(shù)c1和λ,使Δ(c1,λ)最小。利用OPTIMUS集成的優(yōu)化平臺(tái),可以方便快速地搜索出最優(yōu)參數(shù)c1和λ。
在OPTIMUS軟件平臺(tái)中,集成了多種優(yōu)化算法,分為兩大類:全局優(yōu)化算法和局部?jī)?yōu)化算法。全局優(yōu)化算法主要包括:差分進(jìn)化、自適應(yīng)進(jìn)化、模擬退火、遺傳算法[12]等;局部?jī)?yōu)化算法主要包括:模式搜索法[13]、序列二次規(guī)劃、非線性二次規(guī)劃、廣義簡(jiǎn)約梯度等。全局優(yōu)化算法往往收斂速度較慢,而局部?jī)?yōu)化算法收斂速度較快,但如果初始點(diǎn)選取不夠準(zhǔn)確,局部?jī)?yōu)化算法將有可能使搜索進(jìn)入局部極值,導(dǎo)致結(jié)果不準(zhǔn)確。本文結(jié)合全局優(yōu)化算法和局部?jī)?yōu)化算法各自的優(yōu)缺點(diǎn),提出了基于全局優(yōu)化算法和局部?jī)?yōu)化算法相結(jié)合的混合優(yōu)化的概念,全局優(yōu)化算法采用遺傳算法,局部?jī)?yōu)化算法采用模式搜索法,將遺傳算法得到的結(jié)果作為模式搜索法的初始點(diǎn),并進(jìn)一步精確搜索得到最優(yōu)解。采用混合優(yōu)化算法,不僅縮短了仿真時(shí)間、降低了錯(cuò)誤率,而且提高了優(yōu)化精度,滿足了全局最優(yōu)性。下面主要介紹在OPTIMUS軟件平臺(tái)下優(yōu)化設(shè)計(jì)的詳細(xì)步驟。
為了應(yīng)用OPTIMUS軟件平臺(tái),通過可視化的建模方式,將仿真程序和工作流集成在OPTIMUS平臺(tái)的框架中,首先建立M 文件global.variable.dat,其中包含有待優(yōu)化的控制參數(shù)c1和λ的值。接著將事先編好的數(shù)學(xué)模型仿真程序?qū)隣PTIMUS平臺(tái),搭建系統(tǒng)工作流。在OPTIMUS 中選擇控件,并將所有控件用連接線進(jìn)行連接,搭建OPTIMUS工作流。工作流控件中包含設(shè)計(jì)變量、輸入輸出文本文件、仿真程序模塊、輸出向量、輸出變量以及連接線。搭建工作流的目的,在于以圖像化的方式描述仿真中計(jì)算的邏輯和數(shù)據(jù)流程,仿真工作流在輸入文件中更新輸入?yún)?shù)、進(jìn)行仿真計(jì)算,從結(jié)果文件中讀入計(jì)算結(jié)果。
在手動(dòng)輸入完成后,OPTIMUS 就開始自動(dòng)化的計(jì)算。先對(duì)c1和λ取值,修改global.variable.dat中相應(yīng)的數(shù)值。之后運(yùn)行M 文件搭建的動(dòng)力學(xué)仿真模塊,待程序運(yùn)行結(jié)束后讀取Output.dat中的數(shù)據(jù),根據(jù)事先設(shè)置的算法,自動(dòng)計(jì)算出c1和λ更新值,回到第一步進(jìn)行數(shù)據(jù)刷新,然后重復(fù)這個(gè)工作流程,直至得到符合條件的解。整個(gè)工作不需要人工干預(yù),由計(jì)算機(jī)獨(dú)立完成長時(shí)間的運(yùn)算,提高了工作效率。
完成建模后要進(jìn)行試驗(yàn)設(shè)計(jì),通過較少的試驗(yàn)次數(shù),比較準(zhǔn)確地反映輸入?yún)?shù)和輸出參數(shù)的關(guān)系,并且通過響應(yīng)面模型使結(jié)果可視化。OPTIMUS提供了多種試驗(yàn)設(shè)計(jì)(Design of Experiment,DOE)方法,包括全參數(shù)、隨機(jī)試驗(yàn)設(shè)計(jì)、部分參數(shù)等。試驗(yàn)設(shè)計(jì)的目的在于分析參數(shù)相關(guān)性,為響應(yīng)面的建立提供樣本點(diǎn),保證在設(shè)計(jì)變量變化時(shí)(特別是邊界上)不發(fā)生干涉、參數(shù)不匹配及命令不執(zhí)行等情況。此外還應(yīng)對(duì)網(wǎng)格尺度、計(jì)算步長、差分格式及數(shù)等進(jìn)行匹配調(diào)試,以便在計(jì)算時(shí)間、精度及穩(wěn)定性等方面得到平衡。
優(yōu)化變量c1可行域?yàn)椋?1,1],λ可行域?yàn)椋?π,π],圖2為拉丁超立方樣本點(diǎn)分布,圖3為參數(shù)相關(guān)性。
圖2 拉丁超立方樣本點(diǎn)分布Fig.2 Latin super cubic sample point distribution
圖3 拉丁超立方參數(shù)相關(guān)性Fig.3 Latin super cubic parameters correlation
通過實(shí)驗(yàn)設(shè)計(jì)分析得出:該飛行器軌跡優(yōu)化的目標(biāo)泛函對(duì)選擇的參數(shù)敏感性非常高,參數(shù)響應(yīng)面存在很多局部的急劇變化,因此選擇一種合適的優(yōu)化算法顯得尤其重要。試驗(yàn)點(diǎn)的選擇應(yīng)針對(duì)不同模型有側(cè)重地進(jìn)行考慮,雖然盡可能地將試驗(yàn)點(diǎn)個(gè)數(shù)減少,但前提是體現(xiàn)整個(gè)參數(shù)可行域內(nèi)的分布情況。
選擇相應(yīng)的設(shè)計(jì)實(shí)驗(yàn)方法可以得到響應(yīng)面模型,來模擬一系列輸入與輸出參數(shù)之間的響應(yīng)關(guān)系。由于該飛行器優(yōu)化模型非線性極強(qiáng),對(duì)于參數(shù)的響應(yīng)面波動(dòng)十分劇烈,所以需要相當(dāng)精確的擬合方法對(duì)其進(jìn)行擬合,以減少后續(xù)分析處理中由于模型誤差引入的不確定性。以插值法中的RBF 神經(jīng)網(wǎng)絡(luò)方法為例,形成的插值響應(yīng)面如圖4所示。
圖4 RBF神經(jīng)網(wǎng)絡(luò)方法插值響應(yīng)面Fig.4 RBF neural network method interpolation response surface
響應(yīng)面模型搭建的出發(fā)點(diǎn)在于利用一系列的擬合算法,在若干試驗(yàn)點(diǎn)的基礎(chǔ)上,對(duì)待優(yōu)化的精確響應(yīng)面模型進(jìn)行離散樣本點(diǎn)擬合,用數(shù)理分析統(tǒng)計(jì)的方法代替實(shí)際飛行器模型。雖然舍棄了小部分精度,但換來的是后續(xù)的高效運(yùn)算,極大地減少了數(shù)據(jù)存儲(chǔ)和運(yùn)算量,基于該擬合模型,就可以將輸入?yún)?shù)和輸出參數(shù)關(guān)聯(lián)起來,從而具有較快的響應(yīng)速度和較好的響應(yīng)精度,降低了運(yùn)算成本和響應(yīng)時(shí)間。
完成了上述的工作流的建立、以及試驗(yàn)設(shè)計(jì)和響應(yīng)面模型的搭建之后,就可以選擇OPTIMUS軟件平臺(tái)內(nèi)部自帶的遺傳算法和模式搜索法對(duì)參數(shù)c1和λ進(jìn)行尋優(yōu)計(jì)算,并最終得到參數(shù)的最優(yōu)值為c1=0.784 7,λ=0.524 6。
將OPTIMUS軟件平臺(tái)優(yōu)化出的參數(shù)c1和λ的值,代入到公式(22)的滾轉(zhuǎn)角最優(yōu)控制律中,并編寫相應(yīng)的MATLAB 再入軌跡仿真程序,即得到滿足性能泛函和過程約束的再入軌跡。
MATLAB 仿真用到的再入飛行器參數(shù)為:參考面積149.388m2,飛行器質(zhì)量37 362.9kg,最大熱流率794 005 W/m2,最大動(dòng)壓15 425N/m2,最大過載5gn。攻角α變化范圍為[0°,45°];馬赫數(shù)Ma變化范圍為[3,25];高度變化范圍為[0km,120km]。初始條件為:高度123.104km,經(jīng)度和緯度分別為-122.498°和-33.262 9°,速度為7 625.15m/s,航跡傾角為-1.249 4°,航跡偏角為46.061 7°。終端條件為:高度30.427km,經(jīng)度和緯度分別-62.346 5°,35.270 6°,速度為908.151 6m/s,縱程為9 871.4km。圖5~8是仿真結(jié)果,對(duì)任意的給定縱程,都有左右兩邊的兩個(gè)最大橫程。圖5是高度-速度空間中的再入軌跡,從圖中可以看到滿足熱流率、動(dòng)壓和過載約束;圖6是對(duì)應(yīng)的最優(yōu)滾轉(zhuǎn)角;圖7是航跡傾角;圖8是星下點(diǎn)軌跡(給定縱程下對(duì)應(yīng)的最大橫程有左右2條)。
圖5 再入軌跡圖Fig.5 Entry trajectory
圖6 滾轉(zhuǎn)角控制曲線Fig.6 Bank angle profile
圖7 航跡傾角變化曲線Fig.7 Flight path angle profile
圖8 星下點(diǎn)軌跡Fig.8 Trajectory of ground track
本文基于OPTIMUS優(yōu)化軟件平臺(tái),提出了臨近空間飛行器的再入軌跡優(yōu)化方案。與傳統(tǒng)的優(yōu)化方法不同,先利用OPTIMUS平臺(tái)優(yōu)化求出初值,并選出最佳優(yōu)化算法,再編寫優(yōu)化程序進(jìn)行求解。OPTIMUS平臺(tái)在軌跡優(yōu)化中的應(yīng)用,可以快速地搭建起整個(gè)數(shù)據(jù)流程,通過應(yīng)用多種不同的優(yōu)化算法,對(duì)同一個(gè)問題進(jìn)行求解分析,可以更加全面地了解軌跡優(yōu)化問題各個(gè)環(huán)節(jié)的細(xì)節(jié),并獲得優(yōu)化問題的初始猜想,然后針對(duì)具體問題設(shè)計(jì)優(yōu)化方案,所開發(fā)的算法程序具有計(jì)算高效、便捷準(zhǔn)確的特點(diǎn)。
(References)
[1]趙漢元,飛行器再入動(dòng)力學(xué)和制導(dǎo)[M].長沙:國防科學(xué)技術(shù)大學(xué)出版社,1997
Zhao Hanyuan.Flight vehicle entry dynamics and guidance[M].Changsha:University of National Defense Science and Technology Press,1997(in Chinese)
[2]Vinh N X.Optimal trajectories in atmospheric flight[M].NewYork:Elsevier Scientific Publishing Company,1981
[3]Lu P,Xue S B.Rapid generation of accurate entry landing footprints[J].Journal of Guidance,Control and Dynamics,2010,33(3):756-767
[4]雍恩米,唐國金,陳磊.基于Gauss偽譜方法的高超聲速飛行器再入軌跡快速優(yōu)化[J].宇航學(xué)報(bào),2008,29(6):1766-1772
Yong Eenmi,Tang Guojin,Chen Lei.Rapid trajectory optimization for hypersonic reentry vehicle via Guass pseudospectral method[J].Journal of Astronautics,2008,29(6):1766-1772(in Chinese)
[5]宗群,田栢苓,竇立謙.基于Gauss偽譜法的臨近空間飛行器上升段軌跡優(yōu)化[J].宇航學(xué)報(bào),2010,31(7):1775-1781
Zong Qun,Tian Bailing,Dou Liqian.Ascent phase trajectory optimization for near space vehicle based on Gauss pseudospectral method[J].Journal of Astronautics,2010,31(7):1775-1781(in Chinese)
[6]John T B.Survey of numerical methods for trajectory optimization[J].Journal of Guidance,Control and Dynamics,1998,21(3):193-207
[7]周克強(qiáng),高曉光,白奕.改進(jìn)的快速擴(kuò)展隨機(jī)樹在航跡規(guī)劃中的應(yīng)用[J].系統(tǒng)工程與電子技術(shù),2006,28(10):1538-1540
Zhou Keqiang,Gao Xiaoguang,Bai Yi.Trajectory planning using improved rapidly exploring random tree[J].System Engineering and Electronics,2006,28(10):1538-1540(in Chinese)
[8]周浩,陳萬春.高超聲速飛行器滑行航跡優(yōu)化[J].北京航空航天大學(xué)學(xué)報(bào),2006,32(5):513-517
Zhou Hao,Chen Wanchun.Optimization of glide trajectory for a hypersonic vehicle[J].Journal of Beijing University of Aeronautics and Astronautics,2006,32(5):513-517(in Chinese)
[9]嚴(yán)俊峰,吳寶元,逯婉若.基于OPTIMUS的渦輪氣動(dòng)優(yōu)化設(shè)計(jì)[J].火箭推進(jìn),2008,34(2):13-17
Yan Junfeng,Wu Baoyuan,Lu Wanruo.Turbo aerodynamics optimization design based on OPTIMUS flat[J].Journal of Rocket Propulsion,2008,34(2):13-17(in Chinese)
[10]周健,姜曉菡,周喻虹.基于OPTIMUS的控制系統(tǒng)參數(shù)優(yōu)化計(jì)算[J].彈箭與制導(dǎo)學(xué)報(bào),2008,28(3):13-18
Zhou Jian,Jiang Xiaohan,Zhou Yuhong.Optimization of parameters for control system based on OPTIMUS[J].Journal of Projectiles,Rockets,Missiles and Guidance,2008,28(3):13-18(in Chinese)
[11]Lu P.Entry trajectory optimization with analytical feedback bank angle law,AIAA-2008-7268[R].Washington:AIAA,2008
[12]袁麟博,章衛(wèi)國,李廣文.一種基于遺傳算法-模式搜索法的無人機(jī)路徑規(guī)劃[J].彈箭與制導(dǎo)學(xué)報(bào),2009,29(3):279-282
Yuan Linbo,Zhang Weiguo,Li Guangwen.A path planning for UAV based on genetic-pattern searching algorithm[J].Journal of Projectiles,Rockets,Missiles and Guidance,2009,29(3):279-282(in Chinese)
[13]張煜東,吳樂南,王水花.基于遺傳算法與模式搜索的混合優(yōu)化算法[J].南京信息工程大學(xué)學(xué)報(bào),2012,4(1):34-39
Zhang Yudong,Wu Lenan,Wang Shuihua.A hybrid optimization method based on genetic algorithm and pattern search[J].Journal of Nanjing University of Information Science & Technology,2012,4(1):34-39(in Chinese)