豐志偉，江增榮，張青斌，葛健全，黃 浩

(1. 國(guó)防科技大學(xué) 空天科學(xué)學(xué)院， 湖南 長(zhǎng)沙 410073； 2.中國(guó)人民解放軍96901部隊(duì)， 北京 100094)

高超聲速滑翔飛行器(Hypersonic Glider Vehicle，HGV)具有高速度、高機(jī)動(dòng)、長(zhǎng)航程等諸多優(yōu)勢(shì)，成為未來(lái)飛行器發(fā)展的一個(gè)重要方向，具有廣闊的軍事應(yīng)用前景[1]。

彈道優(yōu)化技術(shù)是高超聲速滑翔飛行器的關(guān)鍵技術(shù)之一。隨著研究的不斷深入，在進(jìn)行高超聲速滑翔飛行器彈道設(shè)計(jì)之中，不僅要考慮飛行器在飛行過(guò)程中受到的氣動(dòng)熱、動(dòng)壓、過(guò)載、控制量、終端約束等復(fù)雜約束條件的限制，還需要考慮面向具體作戰(zhàn)任務(wù)的要求，如航路點(diǎn)、禁飛區(qū)等。面對(duì)諸多性能指標(biāo)，設(shè)計(jì)者往往不能追求單個(gè)指標(biāo)最優(yōu)，而是需要同時(shí)考慮多個(gè)相關(guān)的性能指標(biāo)，利用多目標(biāo)優(yōu)化方法進(jìn)行求解，獲得Pareto最優(yōu)解，從而對(duì)各個(gè)指標(biāo)進(jìn)行權(quán)衡和折中[2]。

對(duì)于高超聲速滑翔飛行器的多目標(biāo)彈道優(yōu)化問(wèn)題而言，工程上常用的方法是利用偏好(權(quán)重)信息，將多個(gè)目標(biāo)聚合為單個(gè)目標(biāo)，并采用成熟的單目標(biāo)優(yōu)化算法進(jìn)行求解。主要優(yōu)化方法包括加權(quán)法、理想點(diǎn)法、物理規(guī)劃法、主要目標(biāo)法等[3]。雍恩米等[4]應(yīng)用物理規(guī)劃法設(shè)計(jì)再入飛行最優(yōu)攻角；謝愈等[2]將物理規(guī)劃法應(yīng)用于高超聲速飛行器多目標(biāo)復(fù)雜約束滑翔彈道優(yōu)化設(shè)計(jì)，結(jié)合高斯偽譜法獲得了滿足設(shè)計(jì)者偏好的滿意解。上述基于偏好的方法每次只能優(yōu)化得到一個(gè)解，無(wú)法得到整個(gè)Pareto前沿(Pareto Front, PF)。為了獲得多個(gè)Pareto最優(yōu)解(Pareto Set, PS)，通常需要采用多目標(biāo)優(yōu)化方法來(lái)進(jìn)行設(shè)計(jì)。主要優(yōu)化方法包括各種多目標(biāo)進(jìn)化算法、多目標(biāo)粒子群算法等[5]。例如，陳剛、王明光等[6-7]將多目標(biāo)遺傳算法NSGA-II成功應(yīng)用于可重復(fù)使用運(yùn)載火箭最優(yōu)再入彈道優(yōu)化設(shè)計(jì)；豐志偉等[8]應(yīng)用分解多目標(biāo)進(jìn)化算法(Multiobjective Optimization Evolutionary Algorithm based on Decomposition, MOEA/D)求解火星再入飛行器一體化多目標(biāo)優(yōu)化設(shè)計(jì)問(wèn)題。

