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

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

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

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

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

相關(guān)學(xué)者從2005年起開始對多彈協(xié)同算法的研究。文獻(xiàn)[1]最早設(shè)計(jì)了固定總時(shí)間約束下的單彈最優(yōu)制導(dǎo)律(ITCG),實(shí)現(xiàn)了多彈同時(shí)到達(dá)目標(biāo)點(diǎn)；文獻(xiàn)[2-3]分別提出了“領(lǐng)彈-從彈”及雙層協(xié)同架構(gòu),后續(xù)的多彈協(xié)同末制導(dǎo)方法多基于這兩種架構(gòu),協(xié)同方式主要集中在時(shí)間協(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)化方法實(shí)現(xiàn)多滑翔飛行器時(shí)間協(xié)同。文獻(xiàn)[5-6]通過調(diào)節(jié)航向角誤差走廊寬度的方式改變側(cè)向運(yùn)動幅度,間接控制再入總時(shí)間；文獻(xiàn)[6]在文獻(xiàn)[5]的基礎(chǔ)上,采用神經(jīng)網(wǎng)絡(luò)擬合了航向角誤差走廊邊界與到達(dá)時(shí)間的關(guān)系,動態(tài)改變走廊邊界進(jìn)行時(shí)間調(diào)節(jié)。文獻(xiàn)[7-8]基于預(yù)測校正方法,直接迭代調(diào)整控制量剖面以實(shí)現(xiàn)時(shí)間控制；文獻(xiàn)[9]利用模型預(yù)測靜態(tài)規(guī)劃(MPSP)算法替代了預(yù)測校正的迭代過程；文獻(xiàn)[10-11]將“領(lǐng)-從”模式及ITCG方法應(yīng)用到了多高超聲速飛行器的末制導(dǎo)段。上述方法存在的問題在于: ① 難以獲取再入總時(shí)間的可調(diào)范圍并確定合理的協(xié)同時(shí)間; ② 標(biāo)準(zhǔn)軌跡方法中航向角誤差走廊的增大會導(dǎo)致終端誤差的增大; ③ 預(yù) 測校正方法一般只能滿足終端位置與協(xié)同時(shí)間需求,很難同時(shí)滿足終端角度約束,迭代求解的收斂性受初值影響較大。

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

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

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

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

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

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

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

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

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

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

(1)

(2)

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

e=1/r-V2/2

(3)

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

(4)

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

1.3 約束條件

1.3.1 控制量約束

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

(5)

式中：αmax、αmax L/D分別為滑翔飛行器所允許的最大攻角及最大升阻比所對應(yīng)的攻角；V1、V2為攻角剖面轉(zhuǎn)折點(diǎ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)點(diǎn),需要對滑翔段結(jié)束點(diǎn)的位置、速度等狀態(tài)量進(jìn)行約束,則當(dāng)飛行器到達(dá)終端能量狀態(tài)ef時(shí),終端狀態(tài)量需滿足：

(13)

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

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

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

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

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

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

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

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

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

(14)

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

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

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

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

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

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

(15)

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

(16)

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

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

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

(17)

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

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

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

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

(18)

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

(19)

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

(20)

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

(21)

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

(22)

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

(23)

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

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

2.3.1 時(shí)間項(xiàng)的凸化

將飛行時(shí)間轉(zhuǎn)化為模型狀態(tài)量的線性表達(dá)式,再入段末端時(shí)刻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)的項(xiàng),

(28)

(29)

(30)

2.3.3 問題的離散化

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

Δe=(ef-e0)/N

(31)

則各個(gè)離散點(diǎn)對應(yīng)的能量ei為

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

(32)

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

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

(33)

(34)

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

2.4 固定時(shí)間約束下的序列凸化算法

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

2.4.1 終端約束的處理

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

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

(35)

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

J=pλe+qλφ

(36)

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

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

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

(37)

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

(38)

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

J=mΔT

(39)

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

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

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

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

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

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

(40)

在序列迭代的開始階段,取較大的n值,此時(shí)將在圖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í)各系數(shù)表達(dá)式為

(42)

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

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

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

2.5 序列迭代求解過程

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

(45)

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

(46)

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

圖4 序列凸化算法流程圖

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

步驟2對動力學(xué)方程及各類約束進(jìn)行凸化、離散化,將參考軌跡代入計(jì)算動力學(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é)果與分析

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

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

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

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

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

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

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

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

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

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

3.2 多飛行器同起點(diǎn)序列式再入

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

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

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

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

表4 仿真初始條件

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

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

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

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

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

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

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

圖11 高度-能量剖面對比

圖12 地面投影軌跡對比

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

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

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

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

圖13 三維軌跡

表6 仿真初始條件

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

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

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

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

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

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

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

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

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

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

圖16 終端高度分布

圖17 再入總時(shí)間分布

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

4 結(jié) 論

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

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

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

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