劉超越， 張成

(北京理工大學(xué) 飛行器動(dòng)力學(xué)與控制教育部重點(diǎn)實(shí)驗(yàn)室, 北京 100081)

0 引言

近年來(lái)，導(dǎo)彈協(xié)同攻擊因其在探測(cè)機(jī)動(dòng)目標(biāo)、抵御威脅、穿透防御系統(tǒng)等方面比單個(gè)導(dǎo)彈系統(tǒng)具有更好的性能而受到越來(lái)越多的關(guān)注[1-2]。協(xié)同軌跡規(guī)劃問(wèn)題可看作是非線性最優(yōu)控制問(wèn)題，目前適用于求解最優(yōu)控制問(wèn)題的偽譜法已經(jīng)廣泛地應(yīng)用于無(wú)人機(jī)的協(xié)同軌跡規(guī)劃問(wèn)題中。Vera等[3]利用自適應(yīng)偽譜法生成無(wú)人機(jī)的協(xié)同軌跡，并實(shí)現(xiàn)了避障，但其協(xié)同飛行時(shí)間是預(yù)先規(guī)定的，局限性較大。白瑞光等[4]利用偽譜法求解了以一種加權(quán)和為性能指標(biāo)的兩階段軌跡規(guī)劃問(wèn)題，得到了3個(gè)無(wú)人機(jī)的最優(yōu)協(xié)同航跡，但該研究對(duì)無(wú)人機(jī)的動(dòng)力學(xué)模型進(jìn)行了較大簡(jiǎn)化，將其飛行速度設(shè)為定值，控制量只改變無(wú)人機(jī)的飛行方向。

與無(wú)人機(jī)相比，導(dǎo)彈的彈體動(dòng)力學(xué)具有大范圍快速時(shí)變的特點(diǎn)：飛行高度可高至幾十公里的高空；固體火箭發(fā)動(dòng)機(jī)的工作使導(dǎo)彈的飛行速度從0到幾馬赫之間迅速變化，彈體動(dòng)壓隨之變化幾百倍；火箭彈在快速穿越大氣層時(shí)，大氣密度、溫度等氣象參數(shù)將劇烈變化；導(dǎo)彈無(wú)法像無(wú)人機(jī)一樣作出懸停、后退等機(jī)動(dòng)飛行。這些特性決定了導(dǎo)彈動(dòng)力學(xué)方程的復(fù)雜性，使導(dǎo)彈協(xié)同軌跡規(guī)劃問(wèn)題的收斂域遠(yuǎn)小于無(wú)人機(jī)。因此，若將傳統(tǒng)的無(wú)人機(jī)協(xié)同軌跡規(guī)劃方法直接應(yīng)用到導(dǎo)彈的協(xié)同軌跡規(guī)劃問(wèn)題中，求解效率很低。王芳等[5]采用偽譜法對(duì)3個(gè)導(dǎo)彈組成的編隊(duì)進(jìn)行了協(xié)同軌跡求解，但其以導(dǎo)彈的軸向加速度和法向加速度為控制變量，對(duì)動(dòng)力學(xué)模型進(jìn)行了較大的簡(jiǎn)化，并未考慮導(dǎo)彈飛行高度變化造成的聲速、大氣密度的變化，有較大的局限性。

以上協(xié)同軌跡規(guī)劃算法均為基于集中式架構(gòu)的求解方式，中心控制節(jié)點(diǎn)同時(shí)聯(lián)立所有飛行器的動(dòng)力學(xué)方程建立最優(yōu)控制模型，所有動(dòng)力學(xué)方程基于同一時(shí)間變量，因此不需要協(xié)調(diào)各飛行器的飛行時(shí)間。但隨著協(xié)同飛行飛行器數(shù)量的增加，最優(yōu)控制模型動(dòng)力學(xué)方程的數(shù)量急劇增加，算法的求解時(shí)間呈指數(shù)式上升，因此，集中式方法研究模型中協(xié)同飛行的飛行器數(shù)量一般不能超過(guò)3個(gè)。為了提高求解導(dǎo)彈協(xié)同軌跡規(guī)劃問(wèn)題的效率，需要基于分布式架構(gòu)進(jìn)行協(xié)同軌跡規(guī)劃問(wèn)題的算法研究。

雙脈沖導(dǎo)彈具有速度快、突襲性好、射程遠(yuǎn)等優(yōu)點(diǎn)。在整個(gè)雙脈沖導(dǎo)彈的飛行過(guò)程中，可以區(qū)分為幾個(gè)飛行階段。由于雙脈沖發(fā)動(dòng)機(jī)需要兩次點(diǎn)火，發(fā)射前必須確定一些關(guān)鍵變量[6-7]，如二級(jí)推力的點(diǎn)火時(shí)間、初始射角，這些變量會(huì)顯著影響彈道性能。程仙壘等[8]采用改進(jìn)的粒子群優(yōu)化算法，研究了動(dòng)壓、過(guò)載以及終端狀態(tài)等多約束條件下的非連續(xù)助推彈道優(yōu)化問(wèn)題。但該研究中只優(yōu)化了滑翔段，上升段都是按照飛行程序俯仰角飛行，并未優(yōu)化全程彈道。關(guān)成啟等[9]研究了高斯偽譜法在雙脈沖導(dǎo)彈多階段彈道優(yōu)化問(wèn)題中的應(yīng)用，但二級(jí)推力的點(diǎn)火時(shí)間、制導(dǎo)啟動(dòng)時(shí)間無(wú)法隨著飛行任務(wù)改變，不能很好地發(fā)揮雙脈沖導(dǎo)彈的優(yōu)勢(shì)，對(duì)不同的飛行任務(wù)適應(yīng)性較弱。明超等[10]以裝備雙脈沖發(fā)動(dòng)機(jī)的導(dǎo)彈作為研究對(duì)象，以二級(jí)脈沖推力點(diǎn)火時(shí)間和飛行攻角為優(yōu)化變量，采用高斯偽譜法對(duì)雙脈沖導(dǎo)彈最大射程進(jìn)行彈道優(yōu)化研究，但只對(duì)導(dǎo)彈進(jìn)行豎直平面內(nèi)的建模，動(dòng)力學(xué)模型簡(jiǎn)單。

