丁洪波，蔡 洪

(國(guó)防科技大學(xué)航天與材料工程學(xué)院，長(zhǎng)沙 410073)

基于特征參數(shù)的吸氣式高超聲速飛行器上升段軌跡優(yōu)化①

丁洪波，蔡 洪

(國(guó)防科技大學(xué)航天與材料工程學(xué)院，長(zhǎng)沙 410073)

針對(duì)吸氣式高超聲速飛行器上升段軌跡優(yōu)化問題，提出并研究了基于特征參數(shù)的軌跡優(yōu)化方法。首先，建立了吸氣式高超聲速飛行器動(dòng)力學(xué)模型，給出了氣動(dòng)力和推力模型。根據(jù)上升段軌跡特性，建立了基于指數(shù)函數(shù)和多項(xiàng)式的控制變量的取值模型。該模型取決于若干特征參數(shù)，從而將一個(gè)求解泛函的最優(yōu)控制問題轉(zhuǎn)化為求解特征參數(shù)的非線性規(guī)劃問題，并采用序列二次規(guī)劃算法求解。針對(duì)初值敏感性，提出了基于遺傳算法的初值選取方法，以及基于物理意義的手動(dòng)選取方法。

吸氣式高超聲速飛行器;上升段;特征參數(shù);軌跡優(yōu)化;初值選取

0 引言

吸氣式高超聲速飛行器(AHFV)是未來最有可能實(shí)現(xiàn)高超聲速單級(jí)入軌和高超聲速巡航飛行的飛行器，它代表了未來航空航天領(lǐng)域的研究發(fā)展方向。這類飛行器所采用的超聲速燃燒沖壓發(fā)動(dòng)機(jī)性能與馬赫數(shù)、高度等飛行狀態(tài)密切相關(guān)，且發(fā)動(dòng)機(jī)余量較小。因此，必須對(duì)上升段飛行軌跡進(jìn)行優(yōu)化。

Buren和Mease［1］建立了單級(jí)入軌飛行器動(dòng)力學(xué)模型，采用了軌跡優(yōu)化軟件OTIS進(jìn)行了優(yōu)化，得到上升段飛行軌跡。Corban和Calise［2］采用能態(tài)近似法和奇異攝動(dòng)方法，在若干假設(shè)簡(jiǎn)化條件下，設(shè)計(jì)了一種次優(yōu)的軌跡優(yōu)化與制導(dǎo)方法。Lu Ping［3］采用動(dòng)態(tài)逆的思想設(shè)計(jì)了空天飛機(jī)的上升段軌跡。該飛行器的飛行任務(wù)是進(jìn)入地球軌道，發(fā)動(dòng)機(jī)采用吸氣式發(fā)動(dòng)機(jī)和火箭發(fā)動(dòng)機(jī)組合方式。吸氣式發(fā)動(dòng)機(jī)主要在稠密大氣中作用。吸氣式發(fā)動(dòng)機(jī)的風(fēng)門系數(shù)采用三次多項(xiàng)式擬合，并采用參數(shù)優(yōu)化的方法求解未知量。Keshmiri［4］建立了簡(jiǎn)化形式的吸氣式飛行器上升段動(dòng)力學(xué)方程，借鑒直接法的思想采用MATLAB工具箱進(jìn)行了軌跡優(yōu)化。Park［5］以及宗群［6］等使用偽譜法進(jìn)行吸氣式超聲速飛行器的上升段軌跡優(yōu)化，設(shè)計(jì)了優(yōu)化策略，并進(jìn)行了縱向平面內(nèi)的仿真計(jì)算。佘智勇等［7］采用間接法將上升段軌跡優(yōu)化問題轉(zhuǎn)化為兩點(diǎn)邊值問題，然后采用伴隨方法進(jìn)行求解。

綜合前人的研究成果，在20世紀(jì)90年代，學(xué)者主要采用間接法，并借鑒奇異攝動(dòng)的思想研究和求解AHFV的軌跡優(yōu)化。近幾年來，隨著軌跡優(yōu)化方法的發(fā)展以及計(jì)算機(jī)硬件水平的提高，直接法特別是偽譜法得到了廣泛應(yīng)用，并取得較好效果。

