樊朋飛,郭云鶴2,凡永華,閆 杰

(1.西北工業(yè)大學(xué) 航天學(xué)院，西安 710072; 2.上海機(jī)電工程研究所，上海 201109)

0 引言

高超聲速滑翔飛行器具有高升阻比氣動外形，能夠在大氣層內(nèi)實現(xiàn)高速、遠(yuǎn)距離滑翔和大范圍橫向機(jī)動。因其落點覆蓋范圍大、機(jī)動能力強，難以進(jìn)行軌跡預(yù)測和攔截，因而具有重大的軍事應(yīng)用前景。平衡滑翔式軌跡是指滑翔飛行器飛行中升力在縱向的投影與離心力、重力近似滿足平衡關(guān)系的一類軌跡。該類軌跡具有良好的熱流率曲線以及穩(wěn)定的操控特性和機(jī)動能力，在軌跡優(yōu)化算法的研究和工程應(yīng)用中受到廣泛關(guān)注[1]。

再入可達(dá)區(qū)域體現(xiàn)了飛行器的打擊覆蓋能力和機(jī)動性，是衡量其作戰(zhàn)性能的重要指標(biāo)，對該指標(biāo)進(jìn)行精確高效的計算對于打擊任務(wù)規(guī)劃、目標(biāo)和航路點選取具有實用參考價值。文獻(xiàn)[2]總結(jié)了計算可達(dá)區(qū)域的四種常用方法：軌跡優(yōu)化法、剖面規(guī)劃法、常值傾側(cè)角法和橢圓近似法。軌跡優(yōu)化法中，常采用偽譜法[3-4]、粒子群算法[5]、遺傳算法[6]等優(yōu)化算法求解可達(dá)域的邊界軌跡。該類算法由于最優(yōu)性能指標(biāo)的保證，具有較高計算精度，但其缺點在于求解再入軌跡這種具有復(fù)雜狀態(tài)約束和動力學(xué)模型的問題時，效率不高。同時，由于高速、高升阻比的特點，再入滑翔軌跡具有天然的“振蕩特性”，平衡滑翔軌跡在數(shù)值軌跡優(yōu)化算法中難以得到。文獻(xiàn)[7]和[8]通過在動力學(xué)方程中加入振蕩抑制反饋項以提升軌跡的阻尼特性，并采用偽譜法對改進(jìn)后的動力學(xué)模型進(jìn)行軌跡優(yōu)化。文獻(xiàn)[9]通過引入無損松弛技術(shù)，得到了具有仿射結(jié)構(gòu)的動力學(xué)模型，以減弱控制量對動力學(xué)方程的耦合作用。在此基礎(chǔ)上，文獻(xiàn)[10]設(shè)計了光滑的標(biāo)稱航跡傾角剖面，并將其與實際傾角的偏差引入最優(yōu)指標(biāo)中，保證求解軌跡的光滑特性。

本文研究采用連續(xù)凸優(yōu)化算法求解平衡滑翔式軌跡的再入可達(dá)區(qū)域問題。為提高優(yōu)化求解的效率，首先引入“準(zhǔn)平衡滑翔假設(shè)”條件對動力學(xué)方程進(jìn)行降階，進(jìn)而以速度變量替代時間作為自變量，使得動力學(xué)方程的階次減小一倍。通過這種方法，不僅大大降低了最優(yōu)軌跡的求解難度，而且消除了高升阻比飛行器浮沉特性對于控制輸入敏感的問題，此時，軌跡的平滑度僅與控制量的平滑程度相關(guān)，而無需引入額外的阻尼項。隨后，對優(yōu)化軌跡的求解問題進(jìn)行了轉(zhuǎn)化，通過將降階的動力學(xué)方程進(jìn)行線性化、離散化處理，并引入二階錐約束條件對控制量的平滑度進(jìn)行約束，將非線性最優(yōu)控制問題轉(zhuǎn)化為可被高效求解的二階錐規(guī)劃(SOCP)問題，通過連續(xù)求解一系列凸優(yōu)化子問題獲得最優(yōu)的可達(dá)域邊界軌跡。以CAV-H高超聲速滑翔飛行器為模型的可達(dá)區(qū)域求解算例驗證了方法的有效性。

