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

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

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

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

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

1. 问题背景

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

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

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

2.3 数学模型

参数

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

变量

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

目标函数：总加工时间最小
$$minmax(t_i)$$
约束条件

1. 每个任务都被服务一次
   1. 机器从虚拟起点出发
      $$\sum _{j\in C_k}{x}_{0jk}=1,\forall k\in M$$
   2. 每个任务被服务一次
      $$\sum _{i\in D_j\cup 0}\sum _{k\in O_j}{x}_{ijk}=1,\forall j\in N/0$$
   3. 任务的流平衡
      $$\sum _{i\in D_j\cup 0}{x}_{ijk}-\sum _{i\in D_j\cup 0}{x}_{jik}=0,\forall k\in M,j\in C_k$$
2. 任务的开工时间满足时序关系
   $$t_j\ge t_i+d_i,\forall j\in N,i\in P_j$$
3. 每个产品的工序C紧跟着工序B结束立马开工
   $$t_j=t_i+d_i,\forall i,j\in A_{ij}$$
4. 同一个设备同时只能加工一道工序
   $$t_j\ge t_i+d_i-B(1-x_{ijk}),\forall i,j\in C_k,k\in M$$
5. 当前产品存在多次工序B，且使用机器与上次相同（工序B是否需要准备时间）
   $$y_j+\sum _{k\in M}{x}_{ijk}=1,\forall i,j\in BB$$
6. 加上准备时间后的加工时长
   $$d_i=w_i+y_i\ast T,\forall i\in N$$
7. 变量取值范围
   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 itertools

import pandas as pd
import numpy as np
import math


def 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'] = 0
        for 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] * 60
                else:
                    route_temp.loc[j, 'work_time'] = float(route_temp.iloc[j, 4][:-1]) * 60
            else:
                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_temp
        else:
            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 result


class Data:
    def __init__(self):
        self.order = []
        self.equ_info = []
        self.order_num = 0
        self.pro_route = []
        self.BC_order = []
        self.equ_order = []
        self.order_equ = []
        self.conf_ = []
        self.BB_order = []
        self.equ_num = 0
        self.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 = 1
                    for 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 += 1
                        else:
                            break
                    pro_route[i] = temp_pro
                else:
                    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_route
        self.BC_order = BC

    def 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] = 1
                    order.append(i)
            equ_order[j] = order
        var_start = np.zeros(len(equ_info))
        var_start[0] = 1
        var_start[1] = 1
        conf_ = np.insert(conf_, len(conf_), np.ones(len(equ_info)), axis=0)
        self.conf_ = conf_
        self.equ_order = equ_order
        order_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_equ
        self.order_order = order_order

    def 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] = arr
                else:
                    self.ready_B.append(df.loc[j, 'index'])
        self.BB_order = BB_order

    def 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 = result
        self.equ_num = len(equ_info)

---

模型求解

import itertools

import pandas as pd
import numpy as np
import gurobipy as grb
from gurobipy import GRB

from Data import Data
'''图形工厂'''
import plotly.figure_factory as ff


if __name__ == '__main__':
    product = pd.read_csv('../data/产品信息.csv')
    # 取前4条数据
    product = product.loc[:3]
    data = Data()
    data.readData(product)

    task_num = data.order_num + 1
    equ_num = data.equ_num
    Ck = data.equ_order
    conf = data.conf_
    pro_ = data.pro_route
    BC = data.BC_order
    w = np.array(data.order['work_time'])
    T = 30
    B = 100000
    equ_info_ = data.equ_info
    OE = data.order_equ
    OO = data.order_order
    BB = data.BB_order
    RB = 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 = key
            model.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.1
    model.Params.NoRelHeurTime = 30
    model.Params.TimeLimit = 180
    # model.Params.Heuristics = 0.7

    # Optimize model
    model.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()

结果展示

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