车间调度丨粒子群算法初探:以算例MK01为例

车间调度系列文章：

1、车间调度的编码、解码，调度方案可视化的探讨
2、多目标优化:浅谈pareto寻优和非支配排序遗传算法-NSGAII的非支配排序及拥挤度
3、柔性车间调度问题:以算例MK01初探数据处理和多个遗传算子
4、车间调度丨粒子群算法初探:以算例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个工件的所有工序的可选机器可加工时间，后面的工序以此类推。

下面的每一行以此类推。

具体的数据处理上一篇推文有介绍，数据格式是txt,处理的代码如下：

import numpy as np class data_deal:def __init__(self,job_num,machine_num):self.job_num=job_numself.machine_num=machine_numdef read(self):f=open('./MK01.txt')f1=f.readlines()c,count=[],0for line in f1:t1=line.strip('\n')if(count>0):cc=[]for j in range(len(t1)):if(t1[j]!=' '):cc.append(int(t1[j]))c.append(cc)count+=1return cdef translate(self,tr1):sigdex,mac,mact,sdx=[],[],[],[]sigal=tr1[0]tr1=tr1[1:len(tr1)+1]index=0for j in range(sigal):sig=tr1[index]sdx.append(sig)sigdex.append(index)index=index+1+2*sigfor ij in range(sigal):del tr1[sigdex[ij]-ij]for ii in range(0,len(tr1)-1,2):mac.append(tr1[ii])mact.append(tr1[ii+1])return mac,mact,sdxdef widthxx(self,strt):widthx=[]for i in range(self.job_num):mac,mact,sdx=self.translate(strt[i])siga=len(mac)widthx.append(siga)width=max(widthx)return widthdef tcaculate(self,strt):width=self.widthxx(strt)Tmachine,Tmachinetime=np.zeros((self.job_num,width)),np.zeros((self.job_num,width))tdx=[]for i in range(self.job_num):mac,mact,sdx=self.translate(strt[i])tdx.append(sdx)siga=len(mac)Tmachine[i,0:siga]=macTmachinetime[i,0:siga]=mactreturn Tmachine,Tmachinetime,tdxdef cacu(self):strt=self.read()Tmachine,Tmachinetime,tdx=self.tcaculate(strt)to,tom,work=0,[],[]for i in range(self.job_num):to+=len(tdx[i])tim=[]for j in range(1,len(tdx[i])+1,1):tim.append(sum(tdx[i][0:j]))work.append(i)tom.append(tim)return Tmachine,Tmachinetime,tdx,work,tom

粒子群算法简介：

网上找的介绍：粒子群优化算法（Particle Swarm Optimization，PSO）属于进化算法的一种，是通过模拟鸟群捕食行为设计的。从随机解出发，通过迭代寻找最优解，通过适应度来评价解的品质。PSO 初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个"极值"来更新自己。第一个就是粒子本身所找到的最优解pbest。另一个极值是整个种群目前找到的最优解，即全局极值gbest。
所有的粒子具有以下两个属性。速度、位置，更新公式如下：
速度：

vi+1=w∗vi*+c1∗ran**d1∗(pbest**i−x**i)+c2∗ran**d2∗(gbest**i−*x**i)

位置：
xi+1=xi+vi+1

w为惯性因子，一般取1；c1、c2为学习因子，一般取2；rand*1、*rand2为两个（0,1）之间的随机数；vi和xi分别表示粒子第i维的速度和位置；pbesti、gbesti分别表示某个粒子最好位置第i维的值、整个种群最好位置第i维的值。

算法求解流程：

本文写了甘特图的画图函数；工序，机器，加工时间编码的生成函数；编码的解码函数。甘特图和解码前面推文有介绍，为了能在粒子群算法使用，下面简述一下编码的生成：

步骤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到4倍编码长度生成55个不重复的数，即0到220生成55个不重复的数。

