超启发式交叉熵算法求解分布式装配柔性作业车间调度问题

罗文冲 钱 斌 胡 蓉张长胜向凤红

(1.昆明理工大学信息工程与自动化学院,云南昆明 650500;2.昆明理工大学云南省人工智能重点实验室,云南昆明 650500)

1 引言

生产调度是智能制造中的重要环节,是企业优化资源配置和提高经济效益的关键.近年来,随着经济全球化的深入,许多制造企业积极地参与进国际和区域间的合作.由于产地原料分布不同,产品的制造往往需要交由不同区域的生产基地加工完成,原有的单工厂模式正在被多工厂模式取代[1].这种多工厂模式下的调度问题也被称为分布式调度问题[2].与此同时,为了应对市场的快速变化和产品需求的频繁多变,对于电子行业、服装业等新型制造企业,已经实现了从“短时间内快速生产相同产品”到“短时间内快速生产不同产品”的转变.如何在最短时间内完成不同订单的原料加工、产品装配已经是当下新型制造企业发展的关键.因此,分布式装配柔性作业车间调度问题(distributed assembly flexible job shop scheduling problem,DAFJSP)被提出并逐渐引起研究者的关注.在DAFJSP中,加工阶段每一个加工工厂都是一个独立的柔性作业车间,多个加工工厂需要协作来完成给定的生产任务;加工完成后的工件经运输到达装配工厂,进一步将工件装配成为产品.显然,DAFJSP属于NP-hard问题,研究其优化求解方法具有重要的学术价值;同时,相较于传统的车间调度问题,DAFJSP更加符合实际的生产状况,其研究成果具有重要的应用价值.

与单工厂调度不同,分布式调度问题的求解是从资源配置的总体层面进行考虑,通过各工厂间的协作以实现经济效益的最大化.分布式柔性作业车间调度问题(distributed flexible job shop scheduling problem,DFJSP)是柔性作业车间调度问题和分布式问题的结合.在DFJSP的研究上,Giovanni等[3]以最小化最大完工时间(makespan)为优化目标,提出了一种改进遗传算法(improved genetic algorithm,IGA)进行求解,在算法中设计了4种变异算子来增强算法的搜索性能.Liu等[4]提出了一种将概率融入编码的混合遗传算法,该算法可以有效地减少染色体的长度从而减小解空间.Li等[5]提出了一种基于混合帕累托解的禁忌搜索算法,用来求解同时带有4个目标的DFJSP.Meng等[6]以最小化makespan为优化目标,提出了4种混合整数线性规划(mixed integer linear programming,MILP)模型和一种约束规划(CP)模型.Wu等[7]以makespan为优化目标,提出了一种改进人工蜂群算法(improved artificial bee colony algorithm,IABC),该算法采用工序排序、工厂选择和机器选择的三维向量编码方案,并结合问题特点针对性地设计了多种策略用于种群初始化.

在实际生产中,制造过程往往由生产和装配2个阶段构成,这类两阶段的集成调度问题也引起了学者们的广泛关注.Hatami等[8]首次提出了分布式装配置换流水线调度问题(distributed assembly permutation floe shop scheduling problem,DAPFSP),以最小化makespan为优化目标建立MILP模型,并提出12种基于不同结构的启发式方法和一种变邻域下降(VND)算法进行求解.Sang等[9]为求解以最小化makespan为优化目标的带总流经时间标准的DAPFSP,提出了一种包含产序列和工件序列的编解码方法,设计了3种入侵杂草优化(invasive weed optimization,IWO)算法进行求解.Deng等[10]提出一种考虑同步性和准时性的3阶段装配集成调度问题(three-stage assembly integrated scheduling problem considering assembly synchronization and delivery punctuality,3sAISPSP),以加工–运输–装配同步性和交货准时性的加权和为优化目标,设计一种混合分布估计算法进行求解.Lei等[11]提出了一种考虑序相关的分布式两阶段混合流水车间调度问题(distributed two-stage hybrid flow shop scheduling problem,DTHFSP),并设计了一种有效的改进蛙跳算法求解该问题.Wu等[12]首次以最小化E/T,生产成本为目标建立了分布式装配柔性作业车间调度问题DAFJSP的MILP模型,采用平衡调度机制来权衡求解目标,并提出了一种改进的差分进化模拟退火算法,有效地解决了DAFJSP.从上述文献可知,由于两阶段集成调度问题求解复杂度较大,现有研究主要根据问题特点构造启发式方法将离散优化问题转换为连续优化问题,简化算法中进化操作算子的操作.此外,由于DAFJSP的计算复杂度较高,且当前研究并没有针对DAFJSP的问题特点设计良好的启发式方法.因此,构造合适的启发式方法,并设计高效算法求解该问题具有重要意义.

超启发式算法(hyper-heuristic algorithm,HHA)是一种新型智能优化算法.该算法通过某种高层策略控制多种低层启发式算法(low-level heuristic,LLH)实现对解空间的搜索[13].近年来,超启发式算法已经成功地应用在调度领域.Park等[14]针对工件动态到达和机器故障的作业车间调度问题,以最小化平均加权延迟时间为优化目标,设计一种结合遗传编程的超启发式算法进行求解,该算法低层设计30种调度规则作为备选集合,高层采用遗传编程算法确定新解.Li等[15]设计了一种自适应遗传算法用于求解以最小化最大模糊完工时间为优化目标的模糊柔性作业车间调度问题,该算法高层用自适应遗传算法对不同的6种启发式规则排列进行优化,低层执行高层产生的规则来生成解并执行搜索.此外,为增加算法的深度搜索能力,部分研究者[16–17]在HHA的每个低层启发式操作中加入模拟退火机制,并通过仿真实验验证可该机制可以有效地提高算法性能.HHA通过高层策略来优化低层启发式操作以实现对问题解空间的细致搜索,相较于其他算法,HHA在搜索深度和搜索效率方面具有较大优势.目前大部分研究者采用GA等智能算法作为高层策略,该类超启发式遗传算法虽能较好求解相关问题,但仍存在一定的不足:算法高层仅对规则组合或简单操作组合通过交叉、变异等操作进行优化,一定程度上限制了算法的性能.

