Gurobi 生产计划调度学习案例（含代码实现） （生产切换、装配计划）

 作者：help@gurobi.cn

一、前言

利用数学规划建立和解决生产计划调度问题，已经是运筹学成熟应用领域之一。各种书籍、案例、最佳实践不胜枚举。

大部分数学规划方法在建立生产计划和调度模型时，假设以某个自然的时间间隔为建模时间单位，例如小时、日、天、周或者班次等。在建模时间单位内统计可用工时、初始库存和原料等，再优化每个时间单位内的加工种类和数量、工时分配、人员安排等。如果有一定次序关系，例如工序或者BOM等，次序隐含在以建模时间单位为最小单元的计划中。用到的设备和模具的产能是否满足一般是以累积工时的方式进行判断，也就是我们所说的“标量加和”，而不是“矢量加和“。这种建模思路比较适合件小量大场合。

二、案例场景

本文中的二个案例，也是生产计划和调度常见的二个场景。一个涉及到产线切换，一个涉及到多模式装配计划。为了方便理解，案例进行了简化，但保留了现实问题中的特点和重点。具体如下

（一）案例一：生产切换

场景：

（1）    二条产线 L1, L2，产能为 L1:24, L2:24

（2）    四种产品：A1, A2, B1, B2, 需求量：A1: 14, A2: 10, B1:12, B2: 12

（3）    四种产品切换的成本是

	A1	A2	B1	B2
A1	0	1	4	4
A2	1	0	4	4
B1	4	4	0	1
B2	4	4	1	0

（4）初始条件：上个班次最后的加工状态是 L1:A1, L2:A2

目标：

如何分配产品到产线上，在完成需求量的同时，切换成本最小

（二）案例二：装配计划

场景：

（1）多个客户装配订单：O1,O2,O3,O4

（2）每个订单有三道工序 [1,2,3]，处理时间不同。

	工序1	工序2	工序3
O1	10	20	30
O2	40	30	20
O3	20	20	30
O4	40	40	30

（3）加工资源有 R1, R2，模具有M1,M2

（4）每道工序需要的加工资源和模具也不相同。这个范例假设每道工序都需要配备一台资源和一个模具才能开始加工。

（5）每道工序的资源和模具允许存在多种选择。例如对于加工资源

	工序	资源	可选？
O1	1	R1	Y
O1	1	R2	Y
O1	2	R1	Y
O1	2	R2	N
（完整数据请参考后续代码）

O1 的工序1既可以选择R1，也可以选择R2，但O1的工序2只能选择R1，不能选择R2，以此类推。模具的选择也类似。

目标：

优化目标是完成全部订单的时间最短，这样可以最大化利用全部资源和模具。

三、建模思路

以上二个案例如果沿用传统的数学规划建模方式，那么可能需要

（1）对计划周期进行离散化，确定建模时间单位。例如以10分钟、20分钟、60分钟等作为时间单位，将所有时长变成建模时间单位的整倍数。

（2）在每个建模时间单位上，判断每个设备或者模具的占用率，确保不超过100%

（3）每道工序或者加工需要确保在多个建模时间单位上连续加工，一般不允许出现中间有中断。

（4）加工次序将变成如何找到适合某道工序的时间起始位置。

但离散化带来的负面影响也是非常明显，会造成：

（1）决策变量急剧增加，而且时间单位选择不当会造成优化结果不准确。在生产计划调度建模中，时间往往是最大的敌人。

（2）同时需要引入大量辅助变量，来判断生产的连续性，不重叠性，产生切换与否等情况。

（3）优化难度增加，严重影响优化效率。

对于此类问题，本文范例中引入了“虚拟时段”（slot）的概念。就是将计划周期上设置很多虚拟时段，这些虚拟时段有几个特点：

（1）每个时段内只能容纳单个任务。这个单个任务可以是一道工序，或者一次产品加工，或者一次资源分配

（2）每个时段的起讫时间是可变的，各不相同，但不重叠，也不必连续。

（3）时段的主要作用是为了给各种任务指定合理的次序。

