基于并行交换的增强粒子群优化算法在气动优化中的应用

王培君, 夏露, 周文硕, 栾伟达

(西北工业大学 航空学院, 陕西 西安 710072)

随着计算流体力学(CFD)技术的不断发展,气动设计中CFD的占比逐渐提升,数值优化方法已成为飞行器气动外形设计中的一个重要组成部分。在各类优化方法中,启发式算法表现出了十分强大的寻优能力。目前较常用的有以粒子群优化(particle swarm optimization,PSO)算法[1]、蚁群算法(ant clony optimization,ACO)、布谷鸟搜索算法(cuckoo search algorithm,CSA)[2]为代表的群集智能算法,和以差分进化算法(differential evolution,DE)[3]、遗传算法(genetic algorithms,GA)为代表的进化算法。其中PSO算法简单方便、收敛速度较快,对许多问题都有较好的优化效果,但是其更新公式会导致粒子被吸引在当前全局最优解附近,搜索空间被大幅压缩,分散度下降从而导致算法陷入局部最优[4]。CSA算法是一种全局寻优能力较强的算法,较强的搜索随机性极大提升了该算法的全局搜索能力,但是CSA也存在收敛速度慢的问题。DE算法主要通过差分变异操作来生成新的变异个体,再通过交叉和选择得到新的种群个体,贪婪选择策略使其有较强的局部搜索能力,但也容易陷入局部最优。由此可见,有的算法不擅长全局搜索,容易陷入局部最优;有的算法不擅长局部搜索,搜索效率低,在单峰问题中表现差。所以如何补足短板,使算法适应更多形式问题,增强算法的鲁棒性是一个难题。

为了提高启发式算法的鲁棒性,已经有了很多改进措施,比如借鉴其他算法优秀思想对算法更新公式进行改进或在一个种群内通过一定策略使用多种更新公式;对种群拓扑结构进行改进[5],改善种群多样性;还有通过多个种群并行,再利用某种信息交换机制来提高全局搜索能力。其中对算法更新公式进行改进和单一种群内使用多个更新公式的方法改善了原算法搜索模式,但是其搜索模式仍有很大相似度,原算法的短板无法得到有效提升。种群拓扑方法提高了算法全局性能,但是由于算法更新公式并未改变,所以鲁棒性提高不显著。而多种群并行中的分层交换[6]的方法将2种不同算法有机结合起来,共享最优信息,互相引导,取长补短,可以较好改进全局搜索能力。李丁等[6]通过分层交换操作将PSO和DE这2种容易陷入局部最优的算法结合起来得到了RADPSO算法,新算法鲁棒性得到了较大提升,体现了分层交换操作的有效性。但是限于其使用的PSO和DE算法在很多情况下无法得到较优解,所以该算法在有些多峰问题中仍会陷入局部最优,其鲁棒性仍有较大提升空间。

为了提高PSO算法鲁棒性,本文选择将PSO种群与其他算法种群进行并行,然后进行分层交互操作,其中选择什么算法种群进行并行至关重要。科研人员通常借鉴全局性较好的算法对粒子群优化算法进行改进。石季英等[7]构建了PSO-CS混合算法,其核心是将PSO算法每轮迭代后个体历史最优位置进行CSA算法位置更新,新算法全局搜索能力虽提升有限,但是说明了通过CSA算法来提升PSO算法全局寻优能力是可行的。鉴于此,如果通过使用分层交换操作来引入CSA算法,那就可以通过改善算法全局寻优能力来达到高鲁棒性。故本文提出了基于并行交换的增强粒子群优化算法(enhanced particle swarm optimization based on parallel excha-nge,EPSOBPE)。该算法首先借鉴差分进化的变异行为对粒子群优化算法进行增强,在不减少全局收敛性的基础上加强快速收敛能力。然后通过让增强后的粒子群与布谷鸟种群进行并行运行,每次迭代后进行分层交换,使二者种群可以有机结合起来,共享最优信息的同时增加种群多样性,以达到增强算法鲁棒性的目的。

1 优化算法

1.1 粒子群优化算法