交叉熵算法(cross-entropy algorithm,CEA)是一种新型概率估计算法,最早由Rubinstein[18]提出,该算法用于估计复杂随机网络中稀有事件发生的概率.CEA基于马尔科夫方法和采样技术,通过统计学习和更新,使最优解出现的概率逐渐增大并最终趋向于1.目前,已有部分学者将CEA应用到组合优化问题中.Wang等[19]提出了一种改进交叉熵算法(ICEA)求解炼钢–连铸问题,验证了CEA具有较强的鲁棒性和较好的全局搜索能力;同时也发现CEA与遗传算法(GA)和蚁群算法(ACO)等常规搜索算法不同,其在寻找最优解的过程中,逐渐增加的概率可以引导搜索过程往最有希望的方向寻找.She等[20]提出一种混合交叉熵算法求解以最小化模糊最大完工时间和模糊总能耗为优化目标的模糊分布式装配流水线低碳调度问题,该算法在用CEA全局搜索的基础上,设计了一种自适应变邻域局部搜索方法以实现对解空间不同区域的有效搜索.CEA具有较好的全局搜索能力,在执行过程中可以保留当前优质解的结构信息,指引算法的搜索导向,进而加速算法的收敛速度.因此,本文采用CEA作为HHA的高层策略,通过CEA学习和积累高层策略域中优势个体的位置信息,继而在生成新个体时能较准确地确定优势个体在解中所放置的具体位置,更合理地引导低层算法的搜索方向,提高搜索效率.

综上,本文研究DAFJSP的建模和求解.在建模方面,针对实际生产活动中广泛存在的分布式生产和装配两阶段过程,建立以最小化makespan为优化目标的DAFJSP模型;在求解方面,由于DAFJSP这类问题解空间巨大且复杂,常规智能算法难以实现有效搜索,故本文根据DAFJSP的特点,构造4种启发式方法得到较优的初始解,有效地提高解的质量,节约搜索时间;同时,设计了一种融合HHA和CEA的超启发式交叉熵算法(hyper-heuristic cross-entropy algorithm,HHCEA)进行求解.首先,在算法的编码和解码阶段,设计了一种基于工序序列、工厂分配和产品序列的三维向量编码规则和一种基于贪婪规则的解码规则,并使用所提的启发式方法有效提高解的质量;然后,使用设计的HHCEA 对DAFJSP 进行求解.HHCEA的高层采用CEA学习和积累优质排列的信息,其中各排列由结合问题特点设计的11种启发式操作(即11种有效的邻域操作)构成;低层为增加在解空间中的搜索深度,将高层确定的每个排列中的启发式操作依次重复执行指定次数并在执行过程中加入基于模拟退火的扰动机制,以此作为一种新的启发式方法实现对问题解空间中不同区域的细致搜索.最后,通过仿真实验和算法比较验证了所提算法的有效性.

2 分布式装配柔性作业车间调度问题

2.1 符号定义

表1 符号表Table 1 Notations

2.2 问题描述

DAFJSP是一个多工厂协同调度问题,主要分为加工阶段和装配阶段,如图1所示.加工阶段由一系列加工工厂组成,每一个加工工厂都是一个柔性加工工厂.装配阶段由一条装配线构成,可以在每个产品的所有工件加工完成后且装配机处于空闲时被执行.

图1 分布式装配柔性作业车间系统示意图Fig.1 DAFJSP’s system diagram

本文研究的DAFJSP可被描述如下:产品集合P={P1,P2,···,PQ}中每一个产品都根据原料清单由一系列工件装配而成.工件集合J={J1,J2,···,Jn}中的工件需要被分配到F个处于不同地理位置的加工工厂中.每一个加工工厂中装备有一系列加工机器M={M1,M2,···,Mm}.每个工件Ji包含Ui道加工工序,工序之间遵循工序约束,每道工序均可在多台加工机器上加工.一个工件的所有工序只能在同一个工厂中加工完成,不能转移到其他工厂.

每个产品Ph中包含ωh个工件.并且每个工件只属于一个产品.一旦产品的工件全部加工完成,它们就立即交付给装配机进行装配.在装配阶段,产品Ph在装配机器MA上的装配时间定为.此外,在0时刻,所有工件都是独立的,所有机器不考虑故障且都是可用的.加工时间都是确定的,设置时间和运输时间被认为包括在加工时间中.工件之间没有加工顺序约束,同一工件的工序之间有加工顺序约束.一道工序只能选择一台机器加工,且只能被加工一次.一台机器在同一时刻只能加工一道工序,且不允许抢占.DAFJSP的目标是确定工厂的任务分配,工序的加工顺序,工序的机器选择和装配阶段的产品顺序,从而使整个制造过程中的makespan最小化.

其中:式(1)表示工件在加工阶段的完工时间是最后一道工序的完工时间;式(2)表示同一工件内的工序必须按顺序加工;式(3)表示工序完工时间等于开始加工时间加上实际加工时间;式(4)表示一个工件只能在一个工厂内加工;式(5)表示一台加工机器同一时刻只能加工一道工序;式(6)表示一道工序只能选择一个工厂中的一台机器进行加工;式(7)表示一道工序只能被所选择的机器加工一次;式(8)表示产品Ph在加工阶段的完工时间为最后1个工件的加工完成时间;式(9)表示产品序列中第1个产品的装配完成时间等于加工完成时间加上实际装配时间;式(10)表示产品Ph(h ＞1)必须在其包含的所有工件加工完成后且装配机器MA处于空闲时才可以进行装配;式(11)表示序列π的makespan等于最后一个产品的装配完成时间;式(12)表示DAFJSP的目标是在所有可行调度中寻找makespan最小的调度π*.

3 超启发式交叉熵算法求解DAFJSP

3.1 编码和解码

3.1.1 高层策略域编码和解码

在HHCEA中,首先在高层策略域利用CEA对11种操作的组合顺序(每种顺序为一个高层种群个体)进行优化,从而得到高层策略域新个体;然后,在低层问题域将每个高层策略域个体作为独立的启发式算法,用于对低层相应个体进行变邻域搜索,同时加入模拟退火机制来避免低层搜索陷入局部较小,从而实现多种启发式算法的并行搜索.HHCEA两层示意图如图2所示.

