ジョブショップスケジューリング

ジョブショップスケジューリング問題では、いくつかのジョブ（仕事）をいくつかの機械で処理します。各ジョブは、決められた順序ですべての機械を 1 回ずつ回り、機械ごとに決められた時間の作業をします。 1 台の機械は一度に 1 つの作業しかできず、始めた作業は途中で止めません。目的は、メイクスパン（すべてのジョブが終わる時刻）を最小にすることです。

このページでは、Fisher と Thompson が 1963 年に示した 6 ジョブ × 6 機械の問題 ft06 を解きます。小さな問題ですが、スケジューリングの手法を比べるベンチマークとして広く使われていて、最小のメイクスパンは 55 であることが知られています。

各ジョブが機械を回る順序と処理時間は次のとおりです。「M3: 1」は、機械 3 で 1 の時間の作業をすることを表します。

ジョブ	1 番目	2 番目	3 番目	4 番目	5 番目	6 番目
J1	M3: 1	M1: 3	M2: 6	M4: 7	M6: 3	M5: 6
J2	M2: 8	M3: 5	M5: 10	M6: 10	M1: 10	M4: 4
J3	M3: 5	M4: 4	M6: 8	M1: 9	M2: 1	M5: 7
J4	M2: 5	M1: 5	M3: 5	M4: 3	M5: 8	M6: 9
J5	M3: 9	M2: 3	M5: 5	M6: 4	M1: 3	M4: 1
J6	M2: 3	M4: 3	M6: 9	M1: 10	M5: 4	M3: 1

定式化

ジョブ $j$ が機械 $m$ で作業を始める時刻を整数変数 $s_{j,m}$ で表します。処理時間を $p_{j,m}$ とすると、その作業は $s_{j,m} + p_{j,m}$ に終わります。

時刻の上限 $T$ には、すべての処理時間の和（この問題では 197）を使います。すべての作業を 1 つずつ順に行っても $T$ までに終わるので、最適なスケジュールは必ずこの範囲にあります。そこで、$s_{j,m}$ の範囲を $0 \le s_{j,m} \le T - p_{j,m}$ とします。

ジョブの順序. ジョブ $j$ が機械 $a$ の次に機械 $b$ を回るとき、機械 $b$ の作業は機械 $a$ の作業が終わってから始めます。

\[s_{j,a} + p_{j,a} \le s_{j,b}\]

機械の重なり. 同じ機械 $m$ を使う 2 つのジョブ $i$ と $j$ は、どちらかが先に終わってから、もう一方が始まります。どちらが先かを表すバイナリ変数 $x_{m,i,j}$（$i$ が先なら 1、$j$ が先なら 0）を使うと、次の 2 つの不等式で書けます。

\[\begin{aligned} s_{i,m} + p_{i,m} - s_{j,m} &\le T\,(1 - x_{m,i,j}) \\ s_{j,m} + p_{j,m} - s_{i,m} &\le T\,x_{m,i,j} \end{aligned}\]

$x_{m,i,j} = 1$ のとき、1 つ目は「$i$ が終わってから $j$ が始まる」という制約になります。 2 つ目の右辺は $T$ になり、$s$ の範囲から左辺は $T$ 以下なので、いつも成り立ちます。 $x_{m,i,j} = 0$ のときは役割が入れ替わります。このように大きな定数を掛けて不要なほうの制約を無効にする書き方は、big-M 法と呼ばれます。

メイクスパン. ジョブ $j$ の最後の作業が終わる時刻を $C_j$ とすると、メイクスパンは $\max_j C_j$ です。 3 つ以上の式の最大値を最小にするには、最大値・最小値の最小化と最大化と同じように、整数変数 $z$ と ReLU を使います。

モデルは次のとおりです。

\[\begin{aligned} \text{最小化}\quad & 10000\,z + \sum_{j,m} s_{j,m} + 20000 \sum_j \mathrm{relu}(C_j - z) \\ \text{制約}\quad & \text{ジョブの順序と機械の重なりの不等式} \end{aligned}\]

重みは次のように決めています。

$\mathrm{relu}(C_j - z)$ の重み 20000 は $z$ の重み 10000 の 2 倍です。$z$ がある $C_j$ より小さいときは、$z$ を 1 増やすとエネルギーが下がるので、最適解では $z$ がちょうど最大値 $\max_j C_j$ になります。
開始時刻の和 $\sum s_{j,m}$ は、すべての作業をできるだけ早く始めさせるための項です。この和は最大でも 6895（36 個の作業の $T - p_{j,m}$ の和）なので、メイクスパンを 1 短くする得（10000）のほうが常に大きくなります。そのため、まずメイクスパンが最小になり、その中で開始時刻の和が最小になります。試した範囲では、この項を加えたほうがソルバーが速く最適値に届きました。
制約は qbpp::cons() で書き、重み 10000000 を付けます。制約を満たすスケジュールの目的関数は $10000 \times 197 + 6895$ 以下で、どの制約を破っても 10000000 以上のエネルギーが加わるので、最適解は必ずすべての制約を満たします。QUBO++ は cons() の不等式を直接扱うので、スラック変数は要りません（非線形関数とネイティブ制約）。

変数は、開始時刻の整数変数 36 個、順序のバイナリ変数 90 個（6 台の機械 × 15 組）、$z$ の 1 個です。

プログラム

#include <algorithm>
#include <iostream>
#include <vector>

#include <qbpp/abs3_solver.hpp>
#include <qbpp/qbpp.hpp>

