混合交叉熵算法求解复杂零等待流水线调度问题

张梓琪 钱 斌 胡 蓉

(1.昆明理工大学机电工程学院，云南昆明 650500;2.昆明理工大学信息工程与自动化学院，云南昆明 650500)

1 引言

生产调度在现代化生产制造系统中起着十分重要的作用，尤其是在德国实施工业4.0战略、美国加速工业互联网建设和中国推进智能制造2025的国际化大背景下，生产调度在智能制造的关键环节上显得至关重要.流水线调度问题(flow shop scheduling problem，FSSP)是制造系统中广泛存在的一类典型生产调度问题.数学上，FSSP属于NP-hard问题.近年来，FSSP在计算机科学和运筹优化领域被广泛关注与研究[1-3].零等待流水线调度问题(no-wait flow-shop scheduling problem，NFSSP)是FSSP的特例.NFSSP作为一类具有实际工程应用背景的重要生产调度问题[4-5]，要求加工过程必须连续进行，不能出现中断或等待.NFSSP通常假设工件加工时间包含设置时间且工件无释放时间，但实际生产制造中难以避免会频繁出现切换生产工况、更换工具和调整设备参数等情况，这将导致存在除加工时间以外的设置时间[6-7].一般而言，设置时间可以划分为两种类型[6]，即序无关的设置时间(sequence independent setup times，SISTs)和序相关的设置时间(sequence dependent setup times，SDSTs).绝大多数情况下，设置时间不仅由当前机器上待加工的工件决定，还依赖于该机器上前一个加工的工件[8].因此，考虑带有序相关的设置时间的情形更加具有现实意义.此外，由于实际生产过程中待加工的工件在进入生产车间之前存在运输过程，即工件的开始加工时间不仅取决于当前可供加工使用的机器是否空闲，还必须等到待加工工件的释放时间满足条件后才能开始加工，故机器的设置时间和工件的释放时间是实际生产中不可回避的两个重要约束.因此，研究带有序相关的设置时间和释放时间的零等待流水线调度问题(no-wait flow-shop scheduling problem with sequence dependent setup times and release times，NFSSP-SDSTs-RTs)具有重要的理论价值和实际意义.为了满足实际生产过程中准时制(just in time，JIT)的生产要求[9]，本文综合考虑生产过程中的提前与延迟，以最小化提前和延迟总代价(total earliness and tardiness，TET)为优化目标.根据调度问题的三元组描述法[10]，优化目标为TET的NFSSP-SDSTs-RTs可归类为Fm/no-wait，STsd，在计算复杂度上，优化目标为TET的一类单机调度问题Tj)已经被证明为NP-hard问题[11]，显然，这一类单机调度问题又可以约归为Fm/no-wait，故本文研究的问题属于NPhard.同时，NFSSP-SDSTs-RTs也是NFSSP中较为复杂的一类问题[4-5]，许多NFSSP及其扩展问题都可以归约为NFSSP-SDSTs-RTs.因此，研究优化目标为TET的NFSSP-SDSTs-RTs问题具有重要的学术意义和工程应用价值.

近年来，针对复杂NFSSP得到了较多的关注与研究[12-23].在问题求解层面，NFSSP-SDSTs-RTs具有非线性和强耦合等诸多特性.在问题的建模方面，一般可以建立两类模型，即数学规划模型和排序模型.Bianco等[12]首次建立了NFSSP-SDSTs-RTs的数学规划模型，并从数学上推导得出该类问题可以等效为带准备时间的非对称旅行商问题(asymmetric travelling salesman problem with ready time，ATSP-RT).针对数学规划模型的求解目前主要有分支定界(branch and bound，B&B)、列生成(column generation)和Benders分解(benders decomposition)等方法.传统运筹学方法可通过商用求解器如GUROBI或CPLEX等实现问题解空间的全部遍历或部分遍历，一般能在较短时间内给出小规模问题的最优解，但对于大规模复杂优化问题，这类方法在可接受时间内获取满意解的难度较大.对于排序模型，通常是由一系列用于确定工件在机器上开始加工时间的计算式组成，问题的约束条件隐式包含在排序编码中，通过解码来实现对约束的处理.排序模型一般采用近似算法、启发式算法和智能优化算法进行求解.近似算法如动态规划算法，需要构造Behrman状态递推方程，并严格证明问题最优解包含在每次迭代的状态空间中方可使用.动态规划算法虽然理论上能保证搜索到问题的最优解，但随着问题规模的递增，其状态空间往往呈指数级增长，导致动态规划算法难以实际应用.启发式算法通常基于特定的调度规则或者问题的特性构造可行解，譬如，Johnson算法、Palmer算法、Nawaz-Enscore-Ham算法、Gupta算法和Campbell-Dudek-Smith算法等[10].这类启发式算法可以在较短时间内获得可行且质量相对较高的问题解，但并不能保证全局意义上的最优性.智能优化算法基于计算智能机制，结合问题的特性并通过模拟各类自然机理来实现对问题解空间的有效搜索，往往在较短时间内就能获得复杂生产调度问题的最优解或满意解.He等[13]指出，智能优化算法的有效性是由调度问题排序模型的解空间特性和解的特点，以及智能优化算法自身机制共同决定.由于离散优化问题排序模型的解空间具有“极为扁平”性，解与解之间具有“无梯度”性，使得智能优化算法在较短时间内仅需搜索排序模型解空间中很小的区域，就能够在目标函数上获得明显多于传统运筹学方法的较优解，从而能够有效地求解各类复杂调度问题.

