改进樽海鞘群算法求解柔性作业车间调度问题

赵文超，郭鹏,2，王海波,2，雷坤

（1.西南交通大学 机械工程学院,四川 成都 610031;2.轨道交通运维技术与装备四川省重点实验室,四川 成都 610031）

先进的调度方法能够提高生产效率和经济效益，其主要任务是在时间和资源的双重约束下，合理、科学地对可用共享资源和生产任务进行分配与管理，尽可能使性能评价指标达到最优。柔性作业车间调度问题(flexible job shop scheduling problem，FJSP)作为车间作业调度的延伸，突破了一道工序只能在一台设备加工的局限，引入了设备指派决策，扩大了可行域的搜索范围，具有更高的理论计算复杂度。FJSP 广泛应用于航空航天、交通运输、智能制造和运筹优化等领域，因此针对FJSP 的研究具有重要理论价值和实际应用价值。

针对柔性作业车间调度优化问题，目前主要的求解算法可以分为两大类：精确算法和近似算法[1]。由于问题的高度复杂性，研究者大多采用智能优化算法求解，诸如模拟退火与禁忌搜索混合算法[2]、遗传-禁忌搜索混合算法[3]、改进迭代贪婪算法[4]等。近年来随着计算机科学技术的不断进步，许多新的群体智能算法不断涌现，如人工免疫算法[5]、混合灰狼优化算法[6]等。也有使用分派规则构建启发式算法的文献[7]。此外针对FJSP的扩展问题，先后有文献讨论顺序相关调整时间[8]、能耗约束[9]、动态事件[10-11]、多时间约束[12]等变种。尽管出现了大量启发式求解策略，但求解效率依然有进一步提升的空间和基于群集智能发展高效调度框架的必要。

樽海鞘群算法在2017 年由澳大利亚学者Mirjalili等[13]提出，该算法通过模拟自然界中樽海鞘群在海洋中的觅食行为，实现对解空间的探索。该算法具有收敛速度快，鲁棒性强且易于实现等显著特点，已成功地应用于光伏系统优化[14]、特征提取[15]、图像处理[16]和生物医学信号处理[17]等问题求解。Jia 等[18]提出了SSACL 算法，使用交叉策略和Lévy 飞行分别改进了樽海鞘群领导者和跟随者的运动方式，在基准函数测试中表现出良好的计算性能。然而，该算法在生产调度领域的研究并不多见，且现有的文献多集中于并行机调度问题的研究。Jouhari 等[19]针对不相关的并行机调度问题，利用樽海鞘群算法进行局部搜索以增强优化器的性能并减少计算时间。Ewees 等[20]提出了一种基于萤火虫算法的改进樽海鞘群算法来求解带调整时间不相关的并行机调度问题。Sun 等[21]提出了将樽海鞘群算法应用于车间调度问题，结合樽海鞘群算法全局搜索框架，基于关键路径设计多个邻域算子进行局部搜索，对可重入作业车间调度问题进行了有效的求解。本文利用Lévy飞行和交叉变异算子对樽海鞘群算法进行改进，并结合局部邻域搜索策略来求解柔性作业车间调度问题。鉴于樽海鞘群算法的优越性，非常有必要研究其在生产调度领域的应用。

本文提出一种融合模拟退火策略的改进樽海鞘群算法求解FJSP。基于设备负载进行种群初始化，利用交叉和变异遗传算子对个体进行改进，采用Lévy 飞行对领导者个体位置进行更新。为防止求解陷入局部最优，融合模拟退火进行邻域搜索，提升算法的全局搜索和局部搜索能力。同时，为验证模型的准确性并为所提算法提供参考依据，采用Gurobi 数学求解器进行精确求解混合整数规划模型。运用改进的樽海鞘群算法求解FJSP，并和其他研究文献进行综合对比，进一步阐明所提算法的有效性。

1 问题描述与模型构建

1.1 问题描述

