车间调度-灰狼算法的应用:以算例MK01为例


车间调度系列文章:

柔性车间调度问题

柔性车间调度问题可描述为:多个工件在多台机器上加工,工件安排加工时严格按照工序的先后顺序,至少有一道工序有多个可加工机器,在某些优化目标下安排生产。
柔性车间调度问题的约束条件如下:

  • (1)同一台机器同一时刻只能加工一个工件;
  • (2)同一工件的同一道工序在同一时刻被加工的机器数是一;
  • (3)任意工序开始加工不能中断;
  • (4)各个工件之间不存在的优先级的差别;
  • (5)同一工件的工序之间存在先后约束,不同工件的工序之间不存在先后约束;
  • (6)所有工件在零时刻都可以被加工。

MK01算例:

10 6 2
6 2 1 5 3 4 3 5 3 3 5 2 1 2 3 4 6 2 3 6 5 2 6 1 1 1 3 1 3 6 6 3 6 4 3
5 1 2 6 1 3 1 1 1 2 2 2 6 4 6 3 6 5 2 6 1 1
5 1 2 6 2 3 4 6 2 3 6 5 2 6 1 1 3 3 4 2 6 6 6 2 1 1 5 5
5 3 6 5 2 6 1 1 1 2 6 1 3 1 3 5 3 3 5 2 1 2 3 4 6 2
6 3 5 3 3 5 2 1 3 6 5 2 6 1 1 1 2 6 2 1 5 3 4 2 2 6 4 6 3 3 4 2 6 6 6
6 2 3 4 6 2 1 1 2 3 3 4 2 6 6 6 1 2 6 3 6 5 2 6 1 1 2 1 3 4 2
5 1 6 1 2 1 3 4 2 3 3 4 2 6 6 6 3 2 6 5 1 1 6 1 3 1
5 2 3 4 6 2 3 3 4 2 6 6 6 3 6 5 2 6 1 1 1 2 6 2 2 6 4 6
6 1 6 1 2 1 1 5 5 3 6 6 3 6 4 3 1 1 2 3 3 4 2 6 6 6 2 2 6 4 6
6 2 3 4 6 2 3 3 4 2 6 6 6 3 5 3 3 5 2 1 1 6 1 2 2 6 4 6 2 1 3 4 2

第一行的10,6是工件数和机器数。

第二行第一个加粗的数字6表示,工件1有6道工序。斜体的2 1 5 3 4,表示工件1的第一道工序有两个可选机器,分别是1和3,加工时间是5和4,后面的3 5 3 3 5 2 1表示工件1的第二道工序有3个可选机器,分别是5,3,2,加工时间是3,5,1,一行就是1个工件的所有工序的可选机器可加工时间,后面的工序以此类推。

下面的每一行以此类推。

布谷鸟算法

灰狼优化算法是澳大利亚格里菲斯大学学者Mirjalili等人在2014年提出的一种优化搜索算法,其启发来自自然界中灰狼的捕猎行为。该算法的收敛性能较强、参数少及容易实现的优点使其一出现就受到学者的大量关注,通过近些年的研究,其已经很好地应用在在车间调度问题的求解、图像分类、参数优化等领域。

灰狼优化算法的优化过程一般分为社会等级分层、猎物跟踪、猎物包围和最后的猎物狩取阶段。灰狼的社会分层有四层,分别对应αβγω狼。α狼是灰狼捕猎行为的决策者,等级最高。β狼是α狼的候选人,平时辅助α狼决策,服从α狼。γ狼服从αβ,是其余狼群的领导者。ω狼是等级最低的阶层。猎物跟踪和狩猎过程都由种群中适应度最高的αβγ三只狼引导完成,每次狩猎完成后狼群中个体位置调整的依据与三只狼的距离及相关参数。

为使算法跳出局部最优,跟踪过程设置随机向量A, |A|>1时搜索代理能够远离猎物,实现全局搜索,同时设置在[0,2]上的随机值构成的向量C是一个随机搜索系数,增加了猎物搜索的随机性,避免了算法的局部最优。