目前，应用智能优化算法求解复杂NFSSP已有大量的研究成果[4].根据所研究的问题可以划分为3种类型，分别是带SISTs的NFSSP(NFSSP-SISTs)[14]、带SDSTs 的NFSSP(NFSSP-SDSTs)[15，22]和NFSSPSDSTs-RTs[12，23].在NFSSP-SISTs的研究方面，Allahverdi等[6-7]针对NFSSP-SISTs进行了全面且系统的综述.Ding等[14]针对优化目标为最小化最大完工时间的NFSSP-SISTs，提出一种带禁忌机制的改进迭代贪婪(tabu mechanism improved iterated greedy，TMIIG)算法进行求解.TMIIG算法在IG算法的基础上加入禁忌重构机制，利用启发式算法生成高质量的初始解，同时设计基于多种邻域结构的局部搜索，进一步提高算法的性能.实验结果表明TMIIG 是求解NFSSPSISTs的有效算法.在NFSSP-SDSTs的研究方面，优化指标通常为最小化最大完工时间.Ruiz和Allahverdi[15]设计了几种启发式算法并提出一种基于Insert邻域结构的迭代局部搜索算法进行求解.Araujoa和Nagano[16]分析问题结构性质并提出一种有效的启发式算法进行求解.Jolai等[17]分别设计了基于种群的模拟退火算法、适应性帝国主义竞争算法以及二者混合的优化算法对柔性NFSSP 进行求解.Allahverdi和Aydilek[18]提出了遗传算法、差分进化算法和模拟退火算法等一系列有效算法进行求解，计算结果表明在相同的计算时间内所提改进模拟退火算法(ISA-2)的性能显著优于其他算法.Nagano等[19]提出一种基于进化聚类机制的混合优化算法进行求解.Samarghandi和ElMekkawy[20]设计一种矩阵编码方式对解序列进行编码，同时提出一种基于该编码方式的的粒子群算法进行求解.Samarghandi[21]提出一种基于遗传算法(GA)的优化算法进行求解，有效地解决了服务器端约束对最大完工时间的影响，该方法改进了标准测试集中部分已知的最优解.此外，除了考虑最小化最大完工时间指标外，Arabameri和Salmasi[22]还进一步针对优化目标为最小化加权的提前和延迟时间的NFSSPSDSTs，设计一种基于禁忌搜索和粒子群优化的混合算法进行求解.在NFSSP-SDSTs-RTs的研究方面，优化指标主要考虑最小化最大完工时间.Bianco等[12]分别提出一种最优添加启发式(best adding heuristic，BAH)算法和一种最优插入启发式(best insertion heuristic，BIH)算法进行求解，计算结果表明BIH算法的性能优于BAH算法.Franca等[23]设计一种混合遗传算法(HGA)进行求解，并在HGA中嵌入递归局部搜索方法进一步增强算法的性能.综上所述，目前针对NFSSP-SDSTs-RTs的研究仍十分有限，故对其开展建模与算法研究具有较大的理论和实际意义.

交叉熵算法(cross entropy algorithm，CEA)是一种估计复杂网络中稀有事件的有效方法[24].CEA通过将优化问题转化为该问题对应解空间上的概率分布估计问题，利用小概率事件信息更新概率模型参数，采样概率模型获得新解，实现对问题解空间的有效搜索，引导搜索方向逐渐逼近全局最优.近年来，CEA及其改进算法已被广泛应用于求解各类组合优化问题.Caserta等[25]针对带设置成本约束的一类背包问题，提出了一种混合CEA进行求解，仿真实验验证了所提算法的性能优于分支定界法.Santosa等[26]针对优化目标为最小化最大完工时间的零等待作业车间调度问题，设计了一种混合交叉熵遗传算法(HCEGA)进行求解.HCEGA利用交叉熵算法引导全局搜索方向，同时采用遗传算法的交叉和变异操作来丰富算法的搜索行为，以实现对解空间不同区域的有效搜索.Wang等[27]针对加工时间不确定的炼钢连铸生产调度问题，提出了一种串级CEA进行求解.Wang等[28]针对可再生能源发电不确定性的电力系统优化调度问题，提出了一种多目标CEA进行求解.因此，本文选用CEA为基本框架，设计求解NFSSP-SDSTs-RTs的有效算法.本文研究的NF-SSP-SDSTs-RTs具有诸多复杂性，不仅要考虑零等待的工艺约束和序相关设置时间与释放时间的加工约束，同时还需要满足准时制的生产要求.问题求解的关键难点在于可行解空间十分庞大且复杂，在这样庞杂的解空间中发现优质解是一种极小概率的事件.因此，非常适合采用CEA进行求解.目前，尚无应用CEA求解NFSSP-SDSTs-RTs的相关研究.