粒子群优化算法[1]是通过群体中的个体互相协助和群体最优信息共享来探索最优解的。粒子群中每个粒子都有一个初始速度和适应度值,个体在迭代过程中根据自身经历过的最佳位置以及种群得到的最优解的位置来生成速度,再根据目前位置和生成的速度来得到新的在空间中的位置[1]。粒子群优化算法的数学描述如(1)式所示

式中,下标i和j分别代表的是第i个粒子与维度的第j维,上标t代表迭代代数。ω为惯性权重因子;c1,c2为学习因子;r1和r2是介于[0,1]之间的随机数。

粒子群算法综合能力较强,但因为更新公式导致粒子被限制在当前全局最优解附近,在很多情况下搜索空间会被压缩,使得搜索过早停滞,最终导致算法在大部分情况下无法获得较好的结果。而布谷鸟算法凭借莱维飞行这种长短步相间的操作使其不会出现搜索空间被压缩的情况,有较强的全局搜索能力。

1.2 布谷鸟搜索算法

布谷鸟搜索算法(CSA)[2]中有2种更新位置的方式,一种是布谷鸟寻找鸟窝下蛋的寻找路径采用了莱维飞行。莱维飞行的主要特点是长短步相间,采用莱维飞行更新位置的公式如下

Xt+1=Xt+α⊗Levy(β)

(3)

式中：α是步长缩放因子；Levy(β)是莱维随机路径；⊗是点乘运算。其中生成服从莱维随机步长的公式使用Mantegna方法,公式如(4)式所示

(4)

式中，u～N(0,σ2),v～N(0,1),其中σ如下所示

(5)

式中：β∈[1,3],本文β取1.5；α是步长缩放因子,在布谷鸟算法中针对α有不同的变种,本文中α取值如(6)式所示

(6)

式中：Xi与Xj是任意不同鸟窝的位置；D为问题自变量维数。

另一种位置更新方式是宿主鸟以一定概率Pa发现外来鸟后重新建窝的位置路径,这个路径可以采用莱维飞行或偏好随机游动方式,本文使用偏好随机游动方式,更新位置公式如(7)式所示

Xt+1=Xt+r⊗Heaviside(Pa-ò)⊗(Xi-Xj)

(7)

式中：r,ò是服从均匀分布的随机数；Heaviside是跳跃函数(x>0,y=1;x<0,y=0)。

布谷鸟搜索算法通过莱维飞行可以较大程度上保证样本的多样性,但是也导致其局部搜索能力相对于PSO和DE较差,收敛速度慢,鲁棒性也较差。

1.3 基于并行交换的增强粒子群优化算法

PSO算法前期收敛速度快,后期种群多样性的消失使其容易陷入局部最优,而CSA算法的莱维飞行使其全局搜索能力较强,但收敛速度慢,缺少活力。二者鲁棒性均较差,但是各有优点,为了获得一个鲁棒性较强并兼具二者优点的算法,提出了基于并行交换的增强粒子群优化算法(EPSOBPE)。

在EPSOBPE中,由于CSA种群分走了一部分粒子数量,需要对粒子群种群进行增强,在保证全局搜索能力的前提下加速收敛,所以利用差分进化中的突变操作对粒子群全局最优解进行学习。

(8)

式中，r1,r2,r3,r4∈[1,N]的整数,N为粒子数量,且r1≠r2≠r3≠r4,g为粒子群全局最优解,F(t)为比例因子。

为了使算法在前期具有一定的全局搜索能力避免快速收敛陷入局部最优,所以前期选择开关T值应稍小,F应稍大,使得早期使用增强学习的粒子数量较少,使用增强学习的粒子得到的扰动更大,后期使用增强学习策略的粒子较多,使用增强学习的粒子得到的扰动较小,加速收敛,得到最优解。T与F更新公式如(9)～(10)式所示

(9)

(10)

式中,Gmax为最大迭代轮数。

采用增强学习策略,既可保证算法收敛,又可防止算法早熟。优化前期,粒子之间距离较远,差分变化较大,引入的扰动更强,粒子不会过早地聚集在g周边,粒子分散度大,探索能力强,对空间搜索更为充分;优化后期,粒子相互靠近,差分变化减小,粒子向g聚集,加速收敛。增强学习策略流程图如图1所示。R表示均匀分布在[0,1]区间的随机数。