（4）虚拟时段的数量需要涵盖可能发生的总任务数量。例如案例一中，每条产线最多加工4种产品，因此虚拟时段数量为4；案例二中，每个资源（理论上）可以被4个订单中的任何工序占用，因此时段数量选择为订单数乘以工序数。当然根据具体场景，可以灵活调整。

虚拟时段的设定，避免了时间离散化。在保留次序选择的同时，也允许时间连续化，降低了变量数量。

四、示范代码

 示范代码是基于Gurobi Python 语言，用户可在后台回复「Gurobi生产调度」获取源文件。

案例一的代码
# -*- coding: utf-8 -*-
"""
@author: help@gurobi.cn
"""

import gurobipy as gp
from gurobipy import *


###############################   输入   ###################################

# 产品
PRODUCTS = ['A1','A2','B1','B2']

# 产线
LINES = ['L1', 'L2']

# 上一次生产中最后的产品
LAST_PRODUCTION = {
    'L1':'A1',
    'L2':'A2',
    }

# 产品需求
DEMAND = {
    'A1':14,
    'A2':10,
    'B1':12,
    'B2':12,
    }

# 切换成本
CHANGEOVER_COST = {
    ('A1', 'A1'):0,
    ('A1', 'A2'):1,
    ('A1', 'B1'):4,
    ('A1', 'B2'):4,
    ('A2', 'A1'):1,
    ('A2', 'A2'):0,
    ('A2', 'B1'):4,
    ('A2', 'B2'):4,
    ('B1', 'A1'):4,
    ('B1', 'A2'):4,
    ('B1', 'B1'):0,
    ('B1', 'B2'):1,
    ('B2', 'A1'):4,
    ('B2', 'A2'):4,
    ('B2', 'B1'):1,
    ('B2', 'B2'):0
    }

# 产线产能
LINE_CAPACITY = 24

# 时段数量
SLOTS = range(len(PRODUCTS))


###############################   模型   ###################################

# 创建模型
model = gp.Model('production_scheduling_with_changeover')

# 变量：每个产线每个产品在每个时段的加工数量
quantity = model.addVars(LINES, SLOTS, PRODUCTS, vtype = GRB.INTEGER, name="qty")

# 变量：根据加工数量来判断每个产线的每个产品的每个时段是否被占用
isBusy = model.addVars(LINES, SLOTS, PRODUCTS, vtype = GRB.BINARY, name="isBusy")

# 变量：每个产线每个时段是否被任何一个产品占用
slotBusy = model.addVars(LINES, SLOTS, vtype = GRB.BINARY, name="slotBusy")

#变量：每个产线每个时段从上时段产品到本时段产品的切换
changeOver =  model.addVars(LINES, SLOTS, PRODUCTS, PRODUCTS, vtype = GRB.BINARY, name="changeOver")


# 1. 满足需求
for product in PRODUCTS:
    model.addConstr(((sum(quantity[line, n, product] for line in LINES for n in SLOTS)) == DEMAND[product]), name = 'meet_demand' + '_' + product)

# 2. 不超过产能
for line in LINES:
    model.addConstr(((sum(quantity[line, n, product] for product in PRODUCTS for n in SLOTS)) <= LINE_CAPACITY), name = 'line_capacity'+'_'+line)


# 3. 建立 isBusy 和 quantity 之间的标志关系
for line in LINES:
    for product in PRODUCTS:
        for n in SLOTS:
            model.addGenConstrIndicator(isBusy[line, n, product], 0, quantity[line, n, product] == 0)
            model.addGenConstrIndicator(isBusy[line, n, product], 1, quantity[line, n, product] >= 1)



#4. 每个产线的每个时段只能允许最多一个产品
for line in LINES:
    for n in SLOTS:
        model.addConstr((sum(isBusy[line, n, product] for product in PRODUCTS) <= 1), name = 'One_Product_Per_Slot' + '_' + line + '_' + str(n))

