作业车间调度问题求解框架：PuLP 求解框架¶

发布于：2021-08-29 | 分类：optimization , mathematics

本文根据作业车间调度问题的数学描述，利用 Python 整数规划框架 PuLP 进行建模，并分别使用默认的CBC求解器及商业求解器Gurobi求解。

PuLP¶

PuLP 是开源基金会 COIN-OR (Computational Infrastructure for Operations Research) 维护的一个线性规划建模工具（LP modeler），内置了 CBC（COIN-OR Branch-and-Cut）求解器，也支持集成各种开源（如GLPK）或者商业（CPLEX、GUROBI等）的线性规划求解器。

作为 Python 库，PuLP 支持 pip 一键安装。

pip install pulp

更多入门材料参考：

数学模型¶

前文给出了作业车间调度问题的数学模型。为了方便接下来的建模描述，改为以工序 开始时间 为决策变量的描述方式：

\begin{align*} &\min\,\, c_{m} \\ \\ &s.t.\quad \begin{cases} s_{ij} \geq 0 & (i,j) \in \mathbb{O} \\ \\ c_m-s_{ij} \geq p_{ij} & (i,j) \in \mathbb{O} \\ \\ s_{ij}-s_{kj} \geq p_{kj} & if \,\, (k,j) \rightarrow (i,j), \,\, i, k \in \mathbb{M}^j \\ \\ s_{ij} - s_{ik} \geq p_{ik} \,\, or \,\, s_{ik}-s_{ij} \geq p_{ij} & (i,j), (i,k) \in \mathbb{O}, \,\, j \neq k \end{cases} \end{align*}

其中两个基本变量：

$s_{ij}$ 表示工序 $(i,j)$ 即作业 $j$ 的某个在机器 $i$ 上进行的工序的开始时间
$p_{ij}$ 表示作业 $j$ 的某个在机器 $i$ 上进行的工序所需的加工时间

四个基本约束：

第一个表明工序的开始时间非负
第二个定义了最大完工时间 $c_m$
第三个是工序序列约束，即同一个作业的相邻工序必须满足先后顺序
第四个是资源约束，即每一台机器在某个时刻只能加工一道工序

我们需要特别关注一下第四式，因为它出现了一个或的约束关系，但线性规划问题的约束式通常都是与也就是且的关系。为了正确实施这个约束，我们需要引入额外的变量 $b_{jk}$ ，表示分配到机器 $m_i$ 的所有工序中任意两者例如工序 $(i,j)$ 和工序 $(i,k)$ 的先后关系：

\begin{align*} b_{jk} = \begin{cases} 0 & if \,\, (i,j) \rightarrow (i,k) \\ \\ 1 & if \,\, (i,k) \rightarrow (i,j) \end{cases} \end{align*}

即， $b_{jk}$ 是一个 $0-1$ 变量，当工序 $(i,j)$ 排在工序 $(i,k)$ 之前取值 0，否则 1。

于是，原来写法中逻辑或的两种情况被等效为 $0-1$ 变量 $b_{jk}$ 的两种可能。进而，改写约束式的思路是合并/统一原来两个分支，一旦确定 $b_{jk}$ 则退化为其中一个分支。这一步还需要一点小技巧，最终约束四被改写为以下两个式子：

\begin{align*} \begin{cases} s_{ij}-s_{ik} \geq p_{ik} - C*(1-b_{jk}) & (i,j), (i,k) \in \mathbb{O}, \,\, j \neq k \\ \\ s_{ik}-s_{ij} \geq p_{ij} - C*b_{jk} & (i,j), (i,k) \in \mathbb{O}, \,\, j \neq k \end{cases} \end{align*}

其中， $C$ 是一个足够大例如满足 $C>\Sigma \, p_{ij}$ 的常数。

验证一下，假设工序 $(i,j)$ 排在工序 $(i,k)$ 之前即 $b_{jk}=0$ ，则第二式反映正确的约束关系，而第一式显然成立即不反应任何有效信息；相反， $b_{jk}=1$ 时，第一式反映目标约束而第二式显然成立。

PuLP 模型¶

使用 PuLP 建模和求解的步骤与 Google OR-Tools 极为类似，但 PuLP 并未提供如 OR-Tools 一般便捷的建立约束的 API。例如，OR-Tools 的区间类型变量配合 NoOverlap 约束可以很方便地实现机器上工序不重叠的约束（即上一节描述的约束四），但 PuLP 需要更基础的数学操作，例如上一节进行的约束转换即为此服务。

PuLP 建模常用的API：

