任洋, 曹林平, 國海峰, 丁達(dá)理

(空軍工程大學(xué) 航空航天工程學(xué)院, 陜西 西安 710038)

0 引言

近年來,物理規(guī)劃方法[1]在解決多目標(biāo)優(yōu)化問題上顯現(xiàn)出較大的優(yōu)勢(shì),特別是它能夠根據(jù)設(shè)計(jì)者的偏好,得出合理的折中解。胡丹等[2]利用物理規(guī)劃解決了數(shù)據(jù)中繼衛(wèi)星鏈路多目標(biāo)配置決策問題,但對(duì)各項(xiàng)指標(biāo)的產(chǎn)生辦法沒有做出解釋。龍瓊等[3]利用物理規(guī)劃對(duì)路徑誘導(dǎo)問題進(jìn)行了定量描述,通過與A*算法的結(jié)合將路徑?jīng)Q策與駕駛者的偏好結(jié)合起來。雍恩米等[4]運(yùn)用物理規(guī)劃解決了高超聲速滑翔飛行器再入軌跡優(yōu)化問題,但其運(yùn)動(dòng)學(xué)模型只考慮了二維平面,沒有考慮地球自轉(zhuǎn),相對(duì)簡單,且采用遺傳算法作為優(yōu)化算法,若要達(dá)到較高尋優(yōu)精度,所需運(yùn)算時(shí)間相對(duì)較長。

本文以高超聲速滑翔飛行器三維再入段軌跡優(yōu)化為研究對(duì)象,運(yùn)用了hp自適應(yīng)偽譜法,在物理規(guī)劃的框架下對(duì)多目標(biāo)再入軌跡進(jìn)行了尋優(yōu);針對(duì)hp自適應(yīng)偽譜法的尋優(yōu)特點(diǎn),對(duì)物理規(guī)劃框架的偏好函數(shù)進(jìn)行了改進(jìn),增大了其適用范圍;對(duì)不同偏好的再入段最優(yōu)軌跡進(jìn)行了仿真對(duì)比,為優(yōu)化方案選擇提供了參考。

1 三自由度再入軌跡數(shù)學(xué)模型

飛行器再入軌跡數(shù)學(xué)模型包括飛行器三自由度動(dòng)力學(xué)模型和氣動(dòng)力模型。

1.1 飛行器動(dòng)力學(xué)模型

考慮地球?yàn)樾D(zhuǎn)圓球的三自由度再入運(yùn)動(dòng)模型為:

(1)

式中,r為地心距;θ為經(jīng)度;φ為緯度;V為速度;γ為航跡角;ψ為航向角;ω,g,β,L,D分別為地球自轉(zhuǎn)角速度、引力加速度、側(cè)傾角、升力和阻力。

1.2 氣動(dòng)力模型

在高超聲速飛行器再入過程中,要通過臨近空間穿過大氣層到達(dá)指定空域,而臨近空間大氣環(huán)境變化復(fù)雜,需對(duì)氣動(dòng)模型進(jìn)行較為精確的擬合。飛行器受到的升力L和阻力D的計(jì)算公式如下:

(2)

式中,S為氣動(dòng)參考面積;ρ為大氣密度;升力系數(shù)CL(α,Ma)和阻力系數(shù)CD(α,Ma)的表達(dá)式為:

(3)

式(3)中的系數(shù)l(j),d(j)通過對(duì)X-33飛行數(shù)據(jù)進(jìn)行擬合得到[5],如表1所示。

表1 l(j),d(j)的數(shù)值Table 1 The value of l(j),d(j)

本文中的大氣密度計(jì)算公式用指數(shù)大氣模型來近似代替:

ρ=ρ0e-h/H

(4)

式中,ρ0為海平面的大氣密度;高度h=r-R0,R0為地球平均半徑;H為標(biāo)量高度系數(shù)。

1.3 約束條件

(1)終端約束

本文中再入終端約束包括高度、經(jīng)度、緯度和速度約束,相關(guān)約束條件為:

h(tf)=hf,θ(tf)=θf

φ(tf)≥φf,V(tf)≥Vf

(2)熱流、過載、動(dòng)壓約束

(5)

1.4 優(yōu)化目標(biāo)

f1=min(tf)

f2=max(δ)=max|φf-φ0|

1.5 軌跡優(yōu)化問題描述

