袁帅鹏, 李铁克, 王柏琳, 张文新, 张卓伦, 余娜娜

(1.北京科技大学 经济管理学院,北京 100083; 2.钢铁生产制造执行系统技术教育部工程研究中心,北京 100083)

0 引言

本文研究了一类具有成组可重入特征的单机调度问题,该问题源于钢铁企业宽厚板生产过程中热轧工序的生产管理需求。热轧是宽厚板生产过程中极其关键的一道工序,其工艺流程如图1所示,来自连铸车间的板坯需按照设定的轧制节奏进入轧机进行轧制,随后进入下游工序加工。与其他常规加工过程不同,板坯的轧制过程具有以下特点:①由于轧制工艺的要求,板坯需依次进行粗轧、精轧两个阶段的轧制,且阶段间需要一定的等待时间,以防止轧制过程中发生不可控的形变;②不同板坯两个轧制阶段的加工时间和等待时间均不相同,为了提高轧机利用率,充分利用粗轧和精轧间的等待时间,企业通常会设计调度策略安排两块板坯作为一组进行轧制,即如果相邻板坯满足成组轧制的条件,会在第一块板坯完成粗轧后,安排第二块板坯进入轧机进行粗轧;待第二块板坯粗轧完成且达到第一块板坯规定的等待时间后,第一块板坯重新进入轧机完成精轧操作,随后第二块板坯再次进入轧机执行精轧操作。若相邻板坯不满足成组轧制的条件,则无法进行成组轧制,只能逐一加工。

图1 宽厚轧制工艺流程

若将板坯视为工件,整个轧制过程可视为一类特殊的单机调度问题,且存在以下调度特征:①工件需分阶段重复进入同一台机器加工,加工过程具有可重入性;②不同阶段间存在固定的等待时间,工件在等待期间可暂时释放机器,达到规定等待时间后必须马上进入机器;③对于相邻加工的两个工件,如果满足成组加工的条件,可以安排其成组加工。综上,该问题可抽象为一类具有可重入和成组加工特征(本文将其定义为成组可重入特征)的单机调度问题,通过确定工件间的调度顺序、以及相邻工件是否成组加工两个子问题,以实现系统生产效率的最大化。

可重入调度在半导体制造、钢铁生产等工业领域具有广泛应用,已有不少学者对该问题进行了研究。FRIHAT等[1]从皮革工业中提炼出一类可重入混合流水车间调度问题,以最大完工时间和拖期成本为目标,分别构建了混合整数线性规划模型和约束规划模型;XU等[2]研究了一类可重入置换流水车间调度问题,以最大完工时间为目标建立了数学模型,并提出一种文化基因算法;朱艳艳等[3]研究了具有不确定加工时间的可重入调度问题,以最小化提前惩罚和拖期惩罚为目标构建了混合整数线性规划模型,提出一种变邻域化学反应优化算法;王凯等[4]研究了不确定服务时间的可重入手术系统调度问题,采用三角模糊函数来表征各项服务的持续时长,并以最小化术后恢复平均完成时间为目标,构建了混合整数线性规划模型,提出一种改进的遗传算法;姚远远等[5]针对可重入混合流水车间调度问题,提出一种改进的多目标灰狼优化算法来优化完工时间和拖期时间。以上文献主要针对流水车间展开,针对考虑阶段间等待过程的单机可重入调度问题,尚未发现相关成果。

另一方面,具有成组加工特征的调度问题同样是国内外学者关注的焦点:COSTA等[6]研究了一类流水车间环境下的成组调度问题,构建了问题的调度模型,提出一种基于并行优化的自适应遗传算法;LIN和YING[7]提出一种基于LKH规则和0-1整数规划模型的两阶段元启发式算法,通过数据实验证实了算法的有效性;袁帅鹏等[8]考虑了一类带有阶段间运输的两阶段流水车间成组调度问题,结合问题特征提出一种协同进化迭代贪婪算法。上述文献为进一步研究提供了较好的基础,但已有成果均假定同组内的工件需连续加工,且工件加工过程不允许中断,这与本文的成组可重入调度特征具有较大差别。

综上,尽管可重入调度和成组调度方面均已存在较多的研究成果,但已有成果与本文所研究的具有成组可重入特征的调度问题存在实质性不同。实际上,此类具有成组可重入特征的调度问题不仅存在于钢铁企业的轧制车间,在其他离散制造行业也同样存在,因此对该调度问题展开研究具有重要意义。

1 问题描述与建模

1.1 问题描述