FJSP 描述如下：给出n个工件的集合J和m台设备的集合M。每个工件i的加工工序都是确定的且不尽相同，其中工件i的第j道工序(Oij)可以在兼容设备子集Mij⊆M的任何设备上加工。每道工序的加工时间与所选择的加工设备有关，每个工件都必须严格按照预定好的顺序进行加工，且每台设备在同一时刻只能加工一个工件。调度的目标是让每道工序选择最合适的加工设备，并确定每台设备上的最佳处理顺序和加工开始结束时间。因此，该组合优化问题分为两个子问题：1) 设备选择问题，即为工序选择合适的加工设备；2)工序排列问题，即确定每台操作设备的最佳加工工序序列，以此确定两个子问题的最优值来使得整个调度系统实现特定的优化目标。本文考虑最小化最大完工时间Cmax，即makespan。

此外，在建模过程中还需满足如下约束条件：

1)工件在到达时间之前不允许加工，所有工件在零时刻都可以被加工；

2)同一工件的工序有加工顺序约束，不同的工件没有加工顺序约束；

3)同一台设备在同一时刻只能加工一个工件；

4)同一工件的同一道工序在同一时刻只能被一台设备加工；

5)工件在某台设备上开始加工，不允许中断；

6)所有加工设备是连续可用的，设备的调试与设置时间可以忽略不计。

1.2 模型构建

针对柔性作业车间调度问题，以最小化最大完工时间Cmax为目标，基于文献[8]的模型构建本问题的整数规划(MIP)模型，各参数定义见表1。

表1 符号与参数Table1 Symbols and parameters

目标函数(1)最小化最大完工时间makespan。式(2)确保每道工序当且仅能指派给可用设备集中的一台进行加工。式(3)保证同一工件相邻工序间的先后关系。式(4)和式(5)防止在同一设备k上进行处理的工序作业时间发生重叠，当且仅当工序Oij和工序Oi′j′都分配给设备k(即zijk=zi′j′k=1)时，两个约束方才起作用。式(4)确保当工序Oi′j′完工后处理工序Oij，则yiji′j′=0。式(5) 表示yiji′j′=1 的情况。对于其他情况，式(4)和式(5)则自动满足。式(6)确定最大完工时间。式(7)和(8)决定两组0、1 变量的取值。

2 樽海鞘群算法

樽海鞘群算法模拟樽海鞘的聚集行为，组成樽海鞘链在海底进行捕食和移动。在SSA 中，樽海鞘链可以分为两类：领导者和追随者，链的头部是领导者，其余的都是跟随者。与其他群体优化算法类似，樽海鞘群的位置定义在多维搜索空间内，寻找多维空间中食物源F的位置是樽海鞘群追随的重要指标。领导者根据F的位置修正自身的位置，追随者跟随领导者位置变化更新自身的位置。

在樽海鞘群算法中，定义每个樽海鞘个体的位置向量X用于在N维空间中搜索，其中N为决策变量的数目。樽海鞘群算法中位置向量X将由维度为D的N个樽海鞘个体组成，其中第i个樽海鞘个体位置表示为樽海鞘个体位置为

因此，种群向量由N×D维矩阵构成，即：

矩阵中的第一个向量为领导者，其余向量为追随者。食物源F的位置是所有樽海鞘个体的目标位置。因此，领导者的位置更新公式为

由式(9)可知，领导者的位置更新与食物源的位置相关。系数c1叫做收敛因子，其定义如下：

式中：t为当前迭代次数；tmax为 最大迭代次数。

参数c2和c3是 [0,1] 区间内均匀生成的随机数，c2决定了在第j维空间领导者更新的移动步长，c3决定的是移动的方向的正反。

根据牛顿运动定律，对跟随者的位置更新进行公式化：

式中：i≥2，表示第j维空间跟随者的位置；v0表示初始速度；t0表示时间；a=(vfinal−v0)/t，vfinal=表示第i−1个追随者在第j维空间跟随者的位置。

由于时间在优化过程中表示为迭代，所以t0=1迭代差值为1，初始速度v0=0，并且追随者的位置更新只与它前一个樽海鞘个体位置相关，所以位置更新公式进而可以表示为

式中：i≥2，表示第i个樽海鞘个体在第j维空间的位置。根据以上个体位置更新方式，可以模拟樽海鞘群的行为机制。

3 改进樽海鞘群算法

