多目标柔性车间调度丨mogv算法:以算例MK01为例

车间调度系列文章：

1、车间调度的编码、解码，调度方案可视化的探讨
2、多目标优化:浅谈pareto寻优和非支配排序遗传算法-NSGAII的非支配排序及拥挤度
3、柔性车间调度问题:以算例MK01初探数据处理和多个遗传算子
4、车间调度丨粒子群算法初探:以算例MK01为例
5、车间调度丨布谷鸟算法改进:以算例MK01为例
6、车间调度丨自适应灰狼算法改进：以算例MK01为例
7、车间调度丨模拟退火算法改进：以算例MK01为例
8、车间调度入门系列资料
9、多目标柔性车间调度丨改进灰狼算法:以算例MK01为例
10、多目标柔性车间调度丨NSGA-II:以算例MK01为例
11、书本算法重现丨遗传算法：以MK01为例
12、书本算法重现丨元胞粒子群算法：以MK01为例
13、车间调度丨遗传算法求解动态调度问题：重调度
14、柔性车间调度问题丨一种贪婪策略的应用:以算例MK02例
15、柔性车间调度问题丨教学优化算法和均匀交叉算子:以算例MK01例
16、多目标柔性车间调度丨mogv算法:以算例MK01为例

柔性车间调度问题

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

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

MK01算例：

算例的工件数和机器数分别是10和6。

excel的第一行和第一列是编号，不用考虑，修改与否也不影响。

第一行第一个数字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个工件的所有工序的可选机器可加工时间，后面的工序以此类推。

数学模型

符号定义：

n	工件总数	makespani	工件i的完工时间
m	机器总数	makespan	最大完工时间
i,h	工件号	Load_max	机器总负荷
j,k	工序号	Load_all	总能耗
z	机器号	Xijz	工序oij是否在机器z上加工，为0-1变量，在z上加工为1
qi	工件i的工序数	Gijhk	工序oij和工序ohk的先后顺序，为0-1变量，ij在前为1
oij	工件i的第j道工序	M	一个大正实数
Mij	工序oij的可选机器	Tijz	工序oij在机器z的加工时间
Sij	工序oij的开工时间
Cij	工序oij的完工时间
Load_z	机器负荷

模型：

目标函数：

柔性作业车间工具

工序编码

步骤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得到如下编码：

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]

job就是一个可行工序编码。

机器和加工时间编码：

参考文献的3种机器编码生成方法：全局选择、局部选择和随机选择。

对于6个加工机器的mk01。

全局选择：依次安排每个工件的加工，每道工序选择最小负荷的机器。

局部选择：依次安排每个工件的加工，每次安排完一个工件加工后，各个机器的负荷清0，每道工序选择最小负荷的机器。

随机选择：依次安排每个工件的加工，每道工序随机选择可加工机器

核心代码

if r<self.p1 or r>1-self.p2: for k in range(len(n_machine)):m=int(n_machine[k])-1index_select.append(m)t=n_time[k]a_global[0,m]+=t               #全局负荷计算a_part[0,m]+=t                 #局部负荷计算if r<self.p1:                       #全局选择select=a_global[:,index_select]idx_select=np.argmin(select[0])else:                               #局部选择select=a_part[:,index_select]idx_select=np.argmin(select[0])m_select=n_machine[idx_select]t_index=n_machine.index(m_select)machine.append(m_select)machine_time.append(n_time[t_index])
else:                                       #否则随机挑选机器index=np.random.randint(0,len(n_time),1)machine.append(n_machine[index[0]])machine_time.append(n_time[index[0]])

一次随机生成的机器和加工时间编码如下：

machine=[3.0, 2.0, 6.0, 1.0, 3.0, 4.0, 2.0, 3.0, 1.0, 4.0, 1.0, 2.0, 6.0, 1.0, 3.0, 1.0, 1.0, 2.0, 3.0, 5.0, 6.0, 2.0, 1.0, 2.0, 1.0, 4.0, 6.0, 6.0, 1.0, 2.0, 2.0, 1.0, 4.0, 6.0, 4.0, 3.0, 5.0, 3.0, 6.0, 2.0, 1.0, 2.0, 4.0, 6.0, 1.0, 4.0, 1.0, 2.0, 4.0, 6.0, 2.0, 5.0, 6.0, 4.0, 1.0]