initial_a= [138, 200, 170, 51, 2, 119, 9, 0, 82, 81, 165, 39, 174, 71, 11, 74, 24, 107, 23, 204, 15, 33, 162, 50, 37, 31, 89, 148, 97, 43, 106, 25, 136, 61, 190, 45, 178, 145, 93, 149, 209, 188, 196, 159, 157, 137, 155, 110, 100, 187, 48, 173, 122, 69, 88]

步骤3：对步骤2得到的编码，进行从小到大的排列，按照从小到大的顺序取出其在编码的位置。

index_work= [7, 4, 6, 14, 20, 18, 16, 31, 25, 21, 24, 11, 29, 35, 50, 23, 3, 33, 53, 13, 15, 9, 8, 54, 26, 38, 28, 48, 30, 17, 47, 5, 52, 32, 45, 0, 37, 27, 39, 46, 44, 43, 22, 10, 2, 51, 12, 36, 49, 41, 34, 42, 1, 19, 40]

initial_a的第8（7+1）个数最小，第5（4+1）个数第二小，依次类推。

步骤4：依据步骤3的得到编码位置，在work中找到对应工序编码，即可得到长度55的工序编码。

job= [1, 0, 1, 2, 3, 3, 3, 5, 4, 4, 4, 2, 5, 6, 9, 4, 0, 6, 9, 2, 2, 1, 1, 9, 4, 7, 5, 8, 5, 3, 8, 0, 9, 5, 8, 0, 6, 5, 7, 8, 8, 8, 4, 1, 0, 9, 2, 6, 9, 7, 6, 7, 0, 3, 7]

work的第8（7+1）个数是1，第5（4+1）个数是0，依次类推。

代码在fjsp里。

机器和加工时间编码：得到工序编码后，依据工序编码取加工机器和加工时间，如上的job的前两个工序是工件2和工件1的第一道工序，对应的可选机器是[2]和[1，3]，加工时间是[6]和[5，4]加工机器上面的表。

工件2的第一道工序的可选机器只有2，只有一个选择，不做讨论，工件1的第一道工序的可选机器是1和3，有两个选择，设计选择方法如下：

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

整个柔性车间作业车间代码如下：

