徐建中，晏 福

(哈尔滨工程大学 经济管理学院，黑龙江 哈尔滨 150001)

0 引言

许多现实复杂问题都可以通过转化为优化问题而进行求解，优化已经成为工程和管理领域不可缺少的有机整体。优化的目的就是要在时间、金钱和能源等有限的前提下使得系统的效率、性能、产量或社会福利达到最大化[1]。因此，在现有资源受限的前提下，找到一种高效可靠的优化方法十分必要。近年来，随机搜索优化算法因为其对现实复杂优化问题求解时，表现出的求解准确度高、速度快，低成本以及强鲁棒性而受到全世界的广泛关注。与传统优化算法相比，随机搜索算法不需要任何目标问题的梯度信息，并且十分简单和易于计算机编程实现[2]。而在众多随机搜索算法中，群智能优化算法(SI)得到了广大科研工作者的研究和改进。在群智能算法中，每个个体(诸如鱼、鸟、萤火虫、蚂蚁等)为了求得生存和发展从而表现出独特的自适应行为，这些行为包括捕食、航行和搜索等，SI正是基于这些个体之间的协同行为来对复杂问题进行求解的，而大量研究成果也表明个体间的相互作用所构成的群体行为对于解决不同领域的复杂优化问题十分有效。例如，受蚂蚁搜索行为启发的蚁群优化算法(ACO)对结构优化问题[4]、交通区域控制问题[5]以及基因组优化问题[6]等都具有良好的应用效果。粒子群优化算法(PSO)对双边规划问题[7]、电力系统[8]、海洋石油存储[9]和图像处理等[10]问题的解决带来了便利。

鲸鱼优化算法(WOA)是一种受自然界座头鲸捕食行为启发的群智能优化算法[3]，它模拟了鲸鱼群搜索，包围，追捕和攻击猎物等自适应行为，从而达到优化复杂问题的目的[11]。它的优化性能要比PSO、灰狼算法(GWO)和蜻蜓算法(DA)更优，且更容易编程实现[11]。它在工程应用领域和图像分割处理中也得到了很好的应用。然而大量的群智能优化算法(包括鲸鱼优化算法)的提出，为一些特定的优化问题选择合适的优化方法带来了极大的挑战，这是因为大多数算法都只适用于广义的优化概念，并不具备对每一个优化问题都拥有先验认识。面对这种情形，组合优化算法则利用其他优化算法的优势来实现不同算法的优势互补，取长补短，从而使得组合算法的优化性能要大大优于任意单一的算法[12]。尽管WOA的全局搜索能力具有一定的随机性和静态的群集行为，但它跳出局部最优的能力有限，导致在求解目标问题的过程中容易陷入局部最优。为了克服这些缺点，本文提出了一种新的基于黄金分割搜索的组合鲸鱼算法(GWOA)。数值仿真结果和工程实例分析结果表明GWOA算法是一种高效、简单、求解性能优异和鲁棒性强的群智能优化算法。

1 鲸鱼优化算法

根据Hof和Van Der Gucht[13]的观点，鲸鱼在其大脑的某些区域有与人类的梭形细胞相似的普通细胞。这些细胞主导着人的判断、情绪和社会行为。换句话说，梭形细胞使我们有别于其他生物。鲸鱼的数量是成年人类的两倍，这是它们聪明的主要原因。事实证明，鲸鱼不但可以进行思考、学习、判断、交流，甚至还能像人类一样富有情绪化，但它的聪明程度要比人类低得多。另一个有趣的现象是鲸鱼的社会行为，即它们大多是成群出现的。他们的一些物种(例如虎鲸)可以在整个生命周期内生活在一个家庭中。座头鲸是最大的须鲸之一。一头成年座头鲸形体能比肩一辆校车。它们最喜欢的猎物是磷虾和小鱼群。图1展示了座头鲸捕猎食物的原理，它是通过沿着圆形或“9”形状的路径形成独特的气泡进行诱捕猎物，这种觅食行为被称为气泡-净喂养方法，它包括两个主要阶段：向上螺旋和双循环[11]。Mirjalili和Lewisa基于以上原理于2016年提出了鲸鱼优化算法(WOA)。