由于在气动优化问题中,最优解一般位于部分维度的边界上,所以有必要对边界处理进行改进,故在边界处理中采用了减幅的处理方法,即当粒子位置超出了边界后,粒子的位置仍然重置于边界x上,但x方向的速度以一个随机大小相反方向返回到边界内,如图2所示,r表示均匀分布在[0,1]区间的随机数。

图1 增强学习策略流程图

图2 速度减幅边界处理

EPSOBPE算法策略是一种种群并行策略,2种算法分别使用改进后的PSO和CSA算法。改进后PSO算法虽然通过增强学习策略增强了局部搜索能力和收敛速度,但是在很多情况下仍不可避免陷入局部最优。CSA算法有着较强的全局搜索能力,但是收敛能力较差,所以通过分层交换策略将2种算法取长补短,共享二者优势信息使算法具有较强鲁棒性,而且既可以快速收敛,也可以保证其全局搜索能力。新算法没有通过只改进更新公式或改进种群拓扑结构来增强算法鲁棒性,而是使用多种群并行与分层交换机制。这点与其他改进算法有较大差别,有效提高了新算法的鲁棒性,使新算法在面对不同类型问题时可以获得更好的解。

EPSOBPE算法的具体步骤如下:

1) 对算法进行初始化,给定相关参数。

2) 将整个大种群等分成2个种群,分别运行改进PSO算法和CSA算法。

3) 在一轮后得到2个种群的函数最优值后进行比较,选出最优者。

4) 更新种群的最优值并进行替换,用最优者替换改进PSO种群中的最优者,或CSA种群中任意一个个体。

5) 对改进PSO种群进行一次分层,即求出改进PSO群体中所有个体的平均适应度值FPSO 1;对CSA种群进行一次分层,即求出CSA群体中所有个体的平均适应度值FCSA 1。

6) 对改进PSO种群进行二次分层,即求出改进PSO群体中适应度值优于FPSO 1的所有个体的平均适应度值FPSO 2;对CSA种群进行二次分层,即求出CSA群体中适应度值优于FCSA 1的所有个体的平均适应度值FCSA 2。

7) 判断FPSO 2是否优于FCSA 2,若满足，用PSO群体中适应度值优于FPSO 2的个体替换CSA群体中适应度值较差的那些个体;若不满足，用CSA群体中适应度值优于FCSA 2的个体替换PSO群体中适应度值较差的那些个体。

8) 判断是否满足停止准则,若满足,则结束,若不满足,则转至步骤2)。

本文在PSO算法的基础上首先进行了增强行为,增强了PSO算法的局部搜索能力;同时修改了边界处理方法。然后通过分层交换策略引入CSA算法,来增强PSO算法的全局寻优能力,最终获得了具有强鲁棒性的EPSOBPE算法。

2 函数测试

为了测试EPSOBPE算法性能,本文对PSO算法,DE算法,CSA算法,EPSOBPE算法和RADPSO算法[6]进行比较,测试函数[8]如表1所示,函数具体设置可参考文献[8]。

表1 测试函数

由于测试函数维度不同,故不同维度算法种群个数也不同,具体设置如表2所示。

表2 算法设置

在函数测试中,每个算法计算30次,结果取30次平均值。PSO算法中C1=C2=2,惯性权重因子ω随进化代数由0.9线性递减到0.4[9]。DE算法交叉因子CR=0.7,缩放因子F=0.5。RADPSO算法参数取上述PSO算法和DE算法的参数。CSA算法中Pa=0.25。EPSOBPE算法的设置同上述PSO算法和CSA算法具有相同设置,缩放因子F从0.7线性递减至0.3,T从0.3线性递增至0.7。每次寻优中算法收敛条件为调用函数次数达到100×D×N。表3展示了不同算法30次运行结果的均值和标准差。

表3 函数测试的算法性能

由表3可见,EPSOBPE算法在3个单峰函数中的寻优效果并没有因为执行粒子群优化算法的个体数减少而变差,在单峰函数中得到的结果都较好,体现了对PSO进行增强是十分有效的。但是对于测试函数f3(x)，EPSOBPE算法并没有获得最好的寻优结果,这是主要是因为测试函数f3(x)具有一定的峡谷特征,如图3所示优化算法很容易在峡谷中陷入停滞。

图3 2-Df3(x)函数示意图

