唐红涛，杨志鹏，刘家毅

(武汉理工大学 机电工程学院，湖北 武汉 430070)

批调度[1](batch scheduling)作为车间调度的一个重要分支，是智能制造领域研究的热点和难点问题之一。不同于传统调度问题中机器每次只能加工一个工件，批处理机在不超过其容量限制的前提下可以同时加工多个工件。而带有工序并行性的批调度问题，约束更加复杂，求解难度更大，在工程上该问题主要存在于半导体企业和铸造企业，对该问题的研究能够更好地指导企业生产，提高生产效率。

针对批调度和并行工序[2]的研究出现了众多研究成果。Ikura等[3]研究了工件大小相同，到达时间和截止日期一致的单机批调度问题，并提出了一种多项式算法来优化最大完工时间。在铸造车间调度问题上，Stawowy等[4]研究了造型车间和制芯车间协同生产与调度的问题,并提出了一种拉格朗日启发式算法来求解不同车间协调生产的问题，但其忽略了造型工序与制芯可并行的情形；胡常伟等[5]研究了铸造车间工件同材质情况下如何进行分批与调度的问题。在工序可并行加工的问题上，刘晓平等[6]通过给并行工序设定优先级来优化最大完工时间。在差异工件批次划分问题上，黄锦钿等[7]采用启发式规则对差异工件进行分批以优化总加权拖期；徐本柱等[8]采用试探法解决了批调度过程中工件分批无规则的问题；王万良等[9]提出了一种基于产品需求量的批次划分方案来解决工件盲目分批的问题，并采用差分进化算法优化总生产成本。这些学者的研究丰富了车间调度理论，其研究成果对铸造车间实际生产有一定的指导意义。

国内外对铸造车间批调度问题已经有了较丰富的研究成果，但由于实际铸造车间生产比较复杂[10]，一般都会进行近似化处理。如文献[4]忽视了工序之间的并行性，且工件组批过程中不存在约束，皆不能较好地应用于国内实际铸造生产车间。国内铸造行业以单件小批量生产模式为主导，具有产品种类与材质多样、生产周期长、工艺约束多、排产紊乱等特点，本文以铸造行业为背景，结合铸造车间实际生产过程中存在不相容工件族约束、沙箱尺寸约束、熔炼重量约束，构建一种混合整数规划模型，并设计一种带有单工序编码与批首次匹配规则解码的改进和声算法对模型求解。最后根据企业实际生产数据进行案例仿真，验证本文研究的有效性。

1 问题描述及数学模型构建

1.1 问题描述

铸造车间按先组批后分配过程进行生产。生产阶段由班组完成，班组类似于批调度中的批处理机，因此本文将班组虚拟化为批处理机，将班组加工时间虚拟化为批处理机加工时间进行研究。如图1所示，按工艺流程主要分为工件组批、沙箱选择、任务分配3个阶段。不同阶段存在不同约束，主要有以下3种。

图1 铸造车间前工段产品工艺流程图Figure 1 The flowchart of products in previous-stage casting workshop

1) 材质约束。在工件组批阶段，同一批工件在同一熔炼炉进行熔炼，因此工件的材质必须相同。

2) 尺寸约束。在沙箱选择阶段，同一批工件的尺寸和不能超出所选沙箱的最大尺寸。

3) 重量约束。在工件组批阶段，同一批工件的重量和不能超出熔炼炉单次最大熔炼重量。

传统批调度问题可以描述为：n个工件被分为k个不同的批次，在m台批处理机(加工资源)上进行加工，每台批处理机同一时间段只能加工一个批次的一道工序，需要确定不同工件加工顺序，并为每个工件选择加工机器。铸造车间调度问题和传统批调度问题比较类似，但在铸造车间每个批次的造型工序和制芯工序可以同时由不同机器加工，有沙箱和批处理机2种加工资源可供选择。因此，铸造车间调度问题可以被描述为具有工序并行性的批调度问题。在该调度问题中，需要解决以下4个子问题：工件分批；批次选择沙箱；确定各工件的加工顺序；确定各工序的加工机器。

1.2 数学模型构建

针对工序并行性的批调度问题，以优化最大完工时间和沙箱空置率为目标，建立调度模型。为方便描述，引入以下假设条件和变量符号进行说明。

假设如下。

1) 单个工件的尺寸未超出所有沙箱中的最大尺寸；

2) 单个工件的重量未超出熔炼炉单次最大熔炼重量；

3) 同类型沙箱个数没有限制；

4) 同一个工件只能被分配到一个工件族且只能被分配一次；