本文研究NFSSP-SDSTs-RTs的建模与求解.首先，在问题建模方面，建立以最小化总提前和延迟为优化目标的NFSSP-SDSTs-RTs的排序模型;其次，在问题求解方面，分析问题零等待性质，设计一种快速评价方法，有效降低评价解的计算复杂度;然后，在算法设计方面，采用交叉熵算法学习并积累优质解的结构特征，建立概率模型对优质解中工件块分布进行合理且有效估计，通过采样概率模型生成新种群，实现对解空间中优质区域的全局搜索.同时，利用带两种搜索策略的快速局部搜索对全局发现的优质区域进行更为细致且深入的搜索.最后，通过大量仿真实验和算法对比分析，验证了所提算法的有效性和鲁棒性.

2 问题描述

2.1 符号定义

本文所涉及的数学符号及定义如表1所示.

表1 符号表Table1 Notations

2.2 问题模型

NFSSP-SDSTs-RTs描述为:n个工件在m台机器上进行加工，每个工件的加工时间、序相关设置时间和释放时间是确定的.每台机器在同一时刻只能加工一个工件，并且每个工件在同一时刻只能由一台机器进行加工.为了满足零等待工艺约束，各工件在机器上必须连续不断地进行加工，任意时刻都不允许中断.所有工件按照相同工序在每台机器上进行加工且无任何等待.每台机器上当前加工工件的完工时间到下一个工件的开始加工时间，这期间内存在依赖于工件加工顺序的序相关设置时间.此外，如果生产线上第一台机器已经准备就绪但待加工工件还未到达，此时该机器只能处于空闲等待状态直到工件的释放时间满足要求为止.NFSSP-SDSTs-RTs的优化目标为最小化总提前和延迟时间.基于以上描述，建立如下模型:

调度目标为在所有可行调度方案集合Π中找到总提前和延迟时间最小的最优调度方案π*，即

其中:式(1)为机器l上工件ji-1与工件ji完工时间之间最小延迟的计算公式，式(2)为第1台机器上工件ji-1与工件ji开工时间之间最小延迟的计算公式，式(3)为工件ji在第1台机器上开工时间的计算公式，式(4)至式(6)分别为工件ji在机器m上的完工时间、提前时间和延迟时间的计算公式，式(7)为优化目标，式(8)表示在所有方案集合Π中找到使得TET最小的最优排序π*.图1为n=3和m=3的示例甘特图.

图1 示例的甘特图Fig.1 The Gantt chart for the example

2.3 快速评价方法

根据NFSSP-SDSTs-RTs排序模型的零等待特征可知，式(2)计算ji-1与工件ji在第1台机器上开工时间之间的最小延迟仅由工件ji-1与工件ji的加工顺序所决定，与工件ji-1与工件ji的释放时间无关[12].本节基于该性质提出计算TET(π)的快速评价方法，以有效减少计算式(7)的复杂度.具体描述如下:

步骤1利用工件jk(k=1，2，···，n)的加工时间计算并保存对应的spjk，其中spjk=

步骤2利用工件ju与工件jv在各机器上的加工时间和对应的序相关设置时间，由式(1)和式(2)计算并保存工件ju与工件jv在第1台机器上开工时间之间的最小延迟，其中ju，jv ∈{1，2，···，n}，jujv;

步骤3对于给定的解序列π=[j1，j2，···，jn]，计算TET(π)时，式(3)和式(4)中涉及到的直接调用步骤1和步骤2中提前计算并保存的数据.

根据上述步骤可知，采用快速评价方法能够较大程度上降低评价解的计算复杂度，使第2.2节式(7)中计算TET(π)的复杂度从O(nm)下降到O(n).基于问题特性的快速评价方法有利于提高算法的搜索效率，使得算法在相同时间内评价解的次数增加，从而能较快探索到更多优质解所在区域.

3 混合交叉熵算法

3.1 解的编码、工件块与相似块

智能优化算法求解基于排序模型的零等待流水线调度问题，通常采用基于工件加工顺序的编码方式来表征问题的解[14，23].种群中的每个个体都对应问题解空间中的一个解，即一个长度为n的工件序列π.例如，π=[2，6，4，8，3，9，1，5，7]表示优先加工工件序列π首位置上的工件2，其次加工位置2上的工件6，最后再加工位置8上的工件7.

工件块定义为解序列中相邻的两个工件.譬如，对于序列π=[2，6，4，8，3，9，1，5，7]，位置1到位置8上的工件块分别为[2，6]，[6，4]，[4，8]，[8，3]，[3，9]，[9，1]，[1，5]和[5，7].在此基础上，相似块(similar block)定义为不同解序列中共同存在的工件块.为便于理解，随机生成一个n=9的NFSSP-SDSTs-RTs实例，其解空间规模为9!=362880.通过穷举所有解，可得最优解序列的TET目标值为1520.图2给出部分解序列和对应的TET目标值，其中π1为该实例的最优解，相应的工件块[2，6]，[8，3]，[1，5]和[9，7]为相似块.

图2 工件块示意图Fig.2 The job blocks for the example

3.2 交叉熵算法

