張利民 李 茜,2 李偉勛

1.中原工學(xué)院 電子信息學(xué)院， 鄭州 451191

2.天津工業(yè)大學(xué) 機(jī)械工程學(xué)院，天津 300387

3.天津職業(yè)技術(shù)師范大學(xué) 理學(xué)院，天津 300222

0 引言

巡航導(dǎo)彈攻擊方向和時間的協(xié)調(diào)對于提高突防概率有著特殊意義[1-3]。近幾年，同時考慮攻擊角度和時間的制導(dǎo)規(guī)律設(shè)計(jì)成為研究熱點(diǎn)[4-6]。例如，2016年，廖沫等[7]分析了戰(zhàn)術(shù)導(dǎo)彈協(xié)同作戰(zhàn)模式和關(guān)鍵技術(shù)，重點(diǎn)針對導(dǎo)彈協(xié)同任務(wù)規(guī)劃技術(shù)進(jìn)行分析和研究，提出三級結(jié)構(gòu)的導(dǎo)彈協(xié)同任務(wù)規(guī)劃體系，并基于V圖和A*算法進(jìn)行了優(yōu)化設(shè)計(jì)。2018年，李文等[8]基于李雅普諾夫方法推導(dǎo)了帶指定攻擊角度的時間控制導(dǎo)引律,使導(dǎo)彈按指定落角和落向以指定攻擊時間打擊目標(biāo)，并借助Dubins最短路徑，彈道軌跡邊界分析了多約束條件下可指定攻擊時間的取值范圍。2019年，吳放等[9]為解決速度變化對導(dǎo)彈攻擊時間和攻擊角度的控制問題，提出一種基于成型理論和非奇異終端滑模理論的控制制導(dǎo)律，并證明了其Lyapunov穩(wěn)定性。2020年，許志等[10]針對導(dǎo)彈速度不可控且變化規(guī)律未知條件下時間與角度控制的問題，提出了一種分段解析迭代的導(dǎo)彈速度預(yù)測及在線航程修正方法。

PSO(Particle Swarm Optimization)是一種基于群智能的全局優(yōu)化方法。算法根據(jù)粒子在解空間中所處的情況進(jìn)行搜索，沒有復(fù)雜的尋優(yōu)操作，易實(shí)現(xiàn)，且需要調(diào)節(jié)的參數(shù)少，收斂速度快，已廣泛應(yīng)用于各種優(yōu)化問題[11-15]。

本文提出以Gauss偽譜法(Gauss pseudospectral method，GPM)思想基礎(chǔ)，直接基于非線性導(dǎo)引運(yùn)動學(xué)模型，將控制時域離散化為有限時間段，在考慮過程約束和禁飛區(qū)約束等約束條件下，將具有指定時間和終端入射角的反艦導(dǎo)彈最優(yōu)航跡規(guī)劃最優(yōu)控制問題轉(zhuǎn)化為一個PSO能求解的非線性規(guī)劃問題，然后采用混沌PSO算法求解該問題，另外，本文還對傳統(tǒng)過載目標(biāo)函數(shù)進(jìn)行了改進(jìn)，降低了末段控制需求。

1 導(dǎo)彈-目標(biāo)相對運(yùn)動方程

導(dǎo)彈-目標(biāo)的相對運(yùn)動關(guān)系如圖1所示。

圖1 導(dǎo)彈-目標(biāo)相對位置圖

本文假定所攻擊的目標(biāo)靜止，導(dǎo)彈速度為固定值V，選擇坐標(biāo)系XOY為地面坐標(biāo)系的水平面。

圖1中，(Xm,Ym)和γm表示導(dǎo)彈當(dāng)前的位置和航向角，A表示其法相加速度, 始終與導(dǎo)彈速度V垂直來改變其方向，(Xt,Yt)和γt表示攻擊的目標(biāo)位置和航向角。根據(jù)圖1，導(dǎo)彈-目標(biāo)的質(zhì)點(diǎn)模型課表示為[16]：

(1)

式中，(x,y)為導(dǎo)彈的位置，t為時間，γ為導(dǎo)彈航向角，A為導(dǎo)彈法向加速度，u為導(dǎo)彈法向歸一化過載。

2 基于改進(jìn)PSO的最優(yōu)航跡規(guī)劃

2.1 混沌PSO算法

PSO算法首先初始化為一群隨機(jī)粒子，再迭代尋找其最優(yōu)解。在迭代過程中，每一個粒子跟蹤兩個極值并更新自己，其中一個極值是粒子本身所找到的最優(yōu)解，即個體極值；另一個是種群當(dāng)前所尋到的最優(yōu)解，即全局極值[11-13]。