公式是网上的截图:
在这里插入图片描述

柔性作业车间工具

本文写了甘特图的画图函数;工序,机器,加工时间编码的生成函数;编码的解码函数。甘特图和解码前面推文有介绍,为了能在布谷鸟算法使用,下面介绍一下编码的生成:

  • 步骤1:按照工件的工序数依次生成工序编码如下:
work = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9]

程序里为方便运算,0表示工件1,依次类推。

  • 步骤2:work长度是55,即编码长度。在在0到1之间生成55个数。
initial_a= [0.7937249945850672, 0.16522396435630604, 0.6600917203383786, 0.6596651342253383, 0.9774747473336962, 0.9883891883783963, 0.7156533671841291, 0.04997505800565305, 0.2578073167205841, 0.1836103950459067, 0.20165471162164583, 0.8038018885972078, 0.41850414257851287, 0.95008978291153, 0.11549641941696531, 0.5453005466735185, 0.5471875726880646, 0.8365196679989834, 0.9380606558784159, 0.9499683399591661, 0.4119930638242061, 0.06769082273628235, 0.1248094943364636, 0.5137257432003801, 0.6318423906858445, 0.9488907758993115, 0.08482854147556895, 0.2571981915489773, 0.9716139935057354, 0.15806881092191782, 0.9624616622251282, 0.16955974676704666, 0.8796549710395518, 0.07510161533653481, 0.15127075291529657, 0.836649773542874, 0.2844859458891926, 0.6843108506156074, 0.8546830769047646, 0.6680995444587688, 0.05461130146431059, 0.028553253256760835, 0.6630167299173496, 0.885336533372363, 0.9962311484426508, 0.7236930059814548, 0.5636537858241203, 0.8528428792173608, 0.6304944808908045, 0.5359472633322336, 0.9163032901730644, 0.8258172304840811, 0.7824920279848068, 0.8950392199000986, 0.05567547868012779]
  • 步骤3:对步骤2得到的编码,进行从小到大的排列,按照从小到大的顺序取出其在编码的位置。
index_work=[41, 7, 40, 54, 21, 33, 26, 14, 22, 34, 29, 1, 31, 9, 10, 27, 8, 36, 20, 12, 23, 49, 15, 16, 46, 48, 24, 3, 2, 42, 39, 37, 6, 45, 52, 0, 11, 51, 17, 35, 47, 38, 32, 43, 53, 50, 18, 25, 19, 13, 30, 28, 4, 5, 44]

initial_a的第42(41+1)个数最小,第7(6+1)个数第二小,依次类推。

  • 步骤4:依据步骤3的得到编码位置,在work中找到对应工序编码,即可得到长度55的工序编码。
job= [7, 1, 7, 9, 4, 6, 4, 2, 4, 6, 5, 0, 5, 1, 1, 5, 1, 6, 3, 2, 4, 9, 2, 3, 8, 8, 4, 0, 0, 7, 7, 6, 1, 8, 9, 0, 2, 9, 3, 6, 8, 7, 5, 8, 9, 9, 3, 4, 3, 2, 5, 5, 0, 0, 8]

work的第42(41+1)个数是7,第7(6+1)个数是1,依次类推。

代码在fjsp.py里。

机器和加工时间编码:

随机生成0到1之间的数,如果该随机数小于pi,选择加工时间最小的机器,否则随机选择加工机器,pi自行调整。

算法运行环境

本代码运行环境是win7系统,python3.5.2,第三方库:

numpy==1.18.5
matplotlib==3.0.3

改进灰狼算法

本文对工序编码用灰狼算法更新,机器编码迭代加入自适应算子:遗传算法的均匀交叉

算法步骤:

  • 步骤1:初始化多个工序,机器,加工时间编码
  • 步骤2:对工序编码用灰狼算法更新
  • 步骤3:对机器编码用均匀交叉进行更新,如果新的编码的完工时间小于原来的编码,则取代原编码,反之则反
  • 步骤4:判断是否达到最大迭代次数,是的话输出结果,否则转到步骤1

