基于α-stable分布的多目标粒子群算法研究及应用

樊华羽, 詹浩, 程诗信, 米百刚

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

随着计算流体力学技术(CFD)和计算机软硬件的快速发展，各类数值优化算法结合先进CFD求解器进行高效的气动力优化设计成为研究的热点[1-2]。经过几十年的发展，在气动优化设计领域先后开发出了基于遗传算法[3]、粒子群算法[4]和差分进化[5]等智能优化算法的设计模式，对包括二维翼型、三维机翼、全机构型以及飞发一体化等问题进行了有效的气动优化设计研究，取得了丰硕的研究成果[6-7]，其中应用最普遍的优化算法之一就是粒子群算法(PSO)。

粒子群算法通过模拟鸟群的觅食行为，将解空间中的设计变量当成一群质量和体积均为0的“鸟”，也就是所谓的粒子，将所求解问题的解当成鸟群寻觅的“食物”，以单个粒子的最优解和整个群体的最优解来引导下一代种群的搜索，以此来趋向最优解。该算法属于进化算法的一种，与遗传算法类似，也是从随机解出发，通过迭代寻找最优解，使用适应度来评价解的品质。相比遗传算法，粒子群算法规则简单，没有“交叉”和“变异”操作，仅是通过追随当前搜索到的最优值来寻找全局最优，容易实现，精度高、收敛快，但是该算法的缺点在于进化后期的早熟现象导致算法容易陷入局部最优。为此，研究人员先后建立了多个改进算法来提高其全局寻优能力。李丁等[8]将粒子群算法和差分进化算法结合，使得粒子群和差分群分层交换，共享最优信息，降低了算法陷入局部最优的风险，改进的算法在对RAE288翼型的减阻优化中取得了良好的效果；陈进等[9]提出一种反双曲余弦函数形式的非线性惯性权重表示法，改进了标准的粒子群算法，并对不同厚度的风力机翼型进行了优化，提高了相应的气动特性；李鑫等[10]采用基于二阶振荡和自然选择的改进粒子群算法对二维SC-0717超临界翼型进行了优化设计，显著提高了升阻比并降低了气动噪声。

尽管改进型算法改善了原算法容易陷入局部最优的问题，但是多数研究仍然局限在单点单目标状态，而气动优化设计是一种高维强非线性问题，考虑多个优化目标才能更好发挥设计的作用。为此，本文针对传统的粒子群算法进行改进，提出了一种基于α-stable分布的多目标粒子群算法。该算法将α-stable分布理论引入传统的粒子群算法，旨在改善粒子群智能算法的全局/局部搜索能力平衡问题。使用多个基准多目标函数对改进的多目标粒子群算法的性能进行了测试，并将改进后的多目标优化算法应用到RAE2822的气动优化设计中，对比分析了改进前后的算法在考虑减阻和力矩绝对值不减小2个目标下的优化效果。

1 基础知识

1.1 多目标优化问题

本文对于一个n维搜索空间S∈Rn,假定fi(x),i=1,…,m是定义于S的目标数量为m的多目标优化向量,那么带约束的最小化连续多目标优化问题(muliti-objective optimization, MO)的通用数学模型可以定义如下

(1)

式中，x=[x1,x2,…,xn]是决策变量向量,变量个数为n,l和p分别是不等式约束和等式约束的数量,gj(x)是第j个不等式约束,hk(x)是第k个等式约束。目标函数就是将决策空间映射到目标空间,即S→Ω∈Rm。

下面给出多目标优化问题中常用的Pareto最优相关定义。

定义1对于任意2个向量u=(u1,…,um)和v=(v1,…,vm),当且仅当:∀i=1,2,…,m,ui≤vi∧∃j=1,2,…,m,uj

定义2有且仅有在搜索空间S中不存在任何向量x∈S使得f(x)f(x*),那么就称决策向量x*∈S为Pareto最优解,或者称为非支配解。

定义3在决策空间S内的Pareto最优解的集合被称为Pareto最优解集,记为P*。Pareto最优解集对应的目标函数值称为Pareto最优前沿,为F*={f(x):x∈P*}。

1.2 基本粒子群算法

粒子群优化算法(particle swarm optimization,PSO)最初是为了模拟鸟群的飞行觅食行为而提出的。其初始化为一群随机粒子,通过迭代寻找最优解。在每一次迭代中,粒子通过跟踪2个极值来更新自己:①粒子本身所找到的最优解,称为个体最优值pbest;②整个种群目前找到的最优解,该极值称为全局最优值gbest。在找到上述2个最优值时,粒子根据(2)式和(3)式来更新自己的速度和位置:

(2)

(3)

式中,c1和c2分别为认知和社会加速系数,取值范围为[0,4],一般取为c1=c2=2;r1和r2为2个在[0,1]内服从均匀分布的随机变量。

1.3 多目标粒子群(CDMOPSO)算法

Coello在2002年首次提出了多目标粒子群(MOPSO)算法[11],其独特的搜索机理使其特别适合处理多目标优化问题。在此基础上,Raqual于2005年提出了CDMOPSO(crowding distance multi-objective PSO)算法[12],CDMOPSO将拥挤距离算子(见图1)[13]引入PSO算法中,用于选择引导粒子运动的最优解以及维护存储非支配解的外部档案,增加了整个种群的多样性,从而避免算法在早期就陷入局部最优。此算法的优点是结构简单,易于实现。其具体流程见图2。

图1 拥挤距离计算示意图

图2 CDMOPSO流程图

2 ASMOPSO算法

2.1 α-stable分布

α-stable分布是一个四参数分布函数族[14-15]是一类广泛应用的随机信号模型,记为S(α,β,γ,σ),其中,α∈(0,2]是稳定性参数,β∈[-1,1]是偏斜参数,γ∈(0,∞)是放大系数,σ∈(-∞,∞)是位移参数, 参数α和β决定了α-stable分布的形状。

α-stable分布是一个非常丰富的分布族,它有3种特殊情况:

2) 当α=1,且β=0时为柯西分布,记为S(1,0,γ,σ);

3) 当α=0.5且β=1时是Levy分布,记为S(0.5,1,γ,σ)。

它的概率密度函数可以用特征函数的连续傅里叶变换来定义

(4)

式中,sgn(t)是t的符号,Φ表示为

(5)

2.2 α-stable变异

在α-stable分布中,稳定性系数α是一个非常重要的参数,它描述了分布的尾迹大小,决定了生成的随机数分布范围。图3给出了不同α值对应的α-stable分布。

图3 α对分布密度的影响

本文的PSO算法中采用α-stable分布产生随机数,对PSO种群中的个体进行变异操作。变异过程中通过动态调整α-stable函数的稳定性系数α来实现变异范围和幅度的变化,变化范围为[1,2],其变化过程见图4。在PSO程序循环的开始阶段,使用较小的α值,此时变异能力强,有助于增强算法的全局搜索能力,避免算法过早陷入局部最优;随着α值逐渐变大,全局变异能力减弱,但是局部搜索能力增强;最终α值增大为2,此时执行的是传统的高斯变异,局部搜索能力最强,有利于提高解的精度。这种动态α-stable变异可以很好地平衡算法的开发和探索能力,使得改进后的多目标粒子群优化算法ASMOPSO算法兼顾精度和全局搜索能力,具有更好的适用性。α-stable变异操作的伪代码见图5。

图4 α值的变化过程

图5α-stable变异伪代码

2.3 本文提出的算法(ASMOPSO)

本文提出的ASMOPSO算法是在CDMOPSO算法的基础上实现的,其详细实现步骤如下:

步骤1 设置粒子群种群POP规模为N,加速系数c1和c2,惯性系数w,外部档案A容量NREP,变异概率pmut,给定最大迭代步数Imax,设置外部档案A=φ;

步骤2 对于种群中的每个粒子(i=1,…,N),随机初始化位置Xi,初始速度设置为Vi=0,计算初始目标向量值f(Xi),初始化粒子的最优位置pi=Xi;

步骤3 令t=0;对种群中所有粒子进行Pareto比较,将非支配粒子储存在外部档案A中;

步骤4t=t+1;计算外部档案A中的非支配解的拥挤距离,将非支配解按照拥挤距离值大小降序排列;

步骤5 对于POP中的每个粒子(i=1,…,N):从外部档案A中顶部一定比例(例如10%)的非支配解中随机选择一个作为全局最优粒子pg;用公式(2)和(3)更新粒子的速度和位置;如果Xi超出搜索空间边界,那么就设置粒子的位置为边界值,同时将速度乘以(-1),使粒子反向移动;使用2.2节的方法对粒子执行α-stable变异操作;计算粒子的目标向量值f(Xi);更新粒子个体最优位置;

步骤6 对种群POP中所有粒子进行Pareto比较,并将新的非支配粒子与外部档案A中非支配粒子进行Pareto比较,将外部档案A中的支配粒子删除;如果外部档案A中的非支配解的数量等于NREP,那么进行下列操作:①计算档案A中的每个非支配粒子的拥挤距离;②将外部档案中的个体按照拥挤距离值降序排列;③从A中最拥挤区域(拥挤距离值最小的10%)中随机选择一个,然后将其用新的非支配解代替;