但是PSO優(yōu)化易陷入局部極值點(diǎn)，可以利用混沌搜索，跳出局部最優(yōu)[17-18]。其基本思想是：采用混沌序列對粒子的速度和位置進(jìn)行初始化，此方法不僅不改變粒子群的隨機(jī)性，還能提高種群的多樣性和粒子的遍歷性；以當(dāng)前群最優(yōu)解為基礎(chǔ)產(chǎn)生混沌序列，并用序列中的最優(yōu)位置代替粒子群中的一個位置，產(chǎn)生局部最優(yōu)解的鄰域點(diǎn)，這樣有助于惰性粒子逃離局部最優(yōu)解。

在混沌搜索的算法中，為了保持種群的多樣性，加強(qiáng)搜索的分散性，在保留一定數(shù)量最優(yōu)微粒的同時，算法根據(jù)群體的最佳位置動態(tài)收縮搜索區(qū)域，并在搜索區(qū)域內(nèi)隨機(jī)產(chǎn)生微粒來代替性能較差的微粒。

2.2 最優(yōu)控制問題轉(zhuǎn)換

由于導(dǎo)彈的航跡優(yōu)化屬于最優(yōu)控制問題，需將其轉(zhuǎn)換為PSO能求解的非線性規(guī)劃問題。借鑒GPM的思想，本文將最優(yōu)控制的狀態(tài)變量和控制變量在一系列Gauss點(diǎn)上離散，并以這些離散點(diǎn)為節(jié)點(diǎn)構(gòu)造Lagrange插值多項(xiàng)式對其進(jìn)行近似。通過對全局插值多項(xiàng)式求導(dǎo)，近似狀態(tài)變量對時間的導(dǎo)數(shù)，將微分方程約束轉(zhuǎn)換為一組代數(shù)約束，性能指標(biāo)中的積分項(xiàng)由Gauss積分計(jì)算。經(jīng)上述變換，可將最優(yōu)控制問題轉(zhuǎn)化為一系列具有代數(shù)約束的參數(shù)優(yōu)化問題，即非線性規(guī)劃問題(NLP)[19-20]。再利用混沌PSO算法對其進(jìn)行求解。

2.3 連續(xù)最優(yōu)控制問題的離散

考慮動力學(xué)方程

(2)

其中，x(t)∈Rn和u(t)∈Rm分別為狀態(tài)量和控制量，函數(shù)f:Rn×Rm×R→Rn。傳統(tǒng)的最優(yōu)控制問題可以寫為如下的標(biāo)準(zhǔn)形式：

(3)

滿足邊界等式約束

φ(x(t0),t0,x(tf),tf)=0

(4)

和過程不等式約束

C(x(t),u(t),t0,tf)≤0

(5)

其中，φ:Rn×R×Rn×R→Rq和C:Rn→Rr。

2.3.1 時域變換

采用GPM需要將航跡規(guī)劃問題的時間區(qū)間[t0,tf]轉(zhuǎn)換到[-1,1]:

(6)

2.3.2 全局插值多項(xiàng)式近似狀態(tài)與控制變量

(7)

以及τ0=-1作為節(jié)點(diǎn)構(gòu)成K+1階Lagrange插值多項(xiàng)式，并以此為基函數(shù)構(gòu)造狀態(tài)變量的近似表達(dá)式,即:

(8)

其中，Lagrange插值基函數(shù)為

(9)

且使節(jié)點(diǎn)處的近似狀態(tài)與實(shí)際狀態(tài)相等，即Xi=X(τi)=x(τi)。同樣，對于控制變量有

u(τk)=U(τk),(k=1,…,K)

(10)

即:

(11)

2.3.3 動力學(xué)微分方程約束

對式(8)求節(jié)點(diǎn)上微分，有

(12)

其中，微分矩陣可以離線確定,即:

(13)

從而動力學(xué)微分方程約束轉(zhuǎn)換為代數(shù)約束:

(14)

其中，k= 1,…,K。

2.3.4 終端狀態(tài)約束

狀態(tài)變量的近似表達(dá)式未包括終端時刻節(jié)點(diǎn)，終端狀態(tài)應(yīng)滿足動力學(xué)方程約束:

(15)

將式(11)離散化并用Gauss積分近似,可得:

X(τf)=X(τ0)+

(16)

其中，ωk為Gauss加權(quán)系數(shù)。

2.3.5 性能指標(biāo)函數(shù)的近似

