车间调度丨粒子群算法初探:以算例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_num
        self.machine_num=machine_num
    def read(self):
        f=open('./MK01.txt')
        f1=f.readlines()
        c,count=[],0
        for 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+=1
        return c
    def translate(self,tr1):
        sigdex,mac,mact,sdx=[],[],[],[]
        sigal=tr1[0]
        tr1=tr1[1:len(tr1)+1]
        index=0
        for j in range(sigal):
            sig=tr1[index]
            sdx.append(sig)
            sigdex.append(index)
            index=index+1+2*sig
        for 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,sdx
    def 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 width
    def 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]=mac
            Tmachinetime[i,0:siga]=mact
        return Tmachine,Tmachinetime,tdx
    def 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]+=1
                if 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_a
    def 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=0
            list_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,tmax
    def 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]+=1
                    plt.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_a
                answer[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 random

class 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])-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 pso_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_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_a
                    pbest[i]=initial_a
                best_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_a
                                answer[i]=C_finish
                                pbest[i]=initial_a
                        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
    
                machine_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_finish
            best_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左右，也没用深究其原因了。但既然初略地算法写出来了，就分享一下。

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