【作业车间调度JSP】通过python调用OR-Tools的区间变量建模求解

---

本文的“问题描述”继承自前文《【作业车间调度JSP】通过python调用PuLP线性规划库求解》

问题描述

JSP（Job Shop Scheduling Problem）问题是指车间作业调度问题，是生产调度领域中的经典问题之一。JSP问题涉及到一组工件（Jobs）和一组机器（Machines），每个工件都有一系列工序需要在不同的机器上完成，每个工序有一个特定的加工时间。通过确定每个工件的工序顺序和在每台机器上的开始时间，以最优化既定的目标。这个问题的复杂性（NP-Hard问题）和实际应用的广泛性使得研究和优化JSP一直是运筹学和制造业领域的热点。

现在考虑这样一个问题，有 $j$ 个工件，$m$ 台机器，每个工件需按给定的顺序在不同的机器上加工（表1），且同一机器加工不同工件的工序时间是不同的（表2）。下表的数据来源于文献¹中的测试数据，其中 $j=6，m=6$。最终优化目标是最小化全部工件的完工时间。

表1：工件的加工顺序

工件编号	加工顺序（机器编号）
1	[3,1,2,4,6,5]
2	[2,3,5,6,1,4]
3	[3,4,6,1,2,5]
4	[2,1,3,4,5,6]
5	[3,2,5,6,1,4]
6	[2,4,6,1,5,3]

表2：工件在不同机器上的加工时间

工件编号	机器1 ~ 机器6
1	[1,3,6,7,3,6]
2	[8,5,10,10,10,4]
3	[5,4,8,9,1,7]
4	[5,5,5,3,8,9]
5	[9,3,5,4,3,1]
6	[3,3,9,10,4,1]

假设条件

对于 JSP 问题常常存在自然且约定俗成的假设，包括：

1. 机器的加工能力是有限的，这里每台机器一次只能做一个工作
2. 工件的工序需按给定顺序进行加工
3. 机器一旦开工则不允许中断

Python调用OR-Tools建模求解

对于JSP问题而言，最基本的约束有两块：（1）工件的工序需按顺序加工，前序工序加工完才能加工后续工序；（2）机器加工的不重叠性，即机器在同一时间只能完成一个加工任务。

特别的，OR-Tools具有关于区间变量的特定约束，即去重叠。因此，本文基于OR-Tools的区间变量进行建模求解，验证区间变量在这类带顺序的模型（VRP/JSP…）中的效果。

1. 引入相关库

OR-Tools的区间变量在它的 CP-SAT 求解器当中，因此这里需要引入 cp_model 模块。
这里我们为了方便存储，引入collections 模块的轻量级数据结构 —— 命名元组“namedtuple”，以及不用判断字典键的默认值字典“defaultdict”作为辅助工具（当然，也可以通过其他方式替代这些工具）。具体如下：

import collections
from ortools.sat.python import cp_model

2. 声明问题的模型

OR-Tools 的 cp_model 在建模之前需要先声明问题的模型。

model = cp_model.CpModel()

3. 创建区间变量

在这里，区间变量指每个任务的加工时间段，即机器被某道工序占用的时间段，当工序在机器上的加工时间固定时，这个区间变量的长度 duration 也固定，进一步地，当工序的开始加工时间固定，则相对应的结束时间也固定。

前文《【作业车间调度JSP】通过python调用PuLP线性规划库求解》的实验数据如下：

# 工序的加工序列
order = {0: [3,1,2,4,6,5],
         1: [2,3,5,6,1,4], 
         2: [3,4,6,1,2,5],
         3: [2,1,3,4,5,6],
         4: [3,2,5,6,1,4],
         5: [2,4,6,1,5,3]}
# 工序在 1~6 机器上的加工时间
process_time = {0: [1,3,6,7,3,6],
         1: [8,5,10,10,10,4], 
         2: [5,4,8,9,1,7],
         3: [5,5,5,3,8,9],
         4: [9,3,5,4,3,1],
         5: [3,3,9,10,4,1]}