交叉熵算法通过建立概率模型对解空间中优质解的结构特征进行合理估计[24].对于本文研究的问题，应用交叉熵算法的关键在于概率模型的确定.根据图2可以看出，不少相似块在各解序列中分布的位置不同.以相似块[8，3]为例，其在π1～π3和π85～π87中出现的位置不同.如果采用已有文献中常用的n×n维概率模型，只能积累优质解中相似块的序信息(即相邻工件前后顺序信息)，但无法准确学习并积累相似块的位置信息.例如，π1～π3和π85～π87中相似块[8，3]的“出现频率”或“出现概率”只能存贮在n×n维概率模型下标为(8，3)的元素中，其中下标本身使得8前3后的序信息可以保留，但无法准确区分相似块[8，3]在各解序列中的具体位置信息，不能够合理地学习优质解的结构特征[29].对于n×n维的概率模型采样生成新解或个体，通常难以引导算法在优质解附近区域执行有效搜索.因此，本文设计了一种n×(n×n)维概率模型P，用于学习并积累优质解序列中的优良模式信息.令Pi=(pi，1，pi，2，···，pi，n×n)为概率模型P中的第i行，且pi，n(u-1)+u=0(u=1，···，n)，同时Pi满足归一化条件，即=1.此时对于任意个体π第i个位置上工件块[ji=u，ji+1=v]的“出现概率”可以存储于元素pi，n(u-1)+v中.例如，对于图2中优质个体π1和π2中的相似块[8，3]，则可以分别存贮在P4的p4，66元素中和P5的p5，66元素中，从而能够同时保留工件块或相似块的序信息和位置信息.当概率模型P的结构确定以后，可以利用如下的交叉熵方法以及选取的优质个体来确定算法每一代概率模型P的取值.

基于交叉熵方法，第2.2节式(7)中的优化目标可以等价表示为

其中:X*对应最优个体π*，γ*是X*在目标约束S下的最优值.假设X服从参数为μ的概率分布f(X;μ)，式(9)中对于X*的求解则可以转化为确定最优分布f(X;μ)的一个概率估计问题

其中:Pμ和Eμ为概率度量和期望，EμI{S(X)≤γ}表示样本中目标优于γ的期望.式(10)中当γ越接近最优值，则概率ξ(γ)越小，此时S(X)≤γ为小概率事件.如果要估计式(10)的期望，必须进行大量的抽样.为了减少抽样，利用重要性抽样方法对Ω上的分布g采样得到N个样本{X1，···，XN}，对ξ进行估计

为了用式(11)中的ˆξ估计ξ，要选择合适概率分布参数μ使式(12)中g的最优分布gopt与分布f的差距最小.通常采用Kullback-Leibler(K-L)距离来衡量两个概率分布之间的相似度，最小化K-L距离等价于最小化交叉熵[22].令概率分布Fp(X)对Fq(X)分别是对应随机变量X的两个分布，则分布Fp对分布Fq的交叉熵定义为Fq分布的自信息对Fp分布的期望

步骤1选择一组随机样本{X1，···，XN}，即初始化种群，并通过式(18)对概率分布P(0)初始化;

步骤2通过轮盘赌采样方式，对第gen次迭代的P(gen-1)进行采样，产生一定数量的可行解，构成本次迭代的一组随机样本{X1，···，XN}.式(11)中样本数N即对应问题解空间中的种群规模popsize;

步骤3针对步骤2中所产生的N个随机样本{X1，···，XN}的{S(X1)，···，S(XN)}进行排序，计算样本的φ分位数，使得P(S(X)≤γgen)=φ.再按式(18)确定本次迭代的概率分布P(gen);

步骤4为了逐代学习并积累优质解的工件块信息，对于第gen+1次迭代的概率分布P(gen)，通过学习速率r调控更新过程如下:

3.3 局部搜索

为了增强算法的搜索能力，需要对CEA全局搜索发现的优质解的邻近区域再进行精细的局部搜索.通常，在常用邻域搜索操作中，交换(Interchange)和插入(Insert)操作产生的新解能最大程度地保持原解的工件块特征，有助于进一步细致的局部搜索.因此，本文对算法全局获取的优质解进行基于这两种操作的局部搜索.具体来说，先对全局搜索得到的优质解进行Interchange操作，扰动产生一个中间解序列;随后再对扰动后的中间解序列进行基于Insert邻域结构的局部搜索.

1.犯罪客体是复杂客体，既侵犯了国家人体器官管理制度，也侵犯了被组织出卖器官者的生命、健康等权益，前者是主要客体，后者是次要客体。依照人体器官移植条例实施的移植不是犯罪，是阻却违法性的医疗行为，指医务工作者出于正当目的，经过当事人本人或其监护人、保护人同意，采取适当方法割伤人的身体的行为。［3］条例强调任何组织或个人不得以任何形式买卖人体器官，不得从事与此有关的活动，器官捐献应当遵循自愿、无偿原则。人体器官不是商品，不能自由买卖，否则，人们的价值观、伦理观将受到严重冲击，社会管理秩序将陷于混乱，衍生杀人、伤害等犯罪，具有极大社会危害性。

3.3.1 邻域结构

3.3.2 快速邻域搜索方法

图3 πn，i，l插入操作示意图Fig.3 Insert operation for πn，i，l

3.3.3 邻域搜索策略

3.3.4 局部搜索流程

局部搜索流程如图4所示.

具体步骤如下:

步骤1令L=0，设置参数KM=5，KK=1.πi0为全局搜索得到的最优解，令πi-t=πi0.随机生成u和v且满足|u-v|＞n/3，对πit执行KM次基于Interchange(πit，u，v)的扰动操作;

