劉哲,陸浩然,鄭偉,*,聞國光,王奕迪,周祥

1. 國防科技大學(xué) 空天科學(xué)學(xué)院,長沙 410073

2. 北京宇航系統(tǒng)工程研究所,北京 100076

3. 北京交通大學(xué) 理學(xué)院,北京 100093

高超聲速滑翔飛行器是一種在臨近空間飛行,最大速度可達(dá)馬赫數(shù)20以上的升力式飛行器,具有速度快、側(cè)向運動能力強、精度高、強突防的特點。然而隨著反導(dǎo)武器防御系統(tǒng)的不斷完善和發(fā)展,單個飛行器的生存能力下降,傳統(tǒng)的多彈同時飽和攻擊方式突防能力有限。因此,高超聲速飛行器集群中的個體可通過精確協(xié)同控制再入總時間及終端狀態(tài)的方式,多角度、多批次抵近目標(biāo),從而實現(xiàn)“探測-打擊-評估”一體化協(xié)同再入,有效提升集群的綜合效能。

相關(guān)學(xué)者從2005年起開始對多彈協(xié)同算法的研究。文獻(xiàn)[1]最早設(shè)計了固定總時間約束下的單彈最優(yōu)制導(dǎo)律(ITCG),實現(xiàn)了多彈同時到達(dá)目標(biāo)點；文獻(xiàn)[2-3]分別提出了“領(lǐng)彈-從彈”及雙層協(xié)同架構(gòu),后續(xù)的多彈協(xié)同末制導(dǎo)方法多基于這兩種架構(gòu),協(xié)同方式主要集中在時間協(xié)同。近年來,文獻(xiàn)[4]首先論證了高超聲速飛行器與常規(guī)導(dǎo)彈協(xié)同的可行性,而后對于高超集群協(xié)同軌跡規(guī)劃與制導(dǎo)的研究逐漸興起。

相比于多彈、多無人機(jī)協(xié)同,滑翔飛行器集群為欠驅(qū)動系統(tǒng),控制能力有限,再入環(huán)境復(fù)雜,多彈協(xié)同方法無法直接應(yīng)用,因此多采用“標(biāo)準(zhǔn)軌跡法[5]”“預(yù)測校正法[6-9]”或優(yōu)化方法實現(xiàn)多滑翔飛行器時間協(xié)同。文獻(xiàn)[5-6]通過調(diào)節(jié)航向角誤差走廊寬度的方式改變側(cè)向運動幅度,間接控制再入總時間；文獻(xiàn)[6]在文獻(xiàn)[5]的基礎(chǔ)上,采用神經(jīng)網(wǎng)絡(luò)擬合了航向角誤差走廊邊界與到達(dá)時間的關(guān)系,動態(tài)改變走廊邊界進(jìn)行時間調(diào)節(jié)。文獻(xiàn)[7-8]基于預(yù)測校正方法,直接迭代調(diào)整控制量剖面以實現(xiàn)時間控制；文獻(xiàn)[9]利用模型預(yù)測靜態(tài)規(guī)劃(MPSP)算法替代了預(yù)測校正的迭代過程；文獻(xiàn)[10-11]將“領(lǐng)-從”模式及ITCG方法應(yīng)用到了多高超聲速飛行器的末制導(dǎo)段。上述方法存在的問題在于: ① 難以獲取再入總時間的可調(diào)范圍并確定合理的協(xié)同時間; ② 標(biāo)準(zhǔn)軌跡方法中航向角誤差走廊的增大會導(dǎo)致終端誤差的增大; ③ 預(yù) 測校正方法一般只能滿足終端位置與協(xié)同時間需求,很難同時滿足終端角度約束,迭代求解的收斂性受初值影響較大。

為滿足終端位置與角度約束,偽譜法等優(yōu)化方法將問題離散化后求解,但計算用時較長,且同時滿足協(xié)同時間約束的文獻(xiàn)相對較少。MPSP算法通過解整個時間段的協(xié)態(tài)向量來更新控制變量解,不需要直接對非線性動力學(xué)模型采用線性化假設(shè),計算效率、數(shù)值精度高,但大多需要對控制量進(jìn)行預(yù)設(shè)計以獲取初始猜想解[12-13]。

