Processing math: 68%

上海交通大学学报, 2024, 58(8): 1271-1281 doi: 10.16183/j.cnki.jsjtu.2023.070

机械与动力工程

基于Benders分解和分枝定界的随机交期批量流流水车间调度

石亚东, 刘冉,, 王铖恺, 吴泽锐

上海交通大学 机械与动力工程学院,上海 200240

Stochastic Due-Date Lot-Streaming Flowshop Scheduling with Benders Decomposition and Branch-and-Bound

SHI Yadong, LIU Ran,, WANG Chengkai, WU Zerui

School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

通讯作者: 刘 冉,副教授,博士生导师;E-mail:Liuran2009@sjtu.edu.cn.

责任编辑: 李博文

收稿日期: 2023-03-2   修回日期: 2023-04-21   接受日期: 2023-05-12  

基金资助: 上海市科委“科技创新行动计划”高新技术领域项目(22511103603)

Received: 2023-03-2   Revised: 2023-04-21   Accepted: 2023-05-12  

作者简介 About authors

石亚东(1999-),硕士生,主要研究方向为运筹优化算法设计.

摘要

针对交期随机的批量流车间调度问题,以最小化工件延期期望之和为目标,推导出工件交期符合3类经典随机分布条件下问题目标的闭式计算表达式.建立考虑换模时间与随机交期的问题数学模型,针对模型高度非线性特征对其线性化.设计一种基于逻辑的Benders分解(LBBD)与分枝定界相结合的优化算法,提出两种有效加速策略提升算法求解效率.数值实验结果验证了算法的有效性,通过随机交期与确定交期结果的比较,验证了考虑随机的必要性.

关键词: 随机交期; 批量流; Benders分解; 分枝定界

Abstract

The lot-streaming flowshop scheduling problem with stochastic due time is addressed in this paper, with the objective of minimizing the sum of expected job delays. Closed-form expressions for the expected delays of jobs are derived under three classical distribution conditions. A mathematical model is then formulated, considering set-up times and stochastic due time. To address the highly nonlinear nature of the model, a linearization is performed. Furthermore, an optimization algorithm is designed using a Logic-based Benders decomposition (LBBD) approach combined with branch-and-bound. Two effective acceleration strategies are introduced to improve the efficiency of the algorithm. The numerical experiments demonstrate the effectiveness of the proposed algorithm, and the necessity of considering stochastic lead times is verified by comparing the results with those obtained from deterministic due time.

Keywords: stochastic due time; lot-streaming; Benders decomposition; branch-and-bound

