陸浩然 邱 薇 孫海亮 鄭 偉

1.國防科技大學(xué)，長沙 410073 2.北京宇航系統(tǒng)工程研究所，北京 100076 3.火箭軍研究院，北京 100089

0 引言

高速飛行器具有速度快、升阻比大、航程遠(yuǎn)、機(jī)動(dòng)性強(qiáng)等特點(diǎn)，已成為各種空間力量發(fā)展的重點(diǎn)方向。與單個(gè)飛行器相比，高速飛行器編隊(duì)在協(xié)同感知、博弈對(duì)抗、多維運(yùn)用等方面具有更為廣闊的應(yīng)用前景，因此研究其編隊(duì)控制問題具有重要意義。

針對(duì)多飛行器編隊(duì)控制問題，國內(nèi)外學(xué)者開展了大量的研究，主要方法有PID控制、最優(yōu)控制、高斯偽譜、滑?？刂频?。韋常柱等[1]建立了導(dǎo)彈的相對(duì)運(yùn)動(dòng)模型，將領(lǐng)彈的運(yùn)動(dòng)狀態(tài)視為輸入擾動(dòng)，采用線性二次(LQR)最優(yōu)控制理論設(shè)計(jì)了導(dǎo)彈編隊(duì)飛行保持控制器。張磊等[2]采用仿射非線性系統(tǒng)最優(yōu)控制理論，設(shè)計(jì)了基于領(lǐng)彈-從彈法的多導(dǎo)彈編隊(duì)控制器，使之能夠在領(lǐng)彈機(jī)動(dòng)地情況下快速、穩(wěn)定地實(shí)現(xiàn)編隊(duì)隊(duì)形的形成和保持。Wen等[3]基于最小化隊(duì)形代價(jià)得到了一種隊(duì)形沖突預(yù)報(bào)和協(xié)調(diào)的算法，并通過仿真分析證明了其有效性。Wang等[4]考慮多種約束的編隊(duì)隊(duì)形最優(yōu)軌跡，利用跟蹤算法實(shí)現(xiàn)了在多種擾動(dòng)條件下的隊(duì)形形成與穩(wěn)定控制。李文等[5]面向多導(dǎo)彈協(xié)同飽和攻擊同一目標(biāo)需求，設(shè)計(jì)了綜合考慮脫靶量、攻擊時(shí)間、攻擊角度和視場角限制的三維協(xié)同制導(dǎo)律。但是，上述研究對(duì)象是彈道式或巡航式導(dǎo)彈，其氣動(dòng)力常被忽略或簡化為常值，與實(shí)際偏差較大，而對(duì)于高速飛行器而言，先后經(jīng)歷寬速域、大空域的復(fù)雜飛行環(huán)境，高動(dòng)態(tài)氣動(dòng)特性是必須重點(diǎn)考慮的因素。

近年來，偽譜法(Pseudospectral Method, PM)引起了學(xué)者們的廣泛關(guān)注，在航空航天飛行器軌跡優(yōu)化領(lǐng)域得到深入研究。王芳等[6]基于高斯偽譜法，提出編隊(duì)形成-保持-攻擊一體化的時(shí)間最優(yōu)控制算法，實(shí)現(xiàn)了導(dǎo)彈編隊(duì)協(xié)同攻擊全過程快速控制。姚寅偉等[7]提出了一種基于“初值軌跡生成+Gauss偽譜法+SQP求解NLP”的高超聲速飛行器多維再入軌跡優(yōu)化方法，既利用了Gauss偽譜法收斂快#精度高的特點(diǎn)，又結(jié)合初值軌跡生成算法，彌補(bǔ)了Gauss偽譜法對(duì)初值敏感的不足。雷虎民等[8]提出一種自適應(yīng)更新的偽譜法，以較小的網(wǎng)格規(guī)模獲得較高的精度，并成功應(yīng)用于最優(yōu)控制問題。

