填充式防護結(jié)構(gòu)彈道極限方程的差異演化優(yōu)化

賈光輝， 姚光樂， 張帥

(1. 北京航空航天大學(xué) 宇航學(xué)院, 北京 100083； 2. 中國空間技術(shù)研究院 總體部， 北京 100094)

隨著中國載人航天工程的發(fā)展，填充式防護結(jié)構(gòu)成為航天器防護的重點[1]。在航天器的防護問題中，彈道極限方程一直是最主要的研究內(nèi)容之一，也是航天器空間碎片撞擊風(fēng)險評估失效判斷的依據(jù)[2]。為此，需要對填充式防護結(jié)構(gòu)彈道極限方程的準確建模進行研究。

文獻[3]為獲得適用于中國航天工程特定填充式防護結(jié)構(gòu)的彈道極限方程，采用窮舉法對不同預(yù)測指標[4]組合下NASA 填充式防護結(jié)構(gòu)彈道極限方程[5]進行雙待定參數(shù)修正，修正后方程的預(yù)測能力更優(yōu)。由于窮舉法具有全空間搜索能力，其結(jié)果為所枚舉變量的最優(yōu)值，因此，文獻[3]所獲得的彈道極限方程的預(yù)測能力基本上為最優(yōu)。但此方法的計算量與待定參數(shù)個數(shù)呈指數(shù)型的關(guān)系，因此不適合待定參數(shù)較多且數(shù)據(jù)量大的彈道極限方程建模。

差異演化(Differential Evolution，DE)算法[6]是一種基于群體差異的算法，具有全局搜索能力，且簡單易用，收斂速度快[7]。文獻[8]通過函數(shù)測試表明，差異演化算法在絕大多數(shù)情況下比遺傳算法等具有更強的全局搜索能力。

除了借鑒NASA 填充式防護結(jié)構(gòu)彈道極限方程之外，建立國內(nèi)填充式防護結(jié)構(gòu)的彈道極限方程具有重要意義。文獻[9]通過量綱理論建立了該彈道極限方程形式，確定該方程中包含的3個待定系數(shù)和8個待定指數(shù)之后，即完成該方程的建模工作。為了求得這些待定參數(shù)的數(shù)值，本文將此方程中的指數(shù)和系數(shù)作為優(yōu)化變量，預(yù)測指標(概率型指標和誤差型指標)[4]作為目標函數(shù)，基于國內(nèi)數(shù)據(jù)，采用差異演化算法對方程中的優(yōu)化變量進行確定。

1 差異演化算法基本原理

差異演化算法的演化過程與遺傳算法類似，也包括種群初始化、變異、交叉和選擇操作。差異演化是基于實數(shù)編碼的演化算法，因此可將個體視為實數(shù)設(shè)計變量組成的向量。差異演化算法首先創(chuàng)建種群，然后對種群的所有個體進行變異操作，接著進行離散交叉生成中間個體，最后中間個體與父個體進行競爭，較優(yōu)的個體被選中，形成數(shù)目不變的下一代新種群。

1.1 種群初始化

在D維空間隨機產(chǎn)生個體數(shù)為NP的種群，其各個個體的變量為

i=1,2,…,NP,j=1,2,…,D

(1)

1.2 變 異

變異是差異演化算法的關(guān)鍵步驟，基于個體向量差進行。設(shè)當前種群中的個體為Xi(t)，t為演化代數(shù)。從當前種群中隨機選擇3個個體Xr1(t)、Xr2(t)和Xr3(t)(r1≠r2≠r3≠i),取后2個個體向量的差(Xr2(t)-Xr3(t))，使用縮放因子F調(diào)節(jié)大小，然后與第1個個體向量Xr1(t)求和，獲得變異后的個體hi(t+1)：

hi(t+1)=Xr1(t)+F(Xr2(t)-Xr3(t))

(2)

式中：F為縮放因子，取值范圍一般為[0,2]。

1.3 交 叉

交叉操作可以增加種群的多樣性，變異后的個體hi(t+1)和種群中當前的演化個體Xi(t)通過離散交叉，生成中間個體vi(t+1)：

vij(t+1)=

(3)

式中：rand(i)為[1,D]之間的隨機整數(shù)；CR為交叉概率，取值范圍一般為[0,1]。

1.4 選 擇

