1機械リリース時刻付き重み付き完了時刻和最小化問題
1機械リリース時刻付き重み付き完了時刻和最小化問題 $1|r_j|\sum w_j C_j$ は，以下に定義される問題である．
単一の機械で $n$ 個のジョブを処理する問題を考える． この機械は一度に1つのジョブしか処理できず，ジョブの処理を 開始したら途中では中断できないものと仮定する． ジョブの集合の添え字を $j=1,2,\ldots,n$，ジョブの集合を $J=\{1,2,\ldots,n\}$ と書く． 各ジョブ$j$ に対する処理時間 $p_j$，重要度を表す重み$w_j$， リリース時刻（ジョブの処理が開始できる最早時刻）$r_j$が与えられたとき， 各ジョブ $j$ の処理完了時刻 $C_j$ の重み付きの 和を最小にするようなジョブを機械にかける順番（スケジュール）を求めよ．
この問題は古くから多くの研究がなされているスケジューリング問題であり，$NP$-困難であることが知られている． ここでは，離接定式化(disjunctive formulation)と呼ばれる単純な定式化を示す．
ジョブ $j$ の開始時刻を表す連続変数 $s_{j}$ と，ジョブ $j$ が他のジョブ $k$ に先行する（前に処理される）とき $1$ になる $0$-$1$ 整数変数 $x_{jk}$ を用いる．
離接定式化は，非常に大きな数 $M$ を用いると，以下のように書ける．
$$ \begin{array}{l l l} minimize & \displaystyle\sum_{j=1}^n w_j s_{j} + \displaystyle\sum_{j=1}^n w_j p_j & \\ s.t. & s_j + p_j -M (1-x_{jk}) \leq s_k & \forall j \neq k \\ & x_{jk} +x _{kj}= 1 & \forall j < k \\ & s_j \geq r_j & \forall j=1,2,\ldots,n \\ & x_{jk} \in \{ 0,1 \} & \forall j \neq k \end{array} $$目的関数は重み付き完了時刻の和 $\sum w_j (s_j+p_j)$ を展開したものであり，それを最小化している． 展開した式の第2項は定数であるので，実際には第1項のみを最小化すれば良い． 最初の制約は，$x_{jk}$ が $1$ のとき，ジョブ$j$ の完了時刻 $s_j+p_j$ より ジョブ$k$ の開始時刻 $s_k$ が後になることを表す． （$x_{jk}$ が $0$ のときには，$M$ が大きな数なので制約にはならない．） 2番目の制約は，ジョブ$j$ がジョブ$k$ に先行するか，ジョブ$k$ がジョブ$j$ に 先行するかの何れかが成立することを表す． これを離接制約(disjunctive constraint)と呼ぶ． これが離接定式化の名前の由来である． 3番目の制約は，ジョブ$j$ の開始時刻 $s_j$ がリリース時刻 $r_j$ 以上であることを表す．
大きな数 $M$ は，実際にはなるべく小さな値に設定すべきである． 制約ごとに異なる $M$ を設定するならば，$r_j +\sum s_j -s_k$， すべての制約で同じ $M$ を用いるならば $\max (r_j) +\sum s_j$ とすれば十分である．
以下では，上の離接定式化を用いて5ジョブの例題を求解する．
def make_data(n):
"""
Data generator for the one machine scheduling problem.
"""
p, r, d, w = {}, {}, {}, {}
J = range(1, n + 1)
for j in J:
p[j] = random.randint(1, 4)
w[j] = random.randint(1, 3)
T = sum(p)
for j in J:
r[j] = random.randint(0, 5)
d[j] = r[j] + random.randint(0, 5)
return J, p, r, d, w
def scheduling_disjunctive(J, p, r, w):
"""
scheduling_disjunctive: model for the one machine total weighted completion time problem
Disjunctive optimization model for the one machine total weighted
completion time problem with release times.
Parameters:
- J: set of jobs
- p[j]: processing time of job j
- r[j]: earliest start time of job j
- w[j]: weighted of job j, the objective is the sum of the weighted completion time
Returns a model, ready to be solved.
"""
model = Model("scheduling: disjunctive")
M = max(r.values()) + sum(p.values()) # big M
s, x = (
{},
{},
) # start time variable, x[j,k] = 1 if job j precedes job k, 0 otherwise
for j in J:
s[j] = model.addVar(lb=r[j], vtype="C", name="s(%s)" % j)
for k in J:
if j != k:
x[j, k] = model.addVar(vtype="B", name="x(%s,%s)" % (j, k))
model.update()
for j in J:
for k in J:
if j != k:
model.addConstr(
s[j] - s[k] + M * x[j, k] <= (M - p[j]), "Bound(%s,%s)" % (j, k)
)
if j < k:
model.addConstr(x[j, k] + x[k, j] == 1, "Disjunctive(%s,%s)" % (j, k))
model.setObjective(quicksum(w[j] * s[j] for j in J), GRB.MINIMIZE)
model.update()
model.__data = s, x
return model
n = 5 # number of jobs
J, p, r, d, w = make_data(n)
model = scheduling_disjunctive(J, p, r, w)
model.optimize()
s, x = model.__data
z = model.ObjVal + sum([w[j] * p[j] for j in J])
print("Opt.value by Disjunctive Formulation:", z)
seq = [j for (t, j) in sorted([(int(s[j].X + 0.5), j) for j in s])]
print("solition=", seq)
1機械総納期遅れ最小化問題
1機械総納期遅れ最小化問題 $1||\sum T_j$ は，以下に定義される問題である．
単一の機械で $n$ 個のジョブを処理する問題を考える． この機械は一度に1つのジョブしか処理できず，ジョブの処理を 開始したら途中では中断できないものと仮定する． ジョブの添え字を $1,2,\ldots,n$ と書く． 各ジョブ$j$ に対する 処理時間$p_j$，納期（ジョブが完了しなければならない時刻）$d_j$， が与えられたとき， 各ジョブ $j$ の納期からの遅れの和（総納期遅れ）を 最小にするようなジョブを機械にかける順番（スケジュール）を求めよ．
この問題は，$NP$-困難であることが知られている． ここでは，大きな数 $M$ を用いない定式化として線形順序付け定式化(linear ordering formulation)を考える．
ジョブ$j$ の納期遅れを表す連続変数 $T_{j}$と，ジョブ$j$ が他のジョブ$k$ に先行する（前に処理される）とき $1$ になる $0$-$1$ 整数変数 $x_{jk}$ を用いる．
上の記号を用いると，線形順序付け定式化は，以下のように書ける．
$$ \begin{array}{l l l} minimize & \displaystyle\sum_{j} w_j T_j & \\ s.t. & x_{jk} +x _{kj}= 1 & \forall j < k \\ & x_{jk}+x_{k \ell} +x_{\ell j} \leq 2 & \forall j \neq k \neq \ell \\ & \displaystyle\sum_{j=1}^n \displaystyle\sum_{k \neq j} p_k x_{kj} + p_j \leq d_j +T_j & \forall j=1,2,\ldots,n \\ & x_{jk} \in \{ 0,1 \} & \forall j \neq k \end{array} $$上の定式化において， 目的関数は，各ジョブ $j$ に対する納期遅れ $T_j$ を重み $w_j$ で乗じたものの和を 最小化することを表す．最初の制約は，ジョブ$j$ がジョブ $k$ の前に処理されるか， その逆にジョブ$k$ がジョブ$j$ の前に処理されるかのいずれかが成立することを表す． （ここまでは離接定式化と同じである．） 2番目の制約は，異なるジョブ$j,k,\ell$ に対して， ジョブ$j$ がジョブ$k$ に先行， ジョブ$k$ がジョブ$\ell$ に先行， ジョブ$\ell$ がジョブ$j$ に先行の3つが同時に成り立つことがないことを表す． この2つの制約によって，ジョブ上に線形順序（ジョブが一列に並べられること）が規定される． これが線形順序付け定式化の名前の由来である． 3番目の制約は，変数 $x$ で規定される順序付けと納期遅れを表す変数 $T$ を繋ぐための 制約であり，ジョブ$j$ の完了時刻が納期 $d_j$ より超えている時間が納期遅れ $T_j$ であることを規定する．
以下では，上の線形順序定式化を用いて5ジョブの例題を求解する．
def scheduling_linear_ordering(J, p, d, w):
"""
scheduling_linear_ordering: model for the one machine total weighted tardiness problem
Model for the one machine total weighted tardiness problem
using the linear ordering formulation
Parameters:
- J: set of jobs
- p[j]: processing time of job j
- d[j]: latest non-tardy time for job j
- w[j]: weighted of job j, the objective is the sum of the weighted completion time
Returns a model, ready to be solved.
"""
model = Model("scheduling: linear ordering")
T, x = {}, {} # tardiness variable, x[j,k] =1 if job j precedes job k, =0 otherwise
for j in J:
T[j] = model.addVar(vtype="C", name="T(%s)" % (j))
for k in J:
if j != k:
x[j, k] = model.addVar(vtype="B", name="x(%s,%s)" % (j, k))
model.update()
for j in J:
model.addConstr(
quicksum(p[k] * x[k, j] for k in J if k != j) - T[j] <= d[j] - p[j],
"Tardiness(%r)" % (j),
)
for k in J:
if k <= j:
continue
model.addConstr(x[j, k] + x[k, j] == 1, "Disjunctive(%s,%s)" % (j, k))
for ell in J:
if ell > k:
model.addConstr(
x[j, k] + x[k, ell] + x[ell, j] <= 2,
"Triangle(%s,%s,%s)" % (j, k, ell),
)
model.setObjective(quicksum(w[j] * T[j] for j in J), GRB.MINIMIZE)
model.update()
model.__data = x, T
return model
n = 5 # number of jobs
J, p, r, d, w = make_data(n)
model = scheduling_linear_ordering(J, p, d, w)
model.optimize()
z = model.ObjVal
x, T = model.__data
for (i, j) in x:
if x[i, j].X > 0.5:
print("x(%s) = %s" % ((i, j), int(x[i, j].X + 0.5)))
for i in T:
print("T(%s) = %s" % (i, int(T[i].X + 0.5)))
print("Opt.value by the linear ordering formulation=", z)
順列フローショップ問題
順列フローショップ問題(permutation flow shop problem) $F|prmu|C_{\max}$は，以下に定義される問題である．
$n$ 個のジョブを$m$台の機械で順番に処理することを考える． この機械は一度に1つのジョブしか処理できず，ジョブの処理を 開始したら途中では中断できないものと仮定する． いま，各機械におけるジョブの処理順序が同一であるとしたとき， 最後の機械 $m$ で最後に処理されるジョブの完了時刻を最小にする処理順を求めよ．
ここでは，この問題に対する位置データ定式化(positional data formulation) と呼ばれる定式化を示す．
いま，$n$個のジョブ（添え字は $j$）と $m$台の機械（添え字は $i$）があり， 各ジョブは，機械 $1,2,\ldots,m$ の順番で処理されるものとする． （これがフローショップと呼ばれる所以である．） ジョブ$j$ の機械$i$ 上での処理時間を $p_{ij}$ とする． ジョブの追い越しがないので，解はジョブの投入順序，いいかえれば順列で表現できる． （これが順列フローショップと呼ばれる所以である．） ジョブを並べたときの順番を $\kappa$ で表すことにし，順番を表す変数 $x_{j \kappa}$ を用いる．
$x_{j \kappa}$は，ジョブ $j$ が $\kappa$ 番目に投入されるとき $1$ になる $0$-$1$ 変数であり， $s_{i \kappa}$は，機械 $i$ の $\kappa$ 番目に並べられてるジョブの開始時刻を表す実数変数であり， $f_{i \kappa}$は，機械 $i$ の $\kappa$ 番目に並べられてるジョブの完了（終了）時刻を表す実数変数である．
この定式化は，ジョブ自身の開始時刻ではなく， 機械にかかるジョブの順番（位置）に対する開始時刻や終了時刻のデータを用いるので， 位置データ定式化と呼ばれる．
定式化は以下のようになる．
$$ \begin{array}{l l l} minimize & f_{m n} & \\ s.t. & \displaystyle\sum_{\kappa} x_{j\kappa} = 1 & \forall j=1,2,\ldots,n \\ & \displaystyle\sum_{j} x_{j\kappa} = 1 & \forall \kappa=1,2,\ldots,n \\ & f_{i \kappa} \leq s_{i,\kappa +1} & \forall i=1,2,\ldots,m, \kappa=1,2,\ldots,n-1 \\ & s_{i \kappa}+\displaystyle\sum_{j} p_{ij} x_{j \kappa} \leq f_{i \kappa} & \forall i=1,2,\ldots,m, \kappa=1,2,\ldots,n \\ & f_{i \kappa} \leq s_{i+1,\kappa} & \forall i=1,2,\ldots,m-1, \kappa=1,2,\ldots,n \\ & x_{j \kappa} \in \{ 0,1 \} & \forall j=1,2,\ldots,n, \kappa=1,2,\ldots,n \\ & s_{i \kappa} \geq 0 & \forall i=1,2,\ldots,m, \kappa=1,2,\ldots,n \\ & f_{i \kappa} \geq 0 & \forall i=1,2,\ldots,m, \kappa=1,2,\ldots,n \end{array} $$上の定式化において， 目的関数は最後の機械($m$)の最後の（$n$番目の）ジョブの完了時刻 を最小化している． 最初の制約は，各ジョブがいずれかの位置（順番）に割り当てられることを表す． 2番目の制約は，各位置（場所）にいずれかのジョブが割り当てられることを表す． 3番目の制約は，$\kappa$番目のジョブの完了時刻より $\kappa+1$番目のジョブの開始時刻が 遅いことを表す． 4番目の制約は，機械$i$ の $\kappa$番目のジョブの開始時刻と完了時刻の関係を 規定する制約であり， ジョブの位置を表す変数 $x$ と開始（終了）時刻を表す変数 $s$ ($f$) を繋ぐ制約である． これは，$x_{j\kappa}$ が $1$ のときに処理時間が $p_{ij}$ になる ことから導かれる． 5番目の制約は，機械$i$ 上で $\kappa$番目に割り当てられたジョブの完了時刻より， 機械$i+1$ 上で $\kappa$番目に割り当てられたジョブの開始時刻が遅いことを表す．
以下では，ジョブ数 $15$，機械数 $10$ の例題を求解する．
def permutation_flow_shop(n, m, p):
""" permutation_flow_shop problem
Parameters:
- n: number of jobs
- m: number of machines
- p[i,j]: processing time of job i on machine j
Returns a model, ready to be solved.
"""
model = Model("permutation flow shop")
x, s, f = {}, {}, {}
for j in range(1, n + 1):
for k in range(1, n + 1):
x[j, k] = model.addVar(vtype="B", name="x(%s,%s)" % (j, k))
for i in range(1, m + 1):
for k in range(1, n + 1):
s[i, k] = model.addVar(vtype="C", name="start(%s,%s)" % (i, k))
f[i, k] = model.addVar(vtype="C", name="finish(%s,%s)" % (i, k))
model.update()
for j in range(1, n + 1):
model.addConstr(
quicksum(x[j, k] for k in range(1, n + 1)) == 1, "Assign1(%s)" % (j)
)
model.addConstr(
quicksum(x[k, j] for k in range(1, n + 1)) == 1, "Assign2(%s)" % (j)
)
for i in range(1, m + 1):
for k in range(1, n + 1):
if k != n:
model.addConstr(f[i, k] <= s[i, k + 1], "FinishStart(%s,%s)" % (i, k))
if i != m:
model.addConstr(f[i, k] <= s[i + 1, k], "Machine(%s,%s)" % (i, k))
model.addConstr(
s[i, k] + quicksum(p[i, j] * x[j, k] for j in range(1, n + 1))
<= f[i, k],
"StartFinish(%s,%s)" % (i, k),
)
model.setObjective(f[m, n], GRB.MINIMIZE)
model.update()
model.__data = x, s, f
return model
def make_data_permutation_flow_shop(n, m):
"""make_data: prepare matrix of m times n random processing times"""
p = {}
for i in range(1, m + 1):
for j in range(1, n + 1):
p[i, j] = random.randint(1, 10)
return p
n = 15
m = 10
p = make_data_permutation_flow_shop(n, m)
model = permutation_flow_shop(n, m, p)
model.optimize()
x, s, f = model.__data
print("Opt.value=", model.ObjVal)
ジョブショップスケジューリング問題
ジョブショップスケジューリング問題 $J| | C_{\max}$ は，以下のように定義される古典的なスケジューリング問題である．
ジョブを $J_1, J_2,\cdots, J_n$， ジョブ $J_j$ に属するオペレーションを $O_{1j},O_{2j},\cdots,O_{m_j j}$， 機械を $M_1,M_2, \cdots, M_m$ とする． オペレーションは $O_{1j},O_{2j},\cdots,O_{m_j j}$ の順で処理されなければならず， オペレーション $O_{ij}$ を処理するには機械 $\mu_{ij}$ を作業時間 $p_{ij}$ だけ占有する． オペレーションが途中で中断できないという仮定の下で， 最後のオペレーションが完了する時刻を最小化する各機械上でのオペレーションの処理順序を求める問題．
この問題は 中規模な問題例でさえ数理最適化ソルバーで解くことは難しい． 実務的な付加条件のほとんどを考慮したスケジューリング最適化システムとしてOptSeqが開発されている． OptSeqについての詳細は，付録1を参照されたい．
ジョブショップスケジューリングのベンチマーク問題例も，OptSeqを使用すると比較的容易に良い解を得ることができる．
例として，以下のサイトから入手できるベンチマーク問題例を読み込んで解いてみよう．
- OR-Lib.: http://people.brunel.ac.uk/~mastjjb/jeb/orlib/jobshopinfo.html
- 多くのジョブショップスケジューリング問題例を含んだサイト: http://jobshop.jjvh.nl/
例としてft06.txtを用いる．このデータは，以下のテキストファイルである．
6 6
2 1 0 3 1 6 3 7 5 3 4 6
1 8 2 5 4 10 5 10 0 10 3 4
2 5 3 4 5 8 0 9 1 1 4 7
1 5 0 5 2 5 3 3 4 8 5 9
2 9 1 3 4 5 5 4 0 3 3 1
1 3 3 3 5 9 0 10 4 4 2 1得られた最大完了時刻 $55$ は最適値である．
from optseq import *
# ジョブショップスケジューリング問題のベンチマーク問題例
fname = "../data/jssp/ft06.txt"
#fname = "../data/jssp/ta01.txt"
f = open(fname, "r")
lines = f.readlines()
f.close()
n, m = map(int, lines[0].split())
print("n,m=", n, m)
model = Model()
act, mode, res = {}, {}, {}
for j in range(m):
res[j] = model.addResource(f"machine[{j}]", capacity=1)
# prepare data as dic
machine, proc_time = {}, {}
for i in range(n):
L = list(map(int, lines[i + 1].split()))
for j in range(m):
machine[i, j] = L[2 * j]
proc_time[i, j] = L[2 * j + 1]
for i, j in proc_time:
act[i, j] = model.addActivity(f"Act[{i}_{j}]")
mode[i, j] = Mode(f"Mode[{i}_{j}]", proc_time[i, j])
mode[i, j].addResource(res[machine[i, j]], 1)
act[i, j].addModes(mode[i, j])
for i in range(n):
for j in range(m - 1):
model.addTemporal(act[i, j], act[i, j + 1])
model.Params.TimeLimit = 1
model.Params.Makespan = True
#model.Params.OutputFlag = True
model.optimize()
import plotly.express as px
import pandas as pd
import datetime as dt
start="2022/1/1"
period="days"
start = pd.to_datetime(start)
def time_convert_long(periods):
if period =="days":
time_ = start + dt.timedelta(days=float(periods))
elif period == "hours":
time_ = start + dt.timedelta(hours=float(periods))
elif period == "minutes":
time_ = start + dt.timedelta(minutes=float(periods))
elif period == "seconds":
time_ = start + dt.timedelta(seconds=float(periods))
else:
raise TypeError("pariod must be 'days' or 'seconds' or minutes' or 'days'")
return time_.strftime("%Y-%m-%d")
L=[]
for i, j in proc_time:
a = act[i,j]
res = None
for r in a.modes[0].requirement:
res = r[0]
break
st = time_convert_long(a.start)
fi = time_convert_long(a.completion)
L.append(dict(Task=a.name, Start=st, Finish=fi, Resource=res, Job= f"Job{i}" ) )
df = pd.DataFrame(L)
fig = px.timeline(df, x_start="Start", x_end="Finish", y="Resource", color="Job", opacity=0.5)
plotly.offline.plot(fig);
fname = "../data/jssp/ft06.txt"
#fname = "../data/jssp/ta01.txt"
#fname = "../data/jssp/ta21.txt"
f = open(fname, "r")
lines = f.readlines()
f.close()
n, m = map(int, lines[0].split())
print("n,m=", n, m)
# prepare data
machine, proc_time = {}, {}
for i in range(n):
L = list(map(int, lines[i + 1].split()))
for j in range(m):
machine[i, j] = L[2 * j]
proc_time[i, j] = L[2 * j + 1]
jobs_data = []
for i in range(n):
row = []
for j in range(m):
row.append((machine[i, j], proc_time[i, j]))
jobs_data.append(row)
import collections
from ortools.sat.python import cp_model
model = cp_model.CpModel()
machines_count = 1 + max(task[0] for job in jobs_data for task in job)
all_machines = range(machines_count)
# Computes horizon dynamically as the sum of all durations.
horizon = sum(task[1] for job in jobs_data for task in job)
# Named tuple to store information about created variables
task_type = collections.namedtuple("task_type", "start end interval")
# Named tuple to manipulate solution information
assigned_task_type = collections.namedtuple(
"assigned_task_type", "start job index duration"
)
# Creates job intervals and add to the corresponding machine lists
all_tasks = {}
machine_to_intervals = collections.defaultdict(list)
for job_id, job in enumerate(jobs_data):
for task_id, task in enumerate(job):
machine = task[0]
duration = task[1]
suffix = "_%i_%i" % (job_id, task_id)
start_var = model.NewIntVar(0, horizon, "start" + suffix)
end_var = model.NewIntVar(0, horizon, "end" + suffix)
interval_var = model.NewIntervalVar(
start_var, duration, end_var, "interval" + suffix
)
all_tasks[job_id, task_id] = task_type(
start=start_var, end=end_var, interval=interval_var
)
machine_to_intervals[machine].append(interval_var)
# Create and add disjunctive constraints
for machine in all_machines:
model.AddNoOverlap(machine_to_intervals[machine])
# Precedences inside a job
for job_id, job in enumerate(jobs_data):
for task_id in range(len(job) - 1):
model.Add(
all_tasks[job_id, task_id + 1].start >= all_tasks[job_id, task_id].end
)
# Makespan objective
obj_var = model.NewIntVar(0, horizon, "makespan")
model.AddMaxEquality(
obj_var,
[all_tasks[job_id, len(job) - 1].end for job_id, job in enumerate(jobs_data)],
)
model.Minimize(obj_var)
# Solve model
solver = cp_model.CpSolver()
solver.parameters.max_time_in_seconds = 360.0
status = solver.Solve(model)
# Create one list of assigned tasks per machine.
assigned_jobs = collections.defaultdict(list)
for job_id, job in enumerate(jobs_data):
for task_id, task in enumerate(job):
machine = task[0]
assigned_jobs[machine].append(
assigned_task_type(
start=solver.Value(all_tasks[job_id, task_id].start),
job=job_id,
index=task_id,
duration=task[1],
)
)
# Create per machine output lines.
output = ""
for machine in all_machines:
# Sort by starting time.
assigned_jobs[machine].sort()
sol_line_tasks = "Machine " + str(machine) + ": "
sol_line = " "
for assigned_task in assigned_jobs[machine]:
name = "job_%i_%i" % (assigned_task.job, assigned_task.index)
# Add spaces to output to align columns.
sol_line_tasks += "%-10s" % name
start = assigned_task.start
duration = assigned_task.duration
sol_tmp = "[%i,%i]" % (start, start + duration)
# Add spaces to output to align columns.
sol_line += "%-10s" % sol_tmp
sol_line += "\n"
sol_line_tasks += "\n"
output += sol_line_tasks
output += sol_line
# Finally print the solution found.
print("Optimal Schedule Length: %i" % solver.ObjectiveValue())
print(output)
数理最適化ソルバーを用いた求解
数理最適化による定式化も可能である． 以下では，次の資源制約付きプロジェクトスケジューリング問題で導入する時刻添え字定式化と， 1機械リリース時刻付き完了時刻和最小化問題で用いた離接定式化による求解を試みる．
ここでは，各ジョブに含まれるオペレーション数が機械の数 $m$ に一致するものと仮定する． オペレーションをジョブと順序の組として $(i,j)$ もしくは $(k,l)$ と書く．
ジョブショップスケジューリング問題は資源制約付きプロジェクトスケジューリング問題の特殊形であるので，時刻添え字定式化は省略する．
離接定式化では， オペレーション $(i,j)$ の開始時刻を表す連続変数 $s_{ij}$ と， オペレーション $(i,j)$ が他のオペレーション $(k,l)$ に先行する（前に処理される）とき $1$ になる $0$-$1$ 整数変数 $x_{ijkl}$ を用いる． また，最大完了時刻を表す実数変数を $z$ とする．
非常に大きな数 $M$ を用いると，離散定式化は以下のように書ける．
$$ \begin{array}{l l l} minimize & z & \\ s.t. & s_{ij} + p_{ij} -M (1-x_{ijkl}) \leq s_{kl} & \forall j \neq k \\ & x_{ijkl} +x _{klij}= 1 & \forall (i,j) < (k,l) \\ & s_{ij} + p_{ij} \leq s_{i,j+1} & \forall i, j=1,\ldots,m-1 \\ & s_{im} \leq z & \forall i \\ & s_{i1} \geq 0 & \forall i \\ & x_{ijkl} \in \{ 0,1 \} & \forall (i,j) \neq (k,l) \end{array} $$fname = "../data/jssp/ft06.txt"
# fname = "../data/jssp/ta21.txt"
f = open(fname, "r")
lines = f.readlines()
f.close()
n, m = map(int, lines[0].split())
print("n,m=", n, m)
machine, proc_time = {}, {}
for i in range(n):
L = list(map(int, lines[i + 1].split()))
for j in range(m):
machine[i, j] = L[2 * j]
proc_time[i, j] = L[2 * j + 1]
from gurobipy import Model, GRB, quicksum
model = Model("scheduling: time-index")
T = 55 # メイクスパン以上の計画期間
s, x = {}, {} # s - start time variable, x=1 if job (i,j) starts at period t
for i, j in proc_time:
s[i, j] = model.addVar(vtype="C", name=f"s({i},{j})")
for t in range(1, T - proc_time[i, j] + 2):
x[i, j, t] = model.addVar(vtype="B", name=f"x({i},{j},{t})")
z = model.addVar(vtype="C", name="z") # completion time
model.update()
for i, j in proc_time:
# job execution constraints
model.addConstr(
quicksum(x[i, j, t] for t in range(1, T - proc_time[i, j] + 2)) == 1
)
# start time constraints
model.addConstr(
quicksum((t - 1) * x[i, j, t] for t in range(2, T - proc_time[i, j] + 2))
== s[i, j]
)
# resource upper bound constraints
for t in range(1, T + 1):
for k in range(m): # machine
model.addConstr(
quicksum(
x[i, j, t_]
for i, j in proc_time
if machine[i, j] == k
for t_ in range(
max(t - proc_time[i, j] + 1, 1), min(t + 1, T - proc_time[i, j] + 2)
)
)
<= 1
)
# time (precedence) constraints, i.e., s[k]-s[j] >= p[j]
for i in range(n):
for j in range(m - 1):
model.addConstr(s[i, j] + proc_time[i, j] <= s[i, j + 1])
model.addConstr(s[i, m - 1] + proc_time[i, m - 1] <= z)
model.setObjective(z, GRB.MINIMIZE)
model.optimize()
from gurobipy import Model, GRB, quicksum
model = Model("scheduling: disjunctive")
M = 3000 # big M
s, x = {}, {} # start time variable, x[i,j,k,l] = 1 if job (i,j) precedes job (k,l), 0 otherwise
for i,j in proc_time:
s[i,j] = model.addVar(vtype="C", name= f"s({i},{j})")
for k,l in proc_time:
if machine[i,j]==machine[k,l] and (i,j) != (k,l):
x[i, j, k, l] = model.addVar(vtype="B", ub=1., name= f"x({i},{j},{k},{l})")
z = model.addVar(vtype="C", name= "z" ) #completion time
model.update()
for i,j in proc_time:
for k,l in proc_time:
if machine[i,j]==machine[k,l] and (i,j) != (k,l):
model.addConstr(
s[i,j] - s[k,l] + M * x[i, j, k, l] <= M - proc_time[i,j]
)
if machine[i,j]==machine[k,l] and (i,j) < (k,l):
model.addConstr(x[i,j,k,l] + x[k,l,i,j] == 1)
for i in range(n):
for j in range(m-1):
model.addConstr(s[i,j] +proc_time[i,j] <= s[i,j+1])
model.addConstr(s[i,m-1]+ proc_time[i,m-1] <=z )
model.setObjective(z, GRB.MINIMIZE)
model.Params.LogFile = "gurobi_jssp.log"
model.Params.TimeLimit = 60
model.Params.OutputFlag = 1
model.optimize()
#パラメータチューニング
# model.Params.TuneTimeLimit = 3600
# model.Params.TuneCriterion = 2 # best feasible solution found.
# model.tune()
import grblogtools as glt
results = glt.parse(model.Params.LogFile)
summary = results.summary()
nodelog_progress = results.progress("nodelog")
nodelog_progress.Incumbent.plot()
nodelog_progress.BestBd.plot();
#summary.T
資源制約付きプロジェクトスケジューリング問題
様々なスケジューリング問題を一般化した，以下の資源制約付き（プロジェクト）スケジューリング問題を考える．
ジョブの集合 $Job$ ， 各時間ごとの使用可能量の上限をもつ資源 $Res$ が与えられている． 資源は，ジョブごとに決められた作業時間の間はジョブによって使用されるが， 作業完了後は，再び使用可能になる． ジョブ間に与えられた時間制約を満たした上で， ジョブの作業開始時刻に対する任意の関数の和を最小化するような， 資源のジョブへの割り振りおよびジョブの作業開始時刻を求める．
時刻を離散化した期の概念を用いた定式化（時刻添え字定式化）を示す．
まず，定式化で用いる集合，入力データ，変数を示す．
集合:
- $Job$: ジョブの集合．添え字は $j,k$
- $Res$: 資源の集合．添え字は $r$
- $Prec$: ジョブ間の時間制約を表す集合． $Job \times Job$ の部分集合で， $(j,k) \in Prec$ のとき， ジョブ $j$ とジョブ $k$ の開始（完了）時刻間に何らかの制約が定義されるものとする
入力データ:
- $T$: 最大の期数．期の添え字は $t=1,2,\ldots,T$． 連続時間との対応づけは，時刻 $t-1$ から $t$ までを期 $t$ と定義する． なお，期を表す添え字は $t$ の他に $s$ も用いる．
- $p_j$: ジョブ $j$ の処理時間．非負の整数を仮定
- $Cost_{jt}$: ジョブ $j$ を期$t$ に開始したときの費用
- $a_{jrt}$: ジョブ $j$ の開始後 $t(=0,1,\ldots,p_j-1)$期経過時の処理に要する資源$r$ の量
- $RUB_{rt}$: 時刻 $t$ における資源 $r$ の使用可能量の上限
変数:
- $x_{jt}$: ジョブ $j$ を時刻 $t$ に開始するとき $1$，それ以外のとき $0$ を表す $0$-$1$変数
以下に資源制約付きスケジューリング問題の定式化を示す．
目的関数: $$ minimize \sum_{j \in Job} \sum_{t=1}^{T-p_j+1} Cost_{jt} x_{jt} $$
ジョブ遂行制約： $$ \sum_{t=1}^{T-p_j+1} x_{jt} =1 \forall j \in Job $$ すべてのジョブは必ず1度処理されなければならないことを表す．
- 資源制約：
ある時刻 $t$ に作業中であるジョブの資源使用量の合計 が，資源使用量の上限を超えないことを規定する． 時刻 $t$ に作業中であるジョブは， 時刻 $t-p_j+1$ から $t$ の間に作業を開始したジョブであるので，上式を得る．
- 時間制約：
ジョブ $j$ の開始時刻 $S_j$ は $\sum_{t=2}^{T-p_j+1} (t-1)x_{jt}$ であるので， 完了時刻 $C_j$ はそれに $p_j$ を加えたものになる． ジョブの対 $(j,k) \in Prec$ に対しては，ジョブ $j,k$ の開始（完了）時刻間に 何らかの制約が付加される．たとえば，ジョブ $j$ の処理が終了するまで，ジョブ $k$ の処理が開始できないことを表す先行制約は，
$$ \sum_{t=2}^{T-p_j+1} (t-1)x_{jt} +p_j \leq \sum_{t=2}^{T-p_k+1} (t-1) x_{kt} \ \ \ \forall (j,k) \in Prec $$と書くことができる． ここで，左辺の式は，ジョブ $j$の完了時刻を表し，右辺の式はジョブ $k$ の開始時刻を表す．
以下に，数理最適化ソルバー Gurobiを用いたモデルを示す． 実際には，Gurobiでは中規模問題例までしか解くことができない． ジョブショップスケジューリング問題と同様に，OptSeqやOR-toolsを用いることによって求解できる．
ただし，実務の問題はさらなる付加条件がつくことが多い．OptSeqは，より一般的なスケジューリング問題をモデル化できるように設計されている．
from gurobipy import Model, quicksum, GRB, multidict
# from mypulp import Model, quicksum, GRB, multidict
def rcs(J, P, R, T, p, c, a, RUB):
"""rcs -- model for the resource constrained scheduling problem
Parameters:
- J: set of jobs
- P: set of precedence constraints between jobs
- R: set of resources
- T: number of periods
- p[j]: processing time of job j
- c[j,t]: cost incurred when job j starts processing on period t.
- a[j,r,t]: resource r usage for job j on period t (after job starts)
- RUB[r,t]: upper bound for resource r on period t
Returns a model, ready to be solved.
"""
model = Model("resource constrained scheduling")
s, x = {}, {} # s - start time variable, x=1 if job j starts on period t
for j in J:
s[j] = model.addVar(vtype="C", name="s(%s)" % j)
for t in range(1, T - p[j] + 2):
x[j, t] = model.addVar(vtype="B", name="x(%s,%s)" % (j, t))
model.update()
for j in J:
# job execution constraints
model.addConstr(
quicksum(x[j, t] for t in range(1, T - p[j] + 2)) == 1,
"ConstrJob(%s,%s)" % (j, t),
)
# start time constraints
model.addConstr(
quicksum((t - 1) * x[j, t] for t in range(2, T - p[j] + 2)) == s[j],
"ConstrJob(%s,%s)" % (j, t),
)
# resource upper bound constraints
for t in range(1, T+1):
for r in R:
model.addConstr(
quicksum(
a[j, r, t - t_] * x[j, t_]
for j in J
for t_ in range(max(t - p[j] + 1, 1), min(t + 1, T - p[j] + 2))
)
<= RUB[r, t],
"ResourceUB(%s,%s)" % (r,t)
)
# time (precedence) constraints, i.e., s[k]-s[j] >= p[j]
for (j, k) in P:
model.addConstr(s[k] - s[j] >= p[j], "Precedence(%s,%s)" % (j, k))
model.setObjective(quicksum(c[j, t] * x[j, t] for (j, t) in x), GRB.MINIMIZE)
model.update()
model.__data = x, s
return model
def make_1r():
J, p = multidict(
{ # jobs, processing times
1: 1,
2: 3,
3: 2,
4: 2,
}
)
P = [(1, 2), (1, 3), (2, 4)]
R = [1]
T = 6
c = {}
for j in J:
for t in range(1, T - p[j] + 2):
c[j, t] = 1 * (t - 1 + p[j])
a = {
(1, 1, 0): 2,
(2, 1, 0): 2,
(2, 1, 1): 1,
(2, 1, 2): 1,
(3, 1, 0): 1,
(3, 1, 1): 1,
(4, 1, 0): 1,
(4, 1, 1): 2,
}
RUB = {(1, 1): 2, (1, 2): 2, (1, 3): 1, (1, 4): 2, (1, 5): 2, (1, 6): 2}
return (J, P, R, T, p, c, a, RUB)
def make_2r():
J, p = multidict(
{ # jobs, processing times
1: 2,
2: 2,
3: 3,
4: 2,
5: 5,
}
)
P = [(1, 2), (1, 3), (2, 4)]
R = [1, 2]
T = 6
c = {}
for j in J:
for t in range(1, T - p[j] + 2):
c[j, t] = 1 * (t - 1 + p[j])
a = {
# resource 1:
(1, 1, 0): 2,
(1, 1, 1): 2,
(2, 1, 0): 1,
(2, 1, 1): 1,
(3, 1, 0): 1,
(3, 1, 1): 1,
(3, 1, 2): 1,
(4, 1, 0): 1,
(4, 1, 1): 1,
(5, 1, 0): 0,
(5, 1, 1): 0,
(5, 1, 2): 1,
(5, 1, 3): 0,
(5, 1, 4): 0,
# resource 2:
(1, 2, 0): 1,
(1, 2, 1): 0,
(2, 2, 0): 1,
(2, 2, 1): 1,
(3, 2, 0): 0,
(3, 2, 1): 0,
(3, 2, 2): 0,
(4, 2, 0): 1,
(4, 2, 1): 2,
(5, 2, 0): 1,
(5, 2, 1): 2,
(5, 2, 2): 1,
(5, 2, 3): 1,
(5, 2, 4): 1,
}
RUB = {
(1, 1): 2,
(1, 2): 2,
(1, 3): 2,
(1, 4): 2,
(1, 5): 2,
(1, 6): 2,
(1, 7): 2,
(2, 1): 2,
(2, 2): 2,
(2, 3): 2,
(2, 4): 2,
(2, 5): 2,
(2, 6): 2,
(2, 7): 2,
}
return (J, P, R, T, p, c, a, RUB)
(J,P,R,T,p,c,a,RUB) = make_2r()
model = rcs(J,P,R,T,p,c,a,RUB)
model.optimize()
x,s = model.__data
print ("Opt.value=",model.ObjVal)
for (j,t) in x:
if x[j,t].X > 0.5:
print (x[j,t].VarName,x[j,t].X)
for j in s:
if s[j].X > 0.:
print (s[j].VarName,s[j].X )