步骤7 如果t=Imax,结束程序,否则转到步骤4。

3 基准函数测试

3.1 评价指标

要对多目标优化算法的性能进行全面的量化分析,需要从它们的收敛性能、解集分布的多样性和宽广性三方面进行比较。本文使用反转世代距离(inverted generation distance,IGD)来综合评价算法的性能,IGD既能评价多目标解集的收敛性,又能在一定程度上评价解集的分布特性。其定义如下

(6)

式中,d(v,P)表示近似解集P中的个体到真实前沿P*中个体v之间的最小欧氏距离;该算法求得的最优解集用P表示。IGD指标值越低,优化结果越理想,算法性能越优异。

3.2 参数设置

本文选择ZDT1、ZDT2、ZDT3、ZDT4、Tanaka以及Srinivas 6个测试函数来评估ASMOPSO算法的性能[16-18],并与CDMOPSO算法进行对比。

表1给出了6个基准函数的基本属性。6个测试函数中,ZDT系列函数属于无约束的函数,而Tanaka和Srinivas属于带有约束的多目标测试函数。ZDT系列测试函数的设计变量维度高,优化难度大,特别是ZDT3函数的Pareto前沿是非连续的,ZDT4函数是一个多峰值函数,拥有多个局部最优前沿。Tanaka的Pareto前沿曲线被分为三段,其中中间段的水平和垂直部分由于解之间的差异比较小,是最难搜索的部分。Srinivas是一个具有2个设计变量和二目标带约束的优化问题,但其具有的2个非线性约束项给优化问题还是带来了一定的困难。CDMOPSO和ASMOPSO算法的初始种群大小均取为100,最大循环次数为200,惯性权值w从0.9线性递减为0.4,加速系数c1=c2=2.0,变异概率为0.4。

表1 基准测试函数

续表1

函数名称维数变量范围函数目标与约束注解ZDT230[0,1]minf1(x)=x1minf2(x)=g1-f1g()2()g(x)=1+9∑ni=2xi(n-1)凹函数 ZDT330[0,1]minf1(x)=x1minf2(x)=g1-f1g-f1g()sin(10πf1)()g(x)=1+9∑ni=2xi(n-1)非连续凸函数 ZDT410x1∈[0,1]xi∈[-5,5]i=2,…,nminf1(x)=x1minf2(x)=g1-f1g()g(x)=1+10(n-1)+∑ni=2(x2i-10cos(4πxi))多模态凸函数 Srinivas2[-20,20]minf1(x)=2+(x1-2)2+(x2-1)2minf2(x)=9x1-(x2-1)2s.t. e1(x)=x21+x22≤225 e2(x)=x1-3x2+10≤0凸函数 Tanaka2(0,π]min f1(x,y)=xmin f2(x,y)=ys.t.-x2-y2+1+0.1cos16arctanxy()≤0 x-12()2+y-12()2-12≤0非凸不连续函数

3.3 函数测试结果分析

对每个基准函数进行30次独立运行之后的IGD值统计结果列于表2,表中AVE表示平均值,STD表示标准差。从表中数据可以看出,本文提出的ASMOPSO在4个ZDT函数上的性能都优于CDMOPSO,尤其在优化难度较大的ZDT3和ZDT4函数上的优势明显;对于2个带约束的函数,CDMOPSO的IGD均值都略小于ASMOPSO,但是差距都非常小。

表2 CDMOPSO算法与ASMOPSO算法的IGD指数比较

图6和图7给出了2种算法任意一次优化搜索到的近似Pareto解与真实Pareto前沿的比较。可以看出改进后的ASMOPSO算法获得非支配解与真实Pareto前沿的偏差很小,而且都均匀地分布在真实前沿上,特别是在求解ZD4测试函数时表现出了非常优异的性能。图9和图10表明CDMOPSO在处理ZDT4这种多模态问题时寻优能力欠佳,同时CDMOPSO生成的ZDT3前沿也没有均匀地分布在真实前沿上,且Pareto解的数量明显偏少。

图6 CDMOPSO测试函数结果

图7 ASMOPSO测试函数结果

4 工程应用

为了进一步验证本文提出的多目标粒子群算法的优化性能,将其应用于RAE2822翼型跨音速优化设计。

4.1 优化模型的建立

以RAE2822翼型在Ma=0.73,α=2.79°,Re=6.5×106时最小化阻力系数和俯仰力矩系数绝对值不减小为目标进行外形优化设计[19]。约束条件为升力系数大于原始翼型的升力系数值,同时保持翼型的最大厚度不减小。由此可得优化的数学模型为

