文献[3]为获得适用于中国航天工程特定填充式防护结构的弹道极限方程，采用穷举法对不同预测指标[4]组合下NASA填充式防护结构弹道极限方程[5]进行双待定参数修正，修正后方程的预测能力更优。由于穷举法具有全空间搜索能力，其结果为所枚举变量的最优值，因此，文献[3]所获得的弹道极限方程的预测能力基本上为最优。但此方法的计算量与待定参数个数呈指数型的关系，因此不适合待定参数较多且数据量大的弹道极限方程建模。
差异演化(Differential Evolution，DE)算法[6]是一种基于群体差异的算法，具有全局搜索能力，且简单易用，收敛速度快[7]。文献[8]通过函数测试表明，差异演化算法在绝大多数情况下比遗传算法等具有更强的全局搜索能力。
除了借鉴NASA填充式防护结构弹道极限方程之外，建立国内填充式防护结构的弹道极限方程具有重要意义。文献[9]通过量纲理论建立了该弹道极限方程形式，确定该方程中包含的3个待定系数和8个待定指数之后，即完成该方程的建模工作。为了求得这些待定参数的数值，本文将此方程中的指数和系数作为优化变量，预测指标(概率型指标和误差型指标)[4]作为目标函数，基于国内数据，采用差异演化算法对方程中的优化变量进行确定。
1 差异演化算法基本原理 差异演化算法的演化过程与遗传算法类似，也包括种群初始化、变异、交叉和选择操作。差异演化是基于实数编码的演化算法，因此可将个体视为实数设计变量组成的向量。差异演化算法首先创建种群，然后对种群的所有个体进行变异操作，接着进行离散交叉生成中间个体，最后中间个体与父个体进行竞争，较优的个体被选中，形成数目不变的下一代新种群。
1.1 种群初始化 在D维空间随机产生个体数为NP的种群，其各个个体的变量为
 | (1) |
式中：xij(0)为初始种群中第i个个体的第j个变量，i为个体在种群中的序号，j为个体中设计变量的序号；NP为种群大小; D为个体中设计变量数；xjL和xjU分别为个体第j个变量的下边界和上边界；rand(0, 1)为在(0, 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)：
 | (2) |
式中：F为缩放因子，取值范围一般为[0, 2]。
1.3 交叉 交叉操作可以增加种群的多样性，变异后的个体hi(t+1)和种群中当前的演化个体Xi(t)通过离散交叉，生成中间个体vi(t+1)：
 | (3) |
式中：rand(i)为[1, D]之间的随机整数；CR为交叉概率，取值范围一般为[0, 1]。
1.4 选择 中间个体vi(t+1)与当前演化个体Xi(t)的适应度进行比较：
 | (4) |
式中：f表示适应度函数。
通过上述变异、交叉和选择，对种群中的每个个体进行循环操作，得到下一代种群，如此演化若干代，获得优化问题的最优解。
2 测试函数演化结果 对基本差异演化算法进行编程实现，以文献[10]中的3个测试函数进行程序验证。
Ackley函数(f1)：
 | (5) |
Rastrigin函数(f2)：
 | (6) |
Griwangk函数(f3)：
 | (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 各函数最优个体的变量绝对值之和与适应度随演化代数变化曲线 Fig. 1 Changing curves of sum optimal individual absolute variables and fitness with generation for each function |
| 图选项 |
由图 1可知，3个函数在种群不断演化过程中，最优个体逐渐到达最优。函数f1和f3在演化代数大于300以后，变量绝对值之和与适应度都基本为零，趋于平稳。函数f2在代数大于1 800以后，变量绝对值之和与适应度都基本为零，趋于平稳。3个函数的优化结果如表 1所示。
表 1 各函数优化结果
| 测试函数 | 理论最优值 | 演化代数 | 优化结果 |
| f1 | 0 | 500 | 0.000 06 |
| f2 | 0 | 2 000 | 0.000 01 |
| f3 | 0 | 500 | 0.000 12 |
表选项
由表 1可以发现，函数f1、f2和f3的优化结果与理论结果误差小于0.000 2，因此本文的差异演化算法程序可得到正确的结果。
3 综合建模方程 3.1 一般形式 针对国内填充式防护结构形式，文献[9]利用量纲理论建立了一种合适的填充式防护结构弹道极限方程形式，称之为综合建模方程，其形式如下所示。
1) 低速区：V≤VL。
 | (8) |