5) 同一个工件族造型或制芯只能被一个批处理机加工，且只能被加工一次；

6) 批次加工过程中不允许被中断，同时也不允许向该批次增加或减少工件。

变量描述如下。

n为待加工工件总个数；J为待加工工件集合J={Ji,1≤i≤n}；Si为工件Ji的尺寸；Wi为工件Ji的重量；l为工件材质种类；ai为工件Ji材质；b为批次的数量；Bmk为造型批次Bk，Bm={Bmk,1≤k≤b}；Bck为制芯批次Bk，Bc={Bck,1≤k≤b}；S Bk为批次Bk的总尺寸；WBk为批次Bk的总重量；Tmij为批处理机Mi选择沙箱Gj的造型加工时间；Tcij为批处理机Mi选择沙箱Gj的制芯加工时间；q为沙箱类型数量；G为沙箱类型集合，G={Gj,1≤Gj≤q}；SG j为第Gj种沙箱的尺寸；m为批处理机的数量；M为批处理机集合，M={Mi,1≤Mi≤m}；Q为熔炼炉单次最大熔炼重量；L为一个足够大的正数；Cmax为最大完工时间；Vmax为沙箱空置率。

本文决策变量如下。

Xik为0-1变量，工件Ji分配到批次Bk为1，否则为0；

Yk j为0-1变量，批次Bk选择沙箱Gj为1，否则为0；

Ymk j为0-1变量，造型批次Bk选择沙箱Gj为1，否则为0；

Yck j为0-1变量，制芯批次Bk选择沙箱Gj为1，否则为0；

Zmik为0-1变量，批处理机Mi造型批次Bk为1，否则为0；

Zcik为0-1变量，批处理机Mi制芯批次Bk为1，否则为0。

本文的数学模型构建如下。

1) 最大完工时间。最大完工时间是指在满足前文约束和假设的条件下，每台批处理机造型和制芯工序结束后所耗费的最大时间。

2) 沙箱空置率。沙箱空置率是指沙箱空置尺寸占沙箱总尺寸的平均值。沙箱空置率越低表明沙箱的有效利用率越高。

式(3)表示一个工件只能被分配到一个批次；式(4)表示一个批次只能选择一种沙箱，并且造型工序和制芯工序选择同一个沙箱；式(5)表示一台批处理机同时只能造型或者制芯一个沙箱；式(6)表示一个批次的总尺寸等于该批次所有工件尺寸之和；式(7)表示一个批次的总重量等于该批次所有工件重量之和；式(8)表示一个批次的尺寸之和不能超出所选沙箱的尺寸；式(9)表示一个批次重量之和不能超出熔炼炉单次最大熔炼重量；式(10)表示只有同种材质工件才可以被分配到一个沙箱。

2 基于IHS-SA算法设计

和声搜索算法(harmony search,HS)是由Geem等[11]提出的一种新型的元启发式算法。近几年诸多学者[12-13]也对其进行改进，并将其应用在车间调度问题上。和声算法本身适用于解决连续性问题，而批调度为离散性问题，因此本文提出一种改进和声算法并将其应用于批调度问题。改进算法由2个阶段组成，全局搜索阶段和局部搜索阶段。全局搜索阶段由和声算法完成，局部搜索阶段由模拟退火(simulated annealing，SA)算法完成。两阶段保证算法有较强全局搜索能力的同时，能大概率跳出局部次优解，最终找到算法的近似全局最优解。

2.1 编码与解码

本文和声算法的每条和声H=[Hj|Hs]基于工件和沙箱产生。其中，Hj中的值表示待加工工件编号，Hs中的值表示对应工件选择的沙箱类型。由于每个工件有2道工序，且可同时加工，故采用单工序编码方式，即单个工序编码代表2道工序。

表2 砂箱类型信息Table 2 Sandbox information

和声解码是将每条和声还原为有意义的解，本文采用批首次匹配[14-15](batch first fit,BFF)规则对每条和声进行解码。首先，从Hj中选择第1个位置的工件建立第1批次；然后，从Hs中选择对应于该位置的沙箱为第1批次的沙箱，取Hj中第2个位置的工件，将该工件在不违背沙箱尺寸约束、熔炼重量约束和工件材质约束的情况下放入前一批次，若违反上述某个约束条件则建立新批次，同时该批次的沙箱类型为Hj中对应位置的沙箱，直到为所有工件分配好批次和沙箱。各批次的顺序即为工件的加工顺序。采用此方法Hs中某些位置的编码会失去作用，但不影响批次的划分。