具有成组可重入特征的单机调度问题描述如下:(1)n个工件J={Ji|i=1,2,…,n}需重复进入同一台机器依次执行两个阶段的加工,加工过程具有可重入特征,记两个阶段的加工时间分别为pi1和pi2;(2)每个工件两个加工阶段间均存在等待时间wi,等待期间工件可暂时释放机器,但到达等待时间后必须立刻进入机器完成第二阶段的加工;(3)由于阶段间存在必要的等待时间,为提高设备利用率,如果相邻的两个工件满足成组加工的条件,则允许两者成组加工,若不满足成组加工的条件,只能逐一加工,本文将其定义为非成组加工。图2给出了成组加工和非成组加工的示意图,显然工件进行成组加工能够充分利用设备的空闲时间,进而缩短制造期。对于先后加工的工件Ji和Jj,进行成组加工的充分条件为:wi≥pj1且pi2≤wj。

图2 成组加工和非成组加工调度示意图

根据三元组表示法,该调度问题可表示为1|rcrc,wj,group|Cmax,其中1表示单机环境,rcrc表示可重入特征,wj表示阶段间等待时间限制,group表示成组特征,Cmax表示目标函数。

1.2 数学模型

为方便描述,表1给出了模型中涉及的符号变量及含义。

表1 符号变量含义

数学模型如下:

(14)

式(1)为目标函数,表示最小化Cmax;约束(2)和(3)表示工件的序列关系约束;约束(4)表示虚拟工件必须被安排在调度序列的首位;约束(5)和(6)分别表示两个工件只有在满足成组加工条件、且存在紧邻关系的前提下才有可能执行成组加工;约束(7)表示对于任一工件,成组加工的次数至多为1;约束(8)和(9)分别为各工件第一阶段、第二阶段的完工时间约束;约束(10)表示相邻加工的两个工件未成组加工时的时序关系;约束(11)和(12)限制了相邻加工的两个工件执行成组加工时的时序关系:后续工件第一阶段的开工时间不小于紧前工件第一阶段的完工时间,且第二阶段的开工时间不小于紧前工件第二阶段的完工时间;约束(13)和(14)分别表示0-1决策变量和连续变量的取值范围。

2 问题性质分析

本节将对问题的复杂性进行分析,同时给出了问题最优解存在的一些性质特征,这些性质特征将为后续算法的设计提供有效支撑。为方便下文描述,首先给出以下定义。

定义1对于问题1|rcrc,wj,group|Cmax的任一调度实例,如果存在工件与其他工件都不满足成组加工的条件,则将该工件记为无成组特征工件;反之,如果存在其他工件能与其进行成组加工,则将该工件定义为有成组特征工件。

定理1问题1|rcrc,wj,group|Cmax具有强NP难特性。

定理2对于问题的1|rcrc,wj,group|Cmax任一实例,如果存在无成组特征工件,那么至少存在一个最优调度方案,其中无成组特征工件是连续加工的,且这些工件加工顺序的变化不会影响Cmax值。

定理3对于任一实例,如果两个工件Ji和Jj满足成组加工的条件,则Ji和Jj进行成组加工比非成组加工情况下Cmax值的减少量sij至多为min(wi+pi2,pj1+wj) 。

定理4对于任一实例,至少存在一个最优调度,相邻工件是否成组加工满足贪婪规则,贪婪规则定义为:从给定调度序列的第1个位置开始,如果当前位置加工的工件与紧后工件满足成组加工的条件,可直接将其与紧后工件进行成组加工,无需考虑紧后工件及后续调度序列的情况。

3 改进的分布估计算法

分布估计算法(Estimation of Distribution Algorithm, EDA)利用概率模型的统计学习取代原有交叉与变异操作,能够有效挖掘出优势基因区块特征,进而产生优秀新个体,已被证实在解决调度类组合优化问题方面具有很强的优势[9,10]。本文提出了一种改进的分布估计算法(Improved Estimation of Distribution Algorithm, IEDA),算法实施过程如下。

3.1 编码解码策略

为提升算法的搜索效率,IEDA采用置换序列的方式进行编码,以此确定工件间的加工顺序,记为X;然后通过图3所示的贪婪规则确定相邻工件是否进行成组加工,其中Y为辅助向量,Y(i)=1表示X中处于位置i的工件与后续工件进行成组加工,Y(i)=0则表示两者不进行成组加工。

图3 基于贪婪规则的解码策略

3.2 初始种群构造

