動画：納期付き機械スケジューリング（Python）

動画「QUBO++ ケーススタディ 納期付き機械スケジューリング（Python）」の手引きです。 章ごとに、要点と動画に出てくるプログラム・出力、詳しい説明のページへのリンクをまとめています。 「▶ 動画の … から再生」をクリックすると、その場面から再生します。

QUBO++ ケーススタディ 納期付き機械スケジューリング（Python）（14:09・YouTube）

この動画で扱うこと：

• 3 台の機械と 12 件の注文、納期を守るスケジュールのうち処理時間の合計が最小のもの
• 各機械は注文を納期の早い順に処理してよいこと（隣どうしの入れ替えで示す）と、決めるのは割り当てだけであること
• 割り当てのバイナリ変数 $y_{j,m}$、割り当ての制約と納期の制約を qbpp.cons() で書くこと（不等式でも補助の変数は要らない）
• 配列の演算（要素ごとの積・qbpp.sum()・qbpp.vector_sum()）と、値が決まっている変数の qbpp.replace() による固定
• メイクスパン（最後の注文が終わる時刻）を整数変数 z で最小にすること

文章での説明は納期付き機械スケジューリングのページにあります。ブラウザ上の機械スケジューリングのデモでも試せます。動画（ケーススタディ）の一覧に戻る

問題

▶ 動画の 0:43 から再生

• ある工場に 3 台の機械と 12 件の注文があります。注文の番号は、あらかじめ納期の早い順に付けてあります（時間の単位は分。– はその機械では処理できないことを表します）。

注文	納期	機械 1	機械 2	機械 3
1	100	100	100	–
2	170	95	–	160
3	190	95	155	155
4	225	85	120	135
5	250	–	75	–
6	300	110	135	175
7	320	90	–	140
8	355	85	95	–
9	360	40	60	80
10	405	65	–	–
11	480	80	130	145
12	500	90	150	–

• 各注文はどれか 1 台の機械で処理します。機械は時刻 0 から、一度に 1 件ずつ、休みなく続けて処理します。どの注文も納期までに終えなければなりません。
• すべての注文をいちばん速い機械に置くと（注文 1 は機械 1 と機械 2 で同じ速さなので機械 1）、12 件のうち 11 件が機械 1 に集まり、番号の順に処理すると 12 件中 10 件が納期に遅れます。

何を決めるか

▶ 動画の 1:58 から再生

• 決めるのは、各注文をどの機械に置くかという割り当てと、各機械での処理の順序です。
• 処理時間の合計は、各注文の、割り当てた機械での処理時間の和なので、順序によりません。
• 納期の早い順に処理してよい理由：機械 1 で注文 10、9、8 をこの順に処理すると（納期は 405、360、355）、3 件とも間に合っています。隣どうしで前の注文のほうが納期が遅い組を入れ替えると、後ろに移った注文は、もともと後ろにあった注文と同じ時刻に終わります。その時刻はもともと後ろにあった注文の納期以前で、その納期は後ろに移った注文の納期より早いので、間に合います。前に移った注文は前より早く終わるので、やはり間に合います。入れ替えを 3 回くり返すと番号の順の 8、9、10 になります。
• つまり、間に合う順序があれば、納期の早い順でも間に合います。そこで各機械は番号の順に処理すると決め、割り当てだけを考えます。
• 割り当ては、処理できる機械の数を掛けた 7776 通りで、そのうちすべての納期を守れるのは 29 通りです。

変数と割り当ての制約

▶ 動画の 3:37 から再生

• 注文 $j$ を機械 $m$ で処理するとき 1、そうでないとき 0 のバイナリ変数 $y_{j,m}$（$j$ は 1〜12、$m$ は 1〜3）で割り当てを表します。処理できない 9 つの組では $y_{j,m} = 0$ です。
• 目的は処理時間の合計 $\sum_{j,m} p_{j,m}\, y_{j,m}$ の最小化です（$p_{j,m}$ は注文 $j$ の機械 $m$ での処理時間）。
• 1 つ目の制約は、各注文をちょうど 1 台の機械で処理すること、つまり各行の変数の和が 1 に等しい $\sum_m y_{j,m} = 1$ で、12 件の注文すべてに付けます。

終了時刻と納期

▶ 動画の 4:33 から再生