本文旨在對(duì)AHFV上升段軌跡進(jìn)行優(yōu)化設(shè)計(jì)，提出一種不同的軌跡優(yōu)化方法。其特點(diǎn)是充分考慮AHFV上升段軌跡的特性，著重從工程應(yīng)用的角度分析控制量的取值特性，在此基礎(chǔ)上建立基于特征參數(shù)法的優(yōu)化模型，并研究了相應(yīng)的優(yōu)化算法。

1 吸氣式飛行器動(dòng)力學(xué)方程

AHFV一般具有升力體外形，采用沖壓發(fā)動(dòng)機(jī)作為推力裝置。因此，它在飛行過程中會(huì)受到推力和氣動(dòng)力的作用。AHFV的縱向運(yùn)動(dòng)方程為［8］

式中u為無量綱化速度;γ為當(dāng)?shù)厮俣葍A角;z為無量綱化地心距;θ為航程;為無量綱化質(zhì)量和為無量綱化阻力和升力為無量綱化推力;α為攻角。

實(shí)際的阻力和升力為D和L，則

式中q為動(dòng)壓;Sref為飛行器有效參考面積。

因?yàn)闅鈩?dòng)系數(shù)與攻角α有關(guān)，因此空氣動(dòng)力主要由α和q決定。

實(shí)際推力大小為

式中為最大推力系數(shù);st為發(fā)動(dòng)機(jī)風(fēng)門系數(shù)，0≤st≤1;Se為引擎進(jìn)氣道面積。

推力大小主要由st和q決定。

綜上所述，AHFV的控制量有2個(gè):一個(gè)是風(fēng)門系數(shù)st，它主要控制推力大小，當(dāng)st=0時(shí)，風(fēng)門關(guān)閉，推力為零，當(dāng)st=1時(shí)，風(fēng)門最大，推力為當(dāng)前狀態(tài)下所能得到的最大推力;另一個(gè)控制量是攻角α，它控制升力和阻力的大小。軌跡優(yōu)化的目的是選取合適的st和α，在滿足約束條件前提下，使性能指標(biāo)最大或最小。

2 基于特征參數(shù)的軌跡優(yōu)化模型

2.1 攻角取值模型

AHFV上升段軌跡與火箭上升段飛行軌跡有類似之處，所以可借鑒火箭上升段大氣層內(nèi)轉(zhuǎn)彎部分軌跡設(shè)計(jì)方法?；鸺仙喂ソ且话悴捎孟率剑?］:

攻角變化規(guī)律參見圖1。開始迅速地達(dá)到負(fù)極值，然后絕對(duì)值開始變小，以指數(shù)速率趨于零，趨于零的速度由參數(shù)a決定。

火箭一般采用垂直起飛方式，需采用負(fù)攻角轉(zhuǎn)彎。AHFV在發(fā)動(dòng)機(jī)工作階段攻角必須大于零，同時(shí)考慮到其初始攻角不一定為零，AHFV攻角模型可采用下述表達(dá)式:

式中tc為時(shí)間的無量綱量。

圖1給出了火箭與AHFV的攻角模型，二者的區(qū)別主要是火箭攻角由零減小，然后又增大，趨近于零;AHFV的攻角則在初始時(shí)刻有一不小于零的值，然后增大至最大值，然后又趨近于初始值。

圖1 火箭與AHFV攻角模型比較Fig.1 Attack angle models for rocket and AHFV

2.2 風(fēng)門系數(shù)取值模型

對(duì)于AHFV而言，風(fēng)門系數(shù)大小可調(diào)整，風(fēng)門系數(shù)可取為三次多項(xiàng)式模型［3］。

假設(shè)初始時(shí)刻t0和終端時(shí)刻tf的油門系數(shù)大小及其變化率分別為cT0、cTf、˙cT0、˙cTf，則

根據(jù)上式，可求出三次多項(xiàng)式的系數(shù)。

式(6)中待確定的參數(shù)為cT0、cTf、˙cT0、˙cTf，以及終端時(shí)刻tf。

需注意的是風(fēng)門系數(shù)的大小是有約束的，其值在［0，1］之間。針對(duì)這個(gè)約束，一般的處理方式是直接采用式(7)［3］。

本文通過對(duì)式(7)的分析，將針對(duì)cT的約束轉(zhuǎn)化為對(duì)待優(yōu)化參數(shù)cT0、cTf、˙cT0、˙cTf的約束，使問題得到簡(jiǎn)化。

式(7)等效于

即cT的最大或最小值均為0和1之間。對(duì)式(5)求導(dǎo)，即

