基于混合化学反应算法的柔性作业车间调度

肖华军，柴子力，张超勇，孟磊磊，任亚平，梅慧文

(华中科技大学 数字制造装备与技术国家重点实验室，湖北 武汉 430074)

0 引言

调度问题通常定义为“把有限的资源在时间上分配给若干个任务，以满足或优化一个或多个目标”。在作业车间调度问题(Job-Shop Scheduling Problem, JSP)中，要加工的工序都必须在指定的机器上加工，然而在实际生产中要加工的工序往往可以在多台机器上进行加工，而且不同的机器加工同一工序所用时间不同，这一类问题称为柔性作业车间调度问题(Flexible Job Shop Scheduling Problem, FJSP)。FJSP需要为工序选择可加工机器，减少机器约束的同时，却会增大搜索空间，是比JSP问题更复杂的NP-hard问题[1]。

自Brucher等[2]在1990年首次提出该问题以来，已经出现很多算法用来解决该问题。求解FJSP的算法主要包括精确算法和近似算法两大类。精确算法包括分支定界法、整数规划法等。Torabietal等[3]在2005年提出了一种整数规划方法用于求解确定性柔性作业车间中的共同周期多产品批量调度问题，需要同时确定机器分配，排序，批量调整和调度决策。Demir等[4]在2013年评估了几种FJSP的数学模型。同年Roshanaei等[5]进一步提出了两种新颖的混合整数线性规划(Mixed Integr Linear Programming, MILP)模型用于解决FJSP。精确算法可以获得精确的优化方案，然而不能有效地解决总工序较多的大规模FJSP。近似方法包括启发式算法和元启发式算法[6]。进化算法是一种有效的元启发式方法，包括遗传算法，进化策略和进化规划。Tay等[7]在2008年将遗传算法与优先分配规则结合使用来解决多目标的FJSP。2014年，赵诗奎等[8]将基于机器空闲时间的邻域搜索策略融入遗传算法以求解作业车间调度优化问题。2016年，张国辉等[9]针对FJSP的NP难特性提出解阈值的指标，采用精英进化策略对其进行求解。群体智能算法是另外一种高效的元启发式算法，包括蚁群优化(Ant Colony Optimization, ACO)算法、粒子群优化(Particle Swarm Optimization, PSO)算法和人工蜂群(Artificial Bee Colony, ABC)算法等。2009年Girish等[10]提出了解决以最大完工时间为目标的FJSP的ACO算法，并且评估了算法的性能。Wang等[11]在2010年提出了求解FJSP的算法。2012年，张静等[12]提出一种新的粒子位置更新方式以改进粒子群算法，使之更有效地求解柔性作业车间批量调度问题。2015年，仲于江等[13]将小生境技术和粒子群算法相结合，求解多目标的FJSP优化问题。群体智能算法全局搜索能力较强，但是局部搜索能力较差，因此将群体智能算法与局部搜索算法结合起来能优势互补，获得更好的搜索能力。局部搜索算法包括禁忌搜索(Tabu Search, TS)算法[14]，模拟退火(Simulated Annealing, SA)算法[15]和变邻域搜索(Variable Neighborhood Search, VNS)算法[16]等。

化学反应优化(Chemical-Reaction Optimization, CRO)算法是一种模仿分子在化学反应中变化情况的群体优化智能算法[17]。CRO算法提出势能、动能和能量中心等概念，避免算法陷入早熟。CRO算法在求解离散优化问题，如二次分配[18]、模糊作业车间调度[19]和网格调度[20]等方面得到有效的应用，并展现了良好的性能。本文针对FJSP的特点提出一种融合CRO算法与TS算法结合的混合算法(CROTS)。混合算法采用有效的编码、解码、初始解产生机制和精英保留策略，设计了CRO算法的4种基本反应；通过交叉算子实现了分子间的碰撞反应，完成全局搜索；通过变异算子实现分子与墙壁的碰撞反应，并释放动能到能量中心；在满足分解反应条件时，将单个分子分解成多个分子；在满足合成反应条件时将多个分子合成单个分子，以此增加分子群的多样性。最后再利用禁忌搜索，让分子群跳出全局最优，使混合算法在分散搜索和集中搜索之间达到更合理的平衡。

