车间调度系列文章：
-
18、书本算法重现丨单目标遗传算法（GA3):以算例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 | 工件号 | Xijz | 工序oij是否在机器z上加工，为0-1变量，在z上加工为1 |
| j,k | 工序号 | Gijhk | 工序oij和工序ohk的先后顺序，为0-1变量，ij在前为1 |
| z | 机器号 | M | 一个大正实数 |
| qi | 工件i的工序数 | Tijz | 工序oij在机器z的加工时间 |
| oij | 工件i的第j道工序 | ||
| Mij | 工序oij的可选机器 | ||
| Sij | 工序oij的开工时间 | ||
| Cij | 工序oij的完工时间 | ||
| Load_z | 机器负荷 |
目标函数：
min makespan=max(makespan_i) makespan_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得到如下编码：
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<p1 or r>1-p2:
for k in range(len(n_machine)):
m=int(n_machine[k])-1
index_select.append(m)
t=n_time[k]
a_global[0,m]+=t #全局负荷计算
a_part[0,m]+=t #局部负荷计算
if r<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=1
del rest[sig][m] #删掉空闲时间段
break
if 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里。
GA3算法
具体的算法参考文献有介绍，介绍一下步骤：
算法步骤：
-
步骤1：固定比例全局选择、局部选择和随机选择3种方式初始多个工序、机器、加工时间编码，并解码
-
步骤2：判断是否达到最大迭代次数，是的话输出结果，否则转到步骤3
-
步骤3：锦标赛选择相同规模的个体
-
步骤4：在交叉概率Pc对选择个体的工序编码进行pox交叉，机器编码进行均匀交叉
-
步骤5：在变异概率Pm下对工序编码进行领域搜索变异，对机器编码进行选择最短加工机器变异
-
步骤6：转到步骤2
锦标赛选择
每次从种群中选择ｋ个个体进行完工时间比较， 将完工时间低的个体插入到交叉池中，如此循环，直到填满交叉池为止。
python里实现比较容易，代码如下：
def select(popsize,number,C_size,T_answer):
next_index=[]
for i in range(popsize):
cab=random.sample(range(number),C_size) #按照C_size生成一组不重复的索引用于锦标赛选择
index,Z=[],[]
for j in range(C_size):
index.append(cab[j]),Z.append(T_answer[cab[j]]);
min_Z=Z.index(min(Z))
min_idx=index[min_Z]
next_index.append(min_idx)
return next_index
其中第一次迭代从父代选择个体，以后的迭代是父代和子代合并选择个体，并且当选择出来的最优解时不优于合并最优解，保存合并最优解
W,M,T=np.vstack((work_job,work_job1)),np.vstack((work_M,work_M1)),np.vstack((work_T,work_T1))
A=answer+answer1
Z_min=A.index(min(A))
next_index=select(popsize,2*popsize,C_size,A)
work_job,work_M,work_T=W[next_index],M[next_index],T[next_index]
answer=np.array(A)[next_index].tolist()
best_index=answer.index(min(answer))
if min(answer)>min(A): #选择出来的最优解时不优于合并最优解，保存合并最优解
work_job[best_index],work_M[best_index],work_T[best_index]=W[Z_min],M[Z_min],T[Z_min]
answer[best_index]=min(A)
机器的均匀交叉
均匀交叉算子的概念比较简单，简单说一下逻辑：假设两个解的工序编码的第一道工序分别选择了机器1和3，随机生成0，1两个数中的一个，如果随机数是1，交换两个解第一道工序的机器选择，否则保持原选择。
代码：
def ma_cross(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
工序的pox交叉
以mk01为例：随机0到9的一个数为6，对应两个进行交叉的工序编码，0到6基因及其位置保持不变，每个编码6到9的基因位置按顺序填入另一个工序编码6到9的基因。
代码：
def job_cross(chrom_L1,chrom_L2): #工序的pox交叉
num=list(set(chrom_L1))
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[0]))-1,np.zeros((1,chrom_L1.shape[0]))-1
sig,svg=[],[]
for i in range(chrom_L1.shape[0]):#固定位置的工序不变
ii,iii=0,0
for j in range(len(jpb_set1)):
if(chrom_L1[i]==jpb_set1[j]):
C1[0,i]=chrom_L1[i]
else:
ii+=1
if(chrom_L2[i]==jpb_set1[j]):
C2[0,i]=chrom_L2[i]
else:
iii+=1
if(ii==len(jpb_set1)):
sig.append(chrom_L1[i])
if(iii==len(jpb_set1)):
svg.append(chrom_L2[i])
signal1,signal2=0,0 #为-1的地方按顺序添加工序编码
for i in range(chrom_L1.shape[0]):
if(C1[0,i]==-1):
C1[0,i]=svg[signal1]
signal1+=1
if(C2[0,i]==-1):
C2[0,i]=sig[signal2]
signal2+=1
return C1[0],C2[0]
最短加工机器变异
比较简单，就是工序选择最短加工时间的机器，当然，如果原工序已经选择了最大短加工时间的机器，该变异不会对编码产生改变。
代码：
def ma_mul(w,m,t,parm_fj,parm_mk):
job_num=parm_fj[0]
Tmachine,Tmachinetime,machines,work,start,lenth=parm_mk[0],parm_mk[1],parm_mk[2],parm_mk[3],parm_mk[4],parm_mk[5]
r=np.random.randint(1,w.shape[0]+1,1)[0] #变异的数量
mul=random.sample(range(w.shape[0]),r) #变异的位置
count=np.zeros((1,job_num),dtype=np.int)
signal=0
sig=0
for i in range(len(machines)):
for j in range(machines[i]):
if signal<len(mul) and sig == mul[signal] :
highs=start[sig][count[0,sig]]
lows=start[sig][count[0,sig]]-lenth[sig][count[0,sig]]
n_machine=Tmachine[sig,lows:highs].tolist()
n_time=Tmachinetime[sig,lows:highs].tolist()
time=min(n_time) #最短加工时间的
index=n_time.index(time)
mac=n_machine[index] #最短加工时间的机器
m[sig]=mac #更新编码
t[sig]=time
signal+=1
sig+=1
return w,m,t
领域搜索变异
参考文献的算法步骤：
步骤1： 在变异染色体中随机选择ｒ个不同基因，并生成其排序的所有邻域。
步骤２ 评价所有邻域的适应值，选出最佳个体作为子代。
可以看出，邻域搜索就是对r个基因进行排队，邻域的个数是r!-1,本文取2到5个工件进行变异，领域是用一种递归交换位置的方法
代码：
def perm(begin,end,A,nums): #领域寻找
if begin == end:
A.append(nums.copy())
else:
for i in range(begin,end+1):
nums[begin],nums[i]=nums[i],nums[begin] #交换对应位置
perm(begin+1,end,A,nums)
nums[begin],nums[i]=nums[i],nums[begin]
return A
def location(job_num,job):
r=np.random.randint(2,6,1)[0] #生成2到5个变异的工件
w=random.sample(range(10),r) #生成对应工件编号
A=[]
for i in w:
loc=job.index(i) #工件出现的第一个位置
A.append(loc)
index=[]
index=perm(0,r-1,index,A) #领域寻找函数
return index
本文设计整个邻域搜索变异的逻辑：新解优于原解就更新编码
代码：
def job_mul(w,m,t,parm_fj,parm_mk):
C_finish,_,_,_,_=fjsp.caculate(w,m,t,parm_fj,parm_mk)
t1=C_finish #原编码的目标计算
job_num=parm_fj[0]
job=w.tolist()
result=location(job_num,job)
mul=np.array(result[0])
for re in result[1:] : #位置索引
w[mul]=w[re] #更新编码
C_finish,_,_,_,_=fjsp.caculate(w,m,t,parm_fj,parm_mk)
#新编码的目标计算
if C_finish>t1 : #如果新编码目标都劣于原目标，不更新编码
w[re]=w[mul]
else: #否则，更新编码，并更新目标
t1=C_finish
return w,m,t
整个算法迭代的核心代码：
for i in range(0,popsize,2): #种群规模下每次选2个个体
W1,M1,T1=work_job[i],work_M[i],work_T[i]
W2,M2,T2=work_job[i+1],work_M[i+1],work_T[i+1]
if np.random.rand()<Pc:
M1,T1,M2,T2=ma_cross(M1,T1,M2,T2) #机器交叉
W1,W2=job_cross(W1,W2) #工序交叉
if np.random.rand()<Pm :
W1,M1,T1=ma_mul(W1,M1,T1,parm_fj,parm_mk) #机器变异
W1,M1,T1=job_mul(W1,M1,T1,parm_fj,parm_mk) #工序变异
W2,M2,T2=ma_mul(W2,M2,T2,parm_fj,parm_mk) #机器变异
W2,M2,T2=job_mul(W2,M2,T2,parm_fj,parm_mk) #工序变异
以上所有代码放在GA3.py里。
结果
代码运行环境
windows系统，python3.6.5,第三方库及版本号如下：
pandas==0.25.3
numpy==1.18.5
matplotlib==3.2.1
xlrd==1.2.0
第三方库需要在安装完python之后，额外安装，以前文章（点击我跳转）有讲述过安装第三方库的解决办法。
命令
import fjsp
import data
import GA3
Tmachine,Tmachinetime,machines,work,start,lenth=data.total(10)
#mk01的相关数据，10是工件数
parm_fj=[10,6,0.3,0.4] #编码解码参数，依次是工件数、机器数 ，全局、局部选择比例
parm_mk=[Tmachine,Tmachinetime,machines,work,start,lenth]
parm_ga=[20,50,2,0.8,0.2] #遗传参数，依次是迭代次数、种群规模、锦标赛个体k，交叉和变异概率
job,machine,machine_time=GA3.ga(parm_ga,parm_mk,parm_fj)
#返回的最优个体
# C_finish,_,_,_,_=fjsp.caculate(job,machine,machine_time,parm_fj,parm_mk)
# print(C_finish)
fjsp.draw(job,machine,machine_time,parm_fj,parm_mk)
运行结果
20次次迭代结果如下：
解的甘特图如下：
小结
算法基本是复现原论文，也有自己的设计，文章的算法也有比较详细的介绍，插入式解码和邻域搜索等一些方法是比较不错的，希望对大家的论文有一些参考价值。
excel数据可更改，工件数、机器数、工件的工序数、工序的可加工机器数等数据对得上就能运行。
参考文献：柔性作业车间调度智能算法及其应用-高亮（第三章)
代码
有5个py文件和一个mk01的excel文件，以及《柔性作业车间调度智能算法及其应用-高亮》电子书。因篇幅问题，代码附在后面。
演示视频：
算法重现第三章丨GA3单目标
完整算法源码+数据：见下方微信公众号：关注后回复：车间调度
# 微信公众号：学长带你飞
# 主要更新方向：1、车辆路径和柔性车间调度问题求解算法
# 2、学术写作技巧
# 3、读书感悟
# @Author : Jack Hao
公众号二维码：
1336
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