滑翔飛行器彈道優(yōu)化問(wèn)題是一個(gè)帶有復(fù)雜約束的最優(yōu)控制問(wèn)題，本文綜合考慮計(jì)算效率和精度，提出一種結(jié)合進(jìn)化算法和配點(diǎn)法的混合、分層多目標(biāo)優(yōu)化策略。為了獲得具有較高精度的全局最優(yōu)解，將多目標(biāo)彈道優(yōu)化問(wèn)題的求解分為兩步：①使用較少節(jié)點(diǎn)對(duì)彈道優(yōu)化問(wèn)題進(jìn)行離散化處理，采用MOEA/D進(jìn)行多目標(biāo)優(yōu)化計(jì)算，獲得Pareto前沿和全局最優(yōu)解集；②以多目標(biāo)最優(yōu)解作為基于二代小波節(jié)點(diǎn)自適應(yīng)配點(diǎn)法[9]的初始解進(jìn)行精細(xì)優(yōu)化，獲得較高精度的彈道。特別是，在MOEA/D和自適應(yīng)配點(diǎn)法中采用了橢球聚合法[10]將多個(gè)目標(biāo)函數(shù)聚合為單個(gè)目標(biāo)函數(shù)。

1 滑翔飛行器多目標(biāo)彈道優(yōu)化模型

1.1 再入動(dòng)力學(xué)方程

以地心距r、經(jīng)度λ和緯度φ描述飛行器位置，以速度大小V、速度傾角θ和速度偏角σ描述飛行器速度，建立的滑翔飛行器三自由度運(yùn)動(dòng)方程[11]可表示為

(1)

其中：X為阻力，Y為升力，L為側(cè)向力，m為飛行器質(zhì)量，g為引力加速度；速度傾角θ為速度矢量與當(dāng)?shù)厮矫娴膴A角，速度偏角σ為速度矢量在當(dāng)?shù)厮矫嫱队芭c正北方向的夾角；ωe為地球自轉(zhuǎn)角速度，ν為側(cè)傾角。由于滑翔飛行器為面對(duì)稱升力體，通常采用傾斜轉(zhuǎn)彎(Bank To Turn, BTT)技術(shù)，在上述運(yùn)動(dòng)方程中假設(shè)側(cè)滑角為零。

1.2 約束條件分析

滑翔飛行器多目標(biāo)彈道設(shè)計(jì)約束條件[12]包括：熱流、動(dòng)壓、過(guò)載、平衡滑翔、禁飛區(qū)、航路點(diǎn)等過(guò)程約束；彈道末端需滿足的終端狀態(tài)約束；對(duì)攻角、側(cè)傾角的控制量約束。

1)氣動(dòng)熱、過(guò)載和動(dòng)壓約束：為了保證飛行器結(jié)構(gòu)安全，對(duì)沿彈道飛行的熱流密度、動(dòng)壓、過(guò)載等進(jìn)行限制。

駐點(diǎn)熱流密度約束為

(2)

過(guò)載約束為

(3)

式中，α為攻角。

動(dòng)壓約束為

(4)

2)擬平衡滑翔約束：為保證彈道平穩(wěn)，沿彈道飛行器可獲得的最大升力須能夠平衡其他力，即

(5)

3)禁飛區(qū)約束：禁飛區(qū)指飛行需要規(guī)避的區(qū)域，包括躲避雷達(dá)探測(cè)等。為便于計(jì)算，假設(shè)禁飛區(qū)為無(wú)限高圓柱，彈道點(diǎn)與禁飛區(qū)中心在地面投影的最短距離滿足如下約束

|R0·arccos[sinφsinφnfz+

cosφcosφnfzcos(λ-λnfz)]|>Rnfz

(6)

其中，λnfz和φnfz分別表示禁飛區(qū)中心經(jīng)緯度，Rnfz表示禁飛區(qū)半徑。

4)控制量約束：由于受執(zhí)行機(jī)構(gòu)能力限制，攻角、側(cè)傾角等物理量的幅值存在約束邊界

(7)

5)終端約束：本文主要開(kāi)展最大射程優(yōu)化，不約束滑翔段終點(diǎn)的位置和速度。

1.3 性能指標(biāo)

按照飛行任務(wù)的不同，可選擇不同的性能指標(biāo)。常用的優(yōu)化目標(biāo)有終端射程最大、飛行時(shí)間最短、總加熱量最小等。如取終端射程最大為優(yōu)化指標(biāo)，則

J=min[-L(tf)]

(8)

當(dāng)選擇總加熱量最小作為性能指標(biāo)時(shí)，則

(9)

2 多目標(biāo)優(yōu)化問(wèn)題的混合求解策略