將傳統(tǒng)的Bolza型性能指標(biāo)函數(shù)中的積分項(xiàng)用Gauss積分來近似，得到近似性能指標(biāo)函數(shù)[18]:

(17)

其中，Φ(X0,t0,Xf,tf).為邊值條件，ωk同前。

根據(jù)上述數(shù)學(xué)變換，最優(yōu)控制問題可轉(zhuǎn)換為：求離散節(jié)點(diǎn)上的狀態(tài)Xi，i=0,…,K和控制變量Uk，k= 1,…,K,以及初末時刻t0、tf(若未給定)，使得性能指標(biāo)(17)最小，并滿足終端狀態(tài)約束(16)，以及原最優(yōu)控制問題的邊界條件(4)和過程約束(5)

從而將原最優(yōu)控制問題轉(zhuǎn)化成非線性規(guī)劃問題(NLP)，即:

minF(y)，y∈RM
s.t.gj(y)≥0,j=1,2,…,p
hj(y)=0,j=1,2,…,l

(18)

其中，y為包含狀態(tài)變量、控制變量和端點(diǎn)時間的設(shè)計(jì)變量，p為不等式約束的個數(shù)，l為等式約束的個數(shù)。

2.4 導(dǎo)彈飛行約束條件

1)初始邊界條件約束

x(t0)=x0，y(t0)=y0，γ(t0)=γ0

(19)

2)終端條件約束

終端狀態(tài)約束條件為：

x(tf)=xf，y(tf)=yf，γ(tf)=γf

(20)

3)控制約束

導(dǎo)彈的法向加速度要滿足一定的約束為：

-Amax≤A≤Amax,Amax>0

(21)

4)戰(zhàn)場禁飛區(qū)約束

戰(zhàn)場禁飛區(qū)通常包括地形、敵方火力等，本文以圓盤表示為：

(x(t)-xq)2+(y(t)-yq)2>Rq,q=1,2,…,N

(22)

式中，(xq,yq)和Rq分別表示第q個禁飛區(qū)的中心和半徑，N表示禁飛區(qū)個數(shù)。

2.5 目標(biāo)最優(yōu)函數(shù)的建立

按照多導(dǎo)彈協(xié)同作戰(zhàn)的性質(zhì)，導(dǎo)彈的攻擊路徑規(guī)劃問題要滿足的要求主要是使能量消耗最小。已發(fā)表的文獻(xiàn)中一般所采用的過載最小性能指標(biāo)為：

(23)

式中，u為式(1)中導(dǎo)彈法向歸一化過載。為了保證末段軌跡的平直，降低控制需求，本文在(23)中引入一個時間加權(quán)函數(shù)f1(t)，其值隨時間增大而增大，以保證導(dǎo)彈末端平直：

(24)

其中，

(25)

2.6 約束處理

由于初始邊界條件約束和終端條件約束，分別在軌跡的初始和末端起重要作用，因此也引入一個性能指標(biāo)：

(26)

其中，

(27)

(28)

且，f2(t)和f3(t)為時變加權(quán)函數(shù)，f2(t)隨著時間增大逐漸減小，即當(dāng)軌跡離初始點(diǎn)越遠(yuǎn)其值越小，而f3(t)逐漸增大，即當(dāng)軌跡離目標(biāo)越近其值越大。

這里取

(29)

(30)

對于變量值域約束，混沌粒子群算法中有對其進(jìn)行處理，所以此處可以忽略其性能指標(biāo)函數(shù)的影響。

當(dāng)有禁飛區(qū)約束時，針對每個禁飛區(qū)的約束，引入一個性能指標(biāo)，

(31)

與之對應(yīng)的距離有關(guān)的變加權(quán)函數(shù)fq(x,y),q=1,…,N，當(dāng)前規(guī)劃位置距離禁飛區(qū)越近其值越大，否則越小，這里取

(32)

這樣，我們得到一個最終的指標(biāo)函數(shù)

J=J1+J2+J3

(33)

2.7 協(xié)同攻擊最優(yōu)航跡規(guī)劃問題

本文所從事的研究是針對地面固定目標(biāo)情形開展的。對于面移動目標(biāo)，可以采用導(dǎo)彈飛行中多次規(guī)劃的滾動時域方法，也是我們后續(xù)的工作，超出了本文的研究范疇。

導(dǎo)彈協(xié)同攻擊靜止目標(biāo)的最優(yōu)航跡規(guī)劃問題可描述為：在滿足動力學(xué)微分方程約束(1)、邊界條件 (19)和(20)、過程不等式約束(21)和(22)下，尋求最優(yōu)控制量u(t)，使性能指標(biāo)(33)最小。

