数字化车间智能排产调度挑战赛（一）—

科大讯飞–数字化车间智能排产调度挑战赛

本系列文章用于记录比赛中模型构建，算法设计，仅用于记录与学习。
系列文章将分为一下几个部分

分析问题，建立数学模型构建，并基于求解器验证
设计启发式规则求解车间调度问题
关键路径+VNS的混合算法求解车间调度问题

这三个部分也是我在解决这个问题过程中，求解方法的一个进化过程。第一次接触车间调度这类问题，涉及的内容也不会很深。下面就开始我们的第一部分内容：分析问题，构建数学模型，小规模样例验证。

1. 问题背景

具体的题目点击这里
这里就简单概述一下。
大赛提供了产品生产数量及产品的工艺路线，我们需求合理安排产品加工使用的机器编号已加工开始时间和结束时间，并需要满足如下约束：

每个设备只能加工某一件产品的一种工序，且一旦开始不能中断。
产品加工顺序必须严格按照工艺路线执行，一道工序结束后才可以开始下一道工序。
每一道工序可在给定设备集合中的任一设备上加工，每道工序只能同时在一个设备上加工。
工序B加工结束后必须立即开始工序C，同时工序B加工时间必须在给定单位加工时间范围内。
1. 工序B对应的机器需求设备准备时间
  1. 产品中第一次出现B
  2. 当前产品存在多次工序B，且其使用机器与上一次使用机器不同；
  3. 当前机器加工的产品与当前机器上一次加工的产品不同；
2. 无需增加准备时间
  1. 当前产品存在多次工序B，且其使用机器与上一次使用机器相同
每一产品中工序C以节为单位加工，且每一产品节与节之间必须连续加工。
这里需要额外说明的是，工序AB是以批为单位生产，工序CD是以节为单位生产，每一个产品编号代表一批，每一批中包含若干节数量。（具体看题目哦）

评价指标
e = 工序 C 在对应设备上的累计加工时间所有产品加工完工时间 e = \frac{工序C在对应设备上的累计加工时间}{所有产品加工完工时间} e=所有产品加工完工时间工序C在对应设备上的累计加工时间

2 分析与建模

2.1 数据分析

大赛提供的样本数据中共有40个产品，36台机器，4中工艺路线，每种工序有对应不同机器类型，且数量也不同。

工序	机器类型	机器数量	总加工时长	平均加工时长
工序A	T-01	2	9600	4800
工序B1	T-03	12	57700	4808
工序B2	T-04	6	19100	3183
工序C	T-02	1	11440	11440
工序D	T-05	15	165120	11008

这里的加工单位为min（分钟）。
我的思路是把所有产品需要加工的工序数看作任务，然后这些任务之间存在一下约束关系，比如紧前关系、准备时间等。所以数据处理前后的表是下面这样的。

index	route_id	route_No	name	equ_type	product_id	work_time
0	A1	1	工序A	T-01	P202201101	240
1	A1	2	工序B	T-03	P202201101	480
2	A1	3	工序C	T-02	P202201101	96
3	A1	4	工序D	T-05	P202201101	480
4	A1	4	工序D	T-05	P202201101	480
5	A1	4	工序D	T-05	P202201101	480
6	A1	4	工序D	T-05	P202201101	480
7	A1	4	工序D	T-05	P202201101	480
8	A1	4	工序D	T-05	P202201101	480
9	A1	4	工序D	T-05	P202201101	480
10	A1	4	工序D	T-05	P202201101	480
11	A1	5	工序B	T-03	P202201101	300
—	—	—	—	—	—	—
24	A1	13	工序D	T-05	P202201101	480
25	A1	13	工序D	T-05	P202201101	480
26	A1	13	工序D	T-05	P202201101	480
27	B1	1	工序A	T-01	P202201102	240
28	B1	2	工序B	T-03	P202201102	200
29	B1	3	工序C	T-02	P202201102	48

你可能会有疑问，为啥一个产品有那么多个连续的工序D但只会出现一个工序C呢？
哎，因为工序CD的生产单位是按“节”生产的，所以相当于直接分开生产；又因为工序C的加工机器只有一台，而且要求连续生产，所以干脆合并到一起就好啦。

2.2 数学建模