同時(shí)，在進(jìn)行空間發(fā)射、導(dǎo)彈試驗(yàn)和武器試驗(yàn)時(shí)，需要進(jìn)行全面的靶場(chǎng)安全分析，以保護(hù)人員、資產(chǎn)和公眾[11]。靶場(chǎng)安全模板工具包[12-13]結(jié)合地理空間信息(如資產(chǎn)位置和人口密度)，可以為制導(dǎo)武器的任務(wù)操作規(guī)劃和安全分析提供傷亡和損害估計(jì)，但無(wú)法通過(guò)主動(dòng)的軌跡規(guī)劃將武器碎片(主要是推力段結(jié)束后分離的發(fā)動(dòng)機(jī)殼體)投送至安全區(qū)域。若能提前主動(dòng)預(yù)測(cè)分離發(fā)動(dòng)機(jī)殼體的落點(diǎn)位置，則可對(duì)安全有更好的保證。

因此，對(duì)于考慮將分離的發(fā)動(dòng)機(jī)殼體投送至安全區(qū)域的雙脈沖導(dǎo)彈協(xié)同軌跡規(guī)劃問(wèn)題研究，必須保證雙脈沖導(dǎo)彈的動(dòng)力學(xué)模型充分反映真實(shí)飛行狀態(tài)，并發(fā)揮其雙脈沖推力特性的優(yōu)勢(shì)，同時(shí)提高多導(dǎo)彈協(xié)同軌跡規(guī)劃的求解效率；其次，須對(duì)分離發(fā)動(dòng)機(jī)殼體的落點(diǎn)位置進(jìn)行主動(dòng)預(yù)測(cè)和落點(diǎn)控制，以提供足夠的安全保證。

上述雙脈沖導(dǎo)彈的軌跡規(guī)劃問(wèn)題是一個(gè)終端時(shí)刻自由且?guī)в新窂郊s束的非線性最優(yōu)控制問(wèn)題。高斯偽譜法是求解最優(yōu)控制問(wèn)題數(shù)值方法的一種，它采用全局插值多項(xiàng)式在一系列拉格朗日高斯(LG)點(diǎn)上近似狀態(tài)量和控制量。與一般的配點(diǎn)法相比，高斯偽譜法能夠以較少的節(jié)點(diǎn)獲得較高的求解精度，同時(shí)，根據(jù)高斯偽譜法共軛變量映射定理[14-15]，采用高斯偽譜法轉(zhuǎn)化得到的非線性規(guī)劃(NLP)問(wèn)題的Karush-Kuhn-Tucker(KKT)條件與原最優(yōu)控制問(wèn)題1階最優(yōu)必要條件的離散形式具有一致性。因此，高斯偽譜方法因其高效、快速收斂和避免顯式數(shù)值積分等優(yōu)點(diǎn)，在求解最優(yōu)控制問(wèn)題方面得到了廣泛的應(yīng)用[16-17]。導(dǎo)彈分離發(fā)動(dòng)機(jī)殼體的落點(diǎn)位置控制問(wèn)題是一個(gè)預(yù)測(cè)控制問(wèn)題。基于神經(jīng)網(wǎng)絡(luò)的機(jī)器學(xué)習(xí)方法，既可以處理復(fù)雜非線性、時(shí)變性和不確定性問(wèn)題，又能夠融合啟發(fā)式經(jīng)驗(yàn)，計(jì)算速度快，已經(jīng)廣泛地應(yīng)用到了預(yù)測(cè)控制領(lǐng)域，實(shí)現(xiàn)了對(duì)制導(dǎo)精度和實(shí)時(shí)性的提升，從而提高制導(dǎo)系統(tǒng)的智能水平，如自主著陸與返航[18]、自主導(dǎo)航[19]、提高軌道預(yù)測(cè)精度[20]。

各導(dǎo)彈的飛行軌跡劃分為發(fā)射段、爬升段、續(xù)航段和制導(dǎo)攻擊段4個(gè)階段。為了發(fā)揮雙脈沖導(dǎo)彈可以根據(jù)飛行任務(wù)確定其二級(jí)推力點(diǎn)火時(shí)間的推力特性，提高協(xié)同飛行的性能，并提高最優(yōu)協(xié)同軌跡的求解效率，提出一種基于分布式并行偽譜- 神經(jīng)網(wǎng)絡(luò)(DGPM-ANN)算法，求解考慮將分離的發(fā)動(dòng)機(jī)殼體投送至安全區(qū)域的雙脈沖導(dǎo)彈協(xié)同軌跡規(guī)劃問(wèn)題。針對(duì)偽譜法對(duì)初始猜測(cè)敏感的問(wèn)題，在算法的迭代過(guò)程中，第1次迭代采用較少的離散點(diǎn)，快速計(jì)算最優(yōu)軌跡；各導(dǎo)彈以得到的最優(yōu)軌跡作為下一次迭代的參考軌跡進(jìn)行求解，直至達(dá)到指定的精度；各導(dǎo)彈并行地獨(dú)立求解各自的優(yōu)化軌跡，并確定各導(dǎo)彈二級(jí)推力的點(diǎn)火時(shí)間；引入飛行時(shí)間下界約束協(xié)調(diào)飛行時(shí)間，以使各導(dǎo)彈的飛行時(shí)間隨迭代收斂一致，保證協(xié)同飛行時(shí)間一致性；采用人工神經(jīng)網(wǎng)絡(luò)對(duì)分離發(fā)動(dòng)機(jī)殼體的射程預(yù)測(cè)函數(shù)進(jìn)行離線擬合，以對(duì)發(fā)動(dòng)機(jī)落點(diǎn)位置進(jìn)行快速、精確地預(yù)測(cè)，形成發(fā)動(dòng)機(jī)落點(diǎn)位置約束。