多目標(biāo)彈道優(yōu)化問(wèn)題的核心是獲得一組較為精確的Pareto全局最優(yōu)解集。為此，本文提出結(jié)合分解多目標(biāo)進(jìn)化算法和配點(diǎn)法的混合求解策略。為了保證兩種多目標(biāo)優(yōu)化方法求解問(wèn)題的一致性，均采用橢球聚合法[10]將多目標(biāo)優(yōu)化問(wèn)題轉(zhuǎn)化為單目標(biāo)優(yōu)化問(wèn)題。

2.1 混合求解策略

混合求解策略優(yōu)化流程如圖1所示。具體計(jì)算過(guò)程如下：

圖1 混合求解策略流程圖Fig.1 Flow chart of hybrid optimization strategy

1)定義多目標(biāo)彈道優(yōu)化問(wèn)題：定義狀態(tài)方程，確定控制量、優(yōu)化目標(biāo)和約束條件。

2)彈道優(yōu)化問(wèn)題的離散：分別采用打靶法(僅離散控制量)和配點(diǎn)法(同時(shí)離散控制量和狀態(tài)量)將彈道優(yōu)化問(wèn)題轉(zhuǎn)化為多目標(biāo)優(yōu)化問(wèn)題。

3)全局多目標(biāo)最優(yōu)解集的構(gòu)造：采用MOEA/D對(duì)控制量離散方法得到的多目標(biāo)優(yōu)化問(wèn)題進(jìn)行優(yōu)化計(jì)算。設(shè)置合理的優(yōu)化參數(shù)，如種群規(guī)模和迭代步數(shù)，獲得全局多目標(biāo)最優(yōu)解。

4)最優(yōu)解集的精細(xì)優(yōu)化：采用非線性規(guī)劃算法對(duì)配點(diǎn)法得到的優(yōu)化問(wèn)題進(jìn)行求解，以MOEA/D得到的控制量和狀態(tài)量作為初始解進(jìn)行優(yōu)化迭代，采用自適應(yīng)配點(diǎn)法進(jìn)行網(wǎng)格自適應(yīng)加密以提高求解精度。

需要說(shuō)明的是：①為了求得全局最優(yōu)解，首先采用進(jìn)化算法求解，考慮計(jì)算效率問(wèn)題，采用打靶法僅離散控制量，以降低設(shè)計(jì)變量的個(gè)數(shù)；②為了獲得較精細(xì)的彈道，利用配點(diǎn)法和橢球聚合法將多目標(biāo)彈道優(yōu)化問(wèn)題轉(zhuǎn)化為單目標(biāo)的非線性規(guī)劃問(wèn)題，利用非線性規(guī)劃算法進(jìn)行求解，并進(jìn)行節(jié)點(diǎn)自適應(yīng)加密。

2.2 橢球聚合法

在上述混合求解策略中，利用非線性規(guī)劃求解多目標(biāo)問(wèn)題時(shí)，需要用聚合公式將多目標(biāo)優(yōu)化問(wèn)題轉(zhuǎn)化為一組單目標(biāo)優(yōu)化問(wèn)題。為了保持目標(biāo)函數(shù)的連續(xù)性，本文采用橢球聚合公式[10]。該方法使用具有特定方向和較大偏心率的橢球(橢圓)與可達(dá)目標(biāo)區(qū)域Ω的最優(yōu)邊界相切，并尋找邊界上的切點(diǎn)，從而在一定條件下給出PF的近似。

對(duì)于極小化問(wèn)題而言，在大多數(shù)情況下，PF是可達(dá)目標(biāo)區(qū)域左下邊界?Ω的一部分。橢球聚合公式使用N個(gè)橢球(橢圓)與上述邊界?Ω相切。通過(guò)位于邊界?Ω上的切點(diǎn)給出PF的近似。該方法要求橢球具有一個(gè)較大的長(zhǎng)軸和若干個(gè)較小的短軸。橢球的中心一般位于可達(dá)目標(biāo)區(qū)域的左下側(cè)，通過(guò)控制橢球長(zhǎng)軸的方向來(lái)獲得不同的邊界點(diǎn)，如圖2所示。在一般情況下，當(dāng)長(zhǎng)軸方向均勻分布時(shí)，橢球聚合法可以產(chǎn)生近似均勻分布的PF。具體的做法是：將長(zhǎng)軸方向與一組均勻分布的權(quán)重矢量相關(guān)聯(lián)。

