柔性车间调度丨模拟退火算法研究:以算例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为例

• 17、柔性车间调度丨复现改进人工蜂群算法论文

• 18、书本算法重现丨单目标遗传算法（GA3):以算例MK01为例

• 19、多目标柔性车间调度丨改进灰狼算法:以算例MK01为例

• 20、多目标柔性车间调度丨NSGA3:以算例MK01为例

• 21、柔性车间调度丨模拟退火算法研究:以算例MK01为例

1、柔性车间调度问题

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

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

2、数学模型

符号定义：

n	工件总数	makespani	工件i的完工时间
m	机器总数	makespan	最大完工时间
i,h	工件号	Load_all	机器总负荷
j,k	工序号	E_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的开工时间	Fz	机器的z的完工时间
Cij	工序oij的完工时间	Pz	机器的z的负载功率
Load_z	机器z负荷	Rz	机器的z的空载功率

模型：

目标函数：

​ makespan=min{max makespani} 完工时间

i=(1,2,…,n)

约束条件：

3、柔性作业车间工具

工序编码

• 步骤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])-1
                        index_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=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]       #更新工序完工时间

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

代码在fjsp.py里。

4、 模拟退火算法

算法介绍：

模拟退火算法(Simulated Annealing，SA)最早的思想是由N. Metropolis 等人于1953年提出。1983 年,S. Kirkpatrick 等成功地将退火思想引入到组合优化领域。它是基于Monte-Carlo迭代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。 从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在求解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。

算法步骤

1、初始温度T0 ， 令 T = T0，随机产生一个初始解 X0 ，并计算对应的目标函数值 E(x0) ，本文是makespan。
2、令T=rT， 其中r取值0到1之间，为温度下降速。
3、对当前解施加随机扰动，在其邻域内产生一个新解 ，并计算对应的完工时间，
4、按照 Metropolis准则判断是否接受新解 ， 若新解的完工时间小于原解，接受新解作为当前解，否则按照概率判断是否接受。
5、在温度 T 下，重复 L 次扰动和接受过程，即执行步骤 3 和 4 ，本问的L就是种群数量。
6、判断温度是否达到终止温度水平，若是则终止算法，否则返回步骤 2
参考博客：https://blog.csdn.net/qq_44865735/article/details/124100778：

核心代码：

                if i != best_idx:          # 保留种群最优个体
                    job1, machine1, machinetime1 = w_new[i], m_new[i], t_new[i]
                    C_max, _, _, _, _ = fjsp.caculate(job1, machine1, machinetime1, self.parm_fj, self.parm_mk)
                    if C_max< answer[i]:   # 如果新的完工时间低于原解
                        answer[i] =C_max   # 接受新解
                        w[i] = job1
                        m[i] = machine1
                        t[i] = machinetime1
                    else:                 # 否则以一定概率接受原解
                        delta_e = C_max- answer[i]
                        if random.uniform(0, 1) < math.exp(-delta_e / self.T):
                            answer[i] = C_max
                            w[i] = job1
                            m[i] = machine1
                            t[i] = machinetime1

            self.T *= self.r      # 降温

本文为保留每次迭代的最优个体，最优个体不与新解比较

工序邻域搜索（工序新解产生）

随机产生2个位置，2个位置之间的编码逆序。

代码：

    def job_new(self,w):
        w_new = []
        index = random.sample(range(len(w[0])), 2)  # 随机生成两个位置
        idx = list(set(index))                      # 位置按先后顺序调整
        for wi in w:
            job1 = wi.copy()
            b = wi[idx[0]:idx[1] + 1] 
            b.reverse()                             # 逆序
            job1[idx[0]:idx[1] + 1] = b
            w_new.append(job1)
        return w_new

加工机器邻域搜索（机器新解产生）

随机产生2个位置，交换相邻2个机器编码2个位置之间的基因

代码：

 def mt_new(self, m, t):
        index = random.sample(range(len(m[0])), 2)
        idx = list(set(index))
        m_new, t_new = [], []
        for i in range(0, len(m), 2):
            m1, m2, t1, t2 = m[i].copy(), m[i + 1].copy(), t[i].copy(), t[i + 1].copy()
            a, b, c, d = m1[idx[0]:idx[1] + 1], m2[idx[0]:idx[1] + 1], t1[idx[0]:idx[1] + 1], t2[idx[0]:idx[1] + 1]
            m1[idx[0]:idx[1] + 1], m2[idx[0]:idx[1] + 1] = b, a  # 两段加工机器编码交换
            m_new.append(m1)
            m_new.append(m2)

            t1[idx[0]:idx[1] + 1], t2[idx[0]:idx[1] + 1] = d, c  # 两段加工时间编码交换
            t_new.append(t1)
            t_new.append(t2)
        return m_new, t_new

5、结果

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

numpy==1.24.3
matplotlib==3.7.1

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

命令

在main.py里

import data
from SA import simul
import fjsp

Tmachine,Tmachinetime,machines,work,start,lenth=data.total(10)
parm_fj=[10,6,0.3,0.4]                  # 编码解码参数，依次是工件数、机器数 ，全局、局部选择比例
parm_mk=[Tmachine,Tmachinetime,machines,work,start,lenth]
parm_sa = [100, 0.95, .01, 40]          # 初始温度，降温速率，终止温度，种群规模

sa = simul(parm_sa, parm_fj, parm_mk)   # 模拟退火模块
job, machine, machinetime, result = sa.simulate()      # 模拟退火算法结果

fjsp.draw(job, machine, machinetime,parm_fj,parm_mk)   # 甘特图
fjsp.change(result)                                    # 完工时间变化图

一次运行结果

部分结果：

第162次迭代，完工时间：46.0
第163次迭代，完工时间：46.0
第164次迭代，完工时间：46.0
第165次迭代，完工时间：46.0
第166次迭代，完工时间：46.0

甘特图：

最大完工时间随迭代次数的变化：

6、代码

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

演示视频：

柔性车间调度丨模拟退火算法研究：以MK01算例为例

完整算法+数据：

# 微信公众号学长带你飞
# 主要更新方向：1、车间调度问题求解算法
#              2、学术写作技巧
#              3、读书感悟
# @Author  : Jack Hao
可加微信：diligent-1618，请备注：学校-专业-昵称。