本文針對(duì)高速飛行器編隊(duì)控制問題，設(shè)計(jì)了一種基于偽譜法的編隊(duì)成形隊(duì)形自適應(yīng)優(yōu)化方法。該方法通過分析高速飛行器的運(yùn)動(dòng)狀態(tài)，利用最低能量原則和最小偏轉(zhuǎn)角度對(duì)編隊(duì)隊(duì)形成形的位置進(jìn)行設(shè)計(jì)，避免了初始狀態(tài)與期望隊(duì)形偏差過大時(shí)問題無解的情況；同時(shí)設(shè)計(jì)出期望的隊(duì)形成形狀態(tài)，將原編隊(duì)控制問題轉(zhuǎn)化為軌跡設(shè)計(jì)問題，采用Radau偽譜法對(duì)軌跡設(shè)計(jì)問題進(jìn)行求解；通過數(shù)值仿真驗(yàn)證了該編隊(duì)設(shè)計(jì)方法的有效性，結(jié)果表明高速飛行器能在一定范圍內(nèi)的初始狀態(tài)下形成穩(wěn)定的空間構(gòu)型。

1 動(dòng)力學(xué)方程及約束

在地心坐標(biāo)系下，高速飛行器運(yùn)動(dòng)方程為：

(1)

式中：r為地心距，是飛行器質(zhì)心到地心的距離；λ和φ分別為飛行器所在位置的地心經(jīng)度、地心緯度；V為飛行器相對(duì)地球的速度大??；γ為速度傾角，是速度矢量與當(dāng)?shù)厮矫娴膴A角；ψ為速度方位角，是速度矢量與正北方向的夾角，順時(shí)針為正；σ為傾側(cè)角，是速度系與航跡系繞速度軸旋轉(zhuǎn)的角度，以逆時(shí)針為正；D和L分別為阻力加速度、升力加速度。

實(shí)際上，阻力、升力滿足以下關(guān)系式：

(L,D)=f(α,β,v,h)

(2)

式中：α為攻角，β為側(cè)滑角，h為飛行高度，均與飛行器當(dāng)前運(yùn)動(dòng)狀態(tài)有關(guān)。需要說明的是，在仿真過程中，上述函數(shù)是通過對(duì)標(biāo)準(zhǔn)氣動(dòng)參數(shù)插值得到的。

高速飛行器在大氣層內(nèi)機(jī)動(dòng)飛行時(shí)，除了需要滿足高度、速度等終端約束外，還必須滿足熱流、動(dòng)壓、過載等過程約束條件。

1)熱流約束

高速飛行器在大氣層內(nèi)飛行時(shí)，與來流發(fā)生劇烈碰撞和摩擦，導(dǎo)致了嚴(yán)重的氣動(dòng)加熱問題。在理論研究中，一般選擇駐點(diǎn)熱流密度來表征氣動(dòng)加熱的嚴(yán)重程度，如式(3)所示：

(3)

為避免駐點(diǎn)處熱流密度超出約束，應(yīng)滿足：

(4)

2)過載約束

高速飛行器在大氣層內(nèi)飛行時(shí)，為避免出現(xiàn)結(jié)構(gòu)性破壞，需對(duì)氣動(dòng)過載加以約束，如式(5)所示：

(5)

式中：n為氣動(dòng)過載，nmax為飛行器能承受的最大氣動(dòng)過載，g0為地球引力加速度。

3)動(dòng)壓約束

動(dòng)壓是飛行器最重要的特征量之一，其極限值主要取決于熱防護(hù)材料強(qiáng)度和氣動(dòng)控制鉸鏈力矩。動(dòng)壓必須限制在一定范圍內(nèi)，以確保表面絕熱材料結(jié)構(gòu)不受破壞，不超過對(duì)應(yīng)于控制氣動(dòng)操作面所要求的最大鉸鏈力矩允許動(dòng)壓。

動(dòng)壓約束q的表達(dá)式為：

(6)

式中：qmax為約束的最大動(dòng)壓。

2 編隊(duì)隊(duì)形快速成形實(shí)現(xiàn)方法