int main() {
  // ft06: 各ジョブが機械を回る順序と、各機械での処理時間
  std::vector<std::vector<int>> route = {{2, 0, 1, 3, 5, 4}, {1, 2, 4, 5, 0, 3},
                                         {2, 3, 5, 0, 1, 4}, {1, 0, 2, 3, 4, 5},
                                         {2, 1, 4, 5, 0, 3}, {1, 3, 5, 0, 4, 2}};
  auto time = qbpp::array({{3, 6, 1, 7, 6, 3},  {10, 8, 5, 4, 10, 10}, {9, 1, 5, 4, 7, 8},
                           {5, 5, 5, 3, 8, 9},  {3, 3, 9, 1, 5, 4},    {10, 3, 1, 3, 4, 9}});
  int n = route.size();
  int machines = 6;

  // T: すべての処理時間の和（1 つずつ順に処理しても T までに終わる）
  int T = 0;
  for (int j = 0; j < n; ++j)
    for (int m = 0; m < machines; ++m) T += time[j][m];

  // s[j][m]: ジョブ j が機械 m で作業を始める時刻（0 以上 T - 処理時間 以下）
  auto s = 0 <= qbpp::int_var("s", n, machines) <= T - time;

  // 各ジョブは決められた順序で機械を回る
  qbpp::Expr f = 0;
  for (int j = 0; j < n; ++j) {
    for (int k = 0; k + 1 < machines; ++k) {
      int a = route[j][k];
      int b = route[j][k + 1];
      f += 10000000 * qbpp::cons(s[j][a] + time[j][a] - s[j][b] <= 0);
    }
  }

  // 1 台の機械では 2 つの作業が重ならない（x[m][i][j] = 1 なら i が先）
  auto x = qbpp::var("x", machines, n, n);
  for (int m = 0; m < machines; ++m) {
    for (int i = 0; i < n; ++i) {
      for (int j = i + 1; j < n; ++j) {
        f += 10000000 * qbpp::cons(s[i][m] + time[i][m] - s[j][m] - T * (1 - x[m][i][j]) <= 0);
        f += 10000000 * qbpp::cons(s[j][m] + time[j][m] - s[i][m] - T * x[m][i][j] <= 0);
      }
    }
  }

  // メイクスパン z（すべてのジョブの終了時刻の最大値）を最小にし、次に開始時刻の和を小さくする
  auto z = 0 <= qbpp::int_var("z") <= T;
  f += 10000 * z + qbpp::sum(s);
  for (int j = 0; j < n; ++j) {
    int last = route[j][machines - 1];
    f += 20000 * qbpp::relu(s[j][last] + time[j][last] - z);
  }
  f.simplify_as_binary();

  auto solver = qbpp::ABS3Solver(f);
  auto sol = solver.search({{"time_limit", 10.0}});

  for (int m = 0; m < machines; ++m) {
    std::vector<std::pair<int, int>> order;
    for (int j = 0; j < n; ++j) order.push_back({sol(s[j][m]), j});
    std::sort(order.begin(), order.end());
    std::cout << "Machine " << m + 1 << ":";
    for (auto [start, j] : order)
      std::cout << " J" << j + 1 << "(" << start << "-" << start + time[j][m] << ")";
    std::cout << std::endl;
  }
  std::cout << "Makespan: " << sol(z) << std::endl;
}

route[j] はジョブ j が機械を回る順序、time[j][m] はジョブ j の機械 m での処理時間です。C++ の添字は 0 から始まるので、表の J1 と M1 はプログラムではジョブ 0 と機械 0 です。
T - time は要素ごとに $T - p_{j,m}$ を並べた整数の配列です。0 <= qbpp::int_var("s", n, machines) <= T - time で、要素ごとに上限が異なる整数変数の配列を作ります（整数変数の「要素ごとに範囲が異なる配列」）。
x は 6 × 6 × 6 のバイナリ変数の配列ですが、使うのは i < j の要素だけです。f に現れない変数はソルバーには渡りません。
ABS3Solver で 10 秒間探索し、機械ごとに作業を開始時刻の順に並べて出力します。

このプログラムの出力は次のとおりです。

Machine 1: J1(1-4) J4(13-18) J3(18-27) J6(28-38) J2(38-48) J5(48-51)
Machine 2: J2(0-8) J4(8-13) J6(13-16) J1(16-22) J5(22-25) J3(27-28)
Machine 3: J1(0-1) J3(1-6) J2(8-13) J5(13-22) J4(22-27) J6(42-43)
Machine 4: J3(6-10) J6(16-19) J4(27-30) J1(30-37) J2(48-52) J5(52-53)
Machine 5: J2(13-23) J5(25-30) J4(30-38) J6(38-42) J1(42-48) J3(48-55)
Machine 6: J3(10-18) J6(19-28) J2(28-38) J1(38-41) J5(41-45) J4(45-54)
Makespan: 55

メイクスパンは 55 で、知られている最小値と一致します。メイクスパンが 55 のスケジュールの中で、開始時刻の和（874）が最小のものはこの 1 つだけです。たとえば機械 1 の J1 が時刻 1 から始まるのは、J1 が最初に機械 3 で時刻 0〜1 の作業をするからです。

大きな問題について

ft06 は小さな問題です。ジョブと機械の数が増えると、順序の変数は機械ごとにジョブの組の数だけ必要になり、探索も難しくなります。処理時間をランダムに決めた 10 ジョブ × 10 機械の問題（最小のメイクスパンは 794）では、同じ書き方で 30 秒探索しても制約を満たすスケジュールが見つからないことがあり、2 分探索しても最小値より 1 割ほど長いスケジュールでした。 QUBO++ でジョブショップスケジューリングを書けることはこのページのとおりですが、大きな問題で最適に近いスケジュールを求めるには、ジョブショップスケジューリング専用の手法のほうが向いています。