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

车间调度系列文章：

1、车间调度的编码、解码，调度方案可视化的探讨
2、多目标优化:浅谈pareto寻优和非支配排序遗传算法-NSGAII的非支配排序及拥挤度
3、柔性车间调度问题:以算例MK01初探数据处理和多个遗传算子
4、车间调度丨粒子群算法初探:以算例MK01为例
5、车间调度丨布谷鸟算法改进:以算例MK01为例
6、车间调度丨灰狼算法改进：以算例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 - aC1 = 2 * r2D_alpha =C1*Alpha-xx1 = x - A1 * D_alphar1 = random.random()r2 = random.random()A2 = 2 * a * r1 - aC2 = 2 * r2D_beta =C2*Beta-xx2 = x - A2 * D_betar1 = random.random()r2 = random.random()A3 = 2 * a * r1 - aC3 = 2 * r2D_delta =C3*Delta-xx3 = x - A3 * D_alphainitial_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文档：

篇幅问题，给出灰狼算法的代码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])-1Ma_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,WCrossdef 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])-1m1[0,i]=Ma_W1[signal1][memory1[0,signal1]] #读取对应工序的加工机器t1[0,i]=Tm_W1[signal1][memory1[0,signal1]]memory1[0,signal1]+=1return m1,t1def 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,TC2def gwo_total(self):global tooj=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_aprint('种群初始的完工时间:%.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 - aC1 = 2 * r2D_alpha =C1*Alpha-xx1 = x - A1 * D_alphar1 = random.random()r2 = random.random()A2 = 2 * a * r1 - aC2 = 2 * r2D_beta =C2*Beta-xx2 = x - A2 * D_betar1 = random.random()r2 = random.random()A3 = 2 * a * r1 - aC3 = 2 * r2D_delta =C3*Delta-xx3 = x - A3 * D_alphainitial_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_awork_M1[i],work_T1[i]=machine_new[0],time_new[0]answer1[i]=C_finishfor 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_finishmachine_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_finishwork_job,work_M,work_T=work_job1,work_M1,work_T1answer=answer1job_init=job_init1result.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],resultho=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

公众号二维码：