表1～3给出了5个工件、2种沙箱和2台批处理机的信息。假设一台熔炼炉单次最大熔炼重量为3，对一条和声H=[Hj|Hs]=[2 4 1 3 5|2 1 2 1 1]进行解码。如图2所示，首先选择第1个工件2建立第1个批次，该批次选择的沙箱类型为2；将工件4加入第1批会违反尺寸约束和材质约束，则以工件4建立第2个批次；将工件1加入第2批会违反材质约束，故以工件1建立第3批；将工件3加入第3批不违反任何约束，说明工件1和工件3可以组成第3批，该批次选择工件1对应位置的沙箱2；工件5加入第3批违反材质约束，则以工件5建立第4批。表4给出了5个工件的分批情况，此时工件3位置对应的沙箱编码未发挥作用，通过编解码解决了工件批次划分和沙箱选择的问题。

表4 工件组批信息Table 4 Workpiece batch information

图2 示例解码图Figure 2 The diagram of case decoding

表1 工件信息Table 1 Workpiece information

2.2 机器分配规则

本文提出2种机器分配规则来解决工序分配和机器选择问题：完工时间最短优先规则(the earliestcompletion time first,ECTF)和机器最早可用优先规则(the earliest available machine first,EAMF)。每个批次的工件有造型和制芯2道工序，但在实际生产中一般优先加工造型工序，因此本文在工序分配时优先分配造型工序。

表3 批处理机加工时间信息Table 3 Batch processing time

式(11)中，TECTF为在ECTF分配规则下每个批次造型和制芯工序总完工时间；Tmxj和Tcyj分别为每个批次造型与制芯时间，根据该规则选择满足造型和制芯完工时间最短的机器，如果有多种机器组合都满足完工时间最短的要求，则随机选择其中一个机器组合。

式(12)中，TEAMF为在EAMF分配规则下每个批次造型和制芯工序总完工时间；Tmxj和Tcyj为每个批次造型与制芯时间，在该规则下首先确定造型工序在所有批处理机中完工最早的机器，其次再确定所有批处理机中制芯工序完工最早的机器，若有多个机器完工时间一致，随机选择其中一个，将2道工序的加工顺序和加工机器作为最终分配方案。

2.3 交换算子和插入算子

插入(insert)算子[16]，如图3所示。在[1,n]之间生成2个随机数，如{2,4}，将位置4的工件插入到位置2，中间各位置工件保持原有次序依次后移，沙箱编码部分采取同样策略。交换(swap)算子如图4所示。随机生成2个位置，在父代中找到2个位置对应的工件，将2个位置的工件互换，将互换后的工件集作为子代，沙箱编码部分采取同样策略。

图3 Insert算子示意图Figure 3 The insert operator

图4 Swap算子示意图Figure 4 The swap operator

2.4 算法流程

为便于下文描述，介绍Pareto支配[17]概念：对于多目标优化问题，若 ∀i，有fi(α)≤fi(β)，且 ∃j，使fj(α)≤fj(β)，则称α支配β(α比β更优)，记做α＜β。Pareto支配数量即种群中受某个解支配的数量，不被种群中任何解所支配的解，其Pareto等级为1，所有Pareto等级为1的解组成Pareto非支配解集，所有非支配解集组成Pareto前沿。

本文定义一种初始化策略提高和声记忆库质量。将工件按照材质分类，并将每种材质的工件按重量升序排列，对应的沙箱采用随机方式产生，在初始化阶段为保证和声记忆库的多样性，按照该策略生成的和声个体占20%，其余的和声采用随机方式生成。由于随机选择的沙箱可能不满足尺寸约束，因此，初始化完成后需要对不满足尺寸约束的沙箱进行纠正，重新选择能满足尺寸约束的最小沙箱类型。

2.4.1 全局搜索阶段

和声搜索算法具有全局搜索能力较强的优点，原始和声算法主要根据和声记忆库取值概率HMCR和音调微调概率PAR产生新解，如式(13)所示。在解决离散问题上，HS算法一般采用最大位置法实现从连续空间到离散空间的映射，但在映射过程中容易失去原始解中包含的有效信息，并且原始和声算法每次迭代只产生一条和声，不能较好地执行全局搜索过程。