1 再入飛行器動力學(xué)建模

忽略地球自轉(zhuǎn)的影響，將地球視為靜止的球體，建立再入滑翔飛行器的質(zhì)心動力學(xué)方程如下：

(1)

(2)

(3)

(4)

(5)

(6)

式中,狀態(tài)量V，γ，ψ，r，θ和φ分別表示飛行器的速度、航跡傾角、航向角、地心距、經(jīng)度和緯度；m為飛行器質(zhì)量；g為重力加速度；σ為傾側(cè)角。D和L分別為阻力和升力，其計算公式為:

(7)

(8)

再入過程中，飛行器受到再入走廊的約束，即滿足：

(9)

(10)

(11)

式(9)～(11)分別表示熱流率約束、法向過載約束和動壓約束。α為攻角；CQ為熱流常數(shù)。

為使軌跡末端滿足一定的交接班條件，需對末端高度和速度進(jìn)行約束，其形式如下：

r(tf)-R0=Hf

V(tf)=Vf

(12)

式中,R0為地球半徑。

2 平衡滑翔軌跡可達(dá)區(qū)域計算

采用軌跡優(yōu)化方法計算飛行器的可達(dá)區(qū)域時，通常將問題轉(zhuǎn)化一系列縱程固定，橫程指標(biāo)最大的最優(yōu)控制問題。首先，計算可達(dá)區(qū)域的遠(yuǎn)界和近界，得到縱程的取值范圍；然后，在縱程區(qū)間內(nèi)依次取值作為末端狀態(tài)約束，并求解最優(yōu)控制使末端橫程指標(biāo)最大。不失一般性，本文取0°經(jīng)度、緯度和90°航向角作為初始位置航向條件，則飛行器的經(jīng)度、緯度值可直接作為縱程和橫程使用，且可達(dá)區(qū)域具有關(guān)于縱軸對稱的特性。該條件計算得到的可達(dá)區(qū)域可經(jīng)球面坐標(biāo)變換[11]后投影至其他初始位置航向條件的情況。

2.1 運動方程簡化

(13)

根據(jù)升阻比關(guān)系，計算阻力：

(14)

則簡化后的動力學(xué)方程為：

(15～18)

考慮到可達(dá)區(qū)域的計算不關(guān)心飛行時間，而速度V的初始和末端狀態(tài)固定，將上述方程組對V求導(dǎo)，消去速度方程和時間變量，得到以V為自變量的運動方程：

(19～21)

式中，下標(biāo)V表示狀態(tài)量相對速度的導(dǎo)數(shù)。此時，方程中狀態(tài)量僅剩ψ、θ和φ；升阻比CL/CD為馬赫數(shù)Ma和攻角α的函數(shù)，由于再入飛行器的攻角通常采用固定剖面，若忽略高空聲速的變化，則CL/CD僅與飛行速度相關(guān)。在方程中的狀態(tài)量和控制量σ確定后，方程中未包含的狀態(tài)信息H和γ也可計算得出。其中H可通過聯(lián)立式(7)和式(13)求解：

(22)

將上式對V求導(dǎo)：

(23)

根據(jù)式(1)與式(4)，相除后又可得：

(24)

聯(lián)立上兩式得：

(25)

式中,Dm=D/m為阻力方向的加速度。

2.2 約束條件轉(zhuǎn)化

再入軌跡的優(yōu)化中，再入走廊通常以V-H或者V-Dm走廊的形式出現(xiàn)。而本文中，H或Dm并不顯含于運動方程(19～21)中，基于 “準(zhǔn)平衡滑翔假設(shè)”，上述形式的約束可通過式(13)或式(14)轉(zhuǎn)化為V-σ形式的約束，即滿足：

|σ|≤σr(V)

(26)

式中,σr表示平衡滑翔狀態(tài)下不超出再入走廊約束的傾側(cè)角最大值。同理，初始和末端高度約束可轉(zhuǎn)化為初始、末端傾側(cè)角約束：

σ0=σ(H0,V0)

σf=σ(Hf,Vf)

(27)

2.3 優(yōu)化問題求解

