一种多尺度协同变异的萤火虫粒子群混合算法及在车间调度中的应用

周艳平,刘永娟

(青岛科技大学 信息科学技术学院, 山东 青岛 266061)

0 引言

作为有效的资源配置与优化手段，生产调度在整个制造体系中扮演关键的角色，应用范围非常广泛，包括工业、商业等各方面。作为整个生产调度的关键部分，车间调度同时也是制造类企业的核心问题，是一个典型的NP-hard问题[1]，一直以来也是组合优化领域内的重点内容，对整个生产流程能否顺利完成有很大的影响。车间调度问题的本质是在符合约束条件的基础上，根据生产指标，给出每个工件的加工顺序、机器、操作方式等的组合问题。流水车间调度作为其中的一类，是流水线生产调度的简化模型，对流水车间调度模型进行研究，有非常重要的意义。

在研究的初期，学者们通过使用精确算法来求解流水车间调度问题，但存在难度大、耗时长等问题，求解效果并不理想。后来研究者们采用近似算法去寻找问题的最优解，近似算法又叫做启发式算法，不断被研究改进的启发式算法成为了调度问题的主要求解办法。

粒子群算法(PSO,particle swarm optimization)是美国的社会心理学家Eberhart和电气工程师Kennedy[2]受达尔文“物竞天择、适者生存”的启发，模拟了鸟群群聚、觅食的群体性行为，于1995年提出的，算法步骤简单、参数少并且易于实现；英国剑桥大学工程系的X.S.Yang[3]在对萤火虫种群的发光行为进行建模的基础上提出了萤火虫算法(FA， firefly algorithm)，算法的提出时间不长，对算法的改进和优化仍然处于初级阶段。但不断被研究改进的萤火虫、粒子群算法已被应用于车间调度、资源调度、组合优化等领域。

J.I.FISTER等人[4]提出了一种模因萤火虫算法(MFFA)，将萤火虫算法与局部搜索启发式算法进行结合，并应用于组合优化问题中，具有一定的实用性。郑捷等[5]采用机器选择策略进行种群生成，并对适应度较优的个体之间进行局部搜索，从而产生了一种改进的离散萤火虫算法，并成功地应用在柔性作业车间调度模型中。刘长平等[6]提出了一种DFA算法，用NEH构造初始解，并根据问题的特性，重新定义了距离公式，通过交叉、变异操作实现位置的更新，并对处在最佳位置的个体进行局部搜索，成功地应用在零空闲置换流水车间调度中。张其亮等[7]结合了粒子群、贪婪算法，通过变异操作来改变粒子停滞或群体陷入局部最优状态，并依据模拟退火接收变异个体，在置换流水车间调度问题的求解质量上有一定的提高。吴琼等[8]结合了粒子群、引力搜索算法，当粒子群进化停滞时，引入引力搜索算法，并采用协同进化的思想，使得两种算法并行优化，并且子种群之间交流，以获得最优粒子，并应用到对作业车间调度问题的求解中。从实验数据来看，以上算法在不同程度上都提高了算法的性能，但在求解精度上还可以继续改进，而且随着车间的日益复杂化，求解算法也应该不断的被改进完善。

单一算法在适应性上是有限的，萤火虫和粒子群算法有一定的相似性，在全局搜索能力方面也有一定的互补性，为更好地发挥两种算法的潜能，取长补短获得更优的解，本文在萤火虫粒子群混合算法的基础上，做出了进一步的改进，提出了一种多尺度协同变异的萤火虫粒子群混合算法(HFPMCV, a hybrid firefly and particle swarm optimization algorithm with multi-scale cooperative variation)，引入多尺度协同变异算子，提高求解精度，混沌初始化种群以提高初始解质量，根据高斯拟合把种群分为两组；平行进化，提高求解速度，函数优化实验验证了改进算法的有效性，流水车间调度问题实验，表明了改进算法的优越性。

1 多尺度协同变异的萤火虫粒子群混合算法

1.1 萤火虫算法