#5. 每个产线上每个产品只能出现在最多一个时段里
for line in LINES:
    for product in PRODUCTS:
        model.addConstr((sum(isBusy[line, n, product] for n in SLOTS) <= 1), name = 'One_Product_Per_line'+'_'+line+'_'+product)


#6. 统计每个时段被占用情况，不允许出现前面时段没有生产，后面时段有生产的情况
for line in LINES:
    for n in SLOTS:
        model.addConstr(slotBusy[line, n] == max_([isBusy[line,n,product] for product in PRODUCTS]), name='slotbusy'+'_'+line+'_'+str(n))

for line in LINES:
    for n in SLOTS[1:]:
        model.addConstr(slotBusy[line, n-1] >= slotBusy[line, n], name='slotbusyincrease'+'_'+line+'_'+str(n))

#7. 统计每个时段的切换情况
for line in LINES:
    for product in PRODUCTS:
        if product == LAST_PRODUCTION[line]:
            model.addConstr(changeOver[line,0,LAST_PRODUCTION[line],product] == 0, name='changeover'+'_'+line+'_'+LAST_PRODUCTION[line]+'_'+product)
        else:
            model.addConstr(changeOver[line,0,LAST_PRODUCTION[line],product] == isBusy[line,0,product], name='changeover'+'_'+line+'_'+LAST_PRODUCTION[line]+'_'+product)

for line in LINES:
    for n in SLOTS[1:]:
        for p1 in PRODUCTS:
            for p2 in PRODUCTS:
                if p1 == p2:
                    model.addConstr(changeOver[line,n,p1,p2] == 0, name='changeover'+'_'+line+'_'+p1+'_'+p2)
                else:
                    model.addConstr(changeOver[line,n,p1,p2] == and_([isBusy[line,n-1,p1],isBusy[line,n,p2]]),name='changeover'+'_'+line+'_'+p1+'_'+p2)


# 目标
obj = sum(changeOver[line,n,p1,p2]*CHANGEOVER_COST[p1, p2] for p1 in PRODUCTS for p2 in PRODUCTS for line in LINES for n in SLOTS)

# 最小化切换成本
model.setObjective(obj, GRB.MINIMIZE)

# 输出到 LP 文件
model.write("changeover.lp")

# 优化
model.optimize()



print('\n\n###############################   输出结果   ######################################\n')
print('总切换成本：'+ '%3d' % model.objval)
print('生产计划：')
tableStr=''
for n in SLOTS:
    tableStr = tableStr + '%18s'% n
print(tableStr)

for line in LINES:
    tableStr1 = line + '%3s'% LAST_PRODUCTION[line] +'    '
    tableStr2 = '         '
    for n in SLOTS:
        for p in PRODUCTS:
            tableStr1 = tableStr1 + '%3s' % p + ' '
            tableStr2 = tableStr2 + '%3d' % quantity[line,n,p].x + ' '
        tableStr1 = tableStr1 + ' | '
        tableStr2 = tableStr2 + ' | '
    print(tableStr1)
    print(tableStr2)

print('切换成本：')
for line in LINES:
    for n in SLOTS:
        for p1 in PRODUCTS:
            for p2 in PRODUCTS:  
                if round(changeOver[line,n,p1,p2].x) == 1:
                    print(line + ' ' + p1 + ' '+ p2 + ' ' + str(CHANGEOVER_COST[p1, p2]))






案例二的代码：
# -*- coding: utf-8 -*-
"""
Make-To-Order 订单装配 APS

@author: help@gurobi.cn


"""

import gurobipy as gp
from gurobipy import *



###############################   输入   ###################################

# 订单
ORDERS = ['O1','O2','O3','O4']

# 工序
STEPS = [1,2,3]

# 设备
RESOURCES = ['R1','R2']

# 模具
MOULDS = ['M1','M2']

# 加工时间,假设为整数
STEP_TIME= {
    ('O1',1): 10,
    ('O1',2): 20,
    ('O1',3): 30,
    ('O2',1): 40,
    ('O2',2): 30,
    ('O2',3): 20,
    ('O3',1): 20,
    ('O3',2): 20,
    ('O3',3): 30,
    ('O4',1): 40,
    ('O4',2): 40,
    ('O4',3): 30,
    }