本節(jié)將最優(yōu)軌跡求解問題由非線性規(guī)劃問題轉(zhuǎn)化為連續(xù)求解的一系列SOCP[12]問題。SOCP是凸優(yōu)化問題的一種特例，其要求性能指標(biāo)函數(shù)為優(yōu)化變量的線性組合，受到線性約束和二階錐約束。由于該類問題可運用對偶內(nèi)點法[13](primal-dual interior-point algorithm)高效求解，因此，將文中所述優(yōu)化問題轉(zhuǎn)化為可被SOCP方法求解的標(biāo)準(zhǔn)形式成為問題的關(guān)鍵。將非線性運動方程寫為如下一般形式：

(28)

首先，對非線性動力學(xué)模型進(jìn)行線性化、離散化處理，使其滿足線性等式約束形式。令{xk(V);σk(V)}表示連續(xù)求解過程中的第k次解，將動力學(xué)方程在其附近一階展開可得：

(29)

其中：

c(xk,σk,V)=f(xk,σk,V)-A(xk,σk,V)xk-

B(xk,σk,V)σk

(30)

隨后，將速度區(qū)間[V0,Vf]等分為N+1個子區(qū)間，則其步長ΔV=(Vf-V0)/N，離散點可表示為{V0,V1,V2,…,VN}，其中，Vi=V0+iΔV。則狀態(tài)量x和控制量σ被離散化為xi=x(Vi)，σi=σ(Vi)。狀態(tài)方程可進(jìn)行數(shù)值積分轉(zhuǎn)化為如下形式：

(31)

(32)

Mz=F

(33)

其中：

(34)

(35)

經(jīng)過線性化和離散化后的最優(yōu)問題可描述為：在離散化的運動方程(33)約束下，求取最優(yōu)控制，使橫程指標(biāo)J=-φ(Vf)最小，并滿足，

1)初始狀態(tài)約束:

x(V0)=x0

(36)

2)末端經(jīng)度約束:

θ(Vf)=θf

(37)

3)再入走廊約束:

|σi|≤σr(Vi)，i=0,…,N

(38)

以及式(27)中的初始和末端傾側(cè)角約束：

σ(V0)=σ0

σ(Vf)=σf

式(37)中，θf∈[θfmin,θfmax]。

可以看出，此時最優(yōu)問題只包含線性的約束條件和性能指標(biāo)，符合SOCP問題的求解要求。為改善最優(yōu)問題的可解性，并且保證最優(yōu)軌跡的平滑特性，對上述最優(yōu)問題的指標(biāo)和約束項進(jìn)行改進(jìn)。其中最優(yōu)指標(biāo)變?yōu)椋?/p>

J=-φ(Vf)+kθθa+kσσvar

(39)

式中,θa、σvar為新增優(yōu)化變量，kθ和kσ為其權(quán)重系數(shù)。θa表示對末端經(jīng)度指標(biāo)的滿足程度，有如下不等式關(guān)系：

|θ(Vf)-θf|≤θa

(40)

σvar用于增加控制量的平滑度，其與離散控制量滿足不等式關(guān)系：

(41)

考慮到線性化處理后，最優(yōu)解的可信域問題，應(yīng)對每次求解結(jié)果的范圍進(jìn)行約束，即需要滿足：

|xk+1(Vi)-xk(Vi)|≤δ，i=0,…,N

(42)

式中,取δ=[10 5 5]T。

(43)

最優(yōu)軌跡求解時，首先給出一組初始猜測值{x0(V);σ0(V)}，并反復(fù)將求解結(jié)果作為下一次求解的輸入值，當(dāng)滿足：

max|xk+1(Vi)-xk(Vi)|≤ε，i=0,…,N

(44)

時，則認(rèn)為連續(xù)凸優(yōu)化過程已收斂至最優(yōu)解。式中：ε為一小量。

3 仿真分析

采用CAV-H的氣動、總體參數(shù)建立再入動力學(xué)模型。仿真算例的初始條件為：H0=56 km，V0=6 000 m/s。

計算程序在MATLAB環(huán)境中編寫實現(xiàn)，運用YALMIP軟件[14]進(jìn)行問題建模并調(diào)用MOSEK軟件包[13]中的SOCP問題求解器進(jìn)行求解。