machine_time=[4.0, 1.0, 2.0, 1.0, 1.0, 3.0, 6.0, 1.0, 2.0, 6.0, 1.0, 6.0, 2.0, 1.0, 4.0, 1.0, 1.0, 6.0, 1.0, 3.0, 2.0, 1.0, 1.0, 6.0, 5.0, 6.0, 6.0, 2.0, 2.0, 6.0, 6.0, 1.0, 2.0, 1.0, 2.0, 4.0, 1.0, 1.0, 2.0, 6.0, 1.0, 6.0, 6.0, 1.0, 1.0, 3.0, 2.0, 6.0, 6.0, 2.0, 6.0, 3.0, 1.0, 6.0, 3.0]

由算例知道工件1有6道工序，所以machine和machine_time的前6个数表示工件1的6道工序依次在机器3.0、2.0、6.0、 1.0、 3.0、4.0加工，加工时间是4.0、1.0、2.0、 1.0、1.0、 3.0。小数点是因为数据类型是浮点型，不影响，为了美观也可变为整型。

同理工件2有5道工序，所以machine和machine_time的第7到11个数表示工件2的5道工序依次在机器 2.0、3.0、1.0、4.0、 1.0加工，加工时间是 6.0, 1.0, 2.0, 6.0, 1.0，后续编码依次类推。

插入式解码：

参考文献的介绍：

简单来说：每次安排工序的加工机器，首先找对应加工机器的空闲时间，尽量把工序安排空闲时间里。
核心代码：

if jobtime[0,svg] >0 :                            #如果工序不是第一道工序if len(rest[sig])>0 :                         #如果空闲时间非空for m in range(len(rest[sig])-1,-1,-1):   #空闲时间从后往前遍历if rest[sig][m][1]<=jobtime[0,svg] :  #如果某个空闲时间段小于等于上一道工序的完工时间break                             #结束遍历else:                                 #否则begin=max(jobtime[0,svg],rest[sig][m][0])  #可开工时间是上一道工序的完工时间和空闲片段开始时间最大值if begin+machine_time[index] <= rest[sig][m][1] : #如果空闲时间段满足要求startime=begin                #更新开工时间signal=1del rest[sig][m]              #删掉空闲时间段breakif signal==0 :                                    #如果不可插入startime=max(jobtime[0,svg],tmm[0,sig])       #开工时间是加工机器结束时间和上一道工序完工时间的最大值if startime>tmm[0,sig] and signal==0:             #如果不可插入且开工时间大于加工机器的完工时间rest[sig].append([tmm[0,sig],startime])       #添加时间段到空闲时间里
if signal==0 :                                    #如果不可插入tmm[0,sig]=startime+machine_time[index]       #更新机器的结束时间
if signal>0 :                                     #如果可插入signal=0                                      #不更新机器结束时间，且可插入信号归零jobtime[0,svg]=startime+machine_time[index]       #更新工序完工时间
load_m[0,sig]+=machine_time[index]                #更新对应机器的负荷

对本文来说，一次插入成功后，就把当前插入的空闲时间段删除，虽然插入后，时间段仍然可能剩余空闲时间，但认为剩余较短，不考虑。

代码在fjsp.py里。

MOGV算法设计

具体的算法参考文献有介绍，介绍一下步骤：
1、算法步骤：

步骤1：固定比例全局选择、局部选择和随机选择3种方式初始多个工序、机器、加工时间编码，并解码
步骤2：nsga2的方法计算拥挤距离，pareto解，初始pareto解就为记忆库
步骤3：动态概率Pc下从记忆库选择一个个体，种群中轮盘赌选择一个个体，否则1-Pc概率下从种群轮盘赌选取两个个体
步骤4：对选择个体的工序编码进行pox交叉，机器编码进行均匀交叉，产生2个新个体
步骤5：对新种群的个体，在变异概率Pm下对工序编码进行领域搜索变异，对机器编码进行选择最短加工机器变异
步骤6：交叉和变异直到新个体数达到种群规模
步骤7：合并记忆库和新种群，nsga2的方法计算拥挤距离，pareto解，新的pareto解为记忆库
步骤8：判断是否达到最大迭代次数，是的话输出结果，否则转到步骤3