图1 座头鲸的气泡净摄食行为

X(t+1)=Xbest(t)-A·|C·Xbest(t)-X(t)|

(1)

其中，t为当前迭代次数，A·|C·Xbest(t)-X(t)|表示鲸鱼捕食的包围步长，Xbest为到目前为止迭代搜索得到的最佳位置向量，A和C的表达式为：

A=2a·rand-a

(2)

C=2·rand

(3)

a=2-2t/tmax

(4)

其中，rand为[0,1]之间的随机数，a表示收敛因子，在迭代过程中，从2线性递减到0。

鲸鱼捕食的双循环模式的数学建模如下：

X(t+1)=|Xbest(t)-X(t)|·

ebl·cos(2πl)+Xbest(t)

(5)

其中，|Xbest(t)-X(t)|表示鲸鱼i与猎物之间的距离，是一个用来定义螺旋线形状的常数，l[-1,1]是之间的随机数。

需要说明的是，鲸鱼在一个收缩的圆圈内同时沿着螺旋形的路径环绕游动。为了对这种行为进行建模，假设鲸鱼有50%的概率在缩小的包围圈内以环绕或螺旋的方式游动捕猎，并以此方式来更新其位置，则该过程的数学模型如下：

(6)

鲸鱼的主要捕食行为除了向上螺旋和双循环游动以外，还有随机捕食行为，其数学描述如下:

X(t+1)=Xrand(t)-A·|C·Xrand(t)-X(t)|

(7)

式中，Xrand(t)为从当前鲸鱼群中随机选择的鲸鱼位置向量。

根据文献[3]，WOA算法的伪代码如算法1所示。

