FJSP：蜣螂优化算法( Dung beetle optimizer, DBO)求解柔性作业车间调度问题（FJSP），提供MATLAB代码

本文主要是介绍FJSP：蜣螂优化算法( Dung beetle optimizer, DBO)求解柔性作业车间调度问题（FJSP），提供MATLAB代码，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

一、柔性作业车间调度问题

柔性作业车间调度问题（Flexible Job Shop Scheduling Problem，FJSP），是一种经典的组合优化问题。在FJSP问题中，有多个作业需要在多个机器上进行加工，每个作业由一系列工序组成，每个工序需要在特定的机器上完成。同时，每个机器一次只能处理一个工序，且每个工序的处理时间可能不同。

FJSP问题的目标是找到一个最优的作业调度方案，使得所有作业的完成时间最小化。这个问题的难点在于需要考虑到多个作业、多个机器和多个工序之间的复杂关系，并且需要在有限的时间内找到最优解。

解决FJSP问题的方法包括启发式算法、精确算法和元启发式算法等。启发式算法通过一系列规则和策略来生成调度方案，常见的方法有遗传算法、模拟退火算法和禁忌搜索算法等。精确算法则通过穷举搜索或者动态规划等方法来找到最优解，但在实际应用中可能面临计算复杂度过高的问题。元启发式算法则结合了多种启发式算法和精确算法的优点，通过组合不同的方法来求解FJSP问题。
FJSP问题的难点主要体现在以下几个方面：

组合爆炸：FJSP问题中，每个工件都有多个工序需要完成，而每个工序都有多个可选的机器可以执行。这导致了组爆炸的问题，可能的调度方案数量非常庞，难以穷举所有可能性。
优化目标多样：FJSP问题通常有多个优化目标，如最小化总加权完成时间、最小化总延迟时间等。这些目标之间可能存在冲突，使找到一个全局最优解变得困难。
资源约束：FJSP问题中，每个机器在同一时间只能执行一个工序，且每个工序需要一定的时间和资源。这些资源约束增加了问题的复杂性，需要在满足约束条件的前提下进行调度。
实时性要求：在实际生产中，FJSP问题通常需要考虑实时性要求，即要求在有限的时间内生成一个可行的调度方案。这增加了问题的难度，需要在有限时间内找到一个较优的解。

柔性作业车间调度问题( FJSP) 的描述如下：n个工件 ${J,J_2,..,J_n\}$ 要在 $m$ 台机器 ${M_1,M_2,..,M_m\}$ 上加工。每个工件包含一道或多道工序，工序顺序是预先确定的，每道工序可以在多台不同加工机器上进行加工，工序的加工时间随加工机器的不同而不同。调度目标是为每道工序选择最合适的机器、确定每台机器上各个工序的最佳加工顺序以及开工时间，使整个系统的某些性能指标达到最优。因此，柔性作业车间调度问题包含两个子问题：确定各工件的加工机器 (机器选择子问题) 和确定各个机器上的加工先后顺序 (工序排序子问题)。

此外，在加工过程中还需要满足下面的约束条件：
(1) 同一台机器同一时刻只能加工一个工件；
(2) 同一工件的同一道工序在同一时刻只能被一台机器加工；
(3) 每个工件的每道工序一旦开始加工不能中断；
(4) 不同工件之间具有相同的优先级；
(5)不同工件的工序之间没有先后约束，同一工件的工序之间有先后约束；
(6)所有工件在零时刻都可以被加工。

1.1符号描述

$n :$ 工件总数；
$m :$ 机器总数；
$i, e :$ 机器序号， $i, e = 1, 2, 3, ..., m$ ;
$j, k :$ 工件序号， $j, k = 1, 2, 3, ..., n;$ $h_j:$ 工件 $j$ 的工序总数；
$h, l :$ 工序序号， $h=1,2,3,...,h_j$ ;
$\Omega_{jh}:$ 工件 $j$ 的第 $h$ 道工序的可选加工机器集；
$m_{jh}:$ 工件 $j$ 的第 $h$ 道工序的可选加工机器数；
$O_{jh}:$ 工件 $j$ 的第 $h$ 道工序；
$M_{ijh}:$ 工件 $j$ 的第 $h$ 道工序在机器 $i$ 上加工；
$p_{ijh}:$ 工件 $j$ 的第 $h$ 道工序在机器 $i$ 上的加工时间；
$s_{jh}:$ 工件 $j$ 的第 $h$ 道工序加工开始时间；
$c_{jh}:$ 工件 $j$ 的第 $h$ 道工序加工完成时间；
$d_j:$ 工件 $j$ 的交货期；
$L$ : 一个足够大的正数；
$C_j$ : 每个工件的完成时间；
$C_{\max}:$ 最大完工时间；
$T_o:\quad T_o=\sum_{j=1}^nh_j$ , 所有工件工序总数；
$x_{ijh}=\begin{cases}1,\text{如果工序}O_{jh}\text{选择机器}i;\\0,\text{否则；}\end{cases}$
$y_{ijhkl}=\begin{cases}1,\text{如果}O_{ijh}\text{先于}O_{ikl}\text{加工};\\0,\text{否则};\end{cases}$