自提出以来，SSA 算法在求解连续优化问题得到了广泛的应用，但在组合优化领域的应用较少。因此，针对FJSP 的特点，提出一种基于离散思想的SSA，引入基于设备负载最小的选择方案来加快算法收敛；采用交叉变异算子调整工序与设备编码，以此扩大解的搜索空间从而加快算法的搜索效率；防止求解陷入局部最优或者无法完全收敛，融合模拟退火策略进行局部领域搜索，进而得到改进离散化樽海鞘群算法(improved discretized salp swarm algorithm，IDSSA)求解FJSP。

3.1 编码

针对FJSP 的特性，IDSSA 编码由设备选择部分(MS)和工序排序部分(OS)两个向量组成，每个向量的长度均等于总工序数之和。设备选择部分从第一个工件的第一道加工工序开始直到最后一个工件的最后一道加工工序依次分配加工设备，数字k为对应工序可选设备集内的第k台设备；工序排序部分中的每一维向量代表待加工工件序号，其出现的次数代表工件的第几道工序。(O11,M3),(O21,M4),(O22,M4),(O12,M1),(O23,M3),(O13,M2),(O31,M3),(O32,M1)是利用二维向量编码的个体。个体的编码方式如图1 所示。

图1 IDSSA 编码示意Fig.1 IDSSA code representation

3.2 种群初始化

种群初始化是樽海鞘群算法的重要步骤，初始解的质量对SSA 的收敛速度和寻优效率至关重要。在生成初始解时，不仅要保证解决方案的多样性而且还要保证解的质量。到目前为止，大多数文献采用随机初始化的方法初始化加工设备，无法保证最短时间的加工设备被选择。因此，最初未选择的设备将永远不会被选择，这使得算法仅赖于设备突变操作来选择这些设备，但突变操作只有很小的发生概率，初始解的多样性也就无法得到保证。FJSP 设备选择的难点在于每一道工序选择加工设备时，对设备而言，设备的负载变化会对加工过程中设备选择的结果产生影响，因此需要考虑加工该工序后的设备负载变化。使用SSA 进行初始化时，由于没有任何先验知识可以参考，本文参考Kacem 等[22]提出的初始种群定位法，该方法在考虑设备的负载前提下为工序分配设备，为每一道工序的可选设备集中寻找时间表中具有最短加工时间的设备，利用该设备处理工序。然后对设备负载进行更新操作，即将前一道工序的加工时间加到同一列的其他项。

3.3 位置更新方式

SSA 中领导者的位置更新方式适用于求解连续函数优化问题，无法直接用于离散问题，因此需要对更新策略进行离散化改进：

式中：Fj表示第j维食物源的位置；系数向量c1是SSA 中最重要的参数，用于平衡算法的局部搜索和全局探索，随着迭代次数的增加逐渐变化，最后趋近于0；r为区间[0，1]的随机数，当r<0.5，樽海鞘群个体xi j向食物源Fj趋近，相反则远离。L为Lévy 飞行步长，可有效避免陷入局部最优，更加快速找到全局最优解，Lévy 飞行步长定义为

式中：µ和υ服从正态分布，其中的 σµ和συ由式(16)计算可得：

式中，Γ是gamma 函数，参数 β是区间[0，2]的随机数，一般 β=1.5。追随者的位置更新策略是让其有针对性的移动，找到更好的适应度位置，从而加快最优解的搜索。追随者的移动是由自身位置和前一个个体的位置综合决定，对前一个个体位置有较强的依赖性。若追随者位置陷入局部最优，容易造成算法搜索停滞。为了更好地平衡算法的开发能力和搜索能力，引入线性递减的惯性权重，平衡先前个体对当前个体的影响。因此，引入惯性权重的跟随者位置更新公式为

式中：ω(s)为初始惯性权重；ω(e)为最大迭代次数时的惯性权重。初始迭代时较大的惯性权重有助于提升算法的搜索能力，迭代后期惯性权重较小时，有助于提升算法的开发能力。

SSA 算法中每个个体经过离散化处理后，利用ROV(ranked order value)规则，分别对每个个体的工序位置向量赋予唯一的ROV 值，然后根据ROV 值即可构造工序编码。个体位置转换为工序排序过程如图2 所示，其中相同的工序编号表示同工件的不同工序，出现的次数代表工件的第几道工序。

图2 个体位置转换为工序排序过程Fig.2 Process of converting individual positions into process sequencing