算法1 WAO算法伪代码Initialize the whale population Xi(i=1,2,…,n)Calculate the fitness of each search agentXbest=the best search agentwhile(t

2 黄金分割搜索改进鲸鱼算法

2.1 基于黄金分割的种群初始化

对群智能优化算法而言，初始种群的好坏在一定程度上决定了算法的优化性能，WOA算法初始种群是以随机方式产生的，而这种随机方法所产生的初始种群并不能保证其多样性[11]。为了选择较好的初始种群，本文利用黄金分割的最优性来提取和过滤搜索空间中的有用信息以保证初始种群的多样性。

黄金分割法是裴波那契法(Fibonacci)的比例近似，并由我国华罗庚教授进行了最优性证明[15]。对于最优化问题:

minf(t),t⊆[LB,UB]

(8)

其中LB表示搜索空间的下限，UB表示搜索空间的上限。黄金分割搜索法对种群进行初始化就是通过不断地缩短[LB,UB]，把初始种群搜索到(8)式的近似最优解附近。

为了缩短区间[LB,UB]，逐步迭代搜索得(8)式最优解t*的近似值，所采取的迭代搜索策略为：在[LB,UB]中任取两个对称点t1和t2(不妨设t1

t1=LB+0.618·(UB-LB)

(9)

t2=UB-0.618·(UB-LB)

(10)

计算f(t1)和f(t2)并比较它们的大小。若f(t1)>f(t2)，则有t*∈[t2,UB]，从而把潜在最优解缩短在区间[t2,UB]内；相反，潜在最优解则缩短在区间[LB,t1]内，依次迭代，直到满足迭代停止条件，其迭代示意图如图2所示。

图2 黄金分割搜索示意图

利用上述不断缩短搜索区间的方法，把初始种群聚集到问题(8)的最优解附近，同时也把搜索空间[LB,UB]分割成n-1个子区间[w1,w2]，[w2,w3]，…，[wn-1,wn]，其具体步骤见算法2。

算法2 基于黄金分割的种群初始化伪代码确定搜索区间[LB,UB],并设定迭代次数n在区间[LB,UB]上选取两个对称点t1和t2,k=0w1=zeros(n-1,dim);w2=zeros(n-1,dim);while(kf(t2) LB=t2;t2=t1;t1=LB+0.68(UB-LB]); else UB=t1;t1=t2;t2=UB-0.68(UB-LB]); end if k= k+1; w1(k,:)=a*ones(1,dim);w2(k,:)=b*ones(1,dim); X=(LB+UB)/2;endwhile

2.2 变区间黄金分割非均匀变异操作

在传统的进化算法中，算子的作用与进化代数是没有直接联系的，因此，当算法演化到一定代数后，由于缺乏多样性操作和局部搜索，传统的进化算子将很难获得收益[16]，本文在文献[16]的基础上给出了一种黄金分割的非均匀变异算子。

设问题的解空间为

[LB,UB]={Xi=(X1,X2,…,Xd)T|lbi≤Xi≤ubi,i=1,2,…,d}

(11)

其中，LB=(lb1,lb2,…,lbd)T，UB=(ub1,ub2,…,ubd)。若个体Xi=(Xi1,Xi2,…,Xid)T被选择参加变异操作，则变异操作方法如下：

步骤1由2.1基于黄金分割种群初始化的迭代搜索得到n-1个子区间[w1,w2]，[w2,w3]，…，[wn-1,wn]。

(12)

其中，ω=1-t/Maxiter，t和Maxiter分别表示当前迭代次数和最大迭代次数，γ∈rand(0,1)，且

Δ(ω,y)=y·(1-rωλ)

(13)

其中，r∈rand(0,1)，λ是决定非均匀变异程度的一个参数，本文取λ=0.618。

综上，本文的变区间黄金分割非均匀变异操作的算法实现如算法3所示。

2.3 GWOA算法的时间复杂度分析和实现过程

综上，通过黄金分割搜索对WOA进行种群初始化和进行变区间非均匀变异操作，进而改进标准WOA算法，依此提出黄金分割搜索改进鲸鱼优化算法(GWOA)。为了直观的评价GWOA算法的性能，必须对该算法的时间复杂度进行分析。

GWOA是将传统的WOA和黄金分割搜索相结合来实现的。其计算时间复杂度取决于算法的实现结构，具体分为4个主要部分：(i)黄金分割搜索粒子适应度计算；(ii)黄金分割搜索粒子位置和边界更新；(iii)WOA算法粒子适应度计算；(iv)WOA算法粒子位置更新。假设黄金分割搜索的迭代次数为t1，WOA算法的迭代次数为t2，搜索代理(种群规模)为N，则(i)和(ii)阶段的计算时间复杂度为O(N2t1)，(iii)和(iv)阶段的计算时间复杂度为o(N2t2)。因此，GWOA算法总的计算时间复杂度为O(N2t1+N2t2)。综上，GWOA的计算时间复杂度要比WOA大O(N2t1)。一般情况t1是要远小于t2的，本文取t1=60和t2=500。即通过增加一定的计算时间复杂度来使得WOA算法的优化性能得到较大提高是科学的。

算法3 变区间黄金分割非均匀变异伪代码调用算法2所产生的黄金分割区间w1和w2和子区间个数n-1产生[0,1]之间的随机数p和λ,计算ω,产生随机数ηfor i=1 to N for j=1 to d if p<0.5 Xji=w1(η,j)-(w2(η,j)-w1(η,j))*(1-rωλ); else Xji=w1(η,j)+(w2(η,j)-w1(η,j))*(1-rωλ); end if end forend for

算法4 GWOA算法伪代码设置种群规模N,维数d,搜索空间的上下限[LB,UB];利用算法2产生初始种群{Xi,i=1,2,…,N}和变区间w1和w2以及区间个数n-1;令t=0;while(tf(Xi) Xi=Xi; end if else if(p≥0.5) 根据式(5)更新当前个体的位置; 根据式(12)进行变区间黄金分割非均匀变异操作; 计算适应度f(X)和f(Xi) if f(Xi)>f(Xi) Xi=Xi; end if end if end for end fort=t+1;end while

3 数值仿真及工程应用

3.1 测试函数及性能指标的选取

表1 测试函数

在目前的相关研究中，使用了许多不同的指标来衡量组合优化算法的性能。而算法性能的优劣，主要是看所提出的算法的有效性，即评估该算法在对所选测试函数进行求解最小值(最大值)时的求解精度和寻优成功率[17]。为了对本文所提出的算法寻优性能进行客观合理的评价，首先选取与文献[3]和文献[11]相同的13个标准测试函数进行仿真试验，对求解结果进行对比分析；其次，为了验证本文算法处理大规模优化问题的能力，选取了与文献[11]相同的15个测试函数，分别就不同维数进行优化运算并与IWOA的优化结果进行对比分析(D=200，500，和1000)，同时，实验参数也与文献[11]相同，即种群规模N=30，最大迭代次数tmax=500，黄金分割搜索的最大迭代次数n=60，通过30次独立运算实验，分别记录他们的平均值，标准差和收敛成功率，所有仿真实验都在Intel(R)Core(TM)i5-7200U CPU@2.50GHz 2.70GHz的计算机上依托MATLAB R2015a实现；最后，将本文的GSWOA与其他最新提出的组合优化算法进行优化性能对比分析。

依据文献[11]和文献[17]，衡量算法性能的两个指标分别是结果精确度(Accurracy，AC)和寻优成功率(Successful ratio, SR)，其定义如下：

AC=f(Xbest)-f(Xept)

(14)

(15)

其中，Xbest为函数经过tmax迭代所求得的最优解，Xopt为函数的全局最优解，z为总的实验次数，z′为算法收敛到问题全局最优解的次数。如算法迭代所得到的结果对应的AC值小于设定的收敛精度(见表1)，则认为该次运算收敛到了全局最优解。

3.2 GWOA与WOA、IWOA和GACO的性能对比分析

为了使本文GWOA算法的优化性能具有直观可比性，首先对文献[3]中的13个标准测试函数(D=30)进行数值仿真计算，其中WOA的仿真结果直接来源于文献[3]，IWO的数值仿真结果直接来源于文献[11]，其对比结果见表2。同时，为了验证本文GWOA算法的高维优化性能，对表1中15个测试函数(D=200，500和1000)进行对应的30次独立运算，IWO的数值仿真结果来源于文献[11]，其仿真结果见表4。

由表2中的对比结果可知，对于函数F1、F2、F3、F4和F10，GWOA的平均值和标准差都要优于IWOA和WOA，而对于F5，GWOA的寻优性能比IWOA和WOA要高30个数量级的精确度，GWOA在函数F7和F8的寻优性能上要略逊IWOA但要优于WOA。GWOA对函数F6、F9、F11、F12和F13的求解获得了全局最优值0，IWOA在函数F6、F9和F11处求得全局最优值0，而WOA仅在函数F9处获得了全局最优值0。

数值仿真结果的均值和标准差是对算法性能评估的有效方法，为了证明结果不是偶然产生的，还必须进行统计检验。本文利用Friedman’s检验[18]和Wilcoxon ranksum[19]检验来对算法的实验结果进行统计分析，从而对算法的综合性能进行评价。Friedman’s检验是一种常用的非参数统计方法，用于发现不同算法的结果之间是否存在显著差异，本文用它来对算法的性能进行排序，优化性能越好的算法得到的值越低，而性能较差的值反而较大。进行Friedman’s检验后，则对各算法进行Wilcoxon ranksum检验，其中，“+”表示较好，“-”表示更差，“0”表示无显著性差异，其统计结果见表3。

表2 GWOA、IWO和WOA对13个低维标准测试函数的数值仿真结果

表3 GWOA、IWOA和WOA对13个测试函数寻优的统计结果

从表3可知，GWOA的Friedman’s Test排序值最低为1.3462，IWOA的排序值为1.7692，WOA的排序值最高为2.8846，由此可知GWOA的寻优性能要优于IWOA和WOA。而Wilcoxon ranksum test结果则表明GWOA对13个标准测试函数相对于IWOA有8个测试函数性能更佳，有两个函数的优化性能稍劣，有1个函数的优化性能相当；GWOA相对于WOA，只有1个函数的优化性能稍劣以及1个函数的优化性能相当，在其余的11个函数上的优化性能要更好。

从表4的实验结果可以看出，当维数D=200时，GWOA除了对函数f5(x)和f6(x)没有获得全局最优解外，其余的函数的收敛精度都达到了表1中关于各函数的收敛精度要求。其中，GWOA在函数f8(x)、f10(x)、f11(x)和f13(x)都达到了理论最优值0。GWOA除了在对函数f4(x)和f6(x)的函数优化性能劣于IWOA外，对其余13个函数的平均值和标准差都要优于IWOA且收敛精度都达到了100%。对于D=500和D=1000，GWOA对各函数的优化性能在D=200的基础上变化不大，与文献[11]关于IWOA对高维优化函数的强鲁棒性具有相似的结论。而对于函数f4(x)，GWOA和IWOA在D=200的优化值都达到了100%的收敛精度且非常接近理论值0，分别为4.75E-172和1.37E-206；对于函数f6(x)，GWOA的寻优成功率在D=200和D=500时为40%低于IWOA的60%，而当D=1000时，GWOA的寻优成功率为30%要高于IWOA的20%。

表4 GWOA和IWOA对15个大规模测试函数的寻优结果比较

为了从统计角度说明GWOA与IWOA的优化性能，GWOA与IWOA的Friedman’s Test和Wilcoxon ranksum test的结果见表5所示，以D=200进行分析。

表5 GWOA、IWOA和WOA对13个测试函数寻优的统计结果(D=200)

由表5可知，GWOA的Friedman’s Test排序值为1.2833低于IWOA的1.7167。Wilcoxon ranksum test统计结果表明，GWOA有10个函数的优化性能要优于IWOA，有两个函数的优化性能要劣于IWOA，有3个函数的优化性能与IWOA的优化性能一致。

此外，为了进一步验证本文提出的黄金分割搜索的有效性，选取了蚁群优化算法(ACO)[20]作为验证算法，即在ACO中嵌入本文提出的黄金分割搜索算法形成GACO算法，用于对文献[3]中的13个标准测试函数进行求解。ACO和GACO的参数设置参考文献[20]，其余参数为：维数D=30，最大迭代次数tmax=500，黄金分割搜索的最大迭代次数t=60，通过10次独立运算实验，分别记录他们的最优值，平均值以及标准差，其优化结果如表6所示。

表6 ACO和GACO对13个低维标准测试函数的数值仿真结果

由表6可知，GACO除了在函数F8的优化性能要劣于ACO外，在其余的12个函数上的优化性能都要优于ACO。这表明本文提出的黄金分割搜索算法能有效改善ACO算法的优化性能。为了进一步的验证本文提出的GWOA的优化性能，由表2中GWOA与表6中GACO对13个测试函数的优化结果可知，GWOA算法在9个测试函数(F1,F2,F3,F4,F7,F9,F10,F12和F13)上的结果要优于GACO，在2个测试函数(F6和F11)上结果相等，在2个函数(F5和F8)上的结果要劣于GACO。以上分析结果表明，本文提出的GWOA算法是有效的，相较于其他对比算法的优化性能更好。

3.3 GWOA在电力负荷调度问题中的应用

发电厂的运行成本控制主要取决于发电机组的燃料成本，并通过最优负荷调度达到最小化。最优负荷调度(Optimal loaddispatch，OLD)问题，它被定义为通过对一组线上发电机组单元的发电成本控制以实现发电成本最小的目的，以满足特定时间段的总电力需求[21]。OLD问题的主要目标就是在满足一定的等式和不等式约束条件下降低发电机的燃料成本。在该问题中，发电机的燃料成本被表示为成本曲线，并通过求解总发电总量等于总功率与发电损失的总和从而计算最小的发电运营成本。

在传统的求解最优负荷调度问题的方法中，因为每个发电机的成本函数是用一个二次函数近似表示，所以必须采用迭代法、梯度法、动态规划等[21]方法来求解。一般来说，这些方法在寻找全局最优解时，通常只能提供局部最优点。此外，传统的方法还需要计算目标问题的导数，并对优化模型的函数的可导性和连续性进行一定的检验。为了克服这些缺点，许多基于自然启发的群智能算法被用来对最优负荷调度问题进行求解。粒子群优化算法是应用于最优负荷调度问题的著名元启发式算法之一[22]。其他方法用于求解最优负荷调度问题的还有诸如：差分进化算法(DE)[23]、人工蜂群算法(ABC)[24]、模拟退火算法(SA)[25]等。

最优负荷调度问题的目标是在满足不同约束条件下，并在满足电力系统的必要负荷需求的同时，将总发电成本最小化。最小化目标函数公式如下：

(16)

其中，F(Pg)表示总的燃料成本(Rs/h)；ai、bi和ci表示第i个发电机的燃料成本系数，单位分别为Rs/MW2，Rs/MW和Rs/h；n为发电机个数；Pgi表示发电机i的发电量。

1)能量平衡约束

所有发电机的总发电量与总的电量需求和相等(不考虑电量耗损的情况)：

(17)

其中，Pd表示电能需求量，单位为MW。

2)发电量约束