萤火虫算法的提出是建立在萤火虫种群中个体的种种社会行为的基础之上的，在实际的应用中，只考虑发光特性，通过发光特性搜索目标附近特定区域，向一定范围内定位较优的萤火虫移动[9]，萤火虫进化和目标函数寻优的最主要的两个关键因素就是荧光亮度和吸引度。亮度表明萤火虫所处方位的好坏，决定了飞行方向，吸引度影响萤火虫的移动距离，利用吸引度和亮度更新、迭代来寻求目标函数最优解。数学描述如下：

在D维搜索空间内萤火虫i的位置向量Xi=(xi1,xi2,xi3,…xiD)，萤火虫荧光亮度公式：

I=I0e-γrij

(1)

式中，I0为自身亮度，γ为光吸收系数，r为距离。

萤火虫之间的相对吸引度正比于亮度：

(2)

式中，β0为其初始吸引度。

萤火虫移动距离：

(3)

α是扰动的步长因子，rand()是一个随机扰动，取值为[-0.5,0.5]范围内的均匀分布，α取值为[0,1]范围。

最亮的萤火虫位置移动：

(4)

中心思想就是，所有个体都飞向比自己更亮的个体四周，飞行的最后结果即所有的个体都围在了亮度最强的个体周围，从而达到算法求解的目的。

萤火虫算法 FA：

1)初始化参数，萤火虫数目n，吸引度β0，步长因子α，迭代次数Generation，搜索精度ε。

2)初始化个体位置，最大荧光亮度I0用萤火虫的目标函数值表示。

3)由1)、2)计算群体中个体的I、β，飞行方向由相对亮度决定。

4)由3)更新个体位置，由4)移动最佳个体。

5)计算个体更新位置后的亮度。

6)满足终止条件则输出最优解，否则转3)步。

1.2 粒子群算法

粒子群算法的核心思想在于利用粒子与种群之间的信息改变粒子的速率和定位，在飞行范围内，粒子随机飞行，飞行过程中通过学习、模仿周围粒子的信息来获取个体、全体最佳位置，综合分析经验后，调整自身的飞行经验以实现更好的飞行，完成种群间个体协作探索，利用速率和定位的不断更新迭代来寻求目标函数的最佳解[10]。数学描述如下：

D维搜索区域内，任意一个粒子都在以一定的初速度随机飞行，粒子i的当前位置向量Xi=(xi1,xi2,xi3,…xiD)，第i个粒子当前速率：

Vi=(vi1,vi2,…viD)，i=1,2,…N

(5)

Pbest=(pi1,pi2,…piD)i=1,2,…,N为个体极值,全局极值为gbest=(pg1,pg2,…pgD)。

粒子进行速度更新：

vid=w*vid+c1r1(pid-xid)+c2r2(pgd-xid)

(6)

ω是惯性权重，平衡全局、局部搜索能力，c1、c2是学习因子，r1、r2是[0,1]里的随机数。速度公式包含三部分，w*vid表示个体之前的速度对现在速度的影响，c1r1(pid-xid)是个体历史最好位置对现在速度的影响，c2r2(pgd-xid)是个体之间的信息传递对现在速度的影响。

粒子进行位置更新：

xid=xid+vid

(7)

c1、c2是学习因子，r1、r2是[0,1]区间的随机数。

步骤如下：

粒子群算法 PSO:

1)设置参数，包括种群规模、个体的速度和位置，并设置其它参数。

2)初始化种群。

3)计算适应度值，选出个体、全局极值。

4)粒子飞行，据式(6)、(7)更新个体的速度、位置。

5)符合终止条件就输出最优解，反之返回3)。

1.3 多尺度协同变异的萤火虫粒子群混合算法

考虑到萤火虫算法和粒子群算法的优缺点，杨小东等提出了萤火虫粒子群混合算法FAPSOMA[11]，本文在此基础上引入多尺度协同变异算子，提出了HFPMCV算法，初期，通过多尺度变异算子快速找到全局最优解空间，不断迭代，靠近最优解领域时，通过小尺度变异算子快速收敛到最优解,以有效的进行变异，并且对种群进行混沌初始化，同时动态自适应策略把种群分成两组，两组种群平行进化以寻找优化问题的精确解。

