CN103440377B - 基于改进并行de算法的飞行器气动外形优化设计方法 - Google Patents

Abstract

本发明是一种基于改进并行DE算法的飞行器气动外形优化设计方法，用于飞行器气动外形优化设计。本发明在模式池中集成不同的DE算法模式，扩展了DE模式选择范围，还新增一种动态参数的模式，该模式可对大范围的未知解空间进行更为有效地搜索。在迭代过程中，经过变异操作和交叉操作产生测试向量；把种群中个体编码还原为对应飞行器外形，将CFD计算分析过程引入个体适应度评价阶段，作为对新一代个体进行选择的依据，使得设计结果分析过程与参数修改过程有效结合，缩短了设计周期；经选择操作产生子代群体。鉴于飞行器气动外形设计过程计算量大的特点，本发明引入了并行计算机制，从而大大降低算法运行时间，提高设计效率。

Description

基于改进并行DE算法的飞行器气动外形优化设计方法

技术领域

本发明涉及飞行器气动外形优化设计技术领域，以及计算机仿真与优化技术领域，具体涉及一种基于改进的并行DE（Differential Evolution，差分进化）算法进行飞行器气动外形优化设计的方法。

背景技术

设计性能优良的气动外形对于整个飞行器设计的重要性是不言而喻的。

传统的飞行器气动外形设计方法以风洞试验或者CFD分析结果为基础，多以依赖经验的“Cut and Try试凑法”为主。当今要求的飞行器性能指标越来越高、指标的越来越多，对这些传统气动设计方式和手段提出了严峻挑战。首先，传统设计方法主要根据经验或大量的风动试验，通过人工试凑与反复修正完成设计。即使引入CFD分析与计算，也多用来对设计方案进行评估，以减少风洞试验量，为人工修正设计方法提供更多参考。其次，人工试凑与修正无法定量掌握其对整个流场及其气动力特性的影响，使得传统设计方式和设计手段难以实现多学科、多目标、多约束气动外形设计的要求。

随着计算机技术以及计算流体力学的发展，在飞行器气动外形设计过程中，越来越多地采用优化设计方法以获得性能更好的飞行器外形。一方面，采用气动外形优化设计可以比采用传统方法获得更好的气动外形，提高飞行器设计质量的同时缩短设计周期，提高效率，从而节约成本；另一方面，随着对飞行器性能要求的增加，飞行器设计过程中所面临的矛盾越来越多，难度也越来越大，利用外形优化设计方法，有助于在这些复杂约束条件下得到良好的设计结果。也正因为如此，气动外形优化设计方法受到了高度重视，新方法也层出不穷。概括地说，这些方法可以分成两类，即基于梯度信息的优化设计方法和基于进化算法的优化设计方法。前者需要准确地确定并有效地利用目标函数随设计变量变化的梯度信息，显然，这在有时候是很不方便的。而基于进化算法的优化设计方法，相对而言求解过程与具体问题的耦合性较弱，因而通用性较强，且设计具有良好的全局优化特性，因而该类方法的应用越来越广泛。

进化算法有很多分支，每个分支又有诸多变形。差分进化(Differential Evolution，DE)算法就是其中的一种基于群体进化的算法。它能够通过种群内个体间的合作与竞争来实现对优化问题的求解，具有记忆个体最优解和种群内信息共享的特点，其本质是一种基于实数编码的具有保优思想的贪婪遗传算法。对于简单的变异操作和一对一的竞争生存策略，差分进化算法用实数编码来实现，这使得遗传操作的复杂性大大降低。和其他进化类算法相比，差分进化算法虽然保留了基于种群的全局搜索策略，但是其特有的记忆能力，使其不需要借助问题的特征信息就可以动态跟踪当前的搜索情况，及时地对搜索策略做出调整。这也是差分进化算法具有较强的全局收敛能力和较高的鲁棒性的主要原因之一。差分进化算法在被提出的第二年就被证明是最快的进化算法。同时它还具有内含并行性，算法采用种群的方式组织搜索，因而可以同时搜索解空间的多个区域，并相互交流信息，这种搜索方式使得它虽然每次只执行与种群规模NP成比例的计算，而实质上己经进行了大约O(N²)次有效搜索，能以较少的计算获得较大的利益。由于其简单易用、稳健性好以及强大的全局搜索能力，使得差分进化算法已在多个领域取得成功。