nsga2的方法可以看以前的推文，相比于原文，减少了记忆库的变领域搜索，因为实现比较麻烦且记忆库已经是比较优秀的解了。

2、轮盘赌选择

轮盘赌法是模拟博彩游戏的轮盘赌，扇形的面积对应它所表示的染色体的适应值的大小，适应度值越大个体被选择的可能性也就越大。轮盘赌法的关键部分是概率和累计概率的计算，具体的步骤如下：

python里实现比较容易，代码如下：

def select(self,num,fit):fit=np.array(fit)idx=np.random.choice(np.arange(self.popsize),size=num,replace=False,p=fit/fit.sum())return idx

返回的是轮盘赌选择出的位置索引，其意思是self.popsize个体中选出num个个体，每个个体的选中概率是fit/fit.sum()，replace=False表示不重复。

本文的fit用的是完工时间分之一，完工时间短，适应度越高，当然，对于多目标问题，可以自行设计适应度计算方式，可以多个目标自行组合等等。

3、工序的pox交叉
以mk01为例：随机0到9的一个数为6，对应两个进行交叉的工序编码，0到6基因及其位置保持不变，每个编码6到9的基因位置按顺序填入另一个工序编码6到9的基因。

核心代码：

def job_cross(self,chrom_L1,chrom_L2):       #工序的pox交叉num=list(set(chrom_L1[0]))np.random.shuffle(num)index=np.random.randint(0,len(num),1)[0]jpb_set1=num[:index+1]                  #固定不变的工件jpb_set2=num[index+1:]                  #按顺序读取的工件C1,C2=np.zeros((1,chrom_L1.shape[1]))-1,np.zeros((1,chrom_L1.shape[1]))-1sig,svg=[],[]for i in range(chrom_L1.shape[1]):#固定位置的工序不变ii,iii=0,0for j in range(len(jpb_set1)):if(chrom_L1[0,i]==jpb_set1[j]):C1[0,i]=chrom_L1[0,i]else:ii+=1if(chrom_L2[0,i]==jpb_set1[j]):C2[0,i]=chrom_L2[0,i]else:iii+=1if(ii==len(jpb_set1)):sig.append(chrom_L1[0,i])if(iii==len(jpb_set1)):svg.append(chrom_L2[0,i])signal1,signal2=0,0             #为-1的地方按顺序添加工序编码for i in range(chrom_L1.shape[1]):if(C1[0,i]==-1):C1[0,i]=svg[signal1]signal1+=1if(C2[0,i]==-1):C2[0,i]=sig[signal2]signal2+=1return C1,C2

4、机器的均匀交叉

均匀交叉算子的概念比较简单，简单说一下逻辑：假设两个解的工序编码的第一道工序分别选择了机器1和3，随机生成0，1两个数中的一个，如果随机数是1，交换两个解第一道工序的机器选择，否则保持原选择。

代码：

def ma_cross(self,m1,t1,m2,t2):  #机器均匀交叉MC1,MC2,TC1,TC2=[],[],[],[]for i in range(m1.shape[0]):     index=np.random.randint(0,2,1)[0]if(index==0):  #为0时继承继承父代的机器选择MC1.append(m1[i]),MC2.append(m2[i]),TC1.append(t1[i]),TC2.append(t2[i]);else:                #为1时另一个父代的加工机器选择MC2.append(m1[i]),MC1.append(m2[i]),TC2.append(t1[i]),TC1.append(t2[i]);return MC1,TC1,MC2,TC2

5、领域搜索变异

参考文献的算法步骤和示意图如下：

步骤1：在变异染色体中随机选择ｒ个不同基因，并生成其排序的所有邻域。

步骤２评价所有邻域的适应值，选出最佳个体作为子代。

可以看出，邻域搜索就是对r个基因进行排队，邻域的个数是r!-1,查阅资料，找到一种递归的方法实现这一过程：固定位置的基因不变放到第一个，依次降低r个基因的个数，降低到固定位置的基因和变动位置的基因只有1个时，可变位置与不变位置基因合并。