1.3.1 改进策略

1)混沌初始化：随机初始化会造成个体位置分布不均，个体间差异较大，阻碍全局寻优，影响求解效果[12]，为了获得质量较好的初始解，对种群进行混沌初始化，将变量映射到混沌变量的取值区域，通过线性变换把混沌序列转化到搜索空间[13]。

由Logistic 映射：

Zn+1=μZn(1-Zn)，n=0,1,2,…

(8)

得混沌序列，μ为控制参量，μ∈[0,4]，当μ=4时为完全混沌状态。

将混沌序列映射为求解空间范围的混沌变量：

pn,i=ai+(bi-ai)Zn,i，i=1,2,…，I

(9)

Zn,i为个体在第i维迭代n次得到的混沌序列值，I为搜索空间维度，pn,i为一个混沌变量，[pn,1,pn,2,…pn,I]代表一个可行解。

因为PSO算法主要更新速度、位置、目标函数值，而FA算法主要更新位置和目标函数值，本文用FA算法中的个体位置和目标值来替代PSO算法中的位置和目标值，并在FA算法中创建速度信息，所以在本改进算法中，创建FA算法的速度序列：

zi(n)=μzi(n-1)(1-zi(n-1))

(10)

zi(n)是[0,1]中符合正态分布的随机数，n是迭代数，μ值为4。

2)动态自适应控制策略：在传统算法中，种群分组采用随机分组、平均分组的方法，没有充分考虑种群中个体的运动轨迹。

HFPMCV算法采用高斯拟合函数[14]进行种群分组，在混沌初始化种群后，计算并按照种群中各个体的适应度值进行高斯拟合，因适应度值是由目标函数得来，因此拟合曲线与目标函数曲线有较强关联度[15]，求解拟合曲线峰值的加权平均数[16]，并以此为根据把整个种群划分为两组族群。

3)多尺度协同变异算子：为了提高求解精度，避免陷入局部最优，引入多尺度协同变异算子[17]，多尺度协同变异利用方差不同的高斯变异机制探索解空间，初期，通过大尺度变异算子快速找到全局最优解空间，不断迭代，靠近最优解领域时，通过小尺度变异算子快速收敛到最优解,以有效的进行变异[18]，设尺度个数为M，首先初始化方差:

(11)

对方差、优化变量设置一样的取值范围，随着迭代次数的增加，方差会一直变化，所有个体生成M个子群，个体个数为P=N/M，迭代次数为K，子种群适应度：

(12)

f是适应度函数，第m个变异算子的标准差：

(13)

(14)

(15)

变异算子的进化过程为一个递归的过程，HFPMCV算法要求的是随着迭代次数的增加变异算子逐渐减小，因此需要规范变异算子的标准差：

(16)

W为待优化变量空间的宽度。

位置更新公式如下：

If(vid

Elsevid=rand*Vmax

(17)

其中：Td是设置好的阈值，阈值Td过大会降低算法的深度局部搜索能力，阈值过小会降低算法的搜索速度，因此，自适应设定阈值，不同维的速度设置不同的阈值，当达到该值时，阈值会自动下降，当速度小于阈值时，个体进行逃逸：

其中:

(18)

ifGd(t)>k1then

Gd(t)=0;Td=Td/k2

(19)

k1、k2是常数，Size为种群大小，Gd(t)是种群在第d维速度的逃逸次数。

4)平行进化：结合FA和PSO的优点，使FA和PSO平行的单独进化[19]，每一代进化结束后，各自保存最优个体，选出较优个体，并与父代最优个体进行比较选出更优个体作为本代最优个体，全部迭代结束后，输出整体最优值，大大提升了求解速度[20]。

1.3.2 算法流程

HFPMCV算法流程如图1所示。

图1 HFPMCV算法流程图

1.3.3 函数优化测试