发明内容

本发明将一般DE算法进行改进，使其与CFD（Computational Fluid Dynamics，计算流体动力学）技术结合，形成了一种新型的飞行器气动外形优化设计方法。

本发明提供的基于改进并行DE算法的飞行器气动外形优化设计方法，引入了模式池的概念，具体包括如下步骤：

步骤一：根据飞行器气动外形的性能指标和设计要求，设定优化变量及变量对应范围；

步骤二：设定差分进化的模式；

步骤三：产生初始种群；

步骤四：对种群个体进行适应度评估；

步骤五：根据选定的模式对种群进行变异操作；

步骤六：根据选定的模式对种群进行交叉操作；

步骤七：对生成的子代个体采用并行化操作进行适应度评估；

步骤八：对子代和父代进行选择操作，将两代中对应的个体进行比较，选取适应度较好的个体作为新一代种群的成员；

步骤九：根据预先设定的终止条件和当前状态决定是否终止进化过程，若终止条件满足，退出进化，转到步骤十；否则，转向步骤五；

步骤十：输出并显示最优结果。

所述的步骤四和步骤七中，进行适应度评估的方法是：首先，获得种群个体对应的优化变量的实际数值，由优化变量的数值通过飞行器几何建模技术获取对应的飞行器几何外形；然后，将获取的飞行器几何外形经CFD分析程序的网格划分、仿真计算后获得量化的性能指标，所获得的性能指标越好，则表示所对应的优化变量的适应度越好。

本发明的飞行器气动外形优化设计方法具有如下优点和积极效果：

（1）本发明引入了模式池（Pattern Pool）的概念，应用时可根据优化维度、收敛性及多样性的需求，选择其中一种模式，也可以在利用各种模式进行计算，对各种方法对比后选择最佳的结果；

（2）在模式池中除集成已有经典有效的DE算法模式外，也发明并集成了一种参数服从Lévy分布的动态参数新模式。与其他参数为常数的模式相比，该模式在搜索未知的大规模区域时，搜索效率更高，经测试取得了不错的效果，尤其适用于高维情形（例如，维数大于100），扩大了应用的选择余地。

（3）针对实际中基于CFD进行气动优化设计过程中计算量大的情况，本发明引入了并行计算机制，从而有效地节省算法运行时间，提高求解效率；

（4）本发明有效集成智能优化算法和CFD分析计算过程，将CFD计算作为DE算法中对种群个体适应度评估的手段，使得设计结果分析过程与参数修改过程有效结合。CFD分析过程不再是简单的对设计结果进行简单评估的工具，而且能够参与到循环迭代求解设计过程中，有效地指导迭代过程。从而缩短设计周期的同时，减少人工量，提高设计质量。

附图说明

图1是本发明的飞行器气动外形优化设计方法的主要步骤流程图；

图2是对比了新模式中的动态参数服从的Lévy分布和布朗运动；图(a)为布朗运动，图(b)为二维下的Lévy flight运动；

图3为种群个体进行适应度评估并行计算过程的伪代码及过程相关图示。

具体实施方式

下面将结合附图对本发明方法作进一步的详细说明。

本发明所基于的DE算法是由Rainer Storn和Kenneth Price为求解切比雪夫多项式于1996年共同提出的一种采用浮点矢量编码在连续空间中进行随机搜索的优化算法。该DE算法的原理简单，受控参数少，实施随机、直接的全局搜索，易于理解和实现。

差分进化的基本步骤包括变异、交叉和选择。在基本模式下，随机选择两个不同的个体矢量相减生成差分矢量，将差分矢量赋予权值之后加到第三个随机选择的个体矢量上，生成变异矢量，该操作称为变异。将变异矢量与目标矢量以一定的方式进行参数混合，生成试验矢量，这一过程称之为交叉。如果试验矢量的适应度优于目标矢量的适应度，则用试验矢量取代目标矢量而形成下一代，该操作称为选择。本发明方法中不同模式的差异在于进行变异操作采取不同的方式。下面具体阐述本发明方法各步骤的具体实现方式。