高速飛行器在編隊(duì)飛行過程中，期望隊(duì)形根據(jù)任務(wù)要求設(shè)計(jì)，當(dāng)編隊(duì)開始時(shí)刻飛行器位置相對(duì)于期望隊(duì)形存在較大偏差時(shí)，如果對(duì)編隊(duì)成形隊(duì)形直接賦值，將會(huì)在初始狀態(tài)和氣動(dòng)力約束下出現(xiàn)無解現(xiàn)象。針對(duì)該問題，本節(jié)分析編隊(duì)控制過程中飛行器機(jī)械能的變化，并以此為原則設(shè)計(jì)編隊(duì)成形時(shí)的隊(duì)形位置和過渡時(shí)間，確保從初始狀態(tài)到終端期望隊(duì)形存在可行解。

2.1 問題轉(zhuǎn)化

高速飛行器在大氣層內(nèi)飛行時(shí)，主要受地球引力和氣動(dòng)力作用，升力控制速度方向，阻力改變速度大小，動(dòng)能轉(zhuǎn)化為熱能散掉，整個(gè)飛行過程中的機(jī)械能逐漸降低。因此，對(duì)于多飛行器編隊(duì)控制，將以最低能量飛行器作為參考，設(shè)計(jì)編隊(duì)成形狀態(tài)。

為確定編隊(duì)形成時(shí)的飛行器與初始時(shí)刻飛行器的對(duì)應(yīng)關(guān)系，設(shè)計(jì)一種以最小偏轉(zhuǎn)角度為原則的飛行器標(biāo)號(hào)識(shí)別方法。假定飛行器k的速度和速度傾角不變，經(jīng)過一段時(shí)間Δt后到達(dá)一虛擬位置，計(jì)算該虛擬位置與起點(diǎn)連線的虛擬方位角ψ1k。對(duì)于一種確定的終端編隊(duì)飛行狀態(tài)，存在n!種不同隊(duì)形和對(duì)應(yīng)關(guān)系的排列，分別計(jì)算該位置相對(duì)于起點(diǎn)的方位角ψk，并計(jì)算方位角偏轉(zhuǎn)角度絕對(duì)值|ψ1k-ψk|的總和。通過優(yōu)化方法尋找最小方位偏轉(zhuǎn)角度的編隊(duì)隊(duì)形，從而確定初始時(shí)刻與編隊(duì)形成時(shí)刻飛行器的對(duì)應(yīng)關(guān)系。

設(shè)計(jì)的隊(duì)形確定方案如下所示：

a)計(jì)算初始時(shí)刻各飛行器的機(jī)械能，找出能量最低飛行器，并按照機(jī)械能大小由低到高依次排列

i={i∈(1,2,…,n)|Ei=min(E1,E2,…,En)}

(7)

b)對(duì)機(jī)械能最低的飛行器，利用偽譜法優(yōu)化T時(shí)間內(nèi)的運(yùn)動(dòng)軌跡和終端狀態(tài)，使下式所示的指標(biāo)J最優(yōu)，其中T是選擇的編隊(duì)成形時(shí)間

(8)

c)若以機(jī)械能最低的飛行器T時(shí)刻狀態(tài)為參照，對(duì)于最終期望的編隊(duì)隊(duì)形，可以得到n種不同的編隊(duì)成形隊(duì)形，以及n!種“成形時(shí)刻—初始時(shí)刻”飛行器的一一對(duì)應(yīng)狀態(tài)，利用下式求解不同情況下的方位角偏轉(zhuǎn)角度

(9)

d)尋找方位角偏轉(zhuǎn)角度最小的情形，確定飛行器間的對(duì)應(yīng)關(guān)系

