jsp framework: pulp solver

docs/2021-08-08-作业车间调度问题求解框架.md

@@ -50,6 +50,8 @@ https://github.com/dothinking/jsp_framework
 
 - [OR-Tools 约束求解器](2021-08-22-作业车间调度问题求解框架：OR-Tools约束求解器.md)
 
+- [PuLP求解器](2021-08-29-作业车间调度问题求解框架：PuLP求解器.md)
+
 - [基于规则指派](2021-08-28-作业车间调度问题求解框架：基于规则指派求解器.md)
 
 

docs/2021-08-29-作业车间调度问题求解框架：PuLP求解器.md

@@ -0,0 +1,304 @@
+---
+categories: [optimization, mathematics]
+tags: [job shop schedule]
+---
+
+# 作业车间调度问题求解框架：PuLP 求解器
+
+---
+
+
+本文根据作业车间调度问题的数学描述，利用 Python 整数规划包 PuLP 进行建模和求解，并集成到 `jsp_framework` 框架中。
+
+## PuLP
+
+PuLP 是开源基金会 [COIN-OR](https://github.com/coin-or) (Computational Infrastructure for Operations Research) 维护的一个线性规划建模工具（LP modeler），内置了 CBC（COIN-OR Branch-and-Cut）求解器，也支持集成各种开源（如GLPK） 或者商业（CPLEX、GUROBI等）的线性规划求解器。
+
+> https://github.com/coin-or/pulp
+
+作为 Python 库，PuLP 支持 pip 一键安装。
+
+```python
+pip install pulp
+```
+
+更多入门材料参考：
+
+- [官方文档](https://coin-or.github.io/pulp/index.html)
+- [PuLP简介](http://fancyerii.github.io/2020/04/18/pulp/)
+
+## 数学模型
+
+[前文](2021-08-08-作业车间调度问题求解框架：问题描述.md) 给出了作业车间调度问题的数学模型。为了方便接下来的建模描述，改为以工序 **开始时间** 为决策变量的描述方式：
+
+\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 类
+    ```python
+    LpProblem(name='NoName', sense=LpMinimize)
+    ```
+
+    其中，`sense` 参数目标函数是求极大值（LpMaximize）还是极小值（LpMinimize）。
+
+    建立约束后调用 `solve(solver=None, **kwargs)` 求解，默认为 CBC 求解器。
+
+
+- LpVariable 类
+
+    ```python
+    LpVariable(name, lowBound=None, upBound=None, cat='Continuous', e=None)
+    ```
+
+    其中，
+
+    - `name` 表示变量名
+    - `lowBound` 和 `upBound` 表示变量范围下界和上界，默认负无穷到正无穷
+    - `cat` 指定变量是离散类型（Integer或Binary），还是连续类型（Continuous）
+
+    特别地，当变量个数特别多且满足批量化处理时，可以使用如下方法建立字典格式的变量集合。
+
+    ```python
+    LpVariable.dicts(name, indexs, lowBound=None, upBound=None, cat='Continuous', indexStart=[])
+    ```
+
+    其中，
+
+    - `name` 指定所有变量的前缀
+    - `indexs` 是一个列表，将会为其中每一个元素创建一个变量，并建立对应关系
+
+
+## 基于 `jsp_framework` 实现
+
+本文参考 [案例](https://github.com/KevinLu43/JSP-by-using-Mathematical-Programming-in-Python/)，基于 `jsp_framework` 进行实现和集成，完整代码参考：
+
+> https://github.com/dothinking/jsp_framework/blob/dev/jsp_fwk/solver/pulp.py
+
+
+### 基本流程
+
+参考 [Python建模](2021-08-14-作业车间调度问题求解框架：Python建模.md)，自定义一个继承自 JSSolver 的 PuLPSolver 类，然后重点实现 `do_solve()` 方法。
+
+```python
+class PuLPSolver(JSSolver):
+
+    ...
+
+    def do_solve(self, problem:JSProblem):
+        '''Solve JSP with PuLP default CBC solver.'''
+        # Initialize an empty solution from problem
+        solution = JSSolution(problem)
+
+        # 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)        
+        problem.update_solution(solution) # update solution
+```
+
+!!! warning "注意"
+
+    当前版本采用了默认的 CBC （COIN-OR Branch-and-Cut）求解器。
+
+
+### PuLP 建模
+
+即上一代码片段 `self.__create_model(solution)` 的具体展开。
+
+- 创建问题
+
+    ```python
+    # create the model
+    model = pulp.LpProblem("min_makespan", pulp.LpMinimize)
+    ```
+
+- 创建变量
+
+    - 基本变量：每一道工序的开始时间
+
+        ```python
+        # (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变量：表征同一机器上任意两道工序的前后顺序
+
+        ```python
+        # (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`得到了机器队列工序的两两组合。
+
+- 创建目标函数：总的加工周期
+
+    ```python
+    # objective: makespan / flowtime
+    s_max = pulp.LpVariable(name='max_start_time', \
+                            lowBound=0, \
+                            upBound=max_time, \
+                            cat='Integer')
+    model += s_max
+    ```
+
+    注意：此处目标函数仅仅是一个普通的变量，后续添加约束保证其为最大的加工周期
+
+
+- 添加约束
+
+    - 最大加工周期约束
+
+    ```python
+    # (1) the max start time
+    for _, ops in solution.job_ops.items():
+        last_op = ops[-1]
+        model += (s_max-variables[last_op]) >= last_op.source.duration
+    ```
+
+    - 作业视角的工序顺序约束
+
+    ```python
+    # (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
+    ```
+
+    - 机器视角的工序顺序约束
+
+    ```python
+    # (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秒。
+
+```python
+# benchmark.py
+from jsp_fwk import (JSProblem, BenchMark)
+from jsp_fwk.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(max_time=300)]
+
+# 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  |  pulp  |     6 x 6     |      55      |   55.0   |   0.0   |  3.5  |
+| 2  |   la01  |  pulp  |     10 x 5    |     666      |  666.0   |   0.0   | 209.0 |
+| 3  |   ft10  |  pulp  |    10 x 10    |     930      |  1034.0  |   11.2  | 301.2 |
+| 4  |  swv01  |  pulp  |    20 x 10    |     1407     |  2941.0  |  109.0  | 301.2 |
+| 5  |   la38  |  pulp  |    15 x 15    |     1196     |  1769.0  |   47.9  | 301.2 |
+| 6  |   ta31  |  pulp  |    30 x 15    |     1764     | 13294.0  |  653.6  | 301.1 |
+| 7  |  swv12  |  pulp  |    50 x 10    | (2972, 3003) |  9527.0  |  218.9  | 309.9 |
+| 8  |   ta42  |  pulp  |    30 x 20    | (1867, 1956) | 23202.0  |  1113.8 | 507.1 |
+| 9  |   ta54  |  pulp  |    50 x 15    |     2839     | 26229.0  |  823.9  | 608.0 |
+| 10 |   ta70  |  pulp  |    50 x 20    |     2995     | 36851.0  |  1130.4 | 625.6 |
++----+---------+--------+---------------+--------------+----------+---------+-------+
+```
+
+计算效果并不理想，只有前两个问题即工序总数 50 以内得到了最优解，其余案例误差较大。由此表明，**PuLP 默认的 CBC 求解器在作业车间调度问题上的性能不如 Google OR-Tools 的约束求解器 CpSolver**。
+
+PuLP 支持其他求解器例如商业求解器 CPLEX 和 GUROBI，有机会可以对比一下结果。