近年來,序列凸化算法以其快速收斂性得到了廣泛應(yīng)用[14-25]。文獻(xiàn)[14]利用序列二階錐規(guī)劃(SOCP)求解了再入飛行器的最速下降軌跡,但未涉及給定協(xié)同時間約束下的規(guī)劃方法；文獻(xiàn)[15]直接將協(xié)同時間加入約束中,采用平面相對運動學(xué)模型求解,未考慮實際的動力學(xué)環(huán)境；文獻(xiàn)[16-17]分別通過減弱模型非線性,避免凸化處理的方式提高算法收斂性。此外,相關(guān)學(xué)者通過在求解過程中調(diào)整信賴域半徑的方式加快收斂速度,主要包括將信賴域作為“硬約束”與“軟約束”兩類[18-22]。文獻(xiàn)[18-20]將信賴域作為硬性約束加入到求解模型中,并根據(jù)各次迭代過程的線性化誤差等指標(biāo)調(diào)整信賴域半徑,但求解給定協(xié)同時間下的規(guī)劃問題時,常出現(xiàn)存在多條最優(yōu)軌跡的情況,導(dǎo)致解會在兩條最優(yōu)軌跡間不斷振蕩,序列迭代算法難以收斂；文獻(xiàn)[21-22]分別將信賴域及線性化誤差加入罰函數(shù)中,從而在保證線性化誤差盡可能小的同時加快收斂速度,但此類方法在處理本文問題時會使得罰函數(shù)中懲罰項過多,各項的系數(shù)大小難以確定,不利于實際應(yīng)用。同時,不合適的信賴域懲罰項更新策略會導(dǎo)致算法收斂性下降,甚至搜索不到最優(yōu)解。

針對以上問題,本文首先給出了“探測-打擊-評估”一體化模式下集群再入的協(xié)同規(guī)劃方案與數(shù)學(xué)模型,進(jìn)一步將協(xié)同軌跡規(guī)劃問題轉(zhuǎn)化為3個子問題進(jìn)行求解,即：飛行器再入總時間可調(diào)范圍的快速獲取、期望協(xié)同時間的確定、給定協(xié)同時間約束下的軌跡規(guī)劃。將飛行路徑角初值剖面定義為二次曲線,并將終端經(jīng)緯度偏差、總時間偏差以及飛行路徑角相對于預(yù)設(shè)剖面的偏離程度加入罰函數(shù),提高了問題的可解性。提出了一種罰函數(shù)與信賴域系數(shù)自適應(yīng)策略,在求解過程中不斷調(diào)整目標(biāo)函數(shù)中時間項所占權(quán)重,并根據(jù)罰函數(shù)值的變化收縮信賴域,從而在飛行路徑角平滑剖面附近逐步搜索得到最優(yōu)解,提高了序列迭代算法求解目標(biāo)函數(shù)中含時間項問題的收斂性。實現(xiàn)了不同再入任務(wù)下多滑翔飛行器時間協(xié)同軌跡的快速規(guī)劃,并與偽譜法及MPSP算法求解結(jié)果進(jìn)行了對比。

1 集群時間協(xié)同再入問題描述

1.1 集群時間協(xié)同再入策略

為實現(xiàn)“探測-打擊-評估”一體化協(xié)同再入,采用兩種時間協(xié)同方案： ① 不同起始點的多個集群確定協(xié)同再入時間,同時到達(dá)目標(biāo)點,實現(xiàn)多方向飽和攻擊； ② 單個集群內(nèi)部各飛行器通過采取不同軌跡形式控制再入時間,從而使同一集群內(nèi)部成員多批次到達(dá)目標(biāo)點。其中前序飛行器負(fù)責(zé)探測,并將態(tài)勢信息向后傳遞,后序飛行器負(fù)責(zé)攻擊及效果評估,實現(xiàn)時域協(xié)同。圖1 給出了“探測-打擊-評估”一體化集群協(xié)同再入模式的示意圖。

圖1 “探測-打擊-評估”一體化協(xié)同再入

兩種協(xié)同再入方式首先均需快速確定集群內(nèi)部各飛行器的再入總時間范圍,確定協(xié)同到達(dá)時間后對各子飛行器軌跡進(jìn)行精確的時間控制。可見,高超聲速集群滑翔段時間協(xié)同對軌跡解算的快速性有較高要求,是一種新的再入軌跡規(guī)劃問題,即在滿足式(1)、式(6)、式(12)、式(13)的約束條件下,優(yōu)化得到再入總時間范圍,進(jìn)而進(jìn)行固定期望協(xié)同時間約束下的軌跡規(guī)劃。

1.2 飛行器數(shù)學(xué)模型

對于高超聲速滑翔飛行器再入集群中的第i個成員,其以能量為自變量的無量綱三自由度動力學(xué)模型為

(1)

(2)

式中：ρ為大氣密度；Sref為參考面積；m為飛行器質(zhì)量；CL和CD分別為升力系數(shù)、阻力系數(shù),其值由當(dāng)前時刻的攻角α及馬赫數(shù)Ma二維插值得到。模型的自變量為無量綱化后的能量：

