【车间调度】基于matlab遗传算法求解混合流水车间调度最优问题【含Matlab源码 901期】

一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【车间调度】基于matlab遗传算法求解混合流水车间调度最优问题【含Matlab源码 901期】

获取代码方式2：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、车间调度简介

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]的机器组合。调度的目的还是总的完工时间最短（也可以是其他目标，比如机器最大负荷最短、总的机器负荷最短）

三、遗传算法简介

1 遗传算法概述
遗传算法（Genetic Algorithm，GA）是进化计算的一部分，是模拟达尔文的遗传选择和自然淘汰的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。该算法简单、通用，鲁棒性强，适于并行处理。

2 遗传算法的特点和应用
遗传算法是一类可用于复杂系统优化的具有鲁棒性的搜索算法，与传统的优化算法相比，具有以下特点：
（1）以决策变量的编码作为运算对象。传统的优化算法往往直接利用决策变量的实际值本身来进行优化计算，但遗传算法是使用决策变量的某种形式的编码作为运算对象。这种对决策变量的编码处理方式，使得我们在优化计算中可借鉴生物学中染色体和基因等概念，可以模仿自然界中生物的遗传和进化激励，也可以很方便地应用遗传操作算子。
（2）直接以适应度作为搜索信息。传统的优化算法不仅需要利用目标函数值，而且搜索过程往往受目标函数的连续性约束，有可能还需要满足“目标函数的导数必须存在”的要求以确定搜索方向。遗传算法仅使用由目标函数值变换来的适应度函数值就可确定进一步的搜索范围，无需目标函数的导数值等其他辅助信息。直接利用目标函数值或个体适应度值也可以将搜索范围集中到适应度较高部分的搜索空间中，从而提高搜索效率。
（3）使用多个点的搜索信息，具有隐含并行性。传统的优化算法往往是从解空间的一个初始点开始最优解的迭代搜索过程。单个点所提供的搜索信息不多，所以搜索效率不高，还有可能陷入局部最优解而停滞；遗传算法从由很多个体组成的初始种群开始最优解的搜索过程，而不是从单个个体开始搜索。对初始群体进行的、选择、交叉、变异等运算，产生出新一代群体，其中包括了许多群体信息。这些信息可以避免搜索一些不必要的点，从而避免陷入局部最优，逐步逼近全局最优解。
（4） 使用概率搜索而非确定性规则。传统的优化算法往往使用确定性的搜索方法，一个搜索点到另一个搜索点的转移有确定的转移方向和转移关系，这种确定性可能使得搜索达不到最优店，限制了算法的应用范围。遗传算法是一种自适应搜索技术，其选择、交叉、变异等运算都是以一种概率方式进行的，增加了搜索过程的灵活性，而且能以较大概率收敛于最优解，具有较好的全局优化求解能力。但，交叉概率、变异概率等参数也会影响算法的搜索结果和搜索效率，所以如何选择遗传算法的参数在其应用中是一个比较重要的问题。
综上，由于遗传算法的整体搜索策略和优化搜索方式在计算时不依赖于梯度信息或其他辅助知识，只需要求解影响搜索方向的目标函数和相应的适应度函数，所以遗传算法提供了一种求解复杂系统问题的通用框架。它不依赖于问题的具体领域，对问题的种类有很强的鲁棒性，所以广泛应用于各种领域，包括：函数优化、组合优化生产调度问题、自动控制
、机器人学、图像处理（图像恢复、图像边缘特征提取…)、人工生命、遗传编程、机器学习。

3 遗传算法的基本流程及实现技术
基本遗传算法（Simple Genetic Algorithms,SGA）只使用选择算子、交叉算子和变异算子这三种遗传算子，进化过程简单，是其他遗传算法的基础。

3.1 遗传算法的基本流程
通过随机方式产生若干由确定长度（长度与待求解问题的精度有关）编码的初始群体；
通过适应度函数对每个个体进行评价，选择适应度值高的个体参与遗传操作，适应度低的个体被淘汰；
经遗传操作（复制、交叉、变异）的个体集合形成新一代种群，直到满足停止准则（进化代数GEN>=?）；
将后代中变现最好的个体作为遗传算法的执行结果。

其中，GEN是当前代数；M是种群规模，i代表种群数量。

3.2 遗传算法的实现技术
基本遗传算法（SGA）由编码、适应度函数、遗传算子（选择、交叉、变异）及运行参数组成。
3.2.1 编码
（1）二进制编码
二进制编码的字符串长度与问题所求解的精度有关。需要保证所求解空间内的每一个个体都可以被编码。
优点：编、解码操作简单，遗传、交叉便于实现
缺点：长度大
（2）其他编码方法
格雷码、浮点数编码、符号编码、多参数编码等
3.2.2 适应度函数
适应度函数要有效反映每一个染色体与问题的最优解染色体之间的差距。
3.2.3选择算子