圖2 旋轉(zhuǎn)橢球方法Fig.2 Rotated ellipsoid method

根據(jù)上述思想，將多目標(biāo)優(yōu)化問(wèn)題處理為N個(gè)單目標(biāo)優(yōu)化問(wèn)題。令λ表示權(quán)重矢量，R(λ)表示將第一個(gè)目標(biāo)函數(shù)F1的坐標(biāo)軸旋轉(zhuǎn)到權(quán)重矢量λ方向的坐標(biāo)變換矩陣。橢球與邊界?Ω的切點(diǎn)是式(10)所示單目標(biāo)優(yōu)化問(wèn)題的全局最優(yōu)解。

s.t.x∈Ω

(10)

其中，

(11)

2.3 分解多目標(biāo)進(jìn)化算法

分解多目標(biāo)進(jìn)化算法MOEA/D由Zhang和Li于2007年提出[13]，該算法將多目標(biāo)優(yōu)化問(wèn)題轉(zhuǎn)化為多個(gè)單目標(biāo)優(yōu)化子問(wèn)題，對(duì)子問(wèn)題采用進(jìn)化算法以合作的方式同時(shí)進(jìn)行優(yōu)化。因此，其具有進(jìn)化機(jī)制簡(jiǎn)單、收斂速度快、PF近似均勻等優(yōu)點(diǎn)。

以橢球聚合公式為例，將原始多目標(biāo)優(yōu)化問(wèn)題分解成N個(gè)標(biāo)量?jī)?yōu)化問(wèn)題。令λ1, …,λN為均勻分布的一組權(quán)重矢量，z*為參考點(diǎn)，則第j個(gè)子問(wèn)題的目標(biāo)函數(shù)為

(12)

因?yàn)間RE關(guān)于λ連續(xù)，所以當(dāng)權(quán)重矢量λi和λj相互接近時(shí)，子問(wèn)題gRE(x|a,R(λi),z*)的優(yōu)化解應(yīng)該接近gRE(x|a,R(λj),z*)的優(yōu)化解。因此與權(quán)重矢量λi相鄰的關(guān)于gRE的信息能夠?qū)ψ訂?wèn)題gRE(x|a,R(λi),z*)的優(yōu)化起到輔助作用。MOEA/D包括初始化、更新和停止三個(gè)過(guò)程，具體步驟見(jiàn)文獻(xiàn)[13]。

在MOEA/D中，通常采用歐式距離測(cè)量?jī)蓚€(gè)權(quán)重矢量之間的鄰近性，并用于確定矢量λi的鄰居。獲得較為精確的參考點(diǎn)z*通常需要求解多個(gè)單目標(biāo)優(yōu)化問(wèn)題，比較耗時(shí)，因此該算法首先進(jìn)行初始化，然后在迭代過(guò)程中進(jìn)行更新。這樣可以使得理想點(diǎn)的更新隨著迭代優(yōu)化而進(jìn)行，節(jié)省了計(jì)算時(shí)間。

2.4 自適應(yīng)配點(diǎn)法

針對(duì)軌跡優(yōu)化問(wèn)題中精度和效率之間的矛盾，豐志偉等曾提出基于二代小波的節(jié)點(diǎn)自適應(yīng)加密方法[9]。該方法采用Runge-Kutta方法將原軌跡優(yōu)化問(wèn)題轉(zhuǎn)化為非線性規(guī)劃問(wèn)題，并采用非線性規(guī)劃算法求解；根據(jù)控制或狀態(tài)函數(shù)的小波系數(shù)幅值確定自適應(yīng)加密的節(jié)點(diǎn)，并進(jìn)行序列優(yōu)化。