• 各機械は番号の順に処理するので、機械 $m$ で注文 $j$ が終わる時刻は、処理時間と変数の積を番号の順に足した $C_{j,m} = \sum_{i \le j} p_{i,m}\, y_{i,m}$ です。
• 各注文の終了時刻が納期 $d_j$ 以下という制約 $C_{j,m} \le d_j$ を、機械と注文の 36 通りの組すべてに付けます。
• 注文 $j$ が機械 $m$ に割り当てられていないときも、制約は正しく働きます。そのときの $C_{j,m}$ は機械 $m$ にある直前の注文が終わる時刻で、直前の注文の納期は $d_j$ 以前だからです。そのため、割り当てに合わせて制約を付けたり外したりしなくて済みます。

\[\begin{aligned} \text{最小化}\quad & \sum_{j,m} p_{j,m}\, y_{j,m} \\ \text{制約}\quad & \sum_{m} y_{j,m} = 1 && \text{（すべての注文 } j\text{）}\\ & C_{j,m} = \sum_{i \le j} p_{i,m}\, y_{i,m} \le d_j && \text{（すべての機械 } m \text{ と注文 } j\text{）} \end{aligned}\]

QUBO++での書き方

▶ 動画の 5:55 から再生

• qbpp.var("y", shape=(n, machines)) で 12 行 3 列の変数の配列 y を作ります。添字は 0 から始まるので、y[j][m] は式の $y_{j+1,m+1}$ を表します（y[0][0] は $y_{1,1}$、y[11][1] は $y_{12,2}$）。
• 制約は式を qbpp.cons() で囲みます。割り当ての制約では、qbpp.vector_sum(y, axis=1) で各行の和を並べた配列を作り、それが 1 に等しいという式を囲むと、12 件の制約がまとめてできます。
• 納期の制約は、終了時刻の式 finish が due[j] 以下という式を囲みます。
• 制約を破ると、重みと違反量の 2 乗を掛けた値がエネルギーに加わります。重みには処理時間の合計よりずっと大きい 100000 を使います。不等式の制約でも、qbpp.cons() なら補助の変数は要りません。

詳しくは制約と総和関数のページをご覧ください。

プログラム

▶ 動画の 7:25 から再生

• due は納期のリスト、time は処理時間の 2 次元のリストで、処理時間が 0 の組はその機械で処理できないことを表します。
• qbpp.array(time) と y の積 work は要素ごとの積で、work[j][m] は機械 m で処理するときは処理時間、処理しないときは 0 です。qbpp.sum(work) が処理時間の合計です。
• 納期の制約は、機械ごとに番号の順に finish に work の要素を足し、足すたびに finish <= due[j] を加えます。
• 値が決まっている変数を固定します。処理できない組の変数は 0、処理できる機械が 1 台だけの注文（注文 5 と 10）はその機械の変数を 1 にします。その値を ml に集め、qbpp.replace() で固定した式 g を作ると、固定した 11 個の変数が消え、ソルバーが探す変数は 25 個になります。
• g を simplify_as_binary() で整理し、EasySolver で 1 秒間探します。固定した変数は g の解 sol にはないので、qbpp.Sol(f).set(sol, ml) で固定した値を戻した解 full_sol を作って読みます。
• 最後のループで、機械ごとに 1 になった変数の注文を表示します。添字に 1 を足して注文と機械の番号に戻し、開始時刻はその機械で前に表示した注文の処理時間の和です。

import pyqbpp as qbpp

# 納期と、機械 1, 2, 3 での処理時間（分。0 はその機械では処理できない）
# 注文は納期の早い順に並べてある
due = [100, 170, 190, 225, 250, 300, 320, 355, 360, 405, 480, 500]
time = [[100, 100, 0], [95, 0, 160], [95, 155, 155], [85, 120, 135],
        [0, 75, 0], [110, 135, 175], [90, 0, 140], [85, 95, 0],
        [40, 60, 80], [65, 0, 0], [80, 130, 145], [90, 150, 0]]
n = len(due)
machines = 3

y = qbpp.var("y", shape=(n, machines))
work = qbpp.array(time) * y
f = qbpp.sum(work)

# 各注文をちょうど 1 台の機械で処理する
f += 100000 * qbpp.cons(qbpp.vector_sum(y, axis=1) == 1)

# 各注文は納期までに終わる
for m in range(machines):
    finish = 0
    for j in range(n):
        finish += work[j][m]
        f += 100000 * qbpp.cons(finish <= due[j])

# 処理できない組は 0 に、処理できる機械が 1 台だけの注文は 1 に固定する
ml = {}
for j in range(n):
    for m in range(machines):
        if time[j][m] == 0:
            ml[y[j][m]] = 0
        elif time[j].count(0) == machines - 1:
            ml[y[j][m]] = 1