PDF (1434KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

石亚东, 刘冉, 王铖恺, 吴泽锐. 基于Benders分解和分枝定界的随机交期批量流流水车间调度[J]. 上海交通大学学报, 2024, 58(8): 1271-1281 doi:10.16183/j.cnki.jsjtu.2023.070

SHI Yadong, LIU Ran, WANG Chengkai, WU Zerui. Stochastic Due-Date Lot-Streaming Flowshop Scheduling with Benders Decomposition and Branch-and-Bound[J]. Journal of Shanghai Jiaotong University, 2024, 58(8): 1271-1281 doi:10.16183/j.cnki.jsjtu.2023.070

置换流水车间调度问题(Permutation Flowshop Scheduling Problem, PFSP)是经典的调度问题.当加工机器数大于2时,PFSP是NP-hard问题[1],求解困难.经典PFSP假设每种工件不可分割,即一种工件在当前机器全部加工完成后,才能转移到下一机器加工,但此假设常会造成机器利用率的损失.在很多现实生产场景中,一种工件可被分为若干个子批,每个子批可以被单独加工,即工件的一个子批完成某个机器上的加工后,就可被转移到下个机器.如该子批分割技术使得同种工件的多个子批可同时在多台设备上加工,则被称为批量流技术.批量流技术运用于PFSP就形成了批量流车间调度问题(Lot Streaming and Scheduling Problem,LSSP),是经典PFSP的拓展.相较于经典PFSP,LSSP更为复杂,不仅需要确定工件的加工顺序,还要决策每个工件的子批分割策略.

目前PFSP受到诸多研究关注,其求解方法包括精确算法[2-3]和启发式算法[4-7],而对LSSP的研究则不够深入.首先,对LSSP的研究主要集中在启发式算法[8-11],部分学者研究了最小化延迟目标下的批量流调度问题[12].相较于启发式算法,求解LSSP的精确算法非常有限,且只能解决一些特殊问题.Alfieri等[13]针对同种产品的两机批量流调度问题提出动态规划算法,实现精确求解.除了求解算法的局限性,目前对LSSP的研究均基于确定性假设,即认为所有的信息是确定的,在此基础上进行分批和调度的优化.但在实际场景中,诸多制造因素存在不确定性,如工件的加工时间、交期等受到很多因素影响而不能简单设定为确定已知信息.虽然文献[14]中研究了随机环境下的PFSP,但其所提方法仅考虑工件的调度问题,不涉及工件批量的决策,因此难以运用于随机LSSP.

针对现有研究不足,引入并研究一类新的LSSP,即不确定交期的批量流流水车间调度问题.对比现有文献,本研究有以下特点.首先,考虑每种工件具有交期限制,即下游采购商或下游车间对每种工件提出“交货时间”的要求,对超出交期完成的工件按照其延期的时间计算延迟,并将全部工件的总延迟作为优化目标.其次,在实际场景中,每种工件的交期具有不确定性,这造成了问题具有随机优化性质,增加了问题求解难度.在求解方法方面,推导了不同交期分布下工件延期期望的闭式表达式,以最小化工件延期期望之和为目标建立数学模型,设计Logic-based Benders分解与分枝定界相结合的精确优化算法,实现了复杂不确定环境中批量流调度问题的求解.最后,针对大规模的问题数据,将禁忌搜索算法嵌入Benders分解框架以实现高效启发式求解.研究将LSSP在运用场景和求解方法方面都进行了拓展和创新.

1 问题描述与数学模型

1.1 问题描述

考虑由K台机器组成生产线,其中机器位置固定且缓冲区容量无限,N种工件需要在K台机器组成的生产线上加工.对于每种工件n∈{1, 2, …, N},其总加工任务量为Qn,Qn可以被分为若干个子批,子批批量规定为最小批量Smin的整数倍,同种工件的子批批量不一定相同.工件n最多可被分为un个子批,un=Qn/min S.每种工件的每个子批批量、子批之间的加工顺序在所有机器上均保持一致,即满足LSSP中的“子批一致性”假设.所有子批按照加工顺序依次经过K台机器,只有在前一个子批全部加工完成并离开机器后,下一子批才可进入该机器.对于每一台机器,同一时间仅可加工一个工件,所有机器在系统时刻为0时处于空闲状态,且不考虑机器故障.假设工件的来料时间都为0,工件i在机器k上的加工时间确定,记为pik.

每台机器上有加工切换的换模时间,即机器k完成工件i后继续加工工件j,其换模时间为sijk.机器在换模期间不能加工工件,且换模时间与前后工件种类i,j相关.在一个工件的子批到达机器后才能启动机器换模,且不考虑工件在机器间的运输时间.s0jk表示当第一个子批属于工件j时,在第k台机器上开始加工前的换模时间.

每个工件n都存在不确定的交期dn,假设交期时间具有随机性且工件交期互相独立.如果一个工件的完工时间ten超过dn,则延迟为(ten-dn)+.记工件n的延期为CPn,所有工件的延期之和为

CT=Nn=1CPn=Nn=1(ten-dn)+

以上问题需进行3类决策:①每种工件应分为多少个子批;②每个子批包含多少批量;③子批的加工顺序,以最小化生产线在随机交期下的延期期望.对于确定的交期,批量流调度对于减少工件延期的意义如图1所示.J11代表工件1的第1个子批,J12代表工件1的第2个子批,以此类推.M1, M2, M3分别代表机器1、机器2、机器3,子批需依次进入每一台机器加工.图中,(a)对应不分批策略,(b)与(c)分别对应分批调度策略1和策略2.对比(a)、(b)和(c)可知,对工件的分批操作有利于减少工件的延期.对比(b)与(c)可知, 不同的分批调度策略对于减少延期的效果不同,且策略1优于策略2.由此可知,批量流技术有利于减少工件延期,且不同分批策略效果不同.

图1

图1   批量流调度甘特图

Fig.1   Gantt chart of lot-streaming scheduling problem


1.2 数学模型

由于CT是随机量,所以模型的优化目标是最小化所有工件的延期期望E(CT).基于以上问题定义,建立随机混合整数规划(Mixed Integer Programming,MIP)模型,模型中参数符号与决策变量的含义如表1表2所示.

表1   参数符号及含义

Tab.1  Symbol and meaning of parameter

参数符号含义
k机器的索引号,k∈{1, 2,…, K}
j, h工件i, g的子批索引号,j∈{1, 2,…, ui}, h∈{1, 2,…, ug}
i, g, n工件的索引号,i, g, n∈{1, 2,…, N}

新窗口打开| 下载CSV


表2   决策变量及含义

Tab.2  Decision variables and meanings

决策变量含义
aij辅助变量,用于顺序记录
rij整数变量,工件i的第j个子批的最小批量个数,可以为0,0≤iN, 1≤jui
tijk连续变量,工件i的第j个子批在机器k上的完工时间, 0≤iN, 1≤jui, 1≤kK
yij0-1变量,rij>0时,yij=1,rij=0时,yij=0. 0≤iN, 1≤jui
zijgh0-1变量,当工件i的第j个子批和工件g的第h个子批均不为0,且工件i的第j个子批为工件g的第h个子批的紧前子批时,zijgh=1;否则,zijgh=0. 0≤iN, 1≤jui, 0≤gN, 1≤hug
Sij整数变量,工件i的第j个子批的工件数量,0≤iN, 1≤jui

新窗口打开| 下载CSV


目标函数为

min E(CT)= Nn=1E(CPn)
(1)

约束为

tijktghk+Sijpik+sgik-(1-zghij)B0iN, 1jui, 0gN 1hug, 1kK, (i, j)(g, h)}
(2)
tijktij(k-1)+Sijpik+sgik-(1-zghij)B0iN, 1jui, 0gN1hug, 1kK, (i, j)(g, h)}
(3)
yijrijyijB, 0≤iN, 1≤jui
(4)
Sij=rijSmin, 0≤iN, 1≤jui
(5)
uij=1Sij=Qi, 0≤i≤N
(6)
tei≥tijK, 1≤i≤N, 1≤j≤ui
(7)
tij0=0, 1≤iN, 1≤jui
(8)
Ii=0uij=1zijgh=ygh, 1≤g≤N, 1≤h≤ug
(9)
Ig=1ugh=1zijgh≤yij, 0≤i≤N, 1≤j≤ui
(10)
Ig=1ugh=1z01gh=1
(11)
aij-agh+BzghijB-1, 0iN 1jui, 0gN, 1jug}
(12)
y01=1
(13)
yij, zijgh{0, 1}, tijk0int Sij0, int rij0, uij00iN, 1jui, 0gN 1hug, 1kK}
(14)

目标函数式(1)表示最小化所有工件的延期E(CT).约束式(2)表示同一机器上前后子批的加工约束.约束式(3)表示同一子批在相邻两台机器上的加工约束.约束式(4)表示yijrij之间的关系.约束式(5)表示每个子批包含的工件数量为最小批量的整数倍.约束式(6)表示同种工件全部子批包含的工件数量之和等于该工件的加工任务量.约束式(7)表示每种工件的完工时间,由该工件的最后一个子批在最后一台机器上的完工时间决定.约束式(8)表示所有工件在零时刻都可加工.约束式(9)与(10)分别约束了每一个子批的紧前与紧后子批最多只有一批.约束式(11)表示虚拟工件0为第一个子批.约束式(12)保证模型得到的加工顺序无闭环.约束式(13)表示虚拟工件是非空子批.约束式(14)表示各变量的取值范围.B是一个很大的正数,只有zijgh=1时约束式(2)、(3)才生效.

2 延期期望的计算与线性化

模型的目标是最小化延期E(CT),由于难以直接表达,本章研究当交期服从3种经典分布时 E(CT) 的闭式表达式及线性化方法.首先引入定理1.

定理1 对于随机交期的LSSP问题, fnFn表示dn的概率密度函数和分布函数,符号An代指limx-xFn(x),当An存在,延期E(CT)计算公式为