e=1/r-V2/2

(3)

則在地心距為r處,飛行器的速度為

(4)

再入過程中控制量為攻角α和傾側(cè)角σ,攻角剖面一般設(shè)置為固定形式,則再入軌跡可由傾側(cè)角剖面σ(e)唯一確定。

1.3 約束條件

1.3.1 控制量約束

在進(jìn)行滑翔段軌跡規(guī)劃時,采用“固定攻角剖面+可控制傾側(cè)角剖面”的控制方式。根據(jù)工程經(jīng)驗,攻角一般定義為線性剖面的形式：

(5)

式中：αmax、αmax L/D分別為滑翔飛行器所允許的最大攻角及最大升阻比所對應(yīng)的攻角；V1、V2為攻角剖面轉(zhuǎn)折點對應(yīng)的速度值。本文限制飛行過程中升力方向一直向上,則傾側(cè)角σ的允許變化范圍可以表示為

σmin≤|σ|≤σmax

(6)

為了避免控制量振蕩對于軌跡求解收斂性的影響,采用文獻(xiàn)[14]中的控制量選取方式,轉(zhuǎn)化后的新控制量可分解為向量形式u=[u1,u2]：

u1=cosσ,u2=sinσ

(7)

則僅需對控制量的分量進(jìn)行約束：

(u1)min≤u1≤(u1)max

(8)

由于分量u1為標(biāo)量,因此式(8)滿足凸約束條件。此外控制量還需滿足：

(9)

1.3.2 過程約束

(10)

式中：KQ=9.43×10-5,標(biāo)準(zhǔn)大氣密度ρ可展開為

ρ=ρ0e-h/hs

(11)

式中:hs為基準(zhǔn)高度。

h≥hQ(e),h≥hq(e),h≥hn(e)

(12)

1.3.3 終端約束

高超聲速飛行器滑翔段結(jié)束后,為使飛行器能夠充分利用自身機(jī)動能力順利到達(dá)預(yù)定目標(biāo)點,需要對滑翔段結(jié)束點的位置、速度等狀態(tài)量進(jìn)行約束,則當(dāng)飛行器到達(dá)終端能量狀態(tài)ef時,終端狀態(tài)量需滿足：

(13)

2 考慮協(xié)同時間約束的改進(jìn)序列凸化方法

由于目前相對缺少對多滑翔飛行器協(xié)同軌跡規(guī)劃算法的研究,協(xié)同時間約束難以滿足,因此本文針對實現(xiàn)“探測-打擊-評估”一體化協(xié)同再入的任務(wù)需求,提出了協(xié)同到達(dá)時間約束下的軌跡規(guī)劃方法。首先,根據(jù)兩種集群時間協(xié)同再入形式,分別提出了協(xié)同時間的確定策略。而后利用序列凸化算法快速解算協(xié)同軌跡。圖2給出了時間協(xié)同軌跡規(guī)劃的具體流程。

圖2 時間協(xié)同軌跡規(guī)劃流程

2.1 時間協(xié)同軌跡規(guī)劃問題

2.1.1 協(xié)同時間的確定策略

首先,通過時間協(xié)同決策模塊綜合飛行器及任務(wù)目標(biāo)點的狀態(tài)信息,利用序列凸化算法快速優(yōu)化得到各飛行器滑翔段的總時間可調(diào)范圍,進(jìn)而由上層時間協(xié)調(diào)模塊根據(jù)任務(wù)場景確定協(xié)同到達(dá)時間或到達(dá)時間序列。接下來,將協(xié)同時間下發(fā)至各飛行器,各飛行器分別進(jìn)行固定期望協(xié)同時間約束下的軌跡規(guī)劃問題求解。不同任務(wù)場景下的協(xié)同時間確定策略為

1) 多個飛行器同起點同時開始滑翔,序列多波次到達(dá)同一目標(biāo)點。

步驟1利用序列凸化算法快速解算該起點到目標(biāo)點的滑翔最長時間tmax以及最短時間tmin。

步驟2按集群內(nèi)部飛行器個數(shù)確定到達(dá)時間間隔Δt,生成到達(dá)時間序列ti：

(14)

式中：n為集群中飛行器總個數(shù)；i為飛行器編號,i=1,2,…,n-1。

步驟3將協(xié)調(diào)時間指令ti下發(fā)至各個飛行器,各子飛行器分別進(jìn)行固定滑翔段總時間約束下的軌跡規(guī)劃。

2) 多個飛行器不同起始點同時開始滑翔,同一時間到達(dá)同一目標(biāo)點。