有了前面的数据分析处理，这里的变量与决策变量比较容易确定，一个是加工顺序，一个是开始时间。
x i j k x_{ijk} xijk：表示加工顺序，机器 k k k加工为任务 i i i后加工任务 j j j。（注意我这里把所有产品工序，处理成单个任务了）
t i t_i ti：任务 i i i的开始加工时间
下面来单独看每条约束
先看前三条
① 每个设备只能加工某一件产品的一种工序，且一旦开始不能中断。
② 产品加工顺序必须严格按照工艺路线执行，一道工序结束后才可以开始下一道工序。
③ 每一道工序可在给定设备集合中的任一设备上加工，每道工序只能同时在一个设备上加工。
比较好理解，就是每个任务仅且只能被加工一次、任务间存在紧前关系、同一个机器的加工顺序。就可以构建如下约束：
∑ j ∈ N , k ∈ M x i j k = 1 t i + w i ≤ t j t j ≥ t i + w i − B ( 1 − x i j k ) \sum_{j\in N,k\in M}x_{ijk} =1 \\ t_i + w_i \le t_j \\ t_j≥t_i+w_i-B(1-x_{ijk} ) j∈N,k∈M∑xijk=1ti+wi≤tjtj≥ti+wi−B(1−xijk)
④ 这个约束可以分成两个部分，一个是工序B结束后工序C立即开始；另一个是工序B准备时间问题。
4.1 工序B结束后工序C立即开始。可以用一个集合 A i j A_{ij} Aij表示工序BC的集合，对于这部分的任务呢存在如下约束
t j = t i + w i , ∀ i , j ∈ A i j t_j=t_i+w_i,∀i,j∈A_{ij} tj=ti+wi,∀i,j∈Aij
4.2 工序B的准备时间。关于需要准备时间的描述给出了三种，不需要准备时间的只有一种，那我们的想法是只关注不需要准备时间的那一种情况，剩余的都需要准备时间。
逻辑也很简单，只有同一个产品同一个加工机器类型的工序B才有可能不需要准备时间。这里我们引入一个集合 B B BB BB，一个变量 y i y_i yi。
B B BB BB：表示同一产品同一加工类型的工序B前同类的任务序号，例如{11:[1], 15:[1,11], 17:[1,11,15]},表示任务11前面存在任务1跟它是同一产品同一类型。
y i y_i yi：表示任务 i i i是否需要准备时间。如果需要则为1，否则为0。
就可以构建如下约束了：
y j + ∑ k ∈ M x i j k = 1 , ∀ i , j ∈ B B y_j+∑_{k∈M}x_{ijk} =1,∀i,j∈BB yj+k∈M∑xijk=1,∀i,j∈BB
同时在引入一个变量 d i d_i di记录实际每个任务的加工时间。
d i = w i + 30 ∗ y i , ∀ i ∈ N d_i=w_i+30*y_i ,∀i∈N di=wi+30∗yi,∀i∈N
到这一步，需要将前面出现过 w i w_i wi的约束都替换成 d i d_i di。
约束⑤其实我们已经在前面做了处理，将每个产品每个步骤的工序C合并。
到这里，模型大体已经构建完成，下面来看看完整的模型参数、变量、目标及约束。

2.3 数学模型

参数

符号	定义
k ∈ M k∈M k∈M	设备序号
h ∈ H h∈H h∈H	产品集合
w i w_i wi	任务i的加工时长
P i P_i Pi	任务i的紧前任务集合
C k C_k Ck	设备k可以加工任务i
A i j A_{ij} Aij	工序C开始时间等于工序B的结束时间
B B BB BB	表示同一产品同一加工类型的工序B前同类的任务序号，例如{11:[1], 15:[1,11], 17:[1,11,15]},表示任务11前面存在任务1跟它是同一产品同一类型。
B i g M BigM BigM	极大值

变量

符号	定义
t i t_i ti	任务i的开始加工时间
x i j k x_{ijk} xijk	设备k加工完i之后加工j，01变量
y i y_i yi	工序是否需要准备时间
d i d_i di	加上准备时间后的加工时长

目标函数：总加工时间最小
m i n m a x ( t i ) min max(t_i) minmax(ti)
约束条件

每个任务都被服务一次
1. 机器从虚拟起点出发
  ∑ j ∈ C k x 0 j k = 1 , ∀ k ∈ M ∑_{j∈C_k}x_{0jk}=1,∀k∈M j∈Ck∑x0jk=1,∀k∈M
2. 每个任务被服务一次
  ∑ i ∈ D j ∪ 0 ∑ k ∈ O j x i j k = 1 , ∀ j ∈ N / 0 ∑_{i∈D_j∪{0}}∑_{k∈O_j}x_{ijk} =1,∀j∈N /{0} i∈Dj∪0∑k∈Oj∑xijk=1,∀j∈N/0