步骤2若L/=1，则对解πi-t执行邻域搜索操作得到解πi-1，即πi-1=TS-SP-FindBestFisrtSkipNInsert(πit，KK);若L=1 则跳转至步骤4;

步骤4若TET(πi-t)≤TET(πi0)，则πi0=πi-t;

步骤5输出更新后的解序列πi0.

3.4 算法流程

根据上述描述，HCEA的具体步骤如下:

步骤1由式(1)和式(2)计算并保存和

步骤2初始化种群并初始化工件块概率分布P(0)，同时设置算法的关键参数popsize，φ，r;

步骤3对种群进行评价并得到πbest;

步骤4通过交叉熵算法对解空间执行全局搜索，采样概率分布产生popsize个个体，更新πbest;

步骤5选择popsize×φ个优质个体组成优质子种群，更新概率分布.对优质子种群执行第3.3.4节中的局部搜索并更新πbest;

步骤6判断是否满足终止条件，若不满足则转到步骤4，否则终止循环.

由以上算法步骤可知，步骤4建立概率分布学习并积累优质解的结构特征，用于在解空间中进行全局搜索，发现优质解所在区域.步骤5针对全局搜索找到的优质区域进行较为细致的局部搜索，进一步提升算法的性能，同时使算法在全局和局部搜索上达到合理平衡.HCEA流程如图5所示.

3.5 算法复杂度分析

HCEA关键操作的计算复杂度分析如下:

1) 快速评价方法(见第2.3节)计算Lji-1，ji和spji的复杂度分别为O[n2×m]和O[n×m].基于快速评价方法(见第2.3节)计算TET(π)的复杂度为O[n]，相应地对种群进行快速评价的复杂度为O[popsize×n].

2) HCEA全局搜索阶段，对种群进行快速排序的复杂度为O[popsize×log popsize]，初始化工件块概率分布的复杂度为O[n×n(n-1)]，采样概率分布产生新种群的复杂度为O[popsize×n(n-1)]，更新概率分布的复杂度为O[n×n(n-1)×popsize×φ].

3) HCEA局部搜索阶段，每执行一次快速邻域搜索TS-SP-FindBestFirstSkipNInsert(π，i)(见第3.3.3节中步骤3)的计算复杂度为O[n×logn].针对邻域NInsert(π)进行快速邻域搜索(见第3.3.3节中步骤2到步骤5)的复杂度为O[(n-KK+1)×n×logn].

综上，并考虑通常n ＞m，第3.4节中HCEA的整体计算复杂度Tcc为

根据以上分析可知，本文第2.3节提出的快速评价方法和第3.3.2节提出的快速邻域搜索方法可以有效降低计算复杂度.同时，由式(21)可看出HCEA的整体计算复杂度取决于工件数n，且最高次仅为3.除了问题规模以外，影响算法计算效率的因素主要还有种群规模popsize和优势个体占比φ，实验部分将对这两个因素进行详细探讨.

4 实验设计与分析

为了验证HCEA求解NFSSP-SDSTs-RTs的有效性，进行大量数值实验与算法对比分析.所有测试实验均通过Embarcadero RAD Studio 10.4编程实现，实验平台Windows 10，Inter(R)Core(TM)i7-8700 CPU 3.20 GHz处理器16 GB内存的PC机.

4.1 实验设置

由于目前还没有适合NFSSP-SDSTs-RTs的标准测试集，根据文献[18]中解决流水线调度问题的方法，针对工件n与机器m，随机生成不同规模的测试问题，n×m范围为{20，30，50，70，100}×{5，10，20}.工件加工时间和序相关设置时间分别服从区间[1，100]和[0，100]的离散均匀分布，工件释放时间为区间[0，150nα]内随机生成的一个整数，其中参数α控制工件的到达速率.此外，参数α的选择直接影响问题解空间的特性与问题求解难度.为了能够充分地测试算法的性能，分别选取7种不同的α值:0，0.2，0.4，0.6，0.8，1.0和1.5，一共生成105个算例进行测试.工件的交货时间由以下步骤产生:

步骤1对测试问题p随机产生一个工件排序;

步骤2计算步骤1产生的排序中各工件ji的完工时间Cp，ji，i=1，2，···，n;

步骤3工件ji的交货时间按式(22)产生:

4.2 性能指标

定义π(α)为在参数α水平下的一个工件排序，π′(α)为各工件按照其释放时间的升序得到的排序，即满足先到先加工原则.TET(π(α))为π(α)的总提前与延迟时间.TETavg(π(α))为算法运行30次对应TET(π(α))的平均值.为了评价各算法对所有测试问题的求解能力，定义两个统计指标，分别为

其中:式(23)中的ARI(α)为TET(π′(α))的平均百分提高率，式(24)中的SD(α)为TET(π(α))的标准偏差.令Sα为参数α取值集合，|Sα|为Sα中元素个数.为评价算法的整体性能，定义两个性能指标，分别为:

显然，性能指标ARI越大越好，SD越小越好.

4.3 参数设置

HCEA的3个关键参数分别为种群规模popsize，优势个体占比φ和学习速率r.本节选取中等规模问题n×m=50×10，其中α取第3.1节中7种不同值，采用实验设计方法(design of experiment，DOE)[31]探讨各参数对算法性能的影响.算法参数选取5个水平值，如表2所示.选择规模为L25(53)的正交实验进行参数整定.算法在每种参数组合下独立重复运行30次，终止条件为4×50×10毫秒.算法在不同参数组合下各算例得到的ARI的均值作为平均响应值(average response value，ARV).显然，ARV的值越大，其所对应参数组合下算法的求解性能就越好.