图2 HHCEA的两层结构Fig.2 Two-layer structure of HHCEA

对于高层策略域,CEA种群中的每个个体都由策略集中的11个低层启发式操作构成,个体长度为11,且同一个体中允许出现相同的低层启发式操作.解码高层策略域中的个体时,对每个低层问题域中的解,按从左到右的顺序依次执行个体中不同位置对应的低层启发式操作,每执行一种低层启发式操作,就将得到的新解与旧解进行对比,若新解优于旧解(即新解的目标函数值小于旧解的目标函数值),则用新解替换旧解且不再执行剩余的低层启发式操作,否则继续执行直到高层策略域个体中的操作全部被执行完.高层策略域中个体的目标函数值等于其对应的低层问题域中个体的目标函数值.图3表示一个高层策略域个体.

图3 高层策略域个体示意图Fig.3 Diagram of an individual in the upper layer

3.1.2 低层问题域编码和解码

在低层问题域中,每个个体就是原问题的一种工序排列序列(即原问题的解).高层策略域的一种算法对应更新一个解,因此低层问题域的种群大小与高层策略域相同.

由DAFJSP的问题特性可知,DAFJSP中加工阶段属于柔性加工,包括工序调度和机器选择两部分内容.在本文中,在加工阶段的解码时采用贪婪策略,以此将工序排列序列中的每一个工序放在能加工该工序的所有机器上进行活动化解码[21],然后选择能在当前完工时间最短的那台机器上加工,形成具体调度.该方法与文献[22]中所用的两阶段编码解码相比,编码序列更少,可在选择较优机器的同时充分利用每台机器上的空闲时间间隔,减少加工时间的浪费,进而提高解的质量.

表2给出了一个3台机器3个工件的柔性作业车间工序加工时间的例子,其中O1,1和O1,2代表工件J1的加工工序,O2,1,O2,2和O2,3代表工件J2的加工工序,O3,1,O3,2和O3,3代表工件J3的加工工序,M1,M2和M3为3台加工机器.此外,*表示该工序无法在这一机器上加工.

使用表2的数据可得一个编码个体实例为πO={3,1,3,2,1,2,3,2}.下面以该编码为例阐述本文中柔性加工阶段解码方法的具体过程.解码时按照顺序先加工O3,1,从表2 中可知O3,1可在M1和M3上加工,此时若在M3上加工,则当前最大完工时间为(9),比在M1上得到的当前最大完工时间短,因此选择在M3上加工O3,1.然后加工O1,1,该工序在3台机器上都可以加工,由于O1,1是M1上加工的第一个工序,所以该机器加工得到的当前最大完工时间为(7),同理可得在M2和M3上加工的当前最大完工时间分别为(10)和(20),因此O1,1选择在机器M1上加工.后面的解码过程以此类推,不再赘述,最终得到的甘特图如图4所示.

表2 柔性作业车间工序加工时间Table 2 Processing time of the operations in flexible job shop

图4 贪婪策略下活动化解码得到的甘特图Fig.4 Gantt chart obtained by active resolving codes under greedy strategy

加工阶段采用贪婪策略下活动化解码方法后,DAFJSP包括工厂分配、工序调度和产品排列等3个子问题.因此,DAFJSP的编码可以由3部分组成,分别是工序序列、工厂序列和产品序列.工序序列和工厂序列的长度均等于所有工件中的所有工序数之和,产品序列的长度等于装配阶段的产品数.对应表1中的示例,其中一个可行编码如图5所示,其中产品1由工件1组成,产品2由工件2和工件3组成.

工序序列中每个工件的工序都由该工件的工件索引号表示,某个工件的工件索引号出现第几次就代表是该工件的第几道工序,如图5工序序列中第3个元素为第2次出现的3,代表工件3的第2道工序O3,2;工厂序列由工厂索引号组成,工序序列中的工序在对应的工厂序列中元素代表的工厂中加工,如图5中工序序列的第2位和第4位元素是2,对应的工厂序列的第4位和第6位元素也是2,意味着工件2在工厂2中加工,而工件1和工件3根据图5所示在工厂1中加工;产品序列由产品索引号组成,代表装配阶段产品的装配顺序,如图5中产品序列意味着在装配阶段先装配产品2,再装配产品1.

图5 编码示意图Fig.5 Encoding diagram

初始化时,为避免解的分布过于集中,高层策略域中的每个个体随机产生,但确保每个个体中的低层启发式操作不出现重复.低层问题域中的种群个体均随机生成.

3.2 启发式方法

针对DAFJSP的性质,Wu等[12]提出了一种基于工序和工件的启发式编码解码方法:1)工序按照加工时间升序排列,优先选择位于排列顺序前面的工序进入工厂;2) 工件按照加工时间升序排列,优先选择位于排列顺序前面的工件进入工厂.

为了方便理解,本文给出一个7个工件,2台机器,2个工厂,3个产品的例子.表3给出了工序的加工时间和产品的装配时间.

表3 部分工件加工时间及产品装配时间Table 3 Processing time of the operations and assembly time of the product for the example

其中,产品装配程序为:N1={3,6},N2={1,2,7},N3={4,5}.πA代表一个产品序列,例如πA={1,3,2}是所给例子的一个可能序列.而产品Ph由ωh个工件组成,表示产品Ph所含工序排列,例如:={3,6,3,6},={1,1,7,2,2,7}和={5,4,4,5}.根据πA将所有产品的工序组成组合到,就能能够得到一个完整的工序序列,例如πO={3,6,3,6,5,4,4,5,1,1,7,2,2,7}.为了更好地解释所提方法,对以下符号进行说明:

TOPj:表示工序j在所有可行机器上的加工时间之和.

TJPi:表示工件i包含的所有工序在所有可行机器上的加工时间之和,TJPi,j=

Rh:产品Ph的初始排列.

Js:工件排列.

Fs:Js中工件对应的工厂序号排列.