3. 任务的流平衡
  ∑ i ∈ D j ∪ 0 x i j k − ∑ i ∈ D j ∪ 0 x j i k = 0 , ∀ k ∈ M , j ∈ C k ∑_{i∈D_j∪{0}}x_{ijk} -∑_{i∈D_j∪{0}}x_{jik} =0 ,∀k∈M,j∈C_k i∈Dj∪0∑xijk−i∈Dj∪0∑xjik=0,∀k∈M,j∈Ck
任务的开工时间满足时序关系
t j ≥ t i + d i , ∀ j ∈ N , i ∈ P j t_j≥t_i+d_i,∀j∈N,i∈P_j tj≥ti+di,∀j∈N,i∈Pj
每个产品的工序C紧跟着工序B结束立马开工
t j = t i + d i , ∀ i , j ∈ A i j t_j=t_i+d_i,∀i,j∈A_{ij} tj=ti+di,∀i,j∈Aij
同一个设备同时只能加工一道工序
t j ≥ t i + d i − B ( 1 − x i j k ) , ∀ i , j ∈ C k , k ∈ M t_j≥t_i+d_i-B(1-x_{ijk} ),∀i,j∈C_k,k∈M tj≥ti+di−B(1−xijk),∀i,j∈Ck,k∈M
当前产品存在多次工序B，且使用机器与上次相同（工序B是否需要准备时间）
y j + ∑ k ∈ M x i j k = 1 , ∀ i , j ∈ B B y_j+∑_{k∈M}x_{ijk} =1,∀i,j∈BB yj+k∈M∑xijk=1,∀i,j∈BB
加上准备时间后的加工时长
d i = w i + y i ∗ T , ∀ i ∈ N d_i=w_i+y_i*T ,∀i∈N di=wi+yi∗T,∀i∈N
变量取值范围
x i j k , y i ∈ { 0 , 1 } t i , d i ≥ 0 x_{ijk},y_i \in \{0,1\} \\ t_i, d_i \ge 0 xijk,yi∈{0,1}ti,di≥0

3. 小规模样例验证

代码分两个部分，一部分是数据处理，另一部分是模型求解，这里用了Gurobi求解器进行求解。