基于二代小波的節(jié)點(diǎn)自適應(yīng)算法輸入?yún)?shù)包括：初始分辨率水平J0(控制初始節(jié)點(diǎn)個(gè)數(shù))；最大分辨率水平Jmax(控制背景節(jié)點(diǎn)個(gè)數(shù)，即最小節(jié)點(diǎn)間隔)；小波系數(shù)幅值的閾值ε。終止條件：前后兩次節(jié)點(diǎn)位置相同或達(dá)到預(yù)定的迭代次數(shù)。計(jì)算流程如下：

1)根據(jù)初始分辨率J0生成均勻節(jié)點(diǎn)GJ0(或{τi})以及該均勻節(jié)點(diǎn)上的狀態(tài)函數(shù)x和控制函數(shù)u的初始猜測(cè)；

2)利用非線性規(guī)劃算法優(yōu)化由原軌跡優(yōu)化問(wèn)題轉(zhuǎn)化得到的非線性規(guī)劃問(wèn)題，得到狀態(tài)函數(shù)x和最優(yōu)控制函數(shù)u；

3)根據(jù)節(jié)點(diǎn)自適應(yīng)算法確定下一步優(yōu)化的自適應(yīng)節(jié)點(diǎn)，若存在多個(gè)控制函數(shù)，則將這些控制函數(shù)分別得到的自適應(yīng)節(jié)點(diǎn)進(jìn)行合并；

4)利用前一次優(yōu)化計(jì)算的解x和u，利用小波插值計(jì)算新的自適應(yīng)節(jié)點(diǎn)上的初始猜測(cè)，并進(jìn)行循環(huán)迭代，直到滿足終止條件。

需要說(shuō)明的是，小波系數(shù)幅值的閾值通?？扇棣?γ(umax-umin)。對(duì)于間斷函數(shù)，γ可取0.005～0.01；對(duì)于連續(xù)函數(shù)，γ可取0.001～0.005。

3 仿真分析

3.1 仿真條件和參數(shù)

以通用航空飛行器[14](Common Aero Vehicle，CAV)氣動(dòng)性能為例進(jìn)行多目標(biāo)彈道優(yōu)化，通過(guò)擬合可得到CAV-H模型的氣動(dòng)系數(shù)

(13)

優(yōu)化指標(biāo)包括射程最大和總加熱量最小?？紤]的約束條件包括駐點(diǎn)熱流密度、動(dòng)壓、過(guò)載以及禁飛區(qū)約束，仿真參數(shù)設(shè)置如表1所示。

利用本文提出的混合優(yōu)化算法，將攻角和側(cè)傾角離散為15個(gè)點(diǎn)，加上終端時(shí)間共31個(gè)設(shè)計(jì)變量。

在MOEA/D優(yōu)化過(guò)程中采用旋轉(zhuǎn)橢球聚合公式，種群規(guī)模(子問(wèn)題個(gè)數(shù))為300，迭代步數(shù)為1000，鄰居規(guī)模為20。采用300個(gè)均勻分布的權(quán)重矢量，每個(gè)權(quán)重矢量的分量為正值，和為1。節(jié)點(diǎn)自適應(yīng)算法中，為了節(jié)省時(shí)間，將優(yōu)化子問(wèn)題數(shù)目取為21，小波系數(shù)幅值的閾值取εi=0.001·(uimax-uimin)。針對(duì)每個(gè)子問(wèn)題，取MOEA/D優(yōu)化結(jié)果中相應(yīng)權(quán)重系數(shù)對(duì)應(yīng)的子問(wèn)題的解(控制量和飛行時(shí)間)進(jìn)行積分得到狀態(tài)變量，并作為配點(diǎn)法中的初始解。采用序列二次規(guī)劃算法SNOPT[15]求解非線性規(guī)劃問(wèn)題。

表1 彈道優(yōu)化約束條件

3.2 彈道約束條件處理

由于MOEA/D無(wú)法直接處理多目標(biāo)彈道優(yōu)化問(wèn)題中存在的復(fù)雜非線性約束，需要對(duì)約束進(jìn)行處理。對(duì)于不等式約束，在聚合目標(biāo)gte(x|λ)中引入懲罰項(xiàng)，MOEA/D算法更新新解時(shí)不僅考慮了解的聚合目標(biāo)，而且還融合了約束違約程度。采用罰函數(shù)的MOEA/D在求解控制-結(jié)構(gòu)一體化設(shè)計(jì)中得到了成功應(yīng)用[16]。