α、β、γ狼和a的变化
每次迭代时,按完工时间从小到大对上一次迭代种群的个体排序,所以前三个个体就是灰狼算法中的αβγ狼。首次迭代对初始的种群排序。核心代码:

index_sort=np.array(answer).argsort()        #返回完工时间从小到大的位置索引
            work_job1,work_M1,work_T1=work_job[index_sort],work_M[index_sort],work_T[index_sort]
            answer1=np.array(answer)[index_sort]
            job_init1=job_init[index_sort]

            Alpha=job_init1[0]   #α狼
            Beta=job_init1[1]    #β狼
            Delta=job_init1[2]   #δ狼
            a = 2*(1-gen/self.generation)
  

a是由2到0的变化,gen表示第gen次迭代,self.generation是总迭代次数。

灰狼算法更新方法

公式上面有介绍,核心代码:

				r1 = random.random()   #灰狼算法解的更新
                r2 = random.random()
                A1 = 2 * a * r1 - a
                C1 = 2 * r2
                D_alpha =C1*Alpha-x
                x1 = x - A1 * D_alpha

                r1 = random.random()
                r2 = random.random()
                A2 = 2 * a * r1 - a
                C2 = 2 * r2
                D_beta =C2*Beta-x
                x2 = x - A2 * D_beta

                r1 = random.random()
                r2 = random.random()
                A3 = 2 * a * r1 - a
                C3 = 2 * r2
                D_delta =C3*Delta-x
                x3 = x - A3 * D_alpha

                initial_a=(x1+x2+x3)/3   #更新公式

具体的代码细节在代码在GWO里。

结果

改进算法的运行命令如下:

ho=pso([10,6,0.5],50,100)    

#第一个中括号是工件数,机器数,选择最短机器的概率
#数50,100分别代表迭代的次数和种群的规模
a,b,c,d=ho.pso_total()  #最后一次迭代的最优解

job,machine,machine_time=np.array([a]),np.array([b]),np.array([c])
to.draw(job,machine,machine_time) #画图

result=np.array(result).reshape(len(result),2)
plt.plot(result[:,0],result[:,1])                   #画完工时间随迭代次数的变化
font1={'weight':'bold','size':22}
plt.xlabel("迭代次数",font1)
plt.title("完工时间变化图",font1)
plt.ylabel("完工时间",font1)
plt.show()

结果如下:

在这里插入图片描述

结果的甘特图如下:

在这里插入图片描述

完工时间随迭代次数的变化图如下:

在这里插入图片描述

从上图可以看出:代码运行结果是Mk01.txt数据集求解结果中,最优的

代码

有3个代码和一个mk01的text文档:

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-CtHKEdZI-1647748892877)(C:\Users\Administrator\Desktop\灰狼算法12.18\5)]

篇幅问题,给出灰狼算法的代码GWO:

from data_solve import data_deal
from fjsp import FJSP
import numpy as np
import random
import math
import matplotlib.pyplot as plt 