"""
用于处理模型数据
1. 将产品信息划分为多个子任务，同时换算加工时间
2. 构建任务-设备资质矩阵
"""
import itertoolsimport pandas as pd
import numpy as np
import mathdef result_process(product_info, route_info):# 将产品信息划分为多个子任务result = []for i in range(len(product_info)):route_id = product_info.iloc[i, 2]route_temp = route_info[route_info['route_id'] == route_id].reset_index(drop=True)route_temp['product_id'] = product_info.iloc[i, 0]route_temp['work_time'] = 0for j in range(len(route_temp)):if route_temp.iloc[j, 4][-1] == 'h':if route_temp.iloc[j, 2] == "工序C":route_temp.loc[j, 'work_time'] = float(route_temp.iloc[j, 4][:-1]) * product_info.iloc[i, 1] * 60else:route_temp.loc[j, 'work_time'] = float(route_temp.iloc[j, 4][:-1]) * 60else:route_temp.loc[j, 'work_time'] = float(route_temp.iloc[j, 4][:3])route_temp.loc[j, 'ready_time'] = float(route_temp.iloc[j, 4][4:7])route_temp = route_temp.fillna(0)route_D = route_temp[route_temp['name'] == '工序D']for j in range(product_info.iloc[i, 1] - 1):route_temp = pd.concat([route_temp, route_D], axis=0)if i == 0:result = route_tempelse:result = pd.concat([result, route_temp], axis=0)result = result.sort_values(['product_id', 'route_No']).reset_index(drop=True)result = result.sort_values(['product_id', 'route_No']).reset_index(drop=True).reset_index()return resultclass Data:def __init__(self):self.order = []self.equ_info = []self.order_num = 0self.pro_route = []self.BC_order = []self.equ_order = []self.order_equ = []self.conf_ = []self.BB_order = []self.equ_num = 0self.order_order = []self.ready_B = []def pro_order(self,result):# 获取每个任务的紧前任务pro_route = {1: [0]}BC = []for i in range(2, len(result)):if result.loc[i, 'product_id'] == result.loc[i - 1, 'product_id']:if result.loc[i, 'route_No'] == result.loc[i - 1, 'route_No'] + 1:temp_pro = []count = 1for j in range(i, 0, -1):if result.loc[i, 'route_No'] == result.loc[i - count, 'route_No'] + 1:temp_pro.append(i - count)count += 1else:breakpro_route[i] = temp_proelse:pro_route[i] = pro_route[i - 1]# 获取B、C的组合if result.loc[i, 'name'] == '工序C' and result.loc[i - 1, 'name'] == '工序B':BC.append([i - 1, i])self.pro_route = pro_routeself.BC_order = BCdef conf_equ_order(self, result, equ_info):# 构建任务与设备的资质矩阵、每个设备的子任务集合、每个子任务的设备集合conf_ = np.zeros((len(result), len(equ_info)))equ_order = {}for j in range(len(equ_info)):order = []for i in range(len(result)):if result.loc[i, 'equ_type'] == equ_info.loc[j, 'equ_type']:conf_[i][j] = 1order.append(i)equ_order[j] = ordervar_start = np.zeros(len(equ_info))var_start[0] = 1var_start[1] = 1conf_ = np.insert(conf_, len(conf_), np.ones(len(equ_info)), axis=0)self.conf_ = conf_self.equ_order = equ_orderorder_equ = {}order_order = {}for i in range(len(result)):order_equ[i] = np.where(conf_[i][:] == 1)[0]order_order[i] = equ_order[order_equ[i][0]]self.order_equ = order_equself.order_order = order_orderdef BB_index(self, result, product_info):BB_order = {}for i in range(len(product_info)):p_id = product_info.loc[i, 'product_id']df = result[(result['product_id'] == p_id) & (result['name'] == '工序B')]df = df.reset_index(drop=True)self.ready_B.append(df.loc[0, 'index'])for j in range(1, len(df)):ind = df.loc[j, 'index']equ_type = df.loc[j, 'equ_type']df1 = df[df['equ_type'] == equ_type]arr = np.array(df1['index'])arr = arr[arr < ind]if len(arr) > 0:BB_order[ind] = arrelse:self.ready_B.append(df.loc[j, 'index'])self.BB_order = BB_orderdef readData(self, product_info):route_info = pd.read_csv('../data/工艺路线.csv', encoding="gbk")equ_info = pd.read_csv('../data/设备信息.csv', encoding="gbk")self.equ_info = equ_info# 0表示D分批，1表示D不分批result = result_process(product_info, route_info)# 获取每个任务的紧前任务self.pro_order(result)# # 构建任务与设备的资质矩阵、每个设备的子任务集合self.conf_equ_order(result, equ_info)# 同一个产品的BB工序序号self.BB_index(result, product_info)self.order_num = len(result)self.order = resultself.equ_num = len(equ_info)

模型求解