1 導(dǎo)彈協(xié)同攻擊問(wèn)題描述

1.1 協(xié)同軌跡規(guī)劃算法

協(xié)同軌跡規(guī)劃算法可以分為集中式、分布式兩種算法架構(gòu)，分別如圖1(a)和圖1(b)所示。在集中式算法架構(gòu)中，中心控制節(jié)點(diǎn)作為計(jì)算中心，需要同時(shí)聯(lián)立所有導(dǎo)彈的動(dòng)力學(xué)方程對(duì)協(xié)同軌跡進(jìn)行規(guī)劃，然后反饋至各導(dǎo)彈進(jìn)行任務(wù)分配，不需要設(shè)計(jì)協(xié)調(diào)變量來(lái)保證飛行時(shí)間的一致性，導(dǎo)彈之間也不需要通信。在分布式算法架構(gòu)中，各導(dǎo)彈均具備全局信息以及自主軌跡規(guī)劃能力，相互之間可以通信，不存在主從關(guān)系，導(dǎo)彈加入或退出編隊(duì)也更加靈活。

圖1 軌跡規(guī)劃算法架構(gòu)Fig.1 Trajectory planning algorithm architecture

本文提出的分布式并行算法結(jié)構(gòu)如圖2所示。分布于各導(dǎo)彈的計(jì)算單元并行地獨(dú)立求解各自的軌跡，實(shí)現(xiàn)底層的軌跡規(guī)劃；在上層，導(dǎo)彈之間通過(guò)提出的飛行時(shí)間一致性協(xié)調(diào)策略，使各導(dǎo)彈的飛行時(shí)間隨迭代收斂一致，以保證協(xié)同飛行時(shí)間一致性，從而實(shí)現(xiàn)上層的協(xié)調(diào)控制。

圖2 分布式并行算法結(jié)構(gòu)Fig.2 Decentralized parallel algorithm structure

1.2 軌跡階段的劃分

雙脈沖導(dǎo)彈的飛行過(guò)程可以根據(jù)其發(fā)動(dòng)機(jī)的推力特性劃分為4個(gè)階段：發(fā)射段、爬升段、續(xù)航段和制導(dǎo)攻擊段，如圖3所示。每個(gè)階段的具體表述如下：

1)發(fā)射段。始于發(fā)射點(diǎn)，結(jié)束于一級(jí)推力的關(guān)機(jī)點(diǎn)，為無(wú)控飛行狀態(tài)。一級(jí)推力可為導(dǎo)彈在飛行初期提供較高的飛行速度。發(fā)射段時(shí)長(zhǎng)固定，為一級(jí)推力的工作時(shí)間，設(shè)為[t0,t1)。

圖3 軌跡分段Fig.3 Division of trajectory

2)爬升段。始于一級(jí)推力的關(guān)機(jī)點(diǎn)，結(jié)束于二級(jí)推力的點(diǎn)火點(diǎn)。一級(jí)推力工作結(jié)束后，導(dǎo)彈繼續(xù)爬升，進(jìn)入爬升段。爬升段連接發(fā)射段與續(xù)航段，目的是將導(dǎo)彈續(xù)航段的開(kāi)始時(shí)刻調(diào)整至最佳時(shí)機(jī)。因此，爬升段的工作時(shí)間不固定，設(shè)為[t1,t2)。

3)續(xù)航段。始于二級(jí)推力的點(diǎn)火點(diǎn)，結(jié)束于二級(jí)推力的關(guān)機(jī)點(diǎn)。在續(xù)航段，二級(jí)推力將導(dǎo)彈的姿態(tài)和速度調(diào)整至最佳狀態(tài)，為適應(yīng)戰(zhàn)場(chǎng)環(huán)境、攻擊目標(biāo)做好充分的準(zhǔn)備。續(xù)航段時(shí)長(zhǎng)固定，為二級(jí)推力的工作時(shí)間，設(shè)為[t2,t3)。續(xù)航段結(jié)束時(shí)，發(fā)動(dòng)機(jī)結(jié)束工作，與彈體分離。

4)制導(dǎo)攻擊段。始于二級(jí)推力的關(guān)機(jī)點(diǎn)，結(jié)束于軌跡終點(diǎn)。在制導(dǎo)攻擊段，導(dǎo)彈不斷修正飛行路徑，在滿足約束的情況下，對(duì)地面目標(biāo)進(jìn)行攻擊。制導(dǎo)攻擊段的工作時(shí)間不固定，設(shè)為[t3,tf]。

1.3 導(dǎo)彈及分離發(fā)動(dòng)機(jī)殼體的動(dòng)力學(xué)方程

協(xié)同飛行的各導(dǎo)彈在地面坐標(biāo)系下的質(zhì)點(diǎn)動(dòng)力學(xué)方程為

(1)

導(dǎo)彈在飛行過(guò)程中受到的空氣阻力、升力和側(cè)向力的計(jì)算公式為

(2)

式中：0.5ρv2為動(dòng)壓；ρ為大氣密度；S為導(dǎo)彈參考橫截面積；CD、CL和CC分別為阻力系數(shù)、升力系數(shù)和側(cè)向力系數(shù)，均為攻角α、側(cè)滑角β和馬赫數(shù)Ma的函數(shù)。