中間個體vi(t+1)與當前演化個體Xi(t)的適應(yīng)度進行比較：

Xi(t+1)=

(4)

式中：f表示適應(yīng)度函數(shù)。

通過上述變異、交叉和選擇，對種群中的每個個體進行循環(huán)操作，得到下一代種群，如此演化若干代，獲得優(yōu)化問題的最優(yōu)解。

2 測試函數(shù)演化結(jié)果

對基本差異演化算法進行編程實現(xiàn)，以文獻[10]中的3個測試函數(shù)進行程序驗證。

Ackley函數(shù)(f1)：

xi∈[-32,32]

(5)

Rastrigin函數(shù)(f2)：

xi∈[-5.12,5.12]

(6)

Griwangk函數(shù)(f3)：

xi∈[-512,512]

(7)

式(5)～式(7)中,n=20，即各函數(shù)的優(yōu)化變量個數(shù)，其中Ackley為單峰函數(shù)，Rastrigin、Griwangk為多峰函數(shù)。優(yōu)化計算時，NP的取值原則參考文獻[2]：5D≤NP≤10D，種群大小NP取100，縮放因子F取0.5，交叉概率CR取0.3。

對3個函數(shù)進行差異演化計算時，在每一代的種群中選擇出適應(yīng)度最小的個體，作為最優(yōu)個體。理論上講，需要分別畫出最優(yōu)個體的各個變量變化曲線，并進行討論。鑒于各個函數(shù)都有20個變量，若逐一對變量畫圖分析，就會過于冗雜，因此只對每個函數(shù)的20個變量絕對值求和，并進行畫圖分析。

在對函數(shù)f1、f2和f3優(yōu)化時，最優(yōu)個體的變量絕對值之和與適應(yīng)度隨演化代數(shù)的變化曲線如圖1(a)～(c)所示。

由圖1可知，3個函數(shù)在種群不斷演化過程中，最優(yōu)個體逐漸到達最優(yōu)。函數(shù)f1和f3在演化代數(shù)大于300以后，變量絕對值之和與適應(yīng)度都基本為零，趨于平穩(wěn)。函數(shù)f2在代數(shù)大于1 800以后，變量絕對值之和與適應(yīng)度都基本為零，趨于平穩(wěn)。3個函數(shù)的優(yōu)化結(jié)果如表1所示。

由表1可以發(fā)現(xiàn)，函數(shù)f1、f2和f3的優(yōu)化結(jié)果與理論結(jié)果誤差小于0.000 2，因此本文的差異演化算法程序可得到正確的結(jié)果。

圖1 各函數(shù)最優(yōu)個體的變量絕對值之和與適應(yīng)度隨演化代數(shù)變化曲線Fig.1 Changing curves of sum optimal individual absolute variables and fitness with generation for each function

測試函數(shù)理論最優(yōu)值演化代數(shù)優(yōu)化結(jié)果f105000.00006f2020000.00001f305000.00012

3 綜合建模方程

3.1 一般形式

針對國內(nèi)填充式防護結(jié)構(gòu)形式，文獻[9]利用量綱理論建立了一種合適的填充式防護結(jié)構(gòu)彈道極限方程形式，稱之為綜合建模方程，其形式如下所示。

1) 低速區(qū)：V≤VL。

(8)

2) 高速區(qū)：V≥VH。

(9)

3) 中速區(qū)：VL

dc=dc,V=VL(VH-V)/(VH-VL)+

dc,V=VH(V-VL)/(VH-VL)

(10)

VL=2.6(cosθ)-0.5

VH=6.5(cosθ)-0.75

上述彈道極限方程共涉及16個物理變量,dc為臨界彈丸直徑,cm；V為彈丸速度,km/s；VL為低速區(qū)與中速區(qū)的分界速度,km/s；VH為高速區(qū)與中速區(qū)的分界速度，km/s；ρp為彈丸密度,g/cm3；tb為前板厚度,cm；σb為前板強度,ksi(1 ksi=0.145 MPa=0.145 N/mm2)；S為板間距,cm；tw為后板厚度,cm；ρw為后板密度,g/cm3；σw為后板強度,ksi；σce為玄武巖填充層的強度,ksi；σca為Kevlar填充層的強度,ksi；mce為玄武巖填充層的面密度,g/cm2；mca為Kevlar填充層的面密度,g/cm2；θ為撞擊角度,rad。此方程的待定參數(shù)共11個，包含3個待定系數(shù)和8個待定指數(shù)，待定系數(shù)為ηi(i=1,5,8)，待定指數(shù)為ηj(j=2,3,4,6,7,9,10,11)，均無量綱。