LpProblem 类
```
LpProblem(name='NoName', sense=LpMinimize)
```
其中，sense 参数目标函数是求极大值（LpMaximize）还是极小值（LpMinimize）。

建立约束后调用 solve(solver=None, **kwargs) 求解，默认为 CBC 求解器。
LpVariable 类
```
LpVariable(name, lowBound=None, upBound=None, cat='Continuous', e=None)
```
其中，
- name 表示变量名
- lowBound 和 upBound 表示变量范围下界和上界，默认负无穷到正无穷
- cat 指定变量是离散类型（Integer或Binary），还是连续类型（Continuous）
特别地，当变量个数特别多且满足批量化处理时，可以使用如下方法建立字典格式的变量集合。
```
LpVariable.dicts(name, indexs, lowBound=None, upBound=None, cat='Continuous', indexStart=[])
```
其中，
- name 指定所有变量的前缀
- indexs 是一个列表，将会为其中每一个元素创建一个变量，并建立对应关系

基于 `jsp_framework` 实现¶

本文参考案例，基于 jsp_framework 进行实现和集成，完整代码参考。

基本流程¶

参考 Python建模，自定义一个继承自 JSSolver 的 PuLPSolver 类，然后重点实现 do_solve() 方法。

class PuLPSolver(JSSolver):

    ...

    def do_solve(self):
        '''Solve JSP with PuLP default CBC solver.'''
        # Initialize an empty solution from problem
        solution = self.init_solution()

        # create model
        model, variables = self.__create_model(solution)

        # The problem is solved using PuLP's choice of Solver
        model.solve(pulp.PULP_CBC_CMD(maxSeconds=self.__max_time, msg=1, fracGap=0))
        if model.status!=pulp.LpStatusOptimal: 
            raise JSPException('No feasible solution found.')

        # assign pulp solution back to JSPSolution
        for op, var in variables.items():
            op.update_start_time(var.varValue)        
        self.update_solution(solution) # update solution

注意

当前版本支持 CBC（COIN-OR Branch-and-Cut）、SCIP （开源）和 Gurobi（商业）求解器。其中 CBC 为 Pulp 默认求解器，开箱即用；其余二者需要自行下载安装和授权。

PuLP 建模¶

即上一代码片段 self.__create_model(solution) 的具体展开。

创建问题

# create the model
model = pulp.LpProblem("min_makespan", pulp.LpMinimize)

创建变量

基本变量：每一道工序的开始时间

# (1) start time of each operation
max_time = sum(op.source.duration for op in solution.ops) # upper bound of variables
variables = pulp.LpVariable.dicts(name='start_time', \
                                indexs=solution.ops, \
                                lowBound=0, \
                                upBound=max_time, \
                                cat='Integer')

其中，max_time是所有工序用时之和，即不重叠任何工序的调度结果。显然，这是任意工序开始时间的上限。并且，这个值可以作为问题描述中提及的 足够大的数。

0-1变量：表征同一机器上任意两道工序的前后顺序

# (2) binary variable, i.e. 0 or 1, indicating the sequence of every two 
# operations assigned in same machine
combinations = []
for _, ops in solution.machine_ops.items():
    combinations.extend(pulp.combination(ops, 2))
bin_vars =  pulp.LpVariable.dicts(name='binary_var', \
                                indexs=combinations, \
                                lowBound=0, \
                                cat='Binary')

其中，pulp.combination得到了机器队列工序的两两组合。

创建目标函数：总的加工周期

# objective: makespan / flowtime
s_max = pulp.LpVariable(name='max_start_time', \
                        lowBound=0, \
                        upBound=max_time, \
                        cat='Integer')
model += s_max

注意：此处目标函数仅仅是一个普通的变量，后续添加约束保证其为最大的加工周期

添加约束

最大加工周期约束

# (1) the max start time
for job in solution.job_ops:
    last_op = job.tail_job_op
    model += (s_max-variables[last_op]) >= last_op.source.duration

作业视角的工序顺序约束

# (2) operation sequence inside a job
pre = None
for op in solution.ops:
    if pre and pre.source.job==op.source.job:
        model += (variables[op]-variables[pre]) >= pre.source.duration
    pre = op

机器视角的工序顺序约束

# (3) no overlap for operations assigned to same machine
for op_a, op_b in combinations:
    model += (variables[op_a]-variables[op_b]) >= (op_b.source.duration - max_time*(1-bin_vars[op_a, op_b]))
    model += (variables[op_b]-variables[op_a]) >= (op_a.source.duration - max_time*bin_vars[op_a, op_b])

计算实例¶