同时该函数的峡谷特征正好与坐标轴相匹配,此时DE算法的变异操作可以更好地进行搜索,EPSOBPE算法由于加入了增强学习策略,在一定程度上增强了此类函数的寻优能力, 但是由于并不是所有粒子都会使用增强学习策略,所以EPSOBPE算法在此测试函数的表现上仅次于DE算法。

在之后5个多峰函数中EPSOBPE算法由于使用了增强PSO与CSA并行并进行分层交换,体现出了极强的全局寻优能力,同时,得到的标准差也较低,具有很强的稳定性。由于EPSOBPE算法在更新公式上与PSO和CSA算法较为相似,分层交换操作并不复杂,所以在计算时间方面与其余算法相似。EPSOBPE算法在不同维度不同种类的函数测试均表现良好,这体现出了EPSOBPE算法相较于其他算法的鲁棒性更强,能够适应不同种类的函数,同时具有更好的寻优能力。

3 翼型气动优化设计

3.1 翼型CFD算例验证

本文CFD二维算例验证采用经典算例RAE 2822跨音速流场数值模拟,该算例被广泛应用于二维CFD计算校核。该算例选取设计状态为自由来流马赫数Ma=0.73,雷诺数Re=6.5×106,迎角α=2.85°,设计升力系数CL=0.8。在CFD计算后与试验数据[10]结果进行对比。二维算例的网格量为3万,远场距离翼型表面40倍弦长。本文二维翼型计算时湍流模型选用k-ωSST湍流模型。

图4展示了数值模拟结果与试验值的压力分布对比。从图4可以看出,数值模拟所得压力分布与试验结果基本吻合,能够较为准确地捕捉激波的位置和强度。这说明本文采用的CFD数值计算方法在求解二维翼型的气动特性中具有较高的数值模拟精度。

图4 翼型压力分布对比

3.2 翼型优化模型

考虑到在工程中,翼型均有后缘厚度,故翼型气动优化设计的初始翼型为修改成钝后缘的RAE2822翼型,设计状态为M∞=0.73,Re=6.5×106。 巡航状态下翼型应有优秀的低阻力系数CD、小低头力矩系数CM特性,且升力系数CL应基本保持不变,同时还应考虑结构布置和机翼油箱容积等问题。故而优化目标及约束如下

(11)

式中：A为翼型围成的面积；tmax为翼型最大厚度；角标0为初始翼型相关参数。

翼型参数化使用扰动CST方法[11-12],设计变量为12个,翼型上下表面各6个。分别采用CSA算法、DE算法、EPSOBPE算法、PSO算法和RADPSO算法进行翼型优化设计,各算法种群个体数均为24,其余参数与函数测试中参数设置相同。

CFD设置与算例验证相同,采用k-ωSST湍流模型,网格数量为3万。二维翼型优化用计算机CPU为Intel(R) Core(TM) i7-10700F CPU,内存为32 GB。由于算法调用CFD求解器次数过多所以使用基于Kriging代理模型[13]的优化方法。

本文使用的基于Kriging代理模型的优化方法是首先使用拉丁超立方抽样方法(Latin hypercube sampling,LHS)生成少量的样本点构建初始代理模型,然后通过使用优化算法来求解加点带来的优化子问题得到新的样本坐标,将新样本加入样本库后再次构建代理模型,进行代理模型的更新,直到满足收敛条件。二维翼型优化问题中初始样本点数为240,使用混合加点方式[14]加点200个,共440次CFD计算。

混合加点方式是先使用最大改善期望(expected improvement,EI)加点准则,当E[I(x)]<0.01min(F(x))时使用最小预测值(minimize prediction,MP)加点准则。

EI加点准则选取改善期望最大的设计点作为加点位置。设计点x处的改善值表示为I(x)=ymin-y(x),改善期望为

(12)

式中，φ为标准正态分布的概率密度函数。对(12)式使用分部积分可得

E[I(x)]=σ(x)[uΦ(u)+φ(u)]

(13)

MP方法直接使用代理模型预测值作为目标值进行寻优,把寻优得到的具有最小预测值的设计点作为加点位置。

3.3 优化结果及其分析

图5和图6给出了优化迭代记录及优化前后几何形状与压力系数分布。

表4给出了优化后翼型的性能参数比较。