1.2约束条件

$C_{1}:s_{jh}+x_{ijh}\times p_{ijh}\leq c_{jh}$

其中： $i=1,\ldots,m;j=1,\ldots,n;$ $h=1,\ldots,h_j$
$C_{2}:c_{jh}\leq s_{j(h+1)}$
其中 $:j=1,\ldots,n;h=1,...,h_j-1$
$C_{3}:c_{jh_j}\leq C_{\max}$
其中： $j = 1, ..., n$
$C_{4}:s_{jh}+p_{ijh}\leq s_{kl}+L(1-y_{ijhkl})$

其中 $:j=0,\ldots,n;k=1,\ldots,n;h=1,\ldots,h_j;l=1,\ldots,h_k;i=1,\ldots,m$
$C_{5}:c_{jh}\leq s_{j(h+1)}+L(1-y_{iklj(h+1)})$

其中 $:j=1,\ldots,n;k=0,\ldots,n;h=1,\ldots,h_j-1;\quad l=1,\ldots,h_k;\quad i=1,\ldots,m$
$h_{1}:\sum_{i=1}^{m_{jh}}x_{ijh}=1$
其中： $h=1,...,h_j;j=1,...,n;$

$h_{2}:\sum_{j=1}^n\sum_{h=1}^{h_j}y_{ijhkl}=x_{ikl}$

其中： $i=1,\ldots,m;k=1,\ldots,n;l=1,\ldots,h_k$
$h_{3}:\sum_{i=1}^n\sum_{i=1}^{n_k}y_{ijhkl}=x_{ijh}$

其中： $i=1,\ldots,m;j=1,\ldots,n;\quad h=1,\ldots,h_k$
$C_{6}:s_{jh}\geq0,c_{jh}\geq0$

其中 $j=0,1,...,n;h=1,...,h_j$

$C_{1}$ 和 $C_{2}$ 表示每一个工件的工序先后顺序约束 ;
$C_{3}$ 表示工件的完工时间的约束,即每一个工件的完工时间不可能超过总的完工时间 ;
$C_{4}$ 和 $C_{5}$ 表示同一时刻同一台机器只能加工一道工序 ;
$h_{1}$ 表示机器约束,即同一时刻同一道工序只能且仅能被一台机器加工;
$h_{2}$ 和 $h_{3}$ 表示存在每一台机器上可以存在循环操作 ;
$C_{6}$ 表示各个参数变量必须是正数。

1.3目标函数

FJSP的目标函数是最大完工时间最小。完工时间是每个工件最后一道工序完成的时间，其中最大的那个时间就是最大完工时间(makespan)。它是衡量调度方案的最根本指标，主要体现车间的生产效率，如下式所示：

$f=\min(\max_{\mathrm{l\leq}j\leq n}(C_j))$

参考文献：
[1]张国辉.柔性作业车间调度方法研究[D].华中科技大学,2009.

二、算法简介

蜣螂优化算法( Dung beetle optimizer, DBO)是由 Jiankai Xue 等于2022 年提出的一种群体智能优化算法。其灵感来源于蜣螂的生物行为过程，具有寻优能力强，收敛速度快的特点。
参考文献：Xue, J., Shen, B. Dung beetle optimizer: a new meta-heuristic algorithm for global optimization. J Supercomput (2022). https://doi.org/10.1007/s11227-022-04959-6

原文链接：https://blog.csdn.net/weixin_46204734/article/details/128138381

三、算法求解FJSP

3.1部分代码

dim=2*sum(operaNumVec);
LB = -jobNum * ones(1, dim);
UB = jobNum * ones(1, dim);
Max_iteration = 100;
SearchAgents_no = 100;
fobj=@(x)fitness(x, MachineNum,jobNum,jobInfo,operaNumVec,candidateMachine);%% 优化算法求解FJSP
[fMin , bestX, Convergence_curve ] = DBO(SearchAgents_no,Max_iteration,LB,UB,dim,fobj);
machineTable=GetMachineTable(bestX, MachineNum,jobNum,jobInfo,operaNumVec,candidateMachine);%% 画收敛曲线图
figure
plot(Convergence_curve,'g-','linewidth',2)
xlabel('迭代次数')
ylabel('最大完工时间')
legend('DBO')
saveas(gca,'1.jpg');