對(duì)于本文中的導(dǎo)彈，阻力系數(shù)可看作是關(guān)于α和β的二次函數(shù)，升力系數(shù)可看作關(guān)于α的線性函數(shù)，側(cè)向力系數(shù)可看作關(guān)于β的線性函數(shù)，表示為

(3)

式中：a1=-0.382;a2=2.415；a3=0.149 8；b1=-0.464;b2=6.466；c1=-0.464;c2=6.466.

標(biāo)準(zhǔn)大氣密度ρ、聲速a表示為

(4)

式中：T為大氣溫度，T=292.6-0.010 39y+4.497×10-7y2-6.639×10-12y3.

當(dāng)續(xù)航段完成后，發(fā)動(dòng)機(jī)與彈體分離，并以無(wú)控狀態(tài)拋落至地面。分離發(fā)動(dòng)機(jī)殼體的動(dòng)力學(xué)方程為

(5)

式中：vb為發(fā)動(dòng)機(jī)速度；mb為發(fā)動(dòng)機(jī)質(zhì)量；Fbx為發(fā)動(dòng)機(jī)阻力，

(6)

ρb=1.225exp (-yb/7 201)，

(7)

Sb為發(fā)動(dòng)機(jī)的參考面積，CbD為發(fā)動(dòng)機(jī)的阻力系數(shù)；θb為彈道傾角；ψb為彈道偏角角；(xb,yb,zb)為分離發(fā)動(dòng)機(jī)殼體的位置坐標(biāo)。分離發(fā)動(dòng)機(jī)殼體的狀態(tài)向量可表示為xb=[vb,θb,ψb,xb,yb,zb]T.

1.4 約束條件

1)導(dǎo)彈在飛行的每個(gè)階段都要滿足相應(yīng)的階段約束，在整個(gè)軌跡的不同階段，階段約束有所不同，導(dǎo)彈的質(zhì)量隨著一級(jí)推力以及二級(jí)推力的工作發(fā)生變化。

發(fā)射段：

(8)

爬升段：

(9)

續(xù)航段：

(10)

制導(dǎo)攻擊段：

(11)

式中：tp為發(fā)動(dòng)機(jī)的工作時(shí)間，tp1、tp2分別為一級(jí)推力發(fā)動(dòng)機(jī)和二級(jí)推力發(fā)動(dòng)機(jī)的工作時(shí)間；P1和P2分別為一級(jí)推力發(fā)動(dòng)機(jī)和二級(jí)推力發(fā)動(dòng)機(jī)的推力；ξ1和ξ2分別為一級(jí)推力發(fā)動(dòng)機(jī)和二級(jí)推力發(fā)動(dòng)機(jī)的質(zhì)量秒流量；umax和umin分別為控制量的可行域上、下限。

2)變量的邊界約束。一般化表示為

Φmin≤Φ(x(t0),t0,x(tf),tf)≤Φmax，

(12)

式中：Φ、Φmax和Φmin分別為變量的一般化表示及其邊界約束的上、下限；tf為飛行時(shí)間。

3)過(guò)載約束。為了確保導(dǎo)彈飛行過(guò)程中的結(jié)構(gòu)安全，需要考慮過(guò)載約束：

(13)

式中：nL以及nLmax和nLmin分別為過(guò)載以及過(guò)載約束的上、下限。

4)分離發(fā)動(dòng)機(jī)殼體的落點(diǎn)約束。為了保證人員及財(cái)產(chǎn)安全，需使分離發(fā)動(dòng)機(jī)殼體掉落至安全位置。本文中安全區(qū)為位于地面的圓形區(qū)域，則對(duì)于分離發(fā)動(dòng)機(jī)殼體有如下約束：

(14)

式中：(xbf,ybf,zbf)為分離發(fā)動(dòng)機(jī)殼體的落點(diǎn)位置坐標(biāo)；(xe,0,ze)為發(fā)動(dòng)機(jī)的落地安全區(qū)中心位置坐標(biāo)；Rs為安全區(qū)半徑。

5)相鄰階段連接約束?？杀硎緸?/p>

(15)

1.5 最優(yōu)控制問(wèn)題框架

導(dǎo)彈的多階段協(xié)同軌跡規(guī)劃問(wèn)題可轉(zhuǎn)化為最優(yōu)控制問(wèn)題P1：

(16)

式中：Ψ(·)為性能指標(biāo)中的終端函數(shù)；L(·)為性能指標(biāo)中的被積函數(shù)；Nm為導(dǎo)彈的個(gè)數(shù)。

2 多階段最優(yōu)控制問(wèn)題離散化

高斯偽譜法將最優(yōu)控制問(wèn)題的狀態(tài)量和控制量在一系列LG點(diǎn)上離散，并以離散點(diǎn)為節(jié)點(diǎn)構(gòu)造拉格朗日插值多項(xiàng)式來(lái)逼近狀態(tài)量和控制量。通過(guò)對(duì)全局插值多項(xiàng)式求導(dǎo)來(lái)近似狀態(tài)量對(duì)時(shí)間的導(dǎo)數(shù)，從而將微分方程約束轉(zhuǎn)換為一組代數(shù)約束，性能指標(biāo)函數(shù)中的積分項(xiàng)由高斯求積公式近似，終端狀態(tài)由初始狀態(tài)加右端函數(shù)在整個(gè)過(guò)程的積分獲得。經(jīng)上述變換，可將最優(yōu)控制問(wèn)題轉(zhuǎn)化為具有一系列代數(shù)約束的NLP問(wèn)題：