表2 主要参数与水平表Table 2 Main parameters and level table

正交表和参数响应值如表3所示，根据正交表统计得到各参数的平均响应值如表4所示，各参数对算法性能影响的响应趋势如图6所示.由表4及图6可以发现，不同的参数组合会对算法性能产生较大的影响.根据以上实验结果可知，HCEA的最佳参数设置为popsize=120，φ=0.2，r=0.2，该参数也将应用于后续的算法的性能测试与比较中.

图6 参数响应图Fig.6 Response diagram of parameters

表3 参数设置的正交表Table 3 Orthogonal table of parameter settings

表4 各参数不同水平下的平均响应值和等级表Table 4 Average response value and level table at different levels of each parameter

4.4 仿真结果与比较

为了验证第3.3节所设计的局部搜索方法的有效性，本节首先针对第3.3.2节提出的快速邻域搜索方法和第3.3.3节提出的邻域搜索策略进行探究.对HCEA及其6种不同的局部搜索变体进行比较，具体如下:

1) HCEAV1:将HCEA的带两种策略的快速搜索方 法TS-SP-FindBestFirstSkipNInsert(π，KK)(KK=1)中的快速邻域搜索方法(见第3.3.2节)、首次改进跳出策略和改进后变邻域搜索策略(见第3.3.3节)移除，局部搜索中仅保留快速评价方法(见第2.3节).

2) HCEAV2:将HCEA的带两种策略的快速搜索方 法TS-SP-FindBestFirstSkipNInsert(π，KK)(KK=1)中的快速评价方法(见第2.3节)、首次改进跳出策略和改进后变邻域搜索策略(见第3.3.3节)移除，局部搜索中仅保留快速邻域搜索方法(见第3.3.2节).

3) HCEAV3:将HCEA的带两种策略的快速搜索方 法TS-SP-FindBestFirstSkipNInsert(π，KK)(KK=1)中的首次改进跳出策略和改进后变邻域搜索策略(见第3.3.3节)移除，局部搜索中保留快速评价方法(见第2.3节)和快速邻域搜索方法(见第3.3.2节).

4) HCEAV4:将HCEA的带两种策略的快速搜索方 法TS-SP-FindBestFirstSkipNInsert(π，KK)(KK=1)中的改进后变邻域搜索策略(见第3.3.3节)移除，局部搜索中保留快速评价方法(见第2.3节)，首次改进跳出策略(见第3.3.3节)和快速邻域搜索方法(见第3.3.2节).

5) HCEAV5:将HCEA的带两种策略的快速搜索方法TS-SP-FindBestFirstSkipNInsert(π，KK)中KK设置为0.6，其中表示对δ向上取整.

6) HCEAV6:将HCEA的带两种策略的快速搜索方法TS-SP-FindBestFirstSkipNInsert(π，KK)中KK设置为「0.3，其中表示对δ向上取整.

对于不同规模问题，HCEA与6种变体算法在运行时间为4×n×m毫秒下的统计结果如表5所示.图7为HCEA与6种变体分别在不同工件数n下的柱状统计图.根据第4.2节可知，ARI指标的值越大越好，故首先从图7中可明显看出HCEAV3得到的ARI值均比HCEAV1和HCEAV2的好，验证了本文所设计的快速评价方法与快速邻域搜索方法能够有效降低计算复杂度，使得算法在相同时间内能够评价更多的解.其次，HCEAV4得到的ARI值比HCEAV3的好，说明首次改进跳出策略能够对子邻域NInsert(π，i)进行有效地搜索，在相对有限的时间内能搜索到更多区域.此外，HCEAV6得到的ARI值比HCEAV5的好，表明参数KK直接影响局部搜索的性能，若其值设置过大，则邻域NInsert(π)的实际搜索范围将缩小，进而无法对子邻域NInsert(π，i)进行全面且深入的搜索.

表5 HCEA与6种变体比较结果Table 5 Comparison results of HCEA and its six variants

图7 HCEA与6种变体在不同工件数n下的柱状图Fig.7 Histogram of HCEA and six variants under different n

为了验证HCEA求解NFSSP-SDSTs-RTs的有效性，本节将HCEA与重要国际期刊中的3种算法进行性能比较，即文献[32]中的有效EDA(effective EDA，EEDA)、文献[33]中的迭代贪婪(iterated greedy，IG)算法以及文献[14]中的基于禁忌机制的增强迭代贪婪(TMIIG)算法.EEDA[32]利用二维矩阵学习工件块信息，建立概率模型对解空间中优质解分布进行合理估计，通过调节学习速率更新概率模型，采样概率模型产生新种群，从宏观角度引导搜索方向，目前已被广泛应用于求解各类车间调度问题.IG算法[33]融合启发式算法的高效性和模拟退火算法的优越性，是近年来求解各类流水车间调度问题公认的有效算法，其在求解基于排序模型的FSSP SDSTs方面相比其他智能优化算法有明显优势.TMIIG算法[14]是对经典IG算法的改进，采用改进NEH启发式方法产生高质量的初始种群，并通过基于多种邻域结构的变邻域搜索方法对局部区域进行细致搜索，是近年来求解NFSSP的有效算法.其中，EEDA的参数设置与文献[32]保持一致，具体取值为:种群规模为popsize=150，优势个体数为η=10，学习速率为α=0.1;IG算法的参数设置与文献[33]保持一致，具体取值为:温度常数为T=0.5，移出工件数为d=4;TMIIG算法的参数设置与文献[14]保持一致，具体取值为:温度常数为T0=0.6，移出工件数为d=6，禁忌表长度为ML=1.对于不同规模问题，各算法在运行时间参数ρ=2，4，6下的统计结果如表6-8所示.为能清楚展现算法的求解性能，表6-8中每个问题对应的最优结果用粗体表示.