j={j∈(1,2,…,n!)|Ei=min(A1,A2,…,An!}

(10)

至此，確定了初始時(shí)刻飛行器的編號(hào)和狀態(tài)以及編隊(duì)成形時(shí)刻飛行器的編號(hào)和狀態(tài)，編隊(duì)成形過渡時(shí)間由T確定，將原編隊(duì)成形問題轉(zhuǎn)化為軌跡優(yōu)化問題。

2.2 軌跡設(shè)計(jì)

本文采用Radau偽譜法，求解各飛行器滿足動(dòng)壓、熱流、過載、終端狀態(tài)等約束的飛行軌跡，實(shí)現(xiàn)由散亂的初始狀態(tài)形成期望的編隊(duì)隊(duì)形。

在偽譜法中，取新的自變量τ∈[-1,+1]對(duì)時(shí)間做單位化，將間隔τ∈[-1,+1]劃分為由K個(gè)網(wǎng)格間隔[Tk-1,Tk]，k=1,…,K組成的網(wǎng)格，其中(T0,T1,…,TK)是網(wǎng)格點(diǎn)。

在每一個(gè)k∈[1,…,K]間隔內(nèi)部，狀態(tài)量使用N+1次拉格朗日多項(xiàng)式逼近：

(11)

對(duì)狀態(tài)作微分可以得到：

(12)

目標(biāo)函數(shù)為：

(13)

動(dòng)力學(xué)方程為：

(14)

式中，i=1,…,Nk；Ui是控制量的近似，控制量由N+1階多項(xiàng)式逼近，D(微分形式)表示如下：

(15)

式中：i=1,…,Nk，j=1,…,Nk+1，k=1,…,K。

路徑約束為：

(16)

式中：i=1,…,Nk。

積分約束為：

(17)

式中：i=1,…,Nk，j=1,…,nq。

通過上述處理將原連續(xù)問題離散逼近，進(jìn)一步求得關(guān)于離散優(yōu)化問題的梯度、雅可比矩陣等信息，通過調(diào)用IPOPT、SNOPT等二階或一階求解器解決上述離散系統(tǒng)非線性優(yōu)化問題。

3 仿真分析

不失一般性，以3個(gè)高速飛行器編隊(duì)初始控制為例驗(yàn)證本文所提出算法的有效性，飛行器初始狀態(tài)如表1所示。

表1 高速飛行器初始運(yùn)動(dòng)狀態(tài)

其中，假設(shè)3個(gè)飛行器的氣動(dòng)參數(shù)相同，其質(zhì)量、參考面積、升阻力系數(shù)參照CAV-H。3個(gè)飛行器位于同一高度上，且編隊(duì)時(shí)的緯度相同，沿不同經(jīng)度成一字形排列，相鄰兩個(gè)飛行器間的距離為10km。

令最低能量飛行器偽譜法求解的飛行時(shí)間T=30s，控制量為攻角、傾側(cè)角，設(shè)置范圍為0～30°，±90°，攻角和傾側(cè)角變化率范圍為±10(°)/s和±40(°)/s，采用GPOPS進(jìn)行離散化處理，路徑約束包括過載、熱流、動(dòng)壓約束，范圍分別為0～10g、0～104kW/m2、0～100kPa，相對(duì)求解精度設(shè)置為1e-4，并調(diào)用IPOPT求解器求解如上非線性規(guī)劃問題，仿真結(jié)果如圖1所示。

圖1 飛行器狀態(tài)變化歷程圖

圖1給出了3個(gè)飛行器狀態(tài)變化歷程，其中實(shí)線代表最低能量飛行器，點(diǎn)劃線為中間能量飛行器，虛線為最高能量飛行器。由結(jié)果可知，在高度上三者的變化相對(duì)較為平緩，其中點(diǎn)劃線所對(duì)應(yīng)飛行器平滑過渡到實(shí)線的高度上，而虛線則相對(duì)而言有一個(gè)躍起的過程。在經(jīng)緯度平面上，虛線相對(duì)實(shí)線有聚攏的趨向，點(diǎn)劃線則相對(duì)實(shí)線有一定的偏離。最終三者能夠形成期望的隊(duì)形，其三維圖像如圖1(c)所示，三角表示編隊(duì)成形時(shí)的位置。在圖1(d)中，可以明顯地看到虛線所對(duì)應(yīng)的飛行器在軌跡末段有較大的機(jī)動(dòng)減速過程，而點(diǎn)劃線所對(duì)應(yīng)飛行器速度的變化集中于過渡前段，這與高度的變化規(guī)律是一致的。

攻角、傾側(cè)角變化歷程如圖2所示。由圖2(a)可知，點(diǎn)劃線對(duì)應(yīng)飛行器的大攻角集中于軌跡前段，虛線對(duì)應(yīng)飛行器集中于軌跡后段，與上述關(guān)于機(jī)動(dòng)的分析一致。由圖2(a)，可以看到在前20s范圍內(nèi)，虛線與點(diǎn)劃線的變化基本上相反。根據(jù)初始時(shí)刻對(duì)應(yīng)飛行器的相對(duì)位置可以給出解釋，即虛線對(duì)應(yīng)飛行器過于遠(yuǎn)離最低能量飛行器，而點(diǎn)劃線飛行器則過于靠近，因此形成了2種不同的傾側(cè)角變化方式。

圖2 攻角和傾側(cè)角變化歷程圖

表2給出了飛行過程中過載、動(dòng)壓、熱流約束的最大值。由計(jì)算結(jié)果可知，過載、動(dòng)壓、熱流均滿足過程約束要求，其中動(dòng)壓、熱流密度相對(duì)較小，而最大過載相則對(duì)較大，在7.09g左右，結(jié)合圖2可知，是由于虛線對(duì)應(yīng)的飛行器在軌跡末段機(jī)動(dòng)所致。

表2 飛行過程約束最大值

在飛行器不同初始狀態(tài)下，對(duì)本文所提方法的魯棒性進(jìn)行驗(yàn)證，認(rèn)為初始時(shí)刻飛行器位置、速度在基本初始狀態(tài)下有如表3所示的偏差分布。

表3 初始運(yùn)動(dòng)狀態(tài)偏差分布

在飛行器不同初始狀態(tài)下，對(duì)本文所提方法的魯棒性進(jìn)行驗(yàn)證，認(rèn)為初始時(shí)刻飛行器位置、速度在基本初始狀態(tài)下有如表3所示的偏差分布。同時(shí)，為更好的表現(xiàn)仿真打靶結(jié)果，定義多個(gè)飛行器系統(tǒng)的隊(duì)形偏差為：

(18)

式中：Hi0，λi0和φi0是由當(dāng)前時(shí)刻最低能量飛行器和理想隊(duì)形為基準(zhǔn)確定的其他飛行器的期望高度、經(jīng)度、緯度，Hi，λi和φi為第i個(gè)飛行器的高度、經(jīng)度和緯度。在達(dá)到期望隊(duì)形后，上述系統(tǒng)高度偏差ΔH、經(jīng)度偏差Δλ和緯度偏差Δφ應(yīng)當(dāng)收斂為0。

100組仿真打靶結(jié)果如圖3所示。

圖3 打靶仿真結(jié)果

可以看到，對(duì)于給定初始狀態(tài)偏差條件下的3個(gè)飛行器，經(jīng)過上述方法后均可形成最終的穩(wěn)定隊(duì)形。圖3(a)中所有曲線均在30s處收斂于0，即滿足最終隊(duì)形的同一高度要求。圖3(b)中，偏差同樣在30s處收斂于0，滿足穩(wěn)定狀態(tài)下的經(jīng)度要求。在圖3(c)中，可以看到絕大部分曲線呈緩慢下降趨勢，在30s處收斂于0，極少部分曲線其偏差變化集中于最后5s以內(nèi)，但總體而言均同樣滿足穩(wěn)定狀態(tài)下的緯度要求。

4 結(jié)論

研究了一種考慮多種路徑約束的高速飛行器編隊(duì)隊(duì)形快速成形設(shè)計(jì)方法，利用最低能量原則和最小偏轉(zhuǎn)角度原則確定了編隊(duì)形成時(shí)期望編隊(duì)隊(duì)形的位置和“成形時(shí)刻—初始時(shí)刻”飛行器的一一對(duì)應(yīng)關(guān)系，避免了初始狀態(tài)與期望隊(duì)形偏差過大時(shí)問題無解的情況，從而將編隊(duì)成形問題轉(zhuǎn)化為軌跡設(shè)計(jì)問題，并基于偽譜法對(duì)多飛行器軌跡設(shè)計(jì)問題進(jìn)行求解。仿真驗(yàn)證了該編隊(duì)控制方法的有效性，結(jié)果表明多個(gè)高速飛行器最終能夠形成期望的穩(wěn)定隊(duì)形，打靶結(jié)果表明該方法具有一定的魯棒性，能夠適應(yīng)一定范圍內(nèi)的初始狀態(tài)。