每台发电机的实际发电量都要控制在其最低和最高的发电范围内：

(18)

本文选取两种不同的电力优化调度案例，并用GWOA来求解，以探索其优化潜力，并与蚁狮优化算法(ALO)[26]，粒子群优化算法(PSO)[22]和萤火虫优化算法(FFA)[27]进行对比，目标函数在发电机组的功率范围内受限，同时也考虑了传输损耗。对每个测试用例执行的迭代是500，在两个测试用例中使用的搜索代理(Population)数量为30。

(1)测试案列1：3台发电机组情况

三个发电机组的输入数据如表7所示[26]。

表7 3台发电机组的输入参数

利用本文的GWOA算法求解3台发电机组得到的数据如表8所示。

表8 GWOA和ALO对3台发电机组系统的最优负荷调度优化结果

由表8的仿真结果可以看出GWOA对三台发电机组的优化调度问题求解比较理想，对于电能总需求为600MW来说，GWOA与ALO的优化结果一致。对于电能总需求为500MW时，GWOA的优化结果的总成本比ALO要稍高，但却使得发电机P1和P2的发电负荷分别降低了0.198MW和1.72MW，P3增加1.47MW的发电负荷。对于电能总需求为400MW时，GWOA的优化结果要优于ALO。

(2)测试案列2：6台发电机组情况