由表6-8可知，HCEA求解不同规模问题的统计结果明显优于其他3种比较算法，这表明加入快速邻域搜索后算法的性能得以进一步提升.对于绝大部分问题，HCEA在ρ=2时获得的结果均比EEDA，IG和TMIIG算法在ρ=4和ρ=6时获得的结果都要好.对于不同规模的问题，HCEA得到的方差SD值也较小，这说明了HCEA具有良好的求解性能.

为了进一步分析HCEA与其他算法的性能差异是否显著，采用Tukey’s HSD方法对表6-8中的4种算法分别在ρ=2，4，6三种取值下求解各种规模的问题所得到的ARI值进行方差分析(ANOVA).图8为4种算法在不同运行时间下的均值及95%置信度下Tukey’s HSD检验的置信区间图.由图8可知，在不同运行时间下HCEA得到的ARI值与其他3种对比算法所得到的ARI值相比存在显著差异，具有统计学意义.此外，从图8中还可看出，随着运行时间的增加HCEA得到的ARI均值无明显的变化，即HCEA在3种不同运行时间下的置信区间明显存在重叠，表明HCEA性能稳定并且能够在较短的时间内获得相对较优的结果.图9为4种算法和3种运行时间交互作用对ARI值的交互作用图.图中的折线无交叉，表明交互作用不显著.

图9 各算法在不同时间参数下的交互效应图Fig.9 Interaction plot between algorithms and time factors

表6 算法比较结果(ρ=2)Table 6 Comparison results of algorithms(ρ=2)

表7 算法比较结果(ρ=4)Table 7 Comparison results of algorithms(ρ=4)

表8 算法比较结果(ρ=6)Table 8 Comparison results of algorithms(ρ=6)

图8 各算法在不同时间参数下的均值及95%置信度下Tukey’s HSD检验的置信区间图Fig.8 Means plot and 95% Tukey’s HSD confidence intervals for the interaction between the algorithms and the time factors

为了分析算法的收敛特性，HCEA和3种对比算法在ρ=2时求解70×20规模问题的收敛曲线如图10所示.从图10中可明显看出，HCEA能够快速收敛到较优解，再次验证了第3.3节中提出的快速邻域搜索的有效性.在比较算法中，IG和TMIIG虽然利用启发式方法产生高质量的初始解，由于算法并未考虑排序模型解的工件块分布(见第3.1节)，在搜索过程中“移除”和“重构”操作破坏了优质解的结构特征，使得算法仅在初始解附近反复探索，难以搜索到全局意义上的最优解[34].EEDA基于二维概率模型从宏观角度引导全局搜索，可一定程度上避免IG和TMIIG 中存在的优质解模式破坏的问题，但二维概率模型不能准确保留工件块在解序列中的具体位置信息.同时EEDA的局部搜索邻域数量有限且搜索方式单一，导致算法实际搜索深度有限.HCEA利用交叉熵算法更加合理地学习并积累优质解的工件块信息，同时采用带两种搜索策略的快速局部搜索，有利于引导算法在多个子邻域中进行更为深入且高效地搜索，丰富搜索行为的同时加速搜索过程，进而能发现存在于复杂解空间中的优质解.因此，HCEA能在上述实验中表现出更好的性能，是一种求解NFSSPSDSTs RTs的有效且鲁棒的算法.

图10 收敛曲线图Fig.10 Convergence curves

4.5 案例分析

钢铁制造业是发展国民经济与国防基建的支柱性产业，也是资源高度集中且能源密集型工业.综合考虑技术指标与工艺要求，钢材的生产从开始加工到最终完工期间需要连续不断地进行处理[18].例如，进入轧钢机压轧的连铸坯或钢坯其温度必须预先加热到特定范围.如果由于某种原因在压轧之前出现了较长时间的等待，则钢坯的温度将急剧下降.一旦钢坯温度下降到压轧的临界温度以下时，就需要再次对钢坯预热.这样反复加热将会消耗大量能源，造成不必要的资源浪费[4].钢铁厂中经过炉外精炼处理得到的钢水仅是钢材生产过程中的中间产物，还需要再经过铸块(ingots)、脱模(unmolding)、加热(reheating)、初轧(rougher rolling)和精轧(finish rolling)等多道关键性工序后才能够成为合格的钢材.铸块即液态钢水浇注到模具中，冷却得到特定的外形和厚度的钢坯.脱模指在传送带上将钢坯从模具内取出的操作.高炉加热指将钢坯加热至800度以上(达到金属的相变温度)，随后快速降温冷却以提升钢坯硬度，然后再次加热至略低于相变温度(回火操作)，以降低钢坯硬度并提高韧性，为压轧作准备.初轧过程指在轧机上把钢坯初步轧制成板坯、方坯或异型坯的过程.精轧过程则是将初轧后的坯件经过辊道送入精轧机，最终轧制成满足用户要求的型材.图11为钢铁厂中钢坯生产过程的工艺流程.