接下来介绍本文提出的4种启发式方法,分别是:

1) Heuristic1(H1).

最短加工时间(shortest processing time,SPT)规则是一种常见的分配规则.在SPT中,加工时间最短的工件优先获得加工.该规则目的在于减少产品存货、平均加工时间和工件平均延迟时间等.首先用SPT规则确定产品装配序列;然后在每个产品内工序排列过程中,本文借鉴了Framinan等[23]关于启发式方法的思想.H1方法具体步骤为:

步骤1将所有产品按照装配时间升序排列得到πA;转至步骤2.

步骤2对于每个产品Ph,h=1,2,···,Q.计算Ph包含的所有工序的TOPj=,a=1,2,···,uh;将TOP从小到大排列得到Rh,h=1,2,···,Q;令h=1;转至步骤3.

步骤3若Rh ＞1,转至步骤3.1;否则转至步骤3.4.

步 骤3.1取Rh前2位Rh(1)和Rh(2),若Rh(1)(2),则转至步骤3.1.1;否则转至步骤3.1.2.

步骤3.1.1将Rh(1)和Rh(2)排序可得到2个部分调度{Rh(1),Rh(2)}和{Rh(2),Rh(1)},评价这个部分调度,并将其最大完工时间较小的一个作为当前调度,记为={Rh(1),Rh(2)},令a=3;转至步骤3.2.

步 骤3.1.2将Rh(1)和Rh(2)放 在的 前2位,得到当前调度={Rh(1),Rh(2)},令a=3,k=a,Tn=0,取出Rh的第a个工序Rh(a);转至步骤3.2.

步 骤3.2若Rh(k-1)Rh(a),则 将Rh(a)插入到的第k位置,得到一个新排列,令k=k-1;否则,令Tn=Tn+1,k=k-1.重复本步骤,直到k=1,转至步骤3.3.

步骤3.3评价将Rh(a)插入后得到的a-Tn个排列,将完工时间最小的排列作为当前调度;转至步骤3.4.

步骤3.4令a=a+1,若a≤ωh,则转至步骤3;否则转至步骤4.

步骤4令h=h+1,若h≤Q,则转至步骤3;否则转至步骤5.

步骤5根据πA,得到完整的工序排列πO=;转至步骤6.

步骤6计算工件时间TJPi=,i=1,2,···,n;将TJPi从小到大排列得到Js,s=1,2,···,n;令b=1,f=1;转至步骤7.

步骤7令Fs(1)=f;令b=b+1,f=f+1;若f ＞F,则f=1;若b≤n,则重复步骤7;否则转至步骤8.

步骤8根据所得的Js,Fs和πO,按照Js中工件对应Fs中的工厂号,将πO中的工序分配到工厂,得到π;输出当前调度πA,πO和π,算法结束.

便于理解,用上述提到的例子说明H1的过程.具体过程为:

步骤1根据SPT规则得到πA={1,3,2}.

步骤2计算产品Ph包含的所有工序的TOPa=,得到Rh.以产品P2为例:对于产品P2中的所有工序,先比较每个工件第1道工序的加工时间和,分别为C1,1=8,C2,1=8,C7,1=10,其中O1,1加工时间和最短,将工件1放在R2的第1位.接下来计算C1,2=11,在C1,2,C2,1和C7,1中选择最小的C7,1,R2第2位选择工件7,依此类推,可得R2={1,7,7,1,2,2}.

步骤3以R2为例,阐述的构成过程.

步骤3.1选择R2中的前2位{1}和{7},因为它们不是同一个工件,所以将{1}和{7}插入到的前2位,可得到2个排列{1,7}和{7,1},计算这2个排列的完工时间分别是5和5,由于两个排列完工时间相同,我们选择第一个排列(每次出现具有相同完工时间的排列时,统一选择编号在前的排列).

步骤3.2选择R2中的第3位{7}插入到中,此时有{1,7,7}和{7,1,7}2种情况,完工时间都是6,选择第一种{1,7,7}作为当前排列.选择R2中的第4位{1}插入到中,此时中的排列一共有:{1,7,7,1},{1,7,1,7}和{1,1,7,7}3种,完工时间分别为9,8和8,选择={1,7,1,7}.将R2中的第5位{2}插入到中,得到{1,7,1,7,2},{1,7,1,2,7},{1,7,2,1,7},{1,2,7,1,7}和{2,1,7,1,7}5 种,完工时间分别为12,11,10,13和12,则={1,7,2,1,7}.最后将R2中的第6位{2}插入到中,得到{1,7,2,1,7,2},{1,7,2,1,2,7},{1,7,2,2,1,7},{1,2,7,2,1,7}和{2,1,7,2,1,7}5种,完工时间分别为13,12,12,16和16,最终={1,7,2,1,2,7}.

步骤4使用同样的方法,可以得到产品P2和产品P3的工序排列为别为={6,6,3,3}和={5,4,5,4},对应的完工时间分别为13和10.

步骤5根据πA={1,3,2},得到最终的工序序列为πO={6,6,3,3,5,4,5,4,1,7,2,1,2,7}.

步骤6计算各个工件所包含工序在所有机器上的加工时间和,分别是C1=19,C2=23,C3=30,C4=18,C5=18,C6=16和C7=14,按从小到大排得到工件排列Js={7,6,4,5,1,2,3}.

步骤7由工件分配规则将工件分配到工厂,得到Fs={1,2,1,2,1,2,1}.

步骤8根据Js,Fs和πO,得到最终的调度为π={2,2,1,1,2,1,2,1,1,1,2,1,2,1},πA={1,3,2}和πO={6,6,3,3,4,4,5,5,1,7,2,1,2,7}.

使用H1方法解码后最终的makespan为38.所给例子按照H1方法解码后的甘特图如图6所示.

图6 所给例子用H1方法得到的甘特图Fig.6 Gant chart obtained by the Heuristic1 for the example

