生産計画の最適化(OR-Tools)

#はじめに

- Googleの数理最適化ツール(OR-Tools)を使用して、生産計画を最適化するプログラムを作成してみた。
- https://developers.google.com/optimization

【追記】
今回の生産計画の仕様では、それぞれの設備は同一で区別できない。
あらかじめ任意の一つのジョブを割付ける設備を固定することで
解が得られるまでの時間が短くなる。(タイブレーク条件)

#生産計画の仕様

- 複数の作業ロット(lot)があり、各ロットは連続するジョブ(job)で構成される。
- 各ジョブには、生産に必要な時間(size)が指定される。
- 各ジョブは、複数の設備のいずれかに割付ける必要がある。
- 同一の設備に割付けられたジョブは、互いに重なってはいけない。
- 同一のロットを構成するジョブは順番に生産する必要がある。
- 全てのジョブが終了するまでの時間を最小化する。

#変数の定義方法

- ジョブ毎に、開始時刻(start_var)と終了時刻(end_var)を表す変数を定義する。
- ジョブ・設備毎に割付けるか否かを表すブール変数(bool_var)を定義して、インターバル変数(interval_var)を定義する。
- ジョブ毎にブール変数の合計が1となる制約を定義する。
- 設備毎にインターバル変数が互いに重ならない制約を定義する。
- 同じロットを構成するジョブに前後の順序制約を定義する。

![RCSP.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/192430/7f928e48-071d-5d2e-08af-442b2c26f551.png)

#データの生成

```python
from ortools.sat.python import cp_model
from collections import defaultdict
from dataclasses import dataclass
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import random

num_machines = 3
num_lots = 5
num_jobs_per_lot = 3
num_jobs = num_lots * num_jobs_per_lot

@dataclass
class job_type:
    id: int
    lot_id: int
    size: int

random.seed(1)
jobs_list = [job_type(j, j//num_jobs_per_lot, random.randint(1, 10)) for j in range(num_jobs)]
horizon = sum([job.size for job in jobs_list])
```

#変数と制約条件の定義

```python
model = cp_model.CpModel()
machine_to_intervals = defaultdict(list)
for job in jobs_list:
    start_var = model.NewIntVar(0, horizon, 'start_' + str(job.id))
    end_var = model.NewIntVar(0, horizon, 'end_' + str(job.id))
    job.start_var = start_var
    job.end_var = end_var
    bool_var_list = []
    for m in range(num_machines):
        suffix = str(m) + '_' + str(job.id)
        bool_var = model.NewBoolVar('bool_' + suffix)
        bool_var_list.append(bool_var)
        interval_var = model.NewOptionalIntervalVar(start_var, job.size, end_var, bool_var, 'interval_' + suffix)
        interval_var.job = job
        interval_var.bool_var = bool_var
        machine_to_intervals[m].append(interval_var)
    model.Add(sum(bool_var_list) == 1)

for m in machine_to_intervals:
    model.AddNoOverlap(machine_to_intervals[m])

for j in range(num_jobs-1):
    if jobs_list[j].lot_id != jobs_list[j+1].lot_id: continue
    model.Add(jobs_list[j].end_var <= jobs_list[j+1].start_var)

model.Add(machine_to_intervals[0][0].bool_var == 1) # tie-break
```

#最適化

```python
span_var = model.NewIntVar(0, horizon, 'span_var')
model.AddMaxEquality(span_var, [j.end_var for j in jobs_list])
model.Minimize(span_var)

solver = cp_model.CpSolver()
solver.Solve(model)
print(solver.StatusName(), solver.ObjectiveValue())
```

手元の環境では、およそ230msで最適値28が得られた。

#結果出力

```python
for m in machine_to_intervals:
    for i in machine_to_intervals[m]:
        if solver.Value(i.bool_var) == 0: continue
        i.job.start = solver.Value(i.job.start_var)
        i.job.machine_id = m

result = [[job.machine_id, job.lot_id, job.id, job.start, job.start+job.size, job.size] for job in jobs_list]
result = pd.DataFrame(result, columns=['machine', 'lot', 'job', 'start', 'end', 'size'])
result.sort_values(['machine', 'start'], inplace=True, ignore_index=True)
print(result)
```

||machine|lot|job|start|end|size|
|---:|---:|---:|---:|---:|---:|---:|
|0|0|4|12|0|1|1|
|1|0|2|6|1|9|8|
|2|0|0|1|9|19|10|
|3|0|1|5|19|27|8|
|4|1|3|9|0|4|4|
|5|1|3|10|4|6|2|
|6|1|0|0|6|9|3|
|7|1|1|4|9|11|2|
|8|1|2|7|11|19|8|
|9|1|0|2|19|21|2|
|10|1|4|14|21|28|7|
|11|2|1|3|0|5|5|
|12|2|3|11|6|14|8|
|13|2|4|13|14|21|7|
|14|2|2|8|21|28|7|

#ガントチャートの表示

```python
span_max = max([job.start+job.size for job in jobs_list])
cmap = plt.cm.get_cmap('hsv', num_lots+1)
fig = plt.figure(figsize=(12, 8))
for m in range(num_machines):
    ax = fig.add_subplot(num_machines, 1, m+1, yticks=[], ylabel=m)
    ax.set_xlim(-1, span_max+1)
    ax.set_ylim(-0.1, 1.1)
    for job in jobs_list:
        if job.machine_id != m: continue
        rectangle = matplotlib.patches.Rectangle((job.start, 0), job.size, 1, fill=False, color=cmap(job.lot_id), hatch='/')
        ax.add_patch(rectangle)
        rx, ry = rectangle.get_xy()
        cx = rx + rectangle.get_width()/2
        cy = ry + rectangle.get_height()/2
        lab = str(job.lot_id) + '-' + str(job.id)
        ax.annotate(lab, (cx, cy), ha='center', va='center')
plt.show()
```

![gant.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/192430/bffb4036-9d8c-ac94-25d9-2905ec93046e.png)

ロット毎に色を変えている。ジョブの処理順序を満たしながら、きっちり詰まっている。

#おわりに

- 問題の規模が小さければ、OR-Toolsを使用して生産計画の最適化問題を解くことが可能。
- 但し、現実の問題では変数の数が数万となることもある。見極めが必要。