以上面的[3,1,2]进行介绍：

位置1的基因不变放到第一位，剩余[1,2],

[1,2]只有两个，分别固定第一个位置和第二个位置基因不变，合并得到[1,2]和[2,1]，然后与外层的基因合并得到[3,1,2],和[3,2,1]

[3,1,2]其他位置的邻域也按照这种方法生成：

代码：

def ways(self,mul):if len(mul)<2:                          #基因个数降低为1取原值return [mul]else :                                  #否则遍历基因result=[]for i in range(len(mul)) :   select=mul[i]                   #不变位置last=mul[:i]+mul[i+1:]          #可变位置for mu in self.ways(last):      #原函数递归，减少基因个数result.append([select]+mu)  #合并可变与不变位置的基因return result[1:]                       #返回除原位置的所有邻域

输入[3,1,2],结果是[[3, 2, 1], [1, 3, 2], [1, 2, 3], [2, 3, 1], [2, 1, 3]]

本文设计整个邻域搜索变异的逻辑：新解没有完全劣于原解就更新编码

代码：

def job_mul(self,w,m,t):r=np.random.randint(2,5,1)[0]                   #2到5随机生成一个数mul=random.sample(range(w.shape[0]),r)          #0到55生成不重复的r个位置C_finish,load_max,load_all,_,_,_,_=self.to.caculate(w,m,t) t1=[C_finish,load_max,load_all]                  #原编码的目标计算res=np.arange(r).tolist()result=self.ways(res)mul=np.array(mul)for re in result :                              #位置索引w[mul]=w[mul[re]]                           #更新编码C_finish,load_max,load_all,_,_,_,_=self.to.caculate(w,m,t)t2=[C_finish,load_max,load_all]          #新编码的目标计算if (np.array(t2)>np.array(t1)).all() :   #如果新编码的3个目标都劣于原目标，不更新编码w[mul[re]]=w[mul]else:                                    #否则，更新编码，并更新目标t1=t2return w,m,t

领域搜索的计算量是巨大的，相当于穷举r个位置的所有排列组合，本文暂时先设为r设为2到5之间，可以自行调整。

6、最短加工机器变异

比较简单，就是工序选择最短加工时间的机器，当然，如果原工序已经选择了最大短加工时间的机器，该变异不会对编码产生改变。

代码：

def ma_mul(self,w,m,t):r=np.random.randint(1,w.shape[0]+1,1)[0]  #变异的数量mul=random.sample(range(w.shape[0]),r)    #变异的位置count=np.zeros((1,self.job_num),dtype=np.int)signal=0sig=0for i in range(len(self.machines)):for j in range(self.machines[i]):if signal<len(mul) and sig == mul[signal] :highs=self.tom[sig][count[0,sig]]lows=self.tom[sig][count[0,sig]]-self.tdx[sig][count[0,sig]]n_machine=self.Tmachine[sig,lows:highs].tolist()n_time=self.Tmachinetime[sig,lows:highs].tolist()time=min(n_time)         #最短加工时间的index=n_time.index(time)即mac=n_machine[index]     #最短加工时间的机器m[sig]=mac               #更新编码t[sig]=timesignal+=1sig+=1return w,m,t

整个算法迭代的核心代码：

Pc=1-gen/self.generation                          #动态交叉概率
fit=[1/answer[k][0] for k in range(len(answer))]  #完工时间分之1做适应度for i in range(0,self.popsize,2):                #种群规模下每次选2个个体if np.random.rand()<Pc:loc1=np.random.randint(len(pareto))      #记忆库pareto随机挑选1个个体loc2=self.select(1,fit)[0]               #迭代种群中轮盘赌选1个个体W1,M1,T1=pareto_job[loc1],pareto_machine[loc1],pareto_time[loc1]W2,M2,T2=work_job[loc2],work_M[loc2],work_T[loc2]else:idx=self.select(2,fit)                   #否则轮盘赌选择2个个体loc1,loc2=idx[0],idx[1]W1,M1,T1=work_job[loc1],work_M[loc1],work_T[loc1]W2,M2,T2=work_job[loc2],work_M[loc2],work_T[loc2]W1,W2=self.job_cross(W1,W2)                 #工序交叉M1,T1,M2,T2=self.ma_cross(M1,T1,M2,T2)      #机器交叉if np.random.rand()<self.Pm :              #变异概率W1,M1,T1=self.ma_mul(W1,M1,T1)         #机器变异W1,M1,T1=self.job_mul(W1,M1,T1)        #工序变异if np.random.rand()<self.Pm :W2,M2,T2=self.ma_mul(W2,M2,T2)W2,M2,T2=self.job_mul(W2,M2,T2)