1 FJSP问题描述

FJSP模型建立如下:n个工件J={J1,J2,…,Jn}要在m台机器M={M1,M2,…,Mm}上完成加工。任意工件Ji包含ni道工序{Oi,1,Oi,2,…,Oi,ni},同一工件上的不同工序必须按照确定的先后顺序依次进行加工。每道工序的可加工机器集为Mi,j={Mi,j,1,Mi,j,2,…,Mi,j,k}，Mi,j⊆M。记工序Oi,j在机器Mi,j,k上加工所需时间为pi,j,k(pi,j,k>0)。调度的任务是确定所有工序的加工顺序和所选用的机器。本文考虑的优化目标为寻找一个可行调度方案使得最大完工时间(Makespan)最小。设Ci是工件Ji的完工时间，则最大完工时间Cmax最小的目标函数为：Cmax=min{max(Ci)},(i=1,2,…,n)。

为了描述方便，定义以下符号：

Ci为工件Ji的完工时间;

n为工件总数;

m为机器总数;

M为总机器集;

ni为工件i的工序总数;

Oi,j为工件i的第j道工序;

Mi,j为工序Oi,j的可选加工机器集;

mi,j为工序Oi,j的可选加工机器数;

Mi,j,k为工序Oi,j在第k台机器上加工;

pi,j,k为工序Oi,j在第k台机器上加工所需时间;

Si,j为工序Oi,j的加工开始时间;

Ci,j为工序Oi,j的加工完成时间;

Ci为工件i的完工时间;

Cmax为最大完工时间;

L为一个足够大的正数。

FJSP模型描述如下：

目标函数Cmax=min{max(Ci)},i=1,2,…,n。

(1)

s.t.

Si,j+xi,j,k×pi,j,k≤Ci,j，

i=1,…,n,j=1,…,ni,k=1,…,mi,j;

(2)

Ci,j≤Si,j+1，i=1,…,n,j=1,…,ni-1;

(3)

Ci,ni≤Cmax，i=1,…,n;

(4)

Sg,h+pg,h,k≤Si,j+L(1-yghijk),

g=1,…,n,h=1,…,ng,

i=1,…,n,j=1,…,ni,k=1,…,m;

(5)

(6)

Si,j≥0,Ci,j≥0，i=1,…,n,j=1,…,ni。

(7)

式(2)和式(3)表示同一工件上的工序先后顺序约束，式(4)表示工件的完工时间小于或等于总完工时间，式(5)表示在同一机器上加工的工序先后顺序约束，式(6)表示每个工序只能在一台机器上加工，式(7)表示开始加工时间和完工时间都得为正数。

表1为n=3,m=4的一个柔性作业车间调度问题时间表。

表1 柔性作业车间调度问题加工时间表

2 化学反应优化算法

CRO算法由LAM等[17]于2010年首次提出,是一种启发式优化算法，它模拟了密闭容器中分子间的碰撞情况。该算法中分子对应问题的解，基元反应对应搜索过程。

分子是执行算法操作的个体，每个分子包含分子结构，势能(PE)和动能(KE)。该算法用分子结构ω来刻画问题的解；势能对应当前分子结构所具有的目标函数值；动能则体现了分子发生反应的能力，动能越大分子发生碰撞的概率越大。本文设定的势能为分子的目标函数值，动能为分子的初始动能InitialKE减去父代分子的势能与子代分子势能的差值。则：

PEω=ObjFunc(ω),

(7)

KEω=InitialKE-(PEω-PEω′)。

(8)

式中ObjFunc指目标函数，基元反应对应CRO算法操作算子。CRO算法定义了分子与墙壁的碰撞反应、分子间的碰撞反应、合成反应以及分解反应4种基本操作。其中：分子与墙壁的碰撞反应和分解反应在算法中为对单个个体的扰动；分子间碰撞反应和合成反应在算法中对应多个个体的相互作用。CRO算法流程图如图1所示。