六个发电机组的输入数据如表9所示[26]。

表9 6台发电机组的输入参数

利用本文的GWOA算法求解6台发电机组得到的数据如表10所示。由表10可知，GWOA对6台发电机组的优化调度结果与ALO算法具有几乎一致的仿真结果，除了总的用电需求为800MW的情形下，GWOA对发电机P2的发电负载分配为16.858与ALO的10有较大差异。对于总的用电需求为600MW时，GWOA的优化性能与ALO几乎一致，而对于总的用电需求为700MW时，GWOA要略优于ALO。

表10 GWOA和ALO对6台发电机组系统的最优负荷调度优化结果

4 结论

通过引入黄金分割搜索算法对鲸鱼优化算法(WOA)的初始种群进行优化，使得初始种群能够一开始就聚集到全局最优解附近，有效的提高了WOA算法的寻优精度和全局寻优能力，同时利用黄金分割初始化种群所得到的变区间进行变区间黄金分割非均匀变异，增加了WOA粒子的多样性和跳出局部最优的能力。通过标准的测试函数进行仿真分析，仿真结果和统计结果表明GWOA相对于IWOA和WOA，优化性能得到了有效改善，同时大规模维数的仿真结果也表明GWOA有较强的鲁棒性。此外，利用本文的黄金分割算法用于改进ACO，并取得了较好的实验结果，验证了改进黄金分割搜索算法的有效性。最后，结合两个不同电力负荷优化调度问题进行工程优化应用，结果表明GWOA在电力负荷调度问题中的优化应用是高效的。下一步的研究应用将会推广到更多的工程和管理实践中，并结合神经网络用于预测和分类。