对每个工序都要建立一个区间变量，而区间变量的基本元素包括，工序的开始加工时间（变量），结束加工时间（变量），加工时长（常量）。为了尽量减少可行解空间大小，开始时间和结束时间的上界尽可能紧凑，这里我们取所有工序的加工时长之和作为时间变量的上界。

# 命名元组：存储变量信息（相当于实例化一道工序）
task_type = collections.namedtuple("task_type", "start end interval")
# 存储所有工序信息
all_tasks = {}
# 将工序信息按可上机器（一台）进行存储
machine_to_intervals = collections.defaultdict(list)

# 时间变量的上界
horizon = sum(p_time for job in process_time for p_time in process_time[job])

for job_id, job in order.items():
    for task_id, machine in enumerate(job):
        duration = process_time[job_id][machine - 1]
        suffix = f"_{job_id}_{task_id}"

        start_var = model.NewIntVar(0, horizon, "start" + suffix)
        end_var = model.NewIntVar(0, horizon, "end" + suffix)
        interval_var = model.NewIntervalVar(
            start=start_var, size=duration, end=end_var, name="interval" + suffix
        )
        all_tasks[job_id, task_id] = task_type(
            start=start_var, end=end_var, interval=interval_var
        )
        machine_to_intervals[machine].append(interval_var)

上述创建变量的代码中，先创建了 start_var 和 end_var 变量，然后将这两个变量传入模型的函数NewIntervalVar以创建区间变量。其中，当区间变量内指定了end后，则默认end变量的值与start变量的值的差距为duration。

然后将创建的变量存储进 all_task 字典，并按工序的可上机器存储到 machine_to_intervals 字典。

注意：OR-Tools的CP-SAT求解器支持的基础变量类型只有整数型变量（包括布尔变量），因此对时间变量时创建为整数变量，可通过系数的方式来控制变量精度。

4. 创建约束条件

在OR-Tools中，有专门针对区间变量的约束AddNoOverlap，即指定若干区间变量之间彼此不重叠，但无先后关系保证，这类约束成为不相容约束（Disjunctive constraints），在JSP问题中体现在，机器在同一时间只能加工一个任务（暂不考虑批量加工），此时需要约束同一机器上的区间变量彼此之间不重叠。

约束1：机器上的加工工序不能重叠（不相容约束）

for machine in machine_to_intervals:
	model.AddNoOverlap(machine_to_intervals[machine])

约束2：每个工件内的待加工工序具有先后顺序。

for job_id, job in order.items():
	for task_id in range(len(job) - 1):
		model.Add(all_tasks[job_id, task_id + 1].start >= all_tasks[job_id, task_id].end)

约束3：最大完工时间等于所有工序完工时间的最大值。

obj_var = model.NewIntVar(0, horizon, "makespan")
model.AddMaxEquality(obj_var, [all_tasks[index].end for index in all_tasks])
model.Minimize(obj_var)

5. 求解模型

OR-Tools的CP-SAT需要模型先创建求解器，再调用Solve()函数进行求解。

solver = cp_model.CpSolver()
status = solver.Solve(model)

CP-SAT求解器的 solve() 方法支持写 callback 规则，即求解器在求解过程中，有条件地触发相应的动作（之后安排文章介绍）。这里的solve()函数求解后会返回一个状态码，状态码相应的含义如下：

状态	状态码 status
未知 cp_model.UNKNOWN	0
模型无效 cp_model.MODEL_INVALID	1
找到可行解 cp_model.FEASIBLE	2
证明了不可行 cp_model.INFEASIBLE	3
达到最优解 cp_model.OPTIMAL	4

打印求解结果：

print(f"Optimal Schedule Length: {solver.ObjectiveValue()}")
print("\nStatistics")
print(f"  - conflicts: {solver.NumConflicts()}")
print(f"  - branches : {solver.NumBranches()}")
print(f"  - wall time: {solver.WallTime()}s")

Optimal Schedule Length: 56.0

Statistics
  - conflicts: 14
  - branches : 87
  - wall time: 0.0580873s

值得注意的是：在前文用PuLP求解相同问题时，求解时间为4.86s，而在CP-SAT中求解时间为0.06s，可知在JSP问题中，CP-SAT的区间变量具有极大的性能优势。