3.2 待定參數(shù)優(yōu)化模型

本文在進行綜合建模方程的待定參數(shù)優(yōu)化時，優(yōu)化的目標函數(shù)對應(yīng)多個預(yù)測指標(總體預(yù)測率、安全預(yù)測率、平均相對誤差平方和)。平均相對誤差平方和為誤差型指標，其大小由預(yù)測錯誤的實驗直徑與方程預(yù)測直徑之間的偏差確定；總體預(yù)測率和安全預(yù)測率為概率型指標，其大小是由各類樣本數(shù)和方程預(yù)測正確的個數(shù)確定?？梢?，表示概率型指標的2個函數(shù)為整數(shù)運算結(jié)果，不連續(xù)。同時，這3個目標函數(shù)值的確定依賴于綜合建模方程和實驗數(shù)據(jù)，不能獲得顯式表達式，因此為隱函數(shù)。所以，綜合建模方程的待定參數(shù)優(yōu)化問題本質(zhì)是隱式、不連續(xù)函數(shù)的多目標優(yōu)化問題，面向解析函數(shù)的優(yōu)化方法都不適用于本問題。

分層序列法是多目標決策的一種方法，其本質(zhì)是按一個原始目標選擇若干可行解，再按另一目標縮小可行解的范圍，直至所有目標分析完成為止[11]。本文使用分層序列法時，令目標的先后順序依次為F1(總體預(yù)測率)、F2(安全預(yù)測率)、F3(平均相對誤差平方和)，即首先選擇總體預(yù)測率最大的可行解，然后在這些可行解中篩選出安全預(yù)測率最大的可行解，最后在篩選出的可行解中，將平均相對誤差平方和最小的解作為最優(yōu)解。

在綜合建模方程形式中，所有待定系數(shù)及待定指數(shù)統(tǒng)稱為待定參數(shù)，將其作為優(yōu)化變量，形成優(yōu)化數(shù)學(xué)模型。

優(yōu)化變量為

x=(x1,x2,…，x11)=(η1,η2,…，η11)

目標函數(shù)為

(11)

約束條件為

(12)

3.3 待定參數(shù)優(yōu)化確定

對綜合建模方程的待定參數(shù)進行優(yōu)化時，優(yōu)化變量一共有11個。若采用窮舉法優(yōu)化，如文獻[3](優(yōu)化效果非常好)的變量范圍與步長取值做法，則枚舉個數(shù)為[(1.5-0.5)/0.01+1]11=1.1×1022?？梢姼F舉法的枚舉個數(shù)過大，導(dǎo)致運算時間過長且運算所占內(nèi)存過大，因而不能完成優(yōu)化計算。采用差異演化算法時，種群規(guī)模NP取60，演化代數(shù)的最大值N取100，則總枚舉個數(shù)為(N+1)NP=6 060。與窮舉法相比，優(yōu)化效率獲得極大提升。

基于文獻[3]所使用的所有國內(nèi)數(shù)據(jù)，即來自于文獻[12]的31個數(shù)據(jù)和文獻[13]的3個數(shù)據(jù)，采用差異演化算法進行優(yōu)化。

在差異演化過程中，將適應(yīng)度的比較改進為3個目標函數(shù)的比較，比較方式為多目標分層序列法。經(jīng)過種群100代的演化，各個優(yōu)化變量由初始值改變到優(yōu)化值，如表2所示。在演化過程中，最優(yōu)個體的變量絕對值之和隨演化代數(shù)變化曲線如圖2所示，最優(yōu)個體的總體預(yù)測率與安全預(yù)測率隨演化代數(shù)的變化曲線如圖3所示，最優(yōu)個體的平均相對誤差平方和隨演化代數(shù)的變化曲線如圖4所示。

表2 優(yōu)化前后變量的數(shù)值比較

圖2 最優(yōu)個體的變量絕對值之和隨演化代數(shù)變化曲線Fig.2 Changing curve of sum of optimal individual absolute variables with algebra

