车间调度基于模拟退火求解车间调度问题

Posted 2021-05-24 Matlab走起

模拟退火算法：

       为了解决局部最优解问题， 1983年，Kirkpatrick等提出了模拟退火算法（SA）能有效的解决局部最优解问题。我们知道在分子和原子的世界中，能量越大，意味着分子和原子越不稳定，当能量越低时，原子越稳定。‘退火’是物理学术语，指对物体加温在冷却的过程。模拟退火算法来源于晶体冷却的过程，如果固体不处于最低能量状态，给固体加热再冷却，随着温度缓慢下降，固体中的原子按照一定形状排列，形成高密度、低能量的有规则晶体，对应于算法中的全局最优解。而如果温度下降过快，可能导致原子缺少足够的时间排列成晶体的结构，结果产生了具有较高能量的非晶体，这就是局部最优解。因此就可以根据退火的过程，给其在增加一点能量，然后在冷却，如果增加能量，跳出了局部最优解，这本次退火就是成功的，下面我们就详细讲讲他是如何在局部最优解跳出来到全局最优解的：

        模拟退火算法包含两个部分即Metropolis算法和退火过程。Metropolis算法就是如何在局部最优解的情况下让其跳出来，是退火的基础。1953年Metropolis提出重要性采样方法，即以概率来接受新状态，而不是使用完全确定的规则，称为Metropolis准则，计算量较低。下面先形象的说一下，然后在因此数学公式：

    假设开始状态在A，随着迭代次数更新到B的局部最优解，这时发现更新到B时，能力比A要低，则说明接近最优解了，因此百分百转移，状态到达B后，发现下一步能量上升了，如果是梯度下降则是不允许继续向前的，而这里会以一定的概率跳出这个坑，这各概率和当前的状态、能量等都有关系，下面会详细说，如果B最终跳出来了到达C，又会继续以一定的概率跳出来，可能有人会迷惑会不会跳回之前的B呢？下面会解释，直到到达D后，就会稳定下来。所以说这个概率的设计是很重要的，下面从数学方面进行解释。

    假设前一个状态为,系统根据某一指标（梯度下降，上节的能量），状态变为,相应的，系统的能量由变为,定义系统由变为的接受概率P为：

                                         

从上式我们可以看到，如果能量减小了，那么这种转移就被接受（概率为1），如果能量增大了，就说明系统偏离全局最优值位置更远了，此时算法不会立刻将其抛弃，而是进行概率操作：首先在区间【0,1】产生一个均匀分布的随机数，如果P，则此种转移接受，否则拒绝转移，进入下一步，往复循环。其中P以能量的变化量和T进行决定概率P的大小，所以这个值是动态的。

退火算法的参数控制

       Metropolis算法是模拟退火算法的基础，但是直接使用Metropolis算法 可能会导致寻优速度太慢，以至于无法实际使用，为了确保在有限的时间收敛，必须设定控制算法收敛的参数，在上面的公式中，可以调节的参数就是T，T如果过大，就会导致退火太快，达到局部最优值就会结束迭代，如果取值较小，则计算时间会增加，实际应用中采用退火温度表，在退火初期采用较大的T值，随着退火的进行，逐步降低，具体如下：

    （1）初始的温度T(0)应选的足够高，使的所有转移状态都被接受。初始温度越高，获得高质量的解的概率越大，耗费的时间越长。

    （2） 退火速率。 最简单的下降方式是指数式下降：

                                                    

                 其中是小于1的正数，一般取值为0.8到0.99之间。使的对每一温度，有足够的转移尝试，指数式下降的收敛速度比较慢，其他下降方式如下：

                                                      

                                                        

      （3）终止温度

                        如果在若干次迭代的情况下每有可以更新的新状态或者达到用户设定的阈值，则退火完成。

模拟退火的步骤：   

               1.模拟退火算法可以分解为解空间、目标函数和初始解三部分。

               2.模拟退火的基本思想:

                    (1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)，每个T值的迭代次数L

                    (2) 对k=1, …, L做第(3)至第6步：

                    (3) 产生新解S′

                    (4) 计算增量ΔT=C(S′)-C(S)，其中C(S)为代价函数

                    (5) 若ΔT<0则接受S′作为新的当前解，否则以概率exp(-ΔT/T)接受S′作为新的当前解.

                    (6) 如果满足终止条件则输出当前解作为最优解，结束程序。

                    终止条件通常取为连续若干个新解都没有被接受时终止算法。

                     (7) T逐渐减少，且T->0，然后转第2步。

模拟退火算法新解的产生和接受可分为如下四个步骤：

第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。

第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。

第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若ΔT<0则接受S′作为新的当前解S，否则以概率exp(-ΔT/T)接受S′作为新的当前解S。

第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。