# 可选设备：1 为可选， 0 为不可选
STEP_RESOURCES= {
    ('O1',1,'R1'): 1,
    ('O1',1,'R2'): 1,
    ('O1',2,'R1'): 1,
    ('O1',2,'R2'): 0,
    ('O1',3,'R1'): 0,
    ('O1',3,'R2'): 1,
    ('O2',1,'R1'): 1,
    ('O2',1,'R2'): 0,
    ('O2',2,'R1'): 1,
    ('O2',2,'R2'): 0,
    ('O2',3,'R1'): 0,
    ('O2',3,'R2'): 1,
    ('O3',1,'R1'): 0,
    ('O3',1,'R2'): 1,
    ('O3',2,'R1'): 0,
    ('O3',2,'R2'): 1,
    ('O3',3,'R1'): 1,
    ('O3',3,'R2'): 0,
    ('O4',1,'R1'): 1,
    ('O4',1,'R2'): 1,
    ('O4',2,'R1'): 0,
    ('O4',2,'R2'): 1,
    ('O4',3,'R1'): 0,
    ('O4',3,'R2'): 1,

    }

# 可选模具：1 为可选，0为不可选
STEP_MOULDS = {
    ('O1',1,'M1'): 1,
    ('O1',1,'M2'): 0,
    ('O1',2,'M1'): 1,
    ('O1',2,'M2'): 0,
    ('O1',3,'M1'): 1,
    ('O1',3,'M2'): 0,
    ('O2',1,'M1'): 0,
    ('O2',1,'M2'): 1,
    ('O2',2,'M1'): 0,
    ('O2',2,'M2'): 1,
    ('O2',3,'M1'): 0,
    ('O2',3,'M2'): 1,
    ('O3',1,'M1'): 0,
    ('O3',1,'M2'): 1,
    ('O3',2,'M1'): 0,
    ('O3',2,'M2'): 1,
    ('O3',3,'M1'): 1,
    ('O3',3,'M2'): 0,
    ('O4',1,'M1'): 1,
    ('O4',1,'M2'): 0,
    ('O4',2,'M1'): 0,
    ('O4',2,'M2'): 1,
    ('O4',3,'M1'): 0,
    ('O4',3,'M2'): 1,

    }


# 各种统计数据
nORDERS = len(ORDERS)
nSTEPS = len(STEPS)
nRESOURCES= len(RESOURCES)
nMOULDS = len(MOULDS)

M = 1000

# 时段数量
SLOTS = range(nORDERS * nSTEPS)


###############################   模型   ###################################

# 创建模型
model = gp.Model('AssemblyAPS')

# 变量： 每个订单每道工序的起始时间和终止时间
start =  model.addVars(ORDERS, STEPS, vtype = GRB.INTEGER, name="start")
end =  model.addVars(ORDERS, STEPS, vtype = GRB.INTEGER, name="end")

# 变量：每个设备/模具 每个时段上对应的每个产品每道工序起始时间和终止时间
startR =  model.addVars(RESOURCES, SLOTS, ORDERS, STEPS, vtype = GRB.INTEGER, name="startR")
startM =  model.addVars(MOULDS, SLOTS, ORDERS, STEPS, vtype = GRB.INTEGER, name="startM")
endR =  model.addVars(RESOURCES, SLOTS, ORDERS, STEPS, vtype = GRB.INTEGER, name="endR")
endM =  model.addVars(MOULDS, SLOTS, ORDERS, STEPS, vtype = GRB.INTEGER, name="endM")

# 变量：每个订单每道工序使用的设备
useResource = model.addVars(RESOURCES, SLOTS, ORDERS, STEPS, vtype = GRB.BINARY, name="useResource")

# 变量：每个订单每道工序使用的模具
useMould = model.addVars(MOULDS, SLOTS, ORDERS, STEPS, vtype = GRB.BINARY, name="useMould")

