import random
import math
import numpy as np
from gurobipy import Model, quicksum, GRB, multidict
# from mypulp import Model, quicksum, GRB, multidict
from scipy.stats import norm
新聞売り子問題
ここでは，新聞売り子モデル（newsboy model）とよばれる古典的な確率的在庫モデルを紹介する． 新聞売り子モデルは，以下のようなシナリオに基づく．
新聞の売り子が，$1$種類の新聞を販売している． 新聞が売れる量（需要量）は， 経験からある程度推測できると仮定し，確率変数として与えられているものとする． いま，売れ残りのときの在庫費用と， 品切れのときの品切れ費用の和が最小になるように仕入れ量を決めるものとする． どれだけの量を仕入れれば良いだろうか?
解析に必要な記号を導入する．
- $h$: 新聞 $1$部が売れ残ったときに課せられる在庫費用．正の値を仮定する．
- $b$: 新聞 $1$部が品切れしたときに課せられる品切れ費用．正の値を仮定する．
- $D$: 新聞の需要量を表す確率変数． 非負で連続な確率変数であり，その分布関数を $F(x)$（微分可能と仮定），密度関数を $f(x)$ と記す． $$ F(x) =\Pr \{D \leq x \} $$ $$ f(x)=\frac{\partial F(x)}{\partial x} $$
仕入れ量が $s$ のときの総費用の期待値 $C(s)$ は， $$ C(s)= \mbox{E} \left[ h [s-D]^+ +b [s-D]^- \right] $$ となる． ここで，$[\cdot]^+$ は $\max\{0,\cdot \}$ を表し，$[\cdot]^-$ は $-\min\{0,\cdot\}$ を表す．
これは， $$ \begin{array}{l l l} C(s) &= &h \displaystyle\int_0^{\infty} \max\{s-x,0\} f(x) dx + b \displaystyle\int_0^{\infty} \max\{x-s,0\} f(x) dx \\ &= & h \displaystyle\int_0^{s}(s-x) f(x) dx + b \int_s^{\infty} (x-s) f(x) dx \end{array} $$
と書ける．変数 $s$ による1階偏微分は，Leibnizの規則を使うと， $$ \frac{ \partial C(s)}{\partial s}= h \int_0^{s} 1 f(x) dx + b \int_s^{\infty} (-1) f(x) dx =h F(s) -b(1-F(s)) $$ となる．2階偏微分は $$ \frac{ \partial^2 C(s)}{\partial s^2}= (h+b) f(s) $$ となり，需要量が非負の確率変数であり，かつ $h$ も $b$ も正であるという仮定の下では，非負の値をとる． よって，$C(s)$ は $s$ に関する凸関数であるので， $C(s)$ を最小にする $s$ は，$\partial C(s)/ \partial s = h F(s) -b(1-F(s))=0$ より， $$ F(s^*) =\frac{b}{b+h} $$ を満たす $s^*$ になる．
ここで，$F(s^*)$ は需要が仕入れ量 $s^*$ を超えない確率， 言い換えれば品切れを起こさない確率を表していることに注意すると， 上式は，$b/(b+h)$ が品切れを起こさない確率と同じになるように $s^*$ を設定すれば 良いことを表している． $b/(b+h)$ を臨界率（ctitical ratio）とよび， $\omega$ と記す．
最適な仕入れ量 $s^*$ は，分布関数の逆関数を計算すれば良い． $$ s^* =F^{-1} (\omega) $$
分布関数の逆関数は，パーセント点関数(percent point function)とよばれ，SciPyでは，確率変数オブジェクトのppfメソッドで計算できる．
例として，平均100，標準偏差 10 の正規分布にしたがう需要をもつ商品に対して，品切れ費用が 100円，在庫費用が10円のときの仕入れ量を求めてみよう．
b = 100
h = 10
omega = b / (b + h)
print("臨界率=", omega)
print("仕入れ量=", norm.ppf(omega, loc=100, scale=10))
経済発注量問題
経済発注量問題 (economic ordering quantity problem)はサイクル在庫を決めるための古典である．
単一品目の経済発注量問題（Harrisのモデル）は以下の仮定に基づく．
- 品目（商品，製品）は一定のスピードで消費されており， その使用量（これを需要量と呼ぶ）は$1$日あたり $d (\in \mathbf{R}_+)$ 単位である． ここで $\mathbf{R}_+$ は，非負の実数全体の集合を表す記号である．
- 品目の品切れは許さない．
- 品目は発注を行うと同時に調達される． 言い換えれば発注リード時間（注文してから品目が到着するまでの時間）は $0$である．
- 発注の際には，発注量によらない固定的な費用（これを発注費用と呼ぶ） $F (\in \mathbf{R}_+)$ 円が課せられる．
- 在庫保管費用は保管されている在庫量に比例してかかり， 品目 $1$個あたりの保管費用は $1$日で $h (\in \mathbf{R}_+)$ 円とする．
- 考慮する期間は無限期間とする．
- 初期在庫は $0$とする．
上の仮定の下で，$1$日あたりの総費用を最小化する発注方策を 求めることが，ここで考える問題の目的である．
容易にわかるように，最適な発注方策は以下の $2$つの性質を満たす．
- 定常: 方策が時間に依存しない．
- 在庫ゼロ発注性: 在庫量が $0$ になったときのみ発注する．
上の性質の下では， 在庫レベルの時間的な経過を表すグラフはノコギリの歯型になり， 最適な発注方策を求める問題は $1$回あたりの発注量 $Q$ 求める問題に帰着される．
発注を行う間隔をサイクル時間とよび $T$と書く． 発注量 $Q$ とサイクル時間 $T$ の間には $$ d= \frac{Q}{T} $$ の関係がある． 需要量（需要の速度）$d$ が一定であるという仮定の下では， 発注量 $Q$ を求めることとサイクル時間 $T$ を求めることは 同じである．
まず，発注を $1$回行う間（$1$周期あたり）の総費用を考えよう． 総費用は発注費用と在庫保管費用の和である． サイクル時間 $T$ 内では発注は $1$回だけであり， 在庫保管費用は在庫レベルの面積に比例する． よって，$1$周期あたりの総費用は $$ F + \frac{hTQ}{2} $$ となる．単位時間（$1$日）あたりの費用は，これを $T$ で除することにより $$ \frac{F}{T} + \frac{hQ}{2} $$ となる． $Q$ を消去し，サイクル時間 $T$ だけの式に変形することによって $$ \frac{F}{T} + \frac{hdT}{2} $$ を得る．
この問題は解析的に（微分を用いて）解くことができ，最適なサイクル時間 $T^*$ は $$ T^* = \sqrt{ \frac{2F}{hd} } $$ 最適値 $f(T^*)$ は $$ \sqrt{ 2 F h d} $$ となることが示される．
複数品目経済発注量問題
複数の品目を小売店に卸している倉庫における経済発注量モデルを考える． この場合には，上のように簡単な公式で解析的には解くことができない．ここでは，非線形関数を区分的線形関数で近似することにより，解くことを考える．
品目の集合を $I$ とする． 倉庫には，置けるものの量に上限があり，これを容量制約と呼ぶ． 品目 $i(\in I)$ の大きさを $w_i$ とし，倉庫に置ける量の上限を $W$ とする． また，各品目は，すべて異なるメーカーに注文するので， 発注費用は各品目を注文するたびにかかるものとし，品目 $i$ の発注費用を $F_i$ とする． 品目 $i$ の在庫保管費用を $h_i$，需要量を $d_i$ としたとき，容量制約を破らない条件の下で， 総費用が最小になる発注方策を導こう．
Harrisのモデルと同様に，発注費用と在庫費用の和を最小化するが，容量制約を表現するための制約が付加される． 非線形の目的関数をもつ数理最適化モデルとして定式化すると， 容量を考慮した複数品目モデルは，以下のように書ける．
$$ \begin{array}{l l l} minimize & \sum_{i \in I} \frac{F_i}{T_i} + \frac{h_i d_i T_i}{2} & \\ s.t. & \sum_{i \in I} w_i d_i T_i \leq W & \\ & T_i > 0 & \forall i \in I \end{array} $$ここでは，区分的線形近似を用いて定式化を行い，求解する．
def eoq(I, F, h, d, w, W, a0, aK, K):
"""eoq -- multi-item capacitated economic ordering quantity model
Parameters:
- I: set of items
- F[i]: ordering cost for item i
- h[i]: holding cost for item i
- d[i]: demand for item i
- w[i]: unit weight for item i
- W: capacity (limit on order quantity)
- a0: lower bound on the cycle time (x axis)
- aK: upper bound on the cycle time (x axis)
- K: number of linear pieces to use in the approximation
Returns a model, ready to be solved.
"""
# construct points for piecewise-linear relation, store in a,b
a, b = {}, {}
delta = float(aK - a0) / K
for i in I:
for k in range(K):
T = a0 + delta * k
a[i, k] = T # abscissa: cycle time
b[i, k] = (
F[i] / T + h[i] * d[i] * T / 2.0
) # ordinate: (convex) cost for this cycle time
model = Model("multi-item, capacitated EOQ")
x, y, c = {}, {}, {}
for i in I:
x[i] = model.addVar(vtype="I", name="x(%s)" % i) # cycle time for item i
c[i] = model.addVar(vtype="C", name="c(%s)" % i) # total cost for item i
for k in range(K):
y[i, k] = model.addVar(ub=1, vtype="C", name="y(%s,%s)" % (i, k))
model.update()
for i in I:
model.addConstr(quicksum(y[i, k] for k in range(K)) == 1)
model.addConstr(quicksum(a[i, k] * y[i, k] for k in range(K)) == x[i])
model.addConstr(quicksum(b[i, k] * y[i, k] for k in range(K)) == c[i])
model.addConstr(quicksum(w[i] * d[i] * x[i] for i in I) <= W)
model.setObjective(quicksum(c[i] for i in I), GRB.MINIMIZE)
model.update()
model.__data = x, y
return model
I, F, h, d, w = multidict(
{1: [300, 10, 10, 20], 2: [300, 10, 30, 40], 3: [300, 10, 50, 10]}
)
W = 2000
K = 1000
a0, aK = 0.1, 10
model = eoq(I, F, h, d, w, W, a0, aK, K)
model.optimize()
x, y = model.__data
EPS = 1.0e-6
for v in x:
if x[v].X >= EPS:
print(x[v].VarName, x[v].X)
print("Obj:", model.ObjVal)
途絶を考慮した新聞売り子問題
調達先の途絶だけでなく，需要も不確実性をもつ場合には，陽的な公式を得ることは難しいが， 以下の簡単な数理最適化モデルを用いることによって，最適な発注量を求めることができる．
数理最適化モデルの定式化で用いる記号は，以下の通り．
- $h$: （品目 $1$ 個あたり，$1$ 期間あたりの）在庫費用
- $b$: （品目 $1$ 個あたり，$1$ 期間あたりの）バックオーダー費用
- $S$: 途絶と需要のシナリオの集合；シナリオを表す添え字を $s$ と記す．$n$ 期連続で途絶するシナリオは，確率 $\pi_n$ で発生するものとし， その際には，$n+1$ 期分のバックオーダーされた需要量が発生するものとする
$d_s$: シナリオ $s$ における品目の需要量；基本となる需要が各期ごとに独立な正規分布 $N(\mu,\sigma^2)$ と仮定した場合には， $n$ 期連続で途絶するシナリオに対しては，平均 $\mu (n+1)$，標準偏差 $\sigma \sqrt{n+1}$ の正規分布にしたがう
$p_s$: シナリオ $s$ の発生確率；$n$ 期連続で途絶するシナリオは，確率 $\pi_n$ で発生させる．$n$ が大きくなると $\pi_n$ は小さくなるので， 有限個のシナリオを生成するために適当な $n$ で打ち切り，各 $n$ に対して定数個のランダムな需要シナリオが等確率で発生するように設定する
$x$: 発注量（基在庫レベル）を表す（即時決定）変数
- $I_s$: シナリオ $s$ における在庫量を表す（リコース）変数
- $B_s$: シナリオ $s$ におけるバックオーダー量を表す（リコース）変数
上で定義した記号を用いることによって，途絶を考慮した新聞売り子モデルは，以下のように定式化できる．
$$ \begin{array}{l l l} minimize & \sum_{s \in S} p_s \left( h I_s+ b B_s \right) & \\ s.t. & x +B_s = d_s +I_s & \forall s \in S \\ & x \geq 0 & \\ & I_s \geq 0 & \forall s \in S \\ & B_s \geq 0 & \forall s \in S \end{array} $$上のモデルを用いることによって，途絶リスクと基在庫レベルの関係に対する洞察を得ることができる．
また，リスク要因を考慮するためには，CVaR(conditional value at risk)最小化モデルが有効である． ここでは，上で考えた途絶を考慮した新聞売り子モデルに対してCVaRを評価尺度としたモデルを考える．
与えられた確率 $0<\beta<1$ に対して，費用が閾値 $y$ を超えない確率が $\beta$ 以上になるような 最小の $y$ を $\beta$-VaRと呼ぶ． $\beta$-VaRは計算しにくいので， 代用品として$\beta$-VaRを超えた条件の下での期待値である $\beta$-CVaRを用いる
シナリオ $s$ における費用を表す変数 $f_s$，$f_s$ が $\alpha$ を超過する量を表す変数を $V_s$ とすると， CVaR最小化モデルは，以下のようになる．
$$ \begin{array}{l l l} minimize & y +\frac{1}{1-\beta} \sum_{s \in S} p_s V_s & \\ s.t. & f_s= h I_s+ b B_s & \forall s \in S \\ & V_s \geq f_s -y & \forall s \in S \\ & x +B_s = d_s +I_s & \forall s \in S \\ & x \geq 0 & \\ & I_s \geq 0 & \forall s \in S \\ & B_s \geq 0 & \forall s \in S \\ & V_s \geq 0 & \forall s \in S \\ \end{array} $$h = 1.0 # holding cost
b = 100.0 # backorder cost
T = 1 # 計画期間
mu = 100 # 需要の平均
sigma = 10 # 需要の標準偏差
risk_ratio = 0.5
beta = 0.8
d = {} # 需要量
prob = {}
random.seed(12)
S = 100 # シナリオの数(0..S-1)
#単一期間の需要を正規分布で生成
prob = {}
ND = norm(mu, sigma)
prob[0] = ND.cdf(0)
d[0, 0] = 0
for i in range(1, S):
d[0, i] = i
prob[i] = ND.pdf(i)
model = Model()
x, B, I, V, f = {}, {}, {}, {}, {}
for s in range(S):
V[s] = model.addVar(vtype="C", name=f"V({s})")
f[s] = model.addVar(vtype="C", name=f"f({s})")
alpha = model.addVar(vtype="C", name="alpha") # VaR
for t in range(T):
x[t] = model.addVar(vtype="C", name= f"x({t})")
for s in range(S):
B[t, s] = model.addVar(vtype="C", name= f"B({t},{s})" )
I[t, s] = model.addVar(vtype="C", name= f"I({t},{s})" )
for s in range(S):
I[-1, s] = 0
B[-1, s] = 0
model.update()
for s in range(S):
model.addConstr(
f[s] == quicksum(h * I[t, s] + b * B[t, s] for t in range(T)),
name= f"f_evaluate({s})",
)
model.addConstr(V[s] >= f[s] - alpha, name= f"V_evaluate({s})")
for t in range(T):
model.addConstr(
I[t - 1, s] + x[t] + B[t, s] == d[t, s] + I[t, s] + B[t - 1, s],
name= f"flow_cons({t},{s})",
)
model.setObjective(
(1 - risk_ratio)
* quicksum(
prob[s] * (h * I[t, s] + b * B[t, s]) for s in range(S) for t in range(T)
)
+ risk_ratio* (alpha + 1 / (1 - beta) * quicksum(prob[s] * V[s] for s in range(S))),
GRB.MINIMIZE,
)
model.optimize()
cost = [
(sum((h * I[t, s].X + b * B[t, s].X) for t in range(T)), prob[s])
for s in range(S)
]
risk = alpha.X + 1 / (1 - beta) * sum(prob[s] * V[s].X for s in range(S))
安全在庫配置問題
ここでは，直列多段階システムにおける，安全在庫配置問題を考える． このモデルは以下の仮定に基づく．
- 単一の品目を供給するための在庫点が $n$ 個直列に並んでいるものとする．第$n$段階は原料の調達を表し，第$1$段階は最終需要地点における消費を表す． 第$i$段階の在庫点は，第$i+1$段階の在庫点から補充を受ける．
- 各段階は，各期における最終需要地点における需要量だけ補充を行うものとする．
第$1$段階で消費される品目の $1$日あたりの需要量は，期待値 $\mu$ をもつ定常な分布をもつ． また，$t$日間における需要の最大値を $D(t)$ とする． たとえば，需要が平均値 $\mu$，標準偏差 $\sigma$ の正規分布にしたがい，意思決定者が品切れする 確率を安全係数 $z(>0)$ で制御していると仮定したときには，$D(t)$ は $$ D(t)= \mu t + z \sigma \sqrt{t} $$ と書ける．
第$i$段階における品目の生産時間は，$T_i$ 日である． $T_i$ には各段階における生産時間の他に，待ち時間および輸送時間も含めて考える．
第$i$段階の在庫点は，第$i-1$段階の発注後，ちょうど $L_i$ 日後に品目の補充を行うことを保証しているものとする．これを（第$i$段階の）保証リード時間(guaranteed lead time, guaranteed service time)と呼ぶ． なお，第$1$段階（最終需要地点）における保証リード時間 $L_1$ は，事前に決められている定数とする．
在庫（保管）費用は保管されている在庫量に比例してかかり， 第$i$段階における在庫費用は，品目 $1$ 個，$1$日あたり $h_i (\in \mathbf{R}_+)$ 円とする．
上の仮定の下で，$1$ 日あたりの在庫費用を最小化するように，各段階における安全在庫レベルと保証リード時間を決めることが，ここで考える問題の目的である．
第$i$段階の在庫点を考える．この地点に在庫を補充するのは，第$i+1$段階の在庫点であり， そのリード時間は $L_{i+1}$ であることが保証されている． したがって，それに $T_i$ を加えたものが，補充の指示を行ってから第$i$段階が生産を完了するまでの時間となる． これを，補充リード時間(replenishment lead time)と呼ぶ． また，第$i$段階は第$i-1$段階に対して，リード時間 $L_i$ で補充することを保証している． したがって，第$i$段階では，補充リード時間から $L_i$ を減じた時間内の最大需要に相当する在庫を保持していれば， 在庫切れの心配がないことになる． 補充リード時間から $L_i$ を減じた時間（$L_{i+1}+T_i-L_i$）を 正味補充時間(net replenishment time)とよび， $x_i$ と記す．
第$i$段階における安全在庫量 $I_i$ は， 正味補充時間内における最大需要量から平均需要量を減じた量であるので， $$ I_i= D (x_i) - x_i \mu $$ となり，これに $h_i$ を乗じたものが安全在庫費用になる．
第$i$段階の安全在庫費用を正味補充時間 $x_i$ の凹関数として $f_i(x_i)$ とすると，直列多段階安全在庫配置問題は，以下のように定式化できる．
$$ \begin{array}{lll} minimize & \sum_{i=1}^n h_i f(x_i) & \\ s.t. & x_i + L_i -L_{i+1} = T_i & \forall i=1,2,\ldots,n \\ & L_{n+1} =0 & \\ & L_i \geq 0 & \forall i=2,3,\ldots,n \\ & x_i \geq 0 & \forall i=1,2,\ldots,n \end{array} $$以下では，区分的線形近似を用いた定式化を示す．
def make_data():
n = 30 # number of stages
z = 1.65 # for 95% service level
sigma = 100 # demand's standard deviation
h = {} # inventory cost
T = {} # production lead time
h[n] = 1
for i in range(n - 1, 0, -1):
h[i] = h[i + 1] + random.randint(30, 50)
K = 0 # number of segments (=sum of processing times)
for i in range(1, n + 1):
T[i] = random.randint(3, 5) # production lead time at stage i
K += T[i]
return z, sigma, h, T, K, n
z, sigma, h, T, K, n = make_data()
def convex_comb_sos(model, a, b):
"""convex_comb_sos -- add piecewise relation with gurobi's SOS constraints
Parameters:
- model: a model where to include the piecewise linear relation
- a[k]: x-coordinate of the k-th point in the piecewise linear relation
- b[k]: y-coordinate of the k-th point in the piecewise linear relation
Returns the model with the piecewise linear relation on added variables x, f, and z.
"""
K = len(a) - 1
z = {}
for k in range(K + 1):
z[k] = model.addVar(
lb=0, ub=1, vtype="C"
) # do not name variables for avoiding clash
x = model.addVar(lb=a[0], ub=a[K], vtype="C")
f = model.addVar(lb=-GRB.INFINITY, vtype="C")
model.update()
model.addConstr(x == quicksum(a[k] * z[k] for k in range(K + 1)))
model.addConstr(f == quicksum(b[k] * z[k] for k in range(K + 1)))
model.addConstr(quicksum(z[k] for k in range(K + 1)) == 1)
model.addSOS(GRB.SOS_TYPE2, [z[k] for k in range(K + 1)])
return x, f, z
def ssa(n, h, K, f, T):
"""ssa -- multi-stage (serial) safety stock allocation model
Parameters:
- n: number of stages
- h[i]: inventory cost on stage i
- K: number of linear segments
- f: (non-linear) cost function
- T[i]: production lead time on stage i
Returns the model with the piecewise linear relation on added variables x, f, and z.
"""
model = Model("safety stock allocation")
# calculate endpoints for linear segments
a, b = {}, {}
for i in range(1, n + 1):
a[i] = [k for k in range(K)]
b[i] = [f(i, k) for k in range(K)]
# x: net replenishment time for stage i
# y: corresponding cost
# s: piecewise linear segment of variable x
x, y, s = {}, {}, {}
L = {} # service time of stage i
for i in range(1, n + 1):
x[i], y[i], s[i] = convex_comb_sos(model, a[i], b[i])
if i == 1:
L[i] = model.addVar(ub=0, vtype="C", name="L[%s]" % i)
else:
L[i] = model.addVar(vtype="C", name="L[%s]" % i)
L[n + 1] = model.addVar(ub=0, vtype="C", name="L[%s]" % (n + 1))
model.update()
for i in range(1, n + 1):
# net replenishment time for each stage i
model.addConstr(x[i] + L[i] == T[i] + L[i + 1])
model.setObjective(quicksum(h[i] * y[i] for i in range(1, n + 1)), GRB.MINIMIZE)
model.update()
model.__data = x, s, L
return model
def make_data():
n = 30 # number of stages
z = 1.65 # for 95% service level
sigma = 100 # demand's standard deviation
h = {} # inventory cost
T = {} # production lead time
h[n] = 1
for i in range(n - 1, 0, -1):
h[i] = h[i + 1] + random.randint(30, 50)
K = 0 # number of segments (=sum of processing times)
for i in range(1, n + 1):
T[i] = random.randint(3, 5) # production lead time at stage i
K += T[i]
return z, sigma, h, T, K, n
random.seed(1)
z, sigma, h, T, K, n = make_data()
def f(i, k):
return sigma * z * math.sqrt(k)
model = ssa(n, h, K, f, T)
model.optimize()
x, s, L = model.__data
print("%10s%10s%10s%10s" % ("Stage", "x", "L", "T"))
for i in range(1, n + 1):
print("%10s%10s%10s%10s" % (i, x[i].X, L[i].X, T[i]))
print("Obj:", model.ObjVal)
適応型在庫最適化問題
需要だけが不確実性を含んでいる場合には， 過去の需要系列の履歴のアフィン関数として発注量を決定する適応型モデルが有効であることが知られている．
ここでは，多期間の確率的在庫モデルを考える．このモデルによって，新聞売り子モデルでは単一期間であったため 考慮されなかった複数期間にまたがる影響を考えることができる． また，このモデルは多段階，複数調達などの拡張モデルの基礎となる．
このモデルでは，発注量をシナリオに依存しない変数（即時決定変数）とし， 在庫量ならびにバックオーダー量（品切れ量）はシナリオに依存する変数（リコース変数）とする． 需要が大きかったり，供給が途絶した場合には，在庫不足のため需要が満たされないことが想定される． この際，顧客が再び品目が到着するまで待ってくれる場合（バックオーダー; backorder）と， 需要が消滅してしまう場合（品切れ，販売機会の逸失; lost sales）に分けて考える必要がある． （もちろん，顧客需要の一部が逸失し，一部がバックオーダーされるという場合もあるが，ここでは両極端の場合のみを考える．）
以下に定式化に必要な記号を，パラメータ（定数）と変数に分けて記述する．
パラメータ:
- $T$: 計画期間数．期を表す添え字を $1,2,\ldots,t,\ldots,T$，添え字の集合を $\Upsilon =\{1,2,\ldots,T \}$ と記す
- $S$: 途絶と需要のシナリオの集合；シナリオを表す添え字を $s$ と記す
- $h$: （品目 $1$ 個あたり，$1$ 期間あたりの）在庫費用
- $b$: （品目 $1$ 個あたり，$1$ 期間あたりの）バックオーダー費用（品切れ費用）
- $M$: 発注量の上限．以下では，これを発注の容量と呼ぶ
- $d_t^s$: シナリオ $s$ における期 $t$ の品目の需要量
- $p_s$: シナリオ $s$ の発生確率
変数:
- $I^s_t$: シナリオ $s$ における期 $t$ の在庫量．より正確に言うと，期 $t$ の期末の在庫量． ここで，$I_0^s, B_0^s$ は定数として与えられているものと仮定する
- $B^s_t$: シナリオ $s$ における期 $t$ のバックオーダー量（品切れ量）
- $x_t$: 期 $t$ における発注量
シナリオ $s$ の期 $t$ における発注量を，過去 $j (=1,2,\ldots,\theta)$ 期の需要の $y_j$ 倍を発注量に加えたアフィン関数とする． $$ X_t^s =\sum_{j=1}^{\min \{t-1,\theta\}} d_{t-j}^s y_j +x_t $$
容量 $M$ が小さい場合には，すべてのシナリオ $s$ に対して $X_t^s \leq M$ を満たすために $y_j$ が $0$ になるため， 適応型モデルは静的発注モデルと同じ性能を示す．そこで，需要量 $d_t^s$ が容量 $M$ を超過した量 $E_t^s$ の $z_j$ 倍を 発注量から減じた以下のアフィン関数を定義する．
$$ X_t^s =\sum_{j=1}^{\min \{t-1,\theta\}} \left( d_{t-j}^s y_j - E_{t-j}^s z_j \right) +x_t $$これは，過去の需要量に対する区分的線形なアフィン関数で発注量を決めていることに他ならない． 以下では，この式を用いた適応型モデルを拡張適応型モデルと呼ぶ．
期待値最小化を目的とした適応型モデルは，以下のように定式化できる．
$$ \begin{array}{l l l} minimize & \sum_{s \in S} p_s \sum_{t \in \Upsilon} \left( h I_t^s+ b B_t^s \right) & \\ s.t. & X_t^s =\sum_{j=1}^{\min \{t-1,\theta\}} \left\{ d_{t-j}^s y_j \right\}+x_t & \forall t \in \Upsilon ;s \in S \\ & I_{t-1}^s +\delta^s_t X_t^s +B_t^s = d_t^s +I_t^s+B_{t-1}^s & \forall t \in \Upsilon ;s \in S \\ & 0 \leq X_t^s \leq M & \forall t \in \Upsilon ;s \in S \\ & x_t \geq 0 & \forall t \in \Upsilon \\ & I_{t}^s \geq 0 & \forall t \in \Upsilon ;s \in S \\ & B_t^s \geq 0 & \forall t \in \Upsilon ;s \in S \end{array} $$品切れ（販売機会逸失）の場合には，上の定式化の2番目の制約式を以下のように変更する． $$ I_{t-1}^s +\delta^s_t x_t +B_t^s = d_t^s +I_t^s \ \ \ \forall t \in \Upsilon ; s \in S $$
リスクを加味したCVaR最小化モデルも，同様に構築できる．
$$ \begin{array}{l l l} minimize & y +\frac{1}{1-\beta} \sum_{s \in S} p_s V_s & \\ s.t. & f_s= \sum_{t \in \Upsilon} \left(h I_t^s+ b B_t^s \right) & \forall s \in S \\ & V_s \geq f_s -y & \forall s \in S \\ & I_{t-1}^s + \delta^s_t x_t +B_t^s = d_t^s +I_t^s +B_{t-1}^s & \forall t \in \Upsilon ;s \in S \\ & 0 \leq x_t \leq M & \forall t \in \Upsilon \\ & I_{t}^s \geq 0 & \forall t \in \Upsilon ;s \in S \\ & B_t^s \geq 0 & \forall t \in \Upsilon ;s \in S \\ & V_s \geq 0 & \forall s \in S \\ \end{array} $$適応型モデルは，需要の予測を加味したモデルととらえることができるが，強化学習モデル と考えることもできる． 過去のデータをもとに将来の需要や途絶のシナリオを作成し， 最適なパラメータベクトル $(y,Y,x)$ を計算（学習）しておき，実際の運用の際には，これらのパラメータを用いて 発注量 $X$ を決める．通常の強化学習と異なる点は，将来の需要や途絶のシナリオを作成しておくことであるが， これが難しい場合には，過去のデータをシナリオとして用いて直接パラメータ $(y,Y,x)$ を推定しても良い．
def adaptive_inventory(T, h, b, M, mu, sigma, S, d, theta):
"""
the adaptive invenntory problem:
theta: threshold for the adaptive policy
"""
prob = 1.0 / float(S)
model = Model("robust inv")
x, X, y, B, I, V, f, z = {}, {}, {}, {}, {}, {}, {}, {}
beta = 0.8
risk_ratio = 0.0
for s in range(S):
V[s] = model.addVar(
obj=risk_ratio * prob / (1 - beta), vtype="C", name=f"V({s})"
)
f[s] = model.addVar(vtype="C", name=f"f({s})")
alpha = model.addVar(vtype="C", name="alpha")
for j in range(1, theta + 1):
y[j] = model.addVar(lb=-GRB.INFINITY, vtype="C", name=f"y({j})")
z[j] = model.addVar(lb=-GRB.INFINITY, vtype="C", name=f"z({j})")
for t in range(T):
x[t] = model.addVar(vtype="C", name=f"x({t})")
for s in range(S):
B[t, s] = model.addVar(vtype="C", name=f"B({t},{s})")
I[t, s] = model.addVar(vtype="C", name=f"I({t},{s})")
X[t, s] = model.addVar(ub=M, vtype="C", name=f"X({t},{s})")
for s in range(S):
I[-1, s] = 0
model.update()
for s in range(S):
model.addConstr(
f[s] == quicksum(h * I[t, s] + b * B[t, s] for t in range(T)),
name=f"f_evaluate({s})",
)
model.addConstr(V[s] >= f[s] - alpha, name=f"V_evaluate({s})")
for t in range(T):
model.addConstr(
X[t, s]
== quicksum(
d[t - j, s] * y[j] + max(d[t - j, s] - M, 0) * z[j]
for j in range(1, theta + 1)
if t - j >= 0
)
+ x[t],
name=f"adaptive({t},{s})",
)
model.addConstr(
I[t - 1, s] + X[t, s] + B[t, s] == d[t, s] + I[t, s],
name=f"flow_cons({t},{s})",
)
model.setObjective(
(1 - risk_ratio)
* quicksum(
prob * (h * I[t, s] + b * B[t, s]) for s in range(S) for t in range(T)
)
+ risk_ratio
* (alpha + 1 / (1 - beta) * quicksum(prob * V[s] for s in range(S))),
GRB.MINIMIZE,
)
model.optimize()
print("Opt.value=", model.ObjVal)
for v in y:
if y[v].X != 0.0:
print(y[v].VarName, y[v].X)
if z[v].X != 0.0:
print(z[v].VarName, z[v].X)
for v in x:
if x[v].X > 0.001:
print(x[v].VarName, x[v].X)
T = 20 # number of periods
h = 10.0 # holding cost
b = 100.0 # backorder cost
M = 150 # capacity
mu = 100
sigma = 10
S = 100 # number of samples (scenarios)
d = {}
random.seed(1)
# Cyclic demand
weekly_ratio = [0.1, 1.0, 1.2, 1.3, 0.9, 1.5, 0.2] # 0 means Sunday
cycle = len(weekly_ratio)
yearly_ratio = [1.0 + np.sin(i) for i in range(13)] # 0 means January
for s in range(S):
epsilon = 0.0
for t in range(T):
epsilon = int(random.gauss(0, sigma)) # normal distribution
d[t, s] = mu * weekly_ratio[t % cycle] + epsilon
adaptive_inventory(T, h, b, M, mu, sigma, S, d, theta=cycle)
複数エシェロン在庫最適化問題
実際の在庫最適化においては，多段階の一般型ネットワークで最適化を行う必要がある． この問題の総称は，複数エシェロン在庫最適化問題(multi-echelon inventory optimization problem)とよばれる． この問題は，在庫のシミュレーションと最適化を融合した解法によって，（近似的にではあるが）解くことができる． また上で述べた安全在庫配置問題も，実際問題においては一般のネットワークで求解する必要がある． これらの要件を考慮した在庫最適化の統合システムとして MESSA (MEta Safety Stock Allocation system: https://www.logopt.com/messa/ )が開発されている．