设备编码个体的位置向量限制在[1，m] 内，其中m为加工设备数，若计算得到的向量里元素值超过此区间，则取边界值。经过位置更新后，设备位置向量为小数，对其进行向上取整。位置更新后产生的设备编码可能为非可行解，需要对非可行解进行修正：1)设备向量对应的工序的可选加工设备集中仅有一台设备(即该工序仅能在一台设备上加工)，则选择该工序对应的设备序号更新设备位置向量；2)位置更新后的某道工序选择的设备不能加工该工序，则在该工序的可选加工集中随机选择一台设备替换该设备。

3.4 交叉操作

考虑到遗传算法有比较强的全局搜索能力[23]，融合交叉和变异算子对标准SSA 算法进行改进，以增加种群的多样性并加快收敛速度。SSA 算法不依赖遗传算子进化而是通过个体的位置移动来寻找问题最优个体。引入POX 交叉算子(如图3所示)在OS 部分操作方式如下：

图3 工序编码的交叉操作Fig.3 Cross operation of process code

1)生成一个从1 到工件数的随机整数R；

2)父代个体P1将工件序号小于或等于R的工件复制给子代个体C1；父代个体P2将工件序号大于R的工件复制给子代个体C2，保留其位置，并将对应的设备号复制到相应位置；

3)复制父代P1未出现在子代C2的工件序号到C2，复制父代P2未出现在C1的工件序号到子代C1，并保留其顺序，同时复制设备号到对应位置。

位置更新过程中从子代C1和子代C2中随机选择一个作为后代，POX 交叉过程使子代保留父代在每台设备的加工次序，并保留部分工件的位置。

交叉算子在MS 部分操作如下：针对设备选择部分，采用两点交叉的方法，对于选定的两个父代个体，随机设置两个交叉点，交换所设定的两个交叉点之间的父代个体都有的工序所对应的设备序号。

3.5 变异操作

变异操作目的是为了增加种群多样性，改善算法的寻优能力。FJSP 的变异操作后需要保证变异后解的可行性。和交叉操作一样，分别基于OS 和MS 进行变异操作。设计变异算子时，引入关键路径的思想。FJSP 的一个可行解可用有向图来表示，其中从有向图起点到终点的最长路径称为关键路径。关键路径直接影响调度方案的最大完工时间。关键路径上的关键工序也称作关键工序，通过对关键工序进行扰动，有可能会减少最大完工时间。图4 为图1 对应的甘特图，其中阴影部分为关键工序。

图4 基于关键路径甘特图Fig.4 Gantt chart based on the critical path

在满足加工优先级约束前提下，寻找关键路径上的所有关键块，随机选择关键块为移动源，选择不相邻的工序为移动点，向前或向后插入移动源，其余工序依次向前或向后改变位置，如图5 所示。

图5 关键工序邻域结构示意Fig.5 Diagram of neighborhood structure of key processes

对OS 部分的变异操作引入基于关键路径的变异算子，将变异位置的选择缩短到关键路径上。具体操作过程为随机选择一道关键工序，满足工序加工优先级约束的前提下，将它插入到紧邻的前一道关键工序的某个位置，同时将相应的设备分配同步前移或者后移。图6 为选中关键工序O31的情况，并将其插入在另外一道关键工序O22之前。

图6 工序编码的关键工序变异操作Fig.6 Key process mutation operation for process coding

对MS 部分的变异操作为随机选择一道关键工序，在选择的关键工序可行加工设备集中选择加工时间最少的设备替换当前加工设备。这种选择方式增加种群的多样性，扩大解的搜索范围。

3.6 局部邻域搜索策略

采用上述编码和位置更新方式，有效地提高算法求解效率和收敛速度，但同时容易陷入局部最优解或无法完全收敛。SA 算法具有较强的邻域搜索能力，将SA 算法与IDSSA 算法相结合可以有效提升算法的全局搜索和局部搜索能力，进而提升算法性能。本文基于设备负载的变异方式产生邻域解，具体流程如下：寻找所有加工设备中工作负载最大的设备，并在设备对应加工工序中随机选择一道工序，在考虑工序的优先级约束前提下将其分配到可选设备集中负载最小的设备。融合SA 算法只对MS 向量进行调整，让OS向量匹配到更合适的设备，从而平衡了设备工作负载。但同时会增加搜索时间，因此本文在初始种群中随机选择部分个体进行局部邻域搜索，增加种群的多样性。