假設(shè)第i個(gè)約束gi(x)≤0,i=1,…,nc的違約度[17]定義為

Wi(x)=max[0,gi(x)],i=1,…,nc

(14)

如果Wi(x)=0，則x可行，否則不可行。將目標(biāo)函數(shù)修改為

F′(x)=[f1(x)+aW,…，fm(x)+aW]T

=[f′1(x),…，f′m(x)]T

(15)

其中

W=max[W1,…,Wnc]

(16)

nc是約束的個(gè)數(shù)，參數(shù)a>0。

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

利用進(jìn)化算法得到的近似PF如圖3所示，從圖中可以看出，最大射程解在圖中以“▽”標(biāo)示，此時(shí)射程為13 799 km，總加熱量為18 173 MJ/m2；最小總加熱量解以“△”標(biāo)示，此時(shí)總加熱量為10 745 MJ/m2，射程為10 930 km。兩個(gè)目標(biāo)具有明顯的沖突性，即射程最大時(shí)，總加熱量也最大，反之亦然。為了進(jìn)一步對(duì)比分析，圖4給出了三條多目標(biāo)優(yōu)化典型彈道(兩條單目標(biāo)最優(yōu)彈道；一條典型折中解，圖3中以“□”標(biāo)示)對(duì)應(yīng)的控制量以及主要彈道參數(shù)。從圖4曲線可以看出，最大射程彈道的攻角為最大升阻比攻角；最小總加熱量彈道攻角較大，且高度在縱平面大幅跳躍。

圖3 最大射程-最小總加熱量?jī)赡繕?biāo)PFFig.3 PF for maximizing down range and minimizing total heating

(a) 高度(a) Height

(b) 速度傾角(b) Velocity dip angle

(c) 駐點(diǎn)熱流密度(c) Heating rate of stagnation point

(d) 攻角(d) Attack angle

(e) 側(cè)傾角(e) Bank angle圖4 多目標(biāo)優(yōu)化彈道典型解(MOEA/D)Fig.4 Typical solutions of multiobjective trajectory optimization (MOEA/D)

采用自適應(yīng)配點(diǎn)法優(yōu)化得到的近似PF如圖5所示，與進(jìn)化算法MOEA/D優(yōu)化得到的近似PF相比，近似PF得到較好的改善。在最大射程解附近，MOEA/D較接近配點(diǎn)法獲得的解，而在最小總加熱量解附近其偏差較大。由此說(shuō)明采用15個(gè)點(diǎn)近似攻角和側(cè)傾角來(lái)優(yōu)化最大射程是足夠的，但是為了優(yōu)化最小總加熱量彈道，需要增加設(shè)計(jì)變量的個(gè)數(shù)。

圖5 兩種方法的PF對(duì)比Fig.5 PF comparison of two methods

自適應(yīng)配點(diǎn)法優(yōu)化得到的攻角和側(cè)傾角曲線如圖6所示。整體規(guī)律與進(jìn)化算法得到的解一致，但是最小總加熱量解攻角曲線更復(fù)雜。

(a) 攻角(a) Attack angle

(b) 側(cè)傾角(b) Bank angle圖6 多目標(biāo)優(yōu)化彈道典型解(自適應(yīng)配點(diǎn)法)Fig.6 Typical solutions of multiobjective trajectory optimization (adaptive collocation method)

4 結(jié)論

針對(duì)滑翔飛行器多目標(biāo)彈道優(yōu)化設(shè)計(jì)問(wèn)題，綜合考慮計(jì)算效率和精度，結(jié)合分解多目標(biāo)進(jìn)化算法和自適應(yīng)配點(diǎn)法提出一種分層、混合優(yōu)化方法。以美國(guó)CAV-H為例進(jìn)行仿真分析，計(jì)算結(jié)果驗(yàn)證了所提混合優(yōu)化方法的有效性，不論是MOEA/D還是配點(diǎn)法，均能夠得到較均勻分布的Pareto前沿。