3 仿真校驗(yàn)

仿真1：無禁飛區(qū)仿真

假設(shè)飛行空域無禁飛區(qū)。利用CLSPSO對4枚導(dǎo)彈協(xié)同攻擊原點(diǎn)處靜止目標(biāo)進(jìn)行航跡優(yōu)化。法向加速度最大值A(chǔ)max=45m/s2，各參數(shù)設(shè)置如表1所示。仿真結(jié)果如圖2～圖4所示。

表1 仿真參數(shù)設(shè)置

圖2 無禁飛區(qū)三導(dǎo)彈協(xié)同攻擊航跡圖

圖3 無禁飛區(qū)三導(dǎo)彈協(xié)同攻擊航向角曲線

圖4 無禁飛區(qū)三導(dǎo)彈協(xié)同攻擊法向加速度曲線

圖2中可看出，4枚導(dǎo)彈在飛行初期各自按使過載最小的圓弧形軌跡飛行，以消耗所規(guī)定的飛行時間，當(dāng)接近目標(biāo)點(diǎn)時，均以設(shè)計(jì)好的入射角度飛向目標(biāo)，其末端角度如圖3所示，末端控制量如圖4所示。可以看出導(dǎo)彈協(xié)同攻擊的位置、角度與時間都精確滿足要求，而且對應(yīng)于性能指標(biāo)J1的要求，末段彈道平直，并使導(dǎo)彈后期飛行軌跡平滑，對應(yīng)于性能指標(biāo)J2的要求，末段過載需求降低，并能精確地打擊目標(biāo)。另外，仿真條件為：Matlab2016，惠普ProOne 400，精度要求小于等于10-3m，GPM結(jié)點(diǎn)至少50個時，規(guī)劃出一條最優(yōu)航跡所需時間為0.809017s，如果采用C語言實(shí)現(xiàn)，效率可以提高2個數(shù)量級以上[14]，優(yōu)化時間短。

仿真2：有禁飛區(qū)仿真

假設(shè)仿真1中的航跡上有四個禁飛區(qū)(單位：m)：

(1)(x(t)+ 7500)2+(y(t)+ 2000)2≤10002

(2)(x(t)+ 6000)2+(y(t)-6000)2≤8002

(3)(x(t)+ 2500)2+(y(t)-2500)2≤8002

(4)(x(t)+4000)2+(y(t)+6000)2≤4802

4個圓盤如圖5中所示。引入這些約束作為過程約束，得到的仿真結(jié)果如圖6～圖8所示。

圖5 禁飛區(qū)位置圖

圖6 導(dǎo)彈飛行航跡圖

圖7 航向角曲線

圖8 法向加速度曲線

從圖7可以看出，所規(guī)劃的最優(yōu)航跡在經(jīng)過禁飛區(qū)時均有一個明顯的航向角變化過程，由于約束和性能指標(biāo)J3的存在，可成功地繞過禁飛區(qū)。而對于導(dǎo)彈3和導(dǎo)彈4，由于初始飛行方向和過載最大值的影響，其在繞過禁飛區(qū)前不能有較大的轉(zhuǎn)彎，因此在尋找最優(yōu)航跡時，最終找到一條離禁飛區(qū)比較近的最優(yōu)飛行航跡；而對應(yīng)于性能指標(biāo)J1的要求，圖6中4枚導(dǎo)彈的末段飛行彈道依然平直；對應(yīng)于性能指標(biāo)J2的要求，從圖8可以看出末端過載需求仍然相對比較低。另外，雖然有禁飛區(qū)的約束影響，但飛行時間沒有變化，而且同條件下的優(yōu)化時間顯示了這種算法的魯棒性。

4 結(jié)論

根據(jù)GPM的基本原理，以導(dǎo)引非線性模型作為動力學(xué)方程，充分考慮了導(dǎo)彈的性能約束和戰(zhàn)場環(huán)境約束，研究了基于混沌粒子群算法的多導(dǎo)彈指定入射角和時間條件下協(xié)同攻擊目標(biāo)的航跡優(yōu)化設(shè)計(jì)。通過引進(jìn)目標(biāo)加權(quán)函數(shù)，保證了末段彈道的平直。仿真結(jié)果表明，本方法所設(shè)計(jì)出的航跡能安全地繞過禁飛區(qū)，末段彈道特性理想，優(yōu)化時間短，在工程上具有一定的應(yīng)用參考價(jià)值。