以上所有代码放在MOGV.py里。

7、pareto结果去重

pareto结果会出现重复的目标，可能是因为多个解对应相同的目标，也可能是nsga2的算子设计有问题，或者其他原因。原因暂且不论，为了美观和简洁，对pareto最后的解进行去重，当然，理论上会舍弃一些解，不需要的可以不用这一步。

去重很简单，相当于用新的列表存储原pareto的解，初始新列表为空

每一次遍历原pareto解的元素，如果新列表不存在该元素就添加，否则不添加

代码：

def pareto_simplyfy(self,pareto):pareto_un=[]index=[]for i in range(len(pareto)):if pareto[i] not in pareto_un:  #新列表不存在原pareto的元素就添加机器pareto_un.append(pareto[i])index.append(i)             #记录位置return pareto_un,index

代码在multi_opt.py里

结果

代码运行环境
windows系统，python3.6.0,第三方库及版本号如下：

numpy==1.18.5
matplotlib==3.2.1

第三方库需要在安装完python之后，额外安装，以前文章→点击我跳转有讲述过安装第三方库的解决办法。

命令


da=data_deal(10,6)                   #工件数，机器数
Tmachine,Tmachinetime,tdx,work,tom,machines=da.cacu()
parm_data=[Tmachine,Tmachinetime,tdx,work,tom,machines]
fj=FJSP(10,6,0.3,0.4,parm_data)      #工件数，机器数，全局与局部选择的概率和mk01的数据，第三种选择概率是1-0.3-0.4
mu=mul_op()                          #多目标模块                 parm_mo=[50,100,1]                    #mogv算法的参数，依次是迭代次数，种群规模，变异概率
mod=[fj,mu]                          #柔性模块和多目标模块
mo=mogv(10,parm_mo,mod,parm_data)    #工件数，算法参数，功能模块，mk01数据
pareto,pareto_job,pareto_machine,pareto_time,fit_every=mo.mog()
print('\n原pareto解：\n',pareto)                        #原pareto结果pareto_un,index=mu.pareto_simplyfy(pareto)
print('\n去重pareto解：\n',pareto_un)                     #去重pareto结果
job_un,machine_un,time_un=pareto_job[index],pareto_machine[index],pareto_time[index]mu.draw_change(fit_every)           #每次迭代过程中pareto解中3个目标的变化
mu.draw_3d(pareto)                  #3维图              sig=0                               #取pareto的第一个解
job,machine,machine_time=job_un[sig],machine_un[sig],time_un[sig]C_finish,load_max,load_all,list_M,list_S,list_W,tmax=fj.caculate(job,machine,machine_time)
fj.draw(job,machine,machine_time)  #画甘特图

运行结果

一次迭代结果如下：

3个目标的最大、平均、最小值随迭代次数的变化图如下：

3维图：

一个解的甘特图如下：

从结果上看，相同算例下，已经和之前文章已经介绍过的灰狼算法持平，是比较优秀的算法了。

结论

本文的工作量较大，涉及的内容较多，有编码、解码、多目标优化、nsga2等等。
算法大部分完全是复现原论文，也有自己的设计，文章的算法也有比较详细的介绍，插入式解码和邻域搜索等一些方法是比较不错的，希望对大家的论文有一些参考价值。
excel数据可更改，工件数、机器数、工件的工序数、工序的可加工机器数等数据对得上就能运行。
参考文献：柔性作业车间调度智能算法及其应用-高亮（第七章)

源代码

有5个py文件和一个mk01的excel文件：

篇幅问题，代码附在后面。

演示视频：

多目标求解柔性车间调度问题丨MOGV算法

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

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

公众号二维码：