1)時(shí)域變換。在每個(gè)飛行階段中，t∈[t0,tf]被劃分成K個(gè)時(shí)間區(qū)間[tk-1,tk)，k=1,2,…,K. 采用高斯偽譜法計(jì)算最優(yōu)控制問(wèn)題時(shí)，需要通過(guò)(17)式將每個(gè)時(shí)間區(qū)間轉(zhuǎn)換到[-1,1]：

(17)

2)狀態(tài)變量與控制變量的插值近似。高斯偽譜法需要在一系列離散點(diǎn)上對(duì)狀態(tài)量和控制量進(jìn)行全局插值多項(xiàng)式逼近。在[tk-1,tk)中，采用拉格朗日插值多項(xiàng)式作為基函數(shù)來(lái)近似狀態(tài)量和控制量：

(18)

(19)

3)微分方程的約束轉(zhuǎn)換為代數(shù)約束。對(duì)全局插值多項(xiàng)式(18)式微分后代入運(yùn)動(dòng)學(xué)微分方程中，并在LG點(diǎn)上離散，可得

(20)

4)約束條件離散化。邊界約束(12)式表示為

(21)

在LG點(diǎn)對(duì)路徑約束(13)式進(jìn)行離散化處理，得

(22)

式中：C、Cmax和Cmin分別為路徑約束的一般化表示及其約束的上、下限。

最后，最優(yōu)控制問(wèn)題P1可用高斯求積公式進(jìn)行離散化表示為P2：

(23)

至此，多階段連續(xù)最優(yōu)控制問(wèn)題被轉(zhuǎn)化為離散參數(shù)的NLP問(wèn)題。為了解決NLP問(wèn)題，本文采用了開(kāi)源的NLP求解器IPOPT. IPOPT采用內(nèi)點(diǎn)算法，可以處理具有大量等式和不等式約束的大規(guī)模稀疏非凸問(wèn)題。

3 射程預(yù)測(cè)函數(shù)和人工神經(jīng)網(wǎng)絡(luò)擬合

本節(jié)基于分離發(fā)動(dòng)機(jī)殼體的動(dòng)力學(xué)方程(5)式，構(gòu)造了分離發(fā)動(dòng)機(jī)殼體的射程預(yù)測(cè)函數(shù)，實(shí)現(xiàn)了發(fā)動(dòng)機(jī)由分離位置到著陸位置的距離預(yù)測(cè)，然后建立數(shù)據(jù)集，利用人工神經(jīng)網(wǎng)絡(luò)對(duì)該函數(shù)進(jìn)行了離線擬合。

3.1 分離發(fā)動(dòng)機(jī)殼體的射程預(yù)測(cè)函數(shù)

以t3時(shí)刻導(dǎo)彈的部分狀態(tài)量組成分離發(fā)動(dòng)機(jī)殼體的積分初值[vb0,θb0,ψb0,0,yb0,0]，采用定步長(zhǎng)數(shù)值積分算法對(duì)(5)式進(jìn)行數(shù)值積分。取積分終止條件為yb≤0 km，得到積分終止時(shí)刻tbf的分離發(fā)動(dòng)機(jī)殼體x軸、z軸坐標(biāo)sbfx和sbfz：

(24)

令射程預(yù)測(cè)函數(shù)的輸入向量為sb0=[vb0,θb0,ψb0,yb0]，輸出向量為sbf1=[sbfx,sbfz]，基于(5)式構(gòu)建分離發(fā)動(dòng)機(jī)殼體射程預(yù)測(cè)函數(shù)，表示為

sbf1=fb(sb0).

(25)

3.2 人工神經(jīng)網(wǎng)絡(luò)訓(xùn)練模型

人工神經(jīng)網(wǎng)絡(luò)可以通過(guò)學(xué)習(xí)樣本集自動(dòng)提取規(guī)則，實(shí)現(xiàn)從輸入到輸出的映射，適用于擬合復(fù)雜非線性函數(shù)，因此，它可以根據(jù)訓(xùn)練樣本預(yù)測(cè)輸出結(jié)果。本文利用人工神經(jīng)網(wǎng)絡(luò)的擬合能力，對(duì)分離發(fā)動(dòng)機(jī)殼體射程預(yù)測(cè)函數(shù)(25)式進(jìn)行函數(shù)逼近，以提高預(yù)測(cè)速度，實(shí)現(xiàn)分離發(fā)動(dòng)機(jī)殼體落點(diǎn)的實(shí)時(shí)預(yù)測(cè)，為導(dǎo)彈的在線軌跡規(guī)劃提供條件。因此，需要訓(xùn)練得到一個(gè)神經(jīng)網(wǎng)絡(luò)擬合函數(shù)sbf2=Fb(sb0)，使得對(duì)于任意的sb0:

‖sbf1-sbf2‖≤δn，

(26)

式中:δn為擬合誤差向量的上限，δn=0.5 km. 則分離發(fā)動(dòng)機(jī)殼體的落點(diǎn)位置約束(14)式可表示為

‖(x(t3)+sbf2(1)-xe,z(t3)+sbf2(2)-ze)‖2≤Rs.

(27)

導(dǎo)彈的多階段協(xié)同軌跡規(guī)劃問(wèn)題可以表示為P3：

(28)

本文采用前饋神經(jīng)網(wǎng)絡(luò)進(jìn)行函數(shù)擬合，包括1個(gè)輸入層、1個(gè)隱含層和1個(gè)輸出層。圖4展示了神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)，含有4個(gè)輸入，隱含層有35個(gè)神經(jīng)元，2個(gè)輸出。

圖4 人工神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.4 Artificial neural network structure