为保证种群的质量和多样性,一个有效的做法是结合问题特征设计启发式规则构造若干高质量的人工解,并将其注入初始种群[11]。对于本文问题,我们提出如下启发式规则来构造人工解:首先将待调度工件按照是否存在成组特征分为两个子集N1和N2,对于N1中的工件,计算两两工件间进行成组加工Cmax时值的减少量sij,得到减少量矩阵S,若不能成组加工则记sij=0;然后按照sij值由大到小的顺序依次选择所对应的工件进行加工,同时将被选择的工件从N1中删除,重复上述操作直至N1为空。对于集合N2中的工件,采用随机方式进行排序。

在构造初始种群时,首先利用上述启发式规则产生第一条染色体,然后采用随机方式生成剩余染色体。需要说明的是,由于N2中工件的调度序列不会对Cmax值造成影响,因此在构造初始种群时,N2中的工件序列不会参与后续的优化操作。

3.3 概率模型设置

概率模型是EDA算法的核心。对于调度类的组合优化问题,工件间的前后顺序和工件在某一基因位上出现的频率是决定染色体优劣的两个主要因素。JARBOUI等[12]设计了一种带有工件区块结构特征的概率模型来指导种群进化,即首先按照适应值大小从当前种群中挑选出给定容量Q的优势染色体集,然后按照式(15)生成子代个体:

(15)

其中δij表示当前优势染色体集中工件Ji被安排在位置j上的概率,ηij表示工件Ji在位置j及其之前出现的次数,其值反映了工件前后作业顺序对调度序列的影响;μi′ij表示工件Ji在位置j且Ji′作为其紧前工件上出现的次数,其值反映了工件块结构的相似性;Qj表示在位置j及其之前尚未参与调度的工件集。该模型较好的反映了工件间调度序列以及基因块的相似性,但并不完全适用本文问题,原因如下:①该模型仅考虑了优势染色体中具有相同基因位的优势区块,这会严格限制生成后代种群的多样性,且对于本文问题,优势区块并不一定必须处于相同的基因位置;②采用乘法运算作为ηij和μi′ij两个因素的综合影响率时,容易导致算法在进化过程中丢失有效信息。为了克服以上缺点,对上述概率模型进行了改进,得到式(16)的概率模型:

(16)

其中δij,ηij和Oj与式(15)相同,λi′i表示优势染色体集中Ji作为Ji′的紧后工件出现的次数,Ji′表示位置j-1处被选中的工件,特别地,当j-1时,我们将λi′i的值定义为0,λi′i值反映了优势染色体集中相似区块的重要性。此外,该概率模型采用加法运算作为ηij和λi′i两个因素的综合影响率。

在构造子代个体时,对于每一个位置j,根据式(16)的概率模型从当前个体尚未调度的工件集中选取工件,进而得到完整的个体。每得到一个新的子代个体,将其与当前种群中最差的个体进行对比,如果优于当前种群最差的个体,则替换之,重复上述过程M次,直至生成M个新的个体。

3.4 局部搜索策略

为提升IEDA的局部搜索能力,本节设计了一种基于最大化减少量的重构策略对每条染色体进行深度搜索:对于任一染色体,首先从子集N1对应的调度序列中随机选取并移除1个工件,记为Ji;然后根据定理3计算Ji与调度序列中其他工件进行成组加工时的最大减少量,并选取能与Ji成组加工且能得到最大减少量的工件,记为Jj,按照最大减少量的配置关系将Ji重新插在工件Jj紧前或紧后的位置,即如果最大减少量是由Ji和Jj成组加工得到,则将Ji插入Jj的紧前位置,如果是由Ji和Jj成组加工得到,则插入紧后的位置,进而得到一条新的染色体;最后判断新染色体是否优于当前染色体,如果优于则替换之。

3.5 算法步骤

本文设计IEDA的实施步骤如下:

步骤1基于启发式规则构造容量为p_size的初始种群;

步骤2采用解码策略计算种群内所有个体的适应值,并按适应值由高到低排序;

步骤3从当前种群中选取前Q条高质量的染色体,形成优势染色体集;

步骤4更新概率模型,并基于概率模型生成M个新个体以更新当前种群;

步骤5使用局部搜索策略对当前种群内的个体进行深度搜索;

步骤6统计算法的运行时间,判断是否达到最大运行时间限制,若未达到,转至步骤2;否则,输出当前最优解及其目标值。

4 问题最优解下界分析

本节通过理论分析给出了问题最优解的两个下界。首先,由于机器在同一时刻能且仅能加工1个工件,因此对于任一问题实例,其Cmax至少大于所有工件需要占用机器的时间,因此我们可得到问题的第一个下界:

(17)

其中,N1表示具有成组特征的工件集,N2表示无成组特征的工件集。