軌跡優(yōu)化問題實(shí)際上是最優(yōu)控制問題,在滿足各種約束條件下,尋求使某個(gè)指標(biāo)最優(yōu)的控制量。高超聲速飛行器再入軌跡優(yōu)化可以歸結(jié)為如下最優(yōu)控制問題:

(6)

式中,狀態(tài)變量X(t)=[r(t),θ(t),φ(t),V(t),γ(t),ψ(t)];控制變量U(t)=[α(t),β(t)];路徑約束為C[X(t),U(t),t]。

偽譜法在求解單目標(biāo)軌跡優(yōu)化問題時(shí)顯現(xiàn)出很大的優(yōu)勢(shì)。本文采用hp自適應(yīng)偽譜法作為軌跡尋優(yōu)方法,hp自適應(yīng)偽譜法的特點(diǎn)是運(yùn)算所需時(shí)間較短,能夠通過決定在某個(gè)網(wǎng)格時(shí)間區(qū)間中的狀態(tài)量和控制量是否滿足特定的偏差來調(diào)整當(dāng)前的計(jì)算方式,以達(dá)到提高算法精度的目的。如果沒有滿足精度要求,則配置點(diǎn)的數(shù)目及分配方式要通過增加多項(xiàng)式的階數(shù)或者更新網(wǎng)格來達(dá)到精度要求,具體算法參見文獻(xiàn)[6-7]。

2 物理規(guī)劃方法概述

2.1 設(shè)計(jì)思路

物理規(guī)劃[1]是Messac提出的一種處理多目標(biāo)優(yōu)化設(shè)計(jì)問題的有效方法,是處理多目標(biāo)優(yōu)化問題的新框架。其基本思路是:通過求解單目標(biāo)優(yōu)化問題,確定區(qū)間邊界,從而建立偏好函數(shù),將不同物理意義的各種設(shè)計(jì)目標(biāo)轉(zhuǎn)換為具有相同數(shù)量級(jí)的無量綱的綜合滿意度目標(biāo),再對(duì)綜合滿意度目標(biāo)進(jìn)行優(yōu)化,尋求使其最優(yōu)的設(shè)計(jì)點(diǎn),即為最優(yōu)解。

2.2 數(shù)學(xué)模型

物理規(guī)劃的決策模型為:

(7)

(8)

(9)

式中，gi(k-1)

(10)

λik=gik-gi(k-1)

(11)

(12)

(13)

(14)

(15)

區(qū)間端點(diǎn)信息確定方法如下:

(16)

構(gòu)造偏好函數(shù)時(shí),需要驗(yàn)證其是否滿足二階導(dǎo)數(shù)嚴(yán)格為正的特性,即對(duì)區(qū)間邊界處的二階導(dǎo)數(shù)進(jìn)行驗(yàn)證,可由式(17)計(jì)算得出。

(17)

其中:

(18)

(19)

由于0≤ξik≤1,所以只需要驗(yàn)證a,b是否為正即可[2]。

2.3 算法改進(jìn)

(20)

由1.4節(jié)中確定的優(yōu)化目標(biāo)可知,偏好僅為Class1-S型和Class2-S型,其優(yōu)化模型的綜合滿意度目標(biāo)為:

(21)

3 算法流程

首先運(yùn)用hp自適應(yīng)偽譜法求解出單目標(biāo)情況下的最優(yōu)軌跡,利用得到的各目標(biāo)參數(shù)確定物理規(guī)劃偏好區(qū)間邊界,建立物理規(guī)劃框架。然后將物理規(guī)劃優(yōu)化模型的代價(jià)函數(shù)作為hp自適應(yīng)偽譜法的尋優(yōu)目標(biāo)進(jìn)行尋優(yōu),即將每次循環(huán)得到的各單目標(biāo)值帶入物理規(guī)劃框架,得到綜合滿意度目標(biāo),對(duì)綜合滿意度目標(biāo)進(jìn)行尋優(yōu),得到最優(yōu)解。這樣就將多個(gè)目標(biāo)軌跡優(yōu)化問題轉(zhuǎn)換為一個(gè)反映設(shè)計(jì)者偏好的單目標(biāo)問題。具體流程如圖1所示。

圖1 算法流程圖Fig.1 Flow chart of the algorithm

4 仿真結(jié)果及分析

4.1 初始條件