模拟退火算法与初始值无关，算法求得的解与初始解状态S(是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。

退火算法程序流程图;

作业车间调度问题描述

作业车间调度问题（Job Shop Scheduling, JSP）是最经典的几个NP-hard问题之一。其应用领域极其广泛，涉及航母调度，机场飞机调度，港口码头货船调度，汽车加工流水线等。

JSP问题描述：一个加工系统有M台机器，要求加工N个作业，其中，作业i包含工序数为Li。令，则L为任务集的总工序数。其中，各工序的加工时间已确定，并且每个作业必须按照工序的先后顺序加工。调度的任务是安排所有作业的加工调度排序，约束条件被满足的同时，使性能指标得到优化。

作业车间调度需要考虑如下约束:

Cons1：每道工序在指定的机器上加工，且必须在其前一道工序加工完成后才能开始加工；

Cons2：某一时刻1台机器只能加工1个作业；

Cons3：每个作业只能在1台机器上加工1次；

Cons4：各作业的工序顺序和加工时间已知，不随加工排序的改变而改变。

问题实例

下面给出作业车间调度问题的一个实例，其中每个工序上标注有一对数值（m,p），其中，m表示当前工序必须在第m台机器上进行加工，p表示第m台机器加工当前工序所需要的加工时间。（注：机器和作业的编号从0开始）
jop0=[(0,3),(1,2),(2,2)]
jop1=[(0,2),(2,1),(1,4)]
jop2=[(1,4),(2,3)]
在这个例子中，作业jop0有3道工序：它的第1道工序上标注有(0,3)，其表示第1道工序必须在第0台机器上进行加工，且需要3个单位的加工时间；它的第2道工序上标注有(1,2)，其表示第2道工序必须在第1台机器上进行加工，且需要2个单位的加工时间；余下的同理。总的来说，这个实例中共有8道工序。
该问题的一个可行解是L=8道工序开始时间的一个排列，且满足问题的约束。下图给出了一个可行解（注：该解不是最优解）的示例：

• clc
  clear
  %=========数据录入，参数调整=================
  swarminitNum=20;%初始生成的粒子数；
   
  MM=[1 2 3 4 5 6 
      6 6 6 6 6 6];%工件、工序数量矩阵，MM第一行表示工件，第二行表示每个工件的工序数；
   
  machineNum=6;   %加工机器数；
   
  initT=1000;          %模拟退火初始温度；
   
  gen=500;        %循环迭代数；
   
  w1=0.35;          %变异率；
   
  changeNum=3;     %变异变换对数；
   
  restrictmatrixM=[3     1     2     4     6     5
                   2     3     5     6     1     4
                   3     4     6     1     2     5
                   2     1     3     4     5     6
                   3     2     5     6     1     4
                   2     4     6     1     5     3];%job-shop机器约束矩阵；
               
  restrictmatrixT=[1     3     6    7      3     6
                   8     5     10   10     10    4
                   5     4     8     9     1     7
                   5     5     5     3     8     9
                   9     3     5     4     3     1
                   3     3     9     10    4     1];%job-shop时间约束矩阵；
               
  %===============PSO算法==========================
  swarminit=cell(1,swarminitNum);
  swarminitLong=sum(MM(2,:));          %所有工序数即粒子长度；
  for i=1:swarminitNum,
      swarminit{i}=randomparticle(MM) ;
  end                                  %随机生成初始粒子群体
  [popu,s] = size(swarminit); 
  trace = ones(1,gen); 
  trace(1) = 10000; % 初始全局最佳适应度设为足够大 
  for i = 1:s,
      bestfit(i) = 10000; % 初始个体历史最佳适应度设为足够大 
  end
  bestpar = swarminit; % 个体历史最佳粒子初始化
  for u=1:swarminitNum,
      fitlist=[0]; 
  end
  T=initT;
  for step = 1:gen,
      for q=1:swarminitNum,
              
         swarminit{j}=cross(bestparticle1,swarminit{j},l4,l3);%粒子交叉；
     end 
  end
  [a,b,c]=timedecode2(bestparticle,restrictmatrixM,restrictmatrixT,machineNum);
  disp(['优化目标: 最小平均流动时间'])
  disp(['粒子数:' int2str(swarminitNum)  '   循环代数:' int2str(gen)])
  disp(['变异率:' num2str(w1)  '   变异变换对数:' int2str(changeNum)])
  disp(['模拟退火初始值:' int2str(initT)  '  模拟退火终值:' int2str(T)])
  disp(['迭代循环值:' int2str(trace)])
  disp(['最小平均流动时间:' int2str(a)  '   最大完工时间:'  int2str(b)  '   最小间隙时间:' int2str(c) ])
  disp(['最优粒子' int2str(bestparticle)])
  pause
  gant(bestparticle,swarminitLong,restrictmatrixM,restrictmatrixT,b)

  

  完整代码添加QQ1575304183