圖1 再入可達(dá)區(qū)域

圖2 傾側(cè)角曲線

圖3 高度曲線

再入可達(dá)區(qū)域的計算結(jié)果如圖1～4所示。其中，圖1為可達(dá)域邊界軌跡的經(jīng)緯度信息，從中可以看出，可達(dá)域形狀呈近似橢圓形分布，在上述條件約束下，可達(dá)域縱向遠(yuǎn)邊界達(dá)到69.6°，側(cè)向邊界在經(jīng)度45°處取得最大，為25.9°。

圖2～4分別為軌跡的傾側(cè)角、高度和阻力加速度信息隨速度的變化曲線。從中可以看出，隨著縱程的減小，軌跡的傾側(cè)角逐漸增大，且軌跡的初段傾向于以最大傾側(cè)角變化率增加至某一較大值，隨后逐漸減小，至末端時迅速收斂至滿足末端傾側(cè)角約束。從高度曲線圖可以看出，最優(yōu)軌跡傾向于在初段以較低的高度飛行，以盡快調(diào)整軌跡偏角，而在中、末端以較高高度飛行以增大航程并滿足末端狀態(tài)約束，這與傾側(cè)角曲線圖的表現(xiàn)是一致的。

圖4 阻力加速度曲線

為研究可達(dá)域范圍隨末端速度的變化規(guī)律，分別計算了末速Vf=2 000 m/s和Vf=1 000 m/s時的可達(dá)區(qū)域，并與圖1中Vf=1 500 m/s時的結(jié)果進(jìn)行了對比，其結(jié)果如圖5所示。從中可以看出，可達(dá)區(qū)域的范圍隨末速減小而顯著增大。其最大縱、橫程比較結(jié)果如表1所示。可以看出，末速度變化相同數(shù)值時，可達(dá)域的最大縱、橫程值增加量隨末速度的減小而減小，這是由于隨著速度減小，相同速度增量帶來的航程收益逐漸變小。

圖5 可達(dá)區(qū)域隨末速度變化圖

末速/(m/s)縱程/(°)橫程/(°)100070.9726.82150069.5825.88200067.3024.19

為驗證本文方法的有效性和優(yōu)越性，基于末速Vf=1 500 m/s的算例，將采用本文方法與文獻(xiàn)[2]中敘述的工程常用的剖面規(guī)劃法進(jìn)行了對比，其結(jié)果如圖6所示。從中可以看出，兩種方法得到的可達(dá)區(qū)域形狀基本相符。在可達(dá)遠(yuǎn)邊界和最大橫程的求解方面，本文方法都優(yōu)于剖面規(guī)劃法。需要指出的是，由于本文采用了“準(zhǔn)平衡滑翔”假設(shè)簡化運動方程，在求取較近射程的最優(yōu)軌跡時，所得橫程邊界小于剖面規(guī)劃法，且可達(dá)區(qū)域的近界亦較遠(yuǎn)。這是因為對于較近射程的軌跡，其狀態(tài)變化比較劇烈，往往不滿足平衡假設(shè)條件，此時引入該約束則限制了軌跡變化的幅度，使得軌跡近界的計算不如采用完整動力學(xué)模型精確。但是由于遠(yuǎn)邊界的計算相對極限近界而言更受關(guān)注，因此本文中的假設(shè)條件仍是合理的。

圖6 兩種方法的可達(dá)區(qū)域?qū)Ρ?/p>

4 結(jié)束語

本文研究了再入滑翔飛行器平衡滑翔軌跡的可達(dá)區(qū)域計算問題?；凇皽?zhǔn)平衡滑翔”假設(shè)簡化了運動模型，并采用連續(xù)凸優(yōu)化方法對再入走廊約束下的可達(dá)區(qū)域進(jìn)行了求解。仿真結(jié)果表明：本文方法能夠降低最優(yōu)軌跡的求解難度，獲得光滑無振蕩的平衡滑翔軌跡；所計算得到的可達(dá)區(qū)域遠(yuǎn)界和橫向邊界都優(yōu)于常規(guī)工程方法，體現(xiàn)了該方法的有效性和優(yōu)越性。