基于灰狼优化算法求解柔性作业车间问题matlab

1 简介

柔性车间作业调度（Flexible Jobshop Scheduling Problem，简称FJSP），属于一个典型的NP难问题。解决作业车间调度问题的方法主要有以下四类：1）运筹学方法，包括整数规划，分支定界方法等；2）启发式规则；3）神经网络方法；4）元启发式算法。有很多知名的智能优化算法都已经成功应用FJSP领域，如蚁群算法、遗传算法、模拟退火算法和粒子群算法等。由NFL[1]定理可知，这个定理在逻辑上证明了没有哪种元启发式最适合解决所有的优化问题。这个领域每年都会提出新的元启发式算法。

GWO是一种基于群体的元启发式算法，它设定每个狼群都有层次结构，狼群中的每只狼都有特定的角色。Mirjalili在文献[2]中，通过对29个连续函数优化问题的测试，结果表明GWO在求解精度和稳定性上要明显优于PSO、DE和GSA等算法。GWO中所有候选狼（解决方案）都被最好的三只狼所吸引，从而更快地聚集到这些好狼位置周围。灰狼算法多用于解决连续函数问题的研究，在FJSP这类非常复杂的离散组合优化问题应用还比较少，本文针对经典的FJSP的特点，通过优化初始种群，采用新的解码机制，证明灰狼算法可以求解FJSP问题。

2 部分代码

%该程序用于解决柔性作业车间调度，m个工件，n道工序，其中n为最大工序数，工件的工序
%数可以少于n，加工机器数为M，每个工件的每道工序具有多个机器可以选择，对应的时间
%不同，其中初始种群的储存方式采用cell数据类型
%Version:1.3
%fileDescription:调度机器可选的柔性作业车间问题，甘特图已完善,GWO,8*8实例
%last edit time:2019-6-7
function GWO_Model_FJSP_1_3_8_8()
count = 600;     %迭代次数
N = 50;          %种群规模
m = 8;             %工件数
n = 4;             %工序数
M = 8;             %机器数
a =2;              %计算A/C协同系数的
plotif = 1;        %控制程序是否进行绘图
s = input(m,n);    %数据输入
[p,TN] = initial_p(m,n,N,s,M);    %生成初始种群50,采用细胞结构，每个元素为8*4
P = machine(n,M);
FIT = zeros(count,1);
aveFIT = zeros(count,1);
X1=randperm(count);       %收敛图形的横坐标X
X=sort(X1);
%------------------------输出最优解的时有用------------------------------
best_fit = 1000;            %改变模型需要修改此参数
best_p = zeros(m,n);
best_TN = zeros(m,n);
Y1p = zeros(m,1);
Y2p = zeros(m,1);
Y3p = zeros(m,1);
minfit3 =  1000000000;
%-------------------------进行迭代--------------------------------------
for i = 1:count
  [fit,Y1,Y2,Y3] = object(p,TN,N,P,m,n);  
  [newp,newTN] = GWO(fit,p,TN,N,m,n,s,a);
   a = a-2/(count-1);        %a的值会线性下降
   if best_fit > min(fit)
      [best_p,best_TN,best_fit,Y1p,Y2p,Y3p]=best(best_fit,best_p,fit,best_TN,Y1p,Y2p,Y3p,p,TN,Y1,Y2,Y3);
   end
   p = newp;
   TN = newTN;
   minfit = min(fit);
   if minfit3>minfit
       minfit3 = minfit;
   end
   FIT(i) = minfit3;    %用于适应度函数的
   aveFIT(i) = mean(fit);      %用于适应度函数的
end
%------------------投射最佳方案数据--------------------------------------
 
   fprintf('最优解：%d\n',best_fit);
   fprintf('工序1 工序2 工序3 工序4\n');
   best_p
   fprintf('时间1 时间2 时间3 时间4\n');
   best_TN
%------------------------收敛曲线----------------------------------------
   if plotif == 1
   figure;
   plot(X,FIT,'r');
   hold on;
   plot(X,aveFIT,'b');
   title('收敛曲线');
   hold on;
   legend('最优解','平均值');