本文提出的改进和声算法每次迭代都产生与原始和声记忆库数量相等的新和声。在迭代过程中算法新解的产生基于当前最优调度方案Nbest(Pareto等级为1，若有多个随机选择其中一个)和未调度工件序列Narray产生，如式(14)所示。生成[0,1]上随机数R1。若R1＜HMCR，从和声记忆库中随机选择一条和声，并从第z(z为和声的位置编号)列中选择工件编码和对应的沙箱编码放入新解中，若新解中包含此工件，则从最优解Nbest中选择第一个位置的工件和对应沙箱放入新解中；若R1≥ HMCR，则选择未调度工件序列Narray中第y(随机整数)列工件放入新解中，此时沙箱部分采用随机方式产生，直到新解产生完毕。每次放入新解中的工件都要从最优解Nbest和未调度工件序列Narray中划去，避免非法解的产生。

本文对新解的扰动采用自适应音调微调概率PAR，如式(15)所示。式中，It为当前迭代次数，max It为总迭代次数。当新解产生后采用上述2种算子对新解进行扰动，如式(16)所示。前期解的质量较差，因此，用i nsert算子进行较大范围扰动，后期解的质量比较稳定，故用swap算子进行小范围扰动，式中N表示某个待扰动的解，Nnew表示扰动后产生的新解。

2.4.2 局部搜索阶段

模拟退火算法是一种贪心算法，在迭代过程中引入了随机接受概率，能以一定的概率接受比当前解更差的解，增强了算法跳出局部最优的能力。本文采用基于Pareto支配数量的机制接受新解，即比较邻域扰动前的解θ和邻域扰动后的解θ*支配其他解的数量差Δf=DominateNum(θ)-DominateNum (θ*)(其中，DominateNum(θ)表示解θ支配当前种群中其余解的总数量)，并以Metropolis准则接受新解(即Rand＜exp(-Δf/t)时接受新解，否则不接受，其中，t为当前温度)，当算法达到指定温度或新解达到指定未改善次数时，结束内循环。本文基于2种邻域结构产生新解。

1) 沙箱变异(mutation)。对于每条和声H=[Hj|Hs]，检查Hj中第k维和第k+1维的工件材质是否相同。若材质相同，检查Hs中沙箱类型是否相同；若不同，则令k+1维的沙箱与第k维相同。

2) 批次合并(combine)。对工件进行预分批，检查是否存在2个批次的沙箱和材质都相同。若存在相同的2个批次，则将2个批次的工件放到相近的位置，如图5所示；若B2和B52个批次的沙箱和材质都相同，则将B2和B52个批次的工件放到相近位置。

图5 批次合并示意图Figure 5 The diagram of batch consolidation

算法整体流程如图6所示。

图6 IHS-SA算法流程图Figure 6 The flowchart of IHS-SA algorithm

以表4中5工件4批次的划分方案为例，在2种机器分配规则下，得到调度甘特图，如图7所示。图中字母下标索引表示同一批次的造型和制芯工序，后缀表示批次信息。如Bm-1表示第1批次造型工序，Bc-2表示第2批次制芯工序。从实例可以看出，工件加工顺序相同的情况下采用ECTF规则总完工时间较短。

图7 2种不同规则下调度甘特图Figure 7 The Gantt chart of scheduling under two different rules

3 案例仿真与应用

本文以浙江某企业铸造车间实际生产状况为研究对象，针对该车间工件材质种类多，计划制定效率低等情形，选择该车间某周期任务来验证本文模型的有效性。该车间某周期任务为生产40个工件，工件材质有3种，工件信息如表5所示。该铸造车间目前有3种类型沙箱，第1种沙箱尺寸为1 m3，第2种类型沙箱尺寸为3 m3，第3种类型沙箱尺寸为5 m3，每种类型沙箱个数足够多；3台可用批处理机(M1，M2，M3)，每台都可加工造型和制芯工序，熔炼炉单次最大熔炼重量为20 000 kg，需要确定工件的加工顺序和工序的加工机器以保证完工时间最短和沙箱空置率最小。机器加工时间如表6所示。

表5 待加工工件信息Table 5 The information of workpiece to be processed

表6 批处理机加工时间(造型/制芯)Table 6 The processing time of batch processor (molding/coring) h

算法的运行环境为2.7 GHz CPU，8 G内存，64位win 7系统计算机，编程环境为Matlab 2016。为验证本文IHS-SA算法的有效性，选择与文献[18]中NSGA-II算法及改进前后的和声算法进行对比。经多次试验确定IHS-SA算法的和声记忆库 (种群) HMS=80，HMCR=0.9，PARmax=0.7，PARmin=0.2，模拟退火设定初始温度tb=3，终止温度te=1，温度衰减系数K=0.9，最大未改善次数maxFail=5；改进IHS算法去除局部搜索阶段，其余参数保持和上面一致，原始和声算法HS保持和IHS算法参数一致，但根据式(11)采用连续到离散空间映射方式编解码；NSGA-II算法初始种群数量PopSize=80，交叉概率Pcross=0.6，变异概率Pmutation=0.1。算法迭代次数max It=100。在2种机器分配规则下,每种算法独立运行30次，得出最大完工时间和沙箱空置率的最优值、平均值及最优值出现次数。由于算法受2种规则影响不大，故采用ECTF规则评测算法性能，如表7所示。在ECTF规则下本文所提算法的最优解及最优解的数量均优于其他3种算法。

