A. Lucas, Front. Phys. (2014) 掲載例題の実装と解説 ー 整数長ジョブスケジューリング問題¶
本サンプルコードでは、論文 A. Lucas, "Ising formulations of many NP problems", Front. Phys. (2014) で紹介されている『整数長ジョブスケジューリング問題』に Fixstars Amplify を用いて取り組みます。同論文に紹介されている他の NP 完全・NP 困難な問題も以下で解説しています（カッコ内は論文内で問題に対応する節番号）。
- グラフの分割問題（2.2節）
- 最大クリーク問題（2.3節）
- 厳密被覆問題（4.1節）
- 集合パッキング問題（4.2節）
- 最小頂点被覆問題（4.3節）
- 充足可能性問題（SAT）（4.4節）
- 最小極大マッチング問題（4.5節）
- グラフ彩色問題（6.1節）
- クリーク被覆問題（6.2節）
- 整数長ジョブスケジューリング問題（6.3節）
- ハミルトン閉路問題（7.1節）
- 有向帰還頂点集合問題（8.3節）
- 最小帰還辺集合問題（8.5節）
- グラフ同型性判定問題（9節）
ジョブ割り当て問題¶
$N$ 個のジョブがあり、それぞれのジョブにかかる時間が分かっているとします。それらのジョブを実行できるマシンが $M$ 台あり、$N$ 個のジョブそれぞれをいずれかのマシンに割り当てます。すべてのジョブが完了するまでの時間が最も早くなる割り当て方を求めることを考えます。
ただし、それぞれのマシンは割り当てられたジョブを直列に実行します。つまり、1 つのマシンが複数のジョブを同時に行うことはできません。また、それぞれのジョブにかかる時間は整数であるものとしておきます。
たとえば、それぞれ 1 時間ずつかかるジョブが 3 つあり、マシンが 2 台のとき、2 つのジョブを片方のマシンに、 1 つのジョブをもう片方のマシンに割り当てると、すべてのジョブが完了するまで 2 時間かかります。また、2 時間未満ですべてのジョブを完了することはできないので、これが最適解となります。
ここでは、Fixstars Amplify を用いてこのジョブ割り当て問題を解くプログラムを作成します。定式化は A. Lucas, Front. Phys. (2014) の 6.3 節のものに沿って行います。
問題の作成¶
まず、例題として使用する問題を作成しておきます。ジョブの数とマシンの数、およびそれぞれのジョブにかかる時間を決定します。
import numpy as np
# マシンの数
M = 3
# ジョブの数
N = 7
# 各ジョブにかかる時間
job_lengths = np.array([7, 5, 3, 2, 2, 2, 2])
assert N == len(job_lengths)
定式化¶
以下、$i$ 番目のジョブにかかる時間を $L_i$ とします。
方針¶
$N\times M$ のバイナリ変数テーブル $q$ を用意し、各ジョブをどのマシンで実行するかを表すことにします。 $i$ 番目のジョブを マシン $j$ で行うとき、$q$ の $i$ 行 $j$ 列が $1$ となるようにします。
たとえば、以下のような割り当て方に対応する $q$ は下の表のようになります。
| ジョブ | マシン |
|---|---|
| ジョブ 0 | マシン 0 |
| ジョブ 1 | マシン 2 |
| ジョブ 2 | マシン 2 |
| ジョブ 3 | マシン 1 |
| ジョブ 4 | マシン 1 |
| ジョブ 5 | マシン 1 |
| ジョブ 6 | マシン 1 |
| $q$ | マシン 0 | マシン 1 | マシン 2 |
|---|---|---|---|
| ジョブ 0 | 1 | 0 | 0 |
| ジョブ 1 | 0 | 0 | 1 |
| ジョブ 2 | 0 | 0 | 1 |
| ジョブ 3 | 0 | 1 | 0 |
| ジョブ 4 | 0 | 1 | 0 |
| ジョブ 5 | 0 | 1 | 0 |
| ジョブ 6 | 0 | 1 | 0 |
また、マシンの実行時間の最大値が何であるかを分かりやすくするために、 マシン $0$ の実行時間が最も長くなるようにジョブを割り当てることにします。
目的関数¶
マシン $0$ の実行時間が他のマシンの実行時間より長くなるようにジョブを割り当てるので、 すべてのジョブを完了するまでにかかる時間はマシン $0$ の実行時間と等しくなります。 したがって、目的関数はマシン $0$ の実行時間、つまりマシン $0$ に割り当てられたジョブにかかる時間の総和とすればよいです。 これはジョブにかかる時間 $L$ を用いて
$$ \sum_{i = 0}^{N - 1} L_i q_{i, 0} $$
と書けます。
制約条件¶
$q$ は以下をみたしている必要があります。
- 条件 1 ： 各ジョブはちょうど 1 つのマシンに割り当てられる。つまり、$q$ の各行には $1$ が $1$ つだけある。
- 条件 2 ： それぞれのマシンについて、そのマシンの実行時間はマシン $0$ の実行時間よりも短い。
条件 1 は、
$$ \sum_{j = 0}^{M-1} q_{i, j} = 1 \quad \text{for} \quad i \in \{0, 1, \ldots, N-1\} $$
で表せます。
また、マシン $j$ の実行時間は $\sum_{i = 0}^{N - 1} L_i q_{i, j}$ で表せるので、条件 2 は
$$ \sum_{i = 0}^{N - 1} L_i q_{i, j} \leq \sum_{i = 0}^{N - 1} L_i q_{i, 0} \quad \text{for} \quad j \in \{1, 2, \ldots, M - 1\} $$
と表現できます。
逆に、条件 1 と条件 2 がみたされているとき、$q$ はジョブの割り当て方を表し、かつ目的関数はジョブが完了するまでの時間と等しくなります。
実装¶
上で作成した問題と定式化を使って、実際に問題を解いてみましょう。
まず、Fixstars Amplify SDK の BinarySymbolGenerator を使って $N\times M$ 個のバイナリ変数 $q$ を作成します。
from amplify import VariableGenerator
gen = VariableGenerator()
q = gen.array("Binary", shape=(N, M))
次に、各マシンの総実行時間を、上記で作成した変数配列 q を用いて表しておきます。マシン $j$ の総実行時間は、$\displaystyle \sum_{i = 0}^{N
- 1}
L_i q_{i, j}$ で表されます。この式の $L$ は各ジョブにかかる時間を表す配列であり、コード上では job_lengths という名前の numpy 配列です。
from amplify import PolyArray, einsum
execution_times: PolyArray = einsum("i,ij->j", job_lengths, q) # type:ignore
execution_times
次に、目的関数を作成します。先の説明の通り、目的関数はマシン $0$ の総実行時間です。
cost = execution_times[0]
条件 1 に対応する制約条件を作成します。条件 1 は「それぞれのジョブはちょうど $1$ つのマシンに割り当てられる」ということを意味し、 $q$ の各行にひとつだけ $1$
があるという制約条件です。one_hot 関数の axis パラメータに 1 を指定することで、二次元配列の各行に対する one-hot
制約を一度に生成できます。
from amplify import one_hot
constraint1 = one_hot(q, axis=1)
条件 2 に対応する制約条件を作成します。条件 2 はマシン $0$
の実行時間は他のマシンの実行時間以上であるという条件でした。配列のそれぞれの要素が取る値に対して制約条件を課したい場合、less_equal などの関数の
axis パラメータに空タプルを指定します。
from amplify import less_equal
constraint2 = less_equal(execution_times[1:] - execution_times[0], 0, axis=())
作成した目的関数と制約条件をまとめて、組合せ最適化モデルを構築します。
model = cost + constraint1 + constraint2
また、今回使用している入力変数の数は $N \times M = 21$ ですが、作成した model には不等式制約が含まれるため、論理模型に変換される際に補助変数が発行されます。その結果、論理変数の数は
$q$ に含まれるバイナリ決定変数の数よりも多くなります。
クライアントを設定し、Fixstars Amplify Annealing Engine (AE) で実行します。
from amplify import FixstarsClient, solve
from datetime import timedelta
client = FixstarsClient()
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # ローカル環境等で使用する場合は、Fixstars Amplify AE のアクセストークンを入力してください。
client.parameters.timeout = timedelta(milliseconds=1000) # タイムアウトは 1000 ms
# ソルバーを定義して実行
result = solve(model, client)
解が見つかったかどうかを確認します。Amplify SDK は制約条件をみたす解を自動でフィルターするので、result
が空でなければ、制約条件をみたす解が見つかったと分かります。
if len(result) == 0:
print("解が見つかりませんでした。")
else:
print("解が見つかりました。")
すべてのジョブが完了するまでの時間は目的関数の値と等しいので、以下のようにして確認できます。
result.best.objective
最後に、結果を可視化します。
import matplotlib.pyplot as plt
values = q.evaluate(result.best.values)
assigned_machines = np.where(values == 1)[1]
# x軸を描画
plt.xticks(range(M), [f"machine {i}" for i in range(M)])
# 描画
bottom = np.zeros(M) # 現在の棒グラフの上端
for i, j in enumerate(assigned_machines):
bar = plt.bar(j, job_lengths[i], bottom=bottom[j])
plt.bar_label(bar, labels=[f"job {i}"], label_type="center")
bottom[j] += job_lengths[i]