步驟1分別優(yōu)化得到各個起點到終端目標(biāo)點的滑翔最長時間tmax,i以及最短時間tmin,i。

步驟2多飛行器的協(xié)同到達(dá)時間t可依據(jù)式(15)確定,以使得多飛行器同時到達(dá)目標(biāo)點。

t∈[tmin,1,tmax,1]∩…∩[tmin,i,tmax,i]

(15)

在本文算例中,不同起點的飛行器滑翔段協(xié)同時間t按式(16)計算：

(16)

步驟3將協(xié)同時間指令下發(fā)至各子飛行器,各子飛行器進(jìn)行固定滑翔段總時間約束下的軌跡協(xié)同規(guī)劃。

2.1.2 協(xié)同時間約束下的軌跡優(yōu)化模型

針對以上任務(wù)場景需求,將協(xié)同軌跡規(guī)劃問題轉(zhuǎn)化為2個子優(yōu)化問題,即P1：再入滑翔段飛行最長、最短總時間優(yōu)化問題；P2：固定總時間約束下的軌跡規(guī)劃問題。結(jié)合各類約束條件,滑翔段軌跡優(yōu)化問題求解模型可以被初步描述為

(17)

式中：tf為滑翔段終端時間；δtf為終端時間相對于期望協(xié)同到達(dá)時間的偏差；R0為地球半徑。

2.2 考慮時間約束的預(yù)設(shè)剖面設(shè)計

將再入過程分為初始下降段及滑翔段,初始下降段傾側(cè)角設(shè)置為常值。定義剩余航程s為飛行器當(dāng)前位置到目標(biāo)點連線在地球表面投影的大圓弧長度,其表達(dá)式為

s=R0arccos[sinφsinφf+cosφcosφfcos(λf-λ)]

(18)

對時間求導(dǎo)可得

(19)

式中：Δψ為航向角誤差。由e=1/r-V2/2可得

(20)

則無量綱化剩余航程s相對于能量e的變化率為

(21)

由式(21)可知,飛行器在縱向平面內(nèi)航程s的變化主要由飛行路徑角γ、地心距r和阻力加速度D所決定。由于以能量為離散自變量進(jìn)行軌跡設(shè)計,各個離散點處e的是固定的,且r近似為1,因此速度v和阻力加速度D的變化形式相對固定,s的大小主要取決于飛行路徑角γ的剖面形式[23]。仿真實驗表明,極值總時間或滿足期望時間約束對應(yīng)的軌跡有時并不唯一,為提高算法收斂性,并且使滑翔段軌跡規(guī)劃結(jié)果盡量平滑,將飛行路徑角γ的初值剖面定義為二次曲線的形式：

(22)

式中：γ0、γf為初始點及目標(biāo)點處的飛行路徑角,則當(dāng)中間能量狀態(tài)處的飛行路徑角γmid確定后,即可確定飛行路徑角全程剖面,積分式(21)得到總航程s。迭代求取γmid,當(dāng)式(23)成立時,可認(rèn)為飛行路徑角初值剖面滿足終端航程要求。

(23)

2.3 時間協(xié)同軌跡優(yōu)化問題的凸化處理

為應(yīng)用SOCP算法,需對目標(biāo)函數(shù)中的時間項進(jìn)行凸化處理。此外,將P0中的連續(xù)非線性動力學(xué)方程線性化,轉(zhuǎn)化為離散型有限維最優(yōu)控制問題進(jìn)行求解。

2.3.1 時間項的凸化

將飛行時間轉(zhuǎn)化為模型狀態(tài)量的線性表達(dá)式,再入段末端時刻tf可表示為

(24)

(25)

將式(2)及的表達(dá)式代入式(25)可得

(26)

式中：ρr為大氣密度ρ相對于地心距r的導(dǎo)數(shù);Tr為T相對于地心距r的導(dǎo)數(shù)。

2.3.2 動力學(xué)約束的凸化

按式(7)重新選取控制量后,系統(tǒng)動力學(xué)模型可整理為線性的凸約束形式：

(27)

式中:w(x,e)為與地球旋轉(zhuǎn)有關(guān)的項,

(28)

(29)

(30)

2.3.3 問題的離散化

軌跡規(guī)劃問題為連續(xù)型無窮維非線性規(guī)劃問題,求解時需要對其進(jìn)行離散化。本文以能量為自變量,采用梯形離散化方法,在滑翔段能量變化區(qū)間[e0,ef]內(nèi)取均勻分布的N+1個離散點{e0,e1,e2,…,eN},相鄰兩個離散點的距離Δe滿足:

Δe=(ef-e0)/N

(31)