另外,任一调度方案的Cmax值至少大于所有工件Cmax值的上限与最大减少量的差值,因此可得到问题最优解的第二个下界:

(18)

5 仿真实验

5.1 实验设计

检验IEDA算法的有效性,本节将 IEDA与其他三种元启发式算法进行对比。实验采用离散均匀分布的方式随机生成实验数据,其中工件规模n={40,60,80,100,150};各工件两个阶段的加工时间和阶段间等待时间均为介于a和b之间的离散均匀分布,其中a的取值范围为[5,50β1],β1∈{0.5,1};b的取值范围为[a+1,⎣a(1+β2)」],β2∈{0.5,1}。对于每种工件规模n、时间参数β1和β2的组合,随机生成10组测试算例,共有200组测试算例。采用相对偏差率(Percentage Relative Difference, PRD)作为算法性能评价指标,PRD计算公式为:

PRD=(Cmax(A)-ref)/ref

(19)

其中,Cmax(A)表示当前算例算法A得到的Cmax值,ref为当前测试算例的参考值,计算方式为两个下界较大的值,即ref=max(LB1,LB2)。

5.2 算法参数设置

IEDA中的关键参数包括种群规模p_size,优势群体集容量Q,每次迭代新个体的生成数量M,本文使用田口实验设计方法探讨了这些参数对算法性能的影响,结果表明,当p_size=30,Q=10,M=10时,IEDA具有最好的求解性能,因此采用该参数设置执行后续实验。

5.3 实验结果

由于没有可直接应用于本文问题的相关求解算法,我们选取了三种针对单机调度问题提出元启发式算法,分别为JARBOUI等[12]中的EDA算法(记为EDA_1)、SALAMA和SRINIVAS[13]中的嵌套变邻域搜索的模拟退火算法(Simulation Annealing, SA)和WANG等[14]中的离散人工蜂群算法(Discrete Artificial Bee Colony,DABC),三种算法均采用原文推荐的参数设置。针对200组测试算例,将四种对比算法的最大运行时间均设置为n秒,计算每组实验下10组算例的PRD均值和最大值,结果如表2所示。

表2 四种算法PRD均值(Mean)、最大值(Max)统计结果

由表2可以看出,对于20组测试实验,IEDA算法的PRD均值全部低于其他三种对比算法,且其总体PRD均值仅为3.16%,明显低于EDA_1、SA和DABC,这说明针对本文问题,IEDA算法的求解质量要优于其他三种算法。此外,20组实验IEDA算法求得的PRD最大值也全部低于其他三种算法,表明IEDA在获得较优解的同时,求解性能也相对稳定。与EDA_1相比,所有测试算例下IEDA算法的总体PRD均值降低了15.28%,PRD最大值降低了17.58%,这证实了IEDA中对概率模型的改进是有效的。此外,20组测实验四种算法求得的PRD最大值仅为 9.19%,非常接近于本文提出的下界值LBmax,这说明本文给出的两个问题最优解下界能够较好的评估算法的性能。

为进一步检验四种算法在收敛速度方面的差异性,以n=100的样本数据为例,随机选取一次实验结果,绘制了四种算法的收敛曲线,结果如图4所示。由该收敛曲线可以直观看出,IEDA具有更好的收敛性。此外,为检验四种算法在不同时间参数配置下求解性能的差异性,根据实验结果汇总了β1和β2取不同值时四种算法PRD均值的差异性,并对其进行95%的置信区间分析,结果如图5所示。由图5可以看出,对于不同的时间参数配置,IEDA算法均能获得较低的PRD均值,这进一步证实了IEDA的有效性。

图4 四种算法收敛曲线对比

图5 不同时间参数PRD均值95%置信区间

综上,可得出结论,针对本文问题,IEDA算法的求解质量要优于其他三种对比算法。

6 结语

本文研究了一类具有成组可重入特征的单机调度问题,该问题具有较为广泛的工业应用背景,但尚未发现相关研究成果。针对此调度问题,首先构建了混合整数线性规划模型,然后证明了问题的复杂性和最优解存在的性质特征,随后提出一种改进的分布估计算法,并给出了问题最优解的两个下界。最后通过数据实验证实了所提算法的有效性。研究成果能够为钢铁生产等制造企业解决相应调度问题提供有效的理论支撑。本文只研究了两个工件进行成组加工的情形,实际工业生产中存在三个甚至多个工件进行成组加工的情况,下一步将进一步研究多个工件成组加工的调度问题,并对问题的性质特征进行分析,进而设计更为高效的求解算法。