3.2.4 交叉算子
交叉运算是指对两个相互配对的染色体按某种方式相互交换其部分基因，从而形成两个新的个体；交叉运算是遗传算法区别于其他进化算法的重要特征，是产生新个体的主要方法。在交叉之前需要将群体中的个体进行配对，一般采取随机配对原则。
常用的交叉方式：
单点交叉
双点交叉（多点交叉，交叉点数越多，个体的结构被破坏的可能性越大，一般不采用多点交叉的方式）
均匀交叉
算术交叉
3.2.5 变异算子
遗传算法中的变异运算是指将个体染色体编码串中的某些基因座上的基因值用该基因座的其他等位基因来替换，从而形成一个新的个体。

就遗传算法运算过程中产生新个体的能力方面来说，交叉运算是产生新个体的主要方法，它决定了遗传算法的全局搜索能力；而变异运算只是产生新个体的辅助方法，但也是必不可少的一个运算步骤，它决定了遗传算法的局部搜索能力。交叉算子与变异算子的共同配合完成了其对搜索空间的全局搜索和局部搜索，从而使遗传算法能以良好的搜索性能完成最优化问题的寻优过程。

3.2.6 运行参数

4 遗传算法的基本原理
4.1 模式定理

4.2 积木块假设
具有低阶、定义长度短，且适应度值高于群体平均适应度值的模式称为基因块或积木块。
积木块假设：个体的基因块通过选择、交叉、变异等遗传算子的作用，能够相互拼接在一起，形成适应度更高的个体编码串。
积木块假设说明了用遗传算法求解各类问题的基本思想，即通过积木块直接相互拼接在一起能够产生更好的解。

四、部分源代码

unction [Zp,Y1p,Y2p,Y3p,Xp,LC1,LC2]=JSPGA(M,N,Pm,T,P)

%----------------------------------------------------------------
%   JSPGA.m 流水车间调度遗传算法
%----------------------------------------------------------------
% 输入参数列表
%    M  遗传进化迭代次数
%    N  种群规模（偶数）
%    Pm 变异概率
%    T  m*n的矩阵，存储m个工件n个工序的加工时间
%    P  1*n的向量，n个工序中，每一个工序所具有的机床数目
%----------------------------------------------------------------
% 输出参数列表
%   Zp   最优的Makespan值
%   Y1p  最优方案中，各工件各工序的开始时刻，可用来绘制甘特图
%   Y2p  最优方案中，各工件各工序的结束时刻
%   Y3p  最优方案中，各工件各工序使用的机器编号
%   Xp   最优决策变量的值，决策变量是一个实数编码的m*n矩阵
%   LC1  收敛曲线1，各代最优个体适应值的记录
%   LC2  收敛曲线2，各代群体平均适应值的记录
% 最后程序将绘制三幅图片：两条收敛曲线和各工件调度的甘特图

%----------------------------------------------------------------
% 第一步：变量初始化
[m,n]=size(T);  %m是总工件数，n是总工序数
Xp=zeros(m,n);  %最优决策变量
LC1=zeros(1,M); %收敛曲线1
LC2=zeros(1,N); %收敛曲线2
%----------------------------------------------------------------
% 第二步：随机产生初始种群farm
farm=cell(1,N); %采用细胞结构存储种群
for k=1:N
    X=zeros(m,n);
    for j=1:n
        for i=1:m
            X(i,j)=1+(P(j)-eps)*rand; %eps默认eps(1)表示matlab做能表达的1的精度即1+0.9*eps=1
        end
    end
    farm{k}=X;