3 求解FJSP的混合算法

3.1 编码和解码

FJSP在给加工工序排序的同时还要给工序分配可加工机器，因此问题的可行解需要确定工序先后顺序的编码和给工序分配机器的编码两部分，如图2所示。

基于机器的编码长度也为l，表示为[g1,g2,…,gl]，其中第i个元素gi为第i道工序的加工机器。因此图2所示的分子代表的工序和机器序列为(O2,1,M1)，(O1,1,M1)，(O3,1,M2)，(O1,2,M3)，(O2,2,M1)，(O1,3,M4)，(O3,2,M2)，(O3,3,M4)，加工时间序列为[3 2 1 4 4 2 2 3]。

解码是编码的逆过程，是将分子转换成调度解的过程。设r为工序编码中任意一工序Oi,j，加工时间为pr。记与工序r在同一工件上的前一道工序为PJ[r],后一道工序为SJ[r](如果它们存在)。记与工序r在同一台机器上加工的前一道工序为PM[r]，后一道工序为SM[r](如果它们存在)。设工序r的最早开始加工时间为SE[r]，最早完工时间为CE[r]，则：

SE[r]=max{CE[PJ[r]],CE[PM[r]]}，

(9)

CE[r]=SE[r]+pr。

(10)

解码之后得到问题的半主动调度解，其甘特图如图3所示。

由图3可知，该解的最大完工时间为makespan=14，若将工序O3,3插入到O1,3之前，而不改变其加工机器，所得的解依然是可行解而且减少了最大完工时间，此时makespan=11优化了目标。因此本文引入插入式贪婪解码算法来获取主动调度解。该解码方法描述如下：在满足约束条件且不干扰其他工序的情况下将当前工序安排在可加工机器上的最早的空闲时间段上。由此解码得到的主动调度解的甘特图如图4所示。

3.2 种群初始化

为保证种群的多样性，本文采用多种初始化方法来产生初始种群。第一种方法采用Kacem等[21]提出的局部化方法，种群中20%的个体通过该种方法产生;第二种初始化方法在给工序选择加工机器时优先选择加工时间短的机器，该过程通过二元选择操作来完成，而工序的加工顺序随机产生即可，种群中80%的个体通过该方法产生。

3.3 交叉算子

分子之间的碰撞是通过两个分子的碰撞反应生成两个新的子代分子的过程。结合FJSP的特征，本文采用交叉操作来实现分子间的碰撞反应。针对FJSP，分子中包含两条编码，对于这两部分编码本文分别采用不同的交叉算子。对基于工序的编码，本文采用优先工序交叉(Precedence Operation Crossover, POX)操作和改进优先工序交叉(Improved Precedence Operation Crossover, IPOX)操作。在交叉过程中各按50%的概率随机选择。两种交叉操作介绍如下。

(1)POX交叉操作

设P1和P2为FJSP的两个父代，通过POX交叉产生两个子代C1和C2。交叉步骤如下：

步骤1将工件集J={J1,J2,…,Jn}随机分成两个集合G1和G2。

步骤2父代P1工序编码中属于集合G1的元素直接保留到C1中，且保持原有位置不变；同样，父代P2工序编码中属于集合G1的元素直接保留到C2中，且保持原有位置不变。

步骤3父代P2工序编码中属于集合G2的元素按顺序依次填补到C1中的空白处，同样父代P1工序编码中属于集合G2的元素按顺序依次填补到C2中的空白处。

POX交叉操作过程如图5所示。

(2)IPOX交叉操作

记P1和P2为分子群中的两个父代分子，C1和C2为父代分子经由IPOX交叉得到的子代分子。IPOX交叉步骤如下：

步骤1将工件集J={J1,J2,…,Jn}随机分成两个集合G1和G2。

步骤2父代P1工序编码中属于集合G1的元素直接保留到C1中，且保持原有位置不变；同样，父代P2工序编码中属于集合G2的元素直接保留到C2中，且保持原有位置不变。