在H1方法中,对于n×m×F ×Q的问题,根据装配时间,使用快速排序将Q个产品按升序排列得到πA的复杂度为O(QlogQ).计算Q个产品内工序排列的复杂度为O((n×m)2).计算所有工件所包含的工序在所有可行机器上的加工时间之和的复杂度为O(n×m2).使用快速排序将n个工件按升序排列的复杂度为O(nlogn).将n个工件分配到F个工厂中的复杂度为O(n×F).按照复杂度渐进法则,H1方法的复杂度为O(QlogQ)+O((n×m)2)+O(n×F).

2) Heuristic2(H2).

H2依据先加工先完成的规则,优先考虑在加工阶段先加工完成的产品.一个产品在加工阶段的完工时间就是该产品在装配阶段的最早开始时间.H2的过程为:1) 按照H1中计算产品Ph内工序序列的方法求出每个产品内的工序序列及最早开始装配时间;2)首先根据最早开始时间升序排列得到产品序列πA;3) 根据πA得到最终的工序序列πO;4)按照H1中的工厂分配方法将工件分配到工厂,得到工厂序列π.

所给例子中,每个产品的最早开始装配时间分别为E1=13,E2=12和E3=8.因此,πA={3,2,1},最终的工序序列为πO={5,4,5,4,1,7,2,1,2,7,6,6,3,3}.最后再将工件分配给工厂,得到Js={7,6,4,5,1,2,3}和π={2,1,2,1,1,1,2,1,2,1,2,2,1,1}.使用H2方法解码后最终的makespan为32.所给例子按照H2方法解码后的甘特图如图7所示.

图7 所给例子用H2方法得到的甘特图Fig.7 Gant chart obtained by the Heuristic2 for the example

在H2方法中,对于n×m×F ×Q的问题,计算Q个产品内工序排列的复杂度为O((n×m)2).根据产品的加工时间,使用快速排序将Q个产品按升序排列得到πA的复杂度为O(QlogQ).计算所有工件所包含的工序在所有可行机器上的加工时间之和的复杂度为O(n×m2).使用快速排序将n个工件按升序排列的复杂度为O(nlogn).将n个工件分配到F个工厂中的复杂度为O(n×F).按照复杂度渐进法则,H2方法的复杂度为O((n×m)2)+O(QlogQ)+O(n×F).

3) Heuristic3(H3).

H3跟H2类似,区别在于H3中的构造过程为:在满足工件约束的前提下,将产品Ph包含的所有工序按照在加工阶段所有机器上的加工时间之和的升序排列,得到.其余过程跟H2一样.

所给例子根据H3方法可得到:={6,6,3,3},={1,7,7,1,2,2}和={5,4,4,5},对应产品的最早开始装配时间分别为E1=13,E2=18和E3=11.因此πA={3,1,2},最终的工序序列为πO={5,4,4,5,6,3,3,6,1,7,7,1,2,2}.最后将工件分配给工厂,得到Js={7,6,4,5,1,2,3}和π={2,1,1,2,2,1,1,2,1,1,1,1,2,2}.使用H3方法解码后最终的makespan为32.

在H3方法中,对于n×m×F ×Q的问题,计算Q个产品内工序排列的复杂度为O(n×m2).根据产品的加工时间,使用快速排序将Q个产品按升序排列得到πA的复杂度为O(QlogQ).计算所有工件所包含的工序在所有可行机器上的加工时间之和的复杂度为O(n×m2).使用快速排序将n个工件按升序排列的复杂度为O(nlogn).将n个工件分配到F个工厂中的复杂度为O(n×F).按照复杂度渐进法则,H3方法的复杂度为O(n×m2)+O(QlogQ)+O(n×F).

4) Heuristic4(H4).

H4也跟H2类似,区别在于H4中内的工序序列随机生成.

所给例子中,根据所给例子中,根据H4方法得到的一个可行解为:={6,3,3,6},={7,2,7,1,2,1}和={5,5,4,4},对应产品的最早开始装配时间分别为E1=12,E2=16和E3=11.因此πA={3,1,2},最终工序序列为πO={5,5,4,4,6,3,3,6,7,2,7,1,2,1}.最后将工件分配给工厂,得到Js={7,6,4,5,1,2,3}和π={2,2,1,1,2,1,1,2,1,2,1,1,2,1}.使用方法H4解码后最终的makespan为32.

在H4方法中,对于n×m×F ×Q的问题,计算Q个产品内工序排列的复杂度为O(n×m).根据产品的加工时间,使用快速排序将Q个产品按升序排列得到πA的复杂度为O(QlogQ).计算所有工件所包含的工序在所有可行机器上的加工时间之和的复杂度为O(n×m2).使用快速排序将n个工件按升序排列的复杂度为O(nlogn).将n个工件分配到F个工厂中的复杂度为O(n×F).按照复杂度渐进法则,H4方法的复杂度为O(n×m2)+O(QlogQ)+O(n×F).

为了验证本文提出的4种启发式方法(H1,H2,H3和H4)的有效性,使用At2中的24个算例(At2的生成方法在第4章中有详细介绍)对5种启发式方法进行测试.由于启发式方法并不能直接得到最优解,所以采用相对百分偏差(relative percentage deviation,RPD)作为评价准则.RPD的计算过程如下:

其中:最优值OPTbest是本文提出的HHCEA和所有对比算法得到的最优解,ALGSOL是每个启发式方法求解算例后的结果.具体测试方法为:分别使用4种启发式方法对At2中的24个算例均解码50次,使用式(13)计算每个解的RPD,最后计算RPD的平均值(即一个算例由一种启发式方法得到的50个解的RPD的平均值).计算结果的区间图如图8所示.由图8可知,H2相比其他3种启发式方法有明显优势.

图8 4种启发式方法计算At2算例的结果在均值的95%置信区间下的区间图Fig.8 Means plot and 95% Tukey’s HSD confidence intervals for four heuristics and At2 instances

本文中,所有低层问题域中的个体均使用H2进行解码.

3.3 交叉熵算法

使用交叉熵算法对高层策略域种群进行优化,其关键在于随机样本的抽取和概率分布参数的确定.令X=(x1,x2,···,xN)为某个概率空间Ω的一个随机变量,其中X为高层策略域中的一个个体.随机变量xi(i=1,2,···,11)表示策略域个体内的第i位按概率Pi=(pi1,pi2,···,pi12)在11个低层启发式操作中选择的一个操作,且对于任何一个Pi均满足=1.高层策略域个体的评价函数可表示为