表7 ECTF规则下实例仿真结果Table 7 The simulation results of example under ECTF rules

本文采用收敛性指标γ、多样性指标Δ和支配性指标Ω[18]评测所提算法性能。由于模型的真实最优Pareto前沿(net optimal Pareto front,NOPF)未知，因此采用4种算法的非支配解集进行非支配比较后的解作为NOPF。各指标计算公式如下。

式(17)中，ds为非支配解集中第s个解与NOPF的欧氏距离；ψ为非支配解集个数。γ越小表明算法的收敛性越好。

式(18)中，dσ和dδ分别为算法非支配解集中极端解与边界点之间的欧氏距离；de为第e个解和第e+1个解之间的欧氏距离；d为所有距离的平均值；φ为个体数目。Δ越小表明算法多样性越好。

式(19)中，|C(ξ)|表示集合ξ中的非支配解数量的交集；Pe为第e种算法的非支配解集。Ω越大表明算法的支配性能越好。

如表8所示，通过计算4种算法的评价指标可以得出，本文改进算法的收敛性、多样性和支配性都比其余3种算法好。IHS算法的收敛性、多样性和支配性都比HS算法更好。HS结果较差的原因有2点。1) 连续空间到离散空间映射时失去部分有效信息；2) 原始和声算法迭代过程中产生的新解个数过少。

表8 ECTF规则下4种算法评价指标对比Table 8 The comparison of algorithm evaluation indexes under ECTF rules

根据实验数据绘制Pareto前沿图，如图8所示。通过对比4种算法的Pareto前沿，IHS-SA算法结果较IHS算法更好，原因在于本文加入局部搜索策略后，算法能够跳出局部最优而趋于全局最优解。NSGA-II算法结果和IHS算法相近，说明2种算法均陷入了局部次优解。

图8 ECTF规则下算法Pareto前沿图Figure 8 The Pareto front graph of algorithm under ECTF rules

根据IHS-SA算法，在ECTF 规则下40个工件被分为17个批次，各批次包含工件如表9所示。图9所示为根据其中一个Pareto非支配解绘制的调度甘特图，甘特图中每一行表示在一台机器上加工，每个矩形表示同一个加工批次，矩形长度表示该批次加工时间，各批次的顺序表示工件的加工顺序，字母下标表示同一批次的造型和制芯工序，后缀为批次信息。如Bm-1表示第1批次造型工序，Bc-2表示第2批次制芯工序，属于该批次的工件10、4、9、2、11、3先加工，甘特图中最大完工时间为66 h，沙箱空置率为 20.329 4%。

表9 ECTF规则下工件批次划分表Table 9 The table of workpiece batch division under ECTF rules

图9 ECTF规则下调度甘特图Figure 9 The scheduling Gantt chart under ECTF rules

在不同的机器分配规则下会产生不同的调度甘特图。在EAMF规则下，40个工件被分为17个批次，如表10所示。根据其中一个Pareto非支配解绘制调度甘特图，如图10所示。甘特图中最大完工时间为77 h，沙箱空置率为15.537 3%。

图10 EAMF规则下调度甘特图Figure 10 The scheduling Gantt chart under EAMF rules

表10 EAMF规则下算法工件批次划分表Table 10 The table of workpiece batch division under EAMF rules

4 总结

本文研究了铸造车间工件多材质、组批多约束、工序可并行的批调度问题，建立相应的数学模型，通过改进和声算法结合模拟退火算法对模型进行求解以优化完工时间和沙箱空置率。针对和声搜索算法提出一种单工序编码方式和BFF解码规则解决工件组批和沙箱选择问题，并提出2种不同机器分配规则解决工序分配和机器选择问题。最后将本文所建立模型进行仿真实验，在不同的分配规则下产生不同的调度甘特图，方便企业根据自身的生产状况选择合适调度方案。

本文仅考虑了铸造车间前工段稳定状态下的调度情形，下一步可以对特殊情况下如何进行重调度和批次划分的问题进行研究。