步骤3父代P2工序编码中属于集合G2的元素按顺序依次填补到C1中的空白处，同样父代P1工序编码中属于集合G1的元素按顺序依次填补到C2中的空白处。

IPOX交叉操作过程如图6所示。

针对基于机器的编码，本文采用多点交叉操作(Multi-point Preservative crossover, MPX)，即多点交叉操作。MPX交叉步骤如下：

步骤1产生随机0,1数组R，数组长度为编码长度l。

步骤2选出P1和P2中与数组R中元素1位置相同的元素，将它们互换。

步骤3P1和P2中其他元素保持不变。

由于交换的是同一道工序的加工机器，保证了生成的子代分子依然是可行解。MPX交叉过程如图7所示。

3.4 变异操作

分子与墙壁之间的碰撞是一个分子通过与墙壁之间的碰撞反应得到一个与其相近的子代分子的过程。结合FJSP问题特征，本文采用变异操作来实现分子与墙壁之间的碰撞。对于基于工序的编码，本文采用两种变异操作，在变异过程中各按50%的概率随机选择。

(1)随机插入变异操作，如图8所示，r1=5，r2=3。具体步骤如下：

步骤1产生两随机数r1和r2，r1∈[0,l],r2∈[0,l]且r1≠r2。

步骤2将编码中位置r1处的元素插入到位置r2之前。

(2)邻域变异操作，具体步骤如下：

步骤1在父代P中选择3个不同的元素，生成它的全排列；

步骤2在生成的全排列中任意选择一种代替原来的的排列，其他位置的元素保持不变。

邻域变异操作如图9所示。针对基于机器的编码，本文采用多点变异操作。具体步骤如下：

步骤1产生随机0,1数组R，数组长度为编码长度l。

步骤2找出数组R中元素为1的位置，在对应工序的可加工机器集中随机选取机器替换之前的机器，其他位置的机器保持不变。

多点变异示意图如图10所示。

3.5 单分子反应

在化学反应算法中，分子与墙壁的碰撞反应和分解反应是针对单个个体进行的操作。

3.5.1 分子与墙壁的轻微碰撞反应

设分子ω为父代解，通过与墙壁碰撞产生子代分子ω′，本文采用上述变异操作来实现。分子与墙壁碰撞反应步骤如下：

步骤1对ω进行变异操作，得到子代ω′。

步骤2判断是否满足与墙壁碰撞条件，若满足，则更新个体和能量中心的能量。

3.5.2 分解反应

步骤1判断是否满足分解反应的条件。本文设定的分解条件为在规定的迭代次数α内分子的目标函数值没有更新。

步骤2产生一随机整数r，r∈[-l,l]，其中l为编码长度。

步骤3若r<0,则将编码中前r个元素移动到编码末端，若r>0，则将编码中最后r个元素移动到编码前端。设r分别为-1和2，其针对工序编码的移动过程如图11所示。

3.6 多分子反应

在CRO算法中，分子间的轻微碰撞反应和合成反应是两种单分子反应。

3.6.1 分子间的轻微碰撞反应

分子间的碰撞反应是两个父代分子通过碰撞产生两个子代分子的过程。本文采用上述交叉操作来实现。其步骤如下：

步骤2判断是否满足分子间碰撞反应的条件，若满足则修改分解后的两个子代分子的动能并将这两子代分子加入到新种群。

3.6.2 合成反应

合成反应是两个分子合成一个分子的过程。本文采用距离保持交叉操作[17]来实现合成反应。基于工序编码的距离保持交叉操作如图12所示，其步骤如下:

步骤1判断是否满足合成反应的条件。本文设定的合成条件为两父代分子ω1和ω2的动能同时小于设定的动能下限β，即KEω1≤β且KEω2≤β。

步骤2随机选取两父代分子，比对两父代分子中各位置的元素。

步骤3将父代分子中同一位置相同的元素复制到子代，随机产生其他位置的元素，得到子代分子。