class pso():
    def __init__(self,parm_fjsp,generation,popsize):
        self.job_num=parm_fjsp[0]                       #工件数
        self.machine_num=parm_fjsp[1]                   #机器数
        self.pi=parm_fjsp[2]
        self.generation=generation                  #迭代次数
        self.popsize = popsize                      # 种群规模
        
    def to_MT(self,W1,M1,T1): #把加工机器编码和加工时间编码转化为对应列表,目的是记录工件的加工时间和加工机器
        Ma_W1,Tm_W1,WCross=[],[],[]
        for i in range(self.job_num):#添加工件个数的空列表
            Ma_W1.append([]),Tm_W1.append([]),WCross.append([]);
        for i in range(W1.shape[1]):
            signal1=int(W1[0,i])-1
            Ma_W1[signal1].append(M1[0,i]),Tm_W1[signal1].append(T1[0,i]); #记录每个工件的加工机器
            index=np.random.randint(0,2,1)[0]
            WCross[signal1].append(index)       #随机生成一个为0或者1的列表,用于后续的机器的均匀交叉
        return Ma_W1,Tm_W1,WCross
    def back_MT(self,W1,Ma_W1,Tm_W1):  #列表返回机器及加工时间编码
        memory1=np.zeros((1,self.job_num),dtype=np.int)
        m1,t1=np.zeros((1,W1.shape[1])),np.zeros((1,W1.shape[1]))
        for i in range(W1.shape[1]):
            signal1=int(W1[0,i])-1
            m1[0,i]=Ma_W1[signal1][memory1[0,signal1]] #读取对应工序的加工机器
            t1[0,i]=Tm_W1[signal1][memory1[0,signal1]]
            memory1[0,signal1]+=1
        return m1,t1
    def mac_cross(self,Ma_W1,Tm_W1,Ma_W2,Tm_W2,WCross):  #机器均匀交叉
        MC1,MC2,TC1,TC2=[],[],[],[]
        for i in range(self.job_num):     
            MC1.append([]),MC2.append([]),TC1.append([]),TC2.append([]);
            for j in range(len(WCross[i])):
                if(WCross[i][j]==0):  #为0时继承另一个父代的加工机器选择
                    MC1[i].append(Ma_W1[i][j]),MC2[i].append(Ma_W2[i][j]),TC1[i].append(Tm_W1[i][j]),TC2[i].append(Tm_W2[i][j]);
                else:                #为1时继承父代的机器选择
                    MC2[i].append(Ma_W1[i][j]),MC1[i].append(Ma_W2[i][j]),TC2[i].append(Tm_W1[i][j]),TC1[i].append(Tm_W2[i][j]);
        return MC1,TC1,MC2,TC2
    def gwo_total(self):
        global to
        oj=data_deal(self.job_num,self.machine_num)
        Tmachine,Tmachinetime,tdx,work,tom=oj.cacu()
        parm_data=[Tmachine,Tmachinetime,tdx,work,tom]
        to=FJSP(self.job_num,self.machine_num,self.pi,parm_data)
        answer,result=[],[]
        job_init=np.zeros((self.popsize,len(work)))
        work_job1,work_M1,work_T1=np.zeros((self.popsize,len(work))),np.zeros((self.popsize,len(work))),np.zeros((self.popsize,len(work)))
        work_job,work_M,work_T=np.zeros((self.popsize,len(work))),np.zeros((self.popsize,len(work))),np.zeros((self.popsize,len(work)))
        for gen in range(self.generation):
            if(gen<1):                      #第一次生成多个可行的工序编码,机器编码,时间编码
                for i in range(self.popsize):
                    job,machine,machine_time,initial_a=to.creat_job()
                    C_finish,_,_,_,_=to.caculate(job,machine,machine_time)
                    answer.append(C_finish)
                    work_job[i],work_M[i],work_T[i]=job[0],machine[0],machine_time[0]
                    job_init[i]=initial_a
                print('种群初始的完工时间:%.0f'%(min(answer)))
                result.append([gen,min(answer)])#记录初始解的最小完工时间

            index_sort=np.array(answer).argsort()        #返回完工时间从小到大的位置索引
            work_job1,work_M1,work_T1=work_job[index_sort],work_M[index_sort],work_T[index_sort]
            answer1=np.array(answer)[index_sort]
            job_init1=job_init[index_sort]

            Alpha=job_init1[0]   #α狼
            Beta=job_init1[1]    #β狼
            Delta=job_init1[2]   #δ狼
            a = 2*(1-gen/self.generation)
            
            for i in range(3,self.popsize):     #用最优位置进行工序编码的更新
                job,machine,machine_time=work_job1[i:i+1],work_M1[i:i+1],work_T1[i:i+1]
                Ma_W1,Tm_W1,WCross=self.to_MT(job,machine,machine_time)
                x=job_init1[i]

                r1 = random.random()   #灰狼算法解的更新
                r2 = random.random()
                A1 = 2 * a * r1 - a
                C1 = 2 * r2
                D_alpha =C1*Alpha-x
                x1 = x - A1 * D_alpha

                r1 = random.random()
                r2 = random.random()
                A2 = 2 * a * r1 - a
                C2 = 2 * r2
                D_beta =C2*Beta-x
                x2 = x - A2 * D_beta

                r1 = random.random()
                r2 = random.random()
                A3 = 2 * a * r1 - a
                C3 = 2 * r2
                D_delta =C3*Delta-x
                x3 = x - A3 * D_alpha

                initial_a=(x1+x2+x3)/3   #更新公式
                index_work=initial_a.argsort()
                job_new=[]
                for j in range(len(work)):
                    job_new.append(work[index_work[j]])
                job_new=np.array(job_new).reshape(1,len(work))
                machine_new,time_new=self.back_MT(job_new,Ma_W1,Tm_W1)
                C_finish,_,_,_,_=to.caculate(job_new,machine_new,time_new)
                
                work_job1[i]=job_new[0]  #更新工序编码
                job_init1[i]=initial_a
                work_M1[i],work_T1[i]=machine_new[0],time_new[0]
                answer1[i]=C_finish
            for i in range(0,self.popsize,2):
                job,machine,machine_time=work_job1[i:i+1],work_M1[i:i+1],work_T1[i:i+1]
                Ma_W1,Tm_W1,WCross=self.to_MT(job,machine,machine_time)
                job1,machine1,machine_time1=work_job1[i+1:i+2],work_M1[i+1:i+2],work_T1[i+1:i+2]
                Ma_W2,Tm_W2,WCross=self.to_MT(job1,machine1,machine_time1)

                MC1,TC1,MC2,TC2=self.mac_cross(Ma_W1,Tm_W1,Ma_W2,Tm_W2,WCross)
                machine_new,time_new=self.back_MT(job,MC1,TC1)
                C_finish,_,_,_,_=to.caculate(job,machine_new,time_new)
                if(C_finish<answer1[i]):      #如果更新后的完工时间大于原解,更新机器和加工时间编码
                    work_M1[i]=machine_new[0]
                    work_T1[i]=time_new[0]
                    answer1[i]=C_finish
                machine_new1,time_new1=self.back_MT(job1,MC2,TC2)
                C_finish,_,_,_,_=to.caculate(job1,machine_new1,time_new1)
                if(C_finish<answer1[i+1]):      #如果更新后的完工时间大于原解,更新机器和加工时间编码
                    work_M1[i+1]=machine_new1[0]
                    work_T1[i+1]=time_new1[0]
                    answer1[i+1]=C_finish
            work_job,work_M,work_T=work_job1,work_M1,work_T1
            answer=answer1
            job_init=job_init1
            result.append([gen+1,min(answer)])#记录每一次迭代的最优个体
            print('灰狼算法第%.0f次迭代的完工时间:%.0f'%(gen+1,min(answer)))
            best_index=answer.tolist().index(min(answer))
            
        return work_job[best_index],work_M[best_index],work_T[best_index],result