在神經(jīng)網(wǎng)絡(luò)的訓(xùn)練階段，為了能夠訓(xùn)練得到足夠精確的擬合函數(shù)，應(yīng)該建立充分的訓(xùn)練數(shù)據(jù)集，多次考慮不同狀態(tài)量初值來(lái)求解射程預(yù)測(cè)函數(shù)。各狀態(tài)量初值是以在飛行期間分離發(fā)動(dòng)機(jī)殼體的每個(gè)狀態(tài)量可以假定的最大值和最小值來(lái)定義的。設(shè)分離發(fā)動(dòng)機(jī)殼體質(zhì)量mb=50 kg，阻力系數(shù)CbD=0.149，參考面積Sb=0.292 m2. 各狀態(tài)量初值的上、下界以及步長(zhǎng)由表1所示。

表1 狀態(tài)量初值邊界及步長(zhǎng)Tab.1 Boundaries and step sizes of initial state values

以表1所示的狀態(tài)量初值按相應(yīng)步長(zhǎng)進(jìn)行網(wǎng)絡(luò)劃分，得到4維初值網(wǎng)格空間。以網(wǎng)格空間內(nèi)的所有節(jié)點(diǎn)對(duì)應(yīng)的值作為求解射程預(yù)測(cè)函數(shù)的初值，求解得到4 290 783個(gè)對(duì)應(yīng)的預(yù)測(cè)射程向量，即產(chǎn)生4 290 783組輸入輸出樣本。在完成預(yù)測(cè)射程向量數(shù)據(jù)采集步驟后，利用人工神經(jīng)網(wǎng)絡(luò)對(duì)采集到的數(shù)據(jù)集進(jìn)行離線訓(xùn)練，生成學(xué)習(xí)模型，并將模型的權(quán)值和閾值表示的參數(shù)存儲(chǔ)在數(shù)據(jù)庫(kù)中，從而得到射程預(yù)測(cè)函數(shù)的人工神經(jīng)網(wǎng)絡(luò)擬合函數(shù)Fb(sb0)。圖5顯示了數(shù)據(jù)訓(xùn)練模型。

圖5 數(shù)據(jù)集訓(xùn)練模型Fig.5 Dataset training model

定義擬合結(jié)果的誤差百分比ex和ez分別為

(29)

為了分析人工神經(jīng)網(wǎng)絡(luò)對(duì)以不同積分步長(zhǎng)求解射程預(yù)測(cè)函數(shù)得到的樣本數(shù)據(jù)集進(jìn)行擬合的精度以及擬合函數(shù)的求解速度，在積分步長(zhǎng)分別為0.6 s、0.8 s和1.0 s的條件下，分別求解射程預(yù)測(cè)函數(shù)得到相應(yīng)的數(shù)據(jù)集，然后利用人工神經(jīng)網(wǎng)絡(luò)分別進(jìn)行離線擬合，擬合結(jié)果參數(shù)如表2所示。從表2中可以看出：隨著積分步長(zhǎng)的增大，射程預(yù)測(cè)函數(shù)的積分節(jié)點(diǎn)減少，平均計(jì)算時(shí)間也相應(yīng)減小，但預(yù)測(cè)射程的精確度也勢(shì)必有所減?。煌瑫r(shí)，3種情況下的神經(jīng)網(wǎng)絡(luò)擬合函數(shù)的平均計(jì)算時(shí)間以及擬合函數(shù)預(yù)測(cè)射程x軸方向分量和z軸方向分量的平均誤差百分比幾乎一致，且擬合函數(shù)的預(yù)測(cè)射程非常接近射程預(yù)測(cè)函數(shù)的預(yù)測(cè)射程。積分步長(zhǎng)分別為0.6 s、0.8 s和1.0 s的情況下，采用神經(jīng)網(wǎng)絡(luò)擬合函數(shù)預(yù)測(cè)射程，預(yù)測(cè)速度分別提高了42倍、30倍和24倍。基于以上分析，本文以0.6 s作為積分步長(zhǎng)生成樣本數(shù)據(jù)集，可以得到較高的射程預(yù)測(cè)函數(shù)求解精度，并可以進(jìn)一步得到較高的神經(jīng)網(wǎng)絡(luò)擬合函數(shù)預(yù)測(cè)精度和求解速度。

表2 不同積分步長(zhǎng)對(duì)應(yīng)的數(shù)據(jù)集人工神經(jīng)網(wǎng)絡(luò)擬合結(jié)果Tab.2 ANN fitting results of data sets corresponding todifferent integral step sizes

樣本數(shù)據(jù)經(jīng)人工神經(jīng)網(wǎng)絡(luò)擬合后，擬合結(jié)果x軸方向分量和z軸方向分量的誤差百分比分別如圖6(a)和圖6(b)所示。經(jīng)驗(yàn)證，擬合結(jié)果滿足誤差約束(26)式。由此可見(jiàn)，采用人工神經(jīng)網(wǎng)絡(luò)擬合預(yù)測(cè)模型，可以提供較高的精度和更快的求解速度。

圖6 訓(xùn)練誤差Fig.6 Training error

4 分布式并行偽譜- 神經(jīng)網(wǎng)絡(luò)算法

求解軌跡規(guī)劃問(wèn)題的流程圖如圖7所示。流程圖左上方的3個(gè)方框表示了離線狀態(tài)下分離發(fā)動(dòng)機(jī)殼體的射程預(yù)測(cè)函數(shù)人工網(wǎng)絡(luò)擬合過(guò)程。右上方兩個(gè)方框表示將多階段軌跡規(guī)劃的最優(yōu)控制問(wèn)題轉(zhuǎn)化為非線性規(guī)劃問(wèn)題的過(guò)程。底部的方框表示用結(jié)合序列高斯偽譜法和人工神經(jīng)網(wǎng)絡(luò)的DGPM-ANN算法對(duì)軌跡規(guī)劃問(wèn)題進(jìn)行求解，從而得到原最優(yōu)控制問(wèn)題的解。