3.7 精英保留机制

在CRO算法中，原本适应度很高的分子也有可能生成新的适应度很低的分子。这一现象保证了分子群的多样性，使算法有很好的全局搜索能力，但是这也降低了算法的收敛速度。因此，在改进的算法中引入了精英保留策略。该策略将父代中最优的Popsize×Pr(Popsize为种群规模，Pr为保留概率)个个体直接保留到子代而不进行分子操作。

另外本文还采用二元锦标赛策略来完成选择操作，即每次从分子群中随机选择两个不同个体，同时产生1个0～1之间的随机值，若随机值小于选择概率Ps(一般Ps=0.8),则选择两者中适应度较高的个体，否则选择另一个适应度低的个体。被选中的个体放到下一代种群中参与分子操作。

3.8 禁忌搜索算法

TS算法是一种高效的局部搜索算法，目前，禁忌搜索在组合优化、生产调度、机器学习、电路设计和神经网络等领域得到广泛应用。TS算法由邻域结构、移动评价策略、禁忌表、禁忌长度、特赦准则和终止准则等组成。

本文采用Gambardella等提出的邻域结构[14]和移动评价策略。本文中禁忌对象用(v,k)来表示，其中v表示被移动的工序，k表示移动之前工序v的加工机器。禁忌长度设置为TM(v,k)=P+Mv，其中P为当前关键路径长度，Mv为工序v的可加工机器数。当得到的新解优于当前最优解时即满足特赦准则。当禁忌搜索算法迭代次数达到最大迭代次数，或者在给定的迭代次数下仍没有得到更优解，则算法终止。

禁忌搜索算法程序框图如图13所示。

3.9 CROTS混合算法框架

本文针对FJSP的特点提出融合CRO算法与TS算法的混合算法(CROTS)。CRO算法负责全局搜索而TS算法进行局部搜索，一方面让问题能求得更优的解，另一方面加快了种群收敛速度，使混合算法在分散搜索和集中搜索之间达到更合理的平衡。本文提出的CROTS混合算法步骤如下：

步骤1设置试验参数，初始化分子种群，并根据式(7)和式(8)分别计算每个分子的势能和动能。

步骤2根据精英保留机制，从分子群中选择分子进行优化操作。

步骤3生成0～1的随机数r，若r>MoleColl则进行单分子反应，即执行步骤5；否则进行多分子反应，即执行步骤6。

步骤4判断是否满足分解条件，若满足则进行分解反应，否则进行分子与墙壁的轻微碰撞反应。

步骤5判断合成条件是否满足，若满足则进行合成反应，否则进行分子间的轻微碰撞反应。

步骤6对分子进行禁忌搜索。

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

步骤8算法结束。

3.10 算法复杂度分析

CROTS算法的关键步骤的时间复杂度如表2所示，其中n为编码长度，m为进行操作的个数。混合算法中，每一代有Popsize个分子，其中有Popsize×MoleColl个分子进行多分子反应，Popsize×(1-MoleColl)个分子进行单分子反应，所有的分子都进行禁忌搜索。因此混合化学反应算法的总时间复杂度为：

Popsize×(1+MoleColl)×O(n2)+Popsize×(1-MoleColl)×O(n)。

可以看出，本文提出的混合化学反应算法的时间复杂度与问题规模、种群大小、多分子反应概率等参数有关。

表2 算法中关键步骤的时间复杂度

4 实验结果与分析

上述CROTS算法用C++语言编程，程序运行环境为P4 CPU，主频3.0 G，内存为2 G。为了验证和比较本文所提算法的性能，选取Brandimarte[22]和Fattahi[23]基准问题对算法进行测试。算法的参数设置和实验结果分析如下。

4.1 参数设置