當(dāng)a3=0時(shí)

對(duì)應(yīng)的極值為2

當(dāng)a3≠0，且-12a3a1≥0 時(shí)，cT在［0，1］之間存在2個(gè)極值。

當(dāng)a3≠0，且-12a3a1≥0 時(shí)，cT在［0，1］之間無極值。

綜上所述，要滿足式(7)，只需有式(13)即可:

其中，cT4=cTf。

2.3 過程約束的處理方法

上升段飛行過程中，需考慮過載N以及動(dòng)壓q等過程約束。以動(dòng)壓約束q為例，本文采用如下方法處理過程約束。

引入變量ζ作為一個(gè)新的狀態(tài)變量，使

用這種方法可將飛行過程中超出動(dòng)壓約束的點(diǎn)記錄下來。于是，如果沒有超過過程約束，則ζ(tf)=0，否則 ζ(tf)＞0。

如果性能指標(biāo)為燃料最省，也即

考慮動(dòng)壓約束后，引入新的性能指標(biāo):

其中，ν為常數(shù)，稱為懲罰因子，可取為大于1的值。

2.4 優(yōu)化模型與求解算法

至此，已將控制量α、st隨時(shí)間的變化規(guī)律用若干具有一定物理意義的特征參數(shù)表示出來。根據(jù)初始狀態(tài)以及控制量的值，即可對(duì)動(dòng)力學(xué)模型積分，得到飛行軌跡。根據(jù)狀態(tài)量和控制量的值，可判斷軌跡是否滿足約束條件，并結(jié)合性能指標(biāo)函數(shù)的值，采用迭代尋優(yōu)的方法來求取最佳的特征參數(shù)。這就是基于特征參數(shù)的積分方法的基本思想。

該方法的優(yōu)化模型如下:

式中f代表動(dòng)力學(xué)模型;x為狀態(tài)變量;φ為終端約束條件。

式(16)描述的優(yōu)化問題的實(shí)質(zhì)是一個(gè)非線性規(guī)劃問題，可采用序列二次規(guī)劃(SQP)算法求解。與其他優(yōu)化算法相比，SQP具有收斂性好、計(jì)算效率高、邊界搜索能力強(qiáng)等特點(diǎn)，且能有效處理含約束問題。但該算法需首先給定參數(shù)的迭代初值，迭代初值的好壞直接影響著算法的性能，如果初始不好，還可能導(dǎo)致得不到優(yōu)化結(jié)果。

3 優(yōu)化初值的求取方法

本章提出2種優(yōu)化初值的求取方法，包括遺傳算法和基于物理意義的手動(dòng)選取方法。

3.1 遺傳算法

遺傳算法(GA)的優(yōu)點(diǎn)在于良好的全局尋優(yōu)性能，且不存在初值敏感問題。這里提出“廣義”可行解的概念，旨在解決遺傳算法搜尋可行解困難的問題。

以二維解空間為例進(jìn)行說明，如圖2所示［10-11］。

圖2 “廣義”可行解的概念Fig.2 Concept of extended feasible solution

圖2中，S1～S5代表了原優(yōu)化問題的可行解集域，如果這些區(qū)域在整個(gè)解空間內(nèi)分布較小且相隔較遠(yuǎn)，遺傳算法的隨機(jī)搜索機(jī)制難以保證在有限次遺傳代數(shù)內(nèi)使優(yōu)化結(jié)果收斂于可行解集域。為了給后續(xù)的SQP算法提供1組較好的迭代初值，可適當(dāng)放寬約束條件，如將等式約束改為不等式約束、將不等式約束的范圍擴(kuò)大，于是原優(yōu)化問題轉(zhuǎn)化為中間優(yōu)化問題，其可行解集域?yàn)镋，E內(nèi)任一點(diǎn)即為原優(yōu)化問題的“廣義”可行解。顯然，S1～S5均包含于E內(nèi)，E的大小由約束條件的放寬程度決定。

但是，利用遺傳算法求解上述優(yōu)化問題存在明顯缺陷:

(1)當(dāng)解空間較大時(shí)，其概率性搜索的性質(zhì)無法保證以足夠的精度收斂到最優(yōu)解;

(2)對(duì)約束的處理不如SQP靈活，內(nèi)含兩點(diǎn)邊值問題，導(dǎo)致解集域不連續(xù)，故可行解的搜尋存在困難。