(7)

式中，X是设计变量向量,Cd是阻力系数,Cm是力矩系数,Cl是升力系数,cmax是翼型的最大厚度,上标*表示初始翼型的计算值。

本文的翼型参数化使用型函数为6阶Bernstein多项式的CST方法。翼型上下表面各取6个设计参数,共12个优化设计变量。图8～11显示了RAE2822初始翼型的CST参数化结果。从图中可以看出,CST参数化后的上翼面外形与原翼型几乎完全重合,而下翼面在x/c=0.1～0.5段误差较大。从数值上看上翼面的最大拟合误差小于0.000 2,下翼面的最大拟合误差小于0.000 5,而在工程中对于翼型最大拟合误差小于0.001(弦长为1时),则认为该参数化方法能精确拟合当前翼型。从图12可知经过CST参数化后的拟合翼型表面的Cp分布基本与原翼型重合,下表面有较小的差异,主要集中在x/c=0.2～0.5这一段,上表面的Cp分布的差异较大,在前缘附近未能很好的捕捉负压鼓包区,在翼型中部,激波出现的比原始翼型稍早,但是差异不是很明显;而表3则直观说明了2种翼型的气动力特性差别非常小,升力系数和阻力系数的相对误差都小于1%,力矩系数误差略大,为1.6%。参数化结果表明本文选用的基于6阶Bernstein多项式的CST方法可以较为精确地描述RAE2822翼型的设计空间。

图8 上翼面外形拟合结果 图9 上翼面外形拟合误差 图10 下翼面外形拟合结果

图11 下翼面外形拟合误差

图12 RAE2822翼型CST参数化后压力系数比较

参数值原始翼形CST-BPn6误差/% Cl0.758 00.754 6-0.45 Cd0.017 70.017 6-0.56 Cm-0.087 3-0.085 9-1.60

4.2 优化结果及分析

采用本文改进的ASMOPSO和CDMOPSO 2种方法对RAE2822翼型进行跨音速多目标优化设计并对比两者的优化效果。2个优化算法的种群规模均为100,两者的外部档案大小为100;惯性系数为0.4～0.9。认知加速系数和社会加速系数C1=C2=2.0。ASMOPSO优化算法中的系数α变化范围为1.0～2.0。2种算法的迭代步数均为100步。

由图13可知,多次优化后,ASMOPSO优化算法得到的非支配解集分布均匀,而且比CDMOPSO的优化结果更接近真实Pareto前沿。所以后续的对比中,仅需要与原始翼型对比即可。由于多目标解集更能客观的反应优化问题的真实内涵,选取点更能符合不同设计人员的需求,那么在均衡操纵性和减阻效果下,如图13所示,取Pareto前缘的中间一点。在此点下的气动数据如表4所示。图14与图15分别为原始翼型与优化后翼型的几何外形比较与Cp分布比较。可以看出优化后翼型,最大厚度位置向后移动,前缘压力峰值和上翼面中部的局部压力变化有明显的减缓。图16与图17的压力云图可以看出优化翼型激波明显减小。因此在此取点下的优化结果具有良好的阻力特性和力矩特性。

图13 RAE2822翼型优化后的Pareto前缘

表4 优化前后的气动数据比较

图14 翼型表面压力分布对比

图15 优化翼型外形对比

图16 原始翼型压力云图

图17 ASMOPSO优化后翼型压力云图

5 结 论

本文针对CDMOPSO算法容易陷入局部最优的缺点,结合α-stable分布理论,建立了改进的ASMOPSO多目标优化算法。使用多个基准函数对比测试了改进前后的算法性能,并将改进算法应用到了RAE2822翼型跨声速多目标优化设计中,得到的结论有:

1) 引入α-stable分布理论在PSO算法中对种群进行变异操作,根据α-stable分布随稳定系数α来的特点,通过改变α值来动态调整优化过程中的变异强度和范围,使得改进后的ASMOPSO算法较原始算法具有更好的局部和全局搜索兼顾能力。

2) 无约束的ZDT系列函数和带约束的Tanaka和Srinivas的测试结果显示,改进后的算法能够有效引导非支配解逼近真实Pareto前沿,尤其是对多峰值ZDT4函数的求解,表明改进后的算法在处理多模态问题时有很强的寻优能力;

3) RAE2822翼型的跨声速减阻和俯仰力矩系数绝对值不减小的综合优化设计中,新的ASMOPSO算法得到的结果明显优于CDMOPSO。