其中:Ω是高层策略域个体的解空间,X*是最优个体,γ*为由目标约束S确定的最优个体所对应的评价值.样本X服从概率参数为χ的概率分布,式(14)中对X*和γ*的求解可转化为确定最优概率分布f(X;γ)的概率估计问题:

其中:Pχ和Eχ分别是对分布f(X;χ)的概率度量和期望;I{S(X)≤γ}是对于Ω上不同χ值得示性函数集合,EχI{S(X)≤γ}是样本中目标优于χ的期望值.对l的估计采用重要采样的方法,根据Ω上可以对S(X)接近χ的分布δ进行采样:

为了估计l,此时需要确定概率参数χ使分布f与δ的最优分布δ*差距最小.衡量两个概率分布的Kullback-Leible(K–L)距离以判断他们之间的相似度,最小化K–L距离等价为最小化交叉熵[24].分布p对q的交叉熵定义为q分布的自信息对p分布的期望:

此时,问题转化为求解概率参数χ使得δ*和f(X;χ)的交叉熵最小化问题:

将(16)带入(18)可得到概率参数的推断:

对于χt-1,根据概率分布f(X;χt-1)生成一组随机样本{Y1,Y2,···,YN},样本的分位数为ρ,则γt需满足:

由于样本X由概率矩阵Pil生成,因此分布f(X;χ)可以有以下形式:

根据式(24)确定的概率参数即可更新高层策略域个体,具体描述如下:

步骤1第k次迭代中高层策略域个体编码的概率分布参数为p(k-1),随机生成一定数量的个体,构成本次迭代的一组随机样本{X1,X2,···,XN}.

步骤2由步骤1中产生的随机样本按照式(23)确定本次迭代的概率参数p(k).通过对样本{X1,X2,···,XN}的目标函数值{S(X1),S(X2),···,S(XN)}进行排序,计算样本的ρ分位数,使P(S(X)≤)=ρ.

步骤3对于第k+1次迭代高层策略域个体编码的概率参数,引入平滑系数α,使得

3.4 低层启发式操作

在组合优化问题中,全局最优解在所有邻域结构下均为局部最优解,而接近全局最优解的相似优质解往往在多个邻域结构中作为局部最优解.采用单一的邻域搜索容易较早达到并陷入该邻域结构的局部最优解,但此解一般不是理想解[25].而采用多种有效邻域操作(对应多种邻域结构)的迭代搜索可在到达多种邻域结构共同的局部最优非劣解前一直持续进行搜索,从而可以更快获得全局最优解.此外,DAFJSP作为一个两阶段组合优化调度问题,需要考虑加工阶段多个并行的柔性作业车间和装配阶段产品顺序优化,这使其解空间规模十分庞大和复杂,大量优质局部最优解分散在巨大的“峡谷型”解空间中靠近底部的区域[26],增加了算法发现优质解的难度.因此,本文在算法低层设计了5类共11种有效邻域操作,并动态混合这11种邻域操作构成启发式搜索以增强算法局部搜索的深度.邻域操作具体为:

1) 基于工序的操作.

①LLH1:全局交换操作.在工序序列中随机选择一个元素r,将其替换为一个不等于r的工件号v;然后在工序序列中所有等于v的元素中随机选择一个替换成r.交换参与改变的这两个元素对应的工厂号.

2) 基于产品序列的全局操作.

②LLH2:产品序列交换操作.从产品序列中随机选择两位并进行交换.

③LLH3:产品序列插入操作.从产品序列中随机选择两位,将位置编号大的产品插入到位置编号小的产品之前.

3) 基于工厂的全局操作.

④LLH4:工厂序列交换操作.在所有工件中随机选择两个不相同的工件m和n,将工件m内所包含的所有工序的工厂号替换为工件n所属的工厂号;同理,将工件n内所包含的所有工序的工厂号替换为工件m原来的工厂号.

⑤LLH5:工厂序列变异操作.在所有工件中随机选择一个工件l和工件l对应的工厂号g,然后从可行范围内随机选择一个不等于g的值z,将工件l内所包含的所有工序的工厂号替换为z.

4) 基于加工阶段关键工厂的局部操作.

从解码方法中可以看出,加工阶段在DAFJSP中所占比重较大,且计算过程更加复杂.因此,缩短加工阶段的完工时间将直接影响全局的完工时间.在加工阶段中,由于工件分布不均匀,本地完工时间最大的工厂是其中的关键工厂,具有更大的优化意义.本文设计了4种针对关键工厂的邻域操作.

⑥LLH6:工厂内交换操作.在关键工厂的加工序列随机一个元素x,将其替换为一个不等于x的工件号y;然后在加工序列中所有等于y的元素中随机选择一个替换成x.

⑦LLH7:工厂内插入操作.在关键工厂的加工序列中随机选择两位,将位置编号大的产品插入到位置编号小的产品之前.

⑧LLH8:工厂内相邻交叉操作.在关键工厂内随机选择一个位置,将与它相邻的前一个或后一个与它交换.

5) 基于关键产品的操作.

从解码方法中可以看出,产品内的工序排列在整体解码中具有重要影响,同时,装配阶段的开始时间取决于第一个产品的完工时间.缩短第一个产品的加工时间将直接缩短全局的完工时间,显然,第一个加工的产品就是整个调度的关键产品.

⑨LLH9:产品内交换操作.在关键产品的加工序列随机一个元素a,将其替换为一个不等于a的工件号b;然后在加工序列中所有等于b的元素中随机选择一个替换成a.

⑩LLH10:产品内插入操作.在关键产品的加工序列中随机选择两位,将位置编号大的产品插入到位置编号小的产品之前.

在解码过程中,如果出现某一产品的开始加工时间大于前一个产品的装配完成时间,则意味着缩短该产品的加工时间将缩短全局的完工时间.对该产品进行产品内交换操作,实施步骤参考LLH9.

