摘要：针对制造企业向多工厂分布式生产结构的转变，提出一个考虑订单运输时间的多目标分布式流水车间调度问题，并构建了一个以最小化最大完工时间与总加权提前和拖期惩罚为目标的混合整数规划模型，进而结合问题特点提出一种多目标头脑风暴优化算法进行求解。为了提高算法的搜索能力，利用支配规则设计了一种自适应聚类策略，并基于问题特征采用两种交叉方法和4种邻域结构。实验结果验证了所提算法的优越性。
关键词：分布式流水车间调度;订单运输时间;头脑风暴优化;多目标优化
论文《考虑订单运输时间的多目标分布式流水车间调度》发表在《计算机集成制造系统》，版权归《计算机集成制造系统》所有。本文来自网络平台，仅供参考。
0 引言
随着市场全球化的不断发展和客户差异化需求的持续增长，为了提高生产效率和产品多样性，企业逐渐采用智能生产系统进行加工制造。分布式制造是一种常见的智能生产模式，旨在通过在分散的地理位置上建立多个工厂进行协同生产，利用工厂间的信息交流、技术共享对有限的资源进行优化组合和合理配置，从而快速完成产品的生产制造[1]。分布式制造系统因其在提高生产效率、降低运营成本方面的重要作用而被广泛应用于各行业领域，如汽车制造、炼钢、手术室调度和食品化工工业[2]。近年来，以分布式制造为背景的分布式流水车间调度问题得到了广泛研究，旨在研究工件在工厂间的分配及工件在各工厂内机器上的加工顺序，以实现指定调度指标最优化。通过对分布式流水车间进行合理调度，帮助决策者制定科学高效的生产计划，采用正确的分工和协作方式促进企业在低成本、低风险下生产出高质量的产品。
受经济全球化、现代通信技术飞速发展和客户多样化需求的影响，企业面临着日益激烈的市场竞争。现代企业之间的竞争本质上是供应链与供应链之间的竞争，如何高效组织和协同供应链中各成员实现供需精准对接，不断提高生产效率、降低库存水平、快速响应多变的市场需求并获取较高的顾客满意度，已经成为供应链管理的重点和热点[4]。生产和分销环节在供应链的作用至关重要，研究表明，有效整合生产和分销流程可以帮助企业将总运营成本降低3%~20%[5]。因此，如何科学高效地组织和管理生产与分销环节成为供应链管理的重要内容。在传统的供应链管理中，生产和分销环节往往由不同的部门独立执行，产品在生产完成后需要先存储到仓库再配送到各地，显然这种方式在一定程度上增加了库存成本和不必要的等待时间，难以快速、及时地响应客户需求[6]。另外，在许多涉及时间敏感性产品或按订单交付的行业领域，如报纸、时装、快餐等，往往需要对生产和分销环节进行协同管理和优化[7]。因此，为了在激烈的市场竞争中获得优势，企业有必要在生产制造产品时考虑分销阶段的要求，如订单运输时间、配送路线等。
近年来，考虑订单运输时间的生产车间调度问题在工业和运筹学领域得到了广泛关注。LIU等[8]提出一个针对易腐产品的、考虑订单运输时间的生产车间调度问题，目标是最小化总加权订单交货时间，其将生产环节建模为单机，同时考虑多个客户和车辆容量等约束，设计了一个改进的大邻域搜索算法对问题进行求解;KERGOSIEN等[9]针对医疗卫生行业提出一个考虑订单运输时间的并行机调度问题，目标为最小化最大总延迟，该研究需要用车辆将完工产品交付给位于不同地理位置的客户，因此分销过程是一个多旅行商问题;KESEN等[10]提出一个考虑订单运输时间的不相关并行机调度问题，需要用一组车辆将完工产品在时间窗内配送给位于不同地理位置的客户，并采用CPLEX对该问题进行求解;MOHAMMADI等[11]探讨了更复杂的生产环境，提出一个考虑订单运输时间的流水车间调度问题，目标为最小化最大返回时间，并设计了一个改进的帝国主义竞争算法对问题进行求解;FU等[12]研究了一个以最小化最大完成时间为目标的、考虑订单运输时间的分布式流水车间调度问题，并结合问题特点提出一个增强的黑寡妇优化算法对问题进行求解。
上述研究只考虑了与时间指标相关的单目标优化问题，而在现实工业系统中，决策者往往需要同时考虑多种影响企业经济效益的因素，如客户服务水平、成本、能耗等。因此，近年来，考虑订单运输时间的多目标生产车间调度问题引起了学者们的广泛关注。FU等[13]针对一个金属包装企业提出一个考虑订单运输时间的多目标不相关并行机调度问题，其中生产阶段目标是最小化总安装成本，而分销阶段的目标是最小化运输成本，同时提出一个两阶段启发式算法对该问题进行求解;MOHAMMADI等[14]针对一家家具制造公司研究了一个考虑订单运输时间的多目标柔性作业车间调度问题，产品在生产阶段加工完成后需要由一组车辆配送，因此配送过程被建模为一个车辆路径优化问题，同时针对该问题提出一个目标为最小化总生产和分销成本、最小化总加权提前和拖期的双目标混合整数规划模型，并采用混合粒子群优化算法对模型进行求解;YAGMUR等[15]研究了一个考虑订单运输时间的多目标流水车间调度问题，并采用模因算法和模拟退火算法实现了总旅行时间最小化和总拖期最小化。
目前，考虑订单运输时间的生产车间调度问题主要具有以下特点：①大多数研究在生产环节只考虑单个工厂，且生产车间主要为单机、并行机、流水车间和作业车间等传统车间调度模型，而对生产环境较复杂的分布式车间调度研究较少;②分销阶段往往考虑单车场车辆路径优化问题，路径决策通常局限在一种理想状态中，而对配送环境比较复杂的多车场车辆路径优化问题关注较少;③由于考虑订单运输时间的生产车间调度问题的高度复杂性，智能优化算法在求解该类问题时取得了较好的性能，成为主流求解方法。基于上述对现有研究工作的总结，本文提出一个考虑订单运输时间的多目标分布式流水车间调度问题。综上，本文的创新点如下：
(1) 问题研究方面：本文提出一个考虑订单运输时间的多目标分布式流水车间调度问题，该问题包括生产阶段和分销阶段两个阶段。生产阶段是一个分布式流水车间调度问题，研究工件在工厂间的分配及其在工厂内机器上的加工顺序;分销阶段是一个多车场车辆路径优化问题，研究工件在车辆间的分配决策及客户的访问顺序。为了清晰地描述该问题，本文构建了一个以最小化最大完工时间与总加权提前和拖期惩罚为目标的多目标混合整数规划模型。
(2) 算法研究方面：本文设计了一个多目标头脑风暴优化(Multi-Objective Brain Storm Optimization, MOBSO)算法对问题进行求解。该算法提出一种自适应聚类策略确定聚类数目，同时采用两种交叉方法和4种邻域结构对解空间进行高效搜索，并采用非支配排序和拥挤度距离方法保存较优的个体作为新一代种群。通过与4种流行的多目标进化算法对一组测试算例进行对比实验表明，所设计的算法可以有效解决所提出的问题。
1 问题描述与数学模型
1.1 问题描述
本文所研究的问题描述如下：考虑将n个工件分配给f个工厂加工，每个工厂是一个包含m台串联机器的流水车间。工件加工完成后，由一组具有容量限制的车辆交付给位于不同地理位置的客户。因此，整个过程可以抽象为生产阶段和分销阶段两个阶段，生产阶段是一个分布式流水车间调度问题，而分销阶段可以描述为一个多车场车辆路径优化问题。为了求解该问题，生产阶段需要考虑工件在工厂间的分配及工件在工厂内机器上的加工顺序，分销阶段需要考虑工件在车辆间的分配和客户的访问顺序。需要说明的是，生产阶段的每个工件都对应分销阶段一个客户的订单所需要的产品，文中术语“工件”“客户”“订单”“产品”具有对应关系。每个订单都有时间窗限制，当车辆到达时间早于规定的到达时间或离开时间晚于规定的离开时间时，将产生提前和拖期惩罚。因此，本文旨在对生产和分销阶段进行合理决策，使最大完工时间与总加权提前和拖期惩罚最小。为了清晰地描述所研究的问题，给出一个问题示意图，并建立了一个混合整数规划模型。表1所示为模型中使用的符号、参数和决策变量。
表1 模型中使用的符号、参数和决策变量
类别 符号 说明
集合 F f 个工厂的集合, F={1,2,dots,f}
C n 个工件的集合, C={0,1,2,dots,n}, 其中0是虚拟工件
M m 台机器的集合, M={1,2,dots,m}
K 每个工厂中 s 辆车的集合, K={1,2,dots,s}
索引 g 工厂的索引符号, g in F
j,k,l 工件的索引符号, k,j,l in C
i 机器的索引符号, i in M
h 车辆的索引符号, h in K
参数 pt_{ij} 工件 j 在机器 i 上的加工时间
sigma_{gj} 工厂 g 和客户 j 之间的通行时间
ho_{kj} 客户 k 和客户 j 之间的通行时间
heta_j 工件 j 的重量
s_j 客户 j 的卸货时间
t'_j 客户 j 要求的最早到达时间
t''_j 客户 j 要求的最晚离开时间
varphi_j 客户 j 的单位提前惩罚权重
ilde{varphi}_j 客户 j 的单位拖期惩罚权重
xi 车辆的最大承载能力
lambda 一个很大的数
决策变量 X_{kjg} 若工件 k 先于工件 j 在工厂 g 中加工,则为1,否则为0
Y_{jg} 若工件 j 在工厂 g 中加工,则为1,否则为0
Z_{kjh} 若车辆 h 从客户 k 处行驶至客户 j 处,则为1,否则为0
C_{ij} 工件 j 在机器 i 上的完工时间
V_h 车辆 h 的配送开始时间
A_j 到达客户 j 的时间
L_j 离开客户 j 的时间
1.2 数学模型
一个可行解必须满足以下条件：①在0时刻，所有机器都可用;②每台机器一次只能加工一个工件;③每个工件每次只能在一台机器上加工;④工件一旦开始加工不得中断;⑤每个客户只允许被访问一次;⑥所有车辆的载重不能超过其最大载重。基于上述符号定义，建立如下模型：
目标函数：
约束条件：
其中：式(1)和式(2)分别为最小化最大完工时间及最小化总加权提前和拖期惩罚;式(3)和式(4)分别为每个工件在工厂中只有一个前置工件和一个后置工件;式(5)规定每个工件只能分配给一个工厂;式(6)保证一个工件最多同时有一个前置工件和一个后置工件;式(7)描述两个相邻工件之间的顺序关系;式(8)表示一个工件一次只能在一台机器上加工;式(9)保证一台机器一次只能加工一个工件;式(10)表示每个工厂的车辆数目是有限的;式(11)和式(12)定义每个客户只能被访问一次;式(13)表示车辆到达客户后立即离开;式(14)表示每辆车的载重不能超过其最大承载能力;式(15)规定车辆最多离开其相应的工厂一次;式(16)定义车辆的出发时间;式(17)和式(18)表示车辆的到达时间;式(19)表示车辆的离开时间;式(20)~式(24)定义了决策变量的取值范围。
2 多目标头脑风暴优化算法
因为分布式流水车间调度问题[3]和多车场车辆路径优化问题[16]已被证明是NP-hard，所以本文所考虑的问题也是NP-hard。为了在可接受的时间内找到上述问题的满意解，有必要利用智能优化算法对其进行求解。头脑风暴优化(Brain Storm Optimization, BSO)算法[17]是一种模拟人类头脑风暴过程产生新思想的智能优化算法，其主要迭代聚类、新一代个体生成和选择3个重要的搜索步骤。首先，聚类操作将种群划分为多个聚类，每个聚类包含一组个体，其中每个聚类中最好的个体称为中心个体，其余称为正常个体;然后，采用新一代个体生成操作对中心个体和正常个体进行全局和局部搜索得到一组新个体，其中每个新个体均与种群中的一个个体相关联;最后，通过对每对关联的个体进行选择操作，保留较好的个体作为下一代种群个体。
BSO算法在概念上理解简单，操作上易于实现，近年来已被广泛应用于求解各领域复杂优化问题，如电力系统、航空领域和车辆路径优化等[18-19]。本文将非支配排序思想与BSO算法进行深度融合，提出MOBSO算法对问题进行求解。
2.1 解的表达与初始化种群
本文采用二维向量进行编码[20]，个体用 x=(pi;omega) 表示，其中 pi=(pi_1,pi_2,dots,pi_n) 为一个工件序列向量，表示n个工件的调度顺序;omega=(omega_1,omega_2,dots,omega_n) 为一个工厂分配向量，表示pi中每个位置上工件对应的工厂索引。需要注意的是，这种编码方法仅给出了生产阶段工件的调度顺序和工厂分配决策，无法获得工件在机器上的加工顺序以及分销阶段车辆的工件分配情况与客户的访问顺序。因此，本文提出以下规则将个体解码为一个完整的方案：①工厂中的工件在机器上的加工顺序与其工厂分配顺序一致;②工厂中的所有工件按照其加工顺序依次选择车辆进行装载，车辆的工件载重量不得超过车辆的最大承载能力;③客户的访问顺序与其车辆所配送工件的装载顺序一致。另外，为了构造一个具有较好多样性的初始种群，本文采用随机的方法产生了Psi个个体作为初始种群。
2.2 聚类操作
BSO算法及其变体通常采用k-means方法[21]将种群划分为不同的聚类，该方法中聚类的数目是固定的。由于参数设置通常具有问题依赖性，本文引入一种新的聚类策略，以自适应地确定每次迭代的聚类数。受非支配排序思想[22]的启发，对所有个体按照其支配关系进行排序，设置具有相同排序值的个体属于一个聚类。在该方法中，个体的排序值越小，性能越好。因此，第1个聚类由最好的个体组成，而最差的个体则被用来构造最后一个聚类。
2.3 新个体生成操作
MOBSO算法对随机从两个聚类中产生的两个父代个体进行交叉操作，以产生新个体。针对解的编码特点，本文采用顺序交叉和两点交叉两种交叉方法，其中顺序交叉方法用于产生新个体的工件序列向量，两点交叉方法用于产生个体的工厂分配向量。图2所示为具有8个工件和3个工厂的两个父代个体进行交叉产生新个体的操作过程。首先，在父代个体x_1的工件序列向量中随机生成两个位置索引p_1和p_2，在工厂分配向量中随机生成两个位置索引p_3和p_4;然后，将x_1中p_1和p_2两端的工件信息直接复制到新个体x'相同的位置，同样将p_3和p_4两端的工厂信息直接复制到x'相同的位置;接着，将x'中缺失的工件信息按照其在父代个体x_2中出现的顺序依次填加到工件序列的空缺位置，同时将个体x_2中p_3和p_4之间的工厂分配信息依次填加到工厂分配决策的空缺位置;最后，通过交换两个父代个体x_1和x_2的信息构建两个新个体。本文根据新个体之间的支配关系，选择性能较好的个体作为最终的新个体。
为了充分利用一个个体的良好信息，本文引入两种搜索方法。第1种方法将模拟退火机制应用于一个随机从聚类中选择的中心个体上，其伪代码如下：
在第2种方法中，本文设计了4种问题依赖的邻域结构[23]用于搜索一个随机从聚类中选择的正常个体的周围邻域。4种邻域结构分别为关键工厂内部交换、关键工厂内部插入、关键工厂外部交换和关键工厂外部插入。对于每个被选中的正常个体，本文随机选择一种邻域结构来产生一个新的邻域解，并根据正常个体和邻域解之间的支配关系存储较好的个体作为新个体。本文定义具有最大完工时间的工厂为关键工厂，且4种邻域结构描述如下：
(1) 关键工厂内部交换：随机交换关键工厂中的两个工件。
(2) 关键工厂内部插入：随机将关键工厂中的一个工件插入另一个工件的前面。
(3) 关键工厂外部交换：随机将关键工厂中的一个工件与其他工厂中的一个工件进行交换。
(4) 关键工厂外部插入：随机选择关键工厂中的一个工件，并将其插入其他工厂。
2.4 选择操作
为了确保算法向有希望的领域加速搜索，并在逼近性和分布性之间取得较好的平衡，本文采用带精英策略的快速非支配排序遗传算法(Non-dominated Sorting Genetic Algorithm II, NSGA-II)中使用的非支配排序和拥挤度距离方法[22]设计了一种新的选择操作，过程为：首先将种群中个体和新产生的个体合并，对所有个体进行非支配排序并计算拥挤度距离;然后优先选择非支配排序值小的个体，当非支配排序值相同时，选择拥挤度距离值大的个体。按照此规则，共选择Psi个个体作为下一代种群。为了清楚地展示MOBSO算法，给出算法的流程图如图3所示。
MOBSO算法伪代码如下：
3 仿真实验与结果分析
3.1 测试算例和比较算法
本文从VRF集[24]中选择6个算例用于获取分布式流水车间调度问题的参数信息，分别为VRF30_5_1, VRF30_10_1, VRF60_5_1, VRF60_10_1, VRF100_20_1, VRF200_20_1，其中工件数量集合为{30,60,90,120}，机器数量集合为{5,10}，本文设置工厂数量集合为{2,3,4}。本文从Gehring&Homberger集[25]中选择一个名为C1_2_1的算例，获取了30,60,90,120个客户的坐标、时间窗、服务时间等信息。通过对选择的算例进行组合，共构建了24个测试算例。为了验证MOBSO算法的性能，本文选择NSGA-II[22]、基于分解的多目标进化算法(Multi-Objective Evolutionary Algorithm based on Decomposition, MOEA/D)[26]、多目标鲸鱼群算法(Multi-Objective Whale Swarm Algorithm, MOWSA)[23]和改进的Jaya算法(Improved Jaya Algorithm, IJA)[27]作为竞争算法进行实验。所有算法均采用C++程序语言编译，运行环境为Intel Core i5 CPU/3.2 GHz主频。
3.2 参数设置和评价指标
由于参数设置会影响算法的搜索性能，本文对MOBSO, NSGA-II, MOEA/D和MOWSA采用一个具有60个工件、5台机器和3个工厂的中等规模算例进行正交实验[28]。通过分析实验结果，算法的参数设置具体为：MOBSO算法的种群规模Psi=80，生成子代的参数为r_s=0.6, r_o=0.4, r_t=0.8;NSGA-II的种群规模为80，变异率为0.15;MOEA/D的种群规模为40，邻域规模为25;MOWSA的种群规模为20，交叉率为0.85，变异率为0.15。由于IJA只有一个参数，本文参考IJA原文[27]设置种群规模为20。
为了公平比较5种算法，本文采用估值次数作为停止条件，设置最大估值次数为150cdot ncdot m，其中n和m分别为工件和机器的数量。本文选用IGD指标[29]和HV指标[30]对实验结果进行分析。
IGD(S^, S)=frac{1}{S^
}sum_{sin S^*}distance(s, S) quad (25)
式中：distance(s, S)表示最优解集S^中的点s到某一特定算法获得的非支配解集S的最近距离;$S^
表示S^中解的数目。需要注意的是，最优解集S^$是合并所有算法获得的非支配解集并剔除支配解得到的。
HV指标计算了某一特定算法获得的非支配解集中的个体与参考点P^在目标空间中围成的超立方体的体积，本文设置P^为(1,1)。
3.3 实验结果与分析
为了验证MOBSO算法在求解所提问题时的性能，本文采用IGD指标和HV指标对5种算法的实验结果进行分析。对于每个测试算例，MOBSO, NSGA-II, MOEA/D, MOWSA, IJA 5种算法分别进行20次独立求解，计算20次运行结果的平均值和方差。
表2 MOBSO算法及其竞争算法关于IGD指标的比较结果
算例 MOBSO NSGA-II MOEA/D MOWSA IJA
f m n 均值 方差 均值 方差 均值 方差 均值 方差 均值 方差
2 5 30 0.1054 0.0006 0.2614 0.0018 0.2398 0.0039 0.139 0.0016 0.2846 0.0022
2 5 60 0.1099 0.0003 0.2410 0.0016 0.2657 0.0029 0.1588 0.0006 0.2903 0.0007
2 5 90 0.0930 0.0004 0.2645 0.000 0.3060 0.004 0.1744 0.0010 0.2865 0.0011
2 5 120 0.0756 0.0004 0.312 0.0024 0.3558 0.0050 0.2017 0.0013 0.3382 0.0039
2 10 30 0.1907 0.0023 0.3717 0.0063 0.4007 0.0373 0.2070 0.0027 0.2838 0.0020
2 10 60 0.1021 0.0006 0.3777 0.0022 0.4247 0.0367 0.2295 0.0017 0.3804 0.0024
2 10 90 0.131 0.0008 0.4975 0.0107 0.7583 0.0660 0.2603 0.003 0.3997 0.0046
2 10 120 0.098 0.0003 0.432 0.0027 0.5788 0.0446 0.2407 0.0016 0.3523 0.0007
3 5 30 0.2341 0.0009 0.1772 0.0046 0.1720 0.006 0.2377 0.0005 0.3237 0.0010
3 5 60 0.101 0.0003 0.187 0.0003 0.2218 0.0009 0.115 0.0003 0.2133 0.0003
3 5 90 0.0767 0.0003 0.180 0.0006 0.2020 0.0015 0.1254 0.0002 0.2003 0.0002
3 5 120 0.082 0.0002 0.1941 0.0005 0.2436 0.0010 0.1433 0.0003 0.2083 0.0004
3 10 30 0.1439 0.0017 0.2441 0.0014 0.236 0.0022 0.1539 0.0008 0.279 0.0023
3 10 60 0.0924 0.0004 0.2120 0.0006 0.255 0.0009 0.1497 0.0004 0.2647 0.0005
3 10 90 0.085 0.0002 0.260 0.0009 0.320 0.0011 0.1687 0.0007 0.276 0.0007
3 10 120 0.0747 0.0001 0.2551 0.0006 0.315 0.0015 0.1698 0.0005 0.2868 0.0002
4 5 30 0.3775 0.0011 0.1558 0.0068 0.2069 0.0129 0.3778 0.0005 0.4694 0.0011
4 5 60 0.1015 0.0002 0.1969 0.0005 0.2308 0.0010 0.1155 0.0003 0.2632 0.0005
4 5 90 0.0872 0.0006 0.2456 0.0009 0.2888 0.0007 0.1326 0.0002 0.305 0.004
4 5 120 0.0957 0.0006 0.2475 0.0003 0.2786 0.0007 0.1447 0.0002 0.280 0.0003
4 10 30 0.2643 0.0017 0.1838 0.0056 0.2591 0.0080 0.2324 0.0005 0.314 0.0007
4 10 60 0.1048 0.0004 0.2239 0.0006 0.2547 0.0011 0.1495 0.0003 0.301 0.0005
4 10 90 0.0712 0.0001 0.2090 0.0005 0.2493 0.0011 0.1535 0.0004 0.252 0.0003
4 10 120 0.0694 0.0001 0.2166 0.0006 0.2703 0.0010 0.1447 0.0003 0.252 0.0005
平均值 0.1237 0.0006 0.2562 0.0022 0.3057 0.0101 0.1803 0.0008 0.294 0.0010
表2所示为MOBSO算法及其竞争算法关于IGD指标的结果。可见，MOBSO算法能够在多数算例上获得较好的结果。具体而言，在24个测试算例中，MOBSO算法有21个算例的均值优于NSGA-II和MOEA/D，有23个算例的均值优于MOWSA，在24个算例上获得的均值均优于IJA。另外，MOBSO算法在多数算例上获得的方差均小于其竞争算法，说明其在解决所研究问题方面具有较好的稳定性。表2所示为所有算例关于IGD指标均值和方差的平均值，其中MOBSO, NSGA-II, MOEA/D, MOWSA, IJA的平均均值分别为0.1237, 0.2562, 0.3057, 0.1803, 0.2961，平均方差分别为0.0006, 0.0022, 0.0101, 0.0008, 0.0011，因此MOBSO算法比其竞争算法具有优势，从而验证MOBSO算法在解决所研究问题上具有较好的性能。
表3 MOBSO算法及其竞争算法关于HV指标的比较结果
算例 MOBSO NSGA-II MOEA/D MOWSA IJA
f m n 均值 方差 均值 方差 均值 方差 均值 方差 均值 方差
2 5 30 0.8018 0.0051 0.5546 0.0048 0.5833 0.0108 0.7650 0.0041 0.503 0.0078
2 5 60 0.8141 0.0035 0.5451 0.0048 0.5324 0.0067 0.6917 0.0039 0.426 0.0045
2 5 90 0.8114 0.0051 0.4794 0.0051 0.4184 0.0120 0.6589 0.0035 0.4216 0.0080
2 5 120 0.8652 0.0030 0.4929 0.0042 0.4195 0.0078 0.6867 0.0024 0.324 0.0100
2 10 30 0.7935 0.0077 0.4362 0.0096 0.4223 0.0217 0.7266 0.0077 0.5235 0.0072
2 10 60 0.8417 0.0039 0.4229 0.0035 0.3787 0.0168 0.6892 0.0059 0.4251 0.0080
2 10 90 0.864 0.0064 0.3203 0.0126 0.1495 0.0158 0.7071 0.0054 0.4461 0.0120
2 10 120 0.8518 0.0058 0.3058 0.0088 0.1756 0.0105 0.6339 0.0080 0.3957 0.0048
3 5 30 0.5572 0.0048 0.6813 0.0150 0.7015 0.0316 0.5551 0.0060 0.3515 0.0042
3 5 60 0.7421 0.0022 0.5409 0.0022 0.4957 0.0041 0.6479 0.0034 0.447 0.0031
3 5 90 0.7346 0.0035 0.5184 0.0033 0.4741 0.0050 0.6268 0.0015 0.465 0.0010
3 5 120 0.7878 0.0036 0.5311 0.0026 0.4353 0.0054 0.6423 0.0027 0.466 0.0030
3 10 30 0.7334 0.0126 0.5786 0.0096 0.5625 0.0102 0.7268 0.0065 0.4594 0.0080
3 10 60 0.7806 0.0027 0.5363 0.0034 0.4680 0.0053 0.6582 0.0021 0.425 0.0031
3 10 90 0.7722 0.0035 0.4569 0.0030 0.3717 0.0029 0.6039 0.0032 0.390 0.0029
3 10 120 0.7651 0.0034 0.4482 0.0034 0.3717 0.0051 0.6130 0.0025 0.393 0.0044
4 5 30 0.4672 0.0012 0.8700 0.0012 0.8743 0.0024 0.4575 0.0015 0.265 0.0016
4 5 60 0.7294 0.0037 0.5371 0.0032 0.485 0.0040 0.6568 0.0050 0.4253 0.0032
4 5 90 0.7251 0.0030 0.4981 0.0027 0.4350 0.0028 0.6218 0.0019 0.4241 0.4241
4 5 120 0.7471 0.0043 0.4798 0.0020 0.4192 0.0026 0.6191 0.0018 0.4156 0.0019
4 10 30 0.5763 0.0035 0.7959 0.0061 0.7987 0.0048 0.5795 0.0025 0.3241 0.0048
4 10 60 0.7690 0.0030 0.5372 0.0031 0.4979 0.0036 0.6767 0.0022 0.3846 0.0044
4 10 90 0.7361 0.0039 0.4764 0.0026 0.4176 0.0050 0.6023 0.0039 0.3899 0.0028
4 10 120 0.7608 0.0032 0.4779 0.0018 0.3987 0.0035 0.5984 0.0036 0.4026 0.0028
平均值 0.7512 0.0043 0.5217 0.0049 0.4703 0.0084 0.6436 0.0038 0.4123 0.0048
表3所示为5种算法的HV指标的比较结果。可见，在求解的24个测试算例中，MOBSO算法在21个算例上优于NSGA-II和MOEA/D，在23个算例上优于MOWSA，在全部算例上优于IJA。另外，通过比较MOBSO算法和其竞争算法的方差可见，MOBSO算法更有潜力获得稳定的非支配解集。通过分析MOBSO，NSGA-II，MOEA/D，MOWSA，IJA关于HV指标的均值和方差的平均值，其平均均值分别为0.7512, 0.5217, 0.4703, 0.6436, 0.4123，平均方差值分别为0.0043, 0.0049, 0.0084, 0.0038, 0.0048，与4种竞争算法相比，MOBSO算法在求解所有算例上更具优势，因此其在求解所研究的问题时具有较好的性能。
为了直观地展示5种算法的实验结果，本文绘制了不同规模算例的箱线图。图4和图5所示分别为MOBSO算法及其竞争算法关于IGD指标和HV指标的箱线图，可见MOBSO算法比其竞争算法获得的IGD指标和HV指标具有明显的优势，意味着MOBSO算法更有潜力发现更好、更集中的结果。图6所示为MOBSO算法及其竞争算法在4个算例中得到的非支配解集的分布情况，可见MOBSO算法获得的非支配解集具有较好的逼近性和分布性。基于上述结果和分析，验证了MOBSO算法在求解所提问题时具有较好的性能。
3.4 统计检验与分析
为了验证MOBSO算法和竞争算法之间的性能差异，本文采用Wilcoxon符号秩检验方法[31]对5种算法获得的实验结果在0.05显著性水平下进行统计检验。表4所示为MOBSO算法与其竞争算法两两比较获得的Wilcoxon符号秩检验结果，其中“w+/w~/w-”表示MOBSO算法有w+, w~, w-个算例分别优于、等价于和差于其竞争算法;“R+”表示MOBSO算法优于其竞争算法的秩和，“R-”表示其竞争对手优于MOBSO算法的秩和。由表可知，MOBSO算法在与其竞争算法进行两两比较时获得的R+值远大于R-值，且在95%的置信度下，IGD指标和HV指标获得的p-值均小于0.05。因此，Wilcoxon符号秩检验结果表明MOBSO算法性能显著优于其他4种竞争算法。
表4 Wilcoxon符号秩检验统计分析结果
指标 MOBSO vs. NSGA-II MOBSO vs. MOEA/D MOBSO vs. MOWSA MOBSO vs. IJA
IGD w+/w~/w-: 21/0/3 w+/w~/w-: 21/0/3 w+/w~/w-: 23/0/1 w+/w~/w-: 24/0/0
R+: 277, R-: 23 R+: 285, R-: 15 R+: 293, R-: 7 R+: 300, R-: 0
z: 3.63, p: 0.0001 z: 3.86, p: 0.000 z: 4.09, p: 0 z: 4.29, p: 0
HV w+/w~/w-: 21/0/3 w+/w~/w-: 21/0/3 w+/w~/w-: 23/0/1 w+/w~/w-: 24/0/0
R+: 272, R-: 28 R+: 275, R-: 25 R+: 298, R-: 2 R+: 300, R-: 0
z: 3.49, p: 0.0002 z: 3.57, p: 0.000 z: 4.23, p: 0 z: 4.29, p: 0
3.5 选择算子对比分析
为了检验所设计的选择算子对MOBSO算法搜索性能的影响，本文设计了一个采用基本BSO算法选择算子的竞争算法，记为MOBSO-S。MOBSO-S的选择操作描述为：设每个生成的子代都与当前种群中的一个个体有一一对应关系，根据两个相关联的个体之间的支配关系，保留性能较好的个体。MOBSO与MOBSO-S的实验结果如表5所示。由IGD指标值可知，MOBSO算法在24个算例中获得的均值均优于MOBSO-S，且在24个算例中有19个算例的方差优于MOBSO-S。由HV指标值可知，MOBSO算法在24个算例上获得的均值和方差均优于MOBSO-S，因此MOBSO算法优于MOBSO-S。另外，通过观察所有算例上的平均均值和平均方差值，MOBSO算法也优于MOBSO-S。因此，可以得出结论，MOBSO算法的性能比MOBSO-S更优越，本文设计的选择算子对MOBSO算法的搜索性能具有积极的作用。
表5 MOBSO算法与MOBSO-S关于IGD指标和HV指标的比较结果
算例 IGD HV
MOBSO 均值 方差 MOBSO-S 均值 方差 MOBSO 均值 方差 MOBSO-S 均值 方差
2 5 30 0.1082 0.0007 0.2453 0.0015 0.9280 0.0020 0.6854 0.0184
2 5 60 0.1131 0.000 0.3076 0.0026 0.9919 0.0013 0.6887 0.0257
2 5 90 0.0909 0.0007 0.2837 0.0027 1.0000 0 0.6094 0.0224
2 5 120 0.1002 0.0006 0.3354 0.0010 1.0000 0 0.6386 0.0231
2 10 30 0.2390 0.0117 0.4371 0.0184 0.4592 0.0346 0.3070 0.0393
2 10 60 0.1409 0.0016 0.3488 0.0023 0.9998 0 0.7447 0.0397
2 10 90 0.2137 0.0036 0.4391 0.0107 1.0000 0 0.6996 0.0236
2 10 120 0.1313 0.0024 0.3320 0.0026 0.9933 0.0009 0.7399 0.0291
3 5 30 0.1166 0.0010 0.2835 0.0006 1.0000 0 0.7183 0.0045
3 5 60 0.0877 0.0002 0.2697 0.0009 1.0000 0 0.6377 0.0094
3 5 90 0.0674 0.0003 0.2315 0.0004 1.0000 0 0.6415 0.0078
3 5 120 0.0724 0.0001 0.2813 0.0007 1.0000 0 0.6297 0.0151
3 10 30 0.1545 0.0031 0.2742 0.0053 0.6143 0.0227 0.5050 0.0335
3 10 60 0.1004 0.000 0.2978 0.0001 1.0000 0 0.6880 0.0184
3 10 90 0.0706 0.0005 0.3508 0.0013 1.000 0 0.7380 0.0049
3 10 120 0.0852 0.000 0.3161 0.0008 1.0000 0 0.7214 0.0137
4 5 30 0.1055 0.0009 0.2601 0.0017 1.0000 0 0.6973 0.0122
4 5 60 0.0890 0.0004 0.2567 0.0011 1.0000 0 0.6385 0.0119
4 5 90 0.0944 0.0010 0.332 0.000 1.000 0 0.6028 0.602
4 5 120 0.1047 0.0014 0.353 0.000 1.0000 0 0.5519 0.5519
4 10 30 0.1163 0.0023 0.2160 0.0020 0.7553 0.0032 0.5812 0.5812
4 10 60 0.1016 0.0007 0.2825 0.0008 1.0000 0 0.6899 0.6899
4 10 90 0.0686 0.0003 0.2951 0.0010 1.000 0 0.7083 0.7083
4 10 120 0.0675 0.0002 0.2732 0.000 1.0000 0 0.7123 0.7123
平均值 0.1100 0.0015 0.3043 0.0025 0.947 0.0027 0.6490 0.6490
3.6 MOBSO与CPLEX对比分析
为了验证算法的寻优性能，本文将MOBSO算法与一个数学规划优化软件CPLEX进行比较。表6和表7所示为12个工件的信息，本文分别利用前8个工件、前10个工件和所有工件构建了3个测试问题，每一问题有两个工厂，每个工厂有2台机器和2辆车。3个测试问题的工件数量分别为8,10,12时，对应车辆的最大载重分别为60,80,100，且每个客户的卸货时间为90。由于CPLEX不能直接求解多目标优化问题，本文采用加权求和方法将所提问题转化为单目标优化问题进行求解。针对每个测试问题，本文分别用3个权重向量，即(1.0,0.0),(0.5,0.5),(0.0,1.0)，构建了3个单目标优化子问题的测试算例。令CPLEX的停止条件为3600s，若CPLEX不能在3600s内获得全局最优解，则输出发现的近似最优解。表8所示为CPLEX和MOBSO算法的实验结果，其中“*”表示CPLEX获得的结果为全局最优解。可见，对于具有8个工件的测试问题，CPLEX可以在给定时间内获得最优解，而MOBSO算法在一个算例上获得了最优解;对于具有10个工件和12个工件的测试问题，CPLEX在给定时间内获得了近似最优解，而MOBSO算法能够在较短的时间内获得更好的结果。因此，可以得出结论，当问题规模较大时，MOBSO算法的求解速度和求解质量均优于CPLEX。另外，通过对比MOBSO算法与数学规划优化软件CPLEX，不但验证了MOBSO算法在求解所研究问题时具有较好的性能，而且证明了模型的合理性和算法的有效性。
表6 数值算例的参数信息
工件 加工时间 (mu_1, mu_2) 时间窗 (t', t'') 重量 ( heta_j) 提前/拖期惩罚 (varphi_j, ilde{varphi}_j)
J1 36, 45 105, 130 20 0.2, 0.4
J2 79, 82 193, 253 20 0.3, 0.2
J3 21, 19 162, 192 10 0.2, 0.3
J4 16, 82 190, 294 10 0.4, 0.3
J5 23, 57 141, 160 20 0.2, 0.3
J6 31, 73 157, 192 10 0.3, 0.2
J7 29, 90 234, 300 10 0.2, 0.3
J8 95, 3 311, 345 40 0.1, 0.2
J9 45, 30 308, 353 20 0.3, 0.2
J10 50, 93 278, 364 20 0.2, 0.4
J11 54, 86 332, 378 20 0.3, 0.5
J12 54, 20 371, 400 10 0.4, 0.5
表7 工厂和客户之间及客户和客户之间的通行时间
F1 F2 J1 J2 J3 J4 J5 J6 J7 J8 J9 J10 J11 J12
F1 999 999 52 70 55 76 94 54 23 45 44 12 44 113
F2 999 999 43 47 53 70 86 60 35 53 18 44 37 46
J1 52 43 999 89 94 113 129 24 33 75 22 48 76 70
J2 70 47 89 999 35 37 46 107 37 45 22 67 53 75
J3 55 53 94 35 999 21 39 105 12 18 16 29 53 70
J4 76 70 113 37 21 999 18 125 43 27 36 93 33 27
J5 94 86 129 46 39 18 999 143 28 39 47 56 35 39
J6 54 60 24 107 105 125 143 999 25 36 48 69 23 36
J7 23 35 33 37 12 43 28 25 999 12 33 44 53 12
J8 45 53 75 45 18 27 39 36 12 999 23 34 107 999
J9 44 18 22 22 16 36 47 48 33 23 999 84 18 23
J10 12 44 48 67 29 93 56 69 44 34 84 999 24 34
J11 94 37 76 53 53 33 35 23 53 107 18 24 999 107
J12 113 46 70 75 70 37 53 35 75 37 22 107 17 37
表8 MOBSO算法与CPLEX实验结果
n 权重 CPLEX MOBSO
目标值 运行时间/s 目标值 运行时间/s
8 (1.0,0.0) 246.000* 25.140 246.000 0.956
(0.5,0.5) 292.050* 2242 300.850 0.829
(0.0,1.0) 322.700* 245 327.900 0.853
10 (1.0,0.0) 307.000 3600 307.000 1.044
(0.5,0.5) 428.800 3600 426.700 1.013
(0.0,1.0) 499.300 3600 498.300 1.020
12 (1.0,0.0) 360.000 3600 360.000 1.288
(0.5,0.5) 593.650 3600 565.750 1.298
(0.0,1.0) 767.800 3600 761.500 1.298
通过分析上述仿真实验和统计检验结果可知，MOBSO算法在求解所研究问题时能够获得较高质量的非支配解集，其优越的性能主要得益于良好的全局搜索和局部搜索方法。在搜索过程中，本文采用顺序交叉和两点交叉方法将算法的搜索方向转移到可能具有更高搜索潜力的区域，从而增强算法的全局搜索能力;然后，本文设计了一种模拟退火机制和4种依赖问题的邻域结构帮助算法优化搜索区域，从而增强算法的局部搜索能力;同时，考虑到聚类数目参数具有问题依赖性，应随算法的进化过程不断发生变化，本文基于支配规则设计了一种自适应聚类策略，以动态地将种群划分到不同的聚类中;为了获得一个高质量的种群，本文将子代种群加入种群，并在所有个体中保留良好的解作为下一代种群。通过迭代执行上述操作，使MOBSO算法具有较强的能力，保证算法的全局搜索过程且避免算法陷入局部最优，因此在24个测试问题上的求解结果比其竞争算法表现得好，MOBSO算法在求解所研究问题时可以看作是一个出色的求解器。
4 结束语
本文研究了一个考虑订单运输时间的多目标分布式流水车间调度问题，建立了一个以最小化最大完工时间与总加权提前和拖期惩罚为目标的多目标混合整数规划模型。鉴于所考虑的问题具有NP-hard的性质，本文提出一种MOBSO算法进行求解。该算法引入一种自适应聚类策略帮助算法在进化过程中确定聚类数目，并采用两种交叉方法、一种模拟退火机制和4种问题依赖的邻域结构来提高算法的搜索能力。为了验证该算法的性能，本文选取NSGA-II, MOEA/D, MOWSA和IJA作为竞争算法，并且对24个测试算例进行仿真实验和统计检验，结果表明MOBSO算法在求解所研究问题时获得了较好的效果。
未来工作将侧重于考虑各种实际加工约束，如能耗、质量、成本等，并进一步探讨机器学习与BSO算法相融合的学习机制。
参考文献
[1] WANG Ling, DENG Jin, WANG Shengyao. Survey on optimization algorithms for distributed shop scheduling[J]. Control and Decision, 2016, 31(1): 1-11 (in Chinese). [王凌，邓瑾，王圣尧. 分布式车间调度优化算法研究综述[J]. 控制与决策，2016, 31(1): 1-11.]
[2] FU Y P, HOU Y S, WANG Z F, et al. Distributed scheduling problems in intelligent manufacturing systems[J]. Tsinghua Science and Technology, 2021, 26(5): 625-645.
[3] FU Y P, TIAN G D, FATHOLLAHI-FARD A M, et al. Stochastic multi-objective modelling and optimization of an energy-conscious distributed permutation flow shop scheduling problem with the total tardiness constraint[J]. Journal of Cleaner Production, 2019, 226: 515-525.
[4] CHEN Jian, LIU Yunhui. Operations management innovation enabled by digitalization and intellectualization: From supply chain to supply chain ecosystem[J]. Journal of Management World, 2021, 37(11): 227-240, 14 (in Chinese). [陈剑, 刘运辉. 数智化使能运营管理变革: 从供应链到供应链生态系统[J]. 管理世界, 2021, 37(11): 227-240, 14.]
[5] CHANDRA P, FISHER M L. Coordination of production and distribution planning[J]. European Journal of Operational Research, 1994, 72(3): 503-517.
[6] YAGMUR E, KESEN S E. A memetic algorithm for joint production and distribution scheduling with due dates[J]. Computers & Industrial Engineering, 2020, 142(2): 106342.
[7] MOONS S, RAMAEKERS K, CARIS A, et al. Integrating production scheduling and vehicle routing decisions at the operational decision level: A review and discussion[J]. Computers & Industrial Engineering, 2017, 104: 224-245.
[8] LIU L, LIU S. Integrated production and distribution problem of perishable products with a minimum total order weighted delivery time[J]. Mathematics, 2020, 8(2): 146.
[9] KERGOSIEN Y, GENDREAU M, BILLAUT J C. A Benders decomposition-based heuristic for a production and outbound distribution scheduling problem with strict delivery constraints[J]. European Journal of Operational Research, 2017, 262(1): 287-298.
[10] KESEN S E, BEKTAS T. Integrated production scheduling and distribution planning with time windows[M]//Lean and Green Supply Chain Management. Berlin, Germany: Springer-Verlag, 2019: 231-252.
[11] MOHAMMADI S, CHERAGHALIKHANI A, RAMEZANIAN R. A joint scheduling of production and distribution operations in a flow shop manufacturing system[J]. Scientia Iranica, 2018, 25(2): 911-930.
[12] FU Y P, HOU Y S, CHEN Z H, et al. Modelling and scheduling integration of distributed production and distribution problems via black widow optimization[J]. Swarm and Evolutionary Computation, 2022, 68: 101015.
[13] FU L L, ALOULOU M A, TRIKI C. Integrated production scheduling and vehicle routing problem with job splitting and delivery time windows[J]. International Journal of Production Research, 2017, 55(20): 5942-5957.
[14] MOHAMMADI S, AL-E-HASHEM S M J M, REKIK Y. An integrated production scheduling and delivery route planning with multi-purpose machines: A case study from a furniture manufacturing company[J]. International Journal of Production Economics, 2020, 219: 347-359.
[15] YAGMUR E, KESEN S E. Multi-trip heterogeneous vehicle routing problem coordinated with production scheduling: Memetic algorithm and simulated annealing approaches[J]. Computers & Industrial Engineering, 2021, 161: 107649.
[16] FAN Houming, YANG Cheng, ZHANG Yueguang, et al. Multi-depot mixed fleet vehicle routing problem with mixed time windows[J]. Computer Integrated Manufacturing Systems, 2023, 29(10): 3529-3546 (in Chinese). [范厚明，杨成，张跃光，等. 混合时间窗下多中心混合车队车辆路径优化研究[J]. 计算机集成制造系统，2023, 29(10): 3529-3546.]
[17] SHI Y H. Brain storm optimization algorithm[C]//Proceedings of International Conference in Swarm Intelligence. Berlin, Germany: Springer-Verlag, 2011: 303-309.
[18] KE L J. A brain storm optimization approach for the cumulative capacitated vehicle routing problem[J]. Memetic Computing, 2018, 10(4): 411-421.
[19] HOU Y S, FU Y P, GAO K Z, et al. Modelling and optimization of integrated distributed flow shop scheduling and distribution problems with time windows[J]. Expert Systems with Applications, 2022, 187: 115827.
[20] GUO X W, LIU S X, ZHOU M C, et al. Dual-objective program and scatter search for the optimization of disassembly sequences subject to multi resource constraints[J]. IEEE Transactions on Automation Science and Engineering, 2018, 15(3): 1091-1103.
[21] HAO J H, LI J Q, DU Y, et al. Solving distributed hybrid flowshop scheduling problems by a hybrid brain storm optimization algorithm[J]. IEEE Access, 2019, 7: 66879-66894.
[22] DEB K, PRATAP A, AGARWAL S, et al. A fast and elitist multiobjective genetic algorithm: NSGA-II[J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2): 182-197.
[23] WANG G C, GAO L, LI X Y, et al. Energy-efficient distributed permutation flow shop scheduling problem using a multi-objective whale swarm algorithm[J]. Swarm and Evolutionary Computation, 2020, 57: 100716.
[24] FERNANDEZ-VIAGAS V, RUIZ R, FRAMINAN J M. A new vision of approximate methods for the permutation flowshop to minimise makespan: State-of-the-art and computational evaluation[J]. European Journal of Operational Research, 2017, 257(3): 707-721.
[25] HOMBERGER J, GEHRING H. A two-phase hybrid meta-heuristic for the vehicle routing problem with time windows[J]. European Journal of Operational Research, 2005, 162(1): 220-238.
[26] ZHANG Q F, LI H. A multiobjective evolutionary algorithm based on decomposition[J]. IEEE Transactions on Evolutionary Computation, 2007, 11(6): 712-731.
[27] PAN Y X, GAO K Z, LI Z W, et al. Solving biobjective distributed flow-shop scheduling problems with lot-streaming using an improved Jaya algorithm[J]. IEEE Transactions on Cybernetics, 2023, 53(6): 3818-3828.
[28] MONTGOMERY D C. Design and analysis of experiments[M]. Hoboken, N. J., USA: John Wiley & Sons, 2017.
[29] FU Y P, ZHOU M C, GUO X W, et al. Scheduling dual-objective stochastic hybrid flow shop with deteriorating jobs via bi-population evolutionary algorithm[J]. IEEE Transactions on Systems, Man, and Cybernetics: Systems, 2020, 50(12): 5037-5048.
[30] FU Y P, LI H B, HUANG M, et al. Bi-objective modelling and optimization for stochastic two-stage open shop scheduling problems in the sharing economy[J]. IEEE Transactions on Engineering Management, 2023, 70(10): 3395-3409.
[31] DERRAC J, GARCIA S, MOLINA D, et al. A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms[J]. Swarm and Evolutionary Computation, 2011, 1(1): 3-18.
