【作业车间调度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调用PuLP建模求解

1. 引入PuLP线性规划库

PuLP 是用 Python 编写的 LP 建模器。PuLP 可以生成 MPS 或 LP 文件并调用其他的求解器进行求解。这里我们介绍用 PuLP 库本身求解该 JSP 问题。

import pulp

2. 实例化问题（框架）

通过 PuLP 库的 LpProblem 类实例化一个名为 MIN_makespan 的问题，且该问题是最小化问题。

problem = pulp.LpProblem('MIN_makespan',pulp.LpMinimize)
# 最大化问题 pulp.LpMaximize

3. 创建决策变量

PuLP 预设的变量的下限是负无穷，上限是正无穷，变量的类型 category 只能是 （“Continuous”，“Integer”，“Binary”）

由于这里的每个工件都会在每个机器上有且仅加工一次，因此可以基于工件视角，也可以基于机器视角创建变量，这里以机器视角创建变量如下：

1. $C_{m,j}$：第 $m$ 台机器加工 $j$ 工件的开始时间；
2. $b_{m,j1,j2}$：在第 $m$ 台机器上，$j1$ 工件是否在 $j2$ 工件之前加工，是为1，否为0；由于每台机器均会加工一次所有的工件，因此 $b_{m,j1,j2}+bm,j2,j1=1$，即 $b_{m,j1,j2}=0$ 时，则说明在第 $m$ 台机器上，$j2$ 工件在 $j1$ 工件之前加工；
3. $C_{max}$：最大完工时间

3.1 单个变量创建

这里通过循环一个个地创建变量。

C = {}
for m in range(6):
	for j in range(6):
		C[m,j] = pulp.LpVariable(name = 'start_time',
                                lowBound = 0,
                                upBound = None,	# 不约束上界
                                cat = 'Continuous')	# 默认是连续变量

3.2 矩阵变量创建

这里我们用矩阵的形式创建变量，通过 LpVariable.dicts 方法创建。

# 创建变量 C
C = pulp.LpVariable.dicts("start_time",
	((m, j) for m in range(6) for j in range(6)),
	lowBound=0, upBound = None, cat='Continuous')
# 创建变量 y
b = pulp.LpVariable.dicts("binary_var",
	((m,j1,j2) for m in range(6) for j1 in range(6) for j2 in range(6)
	if j1 != j2), lowBound=0, upBound = None, cat='Binary')
# 创建变量 C_max
Cmax = pulp.LpVariable('Cmax',lowBound = 0, cat='Continuous')

4. 创建约束条件

在约束条件中，需要考虑到工件工序的加工顺序，以及同一台设备的加工工序不能重叠。因此需要引入超参数数据，如下：

# 工序的加工序列
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]}

约束1：每个工件需要按照既定的加工顺序进行加工，即工件的后续工序的开工时间不早于前序工序的开工时间。

for j in range(6):
    for m_index in range(5):
       seq_pre = order[j][m_index]-1
       seq_back = order[j][m_index+1]-1
       model += (C[seq_back,j] - C[seq_pre,j]) >= process_time[j][seq_pre]

约束2：机器上的加工工序不能重叠。

for m in range(6):
    for j1 in range(6):
        for j2 in range(6):
            if j1!=j2:
                model += (C[m,j2] - C[m,j1]) >= (process_time[j1][m] - 100*(1-b[m,j1,j2]))
                model += (C[m,j1] - C[m,j2]) >= (process_time[j2][m] - 100*b[m,j1,j2])

约束3：最大完工时间不小于每个工序的完工时间。

for m in range(6):
    for j in range(6):
       model += (Cmax - C[m,j]) >= process_time[j][m]

5. 添加目标函数

上述约束条件中，model 添加的都是逻辑表达式（返回布尔变量），而目标函数是直接添加变量完成的。本问题的目标为最小化全部工件的完工时间 Makespan，代码如下：

# 目标函数
model += Cmax

6. 求解模型

solve() 方法支持调用其他的求解器。print(LpStatus[status]) 可以打印模型求解状态，print(status) 可以打印求解的状态码。

状态 LpStatus[status]	状态码 status
“Optimal”	1
“Not Solve”	0
“Infeasible”	-1
“Unbounded”	-2
“Undefined”	-3

# 求解模型
model.solve()
pulp.LpStatus[model.status]		

求解日志：

Result - Optimal solution found

Objective value:                56.00000000
Enumerated nodes:               351
Total iterations:               33836
Time (CPU seconds):             4.85
Time (Wallclock seconds):       4.85

Option for printingOptions changed from normal to all
Total time (CPU seconds):       4.86   (Wallclock seconds):       4.86

返回变量及目标函数值：

## 打印变量值
for var in C:
    var_value = C[var].varValue
    print( var[0],'-',var[1] ,var_value)
## 打印目标函数值
total_cost = pulp.value(model.objective)
print ('C_max:',total_cost)

返回结果值：

0 - 0 24.0
0 - 1 34.0
0 - 2 29.0
0 - 3 16.0
0 - 4 42.0
0 - 5 21.0
1 - 0 25.0
1 - 1 3.0
1 - 2 34.0
1 - 3 11.0
1 - 4 31.0
1 - 5 0.0
2 - 0 18.0
2 - 1 8.0
2 - 2 0.0
2 - 3 30.0
2 - 4 24.0
2 - 5 35.0
3 - 0 28.0
3 - 1 42.0
3 - 2 13.0
3 - 3 35.0
3 - 4 52.0
3 - 5 3.0
4 - 0 47.0
4 - 1 19.0
4 - 2 46.0
4 - 3 38.0
4 - 4 35.0
4 - 5 29.0
5 - 0 39.0
5 - 1 29.0
5 - 2 22.0
5 - 3 46.0
5 - 4 38.0
5 - 5 13.0
min cost: 56.0

基于 plotly 展示甘特图

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

1. 引入 plotly 库及相关库

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

import datetime
import plotly.figure_factory as ff

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

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

j_record={}
for var in C:
    var_value = C[var].varValue
    start_time=str(datetime.timedelta(seconds=var_value))
    end_time=str(datetime.timedelta(seconds=var_value+process_time[var[1]][var[0]]))
    j_record[(var[0],var[1])]=[start_time,end_time]

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

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

gantt_data = []
for m in range(6):
    for j in range(6):
        gantt_data.append(
            dict(Task=f"Machine_{m}", Start=f"2018-11-26 {str(j_record[m, j][0])}", Finish=f"2018-11-26 {str(j_record[m, j][1])}", Resource=f"Job_{j}")
        )

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. ↩︎