【车间调度】基于matlab差分进化算法求解作业车间调度问题【含Matlab源码 1743期】

一、车间调度简介

1 车间调度定义
车间调度是指根据产品制造的合理需求分配加工车间顺序，从而达到合理利用产品制造资源、提高企业经济效益的目的。车间调度问题从数学上可以描述为有n个待加工的零件要在m台机器上加工。问题需要满足的条件包括每个零件的各道工序使用每台机器不多于1次，每个零件都按照一定的顺序进行加工。

2 传统作业车间调度
传统作业车间带调度实例

有若干工件，每个工件有若干工序，有多个加工机器，但是每道工序只能在一台机器上加工。对应到上面表格中的实例就是，两个工件，工件J1有三道工序，工序Q11只能在M3上加工，加工时间是5小时。
约束是对于一个工件来说，工序的相对顺序不能变。O11->O12->O13。每时刻，每个工件只能在一台机器上加工；每个机器上只能有一个工件。
调度的任务则是安排出工序的加工顺序，加工顺序确定了，因为每道工序只有一台机器可用，加工的机器也就确定了。
调度的目的是总的完工时间最短（也可以是其他目标）。举个例子，比如确定了O21->O22->O11->O23->O12->O13的加工顺序之后，我们就可以根据加工机器的约束，计算出总的加工时间。
M2加工O21消耗6小时，工件J2当前加工时间6小时。
M1加工O22消耗9小时，工件J2当前加工时间6+9=15小时。
M3加工O11消耗5小时，工件J1当前加工时间5小时。
M4加工O23消耗7小时，工件J2加工时间15+7=22小时。
M1加工O12消耗11小时，但是要等M1加工完O22之后才开始加工O12，所以工件J1的当前加工时间为max(5,9)+11=20小时。
M5加工O13消耗8小时，工件J2加工时间20+8=28小时。
总的完工时间就是max(22,28)=28小时。

2 柔性作业车间调度
柔性作业车间带调度实例（参考自高亮老师论文
《改进遗传算法求解柔性作业车间调度问题》——机械工程学报）

相比于传统作业车间调度，柔性作业车间调度放宽了对加工机器的约束，更符合现实生产情况，每个工序可选加工机器变成了多个，可以由多个加工机器中的一个加工。比如上表中的实例，J1的O12工序可以选择M2和M4加工，加工时间分别是8小时和4小时，但是并不一定选择M4加工，最后得出来的总的完工时间就更短，所以，需要调度算法求解优化。

相比于传统作业车间，柔性车间作业调度的调度任务不仅要确定工序的加工顺序，而且需要确定每道工序的机器分配。比如，确定了O21->O22->O11->O23->O12->O13的加工顺序，我们并不能相应工序的加工机器，所以还应该确定对应的[M1、M3、M5]->[M1、M2、M3]->[M1、M2、M3、M4、M5]->[M2、M3、M4、M5]->[M2、M4]->[M1、M3、M4、M5]的机器组合。调度的目的还是总的完工时间最短（也可以是其他目标，比如机器最大负荷最短、总的机器负荷最短）

二、差分进化算法简介

如同所有的优化算法一样, 差分进化算法基于种群的进化算法。差分进化算法主要的参数主要有种群规模NP, 解空间的维数D, 缩放因子F和交叉概率Cr。D维矢量XGi=[xGi,1, xGi,2, …, xi, DG], i=1, 2, …, Np, 表示G代第i个个体。变异和交叉操作在每一代中产生新的种群, 然后选择操作实现种群优化。基本差分进化算法过程如下。

1 变异
通过变异操作, 新产生的个体由初始种群中任意选取的个体加上任意选取的两个个体的差分量乘以缩放因子。根据如下所示, 缩放因子我们采用自适应缩放因子, 如公式 (1) 所示。

VG+1i=xGr1+F (xGr2-xGr3) (1)

r1, r2, r3∈{1, 2, …, Np}任意选择且互不相等, 并且和i也不相等。缩放因子F控制着差分量 (xr2G-xr3G) 的大小。

2 交叉
交叉操作用来增加种群多样性, UG+1i,j表示第i个个体交叉个体第j个分量, rand ( ) 表示在[0, 1]区间任取的随机数, Cr是交叉概率, 而条件j=rand (0, D) 是为了确保试验个体中至少有一位是由VG+1i贡献的, 增加试验个体的多样性。交叉时我们采用分组交叉思想, 增强算法的局部搜索能力。