E(CT)=Nn=1An+Nn=1ten- Fn(x)dx

证明:

CPn=(ten-dn)+E(CPn)=ten- (ten-x)fn(x)dx=ten- (ten-x)dFn(x)=
[Fn(x)(ten-x)]|ten-+ten- Fn(x)dx=-limx-(Fn(x)ten-Fn(x)x)+ten- Fn(x)dx=
An+ten- Fn(x)dx

由于工件交期相互独立,所以有

E(CT)=Nn=1An+Nn=1ten- Fn(x)dx

进一步对于均匀分布、指数分布、正态分布,An=limx-xFn(x)=0,证明见附录A,因此Nn=1An可从E(CT)表达式中消除.

2.1 3类经典分布下延期期望的计算

基于定理1,定理2~4推导了3类经典分布下 E(CT) 的闭式表达式.

定理2 当交期dn相互独立,且服从在区间[an, bn]的均匀分布,dn~U(an, bn)时,E(CT)的计算表达式为

E(CT)=Nn=1E(CPn)
E(CPn)={0,tenan(ten-an)22(bn-an),antenbnten-an+bn2,tenbn

证明 由于dn~U(an, bn),则

Fn(x)={0,xanx-anbn-an,anxbn1,xbn
E(CPn)=ten- Fn(x)dx=
{0,tenan(ten-an)22(bn-an),antenbnten-an+bn2,tenbn

定理3 当交期dn服从指数分布,dn-an~E(λn),并相互独立.E(CT)的表达式为

E(CT)=Nn=1E(CPn)

其中:

E(CPn)=ten- Fn(x)dx={0,tenanten-an+1λn[exp(-λn(ten-an))-1],tenan

证明 由于dn-an~E(λn),则

Fn(x)={0,xan1-exp(-λn(x-an)),xan
E(CPn)=ten- Fn(x)dx={0,tenanten-an+1λn[exp(-λn(ten-an))-1],tenan

定理4 当交期dn服从正态分布,即dn~N(μn, σn),并相互独立.E(CT)的表达式为

E(CT)=Nn=1E(CPn)

其中:

E(CPn)=(ten-μn)Φ(ten-μnσn)+σn2πexp[-(ten-μn)22σn]

Φ(ten-μnσn)为标准正态分布的分布函数.

证明

E(CPn)=ten- (ten-x)fn(x)dx=ten- (ten-x)1σn2πexp[-(x-μn)22σn]dx

y=x-μn,则

E(CPn)=ten-μn- (ten-μn-y1σn2πexp[-y22σn]dy=(ten-μn)ten-μn- 12πσnexp[-y22σn]dy+
σn2π+(ten-μn)2 exp[-y22σn]dy22σ2n=(ten-μn)Φ(ten-μnσn)+σn2πexp[-(ten-μn)22σn]

2.2 模型线性化

在2.1节的3种分布下,目标函数都是非线性函数,模型求解困难,考虑对目标函数线性化处理,先引入定理5.

定理5 延期E(CPn)是单调递增的凸函数.

证明 见附录A.

基于定理5,E(CPn)可以基于分离规划思想进行线性化[15],以均匀分布为例,若对区间[an, bn]线形化,则目标函数中的E(CPn)可以由G条割线近似.设第t条割线的斜率为kt,截距为bt,E(CPn)表达式可写为

E(CPn)={0,tenanktten+bt,xttenxt+1(t=1, 2, , G)ten-a+b2,tenbn

E(CPn)是单调递增的凸函数,可直接在原模型的基础上增加以下约束:

E(CPn)≥atten+bt (t=1, 2, …, G)
(15)
E(CPn)≥ ten- an+bn2
(16)

对于其他两种经典分布,可以得到类似的线性化数学模型.

3 求解算法

流水车间批量流调度是典型NP-hard问题,本研究模型还考虑到顺序相关的换模时间与随机交期,进一步增加了模型求解难度,小规模算例也包含了大量的约束与变量,使用Gurobi等求解器难以在有限时间内完成求解.Benders分解是求解大规模优化问题的经典算法,通过分离容易处理的变量和复杂变量,将问题分解为主问题和子问题,对变量分别求解.经典Benders分解基于线性对偶,要求子问题必须是线性规划问题,限制了应用范围.LBBD算法是其拓展,使用逻辑意义上的割替代子问题线性对偶产生的割,不再要求子问题是线性规划,提升了算法的适用范围.本模型包括两类决策变量:子批分割决策变量和加工顺序决策变量.Benders分解算法将两类决策变量分配到主问题和子问题中分别求解,不需要同时考虑所有决策变量和约束,可大幅降低问题复杂度,提高求解效率.单独考虑加工顺序时,子问题可转化为经典的流水车间调度问题,提出一种分枝定界算法对其求解.考虑到两类决策变量均为整数变量或0-1变量,故设计LBBD算法求解模型.

3.1 批量流问题Benders分解

首先,去除原模型中与加工顺序有关的约束,只保留子批分割方案相关的约束,形成Benders分解的主问题MP,模型如下.

目标函数

min E(CT)=Nn=1E(CPn)

约束如下:

式(4)~(8), 式(16)

tijktij(k-1)+Sijpik+ming sgik0iN, 1jui,1kK}
(17)
teitij0+uij=1Sijpik1iN, 1kK}
(18)

式中:mingsgik为机器k上切换到工件i所需的最小换模时间.约束(17)为同一子批在相邻机器的完工时间约束.约束(18)为任意一种工件的完工时间不小于该工件的所有子批在任意一台机器上的加工时间之和.这两个约束是子问题的松弛,目的是加速算法的收敛.求解主问题得到批量分配方案,即一组rij,Sijyij的值.子问题SP基于这组参数,对子批的加工顺序进行决策,数学模型如下:

目标函数:

min E(CT)=Nn=1E(CPn)

约束如下:

式(16)~(17)

tijktghk+-Sijpik+sgik-(1-zghij)B0iN, 1jui, 0gN1hug, 1kK, (i, j)(g, h)}
(19)
tijktij(k-1)+-Sijpik+sgik-(1-zghij)B0iN, 1jui, 0gN1hug, 1kK, (i, j)(g, h)}
(20)
Ii=0uij=1zijgh= -ygh, 1≤g≤N, 1≤h≤ug
(21)
Ig=1ugh=1zijgh-yij, 0≤i≤N, 1≤j≤ui
(22)
Ig=1ugh=1z01gh=1
(23)
uij-ugh+BzghijB-10iN, 1jui0gN, 1jug}
(24)
y01=1
(25)
zijgh{0, 1}, tijk0, int uij00iN, 1jui, 0gN1hug, 0kK}
(26)

式中:I为工件集合.约束中的-Sij-ygh为主问题给出的最优解中Sijygh的取值,约束含义与1.2节中相同.针对一组给定的批量分配方案,SP是一个排序问题,一定存在可行解,因此算法中只包括最优割.提出最优割之前,在模型中增加一个指示变量rijw,表示工件i的第j个子批批量是否为w,是则取1,否则取0.相应地,为表示指示变量rijwrij的关系,向模型主子问题中添加如下约束:

uiw=0rijw=1, 0≤i≤N, 1≤j≤ui
(27)
rij= uiw=0wrijw, 0≤i≤N, 1≤j≤ui
(28)

提出Benders割:

CTCSPT {1-[ Ii, j, w|rMPijw=1(1-rijw)]}
(29)

式中:CSPT为子问题给出的最优值;rMPijw为本轮迭代中rijw的取值.通过该割,可确保下一轮迭代中,若MP求解所得各变量取值与本轮相同,目标值CT不小于CSPT.

3.2 Benders子问题的分枝定界算法

以上Benders分解的子问题是典型的排序问题,也具有NP-hard性质,提出求解子问题的分枝定界高效算法.

3.2.1 分枝定界下界设计

模型中,由于交期是随机变量,首先引入随机序的定义, 并提出假设1.

定义1 文献[16]中,一个随机变量X在随机序上小于等于另一个随机变量Y,是指对应的分布函数满足FX(t)≥FY(t),∀t∈R,记为XstY.

假设1 任意两个工件的交期di,需满足随机序上的大小关系,即∀i, g∈1, distddistd.

在分枝定界搜索树中,根节点代表空的排程,其余节点代表部分的子批排程序列σ=(σ(1), …, σ(s)), 1≤sl.其中,元素σ(s)代表每一台机器上第s个位置被加工的子批,s为已排程的子批数量,l为子批的总数量.σ(1)=(2, 3),表示每台机器上第1个被加工的是第2种工件的第3个子批.已排程的子批序列σ,与未排程的l-s个子批构成的任一序列-σ,共同构成一个完整的加工序列σ-σ.由于所有机器上的子批加工顺序相同,故后文用“第t个位置的子批”代指完整序列σ-σ中的第t个元素.将未排程的任意一个子批(i, j)放在已排程序列σ的后一位,构成下一层节点σ(i, j).

工件的延期取决于完工时间,而工件的完工时间取决于该种工件最后一个子批的完工时间.因此在某个节点,工件可以被分为两类:一类是存在未排程子批的工件,这类工件的集合记为J,计算此类工件延期的下界gJLB,需得到其最后一个子批在最后一台机器上完工时间的下界;另一类是所有子批均已排程的工件,这类工件的完工时间是确定的,延期容易求得,记为gN\JLB.考虑上述两类工件,则每个节点延期下界gLB=gJLB+gN\JLB.

推导第t个位置(ts+1)的子批在最后一台机器上的完工时间下界RtK,基于Chung等[2]的研究,发现Rtk递推公式如下.

步骤1-1 对于第1台机器,令Es+1, 1=L1(σ),Rs, 1=L1(σ).其中,Etk是第t个位置的子批在第k个机器上的开工时间下界,Lk(σ)是已排程的子批序列σ在机器k上的完工时间.

p[h]k表示{pijk|(i, j)∈U}从小到大排列的第h个值,U为未排程的子批集合,则第t个位置子批在第一台机器上的完工时间下界Rt1=Rs, 1+t-sh=1p[h]1, t=s+1, …, l,开工时间下界Et1=Rt-1, 1, t=s+2, …, l.

步骤1-2 对于第2台机器到第K台机器,第s+1个位置的子批的开工时间取决于本台机器已排程序列的完工时间或之前机器上已排程序列的完工时间Es+1, k=max{Lk(σ), maxu, k(Es+1, u+q[1], u, k-1)},其中,q[1], u, k-1是{qi, j, u, k-1|(i, j)∈U}中最小值,qijuv=vk=upijk, 0≤uvK,即子批 (i, j) 从机器u加工到机器v所用时间.得到Es+1, k后,有Rsk=Es+1, k.对于第t个位置子批在第k台机器上的完工时间下界,与第1台机器类似,Rtk=Rsk+t-sh=1p[h]k, t=s+1, …, l.

基于上述递推公式,可得第t个位置的子批在第K台机器的完工时间下界RtK,计算集合J中工件的延期下界.例如,若工件i的最后一个子批排在第t个位置, 则RtK是其在最后一台机器上完工时间的下界,因此E(RtK-di)+是工件i延期的下界.

基于以上递推公式,介绍下界gLB的具体计算步骤.

步骤2-1 初始化.首先对每个子批的换模时间和加工时间进行合并.对每一个pijk,将切换到工件i所需的最小换模时间并入加工时间,即pijk=Sijpik+mingsgik.在每台机器上,将未排程的子批加工时间升序排序得到p[1]k, p[2]k, …, p[l-s]k, 1≤kK.

步骤2-2 根据已排程的子批序列σ,更新φ(i),Lk(σ),gN\JLB.其中,φ(i)为集合J中工件i未排程的子批数量.由于σ中子批顺序已知,所以Lk(σ)和gN\JLB容易计算.由步骤1中递推公式计算RtK, s+1≤tl.

步骤2-3 对{φ(i)|iJ}升序排列为φ[1],φ[2], …, φ[|J|], 令ϕh=hi=1φ[i], 1≤h≤|J|.

步骤2-4 对于J中工件,其延期下界为gJLB=|J|h=1E(Rϕh, K-d[h])+, 其中,d[h]是在随机序标准下{di|iJ}从小到大排列的第h个随机变量,此延期下界的证明见附录.最终得到所有工件延期E(CT)的下界gLB=gJLB+gN\JLB.

总结上述步骤,即对于所有未排程子批,同种工件的子批连续排列,工件之间按照交期随机序升序排列,获得节点延期期望的下界.

3.2.2 初始解生成

在分枝定界算法中,通过禁忌搜索获取初始解,作为问题的上界帮助剪枝.首先,禁忌搜索的初始解s0由NEH算法[17]产生.获得解s0后,构造当前解的邻域集合,邻域结构包括子批交换和子批插入.子批交换是指任意两个位置的子批交换位置;子批插入是将某个位置的子批取出,重新插入到另一个子批之前.获得邻域集合后,将本次迭代的邻域最优解作为下一次迭代的搜索起点.为避免算法陷入局部最优,设计禁忌规则如下:若本轮最优邻域解由交换或插入操作获得,则将执行交换操作的子批对或插入操作的子批加入禁忌表中,并在之后的θ轮迭代中禁止对禁忌表中的子批对或子批操作.

对于解的评估,将工件完工时间带入线性化之后的目标函数中,即可计算对应排程的延期.若某邻域解由被禁忌的操作产生,但其目标函数值优于当前最优解,实行赦免操作,即取消对该解的禁忌.本节算法的终止条件为一定的迭代次数,算法运行至设定迭代次数后输出当前最优解.

3.2.3 分枝定界算法流程

整个分枝定界算法流程如下.

步骤3-1 初始化.设置初始层高r=0, 已排程子批序列σr=⌀,未排程子批集合Ur包含全部子批.对于所有(i, j)∈Ur,令节点下界gr(ij)=0.由前文禁忌搜索获得上界Z*.

步骤3-2 下界计算.对于{(i, j)∈Ur|gr(ij)<+∞},计算将子批(i, j)加入序列σr后得到节点下界gr(ij).如果 min(i, j)Urgr(ij)≥Z*,进入步骤3-4,否则进入步骤3-3.

步骤3-3 分枝.如果层高r<l-1,令(i0, j0)=argmin(i, j)Urgr(ij),r=r+1,σr=σr-1(i0, j0),Ur=Ur-1-(i0, j0),gr-1(i0j0)=+∞,进入步骤3-2;如果层高rl-1,此时Ur只包含元素(i0, j0), 令σ*=σr(i0, j0),Z*=gr(i0j0),进入步骤3-4.

步骤3-4 搜索.①令r=r-1,如果r<0则输出最优解σ*与对应目标值Z*;②如果存在 min(i, j)Urgr(ij)≥Z*,则重复步骤3-4,否则进入步骤3-3.

3.3 Benders算法加速策略

3.3.1 两阶段加速策略

在Benders算法中,如果主问题能够通过子问题的松弛来加强,那么收敛的速度会显著加快.在模型中,子问题的松弛是指在MP的解确定后工件延期的下界.显然,对MP确定的部分子批排程得到的是子问题的一个下界.在第一阶段,确定交期均值最小的几种工件为紧要工件,在子问题中只对紧要工件的子批排程,将生成的最优割加入主问题中,直至收敛.第二阶段,在子问题中对所有子批排程,继续进行Benders迭代,直至收敛.步骤如图2所示.

图2

图2   两阶段加速策略流程图

Fig.2   Flow chart of two-stage acceleration strategy


3.3.2 重复枝去除策略

在上文提到的分枝定界算法中,每个节点代表部分的排程,分枝过程就是将未排程的子批分别作为已排程序列的下一个子批.然而在未排程的子批中,可能存在一些同种工件相同批量的子批,这样本质相同的子批对应着完全相同的子树.在分批过程中,对于重复枝的剪枝,就是在分枝过程中检查是否存在本质相同的子批被分成多枝.重复枝的去除可以减少分枝定界过程中的冗余,减少求解时间与问题规模,提高计算效率.

4 数值实验

依据文献[18]中的标准算例,对其工件交期及同种工件之间换模时间修订,产生了数值实验的数据1(1数据已上传至 https://www.aliyundrive.com/s/k5FEzaMnJks).算例用N-J-K符号表示,即在K台机器上加工N种工件,每种工件最多可被分为J个子批.算例中工件交期满足不同的分布类型,用不同的字母表示,U为均匀分布, E为指数分布, N为正态分布.实验均在3.1 GHz CPU和32 GB内存机器运行.

4.1 随机性分析

为验证随机性对实验结果的影响,调整工件交期的标准差,观察结果的变化.由于指数分布中标准差与均值是相等的关系,不能调整工件交期的标准差.以均匀分布和正态分布为例展示随机性的影响.8个算例只有交期标准差与分布类型不同,其他信息相同,计算结果如表3所示.由表可知,当工件交期服从同种分布时,随着交期方差的变化,目标函数值、分批数均有波动,子批的加工顺序也有明显变化.当工件交期服从期望标准差相同的不同分布时,目标函数值、工件分批与排序方案也存在较大差异.当标准差为0,即交期是确定值时,求解结果为235.8,分批数为14,和考虑随机性质下的问题解有明显差别.可见考虑随机性的影响十分必要.

表3   交货时间随机性影响分析

Tab.3  Influence of stochastic due time

算例分布类型交期标准差结果分批数
1N20239.814
2N50252.913
3N100282.613
4N150334.612
5U20241.414
6U50254.215
7U100282.410
8U150331.111

新窗口打开| 下载CSV


4.2 精确算法求解实验

为验证设计的LBBD算法效果,采用15个算例进行数值实验,LBBD结合分枝定界算法的求解结果和LBBD结合Gurobi的求解结果如表4所示,即直接用Gurobi求解Benders分解的子问题而非本文提出的分枝定界算法.全部算例的规模为7-3-5,设置最大求解时间为 3 600 s.除此以外,为避免上下界接近时的拖尾效应,设置算法另一停止条件为上下界相对差距Gapb小于等于1%.考虑到Gurobi求解子问题速度较慢,设定其求解时间最长为 50 s.其中,上界Gap由(BUg-BUb)/BUg求得.

表4   精确求解数值实验结果

Tab.4  Numerical experiment results of exact algorithm

算例分布
类型
LBBD结合分枝定界LBBD结合Gurobi上界
Gap/%
上界,BUb下界,BLbGapb/%时间/s上界,BUb下界,BLbGapb/%时间/s
9N73.973.901 57589.943.4051.60360017.73
10N282.5282.50963306.338.5087.4036007.77
11N397.8397.801 378397.967.1083.1036000.03
12N338.7338.701 621380.264.0083.20360010.92
13N524.0524.001 940639.144.7093.00360018.01
14E400.3379.25.203 600416.5113.3072.8036003.90
15E471.8471.80163478.2195.8059.1036001.34
16E424.4420.31.001 280434.3150.7065.3036002.28
17E163.1161.90.701 228181.052.3071.1036009.87
18E582.5578.50.702 847674.6175.4074.00360013.65
19U77.477.402 880126.714.1088.90360038.91
20U96.296.20422110.525.7076.70360012.94
21U464.1464.102 078541.393.0082.80360014.26
22U47.447.4015247.80.0699.9036000.84
23U205.9205.90948236.712.5594.70360013.01

新窗口打开| 下载CSV


表4可知,首先,使用分枝定界算法求解子问题的LBBD算法,针对不同类型的算例、不同的分布类型,绝大部分算例能够在 3 600 s 内求得最优,不能求得最优的算例Gap也在6%以内,说明本算法适用于复杂的随机环境.同时,使用Gurobi求解子问题的LBBD算法在 3 600 s 难以求得最优,且上下界差距大,最小的Gap也在50%以上.其次,比较两种算法求得的上界,结合分枝定界的LBBD算法在所有算例中都取得了更好的结果,在15个算例中, 有8个算例的上界Gap达到10%以上.需要说明的是,7-3-5的问题规模已经很大,其中至少包含了对21个子批的排序问题,Gurobi在求解单次子问题时都难以求得最优.进一步的实验说明,使用Gurobi直接求解1.2节中随机MIP模型,在 3 600 s 内能最优求解的最大规模是4-3-5,与LBBD算法的求解规模有明显差距.综上所述,设计的LBBD算法与分枝定界算法具有较高的稳定性和求解效率.

4.3 启发式算法求解实验

对于大规模的算例,采用禁忌搜索求解子问题,并进行数值实验.为验证算法的有效性,将启发式算法的结果与4种企业实际策略、Gurobi、遗传算法[9]进行比较. 企业策略包括全分批最小交期策略、全分批最小整体松弛策略、全分批最小加工时间策略、不分批最小交期策略, 将4种企业实际策略编号为策略1~策略4.其中,全分批是指所有工件全部分为最小子批.遗传算法的染色体编码规则与文献[9]中相同,将延期期望作为适应度值,种群规模为 1 000,交叉概率为0.6,变异概率为0.4,设定算法终止条件为 3 600 s 运行时间.相应地,将LBBD启发式算法和Gurobi的最大求解时间也设置为 3 600 s,实验结果如表5所示.用BUmin表示4种企业实际策略中的最优结果,BUl为LBBD启发式算法的结果,BUg为Gurobi的结果,BUc为遗传算法[9]的结果,则Gapmin=(BUmin-BUl)/BUmin,Gapg=(BUg-BUl)/BUg,Gapc=(BUc-BUl)/BUc.

表5   启发式算法数值实验结果

Tab.5  Numerical experiment results of heuristic algorithm

算例分布LBBD企业实际方案Gapmin/
%
Gurobi遗传算法
求解
结果
时间/s策略1策略2策略3策略4时间/s求解
结果
Gapg/
%
时间/s求解
结果
Gapc/
%
时间/s
24N1 7003 6003 4503 6385 1884 986<550.722 06621.553 6002 90541.483 600
25N1 0223 6003 9814 2835 8503 525<571.011 1098.513 6002 48458.863 600
26N1 6083 6004 2564 3907 1815 437<562.222 06328.303 6002 84343.443 600
27N5023 6002 7712 6392 3542 986<578.6858817.193 6002 21177.313 600
28N2 0393 6005 7515 5457 0417 408<563.232 30813.193 6002 52619.283 600
29N2 8043 6005 5325 4475 1116 162<545.143 0458.593 6004 41236.453 600
30E2 0193 6002 7522 8844 0274 813<526.642 2069.263 6002 98132.273 600
31E9993 6001 8732 0002 2692 876<546.661 25926.033 6001 61738.223 600
32E8073 6001 0231 0511 3482 348<521.111 01625.903 6001 49846.133 600
33U1 3013 6002 7562 7564 0665 757<552.791 86543.353 6002 68651.563 600
34U2 1483 6003 8253 9435 1197 178<543.842 48615.743 6003 64541.073 600
35U1 0023 6001 8811 9622 1384 807<546.731 0827.983 6002 73863.403 600

新窗口打开| 下载CSV


表5可知,设计LBBD启发式算法明显优于企业实际策略、Gurobi、遗传算法[9].其中,企业实际策略虽然运行时间短,但结果较差,与LBBD启发式算法最小的差距Gapmin也达到45%.除此以外,LBBD启发式算法在与Gurobi和遗传算法求解时间相同的情况下,解的质量有大幅提升.在12个算例中,有8个算例与Gurobi给出的解相比,结果有10%以上的提高;在所有算例中,相比遗传算法解的质量均有19%以上的提高.这是由于将禁忌搜索嵌入Benders框架中,主问题能更好地考虑多种批量分配方案,而在遗传算法中,批量分配方案则是随机变化的,不容易搜索到好的分配方案,且解的质量具有较强的随机性.

表5中,算例24~29中交期服从正态分布,算例30~32中交期服从指数分布,算例33~35中交期服从均匀分布.在交期服从不同分布类型的同时,算例还选用了多种类型的数据,这些数据的差异体现在工件加工时间与换模时间的比例不同.实验结果证明,针对不同类型的数据和不同的交期分布类型,LBBD启发式算法的求解效果都具有很高的质量,体现出其具有较好的稳定性.

综上所述,表5中结果显示出LBBD启发式算法的良好优化效果,提出的研究方法对实际采用的策略有很大的提升空间,具有较好的实际应用潜力.

5 结语

针对工件交期随机且考虑换模时间的流水车间批量流调度问题,推导出交期服从均匀、指数、正态分布下延期期望的闭式表达式.基于此建立随机混合整数规划模型,针对模型高度非线性特征对其线性化,并设计Benders分解与分枝定界结合的精确算法,通过两阶段、重复枝去除等策略进行算法加速,以提升算法求解效率.通过与Gurobi的结果比较,表明本文精确算法的可行高效.针对大规模的算例,采用启发式算法,通过与4种企业实际加工策略及遗传算法比较,证明结合启发式的Benders算法在求解大规模问题上具有优越性.

附录见本刊网络版(xuebao.sjtu.edu.cn/article/2024/1006-2467/1006-2467-58-08-1271.shtml)

附录A 定理及引理证明

定理5 工件延期期望E(CPn)是单调递增的凸函数.

证明

 E(CPn)=ten-(ten-x)fn(x)dx,ten

引理1 对于均匀分布、指数分布、正态分布,An=limx-xFn(x)=0.

证明

(1) 对于均匀分布,令x =tn~U(an, bn),则∀xan, Fn(x)=0, 因此limx-xFn(x)=0.

(2) 对于指数分布,令x-an =dn-an~E(λn),则∀xan, Fn(x)=0, 因此limx-xFn(x)=0.

(3) 对于正态分布x=dn~N(μn, ), Fn(x)=fn(t)dt, 其中fn(x)=.故|xFn(x)|= \underline{\underline{\text{洛必达法则}}} =x2=0.

定理6 对于集合J中工件,其延期下界gLBJ=h=1JE(Rϕh,K-d[h])+.

证明 首先引入引理2~4如下.

引理2 对于所有未排程子批构成的任一序列σ-,存在一个序列σ-',其中同种工件的未排程子批是连续排列的,序列σ-'对应的延期下界gLBJ(σ-')≤gLBJ(σ-).

证明 考虑两种工件存在未排程的子批,记工件分别为A和B,φ(A)=1, φ(B)=2,其中,φ(i)表示集合J中工件i未排程的子批数量.未排程子批的某种排列ABA对应的延期下界记为gJLB(ABA),排列BAA对应的延期下界记为gLBJ(BAA), gLBJ(ABA)=E(R2K-dB)++E(R3K-dA)+, gLBJ(BAA)=E(R1K-dB)++E(R3K-dA)+.由于RcKc的方向单调递增,显然gLBJ(ABA)>gLBJ(BAA).任何存在不同种工件子批交叉的排列,如ABA,都可以通过交换两个子批位置的操作,得到同种工件子批连续的排列.如BAA,在一种工件延期不变的情况下降低另一种工件的延期,从而降低未完全排程工件的延期.以上证明拓展到3个以及以上的工件显然成立.

由于J中工件的未排程子批存在多种排列,由引理2得到的下界(对于所有未排程的子批,同种工件的子批连续排列)未必是节点本身的下界.考虑到任一工件i的延期为E(RtK-di)+,包含工件最后一个子批的位置s、工件的交期di两项,引理3和4分别对这两项分析下界构造.

基于引理2,将同种工件的未排程子批连续排列后,给定J中工件的任一排序ρ,未排程子批的排序也被确定.例如ρ(1)=3,则表示在已排程子批序列σ后,先连续排列第3种工件的所有未排程子批.φ(i)升序排列为φ[1],φ[2], …, φ[|J|], 令ϕh=i=1hφ[i],表示前hφ(i)的取值之和,可以得到引理3如下.

引理3 给定集合J中工件的任一排序ρ,h=1JE(Rϕh, K-dρ(h))+是此排序下延期的一个下界.

证明 由于Rckc的方向单调递增,所以引理3成立.例如有3种工件A、B、C未完全排程,其中φ(A)=1, φ(B)=2, φ(C)=3, 相应地,ϕ1=1, ϕ2=3, ϕ3=6.对于U中子批,某种连续排列的序列是CCCBBA,对应的延期下界:

T(CCCBBA)=E(R3K-dC)++E(R5K-dB)++E(R6K-dA)+E(R1K-dC)++E(R3K-dB)++E(R6K-dA)+=E(Rφ1, K-dC)++E(Rϕ2, K-dB)++E(Rϕ3, K-dA)+.

引理3用下标ϕh固定了延期表达式中RtK的数值,进一步考虑工件本身的排序,提出引理4如下.

引理4 对于任意给定的完工时间C1C2,若任意两个工件的随机交期满足d1std2,有E(C1-d1)++E(C2-d2)+E(C1-d2)++E(C2-d1)+.

证明f1(t)=E(t-d1)+, f2(t)=E(t-d2)+.

F1, F2分别为d1, d2 的分布函数.d1std2,可得对于任意t,F1(t)≥F2(t).第3章已证得f'1(t)=F1(t),f'2(t)=F2(t),因此对于任意给定的完工时间C1C2,f1(C1)-f1(C2)≤f2(C1)-f2(C2),E(C1-d1)++E(C2-d2)+E(C1-d2)++E(C2-d1)+得证.

由引理2和3,在推导gLB时,对于所有未排程的子批,同种工件的子批连续排列,给定集合J中工件的任一排序ρ,h=1JE(Rϕh, K-dρ(h))+是此排序下延期的一个下界.进一步考虑工件排序,由引理4,若dρ(i)stdρ(i+1),则E(Rϕi, K-dρ(i))++E(Rϕi+1, K-dρ(i+1))+E(Rϕi, K-dρ(i+1))++E(Rϕi+1, K-dρ(i))+.因此按照工件交期的随机序升序排列,可以得到延期下界h=1JE(Rϕh, K-d[h])+h=1JE(Rϕh, K-dρ(h))+.综上,对于J中工件,其延期下界gLBJ=h=1JE(Rϕh, K-d[h])+.

参考文献

GAREY M R, JOHNSON D S, SETHI R.

The complexity of flowshop and jobshop scheduling

[J]. Mathematics of Operations Research, 1976, 1(2): 117-129.

[本文引用: 1]

CHUNG C S, FLYNN J, KIRCA Ö.

A branch and bound algorithm to minimize the total tardiness for m-machine permutation flowshop problems

[J]. European Journal of Operational Research, 2006, 174(1): 1-10.

[本文引用: 2]

GMYS J, MEZMAZ M, MELAB N, et al.

A computationally efficient Branch-and-Bound algorithm for the permutation flow-shop scheduling problem

[J]. European Journal of Operational Research, 2020, 284(3): 814-833.

[本文引用: 1]

RAD S F.

A new ant algorithmic approach for solving PFSP

[J]. Iranian Journal of Science and Technology, Transactions A: Science, 2022, 46(1): 181-188.

[本文引用: 1]

黎阳, 李新宇, 牟健慧.

基于改进模拟退火算法的大规模置换流水车间调度

[J]. 计算机集成制造系统, 2020, 26(2): 366-375.

[本文引用: 1]

LI Yang, LI Xinyu, MOU Jianhui.

Large-scale permutation flowshop scheduling method based on improved simulated annealing algorithm

[J]. Computer Integrated Manufacturing Systems, 2020, 26(2): 366-375.

[本文引用: 1]

赵芮, 顾幸生.

求解零空闲流水车间调度问题的离散正弦优化算法

[J]. 上海交通大学学报, 2020, 54(12): 1291-1299.

DOI:10.16183/j.cnki.jsjtu.2019.321      [本文引用: 1]

针对以最小化最大完工时间(makespan)为目标的零空闲流水车间调度问题(NIFSP),提出一种离散正弦优化算法(DSOA)进行求解.受正弦波形的启发,原始的正弦优化算法(SOA)是一种利用正弦函数对个体位置进行更新的全局优化算法.首先,重新定义了适应组合优化问题的位置更新策略,采用一种去除工件数大小可变的迭代贪婪算法来对个体位置进行更新,以提高算法的探索能力.其次,采用了交叉操作和保留精英解的选择策略,避免算法陷入局部最优.最后,为了提高局部搜索的开发能力和算法精度,引入了一种基于插入的局部搜索方法,以便于在当前最优解的周围寻找更好的解.此外,基于Taillard基准,给出了算法性能比较的仿真结果,实验结果验证了所提出的DSOA算法求解NIFSP的有效性.

ZHAO Rui, GU Xingsheng.

A discrete sine optimization algorithm for No-idle flow-shop scheduling problem

[J]. Journal of Shanghai Jiao Tong University, 2020, 54(12): 1291-1299.

[本文引用: 1]

高丽, 周炳海, 杨学良, .

基于多规则资源分配的柔性作业车间调度问题多目标集成优化方法

[J]. 上海交通大学学报, 2015, 49(8): 1191-1198.

[本文引用: 1]

GAO Li, ZHOU Binghai, YANG Xueliang, et al.

A multi-objective integrated optimization method for FJSP based on multi-rule resource allocation

[J]. Journal of Shanghai Jiao Tong University, 2015, 49(8): 1191-1198.

[本文引用: 1]

汤洪涛, 王丹南, 邵益平, .

基于改进候鸟迁徙优化的多目标批量流混合流水车间调度

[J]. 上海交通大学学报, 2022, 56(2): 201-213.

DOI:10.16183/j.cnki.jsjtu.2020.435      [本文引用: 1]

针对2+1+1型混合流水车间,研究了多目标不相等批量流混合流水车间调度问题,提出一种基于变邻域搜索的自适应候鸟迁徙优化(AMBO)算法,实现了最小化完工时间与最小平均在制品数量的多目标优化.相比原始候鸟迁徙算法,AMBO算法引入变邻域搜索策略,实现每个算子的权重随迭代次数自适应调整,并提出了时间窗算子,以提升交换算子搜索性能和收敛速度.对随机生成不同规模的订单进行算例研究,结果表明AMBO算法比候鸟迁徙优化算法、遗传算法具有更高的求解质量和收敛性能,从而验证了AMBO算法的有效性.

TANG Hongtao, WANG Dannan, SHAO Yiping, et al.

A modified migrating birds optimization for multi-objective lot streaming hybrid flowshop scheduling

[J]. Journal of Shanghai Jiao Tong University, 2022, 56(2): 201-213.

[本文引用: 1]

CHEN T L, CHENG C Y, CHOU Y H.

Multi-objective genetic algorithm for energy-efficient hybrid flow shop scheduling with lot streaming

[J]. Annals of Operations Research, 2020, 290(1): 813-836.

[本文引用: 5]

PAN Y X, GAO K Z, LI Z W, et al.

Improved meta-heuristics for solving distributed lot-streaming permutation flow shop scheduling problems

[J]. IEEE Transactions on Automation Science and Engineering, 2023, 20(1): 361-371.

[本文引用: 1]

SANG H Y, PAN Q K, DUAN P Y, et al.

An effective discrete invasive weed optimization algorithm for lot-streaming flowshop scheduling problems

[J]. Journal of Intelligent Manufacturing, 2018, 29(6): 1337-1349.

[本文引用: 1]

TSENG C T, LIAO C J.

A discrete particle swarm optimization for lot-streaming flowshop scheduling problem

[J]. European Journal of Operational Research, 2008, 191(2): 360-373.

[本文引用: 1]

ALFIERI A, GLASS C, VAN DE VELDE S.

Two-machine lot streaming with attached setup times

[J]. IIE Transactions, 2012, 44(8): 695-710.

[本文引用: 1]

VILLARINHO P A, PANADERO J, PESSOA L S, et al.

A simheuristic algorithm for the stochastic permutation flowshop problem with delivery dates and cumulative payoffs

[J]. International Transactions in Operational Research, 2021, 28(2): 716-737.

[本文引用: 1]

LIU R, FAN X Y, WU Z R, et al.

The physician scheduling of fever clinic in the COVID-19 pandemic

[J]. IEEE Transactions on Automation Science and Engineering, 2022, 19(2): 709-723.

[本文引用: 1]

BELZUNCE F, MARTÍNEZ-RIQUELME C, MULERO J.

An introduction to stochastic orders

[M]. New York, American: Academic Press, 2015.

[本文引用: 1]

NAWAZ M, ENSCORE E E, HAM I.

A heuristic algorithm for the m-machine, n-job flow-shop sequencing problem

[J]. Omega, 1983, 11(1): 91-95.

[本文引用: 1]

RUIZ R, STÜTZLE T.

An Iterated Greedy heuristic for the sequence dependent setup times flowshop problem with makespan and weighted tardiness objectives

[J]. European Journal of Operational Research, 2008, 187(3): 1143-1159.

[本文引用: 1]

/