步骤一：设定优化变量及变量范围；

飞行器的气动阻力、升阻比等性能指标在很大程度上依赖于在飞行器设计过程中所选择的核心几何参数。要根据设计要求综合考虑气动力/气动热、飞行轨迹、运载条件、雷达散射截面等各方面，选择合适的设计几何变量及对应的范围。一般而言，各个性能指标对参数要求是矛盾的，这就需要对其进行折中。一般选择主要指标作为目标函数，将其余的设计指标作为约束进行处理。

步骤二：选择合适的模式；

本发明提供多种模式可供选择，也可以向该模式池中添加新模式以对其进行扩展。各个模式的求解性能因问题规模略有差异，可以根据优化变量的维数大小，在尝试后根据经验和喜好进行选择。

模式以DE/a/b/c的形式给出，其中DE是差分进化（Differential Evolution）的缩写，a指当前被变异的向量是种群中“随机的”（rand）还是“最佳的”（best）个体，“随机的”能够增加多样性，“最佳的”则能够促进收敛；b是所利用的差向量的个数，取值一般限于1和2；c指示交叉程序的操作方法，一般采用二项交叉操作（binomial crossover），简写为bin，也有采用指数型交叉操作的，简写为exp。基本DE算法可以表示为DE/rand/1/bin。图1中还列出了几种其它模式，例如DE/rand/1/exp，DE/best/2/bin等。

步骤三：产生初始种群；

进行本发明方法需要设置合理的初始种群，种群中的个体代表满足约束的求解值。所谓“合理”，就是不违反设定的性能约束条件。在这些个体的基础上，通过不断地迭代进化，最终得到比较好的求解结果。最终求解结果的优劣，在一定程度上依赖于初始值。初始种群比较好，进化过程一般会比较快、顺利，结果也会接近最优解；否则，有可能进化过程会陷入局部最优值而停滞。

在无特殊要求的情况下，可假设初始种群个体的相关参数在可行范围内随机产生。

种群所含种群规模NP是算法的一个重要参数，一般介于5D与10D之间，D为问题空间的维度，但不能少于4，否则无法进行变异操作。NP越大，种群多样性越强，获得最优解概率越大，但是计算时间更长。折中考虑，一般种群规模NP可在20-50取值。

步骤四：对种群进行评估；

对种群个体进行适应性评估可以确定在选定的标准下个体相对于整体的性能优劣。常规评估依据的是个体对应的适应度函数的值，本发明方法没有对应的简单函数可用，采用的是一套CFD计算分析方法，见步骤七。

步骤五：根据选定的模式对种群进行变异操作。

差分进化的变异操作，概括来讲，就是将种群中两个成员之间的加权差向量加到第三个成员上来，以此产生变异向量。由于选择差向量的个数不同、被加向量的选择方式不同、权重不同，产生变异向量的方式千变万化，这也是差分算法拥有众多变种的原因。基本的差分算法利用随机选择的两个向量产生差分向量，将其叠加到随机选择的第三个向量上，即

v_i,G+1＝x_r3,G+F(x_r1,G-x_r2,G)

其中，v_i,G+1表示第G+1代的第i个个体（i=1，2，……，NP，NP为种群规模），x_r1,G,x_r2,G,x_r3,G表示从当前代中随机选择的三个不同个体，F为变异因子。

变异因子F是控制种群多样性和收敛性的重要参数，一般在[0,2]之间取值。大部分模式下，参数F都保持为常数。变异因子F值较小时，种群的差异度减小，进化过程不易跳出局部极值导致种群过早收敛。变异因子F值较大时，虽然容易跳出局部极值，但是收敛速度会减慢，一般可选F的取值范围为0.3～0.6。