考慮到SQP算法具有很強(qiáng)的局部尋優(yōu)能力，在一定程度上彌補(bǔ)了遺傳算法精度不高的缺點(diǎn)。因此，可采用“GA+SQP”的組合優(yōu)化算法求解軌跡優(yōu)化問題，即首先利用遺傳算法獲得“廣義”的最佳可行解，然后將此可行解作為SQP算法的迭代初值進(jìn)行尋優(yōu)計(jì)算。

3.2 基于物理意義的手動(dòng)選取方法

在第2章對(duì)攻角模型及其表達(dá)式中各個(gè)參數(shù)含義進(jìn)行了分析，可看到各個(gè)參數(shù)具有明確的含義?？筛鶕?jù)各個(gè)參數(shù)對(duì)攻角模型的影響以及攻角對(duì)飛行軌跡的影響，采用手動(dòng)方法(MM)來確定各個(gè)參數(shù)的初值。

下面分析各攻角模型中各個(gè)參數(shù)對(duì)軌跡的影響，見圖3。

圖3 α0對(duì)飛行軌跡的影響Fig.3 Effect of α0on trajectory

從圖3可看到，增大α0會(huì)導(dǎo)致攻角增大。升力以及沿軸向的推力也較大，導(dǎo)致速度傾角增大，所以高度也會(huì)增大。

隨著高度的增加，空氣逐漸變稀薄，動(dòng)壓減小，導(dǎo)致發(fā)動(dòng)機(jī)推力減小，推力減小到一定程度后，加速度將小于零，速度將減小。所以，當(dāng)攻角較大時(shí)，速度會(huì)有一個(gè)先增加后減小的過程。

同理，速度傾角也有一個(gè)先增大后減小的過程。當(dāng)速度傾角減到小于零時(shí)，高度會(huì)減小。圖4給出了αm取值為 -2°、0°、2°時(shí)的攻角以及軌跡變化。

圖4 αm對(duì)飛行軌跡的影響Fig.4 Effect of αmon trajectory

由圖4可見，αm正、負(fù)不同時(shí)，攻角的變化規(guī)律相反。增大αm會(huì)增大攻角的峰值，相應(yīng)的，速度傾角和飛行高度均會(huì)增大，而飛行速度會(huì)減小。

從圖5可看到，增大a會(huì)縮短攻角由極值變化到α0的時(shí)間，a越大，α處于最小值的時(shí)間越長(zhǎng)。當(dāng)a＞1時(shí)，基本可認(rèn)為攻角一直取α0。所以，基本可認(rèn)為a的取值范圍為［0，1］。

采用類似方法，也可分析風(fēng)門系數(shù)表達(dá)式中4個(gè)待定參數(shù)對(duì)飛行軌跡的影響，從而可通過手動(dòng)調(diào)整的方式逐步給出各個(gè)參數(shù)的初始值。

圖5 a對(duì)飛行軌跡的影響Fig.5 Effect of a on trajectory

4 仿真算例

AHFV模型參見文獻(xiàn)［12］，初始高度為15 km，初始馬赫數(shù)為4，航跡角為15°，要求終端高度為30 km，終端馬赫數(shù)為6，終端航跡角為0°。

在求取初值的過程中，可將相關(guān)約束項(xiàng)放寬為

表1和表2分別給出了由GA和MM方法得到的初始值，以及用這些初始采用SQP算法得到的優(yōu)化結(jié)果。從表1和表2可看到，2種方法所給出的初始不盡相同，而以這2個(gè)初始得到的最優(yōu)優(yōu)化結(jié)果也是不同的，這也從另一個(gè)側(cè)面說明了SQP算法對(duì)初始值的敏感性。

表1 “GA+SQP”算法優(yōu)化結(jié)果Table 1 Optimal results of the‘GA+SQP’method

表2 “MM+SQP”算法優(yōu)化結(jié)果Table 2 Optimal results of the‘MM+SQP’method

圖6和圖7給出了2種方法得到的軌跡以及相應(yīng)的控制量。從圖6和圖7可看到，2種方法得到的控制量是不相同的，軌跡變化過程也不相同，但2種方法得到的軌跡均滿足約束條件。從性能指標(biāo)上看，二者差別也不大，采用GA+SQP算法結(jié)果略優(yōu)于MM+SQP算法的結(jié)果。