取飛行器質(zhì)量m=907 kg，氣動(dòng)參考面積S=0.48 m2，初始再入點(diǎn)高度h=80 km，初始經(jīng)度θ0=0°，初始緯度φ0=0°，初始速度V0=7000 m/s，初始航跡角γ0=0°，初始航向角ψ0=0°，設(shè)定最終高度hf=20 km，最終經(jīng)度θf=78°,最終緯度φf≥10°,最終速度Vf≥2000 m/s。

狀態(tài)變量及控制變量的約束條件如下:

4.2 偏好結(jié)構(gòu)的確定

為確定物理規(guī)劃問題的偏好結(jié)構(gòu),首先求解和分析單目標(biāo)最優(yōu)軌跡。運(yùn)用hp自適應(yīng)偽譜法,分別以各單目標(biāo)為優(yōu)化目標(biāo)尋優(yōu),可得到當(dāng)滿足某單目標(biāo)最優(yōu)時(shí)其余目標(biāo)的數(shù)據(jù)如表2所示。

表2 單目標(biāo)最優(yōu)時(shí)其余目標(biāo)對(duì)應(yīng)值Table 2 Data corresponding to the rest when single objective is optimum

由于再入過程中各中間變量存在著較強(qiáng)的耦合關(guān)系,不存在使所有目標(biāo)同時(shí)最優(yōu)的設(shè)計(jì)點(diǎn)。通過對(duì)其整體分析可以看出,追求單個(gè)目標(biāo)最優(yōu)時(shí),容易出現(xiàn)其他指標(biāo)太過局限的情況,實(shí)用性不大。

根據(jù)相關(guān)約束的邊界條件及表2計(jì)算結(jié)果,確定三組區(qū)間邊界不同的偏好結(jié)構(gòu),記為偏好1,2和3,分別對(duì)應(yīng)于期望橫程較大、再入時(shí)間較短和總氣動(dòng)加熱較小的三組偏好,如表3所示。

表3 多目標(biāo)軌跡優(yōu)化的偏好結(jié)構(gòu)Table 3 Preference structure of multi-objective trajectory optimization

4.3 優(yōu)化結(jié)果及分析

為作對(duì)比,本文分別采用hp自適應(yīng)偽譜法和文獻(xiàn)[4]中算法對(duì)該問題進(jìn)行尋優(yōu),設(shè)定精度為10-5。

(1)優(yōu)化結(jié)果對(duì)比

根據(jù)表3中的三組偏好結(jié)構(gòu),采用兩種算法,得到的三組目標(biāo)函數(shù)分別對(duì)應(yīng)橫程較大、再入時(shí)間較短、總氣動(dòng)加熱較小。

本文算法所需運(yùn)算時(shí)間分別為12.32 s,12.20 s,11.79 s。三組偏好函數(shù)為:

f1=(1854 s，26.0°，6.48×105kJ/m2，844 kW/m2)

f2=(1601 s，19.0°，6.49×105kJ/m2，794 kW/m2)

f3=(1660 s，21.1°，6.05×105kJ/m2，859 kW/m2) 文獻(xiàn)[4]算法所需運(yùn)算時(shí)間分別為66.47 s,69.55 s,56.06 s。三組偏好函數(shù)為:

f1=(1842 s，25.9°，6.48×105kJ/m2，846 kW/m2)

f2=(1596 s，18.9°，6.48×105kJ/m2，789 kW/m2)

f3=(1662 s，21.1°，6.07×105kJ/m2，861 kW/m2)

對(duì)比可知,兩種算法所得各組目標(biāo)函數(shù)均能體現(xiàn)出設(shè)計(jì)者的偏好,且各項(xiàng)設(shè)計(jì)指標(biāo)都在可接受的范圍內(nèi),是各目標(biāo)之間相互妥協(xié)折中的反映,避免了因追求單個(gè)目標(biāo)最優(yōu)而出現(xiàn)的其他指標(biāo)太過局限的情況。但本文算法在所需尋優(yōu)時(shí)間上有明顯優(yōu)勢(shì),更接近于在線實(shí)時(shí)規(guī)劃；驗(yàn)證了改進(jìn)的物理規(guī)劃方法的有效性,使其可在不改變約束邊界值的情況下,增大適用范圍。

(2)再入過程分析

對(duì)本文算法所得優(yōu)化結(jié)果進(jìn)行分析,最優(yōu)軌跡及相應(yīng)的控制變量、過程約束量如圖2～圖9所示。通過對(duì)比分析可以發(fā)現(xiàn):