圖3 最優(yōu)個體的預(yù)測率隨演化代數(shù)變化曲線Fig.3 Changing curves of predicted rate of optimal individual with algebra

由圖2可知，在演化代數(shù)大于80以后，最優(yōu)個體的變量絕對值之和不再發(fā)生變化。由圖3可知，在演化代數(shù)大于30以后，最優(yōu)個體的2個概率型指標不再發(fā)生改變，總體預(yù)測率達到82.35%，安全預(yù)測率達到100%。由圖4可知，在演化代數(shù)大于80以后，最優(yōu)個體的平均相對誤差平方和的變化量在0.001以內(nèi)，其最終值為0.001 3。可見，采用差異演化算法求解11個待定參數(shù)的彈道極限方程建模問題，可以獲得滿意的結(jié)果。同時可以看出，差異演化代數(shù)到達80時，基本上已經(jīng)實現(xiàn)了綜合建模方程所有待定參數(shù)的優(yōu)化確定，表明了此種算法的高效性。

優(yōu)化后的方程是否具有適用性，需要運用其他來源的實驗數(shù)據(jù)進行檢測。本文選用文獻[14]的12個實驗數(shù)據(jù)(4個數(shù)據(jù)為全玄武巖填充層結(jié)構(gòu)，8個數(shù)據(jù)為玄武巖/Kevlar填充層結(jié)構(gòu))和文獻[15]的37個實驗數(shù)據(jù)(全為玄武巖/Kevlar填充層結(jié)構(gòu))。方程預(yù)測結(jié)果為：總體預(yù)測率為83.67%，相對于82.35%提升了1.32%，安全預(yù)測率為95.92%，相對于100%降低了4.08%，平均相對誤差平方和為0.008 6，相對于0.001 3增加了0.007 3。可見，采用差異演化算法獲得的方程具有一定普適性。

進一步細分對這49個實驗數(shù)據(jù)的分類預(yù)測能力，圖5繪出了參數(shù)確定后的綜合建模方程對此實驗數(shù)據(jù)的預(yù)測效果?？梢?，后板失效的數(shù)據(jù)一共有33個，其中有2個數(shù)據(jù)預(yù)測錯誤，則后板失效數(shù)據(jù)的預(yù)測正確率為93.94%；后板未失效的數(shù)據(jù)一共有16個，其中有6個數(shù)據(jù)預(yù)測錯誤，則未失效數(shù)據(jù)的預(yù)測正確率為62.5%。雖然彈道極限方程對未失效數(shù)據(jù)的預(yù)測率較低，但由圖5可見，預(yù)測錯誤的未失效數(shù)據(jù)距離直線1非常近，表明優(yōu)化后的方程誤差非常小。

圖4 最優(yōu)個體的平均相對誤差平方和隨演化代數(shù)變化曲線Fig.4 Changing curve of average sum of squared relative errors of optimal individual with algebra

圖5 優(yōu)化后方程對實驗數(shù)據(jù)的預(yù)測效果Fig.5 Experimental data prediction results using optimized equation

4 結(jié) 論

1) 本文對差異演化算法進行了編程實現(xiàn)，通過對3個包含20個變量的典型數(shù)學(xué)函數(shù)的優(yōu)化驗證，優(yōu)化結(jié)果與理論結(jié)果在0.000 2以內(nèi)，表明本文所編差異演化優(yōu)化程序結(jié)果可信。

2) 采用差異演化算法，基于文獻[3]中使用的國內(nèi)實驗數(shù)據(jù)，對綜合建模方程的11個待定參數(shù)進行3個目標函數(shù)的優(yōu)化確定，獲得了綜合建模方程各個待定參數(shù)的優(yōu)化值。

3) 待定參數(shù)確定后的綜合建模方程，對優(yōu)化過程中所采用的實驗數(shù)據(jù)預(yù)測情況為：總體預(yù)測率為82.35%，安全預(yù)測率為100%，平均相對誤差平方和為0.001 3；對其他來源的49個實驗數(shù)據(jù)預(yù)測情況為：總體預(yù)測率為83.67%，安全預(yù)測率為95.92%，平均相對誤差平方和為0.008 6。經(jīng)比較發(fā)現(xiàn)，概率型指標的偏差在5%以內(nèi)，誤差型指標的偏差在0.01以內(nèi)，因此，所得綜合建模方程的適用性得到驗證。