# 变量：每个设备/模具在每个时段的起始时间和终止时间
rSlotStartTime = model.addVars(RESOURCES, SLOTS, name="rSlotStartTime")
mSlotStartTime = model.addVars(MOULDS, SLOTS, name="mSlotStartTime")
rSlotEndTime = model.addVars(RESOURCES, SLOTS, name="rSlotEndTime")
mSlotEndTime = model.addVars(MOULDS, SLOTS, name="mSlotEndTime")

# 总时长
timeSpan = model.addVar(vtype = GRB.INTEGER, name="timeSpan")

#1. 订单每道工序起始时间不能早于前道工序结束时间
for step in STEPS:
    for order in ORDERS:
        model.addConstr(end[order, step] == start[order, step] + STEP_TIME[order, step])

for step in STEPS[1:]:
    for order in ORDERS:
        model.addConstr(start[order,step] >= end[order, step-1])

#2. 满足设备需求
for order in ORDERS:
    for step in STEPS:
        model.addConstr(sum(useResource[resource, slot, order, step] for resource in RESOURCES for slot in SLOTS) == 1 )
        model.addConstrs(sum(useResource[resource, slot, order, step] for slot in SLOTS) <= STEP_RESOURCES[order, step, resource] for resource in RESOURCES )

#3. 满足模具需求
for order in ORDERS:
    for step in STEPS:
        model.addConstr(sum(useMould[mould, slot, order, step] for mould in MOULDS for slot in SLOTS) == 1 )
        model.addConstrs(sum(useMould[mould, slot, order, step] for slot in SLOTS) <= STEP_MOULDS[order, step, mould] for mould in MOULDS )

#4. 每个设备/模具每个时段中只能分配被一个产品的一道工序占用
model.addConstrs(sum(useResource[resource, slot, order, step] for order in ORDERS for step in STEPS) <= 1 for resource in RESOURCES for slot in SLOTS)
model.addConstrs(sum(useMould[mould, slot, order, step] for order in ORDERS for step in STEPS) <= 1 for mould in MOULDS for slot in SLOTS)


#5. 每个设备每个时段在每个产品和工序的起始时间终止时间
for resource in RESOURCES:
    for slot in SLOTS:
        for order in ORDERS:
            for step in STEPS:
                model.addConstr(startR[resource, slot, order, step] <= start[order, step] + (1-useResource[resource, slot, order, step])* M)      
                model.addConstr(startR[resource, slot, order, step] >= start[order, step] - (1-useResource[resource, slot, order, step])* M) 
                model.addConstr(startR[resource, slot, order, step] <= useResource[resource, slot, order, step]* M) 
                model.addConstr(endR[resource, slot, order, step] <= start[order, step] + STEP_TIME[order, step] + (1-useResource[resource, slot, order, step])* M)      
                model.addConstr(endR[resource, slot, order, step] >= start[order, step] + STEP_TIME[order, step] - (1-useResource[resource, slot, order, step])* M) 
                model.addConstr(endR[resource, slot, order, step] <= useResource[resource, slot, order, step]* M)

#6. 每个设备每个时段的起始时间和终止时间
for resource in RESOURCES:
    for slot in SLOTS:
        model.addConstr(rSlotStartTime[resource, slot] == sum(startR[resource, slot, order, step] for order in ORDERS for step in STEPS))
        model.addConstr(rSlotEndTime[resource, slot] == sum(endR[resource, slot, order, step] for order in ORDERS for step in STEPS))

#7. 每个模具每个时段在每个产品和工序的起始时间终止时间
for mould in MOULDS:
    for slot in SLOTS:
        for order in ORDERS:
            for step in STEPS:
                model.addConstr(startM[mould, slot, order, step] <= start[order, step] + (1-useMould[mould, slot, order, step])* M)      
                model.addConstr(startM[mould, slot, order, step] >= start[order, step] - (1-useMould[mould, slot, order, step])* M) 
                model.addConstr(startM[mould, slot, order, step] <= useMould[mould, slot, order, step]* M)
                model.addConstr(endM[mould, slot, order, step] <= start[order, step] + STEP_TIME[order, step] + (1-useMould[mould, slot, order, step])* M)      
                model.addConstr(endM[mould, slot, order, step] >= start[order, step] + STEP_TIME[order, step] - (1-useMould[mould, slot, order, step])* M) 
                model.addConstr(endM[mould, slot, order, step] <= useMould[mould, slot, order, step]* M)