为了测试和解释HFPMCV算法的性能和优越性，将HFPMCV算法与FA算法、FAPSOMA算法进行了对比，分别采用HFPMCV算法以及FA算法、FAPSOMA算法对以下函数[21]进行寻优：

xi∈[-5.12,5.12]，i=1,2,…,30

(20)

是一个多峰函数，在30维可行域中有230-1个局部求解，最优值是0，对3个算法中的各参数设定一样的值以提高计算有效性和结果精确度，迭代500次，种群规模是100，维度是30等。

3种算法的性能对比如图2所示，图中可以看出，迭代次数相同时，HFPMCV算法收敛更快，精度更高，求解结果相同时，FA、FAPSOMA需要迭代更多次才能得到最优解，通过函数测试验证了HFPMCV算法的优良性能。

图2 3种算法性能比较图

2 HFPMCV算法求解流水车间调度问题

2.1 流水车间调度问题模型

FSP问题[22]:n个工件在m台机器上加工，i工件工序数为Li，总工序数L，工序加工时间固定，并按工序加工。在约束条件下安排工件加工顺序以优化指标。

具体表述如下：

(21)

凡事都有一个标准，评价调度方案好坏的标准就是是否达到指标，如客户满意度指标、完工时间指标、性能指标等。本章选取的指标为最大完工时间，以最大完工时间中的最小值作为优化目标[23]，目标函数如下：

min{maxTi,j},i=1,…,n;j=1,…,m

(22)

maxTi,j为所有工件中最后完成加工的结束时间；min{maxTi,j}为最大完工时间中的最小值。

以流水车间调度问题为模型，公式(21)中的约束如下：

1)一个机器不可同时加工多个工件。

2)每个待加工工件的加工工序是固定的，只有完成了前面工序的加工工作，才可以进行之后工序的加工工作。

3)工件零等待，完成当前工序加工的工件必须要立刻被运输到下台机器，准备加工。

4)当前机器没有需要加工的工件时可以进行等待，直到工件运输到该机器。

5)工件、机器已知，明确机器上加工工件的时间。

2.2 HFPMCV算法求解流水车间调度问题步骤

HFPMCV算法求解流水车间调度问题的求解步骤如下：

输入：种群规模popsize、迭代次数iteration、搜索范围bound、维度vardim

输出：工件加工最优次序、最小完工时间

步骤1：根据流水车间调度问题，采用基于工件的编码方式，对于整个种群，将位置向量xi[D]，由公式(8)生成向量xi+1[D],xi+2[D],…,xN[D]，N为种群数，通过公式(9)将xi[D],xi+1[D],xi+2[D],…,xN[D]映射到位置变量的取值空间，获取最初位置X1[D],X2[D],…,XN[D]，生成n个类似[pn,1,pn,2,…pn,I]的可行解，即采用混沌映射的方式生成一个可以进化的初始种群。

步骤2：根据动态自适应控制策略，计算并依据初始适应度值进行高斯拟合，把整个种群划分为两组族群。

步骤3：根据初始适应度值，选出两族群各自的最优个体，通过比较，选出较优个体作为本代最优个体，即第一代最优加工次序。

步骤4：开始平行进化，FA和PSO分别搜索最优值：PSO加入多尺度协同变异算子，据式(17)更新速度、位置，重新计算适应度值，FA加入多尺度协同变异算子，据式(1)、(2)更新亮度、吸引度，据式(17)更新位置、速度，重新计算适应度值，FA和PSO各自保存本代最优个体。根据式(13)～(15)分别更新两族群中多尺度协同变异算子的方差，根据变异次数更新阈值。

步骤5：判断两者保存的最优值哪一个更优，比较更优个体与上一代最优个体，选出较优者作为本代最优个体，即最优方案。

步骤6：达到最大迭代次数则输出最优加工次序和相对应的最优加工时间，否则转步骤4，开始下一次并行迭代。

2.3 HFPMCV算法求解流水车间调度问题复杂度分析