UG+1i,j={VG+1i,[rand()≤CR或j=rand(0,D)]xGi,j,其他 (2)

3 选择
选择操作是为了决定交叉个体是否保留到G+1代。对于最小化问题, 交叉个体UG+1i,j和原始个体通过贪婪策略进行比较, 目标函数更小的个体被保留作为新产生的种群个体。选择操作公式如下所示

这里的f (x) 代表目标函数, XiG+1代表新产生的种群个体。

文中算法改进主要有2方面。算法的参数值在算法开始时已经被设定, 随着迭代次数的改变, 可能不适应种群变化的要求。差分算法的控制参数主要有缩放因子 (F) 和交叉概率 (Cr) 。缩放因子 (F) 的大小决定了算法的收敛速度和种群的多样性, 在进化过程中发挥着重要的作用。选择合适的F能够平衡算法全局搜索和局部搜索。较小的F有助于提高算法的局部搜索能力, 但同时也增加了算法陷入局部最优, 出现早熟的可能性;而较大的F有利于增加种群的多样性, 提高全局寻优能力, 但会降低算法的收敛速度。我们提出了一种自适应的缩放因F, 使缩放因子能够随着迭代次数的改变, 适应当前种群的变化情况, 使得最优解可以保存。颜学峰, 余娟等也提出了自适应差分进化算法[7]748, 主要改进了算法的全局寻优能力。本文提出的自适应的缩放因子可以从以下公式得出
F=0.8[ (Gmax-G+1) /Gmax] (4)
Gmax是种群最大迭代次数, G是种群当前迭代次数。同时为了能够更早的找到最优解, 减少算法的计算量。在交叉过程中, 我们采用种群分组交叉思想。经种群变异产生的子代个体, 通过比较其适应度大小, 适应度更好的个体直接进入选择操作, 而剩下的个体根据交叉概率选择性的进入选择操作。通过分组交叉, 可以大大提高种群搜索到最优解的概率, 同时增加算法的健壮性。

三、Job Shop问题描述

Job Shop研究n个工件在m台机器上加工, 每个工件加工工序不同, 每个工序需要由特定的机器经过一定的加工时间才能完成。已知每个工件每道工序的加工时间, 每道工序加工机器号, 使得某些加工性能指标达到最优。如果每个工件在每台机器上加工顺序相同, Job Shop转化为Flow Shop调度问题, 如果每台机器上每个工件的加工次序也相同, 那么Flow Shop转化为置换Flow Shop (PFSP) 调度问题。

1 问题编码
差分进化算法是一种在连续空间的随机搜索算法, 一般采用的是实数编码方式。而作业车间研究的是确定每个工件在每台机器上的加工顺序, 使得莫项指标达到最优。因此算法解空间内的个体无法直接代替工序进行编码。首先我们要将解空间内的实数转化为可以表示工序的整数。对于n个工件, m台机器上的JSP问题, 为了便于理解, 我们这里先讨论3个工件在2台机器上的编码。

第1步:初始化种群个体时, 因为要确定每个工件在每台机器上的加工工序, 所以解空间的维数D=n×m。Xi= (xi, 1, xi, 2, …, xi, D) 代表初始个体。
第2步:初始个体Xi= (xi, 1, xi, 2, …, xi, D) 按照降序进行排列, 中间排序结果用ϕi, j表示。
第3步:降序排列的中间结果ϕi, j按照公式 (5) 转化为工序。
πi, j=ϕi, j%n+1, j=1, 2, …, D;D=n×m (5)
初始个体Xi= (3.2, 1.2, 6.5, 5.8, 0.6, 9.1) 时, 当j=4时, xi, 4=5.8, 中间排序个体为φi, 4=3, 则根据公式 (5) 求的工序为πi, 4=1。

2 适应度函数
JSP问题确定完每个工件在每台机器上的工序后, 就要利用求的的工序计算完工时间。令机器顺序举证为M, 加工时间举证为T, Ti, j表示工件i在机器j上的加工时间。y=πi, j表示当前工件i的第j道工序。p (i) 表示机器i上安排的工件数。O (y) 为当前工件加工的工序数, MOi, j安排在机器i上第j个加工的工件号, N (m) 表示机器m上已经安排的工件数。

步骤1:采用L-O-V原则将初始个体Xi= (xi, 1, xi, 2, …, xi, D) 转化为加工工序集Qi, j= (Qi, 1, Qi, 2, …, Qi, D) 。p (i) =0。
步骤2:令k=1, 对于工件y=Qi, j, O (y) =O (y) +1, m=M[y, O (y) ], N (m) =N (m) +1, MO (m, N) =y。
步骤3:求每道工序的完工时间, 要考虑当前工件是否加工第1道工序, 以及工件所在加工机器是否有紧前工序。因此实际分析时要分4种情况考虑。Sy, O (y) =max{Cy, O (y) -1, Cr, k}, Cy, O (y) =Sy, O (y) +Ty, O (y) , r和k可以由以下程序所得:OPT_NUM表示工序集。

步骤4:y=y+1, N=N+1。
步骤5:k=k+1, 如果k≤D, 转到步骤3。
步骤6:Cmax=max1≤i≤n{Ci,m}作为该种群个体的适应值。

四、部分源代码

%% 定义参数
close all;
clc;
Com_Num=6;  %设备数
Mec_N=6;   %工序数
D=Mec_N*Com_Num; %变量维数
F0=0.7;     %变异系数
CR=0.65;     %交叉系数
NP=50;      %种群规模
Ob=zeros(1,NP);
P=0;PreX=exprnd(3,1,D);
BestX=PreX;
x1=zeros(NP,D);
v=zeros(NP,D);
u=zeros(NP,D);
[~,~,result]=JSP(Mec_N,Com_Num,Time,BestX);
trace(1)=result;
%%  每迭代一次退火一次，直到满足迭代终止条件
for i=1:D
%%  变异操作for ii=1:NPx1(ii,:)=PreX+exprnd(3,1,D);endfor m=1:NPr1=randi(NP);while (r1==m)r1=randi(NP);endr2=randi(NP);while (r2==m)||(r2==r1)endr3=randi(NP);while (r3==m)||(r3==r1)||(r2==r3)r3=randi(NP);endend
%% 交叉操作
r=randi(D);for n=1:D cr=rand(1);if (cr<=CR)||(n==r)u(:,n)=v(:,n);elseu(:,n)=x1(:,n);endendend
figure(1)
plot(trace,'k')
xlabel('G','FontSize',15);
ylabel('Time/s','FontSize',15);
set(gca,'FontSize',13)
figure(2)
a=trace(end);
[end_t,star_t,result4]=JSP(Mec_N,Com_Num,Time,BestX);
plot_gant(Mec_N,Com_Num,end_t,star_t,a)
xlabel('processing time（s）','FontSize',15);
ylabel('Machine','FontSize',15);
set(gca,'FontSize',13)
function plot_gant(Mec_N,Com_Num,star_t,end_t,Best)
% arti=repelem(1:Mec_N,Com_Num);
% x=[150.059399274192  196.260024970295    186.984283323539    212.747530687756    167.363571486705    163.902458562661    107.067493817499    165.250411038459    135.927777433570    165.931695634301    174.450985125151    183.808390673624    189.716304529167    154.962440055274    190.531092055579    132.252304329349    158.260502776165    172.839860476166    178.293493991184    121.502924285128    160.701217067862    144.545455498799    160.797337279283    208.410521459539    161.667493231441    160.941542777386    186.037820807237    142.248355778307    185.953155425369    181.018390705916    168.249162390406    158.793720812979    184.706038367871    194.081824712592    163.633220939848    149.797816925022];
% Order12=[x',arti'];
% Order12=sortrows(Order12,1);
% Order12=Order12';
% x1=Order12(2,:);
% re_t1=zeros(1,Mec_N);
% %% 解码
% Order=x1;
% Com_Order=ones(1,Mec_N);
%   Time1=Time;
% for ii=2:2:2*Com_Num
%   Time1(:,ii)=0;
% end
%   a=zeros(Mec_N,1);
%   b=zeros(size(Order,2),1);
%   end_t=[a,Time1];%工件加工结束时间
%   star_t=[a,Time1];%工件加工开始时间
%   re_t=zeros(1,Mec_N);
%   re_t=[0:Mec_N-1;re_t];%设备释放时间
%   Orders=[Order',b,b,b,b,b];%Orders=[工件 设备 加工开始时间 加工结束时间  属于类别（3） 加工时间]
%   Orders(:,2)=-1;
%   for i=1:Mec_N*Com_Num
%       times=100;
%       [x,y]=find(re_t(1,:)==Time(Order(i),2*Com_Order(1,Order(i))-1));
%         relese=re_t(2,y);%提释放时间
%       mechine_loca=Orders(Orders(:,2)==re_t(1,y),:);
%       if isempty(mechine_loca)%设备未调用
%           end_t(Order(i),2*Com_Order(1,Order(i))+1)=max([relese,end_t(Order(i),2*Com_Order(1,Order(i))-1),re_t1(1,y)])+Time(Order(i),2*Com_Order(1,Order(i)));
%           star_t(Order(i),2*Com_Order(1,Order(i))+1)=max([relese,end_t(Order(i),2*Com_Order(1,Order(i))-1),re_t1(1,y)]);
%           Orders(i,2)=Time(Order(i),2*Com_Order(1,Order(i))-1);%所在设备
%           Orders(i,3)=star_t(Order(i),2*Com_Order(1,Order(i))+1);%开始时间
%           Orders(i,4)=end_t(Order(i),2*Com_Order(1,Order(i))+1);%结束时间
%           if relese<end_t(Order(i),2*Com_Order(1,Order(i))+1)
%           re_t(2,y)=end_t(Order(i),2*Com_Order(1,Order(i))+1);%释放时间
%           end
%           Orders(i,6)=Time(Order(i),2*Com_Order(1,Order(i)));%测试用
%           Com_Order(1,Order(i))= Com_Order(1,Order(i))+1;%测试用
%           Orders(i,5)=1;%测试用
%       end

五、运行结果

六、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.
[3]周萧,王万良,徐新黎.解决作业车间调度问题的混合差分进化算法[J].轻工机械. 2010,28(05)