圖7 軌跡規(guī)劃求解流程圖Fig.7 Flowchart of trajectory planning

本文提出的DGPM-ANN算法是一種基于迭代策略解耦形式的協(xié)同軌跡求解方法。在第1次迭代中，采用較少的離散點(diǎn)，快速計(jì)算最優(yōu)軌跡；以得到的最優(yōu)軌跡作為下一次迭代的參考軌跡進(jìn)行求解，直至達(dá)到指定的精度；每一輪次的迭代中各導(dǎo)彈并行地獨(dú)立求解軌跡規(guī)劃問(wèn)題，并利用神經(jīng)網(wǎng)絡(luò)擬合函數(shù)對(duì)分離發(fā)動(dòng)機(jī)殼體的落點(diǎn)位置進(jìn)行計(jì)算，以對(duì)發(fā)動(dòng)機(jī)落點(diǎn)位置進(jìn)行約束；在每次迭代中引入飛行時(shí)間下界約束來(lái)協(xié)調(diào)各導(dǎo)彈的飛行時(shí)間，以使所有導(dǎo)彈的飛行時(shí)間隨迭代趨于一致；所有導(dǎo)彈的飛行時(shí)間達(dá)到一致后，DGPM-ANN算法完成。

由于在分布式算法的每一次迭代中，各導(dǎo)彈均獨(dú)立求解最優(yōu)軌跡，所以不同導(dǎo)彈的飛行時(shí)間可能不同。為了保證不同導(dǎo)彈飛行時(shí)間的一致性，迭代收斂時(shí)必須滿足以下條件：

(30)

(31)

圖8 DGPM-ANN算法過(guò)程

5 數(shù)值仿真

為了驗(yàn)證所提DGPM-ANN算法的高效性和有效性，通過(guò)兩個(gè)仿真算例，利用一種集中式偽譜- 人工神經(jīng)網(wǎng)絡(luò)(CGPM-ANN)算法與本文提出的DGPM-ANN算法分別求解雙脈沖導(dǎo)彈的多階段協(xié)同軌跡規(guī)劃問(wèn)題，然后將兩種算法的仿真結(jié)果進(jìn)行比較。CGPM-ANN算法同時(shí)對(duì)所有導(dǎo)彈進(jìn)行協(xié)同軌跡規(guī)劃問(wèn)題求解，并采用與CGPM-ANN算法中相同的神經(jīng)網(wǎng)絡(luò)擬合函數(shù)，以續(xù)航段末端狀態(tài)為分離發(fā)動(dòng)機(jī)殼體的初始狀態(tài)，預(yù)測(cè)發(fā)動(dòng)機(jī)落點(diǎn)。采用通用偽譜最優(yōu)控制軟件GPOPS實(shí)現(xiàn)算法，GPOPS軟件使用默認(rèn)的NLP求解程序IPOPT，仿真在一臺(tái)使用英特爾i7-4790@3.6GHz處理器的計(jì)算機(jī)上運(yùn)行。仿真參數(shù)和約束參數(shù)分別如表3和表4所示，表中：αmin和αmax分別為攻角的可行域上、下限；βmin和βmax分別為側(cè)滑角的可行域上、下限；θ0min和θ0max分別為初始彈道傾角的可行域上、下限，各導(dǎo)彈根據(jù)任務(wù)分別確定其最優(yōu)初始彈道傾角θ0. 建立仿真場(chǎng)景如下：4枚導(dǎo)彈同時(shí)從表5所示的不同地點(diǎn)發(fā)射，對(duì)坐標(biāo)位置為(120 km,0 km,60 km)的地面目標(biāo)進(jìn)行協(xié)同攻擊。DGPM-ANN算法與CGPM-ANN算法中各導(dǎo)彈的每個(gè)階段離散點(diǎn)個(gè)數(shù)均為40，GPOPS軟件采用默認(rèn)參數(shù)。

表3 仿真參數(shù)表Tab.3 Simulation parameters

表4 約束參數(shù)Tab.4 Constraint parameters

表5 發(fā)射點(diǎn)位置Tab.5 Launch positions km

圖9 算例1最優(yōu)軌跡Fig.9 Optimal trajectory in Case 1

算例1以飛行時(shí)間最短為仿真目標(biāo)，DGPM-ANN算法和CGPM-ANN算法得到的最優(yōu)軌跡分別如圖9(a)～圖9(d)所示，相鄰階段的連接節(jié)點(diǎn)用黑點(diǎn)標(biāo)記。導(dǎo)彈軌跡用實(shí)線表示；以求解得到的最優(yōu)軌跡導(dǎo)彈續(xù)航段末端狀態(tài)為分離發(fā)動(dòng)機(jī)殼體初始狀態(tài)，根據(jù)(5)式進(jìn)行數(shù)值積分，得到分離的發(fā)動(dòng)機(jī)殼體軌跡，由虛線表示；發(fā)動(dòng)機(jī)落點(diǎn)安全區(qū)由黑色圓圈表示。由圖9可知，DGPM-ANN算法與CGPM-ANN算法均實(shí)現(xiàn)了導(dǎo)彈對(duì)指定地面目標(biāo)的協(xié)同攻擊，且各分離發(fā)動(dòng)機(jī)殼體均落到了安全區(qū)域。

DGPM-ANN算法得到的算例1最優(yōu)飛行參數(shù)如圖10所示。由圖10可知，為了使飛行時(shí)間達(dá)到最短，發(fā)動(dòng)機(jī)分離后，各導(dǎo)彈攻角迅速減小，使導(dǎo)彈獲得負(fù)升力，同時(shí)在側(cè)滑角的控制下，各導(dǎo)彈的飛行方向立即指向目標(biāo)位置，制導(dǎo)攻擊段導(dǎo)彈無(wú)明顯側(cè)向機(jī)動(dòng)，以快速攻擊地面目標(biāo)。整個(gè)飛行過(guò)程中的過(guò)載滿足約束。