3.7 IDSSA 流程

IDSSA 算法流程具体步骤如下：

1)初始化参数，根据3.2 节的初始定位法生成初始解，并随机化食物源F的位置。

2)计算个体的目标函数值，将算法最优适应值为迭代过程中最佳的计算结果。

3)根据改进的位置更新公式分别对樽海鞘领导者和跟随者的位置进行更新。引入惯性权重对追随者个体位置进行优化，并记录当前种群中最优食物源的位置。

4)利用交叉、变异算子分别对个体的工序向量和设备向量进行交叉操作，增加可行解的搜索范围。

5)更新算法相关参数。

6)判断是否达到最优目标值，若新解的目标函数值优于当前解，则更新最优解并将其输出，否则采用SA 算法，随机选择部分个体进行局部邻域搜索。

7)判断是否满足终止条件，即是否达到最大迭代次数。若满足，则跳转到8)，反之，执行2)。

8)算法结束。

4 计算结果与分析

为了测试IDSSA 算法在解决全局优化问题中的效果，将IDSSA 算法和Jia 等[18]提出的SSACL算法、标准SSA 算法进行计算对比。IDSSA 算法利用Lévy 飞行改进领导者更新方式，在追随者更新公式上引入自适应惯性权重，利用交叉、变异策略增加种群多样性，并加入邻域搜索进一步改善算法的局部寻优能力。

采用CEC2005 的10个基准测试函数，选取的测试函数包含单峰和多峰两种类型函数，其中单峰函数在定义上下限区间内只有一个全局最优，因此可以测试算法的开拓能力；多峰函数在定义区间含有多个全局最优和局部最优解，可以测试算法的全局搜索能力。函数的维度、定义域、理论最优解和类型如表2 所示。

表2 基准函数Table2 Benchmark functions

为了对比的公平性，将算法的参数和SSACL算法设置一致：种群大小为30，迭代次数为500。为了避免随机误差的影响，每个测试函数独立运行30 次。计算对比数据如表3 所示，表中每个测试函数的适应度均值和标准差分别反映了不同算法的收敛精度和稳定性。

对于表3 中的测试函数，SSACL 算法和IDSSA算法比标准SSA 算法寻优结果都要好。在求解精度上，IDSSA 算法有4个基准函数计算结果表现最好；就算法测试函数的平均值和标准差而言，IDSSA 算法也表现出良好的计算优势，表明对标准SSA 算法进行改进有效地提升了算法性能。

表3 基准函数优化结果对比Table3 Comparison of optimization results of benchmark functions

在验证IDSSA 求解FJSP 的改进效果和有效性方面，本文采用国际通用的10个Brandimarte基准算例、Fattahi 提出的20个基准算例进行对比分析。Brandimarte 算例中，工件数量n从10 到20，设备数量m从4 到15 进行组合选取，每组算例中的工序数从5 到15；Fattahi 算例中，工件数量n从2 到12，设备数量m从2 到8 进行组合选取，每组算例中的工序数从2 到4。使用Python3.6实现算法。计算机硬件配置为英特尔i56400 CPU(2.70 GHz)、8 GB RAM，在Windows 10 的操作系统下运行程序。

4.1 算法参数选取

参数设置对算法性能有很大的影响，参数选取标准是通过解的质量和算法运行时间来衡量，文中利用Brandimate 提出的算例Mk01 结果来确定算法参数。以解的平均偏差和平均运行时间为基准，经过多次计算从而确定相关参数。本文所提的IDSSA 主要包括 β、T0、µ、ω(s)、ω(e)和Tlim等6个参数。基于计算结果，参数 β设置为1.5；初始温度T0和冷却系数 µ分别被设置为100 和0.95；利用Gurobi 求解限制时间为Tlim=1 200 s；经过正交实验后确定 ω(s)=0.9，ω(e)=0.4 时算法具有最佳性能。

4.2 结果分析