ho=pso([10,6,0.5],50,100)    #第一个中括号是工件数,机器数,选择最短机器的概率
#数50,100分别代表迭代的次数和种群的规模

a,b,c,result=ho.gwo_total()  #最后一次迭代的最优解

job,machine,machine_time=np.array([a]),np.array([b]),np.array([c])
to.draw(job,machine,machine_time) #画图

result=np.array(result).reshape(len(result),2)
plt.plot(result[:,0],result[:,1])                   #画完工时间随迭代次数的变化
font1={'weight':'bold','size':22}#汉字字体大小,可以修改
plt.xlabel("迭代次数",font1)
plt.title("完工时间变化图",font1)
plt.ylabel("完工时间",font1)
plt.show()

完整算法源码+数据:见下方微信公众号:关注后回复:车间调度

# 微信公众号:学长带你飞
# 主要更新方向:1、车辆路径问题求解算法
#              2、学术写作技巧
#              3、读书感悟
# @Author  : Jack Hao

公众号二维码:

  • 1
    点赞
  • 12
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
以下是Matlab版本的多目标灰狼算法,可以对车间调度MK01数据集进行多目标求解: ``` function [best_f, best_x] = MOGWO(nPop, MaxIt, lb, ub, nObj, CostFunction) % nPop: 种群大小 % MaxIt: 最大迭代次数 % lb: 变量下界 % ub: 变量上界 % nObj: 目标函数个数 % CostFunction: 目标函数句柄 % 初始化种群 pop = repmat(struct('x', [], 'f', []), nPop, 1); for i = 1:nPop pop(i).x = unifrnd(lb, ub); pop(i).f = feval(CostFunction, pop(i).x); end % 初始化最优解 best_f = zeros(MaxIt, nObj); best_x = zeros(MaxIt, numel(lb)); [best_f(1, :), idx] = min([pop(:).f]); best_x(1, :) = pop(idx).x; % 初始化参数 a = 2; % 收缩参数 A = repmat(struct('alpha', [], 'beta', [], 'delta', []), nPop, 1); for i = 1:nPop A(i).alpha = pop(i).x; A(i).beta = pop(i).x; A(i).delta = pop(i).x; end % 迭代 for it = 2:MaxIt % 更新alpha、beta、delta for i = 1:nPop for j = 1:nPop if pop(j).f(1) < pop(i).f(1) % 更新alpha A(i).alpha = pop(j).x; elseif pop(j).f(1) > pop(i).f(1) && pop(j).f(1) < pop(i).f(1) % 更新beta A(i).beta = pop(j).x; elseif pop(j).f(1) > pop(i).f(1) && pop(j).f(1) > pop(i).f(1) && pop(j).f(1) < pop(i).f(1) % 更新delta A(i).delta = pop(j).x; end end end % 更新位置 for i = 1:nPop r1 = rand(size(lb)); r2 = rand(size(lb)); A(i).alpha = max(min(A(i).alpha + a .* (2 .* r1 - 1) .* abs(A(i).alpha - pop(i).x), ub), lb); A(i).beta = max(min(A(i).beta + a .* (2 .* r2 - 1) .* abs(A(i).beta - pop(i).x), ub), lb); A(i).delta = max(min(A(i).delta + a .* (2 .* rand(size(lb)) - 1) .* abs(A(i).delta - pop(i).x), ub), lb); pop(i).x = (A(i).alpha + A(i).beta + A(i).delta) ./ 3; pop(i).f = feval(CostFunction, pop(i).x); end % 更新最优解 [best_f(it, :), idx] = min([pop(:).f]); best_x(it, :) = pop(idx).x; end end ``` 下面是如何使用该函数求解车间调度MK01数据集的示例: ``` load mk01.mat; % 加载数据集 nPop = 50; % 种群大小 MaxIt = 100; % 最大迭代次数 lb = zeros(1, size(jobs, 1)); % 变量下界 ub = ones(1, size(jobs, 1)); % 变量上界 nObj = 2; % 目标函数个数 CostFunction = @(x) MOGWO_Cost(x, jobs); % 目标函数句柄 [best_f, best_x] = MOGWO(nPop, MaxIt, lb, ub, nObj, CostFunction); % 求解 ``` 其中,`MOGWO_Cost`函数是对MK01数据集进行求解的子函数,代码如下: ``` function [f, g] = MOGWO_Cost(x, jobs) % x: 决策变量向量,每个元素表示一道工序在机器上的执行时间比例 % jobs: 工件信息 % 计算目标函数值 f1 = max(sum(jobs(:, 1:2) .* x', 2)); % 最大完工时间 f2 = sum(sum(jobs(:, 3:4) .* repmat(x', size(jobs, 1), 1))); % 总加权延迟时间 f = [f1, f2]; % 计算约束函数值 g = []; end ``` 其中,`jobs`变量是MK01数据集,每一行表示一道工序的信息,包括所属工件、所在机器、执行时间、加权延迟时间。对于每个决策变量向量`x`,通过计算工序在机器上的执行时间比例,可以得到每个工件的完工时间以及总加权延迟时间,这两个值就是目标函数。由于MK01数据集没有约束条件,因此`g`为空。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值