填充式防護結構彈道極限方程的差異演化優化

賈光輝， 姚光樂， 張帥

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

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

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

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

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

1 差異演化算法基本原理

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

1.1 種群初始化

在D維空間隨機產生個體數為NP的種群，其各個個體的變量為

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

(1)

1.2 變 異

變異是差異演化算法的關鍵步驟，基于個體向量差進行。設當前種群中的個體為Xi(t)，t為演化代數。從當前種群中隨機選擇3個個體Xr1(t)、Xr2(t)和Xr3(t)(r1≠r2≠r3≠i),取后2個個體向量的差(Xr2(t)-Xr3(t))，使用縮放因子F調節大小，然后與第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]之間的隨機整數；CR為交叉概率，取值范圍一般為[0,1]。

1.4 選 擇

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

Xi(t+1)=

(4)

式中：f表示適應度函數。

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

2 測試函數演化結果

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

Ackley函數(f1)：

xi∈[-32,32]

(5)

Rastrigin函數(f2)：

xi∈[-5.12,5.12]

(6)

Griwangk函數(f3)：

xi∈[-512,512]

(7)

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

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

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

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

由表1可以發現，函數f1、f2和f3的優化結果與理論結果誤差小于0.000 2，因此本文的差異演化算法程序可得到正確的結果。

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

測試函數理論最優值演化代數優化結果f105000.00006f2020000.00001f305000.00012

3 綜合建模方程

3.1 一般形式

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

1) 低速區：V≤VL。

(8)

2) 高速區：V≥VH。

(9)

3) 中速區：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為低速區與中速區的分界速度,km/s；VH為高速區與中速區的分界速度，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。此方程的待定參數共11個，包含3個待定系數和8個待定指數，待定系數為ηi(i=1,5,8)，待定指數為ηj(j=2,3,4,6,7,9,10,11)，均無量綱。

3.2 待定參數優化模型

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

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

在綜合建模方程形式中，所有待定系數及待定指數統稱為待定參數，將其作為優化變量，形成優化數學模型。

優化變量為

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

目標函數為

(11)

約束條件為

(12)

3.3 待定參數優化確定

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

基于文獻[3]所使用的所有國內數據，即來自于文獻[12]的31個數據和文獻[13]的3個數據，采用差異演化算法進行優化。

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

表2 優化前后變量的數值比較

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

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

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

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

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

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

圖5 優化后方程對實驗數據的預測效果Fig.5 Experimental data prediction results using optimized equation

4 結 論

1) 本文對差異演化算法進行了編程實現，通過對3個包含20個變量的典型數學函數的優化驗證，優化結果與理論結果在0.000 2以內，表明本文所編差異演化優化程序結果可信。

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

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