其它的模式，如：DE/best/1/bin，对应的变异向量v_i,G+1＝x_best,G+F(x_r1,G-x_r2,G)；DE/rand/2/bin，对应的变异向量v_i,G+1＝x_r5,G+F(x_r1,G-x_r2,G+x_r3,G-x_r4,G)；DE/best/2/bin，对应的变异向量v_i,G+1＝x_best,G+F(x_r1,G-x_r2,G+x_r3,G-x_r4,G)；DE/rand-to-best/1/bin，对应的变异向量v_i,G+1＝βx_best,G+(1-β)x_r1,G+F(x_r2,G-x_r3,G)，β∈(0,1)等等。其中，x_best,G表示当前第G代的最优个体，x_r1,G、x_r2,G、x_r3,G、x_r4,G和x_r5,G均是从当前G代随机选择的五个不同个体。

除了这些常规变异模式，本发明模式池中也提供了一种动态参数新模式。该模式的变异策略基于Lévy flights。Lévy flights是一种带有长程跳跃的随机游走过程，游走的步长服从Lévy分布。

Lévy分布的概率密度函数L_α,γ(z)为：

$L_{\mathrm{\α},\mathrm{\γ}}\left(z\right)=\frac{1}{\mathrm{\π}}{\∫}_0^{\∞}\exp \left({-\mathrm{\γq}}^{\mathrm{\α}}\right)\cos \left(\mathrm{qz}\right)\mathrm{dq}$

式中γ＞0是范围参数，0＜α＜2用于控制分布形状。

为应对CFD优化面临的大规模搜索空间，本发明引入的新模式将Lévy分布应用到变异因子上。具体地说，第G+1代的变异向量v_i,G+1由以下公式生成

v_i,G+1＝x_i,G+F(x_best,G-x_i,G)+l(x_r1,G-x_r2,G)

即在当代对应个体的基础上叠加两个差分向量，第一个差分向量(x_best,G-x_i,G)由第G代的最优个体和当前个体生成，第二个差分向量(x_r1,G-x_r2,G)由第G代随机选择的两个不同个体生成。通常使用第一个差分向量有助于增加收敛速度，它对应的变异因子F设定为常数；第二个差分向量用来增加多样性，它对应的变异因子l由Lévy分布产生。

实际应用产生Lévy分布服从的随机数时，常采用Mantegna提出的近似方法，即变异因子l由下式产生：

$l=\frac{x}{{\left|y\right|}^{1/\mathrm{\α}}}$

式中，x和y相互独立且都服从正态分布， 一般取

${\mathrm{\σ}}_x=\left[\frac{\mathrm{\Γ}(1+\mathrm{\α})\sin (\mathrm{\π\α}/2)}{\mathrm{\Γ}((1+\mathrm{\α})/2){\mathrm{\α}2}^{(\mathrm{\α}-1)/2}}\right],$ σ_y＝1

其中，Γ(.)为Gamma函数。实验证实，参数α＝1.5，F＝0.3时能取得不错的效果。

图2中(a)的布朗运动和(b)中的Lévy flight运动，运动均从原点开始，运动步数为1000，运动方向满足均匀分布。布朗运动的步长满足标准正态分布，而Lévy flight运动的步长满足Lévy分布，参数α＝1.5，γ=1。由图可以看出，Lévy flight运动类似于布朗运动。但是布朗运动的步长分布属于正态分布，步长取很大值的概率较低；而Lévy分布属于重尾分布，重尾分布意味着可以以更大的概率获得很大的值。因此，Lévy flights在未知的大规模区域内进行更加有效地搜索。具体将其应用到智能优化算法中，Lévy flights能够扩大搜索范围、增加种群多样性，避免陷入局部最优点。

步骤六：根据模式池提供的模式对种群进行交叉操作；

进行交叉操作是为了更加随机化进化过程，增添多样性。交叉操作通过交叉因子CR控制变异向量v_i,G+1和种群个体x_i,G的基因进行交换，进而生成试验向量u_i,G+1。假设种群中的个体维数为D，则试验向量的第j个分量（j=1,2,…,D）表达式为：