則各個離散點對應(yīng)的能量ei為

ei=e0+i·Δei=0,1,…,N

(32)

離散后的狀態(tài)量和控制量可分別表示為(x0,x1,…,xN)和{u0,u1,…,uN}。第i個離散點處控制量為ui=[cosσi,sinσi],第i個離散點處附近對應(yīng)的動力學(xué)約束可以被線性化為

(Aixi+Βiui+bi)]=0i=1,2,…,N

(33)

(34)

式中：ri為第i個離散點處的地心距;Tri,τi為線性化系數(shù),由第k-1次序列迭代的解代入得到。

2.4 固定時間約束下的序列凸化算法

若僅以式(34)作為目標(biāo)函數(shù),則會導(dǎo)致控制量求解結(jié)果不滿足約束條件[14],因此引入復(fù)合型目標(biāo)函數(shù),將終端經(jīng)緯度偏差及時間項加入罰函數(shù),提高問題求解的可行性。此外,含時間約束的優(yōu)化問題可能存在多個解,致使序列迭代難以收斂,通過罰函數(shù)及信賴域系數(shù)的自適應(yīng)調(diào)整避免該問題,提高算法的收斂性。

2.4.1 終端約束的處理

考慮到終端約束過于嚴(yán)格會導(dǎo)致沒有可行解,適當(dāng)放寬終端經(jīng)緯度的限制,在終端位置約束不能滿足的情況下接受最接近目標(biāo)點的解。定義誤差變量λe,φe,則每一次迭代求解結(jié)果中終端經(jīng)緯度均滿足：

|λ(ef)-λf|≤λe, |φ(ef)-φf|≤φe

(35)

進(jìn)一步將λe,φe加入目標(biāo)函數(shù),使求解得到的軌跡終點與目標(biāo)點經(jīng)緯度偏差最小:

J=pλe+qλφ

(36)

式中:p、q為自適應(yīng)系數(shù),取值由2.4.3節(jié)給出。

2.4.2 目標(biāo)函數(shù)中時間項的處理

針對問題P1,直接將線性化后的時間項加入到目標(biāo)函數(shù)中,即

(37)

針對問題P2,在求解模型中加入期望總飛行時間約束：

(38)

式中：Tri、τi第i個離散點處的時間項系數(shù)；Tc為協(xié)同時間；ΔT為當(dāng)前迭代過程中的時間誤差,將ΔT加入到目標(biāo)函數(shù)中,即

J=mΔT

(39)

式中:m為自適應(yīng)調(diào)整的系數(shù),其取值策略由2.4.3 節(jié)給出。當(dāng)ΔT趨于0時,序列凸化求解得到的軌跡將滿足協(xié)同時間的要求。

2.4.3 罰函數(shù)系數(shù)及信賴域自適應(yīng)調(diào)整策略

在實際求解過程中,滑翔段最長、最短總時間及固定總時間等指標(biāo)對應(yīng)的軌跡形式并不唯一,存在多個解對應(yīng)的目標(biāo)函數(shù)值相同的情況,從而導(dǎo)致在序列迭代過程的末段,求解得到的結(jié)果在2條不同的軌跡形式間不斷振蕩,難以收斂。因此本文采用一種罰函數(shù)系數(shù)及信賴域自適應(yīng)調(diào)整策略,即將飛行路徑角相對于預(yù)設(shè)平滑剖面的偏離程度加入罰函數(shù),而后不斷降低其在罰函數(shù)中的權(quán)值系數(shù),同時收縮信賴域,以保證求解過程能夠收斂,并且得到的軌跡較為平滑。飛行路徑角剖面搜索范圍的變化趨勢如圖3所示。

圖3 飛行路徑角剖面搜索范圍變化示意圖

引入式(22)對應(yīng)的飛行路徑角預(yù)設(shè)剖面作為求解過程中的軟約束,即在第i個離散點處飛行路徑角γi需滿足：

|γi-γ0,i|≤Δγi,i=1,2,…,N

(40)

在序列迭代的開始階段,取較大的n值,此時將在圖2中的“Ⅰ”區(qū)域內(nèi)搜索可行飛行路徑角剖面,求解的主要指標(biāo)為軌跡的平滑性。隨著序列迭代的進(jìn)行,逐步降低n值,飛行路徑角預(yù)設(shè)剖面對軌跡求解的影響逐漸減弱,搜索區(qū)域?qū)U(kuò)大至區(qū)域“Ⅱ”,有利于在平滑的初始軌跡附近搜索得到最優(yōu)解。則問題P1、P2的目標(biāo)函數(shù)可表示為

(41)