圖6 2種方法得到的優(yōu)化軌跡Fig.6 Optimal trajectories of the two methods

圖7 2種方法得到的控制量Fig.7 Optimal controls of the two methods

5 結(jié)論

(1)對(duì)吸氣式高超聲速飛行器上升段軌跡進(jìn)行了優(yōu)化設(shè)計(jì)。針對(duì)攻角和風(fēng)門系數(shù)，分別建立了基于特征參數(shù)的取值模型，從而得到以這些特征參數(shù)為優(yōu)化變量的參數(shù)優(yōu)化模型。

(2)著重分析了優(yōu)化初值的選取方法，采用GA+SQP算法和MM+SQP分別進(jìn)行了仿真計(jì)算，得到滿足條件的優(yōu)化結(jié)果，說明了方法的有效性。

［1］Mark A Van Buren，Kenneth D Mease.Minimum-fuel ascent to orbit using air-breathing proplusion［R］.NASA Report，1992:54-58.

［2］Corban J E，Calise A J，F(xiàn)landro G A.Optimal guidance and propulsion control for transatmospheric vehicle［C］//AIAA Guidance，Navagation and Contrl Conference.August 14-16，1989:123-130.

［3］Lu Ping.An inverse dynamic approach for trajectory optimization and guidance for an aerospace plane［R］.AIAA 92-4331-PC:217-223.

［4］Shahriar Keshmiri，Richard Colgren.Trajectory optimization for a generic hypersonic vehicle［C］//AIAA/AHI Space Planes and Hypersonic Systems and Technologies Conference，2006:1-8.

［5］宗群，田苓，竇立謙.基于 Gauss偽譜法的臨近空間飛行器上升段軌跡優(yōu)化［J］.宇航學(xué)報(bào)，2010(7):1776-1781.

［6］佘智勇，馬廣富，沈作軍.基于間接伴隨法的大氣層內(nèi)高超聲速飛行器最優(yōu)上升軌跡研究［J］.宇航學(xué)報(bào)，2010(8):1951-1957.

［7］Jung-Woo Park，Min-Jea Tahk，Hong-Gye Sung.Trajectory optimization for a supersonic air-breathing missile system using pseudo-spectral method［J］.Int.J.of Aeronautical ＆Space Science，2009，10:135-138.

［8］阮春榮.大氣中飛行的最優(yōu)軌跡［M］.北京:宇航出版社，1987:30-33.

［9］賈沛然，陳克俊，何力.遠(yuǎn)程火箭彈道學(xué)［M］.長(zhǎng)沙:國(guó)防科技大學(xué)出版社，1993:217-230.

［10］陳國(guó)良，王熙法.遺傳算法及其應(yīng)用［M］.北京:北京人民郵電出版社，1996.

［11］胡正東，沈薈萍，張士峰，等.SGKW軌道在線生成技術(shù)研究［J］.固體火箭技術(shù)，2009(7):191-195.

［12］Chuang C H，Hitoshi Morimoto.Optima periodic cruise for a hypersonic vehicle with constraints［R］.AIAA 96-3876:321-325.

A characteristic parameter trajectory optimization method for the ascent phase of the air-breathing hypersonic flight vehicle

DING Hong-bo，CAI Hong
(College of Aerospace and Material Engineering ，National University of Defence Technology，Changsha 410073，China)

A Characteristic Parameter Trajectory Optimization Method(CPTOM)for the ascent phase of the air-breathing vehicle was proposed and studied.The expression for aerodynamic force and thrust were specially studied，and the motion equations were founded.According to the characteristics of the ascent motion of air-breathing hypersonic flight vehicle(AHFV)，the control variable models based on exponential and polynomial functions were established.The models were determined by several characteristic parameters to be optimized.So the optimal control problem was translated to a nonlinear programming problem which can be calculated by the Sequence Quadric Programming(SQP).Because of the sensitivity of the initial value for optimization，the Genetic Arithmetic(GA)and Manual Method(MM)were proposed for initial value guess.

air-breathing hypersonic flight vehicle;ascent phase;characteristic parameters;trajectory optimization;initial value guess

V474.28

A

1006-2793(2012)03-0296-05

2011-10-14;

2011-11-21。

丁洪波(1982—)，男，博士生，研究方向?yàn)檐壽E優(yōu)化與制導(dǎo)。E-mail:whhplz2000@sohu.com

book=35,ebook=118

(編輯:崔賢彬)