問題: n個のジョブを1つの機械で任意順に処理する。同時に取り組めるジョブは1つであり、1つのジョブに着手したら途中で他のジョブに切り替えることはできない。ジョブiを始めてから完了するまでの時間はpi (>0)であり、ジョブiを始めるのは時刻ri (≥0)以降でなければならない。また、各ジョブには重み wi>0が与えられており、ジョブの処理を完了する時刻の重み付き総和を最小化したい。つまり、ジョブiを時刻tiから始めるとした時、∑iwi(ti+pi)を最小化したい。最適なt1,...,tnを求めよ。
ジョブスケジューリング問題の一種であり、分類としては1∣rj∣∑wjCjという表記で表される。NP-hardである。
これをMIPとして定式化したい。自然な考え方としては、ジョブiにジョブj以前に取り組むことに対応するバイナリ変数xi,jを使った定式化が考えられる。目的関数は定数項を無視すると∑iwitiであり、制約は以下のようになる:
- 任意のジョブjについて tj≥rj
- 任意の異なるジョブi,jについて tj≥ti+pi−Mxj,i
- 任意の異なるジョブi,jについて xi,j+xj,i=1
- 任意の異なるジョブi,j,kについて xi,j+xj,k≤1+xi,k
- 任意の異なるジョブi,jについて xi,j∈{0,1}
- 任意のジョブjについて tj∈R
これで定式化になっているが、NemhauserとSavelsberghによるとbig Mを使わない定式化が可能らしい。一般性を失わずにr1≤r2≤...≤rnと仮定するとき、tj≥rjとtj≥ti+pi−Mxj,iの代わりに
tj≥rixi,j+k<i,k=j∑pk(xi,k+xk,j−1)+k≥i,k=j∑pkxk,j(*)
とbig Mを使わずに書けるというのが論文の内容である。ただし、便宜上xi,i=1とする。上のbig Mを使った定式化をFormM、(∗)を使った定式化をForm*とする。Form*の正当性を確かめたい。
ジョブを処理する順番を固定して考えることにする。具体的には、次のことを示す:
定理: 1,...,nの順列π=(π(1),...,π(n))が与えられたとする。冒頭の問題に「ジョブをπの順で処理する」という条件を追加したものの最適解の集合は、Form*のすべてのxi,jを順列に従った値を取る定数としたLPの最適解の集合と一致する。
これが示されれば、
冒頭の問題の最適解の集合
= ⋃π:順列(ジョブをπの順で処理するという条件の下での冒頭の問題の最適解の集合) の要素で目的関数を最小化するものの集合
= ⋃π:順列(Form*のすべてのxi,jを順列に従った値を取る定数としたLPの最適解の集合) の要素で目的関数を最小化するものの集合
= Form*の最適解の集合
となってForm*は正しい定式化であると言える。
定理の証明
以降、ジョブjをpos(j)番目に処理するとする。pos(π(k))=kである。
まず、ジョブをπ(1),...,π(n)の順で処理するという条件下でのスケジューリング問題の最適解(t1,...,tn)=(T1,...,Tn)の性質を考える。明らかに、可能な最も早い時刻にジョブを順次処理していくのが最適だから、Tπ(k+1)=max(rπ(k+1),Tπ(k)+pπ(k))である。これを利用してTjをrとpのみで表すことを考える。
まず、Tπ(1)=rπ(1)である。また、π(1)以外の任意のジョブjについては、jより前に行うジョブでリリース時刻ちょうどに処理を始めるものが必ず存在する。その中で最後に処理するジョブをj′とするとき
Tj=max(rj,rj′+k:pos(j′)≤pos(k)<pos(j)∑pk)
である。したがって、次の3つの命題が示せればよい。以下、(∗)の右辺を Ri,j:=rixi,j+∑k<i,k=jpk(xi,k+xk,j−1)+∑k≥i,k=jpkxk,j とする。ただし、ここではxi,jは与えられた順列πに応じて値が定まっている定数である。
命題1: ジョブj=π(1)が与えられたとする。jより前に行いかつリリース時刻ちょうどに処理を始めるジョブで最後に処理するものをj′とするとき、Rj′,j=rj′+∑k:pos(j′)≤pos(k)<pos(j)pkである。
命題2: 任意のジョブjについて、Tj=rjならばRj,j=rjである。
命題3: 任意のジョブi,jについてRi,j≤Tjである。
命題1と2はForm*がtjの下限に関するタイトな制約をすべて含むと主張している。命題3はそれ以外の制約が冗長であると主張している。
命題1の証明
Rj′,j=rj′+∑k<j′,k=jpk(xj′,k+xk,j−1)+∑k≥j′,k=jpkxk,jである。任意のジョブλ=jについて、pλがRj′,jにどう現れるか調べる:
- pos(λ)<pos(j′)の場合
- λ<j′ならpλの係数はxj′,k+xk,j−1=0
- λ≥j′であることはない。なぜならジョブj′はrj′ちょうどに開始するため、pos(λ)<pos(j′)ならrλ<rj′でなければならず、つまりλ<j′だから。
- pos(j′)≤pos(λ)<pos(j)の場合
- λ<j′ならpλの係数はxj′,λ+xλ,j−1=1
- λ≥j′ならpλの係数はxλ,j=1
- pos(j)<pos(λ)の場合
- λ<j′ならpλの係数はxj′,λ+xλ,j−1=0
- λ≥j′ならpλの係数はxλ,j=0
したがって、Rj′,j=rj′+∑λ:pos(j′)≤pos(λ)<pos(j)pλ。
命題2の証明
Rj,j=rjxj,j+∑k<jpk(xj,k+xk,j−1)+∑k>jpkxk,jである。
k=jなら常に xj,k+xk,j−1=0 なので ∑k<jpk(xj,k+xk,j−1)=0。次に∑k>jpkxk,jについて考える。k>jならばrk≥rjだから、jをrjちょうどに開始するという仮定の下ではkは必ずjの後に行うことになる。つまり、xk,j=0。従って総和もゼロであり、Rj,j=rjxj,j=rj。
命題3の証明
pos(i)≤pos(j)の場合とpos(j)<pos(i)の場合にわけて検討する。
まず、pos(i)≤pos(j)の場合を扱う。任意のジョブλ=jについて、命題1の証明と同じようにpλがRi,jにどう現れるか調べる。
- pos(λ)<pos(i)の場合
- λ<iならpλの係数は xi,λ+xλ,j−1=0
- λ≥iならpλの係数は xλ,j=1 ――(A)
- pos(i)≤pos(λ)<pos(j)の場合、命題1の証明と同様にpλの係数は1である
- pos(j)<pos(λ)の場合、命題1の証明と同様にpλの係数は0である
(A)に該当するλが存在しなければ Ri,j=ri+∑λ:pos(i)≤pos(λ)<pos(j)pλ だから、次のことが成り立つ。
補題1: pos(i)≤pos(j)とする。(A)に該当するジョブが存在しなければRi,j≤Tjである。
また、(A)に該当するλが存在する場合は実は次のことが成り立つ:
補題2: pos(i)≤pos(j)とし、(A)に該当するジョブが存在するとする。そのようなジョブの中でもっとも早く行うものλに注目する。言い換えると、λはλ≥iであるジョブの中でもっとも早く行うとする。このとき、Ri,j≤Rλ,jである。
証明: pos(i)≤pos(j)に注意すると、Rλ,j,Ri,jは以下のように表せる:
Rλ,j=rλ+k<λ,pos(λ)<pos(k)<pos(j)∑pk+k≥λ,pos(k)<pos(j)∑pk
Ri,j=ri+k<i,pos(i)<pos(k)<pos(j)∑pk+k≥i,pos(k)<pos(j)∑pk
Ri,jにどのようなpkが現れるか調べるために、総和に登場する添字kの集合を分割する。次の1～4のうち、1が一つ目の総和のkの集合で、2～4が二つ目の総和のkの集合の分割である。
- {k∈[n]:k<i,pos(i)<pos(k)<pos(j)}
- {k∈[n]:i≤k<λ,pos(k)<pos(λ)}
- {k∈[n]:i≤k<λ,pos(λ)<pos(k)<pos(j)}
- {k∈[n]:k≥λ,pos(k)<pos(j)}
i≤λであるλの中でもっとも早く行うものを選んでいるので2は空である。また、1と3の合併は{k∈[n]:k<λ,pos(λ)<pos(k)<pos(j)}に包含され、後者はRλ,jの一つ目の総和の添字集合に等しい。また、4はRλ,jの二つ目の総和の添字集合に等しい。従って、Ri,j≤Rλ,jである。□
次に、pos(j)<pos(i)の場合を検討する。Ri,j=∑k<i,k=jpk(xi,k+xk,j−1)+∑k≥i,k=jpkxk,jである。やはり任意のジョブλ=jについて、pλがRi,jにどう現れるか調べる。
- pos(λ)<pos(j)の場合
- λ<i ならpλの係数はxi,λ+xλ,j−1=0
- λ≥i ならpλの係数はxλ,j=1 ――(B)
- pos(j)<pos(λ)<pos(i)の場合
- λ<iならpλの係数はxi,λ+xλ,j−1=−1
- λ≥iならpλの係数はxλ,j=0
- pos(i)≤pos(λ)の場合
- λ<iならpλの係数はxi,λ+xλ,j−1=0
- λ≥iならpλの係数はxλ,j=0
(B)に該当するλが存在しなければ、Ri,j≤0≤Tjだから、次のことが成り立つ:
補題3: pos(j)<pos(i)とする。(B)に該当するλが存在しなければRi,j≤Tjである。
また、(B)に該当するλが存在する場合は次のことが成り立つ:
補題4: pos(j)<pos(i)とする。(B)に該当するジョブが存在するとする。そのようなジョブの中でもっとも早く行うものλに注目する。言い換えると、λはλ≥iであるジョブの中でもっとも早く行うとする。このとき、Ri,j≤Rλ,jである。
証明: pos(λ)<pos(j)<pos(i)に注意すると、Rλ,j,Ri,jは以下のように評価できる:
Rλ,j=k<λ,pos(λ)<pos(k)<pos(j)∑pk+k≥λ,pos(k)<pos(j)∑pk
Ri,j=−k<i,pos(j)<pos(k)<pos(i)∑pk+k≥i,pos(k)<pos(j)∑pk≤k≥i,pos(k)<pos(j)∑pk
Ri,jの上界∑k≥i,pos(k)<pos(j)pkにどのようなpkが現れるか調べるために、添字kの集合を3つに分割する:
- {k∈[n]:i≤k<λ,pos(k)<pos(λ)}
- {k∈[n]:i≤k<λ,pos(λ)<pos(k)<pos(j)}
- {k∈[n]:k≥λ,pos(k)<pos(j)}
i≤λであるジョブλの中でもっとも早く行うものを選んでいるので1は空である。また、2はRλ,jの一つ目の総和の添字集合に包含され、3は二つ目の総和の添字集合に等しい。従って、Ri,j≤Rλ,jである。□
補題1～4を適用することで、命題3が示せる: 補題1と3より、(A)や(B)に該当するようなλが存在しない(i,j)に対してはRi,j≤Tjである。また、(A)や(B)に該当するようなλが存在する場合、どちらの場合もi<λだから、補題2または4を繰り返し適用することで、i<λ1<λ2<..., Ri,j≤Rλ1,j≤Rλ2,j≤... という列が得られる。高々有限回の適用で(A)や(B)に該当する項がないものが見つかり、補題1または3によってRi,j≤Tjが従う。
これでForm*の正当性が示された。
上のNemhauserとSavelsberghの論文ではLP双対を経由した非常に短い「証明」が書かれている……が、読んでも双対を取ることで問題が自明になっていると思えなかった。双対に注目しても、ギャップを埋めようとすればRi,jのうちどれがどれより大きいのかという細かい議論を避けられないように思える。
実験
big Mを避けることで実際に速くなるかどうか知りたい。FormMとForm*をSCIPに与えて、ランダムインスタンスで比較した。
import math
import random
import statistics
from dataclasses import dataclass
from itertools import pairwise, product
from typing import Any
import pandas as pd
import pulp
@dataclass
class SchedulingInstance:
n: int
release_times: list[float]
weights: list[float]
processing_times: list[float]
@dataclass
class SchedulingSolution:
n: int
start_times: list[float]
model: pulp.LpProblem
def validate_solution(
instance: SchedulingInstance, solution: SchedulingSolution
) -> None:
n = instance.n
release_times = instance.release_times
processing_times = instance.processing_times
start_times = solution.start_times
eps = 1e-8
for i in range(n):
assert release_times[i] <= start_times[i] + eps
for i1, i2 in pairwise(sorted(range(n), key=lambda i: start_times[i])):
assert start_times[i1] + processing_times[i1] <= start_times[i2] + eps
assert math.isclose(
max(start_times[i1] + processing_times[i1], release_times[i2]),
start_times[i2],
rel_tol=eps,
abs_tol=eps,
)
def make_random_instance(n: int, mag: float, seed: int) -> SchedulingInstance:
rng = random.Random(seed)
release_times = sorted((rng.randrange(0, mag) for _ in range(n)))
weights = [rng.randrange(1, mag) for _ in range(n)]
processing_times = [rng.randrange(1, mag) for _ in range(n)]
return SchedulingInstance(n, release_times, weights, processing_times)
def solve(
instance: SchedulingInstance, is_big_m: bool, time_out_sec: float = 600.0
) -> SchedulingSolution | None:
n = instance.n
release_times = instance.release_times
weights = instance.weights
processing_times = instance.processing_times
model = pulp.LpProblem("single_machine_scheduling", pulp.LpMinimize)
start_times = [pulp.LpVariable(f"start_time_{i}") for i in range(n)]
is_priors = [
[pulp.LpVariable(f"is_prior_{i}_{j}", cat=pulp.LpBinary) for i in range(n)]
for j in range(n)
]
model += pulp.lpSum(start_times[i] * weights[i] for i in range(n))
def is_prior(i: int, j: int) -> pulp.LpAffineExpression:
assert i != j
if i > j:
return pulp.LpAffineExpression(1 - is_priors[j][i])
else:
return pulp.LpAffineExpression(is_priors[i][j])
for i, j, k in product(range(n), range(n), range(n)):
if i == j or j == k or k == i:
continue
model += is_prior(i, j) + is_prior(j, k) <= 1 + is_prior(i, k)
for i in range(n):
model += start_times[i] >= release_times[i]
if is_big_m:
big_m = 2 * (release_times[-1] + sum(processing_times)) + 1
for i, j in product(range(n), range(n)):
if i == j:
continue
model += start_times[j] >= start_times[i] + processing_times[
i
] - big_m * is_prior(j, i)
else:
for i, j in product(range(n), range(n)):
if i == j:
continue
model += start_times[j] >= release_times[i] * is_prior(i, j) + pulp.lpSum(
processing_times[k] * (is_prior(i, k) + is_prior(k, j) - 1)
for k in range(i)
if k != j
) + pulp.lpSum(
processing_times[k] * is_prior(k, j) for k in range(i, n) if k != j
)
solver = pulp.SCIP_CMD(timeLimit=time_out_sec)
model.solve(solver)
print(f"{model.solutionTime=}")
if model.status == pulp.LpStatusNotSolved:
return None
assert model.status == pulp.LpStatusOptimal, pulp.LpStatus[model.status]
solution_start_times = [pulp.value(start_times[i]) for i in range(n)]
sol = SchedulingSolution(n, solution_start_times, model)
validate_solution(instance, sol)
return sol
def time_with_param(n: int, mag: float, is_big_m: bool, n_iter: int) -> dict[str, Any]:
times: list[float] = []
objs: list[float] = []
for seed in range(n_iter):
instance = make_random_instance(n, mag, seed)
sol = solve(instance, is_big_m, time_out_sec=600.0)
times.append(sol.model.solutionTime if sol else math.inf)
objs.append(pulp.value(sol.model.objective) if sol else math.nan)
print(objs)
return (
{"n": n, "method": "big M" if is_big_m else "no M"}
| {f"seed {i}": times[i] for i in range(n_iter)}
| {"median": statistics.median(times)}
)
df = pd.DataFrame(
[
time_with_param(10, 10, True, 5),
time_with_param(15, 15, True, 5),
time_with_param(10, 10, False, 5),
time_with_param(15, 15, False, 5),
time_with_param(20, 20, False, 5),
time_with_param(30, 30, False, 5),
time_with_param(50, 50, False, 5),
],
)
print(df.to_markdown(index=False, floatfmt=".2f"))
結果は次のようになった。infは600秒以内に解けなかったことを表す。
| n |
method |
seed 0 |
seed 1 |
seed 2 |
seed 3 |
seed 4 |
median |
| 10 |
big M |
2.26 |
0.94 |
3.95 |
1.30 |
9.83 |
2.26 |
| 15 |
big M |
inf |
298.36 |
inf |
233.59 |
inf |
inf |
| 10 |
no M |
0.21 |
3.47 |
2.25 |
1.20 |
0.07 |
1.20 |
| 15 |
no M |
4.04 |
4.07 |
2.56 |
0.70 |
4.65 |
4.04 |
| 20 |
no M |
3.29 |
0.32 |
3.60 |
0.36 |
3.41 |
3.29 |
| 30 |
no M |
12.25 |
4.03 |
5.63 |
5.40 |
14.96 |
5.63 |
| 50 |
no M |
20.19 |
65.13 |
inf |
inf |
37.51 |
65.13 |
big Mを避けた定式化のほうが速く解けると言える。しかし、この定式化であっても、O(2nn)の自明なDPで扱える規模(n<30)を超えてSCIPで新たに解けるようになる範囲は+10～20程度のようだ。ただし、上の定式化は最低限の制約しか使っておらず、冗長な制約（妥当不等式）を追加していくなどによってさらに速くなる可能性はある。妥当不等式に関しては、NemhauserとSavelsberghの論文にも議論がある。
その他
NemhauserとSavelsberghのテクニックは、機械が1つであれば冒頭の問題設定以外にも様々なリリース時刻付き問題に使える。例えば、ジョブに依存関係がある1∣prec;rj∣∑wjCjや、納期に対する遅れ時間maxi(Ti−di)を最小化する1∣rj∣Lmaxなど。
Khowalaが単一マシンの様々なジョブスケジューリング問題と定式化のベンチマークを行っていて、ここで扱ったような処理順序を表すバイナリ変数を導入するやり方以外の定式化も知ることができる。
この問題に関しては、ジョブの順番に制約がなくて順番を決めれば最適解が自明に求まるという性質上、局所探索が明らかに簡単かつ有効ではある。実用的にはサイズがかなり大きくなっても困らないタイプの問題と言えるかもしれない。（実際はだいたい面倒な制約が付くだろうけど……）
文献
- Brucker, Peter, Jan Karel Lenstra and Alexander H. G. Rinnooy Kan. “Complexity of machine scheduling problems.” (1975).
- Khowala, Ketan. “Single Machine Scheduling: Comparison of MIP Formulations and Heuristics for Interfering Job Sets.” (2012).
- Nemhauser, George L. and Mwp Martin Savelsbergh. “A Cutting Plane Algorithm for the Single Machine Scheduling Problem with Release Times.” (1992).
