多目标柔性车间调度丨NSGA-II:以算例MK01为例

车间调度系列文章：

1、车间调度的编码、解码，调度方案可视化的探讨
2、多目标优化:浅谈pareto寻优和非支配排序遗传算法-NSGAII的非支配排序及拥挤度
3、柔性车间调度问题:以算例MK01初探数据处理和多个遗传算子
4、车间调度丨粒子群算法初探:以算例MK01为例
5、车间调度丨布谷鸟算法改进:以算例MK01为例
6、车间调度丨自适应灰狼算法改进：以算例MK01为例
7、车间调度丨模拟退火算法改进：以算例MK01为例
8、车间调度入门系列资料
9、多目标柔性车间调度丨改进灰狼算法:以算例MK01为例
10、多目标柔性车间调度丨NSGA-II:以算例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个工件的所有工序的可选机器可加工时间，后面的工序以此类推。

下面的每一行以此类推。本系列算例MK01.txt数据文件可在gitee仓库下载：

https://gitee.com/XZDNF-1618/data.git

数学模型

假设条件：

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

数学模型

字符说明：

目标函数：

约束条件：

柔性作业车间工具

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

步骤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就是一个可行工序编码。

代码在fjsp.py里。

机器和加工时间编码：

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

非支配排序遗传算法-NSGA-II的非支配排序及拥挤度

就当前很多组合优化问题，如车间调度，路径优化等问题。运用非精确算法如粒子群算法、遗传算法、蚁群算法解决问题时，我们总需要经过多次迭代才能找到问题的解。每一次迭代时，总会挑选一个或多个优秀的个体指导着目标优化的方向。对于单目标问题，个体的优劣判断比较容易，即目标值的大小。而多目标问题，因为涉及多个目标，很难判断个体的优劣。就遗传算法解决多目标车间调度问题来说，初始化一定规模的种群，但因为我们无法判断个体的优劣，适应度就无法设计，算法最开始的选择操作都无法进行。当然，pareto解肯定是种群的优秀个体，但其规模，分布情况等，只选择其作为向下迭代的种群是不可行的。

非支配排序遗传算法-NSGAII是2000年提出的，在遗传算法的基础上对选择再生方法进行改进：将每个个体按照他们的支配与非支配关系进行再分层，提出了拥挤度和拥挤度比较算子，代替了需要指定共享半径的适应度共享策略，并在快速排序后的同级比较中作为胜出标准，使准Pareto域中的个体能扩展到整个Pareto域，并均匀分布，保持了种群的多样性。最后再做遗传操作，从而达到求解多目标的问题。

简单来说，就是对多个解按照支配关系进行分层，每层又按照拥挤度大小进行排列，这样就能判断每个个体的优劣，其核心是非支配型排序和拥挤度计算。

非支配型排序方法

n(i) 为在种群中支配个体 i 的解个体的数量。S(i) 为被个体 i 所支配的解个体的集合。

首先，找到种群中所有 n(i)=0 的个体，将它们存入当前集合F(1)；
然后对于当前集合 F(1) 中的每个个体 j，考察它所支配的个体集 S(j)，将集合 S(j) 中的每个个体 k 的 n(k) 减去1，即支配个体 k 的解个体数减1(因为支配个体 k 的个体 j 已经存入当前集 F(1) )；
如 n(k)-1=0则将个体 k 存入另一个集H。最后，将 F(1) 作为第一级非支配个体集合，并赋予该集合内个体一个相同的非支配序 i(rank)，然后继续对 H 作上述分级操作并赋予相应的非支配序，直到所有的个体都被分级。其计算复杂度为O(mN^{2})，m为目标函数个数，N为种群大小。

简单来说，假设有两个解p,q，如果p支配q，则q放到p支配的集合S§ 即里，且n(q)在原来的等级基础上上升一个等级。

分层的逻辑:层级为0的个体是最优的一层，把层级为0的个体取出来之后，被该层支配的所有个体等级减1，以此类推，不断取出等级为0的个体，最终所有剩余个体的等级都会为0，分层结束。

拥挤度计算

密度估计：对同一个前沿面的解集合按各个目标分量大小排序，计算每个解在该分量下的两侧点的距离差值，而后进行累加各个分量上的距离作为拥挤系数。

经过快速非支配排序和拥挤度计算，种群中的每一个个体都得到了两个属性：非支配序rank和拥挤度crowd。利用这两个属性，我们可以区分种群中间任意两个个体间的支配和非支配关系。定义拥挤度比较算子，当且仅当i_rank>j_rank或i_rank=j_rank且i_crowd>j_crowd，有个体i优于个体j。为了简化，直接把个体优劣的比较写进拥挤度程序里。为了更好排序，把每层两端的点的拥挤度设为100000和100001，一层只有一个个体时，因为不存在同级的比较，拥挤度设为1。

具体的介绍自行查阅资料，前面的推文也讲到。

非支配排序遗传算法设计

简单来说：每次迭代以上一次迭代得到种群为基础，工序编码pox交叉后根据拥挤度比较淘汰一半个体，机器编码均匀交叉根据拥挤度淘汰一半个体，第一次迭代以初始种群为基础。
算法步骤：