混合化学反应算法的参数主要有：分子种群规模Popsize，多分子碰撞反应概率MoleColl，初始动能InitialKE，保留概率Pr，选择概率Ps，动能下限β，分子没有改善的最大迭代次数α，最大迭代次数maxGen，其中对算法性能影响较大的有MoleColl，InitialKE，β和α。为了考察参数对算法性能的影响，本文采用试验设计方法进行分析，选取Brandimarte[22]基准问题中的MK06案例进行测试，每个试验参数设置4个不同的水平值，如表3所示。选用正交表L16(45)的前四列，算法中分子种群规模Popsize=500，保留概率Pr=0.01，选择概率Ps=0.8，最大迭代次数maxGen=200。为了评估算法的性能，计算算法所得最优解ALG与参考最优值OPT的相对百分偏差RPD=(ALG-OPT)/OPT×100%。算法在每组参数下独立运行20次，以20次实验所得RPD的均值作为影响值(RV)，统计结果如表4所示，使用Minitab 16进行数据分析得到参数的水平值趋势如图14所示。

表3 CROTS算法中参数的水平值

表4 正交表和RV值

续表4

由图14可见，设置不同的参数，CROTS算法的性能也有所不同，根据CROTS算法的设计思想可知，MoleColl，InitialKE，β和α这4个参数决定了算法的收敛速度和跳出局部最优的能力。在算法参数选择时应根据问题对这4个参数取合理的值，根据试验结果本文设定这4个参数的值,如表5所示。

表5 CROTS算法的参数的参考值

4.2 基准实例测试

表6和表7分别显示了Brandimarte和Fattahi实例测试结果的比较，其中对于Brandimarte实例分别给出了与Mastrolilli等[14]采用的禁忌搜索算法;Gao等[16]的遗传算法与变邻域搜索结合的混合算法;Li等[24]的遗传算法与禁忌搜索结合的混合算法HA,以及Wang等[11]的人工蜂群算法的比较。对于Fattahi实例，表8中的AIA算法和HHS算法来自于Yuan[25]等，M2算法来自于Demir和Isleyen[26]。此外，图15～图18分别显示了实例Mk06、MK10、SFJS10和MFJS10的甘特图，表8和表9显示了不同算法对不同案例相对于最优解的RPD值。

表6 对于Brandimarte基准实例测试结果的比较

表7 对于Fattahi基准实例测试结果的比较

表8 对于Brandimarte基准实例不同算法的RPD

表9 对于Fattahi基准实例不同算法的RPD

续表9

表6和表7中带*的数字表示该问题的当前最优解。由表6和表8可知，对于Brandimarte基准实例，本文所求的结果优于M.G.的TS算法、Gao的VNS算法以及Wang的ABC算法，与Li的HA算法所求结果一致，本文效率稍高。这是因为M.G.的TS算法没有进行全局搜索，而本文采用了化学反应算法进行全局搜索，提高了算法性能。本文设计的化学反应算法既有两个分子的相互作用又有分子与墙壁的作用，在传统的交叉、变异操作上又加入了一个分子分解为两个分子的操作，增强了算法的多样性，提高了算法性能。本文算法中采用的交叉、变异等算子都是根据柔性作业车间调度问题的特征设计，本文采用邻域结构也被证明是当前效果较好的，故本文能求出案例的当前最优解。由表7和表9可知，对于Fattahi基准实例，本文所提算法也能求得当前最优解，证明了本文算法的有效性。

5 结束语

本文通过研究FJSP提出了一种结合CRO和TS两种算法的混合算法，给出了算法的基本框架，并且扩展了CRO的4种基本反应，设计了有效的交叉、变异算子，采用了高效的邻域函数。混合算法融合了CRO的全局搜索能力和TS的局部搜索能力。在解码过程中采用了插入式贪婪解码算法以生成主动调度解，使算法更加高效。此外，在算法中引入了精英保留机制以保留较优个体，从而提高算法的收敛速度。最后，使用基准实例测试了该算法，并与其他算法的结果进行了比较，证实了所提算法的有效性。

未来将进一步提高CRO算法的性能，深入分析FJSP问题的结构特征提出更好的邻域结构函数以便减少算法的无效搜索，提高算法的收敛速度。另外，还考虑将所提出的混合算法用于求解多目标的FJSP，以期获得更好的结果。