$u_{\mathrm{ji},G+1}=\left\{\begin{array}{c}{v}_{\mathrm{ji},G+1}\mathrm{if}(\mathrm{rand}\left(j\right)\≤\mathrm{CR})\mathrm{or\; j}=\mathrm{rnbr}\left(i\right)\\ x_{\mathrm{ji},G}\mathrm{if}(\mathrm{rand}\left(j\right)>\mathrm{CR})\mathrm{and\; j}\≠\mathrm{rnbr}\left(i\right)\end{array}\right.$

(i＝1,2,…,NP,j＝1,2,…D)

式中，rand(j)为取值在0～1的随机数；rnbr(i)从序列1,2,……,D中随机取值，用以保证每个u_i,G+1至少有一个分量从v_i,G+1获取；CR为交叉因子。

交叉因子CR可控制个体参数的各维对交叉的参与程度，以及全局与局部搜索能力的平衡，一般在[0,1]之间。交叉因子CR越小，种群多样性减小，容易过早收敛。CR越大，收敛速度越大。但过大可能导致收敛变慢，因为扰动大于了种群差异度。一般CR的取值选在[0.6～0.9]之间。

步骤七：对生成的子代个体采用并行化操作进行适应度评估，评估子步骤包括：生成参数化外形；生成网格；CFD计算；适应度计算。

首先获得种群各个个体对应的设计变量的实际数值，然后由每组设计参数通过诸如NURBS曲线拟合的CAD方法、二次曲线建模或者其他飞行器几何建模技术获得对应的飞行器几何外形，作为CFD分析评估性能模块的输入。然后经CFD分析程序网格划分、仿真计算后获得量化的性能指标，所获得的性能指标越好，则表示所对应的优化变量的适应度越好，这样就完成了种群个体的适应度评估。必要的时候可以将最优个体对应的分析结果进行可视化显示。

鉴于CFD分析计算过程的计算量比较大，况且现在每次迭代要对多个个体进行CFD分析，故一般在此阶段添加并行计算策略，节省时间。

对各个种群个体的适应度评估过程相互独立，Matalb中支持并行计算的Parfor（ParallelFOR-loop，并行for循环）结构是很好的选择。进入评估阶段，使用Parfor结构将不同个体的评估过程分配到多个Matlab工作节点（Worker）上运行。Matlab工作节点作为一个逻辑计算单元供parfor结构进行任务调配，parfor结构检测Matlab工作节点的存在状态，并且自动实现Matlab客户端与工作节点之间的数据交互。而对每个个体的评估均包括生成参数化外形、划分网格、CFD计算和适应度计算的过程。并行计算的伪代码及过程图示见图3。

图3左侧显示的是一个利用Matlab中并行计算工具箱的parfor结构进行并行处理的伪代码例子。首先，为了能进行并行执行求解过程，需要利用matlabpool open命令打开一个Matlab工作节点池；相应地，结束时须利用matlabpool close命令将节点池关闭。循环的主体除了将常规for循环中的for关键字替换成parfor，其余部分没有任何差别，也即随机选择若干种群个体，进行变异操作得到变异向量v_i,G+1（对应图中字符），然后进行交叉操作得到试验向量u_i,G+1（对应图中字符）；通过生成参数化外形、网格划分、CFD计算等过程计算试验向量u_i,G+1的适应度f(u_i,G+1)，将其与原种群对应个体的适应度f(x_i,G)进行比较，选择适应度较好的个体作为新种群的对应个体x_i,G+1（图3中设定的适应度值越小，个体越优）。图3右侧是并行处理过程的形象化描述。每一次种群更新过程进行到变异、交叉和适应度计算等最耗时的部分时，Matlab客户端通过parfor结构将其分配给多个工作节点执行；各节点执行后的结果会自动进行汇总得到最终计算结果。

步骤八：对子代和父代进行选择操作。

到该阶段，本发明有两组一一对应种群个体：当代种群个体（x_1,G，x_2,G，…x_NP,G）和试验向量个体（u_1,G+1,u_2,G+1,…,u_NP,G+1）。根据他们的适应度，决定哪个个体会被选择，成为下一代的成员。具体而言，就是将试验向量种群的每个个体u_i,G+1分别和当代种群中对应的个体x_i,G进行比较，选取适应度较高的那一个作为新一代的成员。

步骤九：判断是否终止优化。

根据预先设定的终止条件和当前状态决定是否终止进化过程，若终止条件满足，退出进化，转到步骤十；否则，转向步骤五。

最大迭代代数t_max一般作为进化过程的终止条件。迭代次数越大，最优解更精确，但同时计算的时间会更长，需要根据具体问题设定。

步骤十：输出并显示最优结果。

从最终种群中选出适应度最高的个体作为结果，并将其进行可视化显示，在其基础上决定后续设计过程。

Claims (4)

1.一种基于改进并行DE算法的飞行器气动外形优化设计方法，其特征在于，该方法包括如下步骤：

步骤一：根据飞行器气动外形的性能指标和设计要求，设定优化变量及变量对应范围；

步骤二：在模式池中选取差分进化的模式；

步骤三：产生初始种群；

步骤四：对种群个体进行适应度评估；

步骤五：根据选定的模式对种群进行变异操作；

步骤六：根据选定的模式对种群进行交叉操作；

步骤七：对生成的子代个体采用并行化操作进行适应度评估；

步骤八：对子代和父代进行选择操作，将两代中对应的个体进行比较，选取适应度较好的个体作为新一代种群的成员；

步骤九：根据预先设定的终止条件和当前状态决定是否终止进化过程，若终止条件满足，退出进化，转到步骤十；否则，转向步骤五；

步骤十：输出并显示最优结果；

所述的步骤四和步骤七中，进行适应度评估的方法是：首先，获得种群个体对应的优化变量的实际数值，由优化变量的数值通过飞行器几何建模技术获取对应的飞行器几何外形；然后，将获取的飞行器几何外形经CFD分析程序的网格划分、仿真计算后获得量化的性能指标，所获得的性能指标越好，则表示所对应的优化变量的适应度越好，CFD表示计算流体动力学；

步骤二中所述的模式池中集成了一种动态参数新模式，在该模式下，第G+1代的变异向量v_i,G+1由以下公式生成：

v_i,G+1＝x_i,G+F(x_best,G-x_i,G)+l(x_r1,G-x_r2,G)

其中，x_i,G表示第G代的第i个个体，x_best,G表示第G代的最优个体，x_r1,G,x_r2,G表示从当前G代中随机选择的两个不同个体，第一个变异因子F为常量，第二个变异因子l服从Lévy分布；

所述的第二个变异因子l，实际实现时由Mantegna提出的近似算法产生，表示为：变量x和y相互独立且都服从正态分布，

其中， ${\mathrm{\σ}}_x=\[\frac{\mathrm{\Γ}(1+\mathrm{\α})sin(\mathrm{\π}\mathrm{\α}/2)}{\mathrm{\Γ}\left(\right(1+\mathrm{\α})/2)\mathrm{\α}{2}^{(\mathrm{\α}-1)/2}}\],$ σ_y＝1，0＜α＜2用于控制分布形状。

2.根据权利要求1所述的飞行器气动外形优化设计方法，其特征在于，所述的第一个变异因子F设置为0.3，Lévy分布中控制分布形状的参数α设置为1.5。

3.根据权利要求1所述的飞行器气动外形优化设计方法，其特征在于，所述的步骤六中，交叉操作通过交叉因子CR控制变异向量v_i,G+1和种群个体x_i,G的基因进行交换，进而生成试验向量u_i,G+1；设种群中的个体维数为D，则试验向量的第j个分量u_ji,G+1为：

$u_{ji,G+1}=\left\{\begin{array}{cc}{v}_{ji,G+1}& if\left(rand\left(j\right)\≤CR\right)orj=rnbr\left(i\right)\\ x_{ji,G}if& if\left(rand\left(j\right)>CR\right)andj\≠rnbr\left(i\right)\end{array}\right.$

(i＝1,2,…,NP,j＝1,2，…D)

其中，rand(j)为取值在0～1的随机数；rnbr(i)为从序列1,2,……,D中的随机取值；交叉因子CR取值在[0.6～0.9]之间，NP表示种群规模。

4.根据权利要求1所述的飞行器气动外形优化设计方法，其特征在于，所述的步骤七中，所述的并行化操作，具体是：利用Matlab中支持并行计算的Parfor结构，将不同种群个体分配到两个以上的Matlab工作节点上进行适应度评估；每个Matlab工作节点作为一个逻辑计算单元供parfor结构进行任务调配，parfor结构检测Matlab工作节点的存在状态，并且自动实现Matlab客户端与工作节点之间的数据交互。