end
counter=0;  %设置迭代计数器
while counter<M     %停止条件为打到最大迭代次数
    %----------------------------------------
    % 第三步：交叉
    newfarm=cell(1,N);  %交叉产生的新种群存在其中
    Ser=randperm(N);    %randperm(N)返回1：N的随机排序向量
    for i=1:2:(N-1)
        A=farm{Ser(i)};     %父代个体
        B=farm{Ser(i+1)};   %父代相邻个体
        Manner=unidrnd(2);  %随机选择交叉方式（横向交叉/纵向交叉）  unidrnd(N)返回1：N的离散随机数
        if Manner==1
            cp=unidrnd(m-1);    %随机选择交叉点
            %双亲双子单点交叉
            a=[A(1:cp,:);B((cp+1):m,:)];    %子代个体
            b=[B(1:cp,:);A((cp+1):m,:)];
        else
            cp=unidrnd(n-1);    %随机选择交叉点
            a=[A(:,1:cp),B(:,(cp+1):n)];    %双亲双子单点交叉
            b=[B(:,1:cp),A(:,(cp+1):n)];
        end
        newfarm{i}=a;   %交叉后的子代存入newfarm
        newfarm{i+1}=b;
    end
    %新旧种群合并
    FARM=[farm,newfarm];
    %----------------------------------------
    % 第四步：选择复制
    FITNESS=zeros(1,2*N);
    fitness=zeros(1,N);
    plotif=0;
    for i=1:(2*N)
        X=FARM{i};
        Z=COST(X,T,P,plotif);   %调用计算费用的子函数，计算适应值
        FITNESS(i)=Z;       
    end
    %选择复制采取两两随机配对竞争的方式，具有保留最优个体的能力
    Ser=randperm(2*N);
    for i=1:N
        f1=FITNESS(Ser(2*i-1));
        f2=FITNESS(Ser(2*i));
        if f1<=f2
            farm{i}=FARM{Ser(2*i-1)};
            fitness(i)=FITNESS(Ser(2*i-1));
        else
            farm{i}=FARM{Ser(2*i)};
            fitness(i)=FITNESS(Ser(2*i));
        end
    end
    %记录最佳个体和收敛曲线
    minfitness=min(fitness);    %找出10个种群中最小的个体
    meanfitness=mean(fitness);  %10个种群的平均适应值
    LC1(counter+1)=minfitness;  %收敛曲线1，各代最优个体适应值的记录
    LC2(counter+1)=meanfitness; %收敛曲线2，各代群体平均适应值的记录
    pos=find(fitness==minfitness); %记录10个种群中最小个体的位置
    Xp=farm{pos(1)};            %将最小个体位置的第一个排序结果记录
    %----------------------------------------
    % 第五步：变异
    for i=1:N
        if Pm>rand %变异概率为Pm
         
            X(I,J)=1+(P(J)-eps)*rand;
            farm{i}=X;
        end
    end
    farm{pos(1)}=Xp;        %变异之后，具有保留最优个体的能力
    counter=counter+1;
end
% 第一步：变量初始化
[m,n]=size(X);
Y1p=zeros(m,n);
Y2p=zeros(m,n);
Y3p=zeros(m,n);
%----------------------------------------------------------------
% 第二步：计算第一道工序的安排
Q1=zeros(m,1);
Q2=zeros(m,1);
R=X(:,1);   %取出第一道工序
Q3=floor(R);    %向下取整即得到各工件在第一道工序使用的机器编号
% 下面计算各工件第一道工序的开始时刻和结束时刻
for i=1:P(1)        %取出机器编号
    pos=find(Q3==i);%取出使用编号为i的机器为其加工的工件编号 find返回向量/矩阵中的位置
    lenpos=length(pos);
    if lenpos>=1
        Q1(pos(1))=0;
        Q2(pos(1))=T(pos(1),1);
        if lenpos>=2
            for j=2:lenpos
                Q1(pos(j))=Q2(pos(j-1));
                Q2(pos(j))=Q2(pos(j-1))+T(pos(j),1);
            end
        end
    end
end
Y1p(:,1)=Q1;
Y2p(:,1)=Q2;
Y3p(:,1)=Q3;
%----------------------------------------------------------------
% 第三步：计算剩余工序的安排
for k=2:n
    R=X(:,k);       %取出第k道工序
    Q3=floor(R);    %向下取整即得到各工件在第k道工序使用的机器编号
    %下面计算各工件第k道工序的开始时刻和结束时刻
    for i=1:P(k)    %取出机器编号
        pos=find(Q3==i);    %取出使用编号为i的机器为其加工的工件编号
        lenpos=length(pos);
        if lenpos>=1
            EndTime=Y2p(pos,k-1);   %取出这些机器在上一个工序中的结束时刻
            POS=zeros(1,lenpos);    %上一个工序完成时间由早到晚的排序
            for jj=1:lenpos
                MinEndTime=min(EndTime);
                ppp=find(EndTime==MinEndTime);
                POS(jj)=ppp(1);
                EndTime(ppp(1))=Inf;
            end
            %根据上一个工序完成时刻的早晚，计算各工件第k道工序的开始时刻和结束时刻
            Q1(pos(POS(1)))=Y2p(pos(POS(1)),k-1);
            Q2(pos(POS(1)))=Q1(pos(POS(1)))+T(pos(POS(1)),k);%前一个工件的结束时刻
            if lenpos>=2
                for j=2:lenpos
                    Q1(pos(POS(j)))=Y2p(pos(POS(j)),k-1);    %预定的开始时刻为上一个工序的结束时刻
                    Q2(pos(POS(j)))=Q1(pos(POS(j)))+T(pos(POS(j)),k);%前一个工件的结束时刻
                    if Q1(pos(POS(j)))<Q2(pos(POS(j-1)))   %如果比前面的工序的结束时刻还早
                       Q1(pos(POS(j)))=Q2(pos(POS(j-1)));
                       Q2(pos(POS(j)))=Q1(pos(POS(j)))+T(pos(POS(j)),k);%前一个工件的结束时刻
                    end
                end
            end
        end
    end

五、运行结果

六、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.

文章来源: qq912100926.blog.csdn.net，作者：海神之光，版权归原作者所有，如需转载，请联系作者。

原文链接：qq912100926.blog.csdn.net/article/details/116993511