G*O[O(m*D)+2*O(M*D)+O(M)+O(M)

整理后为：

O(G*m*D)+7*O(G*M*D)+2*O(G*M)可近似看作O(G*M*D)，即复杂度与基本算法仍处于一个数量级。

2.4 仿真实验

2.4.1 实验准备

在Windows10操作系统下的PyCharm软件上运行程序，采用python语言编程。种群规模为50，迭代800次，保证几个算法参数一致，基于工件编码，选择TA类车间调度数据集[24]，用20*5、50*5及100*5这3个数据集为基础，分别选择每一个数据集下面的一个典型问题，来进行流水车间调度问题的实验以验证算法。

2.4.2 实验结果与分析

为了测试和解释HFPMCV算法的性能和优越性，分别使用FA、FAPSOMA、HFPMCV对3种不同规模的调度集进行对比实验，并以最小化最大完工时间为指标，为避免实验的偶然性，3种算法在对3个不同规模的数据集进行求解的时候，分别运行10次，把10次结果的平均值和最优值作为结果参考指标。实验结果如表1～3所示。

表1 3种算法在规模为20*5下的输出结果

表2 3种算法在规模为50*5下的输出结果

表3 3种算法在规模为100*5下的输出结果

仿真结果如表3所示，对于同一个规模的数据集，单从最优值来看，例如50*5规模，FA的最优值为2 886，FAPSOMA的最优值为2 837，HFPMCV算法最优值为2 739，即HFPMCV算法可以得到比FA、FAPSOMA算法更佳的完工时间，FAPSOMA算法相比FA算法所得到的完工时间更佳。对于同一规模，单从十次结果的平均值来看，例如100*5规模，HFPMCV算法的平均值为5 507，FAPSOMA的平均值为5 565.9，FA的平均值为5 628.4，HFPMCV算法的平均值小于FA、FAPSOMA算法，FAPSOMA算法相比FA算法得到的平均值更小。综合两个指标来看，HFPMCV算法在求解流水车间调度问题方面与FA、FAPSOMA算法相比，性能更优、强度更高、求解结果更加精准。

甘特图，别名条状图，较为简单、直观，以图示的方式来向人们展示工程进度安排，甘特图上面标识着工件加工顺序，每一个工件在每一台机器上的开始加工时间、结束加工时间，整个优化方案的用时等。通过生成的甘特图可以直观的看到算法求解流水车间调度的结果，即生成的最优加工顺序以及最小完工时间。

图3展示了HFPMCV算法求解一次20*5的车间调度问题时所得的甘特图，X轴表示加工时间，Y轴表示机器号，J代表工件号，st代表工件开始加工时间，ft代表工件结束加工时间，如工件8在第5台机器上面的开始加工时间为138，结束加工时间为207，加工用时69，所有的工件在第一台机器上全部加工完成的时间为1 121，所有的工件在第五台机器上全部加工完成的时间为1 305，从图中可以看出工件的加工顺序为[8,7,16,14,5,13,10,11,1,2,15,12,4,17,3,0,18,9,6,19]，最小完工时间为1 305。

图3 HFPMCV算法求解车间调度问题的甘特图

3 结束语

本文针对流水车间调度模型求解问题，并考虑到单种算法对全局最优解的收敛局限性，并且结合前人提出的萤火虫粒子群混合算法，在萤火虫粒子群混合算法的基础上，提出了HFPMCV算法，通过引入多尺度协同变异算子、混沌初始化、动态自适应策略以及平行进化模式，提高了算法的求解速率和精度，函数优化实验验证了HFPMCV算法的可行、有效，用HFPMCV、FA、FAPSOMA算法解决以最小化最大完工时限为指标的流水车间调度问题，对比实验结果显示HFPMCV算法得到的完工时间均小于FA、FAPSOMA算法，说明了HFPMCV算法速率更快、精度更高，与基本算法相比较，在求解流水车间调度模型时，更具有一定的实用价值。但萤火虫位置更新公式仍需要优化，以改善算法时间复杂度，获得更优的求解结果。接下来将着重于研究HFPMCV算法在其他车间调度问题中的应用。