通过求解器的方法 solver.Value( ) 获取变量信息，打印出最终的方案信息：

assigned_jobs = collections.defaultdict(list)
    if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
        print("Solution:")
        for job_id, job in order.items():
            for task_id, machine in enumerate(job):
                assigned_jobs[machine].append((job_id, task_id))

        for machine_ in assigned_jobs:
            print(f"机器{machine_}加工工序的开始时间：")
            for task_ in assigned_jobs[machine_]:
                print(f"工件 {task_[0]} 的 {task_[1]} 工序 {solver.Value(all_tasks[task_].start)}")

        # 完工时间
        print(f"Optimal Schedule Length: {solver.ObjectiveValue()}")
    else:
        print("No solution found.")

返回结果值：

器3加工工序的开始时间：
工件 0 的 0 工序 18
工件 1 的 1 工序 8
工件 2 的 0 工序 0
工件 3 的 2 工序 29
工件 4 的 0 工序 24
工件 5 的 5 工序 35
机器1加工工序的开始时间：
工件 0 的 1 工序 24
工件 1 的 4 工序 34
工件 2 的 3 工序 29
工件 3 的 1 工序 18
工件 4 的 4 工序 42
工件 5 的 3 工序 25
机器2加工工序的开始时间：
工件 0 的 2 工序 25
工件 1 的 0 工序 3
工件 2 的 4 工序 35
工件 3 的 0 工序 12
工件 4 的 1 工序 29
工件 5 的 0 工序 0
机器4加工工序的开始时间：
工件 0 的 3 工序 28
工件 1 的 5 工序 42
工件 2 的 1 工序 13
工件 3 的 3 工序 35
工件 4 的 5 工序 52
工件 5 的 1 工序 3
机器6加工工序的开始时间：
工件 0 的 4 工序 36
工件 1 的 3 工序 29
工件 2 的 2 工序 22
工件 3 的 5 工序 46
工件 4 的 3 工序 35
工件 5 的 2 工序 13
机器5加工工序的开始时间：
工件 0 的 5 工序 47
工件 1 的 2 工序 18
工件 2 的 5 工序 46
工件 3 的 4 工序 38
工件 4 的 2 工序 32
工件 5 的 4 工序 28
Optimal Schedule Length: 56.0

基于 plotly 展示甘特图

甘特图是查看车间调度问题方案的最直观形式，基于上面的求解结果画甘特图的步骤如下：

1. 引入 plotly 库及相关库

这里引入了 datetime 库是为了将上述模型的变量值（秒）转化为日期形式：年-月-日 时:分:秒，以方便生成甘特图。

import datetime
import plotly.figure_factory as ff

2. 取出每个工序的开始时间和结束时间

每个工序的开始时间通过模型的变量值取出，而结束时间通过开始时间加上相应的加工时间。datetime.timedelta() 将变量值（秒）转化为时间形式（时:分:秒）。

j_record = {}
for job_id, job in order.items():
    for task_id, machine in enumerate(job):
        start_ = solver.Value(all_tasks[job_id, task_id].start)
        start_time=str(datetime.timedelta(seconds=start_))
        end_time=str(datetime.timedelta(seconds=start_+process_time[job_id][machine-1]))
        j_record[(machine-1, job_id)]=[start_time,end_time]

3. 设置甘特图元素的信息

上一步骤从模型中取出每个工序的开始时间和结束时间，这一步骤将配置甘特图元素信息，即这些工序（条形图）隶属于哪个机器，以及哪些工序是同一 Job。

gantt_data = []
for machine_ in range(6):
    for job_ in range(6):
        gantt_data.append(
            dict(Task=f"Machine_{machine_}", Start=f"2023-12-01 {str(j_record[machine_, job_][0])}", 
                 Finish=f"2023-12-01 {str(j_record[machine_, job_][1])}", Resource=f"Job_{job_}")
        )

fig = ff.create_gantt(gantt_data, index_col='Resource', group_tasks=True, show_colorbar=True, title="Gantt Chart of JSP")

fig.show()

结果图展示：

---

1. J.F.Muth & G.L.Thompson. Industrial Scheduling. Prentice Hall, Englewood Cliffs, New Jersey,1963. ↩︎