图11 钢板加工工艺流程Fig.11 The process flow chart of steel slab

为了满足实际的生产要求，需要考虑如下因素:1)钢坯的加工工艺路线相同且每道工序在一台加工设备上进行处理;2)钢坯的加工工序之间满足零等待约束[4-18];3)由于加工件的更换，加工设备需要一定的设置时间进行更换夹具或装卸工件等必要操作，该设置时间取决于加工件的特征[6-7];4)加工件通过中间传送装置运输到加工设备上需要一定的时间，即工件存在释放时间[8];5)钢铁厂根据不同客户的订单要求进行生产，为了提供优质服务，在生产过程中需要满足准时制生产要求[9].因此，钢材加工成形的调度问题可建模为NFSSP-SDSTs RTs.

本节在第4.1-4.4节的基础上，采用云南某钢铁公司实际生产过程的数据，经过预处理后进一步对HCEA的性能进行验证.目前该厂车间内的生产调度方案是由公司里有经验的调度员按照常规的先到先加工原则进行计划排程，即根据工件的释放时间从小到大顺序依次进行加工处理.为了按期生产客户订单所需要的钢材，该厂车间内现有一条生产流水线，正常运转的设备共有5台，即机器数m=5.通过对客户订单汇总整理，需要加工的坯件数量n的数值分别从10到200.针对该厂的实际生产需求，分别使用第4.4 节中的4 种算法:EEDA[32]、IG[33]、TMIIG[14]和HCEA进行求解，各算法分别独立重复30次，运行时间设置为4×n×m毫秒，求解结果如表9所示.由表9可看出，与IG相比，TMIIG和EEDA是非常具有竞争潜力的算法，HCEA在所有实际问题上的平均表现都明显优于EEDA、IG和TMIIG这3种算法，HCEA得到的平均值优于其他算法的平均值，再次验证了HCEA具有良好的求解性能.为了更加直观的进行比较，4种算法的ARI值的箱线图如图12所示.因为篇幅所限，仅列出4个实例对应的结果.箱线图越窄说明方差越小，算法稳定性越高.在绝大多数实例中，HCEA对应的箱线图的宽度最窄，表明鲁棒性好.HCEA的ARI值中线优于其他3种算法，说明该算法得到的解质量较好，HCEA能在相对有限的时间内获得较为满意的调度方案.图13和图14分别给出了HCEA解决实际钢铁厂中钢材生产调度问题Case 10-5和Case 20-5的甘特图.图13和图14中Case 10-5和Case 20-5 最优调度的最大完工时间分别为15143和23341，订单的总提前和延迟时间分别为11501和32874.

表9 实例仿真结果Table 9 Comparison results of algorithms for the real world instances

图12 HCEA与3种对比算法的箱线图Fig.12 Box plot of HCEA and three compared algorithms

图13 Case 105调度方案的甘特图Fig.13 The Gantt chart for Case 105

图14 Case 205调度方案的甘特图Fig.14 The Gantt chart for Case 20-5

综上可知，HCEA能够建立有效的概率模型学习优质解的结构信息，推动种群到达较多的优质解区域.通过基于Interchange和Insert的快速局部搜索方法，对全局搜索获得的优质区域进行较为细致的局部搜索，使算法在全局搜索和局部搜索之间达到了较好的平衡，从而增强了算法的整体搜索能力.通过不同测试问题上的大量实验和比较分析，验证了HCEA是求解NFSSP-SDSTs-RTs的有效算法.

5 结论

本文在传统零等待流水线调度问题的基础上，进一步考虑实际生产中普遍存在的机器设置时间和工件释放时间以及准时制的生产要求，研究工作更加具有实际意义和理论价值.本文以同时最小化总提前和延迟时间为优化目标，建立带有序相关设置时间和释放时间的零等待流水线调度问题(NFSSP-SDSTs-RTs)的排序模型，并提出一种混合交叉熵算法(HCEA)进行求解.主要贡献包括:1)分析问题的性质，设计一种快速评价方法，有效降低评价解的计算复杂度.2)利用交叉熵算法对优质解中工件块的分布进行有效地估计，学习并积累优质解的结构特征信息，通过合理的采样与更新方法，实现对解空间中优质区域的全局搜索.3)提出基于Interchange和Insert操作的局部搜索方法，对全局发现的优质区域进行较为深入的搜索.4)设计基于两种邻域搜索策略的快速局部搜索方法，有效降低了搜索的计算复杂度同时也提高了搜索效率，实现对解空间中更多优质区域的细致搜索.大量实验和算法对比验证了本文所提HCEA是求解NFSSP-SDSTs-RTs的一种有效且鲁棒的算法.

分布式调度是生产调度领域近年来研究的热点，节能降耗是实现低碳制造备受关注的焦点.未来的研究工作针对问题的特性，设计高效智能优化算法，将HCEA扩展应用于求解更加复杂且实际的分布式低碳流水线调度问题.