为了更好评估融合交叉与变异算子的IDSSA算法的有效性和稳定性，选取Brandimate 提出的算例Mk01 进行测算。收敛曲线如图7(a)所示，通过观察IDSSA 算法与不含交叉变异操作的改进SSA 算法的收敛曲线，发现引入交叉与变异算子的I D S S A 算法收敛曲线下降更快。较之SSA 算法，IDSSA 算法收敛速度和稳定性得到进一步改善。同时，为验证局部邻域搜索策略的有效性，同样选取算例Mk01 进行测算。收敛曲线如图7(b)所示，将IDSSA 算法与不含局部搜索策略的改进SSA 算法进行对比发现，IDSSA 算法表现出了更高的搜索精度。

图7 IDSSA 与SSA 的收敛曲线比较Fig.7 Comparison of IDSSA and SSA convergence curves

为验证所提MIP 模型的正确性并给所提算法运算结果提供对比参考，采用Brandimarte 基准算例[24]进行测试，运用Gurobi 8.1.1 进行精确求解。为了达到性能评价的目的，在本文所提算法和HGWO 算法[6]、Heuristic 算法[7]以及后续对比算法的参数设置上选择了最佳的、相同的参数配置以保证了算法的可比性。由于问题的难求解性，处理MIP 模型时将Gurobi 时间限定为1 200 s，若求解时间未超过1 200 s，则输出解为最优解。除此之外，还采用Fattahi 基准算例进行对比，并与文献中的算法进行对比，包括HTS/SA 算法[2]、MIG算法[4]、AIA 算法[5]。

比较本文提出算法与其他文献算法相较当前最优解的相对百分偏差RPD(relative percent deviation)，计算公式为

其中对于每个算例，UB表示当前算例目前最佳目标函数值，Cmax表示当前算例求解的最大完工时间。

表4 为Brandimarte 基准算例测试结果。表4中n×m表示对应算例的工件数和设备数。ARPD表示当前算法的相对百分偏差平均值(%)。CPU表示算法运行时间，s。为了排除误差的影响，更加准确地对比所提算法的可行性和有效性，将每个算例运行10 次。可以看出所提出的MIP 在求解较小规模算例得到了令人满意的解，在求解较大规模算例陷入局部最优，不能在有效的时间内得到较优解。本文所提出的IDSSA 和其他参考文献提出的算法相比，计算结果中9个算例优于Heuristic，相比HGWO 存在较大的优势。相较当前最优解的平均偏差来看，IDSSA 和其他算法相比偏差值最小，即从整体上来看，所提出的IDSSA 优于MIP 和参考文献[6-7]所提出的其他两种算法。

表4 Brandimate 算例计算结果对比Table4 Comparison of calculation results of Brandimate instances

续表4

表5 给出了本文算法在Fattahi 测试集上的性能，对比方法有HTS/SA[2]、AIA[4]、MIG[5]等。在较小规模算例上，表中列举算法求解性能差异不大，而在较大规模算例(即MFJS9 和MFJS10)上，本文算法显示出良好的求解性能。图8 为IDSSA算法求解问题MFJS9 输出结果的甘特图。

图8 MFJS9 调度甘特图Fig.8 MFJS9 scheduling Gantt chart

表5 Fattahi 算例计算结果对比Table5 Comparison of calculation results of Fattahi instances

5 结束语

本文在标准樽海鞘群算法的基础上，提出了一种融合模拟退火算法的改进樽海鞘群算法。基于Lévy 飞行对领导者的位置更新方式进行离散化改进，引入自适应惯性权重对追随者的位置进行移动，平衡算法的开发能力和搜索能力。建立问题的混合整数规划模型并利用标准算例检验算法的性能，结果表明：改进的樽海鞘群算法在求解单目标柔性业车间问题取得了比较好的结果，算法的鲁棒性和有效性得到了验证。同时，本文为连续化SSA 算法进行离散化求解工程实际问题提供了一种可行的改进方式。

下一步还将开展算法的计算时间复杂度理论分析及算法的收敛特性理论证明相关工作。同时，深入分析FJSP 的问题特征，在实际的生产作业中，存在着物料组成、工人技能、设备资源受限等诸多约束条件，考虑将改进的樽海鞘群算法应用到更加复杂的工程实践问题当中，进一步验证算法的性能。