class FJSP():def __init__(self,job_num,machine_num,pi,parm_data):self.job_num=job_num                #工件数self.machine_num=machine_num        #机器数self.pi=pi                  #随机挑选机器的概率self.Tmachine,self.Tmachinetime,self.tdx,self.work,self.tom=parm_data[0],parm_data[1],parm_data[2],parm_data[3],parm_data[4]def axis(self):index=['M1','M2','M3','M4','M5','M6','M7','M8','M9','M10','M11','M12','M13','M14','M15','M16','M17','M18','M19','M20']scale_ls,index_ls=[],[]   for i in range(self.machine_num):scale_ls.append(i+1)index_ls.append(index[i])return index_ls,scale_ls  #返回坐标轴信息，按照工件数返回，最多画20个机器，需要在后面添加def creat_job(self):initial_a=random.sample(range(4*len(self.work)),len(self.work))index_work=np.array(initial_a).argsort()job=[]for i in range(len(self.work)):job.append(self.work[index_work[i]])job=np.array(job).reshape(1,len(self.work))ccount=np.zeros((1,self.job_num),dtype=np.int)machine=np.ones((1,job.shape[1]))machine_time=np.ones((1,job.shape[1]))    #初始化矩阵for i in range(job.shape[1]):oper=int(job[0,i])highs=self.tom[oper][ccount[0,oper]]lows=self.tom[oper][ccount[0,oper]]-self.tdx[oper][ccount[0,oper]]n_machine=self.Tmachine[oper,lows:highs]n_time=self.Tmachinetime[oper,lows:highs]ccount[0,oper]+=1if np.random.rand()>self.pi:                #选取最小加工时间机器     machine_time[0,i]=min(n_time)index=np.argwhere(n_time==machine_time[0,i])machine[0,i]=n_machine[index[0,0]]else:                                       #否则随机挑选机器                                index=np.random.randint(0,len(n_time),1)machine[0,i]=n_machine[index[0]]machine_time[0,i]=n_time[index[0]]return job,machine,machine_time,initial_adef caculate(self,job,machine,machine_time):jobtime=np.zeros((1,self.job_num))        tmm=np.zeros((1,self.machine_num))              tmmw=np.zeros((1,self.machine_num))         startime=0list_M,list_S,list_W=[],[],[]for i in range(job.shape[1]):svg,sig=int(job[0,i]),int(machine[0,i])-1  if(jobtime[0,svg]>0):                               startime=max(jobtime[0,svg],tmm[0,sig])     tmm[0,sig]=startime+machine_time[0,i]jobtime[0,svg]=startime+machine_time[0,i]if(jobtime[0,svg]==0):                          startime=tmm[0,sig]tmm[0,sig]=startime+machine_time[0,i]jobtime[0,svg]=startime+machine_time[0,i]tmmw[0,sig]+=machine_time[0,i]list_M.append(machine[0,i])list_S.append(startime)list_W.append(machine_time[0,i])tmax=np.argmax(tmm[0])+1        #结束最晚的机器C_finish=max(tmm[0])            #最晚完工时间return C_finish,list_M,list_S,list_W,tmaxdef draw(self,job,machine,machine_time):#画图C_finish,list_M,list_S,list_W,tmax=self.caculate(job,machine,machine_time)    figure,ax=plt.subplots()count=np.zeros((1,self.job_num))for i in range(job.shape[1]):  #每一道工序画一个小框count[0,int(job[0,i])-1]+=1plt.bar(x=list_S[i], bottom=list_M[i], height=0.5, width=list_W[i], orientation="horizontal",color='white',edgecolor='black')plt.text(list_S[i]+list_W[i]/32,list_M[i], '%.0f' % (job[0,i]+1),color='black',fontsize=10,weight='bold')#12是矩形框里字体的大小，可修改plt.plot([C_finish,C_finish],[0,tmax],c='black',linestyle='-.',label='完工时间=%.1f'% (C_finish))#用虚线画出最晚完工时间font1={'weight':'bold','size':22}#汉字字体大小，可以修改plt.xlabel("加工时间",font1)plt.title("甘特图",font1)plt.ylabel("机器",font1)scale_ls,index_ls=self.axis()plt.yticks(index_ls,scale_ls)plt.axis([0,C_finish*1.1,0,self.machine_num+1])plt.tick_params(labelsize = 22)#坐标轴刻度字体大小，可以修改labels=ax.get_xticklabels()[label.set_fontname('Times New Roman')for label in labels]plt.legend(prop={'family' : ['STSong'], 'size'   : 16})#标签字体大小，可以修改plt.xlabel("加工时间",font1)plt.show()

粒子群算法数据导入

本文对工序编码用粒子群更新，局部最优pbesti用第i个解目前找到的最优编码代替，全局最优gbesti是每次迭代种群的最优解。机器编码用前面介绍的机器编码交叉。

工序编码迭代的核心代码如下：

 for i in range(self.popsize):job,machine,machine_time=work_job[i:i+1],work_M[i:i+1],work_T[i:i+1]Ma_W1,Tm_W1,WCross=self.to_MT(job,machine,machine_time)x=job_initial[i]v[i] = self.W *v[i]  + self.C1 * random.random() * (pbest[i] - x) + self.C2 * random.random() * (gbest - x)initial_a=x+v[i]index_work=initial_a.argsort()job=[]for j in range(len(work)):job.append(work[index_work[j]])job=np.array(job).reshape(1,len(work))machine_new,time_new=self.back_MT(job,Ma_W1,Tm_W1)C_finish,_,_,_,_=to.caculate(job,machine_new,time_new)work_job[i]=job[0]job_initial[i]=initial_aanswer[i]=C_finish

简单来说，每次更新位置initial_a后，在用work读取工序编码。

机器编码的核心代码如下：