步骤1：初始化多个工序，机器，加工时间编码，解码得到加工时间、机器负荷和能耗
步骤2：工序编码进行pox交叉，根据加工时间、机器负荷和能耗，进行非支配排序和拥挤度的计算，淘汰一半个体
步骤3：机器编码进行均匀交叉，根据加工时间、机器负荷和能耗，进行非支配排序和拥挤度的计算，淘汰一半个体
步骤4：判断是否达到最大迭代次数，是的话输出结果，否则转到步骤2

pox交叉
以mk01为例：初始工件编号为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

均匀交叉

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

核心代码：

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

具体的代码细节在代码在nsga_2.py里。

非支配排序和拥挤度

具体的介绍自行查阅资料，前面的推文也讲到。这里贴一下代码，在mult_opt.py里。

代码：

class mul_op():def divide(self,answer):S=[[] for i in range(len(answer))]front = [[]]n=[0 for i in range(len(answer))]rank = [0 for i in range(len(answer))]for p in range(len(answer)):for q in range(len(answer)):# 如果p支配qif (np.array(answer[p])<=np.array(answer[q])).all() and (answer[p]!=answer[q]): if q not in S[p]:S[p].append(q)  # 同时如果q不属于sp将其添加到sp中# 如果q支配pelif (np.array(answer[p])>=np.array(answer[q])).all() and (answer[p]!=answer[q]):n[p] = n[p] + 1  # 则将np+1if n[p]==0:rank[p] = 0if p not in front[0]:front[0].append(p)i = 0while(front[i] != []):Q=[]for p in front[i]:for q in S[p]:n[q] =n[q] - 1  # 则将fk中所有给对应的个体np-1if( n[q]==0):   # 如果nq==0rank[q]=i+1if q not in Q:Q.append(q)i = i+1front.append(Q)del front[len(front)-1]return frontdef dis(self,answer):crowder,crowd=[],[]front=self.divide(answer)for i in range(len(front)):x=[answer[front[i][j]][0] for j in range(len(front[i]))] #取三个目标函数的各个目标y=[answer[front[i][j]][1] for j in range(len(front[i]))]z=[answer[front[i][j]][2] for j in range(len(front[i]))]sig=front[i]clo=[[] for j in range(len(front[i]))]if(len(sig)>1):    #每层的个体大于1个做拥挤度计算x_index,y_index,z_index=np.array(x).argsort(),np.array(y).argsort(),np.array(z).argsort()x.sort(),y.sort();z.sort()dis1,dis2,dis3=[],[],[]dis1.append(100000);dis2.append(100000);dis3.append(100000)if(len(sig)>2):    #大于2个做中间个体的拥挤度计算for k in range(1,len(sig)-1):distance1,distance2,distance3=(x[k+1]-x[k-1])/(x[-1]-x[0]),(y[k+1]-y[k-1])/(y[-1]-y[0]),(z[k+1]-z[k-1])/(z[-1]-z[0])dis1.append(distance1);dis2.append(distance2);dis3.append(distance3)dis1.append(100001);dis2.append(100001);dis3.append(100001)crow=[]x_index=x_index.tolist()y_index=y_index.tolist()z_index=z_index.tolist()for m in range(len(sig)):index1,index2,index3=x_index.index(m),y_index.index(m),z_index.index(m)cro=dis1[index1]+dis2[index2]+dis3[index3]crow.append(cro)crowd.append(crow)index=np.array(crow).argsort()for n in range(len(index)):     #拥挤度排列并取出clo[n]=sig[index[n]]for o in range(len(clo)-1,-1,-1):crowder.append(clo[o])else:crowder.append(front[i][0])crowd.append([1])return front,crowd,crowder

结果

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

numpy==1.18.5
matplotlib==3.2.1

第三方库需要在安装完python之后，额外安装，以前文章有讲述过安装第三方库的解决办法。
功率设计

机器的负载和空载功率设计如下：

p1=[2,1.6,1.8,2.4,2.4,4.1]
p2=[0.6,0.6,0.5,0.4,0.4,0.6]

具体的解码在fjsp.py文件里。

主函数

设计主函数如下：

oj=data_deal(10,6)               #工件数，机器数
Tmachine,Tmachinetime,tdx,work,tom=oj.cacu()
parm_data=[Tmachine,Tmachinetime,tdx,work,tom]
to=FJSP(10,6,0.5,parm_data)      #工件数，机器数，选择最短机器的概率和mk01的数据
oh=mul_op()
ho=nsga_II(50,100,to,oh,work,10)     #数50,100,10分别代表迭代的次数、种群的规模、工件数
#to是柔性车间模块，oh是多目标模块pareto,pareto_job,pareto_machine,pareto_time,fit_every=ho.nsga_total()  #最后一次迭代的最优解
print(pareto)
oh.draw_change(fit_every)        #每次迭代过程中pareto解中3个目标的变化
oh.draw_3d(pareto)               #结果的3维图sig=0
job,machine=np.array([pareto_job[sig]]),np.array([pareto_machine[sig]])
machine_time=np.array([pareto_time[sig]])
# C_finish,Twork,E_all,list_M,list_S,list_W,tmax=to.caculate(job,machine,machine_time)
to.draw(job,machine,machine_time)#画pareto解的第一个解的甘特图

运行结果

结果如下：

帕累托解中完工时间、负荷、能耗随迭代次数的变化图如下：

帕累托解的3维图如下：

帕累托解的第一个解的甘特图如下：

代码

有5个py文件和一个mk01的text文档：

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

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

公众号二维码：