式中:p、q、m、n為自適應(yīng)調(diào)整的系數(shù),第k+1次迭代求解時各系數(shù)表達(dá)式為

(42)

由于動力學(xué)方程及各類約束的線性化過程基于小擾動假設(shè),因此只有待求解變量在原參考點附近取值,即相鄰兩次求解結(jié)果滿足式(43)所示的信賴域約束時,離散后的求解模型才是對原非線性問題的有效近似。若信賴域半徑過小,則求解過程中搜索的范圍有限,可能收斂不到最優(yōu)解；若信賴域半徑過大,雖然易于找到最優(yōu)解,但對于存在多個最優(yōu)解的情況,序列凸化算法難以收斂,因此采用一種信賴域自適應(yīng)收縮策略,定義第次迭代求解時的信賴域約束為

|x(k+1)-x(k)|≤η(k)δ

(43)

式中:x=[r,λ,φ,γ,ψ]T;δ為各狀態(tài)量的初始信賴域[δr,δλ,δφ,δγ,δψ]T;η(k)為隨序列迭代次數(shù)改變的系數(shù),其表達(dá)式為

(44)

式中:ω=1-|(J(k-1)-J(k-2))/J(k-2)|;ε為預(yù)先設(shè)定的閾值,0<ε<1/2;η(k)的選取由收斂系數(shù)ξ及ω共同決定。系數(shù)ω反映了代價函數(shù)值J的變化,當(dāng)ω小于閾值ε時,認(rèn)為代價函數(shù)J的值下降不明顯,解已經(jīng)趨于最優(yōu)解附近,為加快收斂速度,避免因存在多個最優(yōu)解導(dǎo)致求解結(jié)果在多條軌跡間振蕩的情況,引入收斂系數(shù)ξ,滿足0<ξ<1。ξ值過大或過小會導(dǎo)致軌跡精度較差或收斂較慢,因此本文取ξ=0.5；當(dāng)ω大于閾值ε時,令ξ=1。

在迭代初始階段,信賴域半徑δ可選取為一個足夠大的數(shù)值,以保證能夠在較大的初始解空間范圍內(nèi)搜索得到近似最優(yōu)解。在后續(xù)迭代過程中,隨著信賴域半徑的減小,最優(yōu)解的搜索空間逐漸減小,易于序列凸化算法的收斂,相比于固定信賴域半徑方法的求解結(jié)果更接近全局最優(yōu)解。

2.5 序列迭代求解過程

在對問題凸化和離散化之后,多高超聲速飛行器滑翔段時間協(xié)同彈道規(guī)劃問題可被分為P1、P2兩個序列二階錐規(guī)劃問題,兩個子問題的求解模型可表示為

(45)

2) 固定總時間約束下的軌跡規(guī)劃問題

(46)

本文所采用的序列二階錐規(guī)劃算法流程如圖4所示,具體的求解步驟如下。

圖4 序列凸化算法流程圖

步驟1令迭代次數(shù)k=0,迭代求取飛行路徑角初值剖面,獲取軌跡的初始解x(0),并將x(0)定義為參考軌跡x(k)。

步驟2對動力學(xué)方程及各類約束進(jìn)行凸化、離散化,將參考軌跡代入計算動力學(xué)約束中的系數(shù)矩陣A(x(k))、B(x(k))、b(x(1))等,得到SOCP形式的問題求解模型。

步驟3利用SOCP算法得到最優(yōu)解x(k+1),更新參考軌跡。

步驟4判斷|x(k+1)-x(k)|<ε是否滿足,其中ε為預(yù)先設(shè)定的常值向量,若收斂條件不滿足,則更新罰函數(shù)系數(shù)m、n、p、q及信賴域半徑η(k-1)δ,轉(zhuǎn)步驟2,若收斂條件滿足,則轉(zhuǎn)步驟5。

步驟5得到優(yōu)化問題的最優(yōu)解為x(k+1),序列迭代結(jié)束。

3 仿真結(jié)果與分析

為驗證所提滑翔段時間協(xié)同軌跡規(guī)劃算法的有效性,以3個滑翔飛行器構(gòu)成的集群為例,基于2種協(xié)同任務(wù)場景進(jìn)行了數(shù)值仿真實驗。首先,針對問題P1,分別運用改進(jìn)序列凸化算法與基于非線性規(guī)劃求解器GPOPS的h-p自適應(yīng)Radau偽譜法,求解了最長、最短再入時間對應(yīng)的軌跡。而后在確定協(xié)同時間的基礎(chǔ)上求解了問題P2,分析了含時間約束的優(yōu)化問題求解過程中,本文算法與傳統(tǒng)序列凸化算法的收斂效果。最后對比了序列凸化算法與偽譜法的軌跡規(guī)劃結(jié)果。