算例2以控制量最省為仿真目標(biāo)，DGPM-ANN算法和CGPM-ANN算法得到的最優(yōu)軌跡分別如圖11(a)～圖11(d)所示，相鄰階段的連接節(jié)點(diǎn)用黑點(diǎn)標(biāo)記。由圖11可知，兩種算法均實(shí)現(xiàn)了導(dǎo)彈對(duì)指定地面目標(biāo)的協(xié)同攻擊，且各分離發(fā)動(dòng)機(jī)殼體均落到了安全區(qū)域。

算例2 DGPM-ANN算法得到的最優(yōu)飛行參數(shù)如圖12所示。由圖12可知：為了使各導(dǎo)彈控制量最省，在飛行階段的前期，導(dǎo)彈只在攻角和側(cè)滑角作用下進(jìn)行較小機(jī)動(dòng)，以保證分離發(fā)動(dòng)機(jī)殼體被拋落至安全區(qū)域；制導(dǎo)攻擊段的后期，導(dǎo)彈進(jìn)入空氣密度較大的高度后，在攻角和側(cè)滑角作用下對(duì)目標(biāo)進(jìn)行協(xié)同攻擊。整個(gè)飛行過(guò)程中的過(guò)載滿足約束。

兩個(gè)算例中各導(dǎo)彈的優(yōu)化參數(shù)，分別如表6和表7所示。由表6和表7可知，兩個(gè)算例中DGPM-ANN算法的最優(yōu)性能指標(biāo)分別為CGPM-ANN算法的99.5%和92.4%，獲得了更好的性能指標(biāo)。這是因?yàn)椴捎肅GPM-ANN算法求解該多階段的最優(yōu)控制問(wèn)題時(shí)，所有導(dǎo)彈的動(dòng)力學(xué)方程被聯(lián)立為大規(guī)模方程組，所有方程共用同一時(shí)間變量，各導(dǎo)彈無(wú)法獨(dú)立地確定其最適合協(xié)同飛行的最優(yōu)二級(jí)推力點(diǎn)火時(shí)間，因此所有導(dǎo)彈的各階段時(shí)間節(jié)點(diǎn)均相同。作為對(duì)比，DGPM-ANN算法求解時(shí)，各導(dǎo)彈可根據(jù)飛行任務(wù)和協(xié)同飛行任務(wù)的性能指標(biāo)獨(dú)立地優(yōu)化其最優(yōu)軌跡，因此各導(dǎo)彈擁有獨(dú)立的時(shí)間變量，可以確定各自更適合協(xié)同飛行任務(wù)的最優(yōu)二級(jí)推力點(diǎn)火時(shí)間，從而更好地實(shí)現(xiàn)各導(dǎo)彈發(fā)動(dòng)機(jī)參數(shù)與彈道參數(shù)的共同優(yōu)化，使協(xié)同飛行獲得更好的性能指標(biāo)。

圖10 算例1 DGPM-ANN算法最優(yōu)飛行參數(shù)Fig.10 Optimal flight parameters in Case 1 of DGPM-ANN algorithm

圖11 算例2最優(yōu)軌跡Fig.11 Optimal trajectory in Case 2

圖12 算例2 DGPM-ANN算法最優(yōu)飛行參數(shù)Fig.12 Optimal flight parameters in Case 2 of DGPM-ANN algorithm

表6 算例1優(yōu)化參數(shù)Tab.6 Optimal parameters in Case 1

表7 算例2優(yōu)化參數(shù)Tab.7 Optimal parameters in Case 2

由表6和表7可知：兩個(gè)算例中，DGPM-ANN算法的求解時(shí)間分別為CGPM-ANN算法的40.21%和26.8%，具有更高的求解效率；DGPM-ANN算法得到的各導(dǎo)彈飛行時(shí)間也滿足飛行時(shí)間一致性約束。

6 結(jié)論

本文建立了雙脈沖導(dǎo)彈多階段協(xié)同對(duì)地攻擊的軌跡規(guī)劃模型，并考慮將分離發(fā)動(dòng)機(jī)殼體投送至安全區(qū)域。通過(guò)以飛行時(shí)間最短以及飛行過(guò)程控制量最省的兩個(gè)算例仿真，得到結(jié)論如下：

1)利用人工神經(jīng)網(wǎng)絡(luò)對(duì)分離發(fā)動(dòng)機(jī)殼體的射程預(yù)測(cè)函數(shù)進(jìn)行離線訓(xùn)練，生成高精度擬合模型，大大提高了預(yù)測(cè)效率。

2)提出了基于高斯偽譜法和人工神經(jīng)網(wǎng)絡(luò)的分布式并行算法，仿真結(jié)果表明，通過(guò)引入飛行時(shí)間下界約束，使各導(dǎo)彈的飛行時(shí)間隨迭代收斂一致，實(shí)現(xiàn)了各導(dǎo)彈攻擊目標(biāo)的飛行時(shí)間一致性。

3)分布式并行算法允許各導(dǎo)彈獨(dú)立求解最優(yōu)軌跡，因而各導(dǎo)彈可以分別確定其二級(jí)推力點(diǎn)火的最優(yōu)時(shí)間，從而得到更好的性能指標(biāo)。

4)分布式并行算法具有較高的求解效率，因此對(duì)求解多導(dǎo)彈協(xié)同攻擊的軌跡規(guī)劃問(wèn)題更有潛力。