①在要求總氣動(dòng)加熱較小時(shí),高超聲速滑翔再入軌跡的振蕩幅值較大,即通過大幅振蕩來減少動(dòng)能和勢(shì)能向熱能的轉(zhuǎn)換,但此時(shí)對(duì)飛行器機(jī)體結(jié)構(gòu)要求較高。

②當(dāng)設(shè)計(jì)者要求較短時(shí)間再入時(shí),飛行器為快速再入到指定高度,在很長一段時(shí)間內(nèi)保持了較大的再入速度;末端迎角快速減小,導(dǎo)致飛行末端動(dòng)壓較大,同時(shí)引起其他變量的一系列大幅調(diào)整;對(duì)比三組偏好可發(fā)現(xiàn),隨著末端速度的減小,末端動(dòng)壓逐漸回歸到較好的范圍內(nèi)。

③初始下降段一般采用較大的迎角減速,對(duì)比圖6可以看出,要求總氣動(dòng)加熱較小時(shí),初始段迎角平均值較其他偏好大,以增大振蕩幅度;末端迎角大幅變小,以較快進(jìn)入指定高度。

④初始段傾側(cè)角較大,以快速改變?cè)偃敕较?增大橫程;隨著速度的減小,傾側(cè)角逐漸減小,這是由于要兼顧縱程和再入時(shí)間,不能因?yàn)樽非髾M程而使縱程過小或者再入時(shí)間太長;末端由于迎角的調(diào)整,傾側(cè)角有相應(yīng)的振蕩變化,使末端再入過程滿足過載、熱流等約束。

圖2 高度變化曲線Fig.2 Altitude variation curve

圖3 速度變化曲線Fig.3 Velocity variation curve

圖4 橫程變化曲線Fig.4 Cross-range variation curve

圖5 迎角變化曲線Fig.5 AOA variation curve

圖6 傾側(cè)角變化曲線Fig.6 Bank angle variation curve

圖7 過載變化曲線Fig.7 Load variation curve

圖8 動(dòng)壓變化曲線Fig.8 Dynamic pressure variation curve

圖9 熱流密度變化曲線Fig.9 Variation curve of heat flow density

5 結(jié)束語

本文對(duì)高超聲速滑翔式再入飛行器三自由度最優(yōu)軌跡進(jìn)行了研究,在物理規(guī)劃的框架下,對(duì)期望橫程較大、再入時(shí)間較短以及總氣動(dòng)加熱較小的三組偏好進(jìn)行了仿真,驗(yàn)證了改進(jìn)的物理規(guī)劃方法與偽譜法結(jié)合求解多目標(biāo)優(yōu)化問題的有效性,得到了反映設(shè)計(jì)者偏好且實(shí)用性更強(qiáng)的三組折中解,所需優(yōu)化時(shí)間相對(duì)較短,具有實(shí)時(shí)應(yīng)用潛力,可為再入軌跡的優(yōu)化選擇提供一定借鑒。

參考文獻(xiàn)：

[1] Messac A.Physical programming-effective optimization for computational design[J].AIAA Journal,1996,34(1):149-158.

[2] 胡丹,閆忠文.基于物理規(guī)劃的數(shù)據(jù)中繼衛(wèi)星鏈路配置決策方法[J].中國空間科學(xué)技術(shù),2009,10(5):54-60.

[3] 龍瓊,胡列格,張蕾,等.基于物理規(guī)劃的路徑誘導(dǎo)方法[J].中南大學(xué)學(xué)報(bào)(自然科學(xué)版),2012,43(8):3287-3293.

[4] 雍恩米,唐國金,陳磊.基于物理規(guī)劃的高超聲速飛行器滑翔式再入軌跡優(yōu)化[J].航空學(xué)報(bào),2008,29(5):1091-1097.

[5] Bollino Kevin P,Michael Ross,Doman David D.Optimal nonlinear feedback guidance for reentry vehicles[R].AIAA-2006-6074,2006.

[6] Darby C L,Hager W W,Rao A V.An hp-adaptive pseudospectral method for solving optimal control problems[J].Optimal Control Applications and Methods,2011,32(4):476-502.

[7] 國海峰,黃長強(qiáng),丁達(dá)理,等.多約束條件下高超聲速導(dǎo)彈再入軌跡優(yōu)化研究[J].彈道學(xué)報(bào),2013,25(1):10-15.