3.1 多飛行器不同起始點協(xié)同再入

表1給出了各飛行器及目標(biāo)點的初始狀態(tài)信息,表2給出了不同算法對于飛行極值時間的求解結(jié)果。

表1 各飛行器及目標(biāo)點初始狀態(tài)

表2 序列凸化(SOCP)與偽譜法(RPM)結(jié)果對比

由表2可知,各飛行器的再入總時間可調(diào)范圍分別為144.5 s,269.6 s,239.6 s。由各飛行器與目標(biāo)點的相對位置可以看出,初始航向角對飛行器的再入總時間極值有著較大影響,飛行器2、3相比于飛行器1的可調(diào)范圍顯著增大,當(dāng)初始速度方向指向目標(biāo)時,飛行器的再入總時間將有更大的可調(diào)范圍。SOCP算法得到的總時間可調(diào)范圍與偽譜法基本一致,偽譜法求得的最短時間稍小于SOCP算法,總時間可調(diào)范圍稍大,但SOCP算法的計算時間要明顯優(yōu)于偽譜法。根據(jù)時間協(xié)同策略可確定3個飛行器的協(xié)同再入總時間為1 430 s,利用本文提出的固定期望總時間約束下的軌跡規(guī)劃方法求解再入?yún)f(xié)同軌跡,飛行器1～3的三維協(xié)同軌跡與其各自的最速下降軌跡地面投影對比如圖5和圖6 所示。表3 給出了各飛行器軌跡數(shù)值積分驗證結(jié)果的終端狀態(tài)精度情況。

圖5 協(xié)同軌跡三維圖

圖6 協(xié)同軌跡地面投影

表3 協(xié)同軌跡積分驗證結(jié)果

由圖5 可知,各飛行器在到達(dá)時間相同的基礎(chǔ)上,能夠?qū)崿F(xiàn)終端飛行路徑角及航向角的一致。圖6 中各飛行器為實現(xiàn)最短時間再入,速度方向基本保持在起始點與目標(biāo)點連線的縱向平面內(nèi),軌跡的地面投影均較為平直。為實現(xiàn)再入時間的一致,各飛行器在初始階段偏離了最速下降方向,通過側(cè)向機(jī)動的方式延緩到達(dá)時間,從而實現(xiàn)時間協(xié)同。根據(jù)表3中的積分驗證結(jié)果,終端最大高度誤差為140 m,最大再入總時間誤差為6.28 s,終端飛行路徑角與航向角誤差均小于1°,滿足精度要求。

3.2 多飛行器同起點序列式再入

假定3個滑翔飛行器從同起點再入,在滑翔段起始點分離,要求序列到達(dá)同一目標(biāo)點,驗證了算法的有效性。表4 給出了仿真初始條件,圖7～圖9給出了各飛行器軌跡對比及熱流、動壓、過載約束滿足情況,圖10以問題P1的求解過程為例,展示了本文方法對于序列凸化算法收斂性的改進(jìn)情況。

圖7 序列到達(dá)軌跡高度-能量剖面對比

圖8 序列到達(dá)軌跡地面投影對比

圖9 滑翔段動壓、過載、熱流變化情況

表4 仿真初始條件

圖7和圖8表明最長時間到達(dá)軌跡呈周期性

進(jìn)一步將本文變系數(shù)方法與傳統(tǒng)序列凸化方法、固定罰函數(shù)及信賴域收縮系數(shù)方法的迭代求解過程進(jìn)行對比。圖10給出了本文方法對于求解最短時間再入問題收斂過程的影響情況。圖中,縱軸為各方法在第次迭代過程中求解得到的滑翔段總時間相對于最優(yōu)解的絕對誤差?？梢钥闯?本文變系數(shù)方法最先收斂至最優(yōu)解；傳統(tǒng)序列凸化方法求解結(jié)果不斷振蕩,說明存在多組解對應(yīng)的目標(biāo)函數(shù)值相同,使得算法難以收斂；罰函數(shù)

圖10 不同方法最短再入時間求解結(jié)果的絕對誤差隨迭代次數(shù)變化情況

及信賴域系數(shù)均固定的方法在開始階段目標(biāo)函數(shù)值收斂趨勢與變系數(shù)方法一致,但在第3次迭代后偏離最優(yōu)解,說明目標(biāo)函數(shù)中的時間項受到了其他指標(biāo)影響,未在序列優(yōu)化后期起主要作用。