为了进一步提高算法在每个低层启发式操作中跳出局部最优的能力,本文在算法中加入了模拟退火机制.具体而言,在算法执行过程中,设置冷却控制表及初始温度T0,每次迭代时根据冷却系数λ来降温,温度变化式为:

其中:Titer+1为下一次迭代过程中所用的温度,Titer为当前温度.温度通过下式来影响低层启发式操作:

式中Δf为每次新解与旧解的目标函数值之差.迭代时,每个新算法中的每种低层启发式操作,都更新对应的解10次,每次得到的新解若优于旧解,则将旧解用新解替换;否则生成一个0到1之间的随机数r,若r小于P则用这个差解将旧解替换.由此,在迭代后期每一个低层启发式操作将会退化为贪婪搜索,从而使算法最终收敛到全局最优值.

3.5 算法流程

根据上述算法描述,整个算法流程图如图9所示.

图9 HHCEA流程图Fig.9 Flow chart of HHCEA

具体步骤如下:

步骤1初始化高层策略域和低层问题域.种群大小均为ps.第1次迭代时,高层策略域随机生成长度为11的不重复个体,确保每个低层启发式操作在开始时均被执行.

步骤2设置模拟退火初始温度T0=200.

步骤3评价低层问题域种群.

步骤4将高层策略域个体中的启发式操作依次对应更新低层问题域个体,根据目标函数计算得到解的适应值.每次更新时,若产生的新解的适应值优于旧解,则用新解替换旧解,并不再执行后续低层启发式操作;否则根据当前温度以一定概率接收差解.如果接受了当前差解,则判断是否达到了该操作的内循环次数10,没达到则继续执行,否则结束该操作.执行完一个高层策略域个体后,该高层个体的适应值即为对应的低层个体的适应值.

步骤5使用交叉熵方法更新高层策略域种群.

步骤6根据冷却系数降温,并判断是否满足终止条件.若不满足则跳转到步骤3,若满足则终止循环.

3.6 算法复杂度分析

智能算法的复杂度计算为执行加、减、乘、除的次数,HHCEA中包含5个参数,分别为种群容量ps、统计样本的ρ分位数、迭代次数Gmax和高层策略域个体长度L,对于问题规模为n×m×F ×Q的问题,所提算法的计算复杂度分析如下:

算法初始生成高层问题域个体时,其计算复杂度为O(L!),考虑种群规模,产生高层策略域种群的计算复杂度为O(ps×L!).记H2方法计算低层问题域个体的计算复杂度为O(Θ)(其中,Θ=O((n×m)2)+O(QlogQ)+O(n×F)),考虑种群规模,低层问题域种群的计算复杂度为ps×Θ.使用CEA更新高层策略域个体的计算复杂度为O((ps×ρ×L2)/100),对种群中个体采用快速排序法排序的计算复杂度为O(pslogps),低层启发式操作排列构成的启发式方法的计算复杂度为O(10×L×Θ).综上,整个算法的计算复杂度为

4 实验结果及分析

4.1 算例生成

当前并没有标准测试算例可用来测试DAFJSP,故本文所有测试算例均采用随机方式生成.由于在目前研究带有装配阶段的两阶段集成调度问题的文献中,产品装配时间的设定差异较大,且产品的装配时间长短对问题的结果具有重要影响.因此本文根据产品的装配时间设计At1,At2和At3三组不同类型的测试算例,每组类型均包含24个测试算例,其中n={10,30,50},m={3,5},F={3,5},Q={3,5}按照n×m×F ×Q的形式组合.在At1中,工件工序的加工时间和产品的装配时间分别取在[10,30]和[10,20]中按均匀分布随机产生的整数;在At2中,工件工序的加工时间和产品的装配时间分别取在[10,30]和[20,40]中按均匀分布随机产生的整数;在At3中,工件工序的加工时间和产品的装配时间分别取在[10,30]和[40,80]中按均匀分布随机产生的整数1本文仿真实验所用测试算例可在https://github.com/LuoWenchong/test-data.git下载..本文中所有算法的测试程序均由Delphi10.3编程实现,Windows10 64位操作系统,处理器为Intel Core i7-10700@2.9Ghz,内存16 GB.每种算法对每个测试问题均在相同时间下独立运行20次.其中,AVG为算法独立运行20次输出最优结果的平均值,BST为算法独立运行20次输出结果的最优值,WST为算法独立运行20次输出结果的最差值.

4.2 参数设定

在HHCEA中,关键参数主要有种群规模ps、统计样本的ρ分位数、高层策略域中概率参数的平滑系数α和模拟退火机制中的冷却系数λ.为了测试参数设置对算法效果的影响,采用DOE方法[27]分析参数对于算法性能的影响.参与测试的4个参数的水平取值如表4所示.

表4 主要参数水平表Table 4 Main parameters and levels

根据各参数的水平设置,使用正交表L16(44)进行实验,共有16组不同参数不同水平的组合.使用At2中的算例30×3×3×5,对每种不同参数组合的HHCEA独立运行20次,运行时间均为n×ms,每组得到20个最大完工时间,取它们的平均值作为响应值,实验结果越小代表这组参数下的算法性能越强,参数设置的正交表如表5所示.

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

根据上述实验,得到4个主要参数的平均响应值和影响力等级如表6所示,其中等级一栏中数值越小代表该参数的影响力排名越高.由表6可得图10的参数响应趋势.

表6 各参数的平均响应值和等级Table 6 Average response values at different levels of each parameter

通过表5–6和图10可知,当HHCEA的参数设置ps=40,ρ=20,α=0.15,λ=0.85时,算法性能表现良好.

图10 各参数影响趋势Fig.10 The influence trend of each parameter

4.3 算法性能对比

对于不同规模的算例,设定每种算法的运行时间均为n×ms,每种算例均独立运行20次.各问题性能指标的最优解用粗体表示.

4.3.1 高层策略有效性验证

