作业车间调度问题求解框架:OR-Tools 约束求解器

发布于:2021-08-22 | 分类:optimization , mathematics


本文根据作业车间调度问题的数学描述,利用 Google 的整数规划求解器 OR-Tools 进行求解,作为 jsp_framework 内置的第一个求解器。

OR-Tools

OR-Tools 是 Google 基于 C++ 开发的针对组合优化问题的求解库,同时支持 Python,C#,或 Java 语言的调用。

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

pip install ortools

数学模型

前文给出了作业车间调度问题的两种典型描述,其中数学模型描述为:

\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 的步骤为:用内置语言建立模型例如定义变量、约束及目标函数,然后使用内置或者第三方求解器求解。

注意

OR-Tools 使用 CP-SAT 约束求解器求解作业车间调度问题。CP-SAT 求解器仅支持整数变量。

上一节数学模型中将工序完工时刻作为求解变量,为了建模方便,这里为每个工序创建两个变量 s_{ij}c_{ij},分别表示开始时刻和结束时刻。OR-Tools 内置了 NewIntVar 整数变量,根据函数签名可知,我们需要指定变量名及变量范围 [lb, ub]

def NewIntVar(self, lb, ub, name): pass

基于这些变量,

  • AddMaxEqualityc_{ij} 的最大值然后最小化得到目标函数
  • 定义变量下限 lb \geq 0 即可满足第一个非负的约束条件
  • 第二个约束条件基于前后工序的 s_{ij}c_{ij} 实现

为了实施第三个约束条件,即同一机器上的任意两个工序不重叠,需要用到 OR-Tools区间变量 NewIntervalVar。针对不可重叠的所有工序/区间,直接施加内置的 NoOverlap 约束。

def NewIntervalVar(self, start, size, end, name): pass

其中,startendsize 分别表示区间的开始、结束位置及区间的大小,并且始终满足 start + size = end。对于本问题,startend 为上一步的 s_{ij}c_{ij}size 为工序加工时间常量。

基于 jsp_framework 实现

OR-Tools 的官方文档提供了一个求解作业车间调度问题的 完整案例,本文将基于本系列的 jsp_framework 实现。参考 Python建模 部分,基本流程为:

自定义一个继承自 JSSolverGoogleORCPSolver 类,然后重点实现 do_solve() 方法。完整代码参考ortools.py。其中,

  • 为了避免求解规模太大而失去响应,人为设定运行时间的上限:

    solver.parameters.max_time_in_seconds = 300.0
  • 为了实时输出当前最优解,定义如下执行回调函数的类:

    class VarArraySolutionPrinter(cp_model.CpSolverSolutionCallback):
        '''Output intermediate solutions.'''
        def __init__(self, variables:dict, solver:JSSolver, solution:JSSolution):
            '''Initialize with variable map: operation step -> OR-Tools variable.'''
            cp_model.CpSolverSolutionCallback.__init__(self)
            self.__variables = variables
            self.__solver = solver
            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.__solver.update_solution(self.__solution)

    相应地,将普通的求解方式:

    status = solver.Solve(model)

    改为:

    solution_printer = VarArraySolutionPrinter(variables, problem, solution)
    status = solver.SolveWithSolutionCallback(model, solution_printer)

计算实例

最后,求解几个标准问题。求解时间上限设定为300秒。

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

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

# solver
solvers = [GoogleORCPSolver(max_time=300)]

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

结果如下表所示:

  • 前5行变量总数300以内,OR-Tools 基本得到了最优解。例4虽然没能在规定时间内求得最优解,但是相对误差0.5%已经很小。

  • 后5行变量总数在[450, 1000]量级,除了例9有些反常外,其余都没能在300秒内求得最优解,并且精度随着问题复杂度(例如变量总数、工序数等)增加而降低。

+----+---------+------------------+-------+--------------+----------+---------+-------+
| ID | Problem |      Solver      | Scale |   Optimum    | Solution | Error % |  Time |
+----+---------+------------------+-------+--------------+----------+---------+-------+
| 1  |   ft06  | GoogleORCPSolver |  6x6  |      55      |    55    |   0.0   |  0.1  |
| 2  |   la01  | GoogleORCPSolver |  10x5 |     666      |   666    |   0.0   |  0.2  |
| 3  |   ft10  | GoogleORCPSolver | 10x10 |     930      |   930    |   0.0   |  3.4  |
| 4  |  swv01  | GoogleORCPSolver | 20x10 |     1407     |   1414   |   0.5   | 300.1 |
| 5  |   la38  | GoogleORCPSolver | 15x15 |     1196     |   1196   |   0.0   | 168.7 |
| 6  |   ta31  | GoogleORCPSolver | 30x15 |     1764     |   1814   |   2.8   | 300.2 |
| 7  |  swv12  | GoogleORCPSolver | 50x10 | (2972, 3003) |   3339   |   11.8  | 300.3 |
| 8  |   ta42  | GoogleORCPSolver | 30x20 | (1867, 1956) |   2130   |   11.4  | 303.5 |
| 9  |   ta54  | GoogleORCPSolver | 50x15 |     2839     |   2862   |   0.8   | 468.8 |
| 10 |   ta70  | GoogleORCPSolver | 50x20 |     2995     |   3283   |   9.6   | 600.3 |
+----+---------+------------------+-------+--------------+----------+---------+-------+

问题规模较小时,能高效求出最优解;但随着问题规模的增大,求解效率急剧下降。例如,上例中即便将计算时间上限增加到600秒,后5个例子的结果也没有显著改善。这是整数规划等精确求解方法应用于实际大规模调度问题的限制所在。