图5 翼型优化迭代记录图6 翼型几何形状与压力系数比较

表4 翼型优化结果比较

由图5可知,几种优化算法的收敛速度类似,但是接近收敛时,EPSOBPE算法的结果更好,体现了EPSOBPE算法有着更好的鲁棒性,所以可以更好地找到每次代理模型中获得最大改善的点或阻力系数最小的点。由图6可见,优化后的翼型上表面顶点后移,最高点有所下移,下表面顶点后移,最低点下移量与上表面最高点下移量相似。压力系数方面,优化后的前缘吸力峰峰值有所增大,为了保持升力,翼型上表面前部和后部的压力系数有所减少,中部消除了激波,压力系数过渡更为平滑,但也导致压力系数绝对值减少。翼型下表面压力系数与初始值相似。在优化后的翼型中,经过PSO优化的翼型与其他算法优化的翼型相差较大,其阻力系数也在优化后的翼型中最大,说明其陷入了局部最优,体现了PSO算法较差的鲁棒性。在其余算法中EPSOBPE算法得到的结果在厚度和面积上更接近初始翼型的值,且结果的阻力系数是最小的,说明EPSOBPE优化后得到的结果更接近最优解。EPSOBPE算法正如上文分析的一样,具有极强的鲁棒性,面对不断变换的代理模型所带来的子优化问题,可以获得更好的解。

4 机翼气动优化

4.1 机翼CFD算例验证

本文三维算例采用的是经典三维验证算例,ONERA M6机翼[15]流场的数值模拟算例。该算例被广泛应用于三维CFD计算校核。该算例的气动设计状态为Ma=0.839 5,迎角α=3.06°,雷诺数Re=1.172×107。该算例的计算网格量为108万。本文三维机翼计算时湍流模型选用k-ωSST湍流模型。

图7展示了0.44倍和0.9倍展向位置处截面下的数值模拟结果与试验值的压力分布对比,其中,y表示截面处位置距离机翼根部的长度,b表示机翼展长。结果表明数值模拟所得压力分布与试验结果基本吻合,并能够较为准确地捕捉激波的位置和强度。这说明本文采用的 CFD 数值计算方法在求解三维机翼的气动特性中具有较高的数值模拟精度。

图7 y/b=0.44与y/b=0.9截面处的压力分布

4.2 机翼优化模型

以M6机翼为初始机翼,分别采用CSA算法、DE算法、EPSOBPE算法、PSO算法和RADPSO算法进行机翼优化设计,各算法种群个体数均为58,其余参数与函数测试中参数设置相同。CFD计算参数设置与算例验证相同,网格与算例验证时相同均为108万。三维机翼优化用计算机CPU为2个Intel(R) Xeon(R) Gold 6248R,内存为256 GB。由于三维流场计算时间较长所以使用基于Kriging代理模型的优化方法。优化方法流程与二维翼型优化流程相同。三维机翼优化初始样本点数量为300,采用混合加点方式,加点数量为300。参数化方法使用自由变形方法(free-form deformation,FFD)[16],本文通过改变翼根与翼梢2个控制截面处的控制框点在Z方向上的位置实现对外形的改变,每个控制截面翼型上下表面各7个控制点,再对翼尖在XOY平面方向内设置1个扭转角的控制变量,综上,该算例模型共使用29个设计变量。

机翼优化设计状态为M∞=0.839 5,Re=11.72×106。优化目标为阻力系数CD最小化,约束为保持升力系数CL不变以及翼根厚度tr与翼尖厚度tt不减小。可得公式如下

(14)

4.3 优化结果及其分析

图8展示了优化过程中样本中最小阻力系数的收敛过程,表5给出了M6机翼最终的优化结果。

图8 机翼优化迭代记录

表5 机翼优化结果对比

由图8可知, EPSOBPE算法在较高维度的优化中可以适应不断变化的子优化问题,提高了优化收敛速度,体现了EPSOBPE算法的鲁棒性较强,最终可以得到与其他算法相比更好的优化结果。由表5可知,相同计算量下通过使用EPSOBPE算法,阻力系数从0.017 654降低到0.013 979,结果比其他算法中最小值仍少1.35个框。EPSOBPE优化后翼根厚度为0.079 0,与初始机翼翼根厚度较为相近,说明EPSOBPE算法得到的结果更接近全局最优解。