为验证HHCEA中选择CEA作为高层策略及在算法中嵌入模拟退火机制的有效性,将HHCEA与采用GA为高层策略的超启发式遗传算法(hyper-heuristic genetic algorithm,HHGA)和没有嵌入模拟退火机制的HHCEA_noSA进行比较.以At2中的24个算例为测试问题,每种算法的其他参数均保持一致,结果如表7所示.由表7可知,HHCEA在绝大部分算例上的测试结果都优于HHGA和HHCEA_noSA,证明了嵌入模拟退火机制的CEA作为高层策略的有效性.实际上,CEA作为基于概率模型的方法,可在一定程度上避免GA这类传统进化方法普遍存在的较优解模式破坏的问题[28].而嵌入模拟退火机制的HHCEA在不同规模的算例中,性能均好于或等于HHCEA_noSA,而嵌入模拟退火机制的HHGA性能甚至略好于HHCEA_noSA,这也验证了这一机制可以有效地避免低层启发式操作陷入局部最优,增强了低层启发式操作的局部搜索能力.

表7 验证HHCEA高层有效性Table 7 Verify upper level’s effectiveness of HHCEA

4.3.2 本文算法有效性验证

由于当前对于DAFJSP的研究较少,没有找到合适的求解单目标DAFJSP的算法.因此,本文选择了求解跟DAFJSP相近的DFJSP的单目标算法IGA[3],IABC[7]和求解DAPFSP的有效算法IWO[9]作为对比实验,验证HHCEA的有效性.

由于IGA和IABC主要用于求解DFJSP,出于公平原则,先将HHCEA同IGA和IABC一起求解DFJSP.在求解DFJSP中,IGA和IABC的参数设置均同所属文献一致,HHCEA采用IGA中所给解码方式.测试算例采用n={10,30,50},m={3,5}和F={3,5}按照n×m×F的形式组合成的12个算例,其中每个工序的加工时间取在[10,30]中随机产生的整数.每种算法的运行时间均为n×m×0.02 s,每种算例均独立运行20次.各算法比较结果如表8所示.由表8可知,HHCEA在不同规模的算例中,性能均好于或等于IGA和IABC,证明了HHCEA在求解DFJSP中同样是一个有效算法.

表8 各算法求解DFJSP结果Table 8 Comparison results of each algorithm for the DFJSP

在求解DAFJSP的算法比较中,对IGA,IABC和IWO进行适量调整.在IGA中的变异阶段,增加操作LLH9作为对产品序列的变异算子.在IABC中,将所有对机器的操作改为随产品序列的操作.在IWO的局部搜索阶段,增加LLH9作为产品序列的操作算子.各算法的解码方式均采用本文所提的H2方法.各算法在不同装配时间下的运算结果如表9所示,表9中数据表示各算法在当前算例下独立运行20次后得到的平均值.由表9可知,HHCEA在绝大部分算例中,性能优于或等于其他3种算法,证明了HHCEA是求解DAFJSP的有效算法.

表9 各算法求解不同装配时间下的DAFJSP结果Table 9 Comparison results of each algorithm for the DAFJSP under different assembly time

图11为HHCEA求解算例10×3×3×3在不同装配时间下的甘特图.由图11可以看出,由于装配阶段作为DAFJSP的最后阶段,不同的装配时间对整个问题影响较大,若装配时间选取过大,会导致装配阶段所占比重太大,从而削弱整个问题的优化过程,不容易体现出算法性能的差异.

图11 HHCEA求解算例10×3×3×3在不同装配时间下的甘特图Fig.11 Gant chart of HHCEA for solving 10×3×3×3 instance under different assembly time

图11为HHCEA求解算例10×3×3×3在不同装配时间下的甘特图.由图11可以看出,由于装配阶段作为DAFJSP的最后阶段,不同的装配时间对整个问题影响较大,若装配时间选取过大,会导致装配阶段所占比重太大,从而削弱整个问题的优化过程,不容易体现出算法性能的差异.

为了进一步验证HHCEA与其他算法差异性的显著与否,根据3种装配时间下的算法全局对比结果,采用Tukey’s HSD方法对表9中HHCEA等4种对比算法的AVG进行多因素方差分析(ANOVA),受控因素为算法种类和不同的装配时间.图12显示了3种算法在不同装配时间下的均值变化线及95%置信度下Tukey’s HSD检验的置信区间.从图12中可以看出,HHCEA在均值水平上与其他算法均有显著差异,证明了HHCEA能够有效求解DAFJSP.

根据图12及表9中的数据可知,在求解复杂度越大的算例中,HHCEA所表现出的竞争力和有效性就越强.实际上,如果算法中用于搜索的邻域操作(如交叉、变异、插入和交换)数量越多,则各邻域的共同局部最优值越小,算法可在到达此共同局部最优值之前一直持续向下搜索.现有的群智能算法(如本文比较的IGA,IABC和IWO)本质上都是采用数个邻域操作,按某种基本固定的顺序反复执行(执行次数为算法进化代数),以实现对问题解空间的搜索.由于这类算法搜索邻域数量有限且搜索方式单一,导致其实际搜索深度有限.HHCEA高层利用概率模型动态控制低层个体来执行包含多种邻域搜索组合的启发式搜索,有利于引导算法在多种不同邻域区域中进行更深入且有效的搜索,进而能发现存在于复杂解空间中的优质解.因此,HHCEA能在上述实验中取得更好的结果.

图12 各算法在不同装配时间下的均值线及95%置信度下的Turkey’HSD检验的置信区间Fig.12 Means plot and 95% Tukey’s HSD confidence intervals for the interaction between the algorithms and the different assembly time

5 结论

本文针对生产过程中广泛存在的分布式装配柔性作业车间调度问题DAFJSP,以最小化makespan为目标建立其问题模型.由于DAFJSP问题复杂度较大,构造了4种启发式算法来减小解空间,其中,所提的H2方法能够有效提高解的质量.为求解DAFJSP,设计了一种超启发式交叉熵算法HHCEA;在HHCEA低层,设计了11种启发式操作,增强算法执行较深入局部搜索的能力;在HHCEA高层采用CEA持续更新种群,进而动态确定低层每个个体所执行的启发式算法,引导算法较快到达解空间中的多个优质区域执行搜索.后续工作将在DAFJSP中考虑加工过程和客户需求的不确定性,进一步研究不确定处理技术的有效算法设计.