動的ロットサイズ決定問題に対する定式化とアルゴリズム
import numpy as np
import string
import networkx as nx
import matplotlib.pyplot as plt
from collections import defaultdict
import random
from gurobipy import Model, quicksum, GRB, multidict
# from mypulp import Model, quicksum, GRB, multidict
単一段階・単一品目動的ロットサイズ決定問題
ここでは，複雑な実際問題の基礎となる単一段階・単一品目モデルを考える． 単一段階・単一品目の動的ロットサイズ決定問題(dynamic lotsizing problem)の基本形は，以下の仮定をもつ．
- 期によって変動する需要量をもつ単一の品目を扱う．
- 品目を生産する際には，生産数量に依存しない固定費用と数量に比例する変動費用がかかる．
- 計画期間はあらかじめ決められており，最初の期における在庫量（初期在庫量）は $0$ とする． この条件は，実際には任意の初期在庫量をもつように変形できるが，以下では議論を簡略化するため， 初期在庫量が $0$ であると仮定する．
- 次の期に持ち越した品目の量に比例して在庫（保管）費用がかかる．
- 生産時間は $0$ とする．これは，生産を行ったその期にすぐに需要を満たすことができることを表す． （生産ではなく）発注を行う場合には，発注すればすぐに商品が届くこと，言い換えればリード時間が $0$ であることに相当する．
- 各期の生産可能量には上限がある．
- 生産固定費用，生産変動費用，ならびに在庫費用の合計を最小にするような生産方策を決める．
問題を明確にするために，単一段階・単一品目のモデルの定式化を示す 以下に定式化に必要な記号を，パラメータ（定数）と変数に分けて記述する．
パラメータ
- $T$: 計画期間数．期を表す添え字を $1,2,\ldots,t,\ldots,T$ と記す
- $f_t$: 期 $t$ において生産を行うために必要な段取り（固定）費用
- $c_t$: 期 $t$ における品目 $1$個あたりの生産変動費用
- $h_t$: 期 $t$ における（品目 $1$ 個あたり，$1$ 期間あたりの）在庫費用
- $d_t$: 期 $t$ における品目の需要量
- $M_t$: 期 $t$ における生産可能量の上限．これを生産の容量とよぶこともある
変数
- $I_t$: 期 $t$ における在庫量．より正確に言うと，期 $t$ の期末の在庫量
- $x_t$: 期 $t$ における生産量
- $y_t$: 期 $t$ に生産を行うとき $1$，それ以外のとき $0$ を表す $0$-$1$ 変数
上の記号を用いると，単一段階・単一品目の動的ロットサイズ決定問題の基本形は，以下のように定式化できる．
$$ \begin{array}{ l l l } minimize & \sum_{t=1}^T \left( f_t y_t +c_t x_t + h_t I_t \right) & \\ s.t. & I_{t-1} +x_t -I_t =d_t & \forall t=1,\ldots,T \\ & x_t \leq M_t y_t & \forall t=1,\ldots,T \\ & I_0 =0 & \\ & x_t,I_t \geq 0 & \forall t=1,\ldots,T \\ & y_t \in \{0,1\} & \forall t=1,\ldots,T \end{array} $$上の定式化で，最初の制約式は、各期における品目の在庫保存式であり， 前期からの繰り越しの在庫量 $I_{t-1}$ に今期の生産量 $x_t$ を加え， 需要量 $d_t$ を減じたものが，来期に持ち越す在庫量 $I_t$ であることを意味する． 2番目の制約式は，生産を行わない期における生産量が $0$ であり， 生産を行う期においては，その上限が $M_t$ 以下であることを保証するための式である． 3番目の制約式は，初期在庫量が $0$ であることを表す．
各期の生産可能量に制約がない場合には，容量制約なしのロットサイズ決定モデルもしくは発案者の名前をとってWagner--Whitinモデルとよばれる．
動的最適化による求解
容量制約なしの単一段階・単一品目のロットサイズ決定モデルは， 動的最適化を適用することによって，効率的に解くことができる．
$1,\ldots,j$ 期までの需要を満たすときの最小費用を $F(j)$ と書く． 初期（境界）条件は，仮想の期 $0$ を導入することによって， $$ F(0)= 0 $$ と書ける． 再帰方程式は，期 $j$ までの最小費用が， 期 $i (<j)$ までの最小費用に，期 $i+1$ に生産を行うことによって期 $i+1$ から $j$ までの需要をまかなう ときの費用 $f_{i+1} + c_{i+1} \left (\sum_{t=i+1}^j d_t \right) + \sum_{s=i+1}^{j-1} h_s \sum_{t=s+1}^j d_t$ を加えたものになることから， $$ F(j) =\min_{i \in \{1,\ldots,j-1 \} } \left\{ F(i)+ f_{i+1} + c_{i+1} \left (\sum_{t=i+1}^j d_t \right) + \sum_{s=i+1}^{j-1} h_s \sum_{t=s+1}^j d_t \right\} $$ と書ける．上の再帰方程式を $j=1,2,\ldots,T$ の順に計算することによって， もとの問題の最適費用 $F(T)$ を得ることができる．
以下に示すのは，高速化を施した動的最適化であり， ほとんどの問題例に対して $O(T)$ 時間で終了する．
# based on Evans' forward fast implementation
T = 5
fixed = [3, 3, 3, 3, 3]
variable = [1, 1, 3, 3, 3]
demand = [5, 7, 3, 6, 4]
h = [1 for i in range(T)]
F = [9999999 for i in range(T)]
for i in range(T):
if i == 0:
cum = fixed[i] + variable[i] * demand[i]
else:
cum = F[i - 1] + fixed[i] + variable[i] * demand[i]
cumh = 0
for j in range(i, T):
if cum < F[j]:
F[j] = cum
if j == (T - 1):
break
cumh += h[j]
cum += (variable[i] + cumh) * demand[j + 1]
if (
fixed[j + 1] + variable[j + 1] * demand[j + 1]
< (variable[i] + cumh) * demand[j + 1]
):
break
print("Optimal Value=", F[T - 1])
単段階多品目動的ロットサイズ決定問題
ここでは，施設配置定式化と呼ばれる強化された定式化を示す．
我々の想定しているロットサイズ決定問題では品切れは許さないので， 期$t$ の需要は，期$t$ もしくはそれ以前の期で生産された量でまかなわれなければならない． そこで各品目$p$ に対して，期$t$ の需要のうち，期$s (s \leq t)$ の生産によってまかなわれた量を表す変数 $X_{st}^p$ を導入する．
需要 $1$単位あたりの変動費用は，生産変動費用 $c_s^p$ と 在庫費用 $\sum_{\ell=s}^{t-1} h_{\ell}^p$ の和になるので， $X_{st}^p$ の係数 $C_{st}^p$ は，以下のように定義される．
$$ C_{st}^p= c_s^p +\sum_{\ell=s}^{t-1} h_{\ell}^p \ \ \ \forall p \in P $$すべての需要は満たされなければならず，また $X_{st}^p >0$ のときには期$s$ で生産しなければならないので，$y_s^p=1$ となる． ここで，$y_t^p$ は前節の標準定式化における段取りの有無を表す $0$-$1$ 変数で， 品目$p$ を期$t$ に生産をするときに $1$，それ以外のとき $0$ を表すことを思い起こされたい． よって，以下の定式化を得る．
$$ \begin{array}{ l l l } minimize & \sum_{s t: s \leq t} \sum_{p \in P} C_{st}^p X_{st}^p + \sum_{t=1}^T \sum_{p \in P} f_t^p y_t^p & \\ s.t. & \sum_{s=1}^t X_{st}^p=d_t^p & \forall p \in P, t=1,2,\ldots,T \\ & \sum_{p \in P} \sum_{j: j \geq t} X_{tj}^p +\sum_{p \in P} g_t^p y_t^p \leq M_t & \forall t=1,2,\ldots,T \\ & X_{st}^p \leq d_t^p y_s^p & \forall p \in P, s=1,2,\ldots,t, t=1,2,\ldots,T \\ & X_{st}^p \geq 0 & \forall p \in P, s=1,2,\ldots,t, t=1,2,\ldots,T \\ & y_t^p \in \{0,1\} & \forall p \in P, t=1,2,\ldots,T \end{array} $$最初の制約は，各期の品目の需要が満たされることを規定する式であり， 2番目の制約は，各期の生産時間の上限を表す． 3番目の式は，品目が段取りを行わない期には生産ができないことを表す．
以下にベンチマーク用のデータ生成のコードと，数理最適化ソルバーによる最適化のためのコードを示す．
問題例の難しさは期数 $T$，品目数 $N$ と引数factorで制御される．factorが $0.75$ のとき緩い制約の簡単な問題例になり，$1.1$ のとききつい制約の難しい問題例となる．
def trigeiro(T, N, factor):
"""
Data generator for the multi-item lot-sizing problem
it uses a simular algorithm for generating the standard benchmarks in:
"Capacitated Lot Sizing with Setup Times" by
William W. Trigeiro, L. Joseph Thomas, John O. McClain
MANAGEMENT SCIENCE
Vol. 35, No. 3, March 1989, pp. 353-366
Parameters:
- T: number of periods
- N: number of products
- factor: value for controlling constraining factor of capacity:
- 0.75: lightly-constrained instances
- 1.10: constrained instances
"""
random.seed(123)
P = range(1, N + 1)
f, g, c, d, h, M = {}, {}, {}, {}, {}, {}
sumT = 0
for t in range(1, T + 1):
for p in P:
# setup times
g[t, p] = 10 * random.randint(1, 5) # 10, 50: trigeiro's values
# set-up costs
f[t, p] = 100 * random.randint(1, 10) # checked from Wolsey's instances
c[t, p] = 0 # variable costs
# demands
d[t, p] = 100 + random.randint(-25, 25) # checked from Wolsey's instances
if t <= 4:
if random.random() < 0.25: # trigeiro's parameter
d[t, p] = 0
sumT += (
g[t, p] + d[t, p]
) # sumT is the total capacity usage in the lot-for-lot solution
h[t, p] = random.randint(
1, 5
) # holding costs; checked from Wolsey's instances
for t in range(1, T + 1):
M[t] = int(float(sumT) / float(T) / factor)
return P, f, g, c, d, h, M
def mils_fl(T, P, f, g, c, d, h, M):
"""
mils_fl: facility location formulation for the multi-item lot-sizing problem
Requires more variables, but gives a better solution because LB is
better than the standard formulation. It can be used as a
heuristic method that is sometimes better than relax-and-fix.
Parameters:
- T: number of periods
- P: set of products
- f[t,p]: set-up costs (on period t, for product p)
- g[t,p]: set-up times
- c[t,p]: variable costs
- d[t,p]: demand values
- h[t,p]: holding costs
- M[t]: resource upper bound on period t
Returns a model, ready to be solved.
"""
Ts = range(1, T + 1)
model = Model("multi-item lotsizing -- facility location formulation")
y, X = {}, {}
for p in P:
for t in Ts:
y[t, p] = model.addVar(vtype="B", name="y(%s,%s)" % (t, p))
for s in range(1, t + 1):
X[s, t, p] = model.addVar(name="X(%s,%s,%s)" % (s, t, p))
model.update()
for t in Ts:
# capacity constraints
model.addConstr(
quicksum(X[t, s, p] for s in range(t, T + 1) for p in P)
+ quicksum(g[t, p] * y[t, p] for p in P)
<= M[t],
"Capacity(%s)" % (t),
)
for p in P:
# demand satisfaction constraints
model.addConstr(
quicksum(X[s, t, p] for s in range(1, t + 1)) == d[t, p],
"Demand(%s,%s)" % (t, p),
)
# connection constraints
for s in range(1, t + 1):
model.addConstr(
X[s, t, p] <= d[t, p] * y[s, p], "Connect(%s,%s,%s)" % (s, t, p)
)
C = {} # variable costs plus holding costs
for p in P:
for s in Ts:
sumC = 0
for t in range(s, T + 1):
C[s, t, p] = c[s, p] + sumC
sumC += h[t, p]
model.setObjective(
quicksum(f[t, p] * y[t, p] for t in Ts for p in P)
+ quicksum(
C[s, t, p] * X[s, t, p] for t in Ts for p in P for s in range(1, t + 1)
),
GRB.MINIMIZE,
)
model.update()
model.__data = y, X
return model
T, N, factor = 15, 16, 1.1 #0.75
P, f, g, c, d, h, M = trigeiro(T, N, factor)
model = mils_fl(T, P, f, g, c, d, h, M)
model.optimize()
print("Opt.value=", model.ObjVal)
y, X = model.__data
標準定式化
集合:
- $\{1..T\}$: 期間の集合
- $P$ : 品目の集合（完成品と部品や原材料を合わせたものを「品目」と定義する）
- $K$ : 生産を行うのに必要な資源（機械，生産ライン，工程などを表す）の集合
- $P_k$ : 資源 $k$ で生産される品目の集合
- $Parent_p$ : 部品展開表における品目（部品または材料）$p$ の親品目の集合．言い換えれば，品目 $p$ から製造される品目の集合
パラメータ:
- $T$: 計画期間数．期を表す添え字を $1,2,\ldots,t,\ldots,T$ と記す
- $f_t^p$ : 期 $t$ に品目 $p$ に対する段取り替え（生産準備）を行うときの費用（段取り費用）
- $g_t^p$ : 期 $t$ に品目 $p$ に対する段取り替え（生産準備）を行うときの時間（段取り時間）
- $c_t^p$ : 期 $t$ における品目 $p$ の生産変動費用
- $h_t^p$ : 期 $t$ から期 $t+1$ に品目 $p$ を持ち越すときの単位あたりの在庫費用
- $d_t^p$ : 期 $t$ における品目 $p$ の需要量
- $\phi_{pq}$ : $q \in Parent_p$ のとき， 品目 $q$ を $1$ 単位製造するのに必要な品目 $p$ の数 （$p$-units）． ここで， $p$-unitsとは，品目 $q$ の $1$単位と混同しないために導入された単位であり， 品目 $p$ の $1$単位を表す．$\phi_{pq}$ は，部品展開表を有向グラフ表現したときには，枝の重みを表す
- $M_t^k$ : 期 $t$ における資源 $k$ の使用可能な生産時間の上限． 定式化では，品目 $1$単位の生産時間を $1$単位時間になるようにスケーリングしてあるものと仮定しているが， プログラム内では単位生産量あたりの生産時間を定義している
- $UB_t^p$ : 期 $t$ における品目 $p$ の生産時間の上限 品目 $p$ を生産する資源が $k$ のとき，資源の使用可能時間の上限 $M_t^k$ から段取り替え時間 $g_t^p$ を減じたものと定義される
変数:
- $x_t^p$（x）: 期 $t$ における品目 $p$ の生産量
- $I_t^p$（inv） : 期 $t$ における品目 $p$ の在庫量
- $y_t^p$（y）: 期 $t$ に品目 $p$ に対する段取りを行うとき $1$， それ以外のとき $0$ を表す $0$-$1$ 変数
上の記号を用いると、多段階ロットサイズ決定モデルは，以下のように定式化できる．
$$ \begin{array}{ l l l } minimize & \sum_{t=1}^T \sum_{p \in P} \left( f_t^p y_t^p + c_t^p x_t^p + h_t^p I_t^p \right) & \\ s.t. & \ \ I_{t-1}^p +x_t^p = d_t^p+ \sum_{q \in Parent_p} \phi_{pq} x_t^q +I_t^p & \forall p \in P, t=1,\ldots,T \\ & \sum_{p \in P_k} x_t^p +\sum_{p \in P_k} g_t^p y_t^p \leq M_t^k & \forall k \in K, t=1,\ldots,T \\ & x_t^p \leq UB_t^p y_t^p & \forall p \in P, t=1,\ldots,T \\ & I_0^p =0 & \forall p \in P \\ & x_t^p,I_t^p \geq 0 & \forall p \in P, t=1,\ldots,T \\ & y_t \in \{0,1\} & \forall t=1,\ldots,T \end{array} $$上の定式化で，最初の制約式は，各期および各品目に対する在庫の保存式を表す． より具体的には，品目 $p$ の期 $t-1$ からの在庫量 $I_{t-1}^p$ と生産量 $x_t^p$ を加えたものが， 期 $t$ における需要量 $d_t^p$，次期への在庫量 $I_t^p$， および他の品目を生産するときに必要な量 $\sum_{q \in Parent_p} \phi_{pq} x_t^q$ の和に等しいことを表す．
2番目の制約は， 各期の生産時間の上限制約を表す． 定式化ではすべての品目の生産時間は， $1$ 単位あたり$1$ 時間になるようにスケーリングしてあると仮定していたが，実際問題のモデル化の際には， 品目 p を $1$ 単位生産されるときに，資源 $r$ を使用する時間を用いた方が汎用性がある．
3番目の式は，段取り替えをしない期は生産できないことを表す．
def mils_standard(T, K, P, f, g, c, d, h, a, M, UB, phi):
"""
mils_standard: standard formulation for the multi-item, multi-stage lot-sizing problem
Parameters:
- T: number of periods
- K: set of resources
- P: set of items
- f[t,p]: set-up costs (on period t, for product p)
- g[t,p]: set-up times
- c[t,p]: variable costs
- d[t,p]: demand values
- h[t,p]: holding costs
- a[t,k,p]: amount of resource k for producing p in period t
- M[t,k]: resource k upper bound on period t
- UB[t,p]: upper bound of production time of product p in period t
- phi[(i,j)]: units of i required to produce a unit of j (j parent of i)
"""
model = Model("multi-stage lotsizing -- standard formulation")
y, x, I = {}, {}, {}
Ts = range(1, T + 1)
for p in P:
for t in Ts:
y[t, p] = model.addVar(vtype="B", name="y(%s,%s)" % (t, p))
x[t, p] = model.addVar(vtype="C", name="x(%s,%s)" % (t, p))
I[t, p] = model.addVar(vtype="C", name="I(%s,%s)" % (t, p))
I[0, p] = model.addVar(name="I(%s,%s)" % (0, p))
model.update()
for t in Ts:
for p in P:
# flow conservation constraints
model.addConstr(
I[t - 1, p] + x[t, p]
== quicksum(phi[p, q] * x[t, q] for (p2, q) in phi if p2 == p)
+ I[t, p]
+ d[t, p],
"FlowCons(%s,%s)" % (t, p),
)
# capacity connection constraints
model.addConstr(x[t, p] <= UB[t, p] * y[t, p], "ConstrUB(%s,%s)" % (t, p))
# time capacity constraints
for k in K:
model.addConstr(
quicksum(a[t, k, p] * x[t, p] + g[t, p] * y[t, p] for p in P)
<= M[t, k],
"TimeUB(%s,%s)" % (t, k),
)
# initial inventory quantities
for p in P:
model.addConstr(I[0, p] == 0, "InventInit(%s)" % (p))
model.setObjective(
quicksum(
f[t, p] * y[t, p] + c[t, p] * x[t, p] + h[t, p] * I[t, p]
for t in Ts
for p in P
),
GRB.MINIMIZE,
)
model.update()
model.__data = y, x, I
return model
def calc_rho(phi):
D = SCMGraph()
D.add_weighted_edges_from([(i, j, phi[i, j]) for (i, j) in phi])
ancestors = D.find_ancestors()
rho = defaultdict(float)
for v in D.up_order():
for w in D.succ[v]:
for q in ancestors[v]:
if q in ancestors[w]:
if w == q:
rho[v, q] = phi[v, w]
else:
rho[v, q] += rho[w, q] * phi[v, w]
return rho
def mils_echelon(T, K, P, f, g, c, d, h, a, M, UB, phi):
"""
mils_echelon: echelon formulation for the multi-item, multi-stage lot-sizing problem
Parameters:
- T: number of periods
- K: set of resources
- P: set of items
- f[t,p]: set-up costs (on period t, for product p)
- g[t,p]: set-up times
- c[t,p]: variable costs
- d[t,p]: demand values
- h[t,p]: holding costs
- a[t,k,p]: amount of resource k for producing p in period t
- M[t,k]: resource k upper bound on period t
- UB[t,p]: upper bound of production time of product p in period t
- phi[(i,j)]: units of i required to produce a unit of j (j parent of i)
"""
rho = calc_rho(
phi
) # rho[(i,j)]: units of i required to produce a unit of j (j ancestor of i)
model = Model("multi-stage lotsizing -- echelon formulation")
y, x, E, H = {}, {}, {}, {}
Ts = range(1, T + 1)
for p in P:
for t in Ts:
y[t, p] = model.addVar(vtype="B", name="y(%s,%s)" % (t, p))
x[t, p] = model.addVar(vtype="C", name="x(%s,%s)" % (t, p))
H[t, p] = h[t, p] - sum([h[t, q] * phi[q, p] for (q, p2) in phi if p2 == p])
E[t, p] = model.addVar(
vtype="C", name="E(%s,%s)" % (t, p)
) # echelon inventory
E[0, p] = model.addVar(vtype="C", name="E(%s,%s)" % (0, p)) # echelon inventory
model.update()
for t in Ts:
for p in P:
# flow conservation constraints
dsum = d[t, p] + sum([rho[p, q] * d[t, q] for (p2, q) in rho if p2 == p])
model.addConstr(
E[t - 1, p] + x[t, p] == E[t, p] + dsum, "FlowCons(%s,%s)" % (t, p)
)
# capacity connection constraints
model.addConstr(x[t, p] <= UB[t, p] * y[t, p], "ConstrUB(%s,%s)" % (t, p))
# time capacity constraints
for k in K:
model.addConstr(
quicksum(a[t, k, p] * x[t, p] + g[t, p] * y[t, p] for p in P)
<= M[t, k],
"TimeUB(%s,%s)" % (t, k),
)
# calculate echelon quantities
for p in P:
model.addConstr(E[0, p] == 0, "EchelonInit(%s)" % (p))
for t in Ts:
model.addConstr(
E[t, p] >= quicksum(phi[p, q] * E[t, q] for (p2, q) in phi if p2 == p),
"EchelonLB(%s,%s)" % (t, p),
)
model.setObjective(
quicksum(
f[t, p] * y[t, p] + c[t, p] * x[t, p] + H[t, p] * E[t, p]
for t in Ts
for p in P
),
GRB.MINIMIZE,
)
model.update()
model.__data = y, x, E
return model
def mk_instance(seed, NPERIODS, NRESOURCES, MAXPATHS, MAXPATHLEN, factor):
random.seed(seed)
np.random.seed(seed)
random_state = np.random.RandomState(42)
G = nx.DiGraph()
G.add_node(0) # super-source node (to be removed in the end)
G.add_node(-1) # super-target node (to be removed in the end)
NPATHS = random.randint(1, MAXPATHS)
for i in range(NPATHS):
# determing inserting points for a new path in the current graph
paths = list(nx.all_simple_paths(G, source=0, target=-1))
if paths:
path = random.choice(paths)
src = random.randint(0, len(path) - 2)
tgt = random.randint(src + 1, len(path) - 1)
# construct path and insert it in the graph
c = f"{i}-"
N = random.randint(1, MAXPATHLEN)
Gi = nx.path_graph(N, create_using=nx.DiGraph)
map = {j: f"{c}{j}" for j in Gi.nodes()}
nx.relabel_nodes(Gi, map, copy=False)
G = nx.compose(Gi, G)
G.add_edge(0, f"{c}0")
G.add_edge(f"{c}{N-1}", -1)
# connect new path to points (src, tgt) on previous paths
if paths:
G.add_edge(path[src], f"{c}0")
G.add_edge(f"{c}{N-1}", path[tgt])
# remove super-source and super-target nodes, relabel the others
G.remove_node(0)
G.remove_node(-1)
mapping = {node: (i + 1) for (i, node) in enumerate(nx.topological_sort(G))}
nx.relabel_nodes(G, mapping, copy=False)
for (i, j) in sorted(G.edges()): # "sorted" is necessary for deterministic behavior
G.edges[i, j]["phi"] = random.randint(1, 5)
# create lot-sizing data
# start by determining the set of products and conversion (phi)
PRODUCTS = sorted(
list(G.nodes())
) # "sorted" is necessary for deterministic behavior
FINAL = [i for i in PRODUCTS if G.out_degree(i) == 0]
phi = nx.get_edge_attributes(G, "phi")
# now, prepare data for the graph generated
T = NPERIODS
NPRODUCTS = len(PRODUCTS)
PERIODS = range(1, NPERIODS + 1)
K = range(1, NRESOURCES + 1)
f, g, c, d, h, UB, a, M = {}, {}, {}, {}, {}, {}, {}, {}
for t in PERIODS:
for p in PRODUCTS:
if p in FINAL:
d[t, p] = random.randint(0, 10)
else:
d[t, p] = 0
# draw resources required and setup times
a = {(t, k, p): random.randint(0, 3) for t in PERIODS for k in K for p in PRODUCTS}
g = {(t, p): random.randint(1, 3) for t in PERIODS for p in PRODUCTS}
# determine total demand, following graph from target to source
# follow path from final products to raw materials
v = set()
for (i, j) in phi:
v.add(i)
v.add(j)
pred, succ = {}, {}
for i in v:
pred[i] = set()
succ[i] = set()
for (i, j) in phi:
pred[j].add(i)
succ[i].add(j)
# set of vertices corresponding to end products:
final = set(i for i in v if len(succ[i]) == 0)
demand = defaultdict(int)
for t in PERIODS:
for p in final:
demand[t, p] = d[t, p]
curr = set(final)
done = set()
while len(curr) > 0:
q = curr.pop()
done.add(q)
prev = pred[q]
for p in prev:
demand[t, p] += phi[p, q] * demand[t, q]
for s in succ[p]:
if s not in done:
break
else: # no break
curr.add(p)
total_demand = {p: sum(demand[t, p] for t in PERIODS) for p in PRODUCTS}
total_M = sum(
a[t, k, p] * demand[t, p] + g[t, p]
for t in PERIODS
for k in K
for p in PRODUCTS
)
M = {
(t, k): int(0.5 + total_M / (NPERIODS * NRESOURCES) / factor)
for t in PERIODS
for k in K
}
UB = {
(t, p): int(0.5 + total_M / (NPERIODS * NPRODUCTS) / factor)
for t in PERIODS
for p in PRODUCTS
}
# draw costs
for t in range(1, NPERIODS + 1):
for p in PRODUCTS:
f[t, p] = 100 + 10 * random.randint(1, 10)
c[t, p] = random.randint(1, 10)
h[t, p] = random.randint(1, 10) / 100
return T, K, PRODUCTS, f, g, c, d, h, a, M, UB, phi
seed = 12
NRESOURCES = 3
NPERIODS = 5
MAXPATHS = 10
MAXPATHLEN = 3
T, K, P, f, g, c, d, h, a, M, UB, phi = mk_instance(
seed, NPERIODS, NRESOURCES, MAXPATHS, MAXPATHLEN, factor=0.5
)
エシェロン在庫モデルで用いるSCMGraphクラスを準備しておく．このクラスを用いて，部品展開表を綺麗に描画することができる．
class SCMGraph(nx.DiGraph):
"""
SCMGraph is a class of directed graph with edge weight that can be any object.
"""
def layout(self):
"""
Compute x,y coordinates for the supply chain
The algorithm is based on a simplified version of Sugiyama's method.
First assign each node to the (minimum number of) layers;
Then compute the y coordinate by computing the means of y-values
of adjacent nodes
return the dictionary of (x,y) positions of the nodes
"""
longest_path = nx.dag_longest_path(self)
LayerLB = {}
pos = {}
MaxLayer = len(longest_path)
candidate = set([i for i in self]) - set(longest_path)
for i in candidate:
LayerLB[i] = 0
Layer = defaultdict(list)
for i, v in enumerate(longest_path):
Layer[i] = [v]
LayerLB[v] = i
for w in self.successors(v):
if w in candidate:
LayerLB[w] = LayerLB[v] + 1
L = list(nx.topological_sort(self))
for v in L:
if v in candidate:
Layer[LayerLB[v]].append(v)
candidate.remove(v)
for w in self.successors(v):
if w in candidate:
LayerLB[w] = max(LayerLB[v] + 1, LayerLB[w])
MaxLayer = len(Layer)
for i in range(MaxLayer + 1):
if i == 0:
j = 0
for v in Layer[i]:
pos[v] = (i, j)
j += 1
else:
tmplist = []
for v in Layer[i]:
sumy = 0.0
j = 0.0
for w in self.predecessors(v):
(ii, jj) = pos[w]
sumy += jj
j += 1.0
if j != 0:
temp = sumy / j
else:
temp = j
tmplist.append((temp, v))
tmplist.sort()
order = [v for (_, v) in tmplist]
j = 0
for v in Layer[i]:
pos[order[j]] = (i, j)
j += 1
return pos
def up_order(self):
"""
Generator fuction in the reverse topological order
generate the order of nodes from to demand points to suppliers
"""
degree0 = []
degree = {}
for v in self:
if self.out_degree(v) == 0:
degree0.append(v)
else:
degree[v] = self.out_degree(v)
while degree0:
v = degree0.pop()
yield v
# print v
for w in self.predecessors(v):
degree[w] -= 1
if degree[w] == 0:
degree0.append(w)
def find_ancestors(self):
"""
find the ancestors based on the BOM graph
The set of ancestors of node i is the set of nodes that are reachable from node i (including i).
"""
ancestors = {v: set([]) for v in self}
for v in self.up_order():
ancestors[v] = ancestors[v] | set([v])
for w in self.successors(v):
ancestors[v] = ancestors[v] | ancestors[w]
return ancestors
D = SCMGraph()
D.add_weighted_edges_from([(i, j, phi[i, j]) for (i, j) in phi])
pos = D.layout()
nx.draw(D, pos=pos, with_labels=True)
上で生成したデータを標準定式化を用いて求解する．
print("\n\nstandard model")
model = mils_standard(T, K, P, f, g, c, d, h, a, M, UB, phi)
model.optimize()
status = model.Status
if status == GRB.Status.UNBOUNDED:
print("The model cannot be solved because it is unbounded")
elif status == GRB.Status.OPTIMAL:
print("Opt.value=", model.ObjVal)
standard = model.ObjVal
y, x, I = model.__data
print("%7s%7s%7s%7s%12s%12s" % ("prod", "t", "dem", "y", "x", "I"))
for p in P:
print("%7s%7s%7s%7s%12s%12s" % (p, 0, "-", "-", "-", round(I[0, p].X, 5)))
for t in range(1, T + 1):
print(
"%7s%7s%7s%7s%12s%12s"
% (
p,
t,
d[t, p],
int(y[t, p].X + 0.5),
round(x[t, p].X, 5),
round(I[t, p].X, 5),
)
)
for k in K:
for t in range(1, T + 1):
print(
"resource %3s used in period %s: %12s / %-9s"
% (
k,
t,
sum(a[t, k, p] * x[t, p].X + g[t, p] * y[t, p].X for p in P),
M[t, k],
)
)
elif status != GRB.Status.INF_OR_UNBD and status != GRB.Status.INFEASIBLE:
print("Optimization was stopped with status", status)
else:
print("The model is infeasible; computing IIS")
model.computeIIS()
print("\nThe following constraint(s) cannot be satisfied:")
for cnstr in model.getConstrs():
if cnstr.IISConstr:
print(cnstr.ConstrName)
エシェロン在庫を用いた定式化
品目間の親子関係だけでなく，先祖（部品展開表を表す有向グラフを辿って到達可能な品目の集合）を導入しておく．
集合:
- $Ancestor_{p}$: 品目 $p$ の先祖の集合． 親子関係を表す有向グラフを辿って到達可能な点に対応する品目から構成される集合． 品目 $p$ 自身は含まないものとする
パラメータ:
- $\rho_{pq}$: $q \in Ancestor_p$ のとき，品目 $q$ を $1$ 単位生産するのに必要な製品 $p$ の量
- $H_t^p$: 期 $t$ における品目 $p$ のエシェロン在庫費用． 品目 $p$ を生産することによって得られた付加価値に対する在庫費用を表す。品目 $p$ を製造するのに必要な品目の集合を $Child_p$ としたとき，以下のように定義される $$ H_t^p =h_t^p -\sum_{q \in Child_p} h_t^q \phi_{qp} $$
変数:
- $E_{t}^p$: 期 $t$ における品目 $p$ のエシェロン在庫量；自分と自分の先祖の品目の在庫量を合わせたものであり， 以下のように定義される $$ E_t^p = I_t^p +\sum_{q \in Ancestor_p} \rho_{pq} I_t^q $$
上の記号を用いると，エシェロン在庫を用いた多段階ロットサイズ決定モデルの定式化は，以下のようになる．
$$ \begin{array}{ l l l } minimize & \sum_{t=1}^T \sum_{p \in P} \left( f_t^p y_t^p + c_t^p x_t^p + H_t^p E_t^p \right) & \\ s.t. & \ \ E_{t-1}^p +x_t^p - E_t^p = d_t^p + \sum_{q \in Ancestor_p} \rho_{pq} d_t^q & \forall p \in P, t=1,\ldots,T \\ & \sum_{p \in P_k} x_t^p +\sum_{p \in P_k} g_t^p y_t^p \leq M_t^k & \forall k \in K, t=1,\ldots,T \\ & E_t^p \geq \sum_{q \in Parent_p} \phi_{pq} E_t^q & \forall p \in P, t=1,\ldots,T \\ & x_t^p \leq UB_t^p y_t^p & \forall p \in P, t=1,\ldots,T \\ & E_0^p =0 & \forall p \in P \\ & x_t^p,E_t^p \geq 0 & \forall p \in P, t=1,\ldots,T \\ & y_t \in \{0,1\} & \forall t=1,\ldots,T \end{array} $$上の定式化で，最初の制約は，各期および各品目に対するエシェロン在庫の保存式を表す． より具体的には，品目 $p$ の期 $t-1$ からのエシェロン在庫量 $E_{t-1}^p$ と生産量 $x_t^p$ を加えたものが， 期 $t$ における品目 $p$ の先祖 $q$ の需要量を品目 $p$ の必要量に換算したものの合計 $\sum_{q \in Ancestor_p} \rho_{pq} d_t^q$ と、 自身の需要量 $d_{t}^p$ と、次期へのエシェロン在庫量 $E_t^p$ の和に等しいことを表す． 2番目の制約は，各期の生産時間の上限制約を表す． 3番目の制約は，各品目のエシェロン在庫量が，その親集合の品目のエシェロン在庫量の合計以上であること（言い換えれば実需要量が負にならないこと）を規定する． 4番目の制約は，段取り替えをしない期は生産できないことを表す．
この定式化は1段階モデルと同じ構造をもつので、強化式を加えたり、施設配置定式化のような強い定式化に変形することが可能である。
実在庫量の上下限制約を付加したい場合には， エシェロン在庫モデルにおける3番目の制約のスラック変数が、実在庫量になっているので、 制約を以下のように書き換えた後で、実在庫量に対する上下限制約を付加すれば良い。
$$ E_t^p + I_t^p = \sum_{q \in Parent_p} \phi_{pq} E_t^q \ \ \ \forall p \in P, t=1,\ldots,T $$なお，問題に応じた強い定式化を用いたロットサイズ最適化システムとしてOptLot (https://www.logopt.com/optlot/) が開発されている．
print("\n\nechelon model")
model = mils_echelon(T, K, P, f, g, c, d, h, a, M, UB, phi)
model.optimize()
status = model.Status
if status == GRB.Status.UNBOUNDED:
print("The model cannot be solved because it is unbounded")
elif status == GRB.Status.OPTIMAL:
print("Opt.value=", model.ObjVal)
y, x, E = model.__data
print("%7s%7s%7s%7s%12s%12s" % ("t", "prod", "dem", "y", "x", "E"))
for p in P:
print("%7s%7s%7s%7s%12s%12s" % ("t", p, "-", "-", "-", round(E[0, p].X, 5)))
for t in range(1, T + 1):
print(
"%7s%7s%7s%7s%12s%12s"
% (
t,
p,
d[t, p],
int(y[t, p].X + 0.5),
round(x[t, p].X, 5),
round(E[t, p].X, 5),
)
)
for k in K:
for t in range(1, T + 1):
print(
"resource %3s used in period %s: %12s / %-9s"
% (
k,
t,
sum(a[t, k, p] * x[t, p].X + g[t, p] * y[t, p].X for p in P),
M[t, k],
)
)
elif status != GRB.Status.INF_OR_UNBD and status != GRB.Status.INFEASIBLE:
print("Optimization was stopped with status", status)
else:
print("The model is infeasible; computing IIS")
model.computeIIS()
print("\nThe following constraint(s) cannot be satisfied:")
for cnstr in model.getConstrs():
if cnstr.IISConstr:
print(cnstr.ConstrName)