图9 y/b=0.99处翼型与压力系数对比图

通过图9可以看出EPSOBPE得到的翼根处翼型的头部前缘半径更小,翼型靠近后缘部分变薄并有向下弯曲的趋势。上表面相较于初始机翼的翼型更加平坦,为了补偿上表面厚度的减少,下表面厚度有所增加。通过图9展示的压力系数对比图可以看出各算法的优化结果均减弱了前缘的吸力峰,同时减弱了上表面的激波,其中EPSOBPE优化后的上表面压力系数过渡更为平缓,阻力更小。下表面靠近前缘处PSO算法和DE算法得到的压力系数与初始翼型压力系数相近,其余3种算法结果不同程度提高了压力系数。在下表面中部所有算法都有较初始翼型更低的压力系数,其中EPSOBPE算法压力系数更小,十分接近上表面同位置的压力系数。

从图10中可以看出不同优化算法得到的翼尖翼型都有不同程度的低头扭转,新的截面形状与原截面形状区别较大。在所有算法中DE算法扭转最大,RADPSO算法扭转最小,优化结果最好的EPSOBPE算法扭转程度位于中间位置。

图10 y/b=0.05处翼型与压力系数对比图

图10展示了翼尖截面的压力系数情况,可以发现各算法结果都不同程度地减小了翼尖的载荷,使得原机翼翼尖的激波得到了有效减弱。翼尖截面优化后前缘部分上下截面的压力差变小,截面中后部的压力系数由初始的几乎无压力差变为有压力差,提高了气动效率。

图11～16展示了优化前后机翼上表面压力系数云图,可以发现经过优化后各算法都减弱了机翼中部激波,DE和EPSOBPE算法结果前缘吸力的激波较其他算法结果更弱,上表面压力分布更为均匀,其他算法机翼上表面压力分布均形成了2个吸力团。DE算法上表面压力分布虽均匀,但是其机翼中部吸力峰后的压力过渡较快,故其阻力较大。EPSOBPE的上表面压力分布均匀,机翼中部吸力峰后的压力过渡较DE更为平缓,故有更低的阻力。

图11 初始机翼上表面云图 图12 PSO优化后机翼上表面 图13 DE优化后机翼上表面压力系数云图 压力系数云图

图14 CSA优化后机翼上表面 图15 RADPSO优化后机翼上表面 图16 EPSOBPE优化后机翼上表面压力系数云图 压力系数云图 压力系数云图

5 结 论

本文针对粒子群优化算法面对多峰问题容易陷入局部最优,存在鲁棒性较差的问题提出了一种基于并行交换的增强粒子群优化算法。新算法首先使用增强学习策略对粒子群算法进行了改进来增强其局部寻优能力,同时考虑到气动优化的特点,对边界处理进行了改进,得到增强后的粒子群优化算法;然后使增强后的粒子群优化算法和布谷鸟搜索算法并行,每次迭代后通过分层交换操作来交换各种群信息,有机结合了增强后粒子群优化算法的局部寻优能力和布谷鸟搜索算法的全局寻优能力,得到了EPSOBPE算法。新算法同时兼具了强局部寻优能力和强全局寻优能力,具有了极强的鲁棒性。在函数测试中,EPSOBPE算法与PSO算法等共5种优化算法在8种测试函数上进行试验对比。结果表明EPSOBPE算法在处理低维和高维、单峰和多峰优化测试函数都有突出表现,验证了新算法的优越性和鲁棒性。在二维气动优化算例中各算法在优化前期表现较为相近,但是在优化后期EPSOBPE算法逐渐与其他算法拉开差距,最终优化效果更好,展现了对气动优化问题有更强适应能力。在三维气动优化算例中,EPSOBPE算法展现出了更强的寻优能力,能够在不断变化的子问题中不断找到全局最优点,最终使机翼阻力系数优化效果明显好于其他算法,进一步验证了EPSOBPE算法的强鲁棒性,可在不同维度和不同种类问题中得到更优解。

目前提升算法鲁棒性、彻底解决算法搜索模式单一的方法除种群并行外仍缺乏其他方法。继续使用种群并行机制,然后发展其他信息交换机制或提出其他解决方法将是今后重要的研究方向。