jsp framework: ortools solver

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

@@ -18,6 +18,8 @@ tags: [job shop schedule]
 
 https://github.com/dothinking/jsp_framework
 
+
+
 ## 建模
 
 - [问题描述](2021-08-08-作业车间调度问题求解框架：问题描述.md)
@@ -28,3 +30,7 @@ https://github.com/dothinking/jsp_framework
 
 
 ## 求解
+
+- [OR-Tools 约束求解器](2021-08-22-作业车间调度问题求解框架：OR-Tools约束求解器.md)
+
+- [基于规则指派](2021-08-22-作业车间调度问题求解框架：基于规则指派求解器.md)

docs/2021-08-22-作业车间调度问题求解框架：OR-Tools约束求解器.md

@@ -0,0 +1,193 @@
+---
+categories: [optimization, mathematics]
+tags: [job shop schedule]
+---
+
+# 作业车间调度问题求解框架：OR-Tools 约束求解器
+
+---
+
+本文根据作业车间调度问题的数学描述，利用 Google 的整数规划求解器 `OR-Tools` 进行求解，作为 `jsp_framework` 内置的第一个求解器。
+
+## OR-Tools
+
+`OR-Tools` 是 Google 基于 C++ 开发的针对组合优化问题的求解库，同时支持 Python，C#，或 Java 语言的调用。
+
+> https://developers.google.cn/optimization/
+
+作为Python库，支持`pip`一键安装。
+
+```python
+pip install ortools
+```
+
+## 数学模型
+
+[前文](2021-08-08-作业车间调度问题求解框架：问题描述.md)给出了作业车间调度问题的两种典型描述，其中数学模型描述为：
+
+\begin{align\*}
+&\min\,\, c_{m} = \max \left \\{ c_{ij} \, | \, (i,j) \in \mathbb{O} \right \\} \\\\
+\\\\
+&s.t.\quad
+\begin{cases}
+c_{ij} \geq p_{ij} & (i,j) \in \mathbb{O} \\\\
+\\\\
+c_{ij}-c_{kj} \geq p_{ij} & if \,\, (k,j) \rightarrow (i,j), \,\, i, k \in \mathbb{M}^j \\\\
+\\\\
+c_{ij} - c_{ik} \geq p_{ij} \,\, or \,\, c_{ik}-c_{ij}  \geq p_{ik} & (i,j), (i,k) \in \mathbb{O}, \,\, j \neq k
+\end{cases}
+\end{align\*}
+
+具体变量说明参考原文，其中两个基本变量：
+
+- $c_{ij}$ 表示工序 $(i,j)$ 即作业 $j$ 的某个在机器 $i$ 上进行的工序的完成时间
+- $p_{ij}$ 表示作业 $j$ 的某个在机器 $i$ 上进行的工序所需的加工时间
+
+三个约束中：
+
+- 第一个表明工序的开始时间非负
+- 第二个是工序序列约束，即同一个作业的相邻工序必须满足先后顺序
+- 第三个是资源约束，即每一台机器在某个时刻只能加工一道工序
+
+
+## `OR-Tools` 模型
+
+使用 `OR-Tools` 的步骤为：用内置语言建立模型例如定义变量、约束及目标函数，然后使用内置或者第三方求解器求解。
+
+!!! warning "注意"
+
+    `OR-Tools` 使用 `CP-SAT` 约束求解器求解作业车间调度问题。`CP-SAT` 求解器仅支持整数变量。
+
+上一节数学模型中将工序完工时刻作为求解变量，为了建模方便，这里为每个工序创建两个变量 $s_{ij}$ 和 $c_{ij}$，分别表示开始时刻和结束时刻。`OR-Tools` 内置了 `NewIntVar` 整数变量，根据函数签名可知，我们需要指定变量名及变量范围 $[lb, ub]$。
+
+```python
+def NewIntVar(self, lb, ub, name): pass
+```
+
+基于这些变量，
+
+- 用`AddMaxEquality`取 $c_{ij}$ 的最大值然后最小化得到目标函数
+- 定义变量下限 $lb \geq 0$ 即可满足第一个非负的约束条件
+- 第二个约束条件基于前后工序的 $s_{ij}$ 和 $c_{ij}$ 实现
+
+为了实施第三个约束条件，即同一机器上的任意两个工序不重叠，需要用到 `OR-Tools` 的 **区间变量** `NewIntervalVar`。针对不可重叠的所有工序/区间，直接施加内置的 `NoOverlap` 约束。
+
+```python
+def NewIntervalVar(self, start, size, end, name): pass
+```
+
+其中，`start`、`end` 及 `size` 分别表示区间的开始、结束位置及区间的大小，并且始终满足 $start + size = end$。对于本问题，`start`、`end` 为上一步的 $s_{ij}$ 和 $c_{ij}$，`size` 为工序加工时间常量。
+
+
+## 基于 `jsp_framework` 实现
+
+`OR-Tools` 的官方文档提供了一个求解作业车间调度问题的 [完整案例](https://developers.google.cn/optimization/scheduling/job_shop)，本文将基于本系列的 `jsp_framework` 实现。
+
+- 基本流程（参考 [Python建模](2021-08-14-作业车间调度问题求解框架：Python建模.md) 部分）：
+
+    自定义一个继承自 `JSSolver` 的 `GoogleORCPSolver` 类，然后重点实现 `do_solve()` 方法。
+
+
+- 完整代码参考：
+
+    > https://github.com/dothinking/jsp_framework/blob/dev/jsp_fwk/solver/ortools.py
+
+
+其中，
+
+- 为了避免求解规模太大而失去响应，人为设定运行时间的上限：
+
+    ```python
+    solver.parameters.max_time_in_seconds = 300.0
+    ```
+
+
+- 为了实时输出当前最优解，定义如下执行回调函数的类：
+
+    ```python
+
+    class VarArraySolutionPrinter(cp_model.CpSolverSolutionCallback):
+        '''Output intermediate solutions.'''
+        def __init__(self, variables:dict, problem:JSProblem, solution:JSSolution):
+            '''Initialize with variable map: operation step -> OR-Tools variable.'''
+            cp_model.CpSolverSolutionCallback.__init__(self)
+            self.__variables = variables
+            self.__problem = problem
+            self.__solution = solution
+
+        def on_solution_callback(self):
+            '''Pass data back to domain class.'''
+            # assign OR-Tools solution back to JSPSolution
+            for op, var in self.__variables.items():
+                op.update_start_time(self.Value(var.start))
+
+            # update solution
+            self.__problem.update_solution(self.__solution)
+    ```
+
+    相应地，将普通的求解方式：
+
+    ```python
+    status = solver.Solve(model)
+    ```
+
+    改为：
+
+    ```python
+    solution_printer = VarArraySolutionPrinter(variables, problem, solution)
+    status = solver.SolveWithSolutionCallback(model, solution_printer)
+    ```
+
+## 计算实例
+
+最后，求解几个标准问题。求解时间上限设定为300秒。
+
+```python
+from jsp_fwk import (JSProblem, JSSolution)
+from jsp_fwk.solver import GoogleORCPSolver
+
+# benchmarks
+names = ['ft06', 'la01', 'ft10', 'swv01', 'la38', \
+        'ta31', 'swv12', 'ta42', 'ta54', 'ta70']
+
+res = []
+res.append('name,num_job,num_machine,optimum,solution,time')
+for name in names:
+    p = JSProblem(benchmark=name)
+    s = GoogleORCPSolver(max_time=300) # set upper limit of solving time
+    s.solve(problem=p, interval=None)
+
+    # wait for termination
+    while True:
+        if not s.is_running: break
+
+    # collect results
+    items = [name, len(p.jobs), len(p.machines), p.optimum, \
+                p.solution.makespan, s.user_time]
+    res.append(','.join(map(str, items)))
+
+print(res)
+```
+
+结果如下表所示：
+
+- 前5行变量总数300以内，`OR-Tools` 基本得到了最优解。例4虽然没能在规定时间内求得最优解，但是相对误差0.5%已经很小。
+
+- 后5行变量总数在$[450, 1000]$量级，都没能在300秒内求得最优解，并且精度随着问题复杂度（例如变量总数、工序数等）增加而急剧降低。
+
+| |name|num of <br> job|num of <br> machine|optimum|solution|err / %|time|
+|---|---|---|---|---|---|---|---|
+|1|ft06|6|6|55|55|0.0|0.8
+|2|la01|10|5|666|666|0.0|6.8
+|3|ft10|10|10|930|930|0.0|20.6
+|4|swv01|20|10|1407|1414|0.5|300.1
+|5|la38|15|15|1196|1196|0|193.5
+|6|ta31|30|15|1764|1891|7.2|300.1
+|7|swv12|50|10|(2972, 3003)|3339|11.8|300.1
+|8|ta42|30|20|(1867, 1956)|2546|33.2|300.1
+|9|ta54|50|15|2839|3192|12.4|299.7
+|10|ta70|50|20|2995|3948|31.8|299.0
+
+
+问题规模较小时，高效求出最优解；但随着问题规模的增大，求解效率急剧下降。这是整数规划等精确求解方法应用于实际大规模调度问题的限制所在。
+