由表5可得,本文方法與偽譜法優(yōu)化得到的到達(dá)時間基本一致,終端狀態(tài)量精度均滿足要求。傳統(tǒng)序列凸化方法未收斂至最優(yōu)解,導(dǎo)致終端存在500 m高度偏差,最短再入總時間同樣出現(xiàn)較大誤差。接下來驗證了本文方法在處理給定期望總時間約束下的軌跡求解的有效性,并與偽譜法、MPSP方法進(jìn)行對比。

表5 不同方法終端精度對比

選取1 500 s作為期望總時間,分別進(jìn)行給定總時間約束條件下的軌跡解算,圖11和圖12給出了3種方法的縱向剖面及地面軌跡投影的結(jié)果對比。其中,序列凸化方法總用時4.24 s,MPSP方法用時10.82 s,偽譜法用時105.02 s。3種方法求解結(jié)果都能夠滿足期望總時間約束,終端精度均處于合理范圍。由規(guī)劃結(jié)果可知,同一滑翔段總時間對應(yīng)的軌跡形式并不唯一,三者優(yōu)化得到的地面軌跡基本一致,但序列凸化方法計算得到的縱向軌跡相比于偽譜法及MPSP方法上下跳躍幅度明顯減小,軌跡更為平滑,且解算用時更短。此外,相比于序列凸化及偽譜法,MPSP方法的收斂性對于控制量初值選取較為敏感。

圖11 高度-能量剖面對比

圖12 地面投影軌跡對比

3.3 算法適應(yīng)性測試

進(jìn)一步,分別在初始點狀態(tài)大范圍變化及小范圍攝動2種情況下,對算法的收斂性及快速性進(jìn)行仿真測試。首先,選取期望再入總時間為1 400 s,目標(biāo)終端狀態(tài)為

xf=[hf,vf,λf,φf]=[20 km,2 km/s,155°,50°]

選取多個不同的起始點狀態(tài),以驗證算法的收斂性和快速性。各次仿真的初始狀態(tài)及算法各次求解用時情況如表6所示。圖13和圖14給出了各次求解得到的軌跡積分驗證結(jié)果。

圖13 三維軌跡

表6 仿真初始條件

各次仿真優(yōu)化求解的平均時間為3.64 s,最大總時間偏差為7 s ,最大終端高度偏差為516 m,最大終端速度偏差為62 m/s,最大終端經(jīng)度偏差為0.63°,最大終端緯度偏差為0.46°。

最后,通過蒙特卡洛方法測試所提固定總時間約束下軌跡規(guī)劃算法對于初始經(jīng)緯度變化的適應(yīng)性。目標(biāo)終端狀態(tài)為

xf=[20 km,2 km/s,155°,

50°,-0.5°±1°,90°±1°]

初始點狀態(tài)的取值上限定義為

x0max=[50 100 m,6 050 m/s,90°,45°,0°,90°]

初始點狀態(tài)的取值下限定義為

x0min=[49 900 m,5 980 m/s,88°,50°,-1°,90°]

在上述初始點范圍內(nèi)隨機(jī)取值并利用序列凸化算法求解滿足再入總時間約束的可行軌跡。共進(jìn)行50次仿真,各次求解得到的軌跡通過數(shù)值積分驗證后的結(jié)果如圖15～圖17所示。

圖15 蒙特卡洛仿真結(jié)果

圖16 終端高度分布

圖17 再入總時間分布

由仿真結(jié)果統(tǒng)計可得,序列凸化算法各次仿真平均CPU用時為3.92 s；終端高度平均偏差為136.05 m,最大偏差為214.59 m；再入總時間平均偏差為2.37 s,最大時間偏差為8.04 s；平均經(jīng)緯度偏差為0.12°,終端航向角及飛行路徑角均滿足誤差要求。

4 結(jié) 論

1) 針對多高超聲速滑翔飛行器集群協(xié)同再入問題,設(shè)計了“探測-打擊-評估”一體化模式下,集群內(nèi)部時間協(xié)同的軌跡規(guī)劃方案。

2) 給出了集群協(xié)同時間的確定策略,解決了給定期望再入?yún)f(xié)同時間約束下的多飛行器軌跡快速規(guī)劃問題。

3) 提出了一種罰函數(shù)與信賴域自適應(yīng)調(diào)整策略,能夠避免存在多個最優(yōu)解情況下的軌跡振蕩問題,提高了序列凸化算法的收斂性。

4) 仿真結(jié)果表明,所提改進(jìn)序列凸化算法能夠適應(yīng)不同的時間協(xié)同再入任務(wù),收斂性受初值影響相對MPSP算法較小,求解結(jié)果相比于偽譜法軌跡更加平滑,計算速度優(yōu)于偽譜法。