2) 高速区：V≥VH。
 | (9) |
3) 中速区：VL < V < VH。
 | (10) |
上述弹道极限方程共涉及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(平均相对误差平方和)，即首先选择总体预测率最大的可行解，然后在这些可行解中筛选出安全预测率最大的可行解，最后在筛选出的可行解中，将平均相对误差平方和最小的解作为最优解。
在综合建模方程形式中，所有待定系数及待定指数统称为待定参数，将其作为优化变量，形成优化数学模型。
优化变量为
 |
目标函数为
 | (11) |
约束条件为
 | (12) |
式中：PTotal、PSafe、Sa2分别为总体预测率、安全预测率、平均相对误差平方和。
3.3 待定参数优化确定 对综合建模方程的待定参数进行优化时，优化变量一共有11个。若采用穷举法优化，如文献[3](优化效果非常好)的变量范围与步长取值做法，则枚举个数为[(1.5-0.5)/0.01+1]11=1.1×1022。可见穷举法的枚举个数过大，导致运算时间过长且运算所占内存过大，因而不能完成优化计算。采用差异演化算法时，种群规模NP取60，演化代数的最大值N取100，则总枚举个数为(N+1)NP=6 060。与穷举法相比，优化效率获得极大提升。
基于文献[3]所使用的所有国内数据，即来自于文献[12]的31个数据和文献[13]的3个数据，采用差异演化算法进行优化。
在差异演化过程中，将适应度的比较改进为3个目标函数的比较，比较方式为多目标分层序列法。经过种群100代的演化，各个优化变量由初始值改变到优化值，如表 2所示。在演化过程中，最优个体的变量绝对值之和随演化代数变化曲线如图 2所示，最优个体的总体预测率与安全预测率随演化代数的变化曲线如图 3所示，最优个体的平均相对误差平方和随演化代数的变化曲线如图 4所示。
表 2 优化前后变量的数值比较 Table 2 Comparison of variable values before and after optimization
| 变量 | 初始值 | 优化值 |
| x1 | 0.6 | 0.457 |
| x2 | 0.3 | 0.4 |
| x3 | 0.2 | 0.1 |
| x4 | -0.3 | -0.311 |
| x5 | 1 | 1.05 |
| x6 | 0.3 | 0.25 |
| x7 | 0.4 | 0.284 |
| x8 | 1.2 | 1.5 |
| x9 | 0.3 | 0.21 |
| x10 | 0.2 | 0.11 |
| x11 | -0.25 | -0.339 |
表选项
 |
| 图 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 |
| 图选项 |
 |
| 图 4 最优个体的平均相对误差平方和随演化代数变化曲线 Fig. 4 Changing curve of average sum of squared relative errors 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非常近，表明优化后的方程误差非常小。
 |
