【车间调度】基于蜣螂优化算法DBO求解零等待流水车间调度问题NWFSP附Matlab代码

✅作者简介：热爱科研的Matlab仿真开发者，擅长毕业设计辅导、数学建模、数据处理、算法改进、程序设计科研仿真。

🍎完整代码获取 定制创新 论文复现私信

🍊个人信条：做科研，博学之、审问之、慎思之、明辨之、笃行之，是为：博学慎思，明辨笃行。

🔥 内容介绍

在现代制造业中，车间调度问题对于提高生产效率、降低成本至关重要。零等待流水车间调度问题（NWFSP）作为其中的一个重要分支，要求工件在各加工阶段之间无等待时间，这对生产流程的紧凑性和协调性提出了更高要求。蜣螂优化算法（DBO）作为一种新兴的智能优化算法，为解决 NWFSP 提供了新的思路和方法。

零等待流水车间调度问题（NWFSP）

问题描述

假设有 n 个工件需要在 m 台机器上按照相同的加工顺序依次加工。每个工件 i 在机器 j 上的加工时间为 pij。在 NWFSP 中，一旦一个工件在某台机器上开始加工，该工件在后续机器上的加工必须立即开始，即工件在机器间转移时不允许有等待时间。目标是确定工件在各机器上的加工顺序，使得完成所有工件加工的总时间（即最大完工时间 Cmax）最小。

蜣螂优化算法（DBO）

算法灵感来源

蜣螂优化算法受到蜣螂滚动粪球这一自然行为的启发。蜣螂在寻找合适的地点掩埋粪球时，会根据自身经验和周围环境信息不断调整滚动方向和力度。在优化问题中，将每个解看作一个蜣螂，解的质量对应粪球的 “吸引力”，算法通过模拟蜣螂的滚动行为来寻找最优解。

算法基本原理

1. 初始化

   ：随机生成一定数量的蜣螂（即初始解），每个蜣螂的位置代表一种工件加工顺序。同时，根据目标函数计算每个蜣螂的适应度值，适应度值反映了该解对应的最大完工时间 Cmax，值越小表示解越优。

2. 滚动阶段

   ：每个蜣螂根据自身位置和其他蜣螂的位置信息，按照一定规则调整自己的位置，模拟蜣螂滚动粪球的过程。例如，蜣螂可能向适应度值更好的蜣螂靠近，同时也会加入一定的随机因素以避免陷入局部最优。这个过程通过位置更新公式来实现，位置更新公式通常涉及当前位置、目标位置（如适应度最优的蜣螂位置）以及一些控制参数。

3. 挖掘阶段

   ：在滚动一定次数后，部分蜣螂进入挖掘阶段。在挖掘阶段，蜣螂对当前位置进行局部搜索，尝试找到更好的解。这有助于在局部范围内进一步优化解的质量。挖掘阶段可以通过对当前解进行微小扰动，并重新计算适应度值来实现。

4. 更新与迭代

   ：根据滚动和挖掘阶段的结果，更新蜣螂的位置和适应度值。重复上述步骤，直到满足终止条件，如达到最大迭代次数或适应度值收敛。

5. 输出最优解

   ：算法终止时，输出适应度值最优的蜣螂位置，即对应最优的工件加工顺序。

基于 DBO 求解 NWFSP 的实现

编码与解码

1. 编码

   ：采用排列编码方式，将 n 个工件的加工顺序直接编码为一个长度为 n 的排列。例如，排列 [3,1,2] 表示第 3 个工件先加工，然后是第 1 个工件，最后是第 2 个工件。

2. 解码

   ：根据编码得到的工件加工顺序，结合各工件在不同机器上的加工时间，计算每个工件在各机器上的开始时间和完工时间，进而得出最大完工时间 Cmax，作为该编码对应的适应度值。

适应度函数

适应度函数直接采用 NWFSP 的目标函数，即最大完工时间 Cmax。对于给定的工件加工顺序编码，通过解码过程计算出 Cmax，Cmax 值越小，适应度越高。

⛳️ 运行结果

📣 部分代码

% -----------------------------------------------------------------------------------------------------------

% Dung Beetle Optimizer: (DBO) (demo)

% Programmed by Jian-kai Xue

% Updated 28 Nov. 2022.

%

% This is a simple demo version only implemented the basic

% idea of the DBO for solving the unconstrained problem.

% The details about DBO are illustratred in the following paper.

% (To cite this article):

% Jiankai Xue & Bo Shen (2022) Dung beetle optimizer: a new meta-heuristic

% algorithm for global optimization. The Journal of Supercomputing, DOI:

% 10.1007/s11227-022-04959-6

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [fMin , bestX, Convergence_curve ] = DBO(pop, M,c,d,dim,fobj )

P_percent = 0.2; % The population size of producers accounts for "P_percent" percent of the total population size

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pNum = round( pop * P_percent ); % The population size of the producers

lb= c.*ones( 1,dim ); % Lower limit/bounds/ a vector

ub= d.*ones( 1,dim ); % Upper limit/bounds/ a vector

%Initialization

for i = 1 : pop

x( i, : ) = lb + (ub - lb) .* rand( 1, dim );

fit( i ) = fobj( x( i, : ) ) ;

end

pFit = fit;

pX = x;

XX=pX;

[ fMin, bestI ] = min( fit ); % fMin denotes the global optimum fitness value

bestX = x( bestI, : ); % bestX denotes the global optimum position corresponding to fMin

% Start updating the solutions.

for t = 1 : M

[fmax,B]=max(fit);

worse= x(B,:);

r2=rand(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1 : pNum

if(r2<0.9)

r1=rand(1);

a=rand(1,1);

if (a>0.1)

a=1;

else

a=-1;

end

x( i , : ) = pX( i , :)+0.3*abs(pX(i , : )-worse)+a*0.1*(XX( i , :)); % Equation (1)

else

aaa= randperm(180,1);

if ( aaa==0 ||aaa==90 ||aaa==180 )

x( i , : ) = pX( i , :);

end

theta= aaa*pi/180;

x( i , : ) = pX( i , :)+tan(theta).*abs(pX(i , : )-XX( i , :)); % Equation (2)

end

x( i , : ) = Bounds( x(i , : ), lb, ub );

fit( i ) = fobj( x(i , : ) );

end

[ fMMin, bestII ] = min( fit ); % fMin denotes the current optimum fitness value

bestXX = x( bestII, : ); % bestXX denotes the current optimum position

R=1-t/M; %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Xnew1 = bestXX.*(1-R);

Xnew2 =bestXX.*(1+R); %%% Equation (3)

Xnew1= Bounds( Xnew1, lb, ub );

Xnew2 = Bounds( Xnew2, lb, ub );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Xnew11 = bestX.*(1-R);

Xnew22 =bestX.*(1+R); %%% Equation (5)

Xnew11= Bounds( Xnew11, lb, ub );

Xnew22 = Bounds( Xnew22, lb, ub );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = ( pNum + 1 ) :12 % Equation (4)

x( i, : )=bestXX+((rand(1,dim)).*(pX( i , : )-Xnew1)+(rand(1,dim)).*(pX( i , : )-Xnew2));

x(i, : ) = Bounds( x(i, : ), Xnew1, Xnew2 );

fit(i ) = fobj( x(i,:) ) ;

end

for i = 13: 19 % Equation (6)

x( i, : )=pX( i , : )+((randn(1)).*(pX( i , : )-Xnew11)+((rand(1,dim)).*(pX( i , : )-Xnew22)));

x(i, : ) = Bounds( x(i, : ),lb, ub);

fit(i ) = fobj( x(i,:) ) ;

end

for j = 20 : pop % Equation (7)

x( j,: )=bestX+randn(1,dim).*((abs(( pX(j,: )-bestXX)))+(abs(( pX(j,: )-bestX))))./2;

x(j, : ) = Bounds( x(j, : ), lb, ub );

fit(j ) = fobj( x(j,:) ) ;

end

% Update the individual's best fitness vlaue and the global best fitness value

XX=pX;

for i = 1 : pop

if ( fit( i ) < pFit( i ) )

pFit( i ) = fit( i );

pX( i, : ) = x( i, : );

end

if( pFit( i ) < fMin )

% fMin= pFit( i );

fMin= pFit( i );

bestX = pX( i, : );

% a(i)=fMin;

end

end

Convergence_curve(t)=fMin;

end

% Application of simple limits/bounds

function s = Bounds( s, Lb, Ub)

% Apply the lower bound vector

temp = s;

I = temp < Lb;

temp(I) = Lb(I);

% Apply the upper bound vector

J = temp > Ub;

temp(J) = Ub(J);

% Update this new move

s = temp;

function S = Boundss( SS, LLb, UUb)

% Apply the lower bound vector

temp = SS;

I = temp < LLb;

temp(I) = LLb(I);

% Apply the upper bound vector

J = temp > UUb;

temp(J) = UUb(J);

% Update this new move

S = temp;

%---------------------------------------------------------------------------------------------------------------------------

🔗 参考文献

[1]刘柄廷.共生生物算法及其在零等待流水车间调度问题中的研究[D].兰州理工大学,2023.

🍅更多免费数学建模和仿真教程关注领取