import itertoolsimport pandas as pd
import numpy as np
import gurobipy as grb
from gurobipy import GRBfrom Data import Data
'''图形工厂'''
import plotly.figure_factory as ffif __name__ == '__main__':product = pd.read_csv('../data/产品信息.csv')# 取前4条数据product = product.loc[:3]data = Data()data.readData(product)task_num = data.order_num + 1equ_num = data.equ_numCk = data.equ_orderconf = data.conf_pro_ = data.pro_routeBC = data.BC_orderw = np.array(data.order['work_time'])T = 30B = 100000equ_info_ = data.equ_infoOE = data.order_equOO = data.order_orderBB = data.BB_orderRB = data.ready_B# 构建模型print("========模型构建========")model = grb.Model()# 创建变量 包含一个虚拟开始任务,放到最后t_i = model.addVars(task_num, lb=0, ub=np.inf, vtype=GRB.INTEGER, name='t')y_i = model.addVars(task_num, vtype=GRB.BINARY, name='y')d_i = model.addVars(task_num, lb=0, ub=np.inf, vtype=GRB.INTEGER, name='d')x_ijk = model.addVars(task_num, task_num, equ_num, vtype=GRB.BINARY, name='x')f = model.addVar(lb=0, ub=np.inf, vtype=GRB.INTEGER, name='f')for i in range(task_num):for j in range(task_num):for k in range(equ_num):if i == j:x_ijk[i, j, k].setAttr('ub', 0)for k in range(equ_num):x_ijk[task_num - 1, task_num - 1, k].setAttr('ub', 1)for i in RB:y_i[i].setAttr('lb', 1)# 目标函数model.setObjective(f, GRB.MINIMIZE)for i in range(task_num):model.addConstr(f >= t_i[i] + d_i[i], name=f'max_time_{i}')# 约束条件# 1.1. 设备从虚拟起点出发for k in range(equ_num):lin = grb.LinExpr()lin.addTerms(1, x_ijk[task_num - 1, task_num - 1, k])for j in Ck[k]:lin.addTerms(1, x_ijk[task_num - 1, j, k])model.addConstr(lin == 1, name=f'equ_start_{k}')# 1.2. 每个任务被服务一次。流只能从与该任务相关的任务流入，只能被该任务相关的设备服务for j in range(task_num - 1):lin = grb.LinExpr()for k in OE[j]:lin.addTerms(1.0, x_ijk[task_num - 1, j, k])for i in OO[j]:lin.addTerms(1, x_ijk[i, j, k])model.addConstr(lin == 1, name=f'fulfill_{j}')# 1.3. 每个任务的流平衡for j in range(task_num - 1):for k in OE[j]:lin = grb.LinExpr()lin.addTerms(1, x_ijk[task_num - 1, j, k])lin.addTerms(-1, x_ijk[j, task_num - 1, k])for i in OO[j]:lin.addTerms(1, x_ijk[i, j, k])lin.addTerms(-1, x_ijk[j, i, k])model.addConstr(lin == 0, name=f'flowBal_{j}_{k}')# 2. 紧前关系for j in pro_.keys():for i in pro_[j]:model.addConstr(t_i[j] >= t_i[i] + d_i[i], name=f'order_{i}_pro_{j}')# 5. 加上准备时间后的加工时长for i in range(task_num - 1):model.addConstr(d_i[i] == w[i] + y_i[i] * T, name=f'act_wt_{i}')# 6. 同一个设备同时只能加工一道工序for key in Ck.keys():order_arr = itertools.permutations(Ck[key], 2)for arr in order_arr:i = arr[0]j = arr[1]k = keymodel.addConstr(t_i[j] >= t_i[i] + d_i[i] - B * (1 - x_ijk[i, j, k]), name=f'equ_{k}_{i}_before_{j}')# # 4. 工序B是否需要准备时间for j in BB.keys():for i in BB[j]:lin = grb.LinExpr()for k in OE[j]:lin.addTerms(1, x_ijk[i, j, k])lin.addTerms(1, y_i[j])model.addConstr(lin == 1, name=f'BB_{i}_{j}')# 3. BC工序for arr in BC:i = arr[0]j = arr[1]model.addConstr(t_i[j] == t_i[i] + d_i[i], name=f'BC_{i}_{j}')# 优化求解model.write('M.lp')model.Params.MIPGap = 0.1model.Params.NoRelHeurTime = 30model.Params.TimeLimit = 180# model.Params.Heuristics = 0.7# Optimize modelmodel.optimize()if model.status == GRB.INFEASIBLE:model.computeIIS()model.write('M.ilp')else:# 输出求解结果# 格式：产品id，工序id，设备名，开始时间，持续时长，结束时间re = []for key in x_ijk.keys():if x_ijk[key].X > 0:if key[1] != task_num - 1:re.append([key[1], equ_info_.loc[key[2], 'equ_name'], t_i[key[1]].X, w[key[1]]])# 在子任务后面加上这些信息re = pd.DataFrame(re)re.columns = ['index', 'equ_name', 'start', 'act_wt']result_ = pd.merge(data.order, re, left_on='index', right_on='index', how='left')result_['end'] = result_['start'] + result_['act_wt']# result_.to_csv('result_init.csv', encoding='utf-8', index=False)# 绘制甘特图gantt = result_[['equ_name', 'product_id', 'route_No', 'start', 'act_wt']]gantt.columns = ['Task', 'product_id', 'route_No', 'Start', 'Finish']gantt['Finish'] = gantt['Start'] + gantt['Finish']gantt['Start'] = gantt['Start'].apply(lambda x: x * 60 * 1000)gantt['Finish'] = gantt['Finish'].apply(lambda x: x * 60 * 1000)gantt = gantt.sort_values(by=['Task'], ascending=0)color = {'P202201101': '#cfff23','P202201102': '#23dcff','P202201103': '#2363ff','P202201104': '#bcb866',}fig = ff.create_gantt(gantt,colors=color,index_col='product_id',show_colorbar=True,  # 显示颜色柱group_tasks=True)fig.show()

结果展示

到此为止，第一部分的问题分析、建模、验证已经结束，由于问题的规模较大，仅用求解器无法在有效时间内找到可行解，更别说好一点的解，所以后续我们还需设计求解算法进行问题求解。
（留坑）