最后，求解几个标准问题。与 Google OR-Tools 一样，求解时间上限设定为300秒。

# benchmark.py
from jsp import (JSProblem, BenchMark)
from jsp.solver import PuLPSolver

# problems
names = ['ft06', 'la01', 'ft10', 'swv01', 'la38', \
        'ta31', 'swv12', 'ta42', 'ta54', 'ta70']
problems = [JSProblem(benchmark=name) for name in names]

# solver
solvers = [PuLPSolver(name=name, solver_name=name, max_time=300) \
    for name in ('cbc', 'gurobi')]

# solve 
benchmark = BenchMark(problems=problems, solvers=solvers, num_threads=5)
benchmark.run(show_info=True)

结果如下表所示：

+----+---------+--------+---------------+--------------+----------+---------+-------+
| ID | Problem | Solver | job x machine |   Optimum    | Solution | Error % |  Time |
+----+---------+--------+---------------+--------------+----------+---------+-------+
| 1  |   ft06  |  cbc   |     6 x 6     |      55      |   55.0   |   0.0   |  3.5  |
| 2  |   la01  |  cbc   |     10 x 5    |     666      |  666.0   |   0.0   | 209.0 |
| 3  |   ft10  |  cbc   |    10 x 10    |     930      |  1034.0  |   11.2  | 301.2 |
| 4  |  swv01  |  cbc   |    20 x 10    |     1407     |  2941.0  |  109.0  | 301.2 |
| 5  |   la38  |  cbc   |    15 x 15    |     1196     |  1769.0  |   47.9  | 301.2 |
| 6  |   ta31  |  cbc   |    30 x 15    |     1764     | 13294.0  |  653.6  | 301.1 |
| 7  |  swv12  |  cbc   |    50 x 10    | (2972, 3003) |  9527.0  |  218.9  | 309.9 |
| 8  |   ta42  |  cbc   |    30 x 20    | (1867, 1956) | 23202.0  |  1113.8 | 507.1 |
| 9  |   ta54  |  cbc   |    50 x 15    |     2839     | 26229.0  |  823.9  | 608.0 |
| 10 |   ta70  |  cbc   |    50 x 20    |     2995     | 36851.0  |  1130.4 | 625.6 |
+----+---------+--------+---------------+--------------+----------+---------+-------+

+----+---------+--------+---------------+--------------+----------+---------+-------+
| ID | Problem | Solver | job x machine |   Optimum    | Solution | Error % |  Time |
+----+---------+--------+---------------+--------------+----------+---------+-------+
| 1  |   ft06  | gurobi |     6 x 6     |      55      |   55.0   |   0.0   |  0.7  |
| 2  |   la01  | gurobi |     10 x 5    |     666      |  666.0   |   0.0   |  5.2  |
| 3  |   ft10  | gurobi |    10 x 10    |     930      |  930.0   |   0.0   |  60.2 |
| 4  |  swv01  | gurobi |    20 x 10    |     1407     |  1757.0  |   24.9  | 301.4 |
| 5  |   la38  | gurobi |    15 x 15    |     1196     |  1196.0  |   0.0   | 300.9 |
| 6  |   ta31  | gurobi |    30 x 15    |     1764     |  2142.0  |   21.4  | 304.5 |
| 7  |  swv12  | gurobi |    50 x 10    | (2972, 3003) | unsolved |   n.a.  | 308.9 |
| 8  |   ta42  | gurobi |    30 x 20    | (1867, 1956) |  2457.0  |   28.5  | 362.8 |
| 9  |   ta54  | gurobi |    50 x 15    |     2839     | unsolved |   n.a.  | 607.6 |
| 10 |   ta70  | gurobi |    50 x 20    |     2995     |  4513.0  |   50.7  | 608.2 |
+----+---------+--------+---------------+--------------+----------+---------+-------+

cbc求解器的计算效果并不理想，只有前两个问题即工序总数 50 以内得到了最优解，其余案例误差较大。Gurobi的结果略有改进，前5题结果尚可，后5题或者无法求解或者误差较大，整体不如前文的OR-Tools 约束求解器。

由此表明，对于作业车间调度问题，混合整数规划求解器（如CBC、Gurobi）的求解性能不如 约束求解器（如Google OR-Tools 的CpSolver）。

Keys	Action
`?`	Open this help
`n`	Next page
`p`	Previous page
`s`	Search

作业车间调度问题求解框架：PuLP 求解框架¶

PuLP¶

数学模型¶

PuLP 模型¶

基于 jsp_framework 实现¶

基本流程¶

PuLP 建模¶

计算实例¶

基于 `jsp_framework` 实现¶