| 图 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以内，因此，所得综合建模方程的适用性得到验证。
参考文献
| [1] | 闫军, 郑世贵. 填充式防护结构填充层撞击特性研究[J].载人航天, 2013, 19(1): 10–14. YAN J, ZHENG S G. Study on impact characteristics of stuffed shielding structure[J].Manned Spaceflight, 2013, 19(1): 10–14.(in Chinese) |
| [2] | CHRISTIANSEN E, LAMBERT M, STOKES H. IADC protection manual[M].Bremen: Inter Agency Debris Committee, 2002. |
| [3] | 贾光辉, 欧阳智江, 蒋辉, 等. 填充式防护结构弹道极限方程的多指标寻优[J].北京航空航天大学学报, 2013, 39(12): 1573–1583. JIA G H, OUYANG Z J, JIANG H, et al. Multiple indicators optimization for stuffed whipple shield ballistic limit equation[J].Journal of Beijing University of Aeronautics and Astronautics, 2013, 39(12): 1573–1583.(in Chinese) |
| [4] | 贾光辉, 欧阳智江, 蒋辉. 撞击极限方程预测指标剖析与实例[J].航空学报, 2013, 34(10): 2364–2371. JIA G H, OUYANG Z J, JIANG H. Analysis and instances of ballistic limit equations' predictive indicators[J].Acta Aeronautica et Astronautica Sinica, 2013, 34(10): 2364–2371.(in Chinese) |
| [5] | CHRISTIANSEN E L, KERR J H. Ballistic limit equations for spacecraft shielding[J].International Journal of Impact Engineering, 2001, 26(1): 93–104. |
| [6] | STORN R, PRICE K. Differential evolution-A simple and efficient adaptive scheme for global optimization over continuous spaces: TR-95-012[R]. Berkeley: University of California, 1995. |
| [7] | 陈良, 戴光明, 张全元, 等. 差异演化及其改进形式的综述[J].计算机工程与设计, 2008, 29(1): 131–134. CHEN L, DAI G M, ZHANG Q Y, et al. Overview of differential evolution algorithm and its improved algorithms[J].Computer Engineering and Design, 2008, 29(1): 131–134.(in Chinese) |
| [8] | STORN R, PRICE K. Differential evolution-A simple and efficient heuristic for global optimization over continuous spaces[J].Journal of Global Optimization, 1997, 11(4): 341–359.DOI:10.1023/A:1008202821328 |
| [9] | 姚光乐, 贾光辉. 填充式防护结构弹道极限方程形式建模[J].空间碎片研究, 2017, 17(1): 29–33. YAO G L, JIA G H. Formal modeling of ballistic limit equations for stuffed whipple shield[J].Space Debris Research, 2017, 17(1): 29–33.(in Chinese) |
| [10] | 袁俊刚, 孙治国, 曲广吉. 差异演化算法的数值模拟研究[J].系统仿真学报, 2007, 19(20): 4646–4648. YUAN J G, SUN Z G, QU G J. Simulation study of differential evolution[J].Journal of System Simulation, 2007, 19(20): 4646–4648.DOI:10.3969/j.issn.1004-731X.2007.20.011(in Chinese) |
| [11] | 汪勇, 徐琼, 张凌, 等. 基于遗传分层序列法的云制造资源优化配置[J].统计与决策, 2016, 20: 80–83. WANG Y, XU Q, ZHANG L, et al. Optimal allocation of cloud manufacturing resources based on genetic stratified sequence method[J].Statistics and Decision, 2016, 20: 80–83.(in Chinese) |
| [12] | 张宝玺. 超高速撞击玄武岩及Kevlar纤维布填充防护结构优化设计[D]. 哈尔滨: 哈尔滨工业大学, 2011. ZHANG B X. The optimal structural design of stuffed shields with basalt and Kevlar fiber clothes on hypervelocity impacting[D]. Harbin: Harbin Institute of Technology, 2011(in Chinese).http://cdmd.cnki.com.cn/Article/CDMD-10213-1012002274.htm |
| [13] | 哈跃. 玄武岩纤维材料及其填充防护结构超高速撞击特性研究[D]. 哈尔滨: 哈尔滨工业大学, 2009. HA Y. Research on hypervelocity impact properties of woven of basalt fiber and its stuffed shielding structure[D]. Harbin: Harbin Institute of Technology, 2009(in Chinese).http://cdmd.cnki.com.cn/article/cdmd-10213-2011015844.htm |
| [14] | 祖士明. 玄武岩及Kevlar纤维填充式防护结构超高速撞击性能研究[D]. 哈尔滨: 哈尔滨工业大学, 2013. ZU S M. Research on hypervelocity impact performance of shielding structure stuffed with basalt and Kevlar fiber[D]. Harbin: Harbin Institute of Technology, 2013(in Chinese).http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=D416804 |
| [15] | 贾古寨, 哈跃, 庞宝君, 等. 玄武岩/Kevlar纤维布填充防护结构撞击极限及损伤特性[J].爆炸与冲击, 2016, 36(4): 433–440. JIA G Z, HA Y, PANG B J, et al. Ballistic limit and damage properties of basalt/Kevlar stuffed shield[J].Explosion and Shock Waves, 2016, 36(4): 433–440.DOI:10.11883/1001-1455(2016)04-0433-08(in Chinese) |