g = qbpp.replace(f, ml)
g.simplify_as_binary()

solver = qbpp.EasySolver(g)
sol = solver.search(time_limit=1.0)
full_sol = qbpp.Sol(f).set(sol, ml)

total = 0
for m in range(machines):
    print(f"Machine {m + 1}:", end="")
    t = 0
    for j in range(n):
        if full_sol(y[j][m]) == 1:
            print(f" #{j + 1}({t}-{t + time[j][m]})", end="")
            t += time[j][m]
            total += time[j][m]
    print()
print("Total processing time:", total)

詳しくは配列演算・置換関数・変数の固定と番兵のページをご覧ください。

実行結果

▶ 動画の 9:36 から再生

出力：

Machine 1: #2(0-95) #3(95-190) #6(190-300) #9(300-340) #10(340-405) #12(405-495)
Machine 2: #1(0-100) #5(100-175) #8(175-270) #11(270-400)
Machine 3: #4(0-135) #7(135-275)
Total processing time: 1170

• 機械 1 は注文 2、3、6、9、10、12 を 495 分までに、機械 2 は注文 1、5、8、11 を、機械 3 は注文 4 と 7 を処理し、どの注文も納期の手前で終わります。
• 処理時間の合計は 1170 分です。すべてを速い機械に置いた 1010 分より長いですが、納期を守っています。
• 7776 通りの割り当てをすべて調べると合計の最小は 1170 分で、最小になる割り当ては 1 つだけ、出力の割り当てと一致します。
• ただし、機械 1 は 495 分まで動き、機械 3 は 275 分で止まっています。

メイクスパン

▶ 動画の 10:38 から再生

• 最後の注文が終わる時刻をメイクスパンと呼びます。各機械が動いている時間（その機械に割り当てた注文の処理時間の合計）を負荷と呼ぶと、メイクスパンはいちばん大きい負荷です。
• 0〜600 の整数の値をとる変数 z を qbpp.var("z", integer=(0, 600)) で作り、どの機械の負荷も z 以下という制約を 3 台の機械に付け、目的関数に 100 * z を加えます。z は最大の負荷まで下がり、さらに下げるには負荷がそろうように割り当てが変わります。
• 重みを 100 にした理由：時間はすべて 5 分の倍数なので、メイクスパンが 5 分短くなるとエネルギーは 500 下がります。一方、処理時間の合計はどう割り当てても 465 分（注文ごとの最も遅い処理時間と最も速い処理時間の差の合計）までしか変わりません。500 のほうが大きいので、まずメイクスパンが最小になり、メイクスパンが同じ割り当ての中で処理時間の合計が最小になります。
• プログラムでは、納期の制約の後に次の行を加えます。work の各列の和が各機械の負荷なので、qbpp.vector_sum(work, axis=0) で負荷の配列を作り、z を引いた値が 0 以下という制約をまとめて加えます。最後に print("Makespan:", full_sol(z)) で z を表示します。

# メイクスパン z: どの機械も z までに終わる
z = qbpp.var("z", integer=(0, 600))
f += 100 * z
f += 100000 * qbpp.cons(qbpp.vector_sum(work, axis=0) - z <= 0)

出力：

Machine 1: #2(0-95) #3(95-190) #6(190-300) #9(300-340) #10(340-405)
Machine 2: #1(0-100) #5(100-175) #8(175-270) #12(270-420)
Machine 3: #4(0-135) #7(135-275) #11(275-420)
Makespan: 420
Total processing time: 1245

• メイクスパンは 420 分になり、3 台の機械の負荷がそろいます。その代わり、処理時間の合計は 1170 分から 1245 分に増えます。
• すべての割り当てを調べると 420 分が最小で、その割り当ては 1 つです。納期を無視すれば 410 分にできるので、納期を守るために 10 分延びています。

詳しくは整数変数のページをご覧ください。

まとめ

▶ 動画の 13:13 から再生

• 各機械は納期の早い順に処理すると決め、割り当てだけを変数にしました。終了時刻は、納期の早い順に処理時間と変数の積を足した式で表せます。
• 割り当ての制約と納期の制約は qbpp.cons() でそのまま書けました。処理時間の合計や各行の和は配列の演算で短く書け、値が決まっている変数は qbpp.replace() で固定しました。
• メイクスパンは整数変数 z と、負荷が z 以下という制約で最小にしました。重みの大きさで、2 つの目的に優先順位を付けられます。