#8. 每个模具每个时段的起始时间和终止时间
for mould in MOULDS:
    for slot in SLOTS:
        model.addConstr(mSlotStartTime[mould, slot] == sum(startM[mould, slot, order, step] for order in ORDERS for step in STEPS))
        model.addConstr(mSlotEndTime[mould, slot] == sum(endM[mould, slot, order, step] for order in ORDERS for step in STEPS))

#9，起点时间限制为0
for resource in RESOURCES:
    model.addConstr(rSlotStartTime[resource, 0] == 0)

for mould in MOULDS:
    model.addConstr(mSlotStartTime[mould, 0] == 0)


#10. 设备和模具的每个时段的起始时间不能早于前个时段的终止时间
for resource in RESOURCES:
    for slot in SLOTS[1:]:
          model.addConstr(rSlotEndTime[resource, slot-1] <= rSlotStartTime[resource, slot])

for mould in MOULDS:
    for slot in SLOTS[1:]:
          model.addConstr(mSlotEndTime[mould, slot-1] <= mSlotStartTime[mould, slot])                        


#11. 定义 timespan 为最晚完成订单的终止时间
model.addConstr(timeSpan == max_([end[order, step] for order in ORDERS for step in STEPS]))

# 最小化完成时间
model.setObjective(timeSpan, GRB.MINIMIZE)

# 输出到 LP 文件
model.write("assemblyAPS.lp")

# 优化
model.optimize()



print('\n\n###############################   输出结果   ######################################\n')
print('总时长：'+ '%3d' % model.objval)

print('订单工序计划：')
for order in ORDERS:
    for step in STEPS:
        string = order + '\t' + str(step) + '%10s'% str(start[order,step].x) +'%10s'% str(end[order,step].x) +'\t'
        stop = 0
        for resource in RESOURCES:
            for slot in SLOTS:
                if abs(useResource[resource, slot, order, step].x -1) <= 0.01:
                    string = string + resource + '\t'
                    stop = 1
                    break
            if stop == 1:
                break
        stop = 0
        for mould in MOULDS:
            for slot in SLOTS:
                if abs(useMould[mould, slot, order, step].x -1) <= 0.01:
                    string = string + mould + '\t'
                    stop = 1
                    break
            if stop == 1:
                break

        print(string)           

print('\n')

print('设备使用计划:')

for resource in RESOURCES:
    for slot in SLOTS:
        string= ''
        if rSlotEndTime[resource, slot].x > rSlotStartTime[resource, slot].x:
            string = resource + '%10s'% str(rSlotStartTime[resource, slot].x) + '%10s'% str(rSlotEndTime[resource, slot].x) + '\t'
            stop = 0
            for order in ORDERS:
                for step in STEPS:
                    if abs(useResource[resource, slot, order, step].x -1) <= 0.01:
                        string = string + order + '\t' + str(step)
                        stop = 1
                        break
                if stop == 1:
                    break
            print(string)

print('\n')
print('模具使用计划:')

for mould in MOULDS:
    for slot in SLOTS:
        string= ''
        if mSlotEndTime[mould, slot].x > mSlotStartTime[mould, slot].x:
            string = mould + '%10s'% str(mSlotStartTime[mould, slot].x) + '%10s'% str(mSlotEndTime[mould, slot].x) + '\t'
            stop = 0
            for order in ORDERS:
                for step in STEPS:
                    if abs(useMould[mould, slot, order, step].x -1) <= 0.01:
                        string = string + order + '\t' + str(step)
                        stop = 1
                        break
                if stop == 1:
                    break
            print(string)