%-------------------------甘特图-----------------------------------------
figure;
w=0.5;       %横条宽度
set(gcf,'color','w');      %图的背景设为白色
for i = 1:m
   for j = 1:n
       color=[1,0.98,0.98;1,0.89,0.71;0.86,0.86,0.86;0.38,0.72,1;1,0,1;0,1,1;0,1,0.49;1,0.87,0.67;0.39,0.58,0.92;0.56,0.73,0.56];
       a = [Y1p(i,j),Y2p(i,j)];
       x=a(1,[1 1 2 2]);      %设置小图框四个点的x坐标
       y=Y3p(i,j)+[-w/2 w/2 w/2 -w/2];   %设置小图框四个点的y坐标
       color = [color(i,1),color(i,2),color(i,3)];
       p=patch('xdata',x,'ydata',y,'facecolor',color,'edgecolor','k');    %facecolor为填充颜色，edgecolor为图框颜色
           text(a(1,1)+0.5,Y3p(i,j),[num2str(i),'-',num2str(j)]);    %显示小图框里的数字位置和数值
   end
end
xlabel('加工时间/s');      %横坐标名称
ylabel('机器');            %纵坐标名称
title({[num2str(m),'*',num2str(M),'的一个最佳调度（最短完工时间为',num2str(best_fit),')']});      %图形名称
axis([0,best_fit+2,0,M+1]);         %x轴，y轴的范围
set(gca,'Box','on');       %显示图形边框
set(gca,'YTick',0:M+1);     %y轴的增长幅度
set(gca,'YTickLabel',{'';num2str((1:M)','M%d');''});  %显示机器号
hold on;
   end
%--------------------------输入数据---------------------------------
function s = input(m,n)      %输入数据
s = cell(m,n);
s{1,1}=[1 2 3 4 5 7 8;5 3 5 3 3 10 9];
s{1,2}=[1 3 4 5 6 7 8;10 5 8 3 9 9 6];
s{1,3}=[2 4 5 6 7 8;10 5 6 2 4 5];
s{1,4}=[0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0];
s{2,1}=[1 2 3 4 5 7;5 7 3 9 8 9];
s{2,2}=[2 3 4 5 6 7 8;8 5 2 6 7 10 9];
s{2,3}=[2 4 5 6 7 8;10 5 6 4 1 7];
s{2,4}=[1 2 3 4 5 6;10 8 9 6 4 7];

s{3,1}=[1 4 5 6 7 8;10 7 6 5 2 4];
s{3,2}=[2 3 4 5 6 7;10 6 4 8 9 10];
s{3,3}=[1 2 3 4 6 8;1 4 5 6 10 7];
s{3,4}=[0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0];
s{4,1}=[1 2 3 4 5 6 7 8;3 1 6 5 9 7 8 4];
s{4,2}=[1 2 3 4 5 6 7 8;12 11 7 8 10 5 6 9];
s{4,3}=[1 2 3 4 5 6 7 8;4 6 2 10 3 9 5 7];
s{4,4}=[0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0];
s{5,1}=[1 2 3 4 5 7;3 6 7 8 9 10];
s{5,2}=[1 3 4 5 6 7;10 7 4 9 8 6];
s{5,3}=[2 3 4 5 6 7;9 8 7 4 2 7];
s{5,4}=[1 2 4 5 6 7 8;11 9 6 7 5 3 6];

s{6,1}=[1 2 3 4 5 6 8;6 7 1 4 6 9 10];
s{6,2}=[1 3 4 5 6 7 8;11 9 9 9 7 6 4];
s{6,3}=[1 2 3 4 5 7;10 5 9 10 11 10];
s{6,4}=[0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0];
s{7,1}=[1 2 3 4 5 7;5 4 2 6 7 10];
s{7,2}=[2 4 5 6 7 8;9 9 11 9 10 5];
s{7,3}=[2 3 4 5 6 8;8 9 3 8 6 10];
s{7,4}=[0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0];
s{8,1}=[1 2 3 4 6 8;2 8 5 9 4 10];
s{8,2}=[1 2 3 4 5 7;7 4 7 8 9 10];
s{8,3}=[1 2 4 5 6 7 8;9 9 8 5 6 7 1];
s{8,4}=[1 3 4 5 6 7;9 3 7 1 5 8];  

3 仿真结果