for i in range(0,self.popsize,2):job,machine,machine_time=work_job[i:i+1],work_M[i:i+1],work_T[i:i+1]Ma_W1,Tm_W1,WCross=self.to_MT(job,machine,machine_time)job1,machine1,machine_time1=work_job[i+1:i+2],work_M[i+1:i+2],work_T[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<answer[i]):      #如果更新后的完工时间低于原解，更新机器和加工时间编码work_M[i]=machine_new[0]work_T[i]=time_new[0]answer[i]=C_finish

简单来说，每次机器编码更新后，如果完工时间低于原编码，更新原编码，否则不更新。所有代码如下：

from data_solve import data_deal
from fjsp import FJSP
import numpy as np
import randomclass pso():def __init__(self,parm_fjsp,generation,popsize,parm_pso):self.job_num=parm_fjsp[0]                       #工件数self.machine_num=parm_fjsp[1]               #机器数self.pi=parm_fjsp[2]self.generation=generation          #迭代次数self.popsize = popsize  # 粒子个数self.W = parm_pso[0]self.C1 = parm_pso[1]self.C2 = parm_pso[2]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 pso_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_initial,pbest=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)))v=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_initial[i]=initial_apbest[i]=initial_abest_index=answer.index(min(answer))gbest=pbest[best_index]for i in range(self.popsize):job,machine,machine_time=work_job[i:i+1],work_M[i:i+1],work_T[i:i+1]Ma_W1,Tm_W1,WCross=self.to_MT(job,machine,machine_time)x=job_initial[i]v[i] = self.W *v[i]  + self.C1 * random.random() * (pbest[i] - x) + self.C2 * random.random() * (gbest - x)initial_a=x+v[i]index_work=initial_a.argsort()job=[]for j in range(len(work)):job.append(work[index_work[j]])job=np.array(job).reshape(1,len(work))machine_new,time_new=self.back_MT(job,Ma_W1,Tm_W1)C_finish,_,_,_,_=to.caculate(job,machine_new,time_new)work_job[i]=job[0]job_initial[i]=initial_aanswer[i]=C_finishpbest[i]=initial_afor i in range(0,self.popsize,2):job,machine,machine_time=work_job[i:i+1],work_M[i:i+1],work_T[i:i+1]Ma_W1,Tm_W1,WCross=self.to_MT(job,machine,machine_time)job1,machine1,machine_time1=work_job[i+1:i+2],work_M[i+1:i+2],work_T[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<answer[i]):      #如果更新后的完工时间大于原解，更新机器和加工时间编码work_M[i]=machine_new[0]work_T[i]=time_new[0]answer[i]=C_finishmachine_new1,time_new1=self.back_MT(job1,MC2,TC2)C_finish,_,_,_,_=to.caculate(job1,machine_new1,time_new1)if(C_finish<answer[i+1]):      #如果更新后的完工时间大于原解，更新机器和加工时间编码work_M[i+1]=machine_new1[0]work_T[i+1]=time_new1[0]answer[i+1]=C_finishbest_index=answer.index(min(answer))gbest=job_initial[best_index]result.append(answer[best_index])print(answer[best_index])return work_job[best_index],work_M[best_index],work_T[best_index],result

结果分析：

运行命令如下：

ho=pso([10,6,0.5],2,20,[1,2,2])    #第一个中括号是工件数，机器数，选择最短机器的概率
#第一个中括号后面的2,20分别代表迭代的次数和粒子的数目
#第二个是粒子群的参数，1,2，2分别是w,c1,c2。
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) #画图

结果的甘特图如下：

结论

本文就MK01，介绍了粒子群算法对其求解的方法。但粒子群算法刚刚涉猎，可能是算法参数设置不对，或者算法设计不合理，也或者是该问题不太适合粒子群算法，求解结果不是太好。MK01现在的最优结果一般是40及以下低一点，但本文的结果只能到达50左右，也没用深究其原因了。但既然初略地算法写出来了，就分享一下。

完整代码，可见下面微信公众号,扫描关注后，回复：路径优化。