在庫最適化と安全在庫配置システム MESSA (MEta Safety Stock Allocation system)
サービスレベルと安全在庫係数
ここでは，小売店に顧客がやってくることによって需要が発生する場合を考える． 顧客は，毎日小売店にやってきて商品を購入するが，その量は気分によって変わるものとする． 一般には，商品が品切れしてしまうことは避けなければならない． 商品を余分に在庫しておく費用が，品切れしてしまう費用より小さい場合には， 商品を余分に在庫して気まぐれな顧客のために準備しておくことが常套手段である． このように，需要の不確実性に対処するためにもつ在庫を安全在庫（safety inventory, safety stock）とよぶ．
いま，簡単のため，毎日の需要量は，定常な分布にしたがっているものとする． ここで「定常」とは，日によって需要の分布が変化しないことを指す． 分布によっては，負の需要量になる可能性があるが，その場合には需要量を $0$ と定義する． たとえば， 需要が正規分布と仮定した場合には，負の部分を除いた正規分布（切断正規分布）を指すものとする．
wb = make_excel_messa()
wb.save("messa-templete.xlsx")
wb = load_workbook("messa-ex2.xlsx")
G, ProcTime, LTUB, z, mu, sigma, h = prepare_opt_for_messa(wb)
#best_cost, best_sol, best_NRT, best_MaxLI, best_MinLT = tabu_search_for_SSA(G, ProcTime, LTUB, z, mu, sigma, h, max_iter = 10, TLLB =1, TLUB =3, seed = 1)
wb = load_workbook("messa-ex1.xlsx")
G, ProcTime, LTUB, z, mu, sigma, h = prepare_opt_for_messa(wb)
best_cost, best_sol, best_NRT, best_MaxLI, best_MinLT = tabu_search_for_SSA(G, ProcTime, LTUB, z, mu, sigma, h, max_iter = 10, TLLB =1, TLUB =3, seed = 1)
wb = messa_for_excel(wb, best_NRT, best_MaxLI, best_MinLT )
wb.save("messa-ex1-result.xlsx")
n_samples =100
n_periods =100
for i in range(n):
if best_NRT[i]>=1:
s, S = approximate_ss(mu[i], sigma[i], best_NRT[i]-1, b[i], h[i], fc[i])
print(s,S)
cost, I = simulate_inventory(n_samples =n_samples, n_periods =n_periods, mu=mu[i], sigma=sigma[i], LT=int(best_NRT[i]), Q=None, R=s, b=b[i], h=h[i], fc=fc[i], S=S)
print(cost)
net = generate_pyvis_for_messa(wb, G, best_sol, best_MaxLI, best_MinLT)
#net.show("example.html")
Wagner-Whitin法による動的ロットサイズ決定問題の求解関数 ww
動的ロットサイズ決定モデルの最適化（モデルについては動的ロットサイズ決定モデルlotsize参照）
需要予測に基づく不確実性を考慮した在庫最適化（forecast）で利用
引数:
- demand: 多期間需要を表すリスト
- fc: 固定費用（定数かリスト）; 関数内で需要と同じ長さのNumPy配列に変換される．
- vc: 変動費用（定数かリスト）; 関数内で需要と同じ長さのNumPy配列に変換される．
- h: 在庫費用（定数かリスト）; 関数内で需要と同じ長さのNumPy配列に変換される．
返値：
- cost: 最適値
- order: 各期の発注量を入れたリスト
fixed= 3
variable= [1,1,3,3,3]
h = 1.
demand=[5,7,3,6,4]
ww(demand, fc =fixed, vc=variable, h=h)
data = np.random.normal(100,10,1000)
fig, best_dist, best_fit_name, best_fit_params = best_distribution(data)
plotly.offline.plot(fig);
print(best_fit_name, best_fit_params)
data = np.random.normal(100,10,1000)
fig, hist_dist = best_histogram(data, nbins = 20)
plotly.offline.plot(fig);
得られた分布に基いてランダムな需要を発生させる．
best_dist.mean(),hist_dist.mean()
得られた分布に対して，在庫管理に必要な統計量を計算する．
print( best_dist.ppf(0.9), hist_dist.ppf(0.9) ) #90%品切れしないための在庫量
print( best_dist.sf(200), hist_dist.sf(200) ) #在庫が200のときの品切れ率
需要分布に対して最適な適合をもつ分布を返す関数 fit_demand
引数：
- demand_df: 需要データフレーム
- cust_selected: 需要の分布を求めたい顧客名のリスト
- prod_selected: 需要の分布を求めたい製品名のリスト
- agg_period: 集約を行う期。例えば１週間を１期とするなら"1w"とする。既定値は1日（"1d"）
- nbins: ヒストグラムのビン数（既定値は50）
- L: リード時間（在庫管理で用いる．L日前からL-1前の需要の合計の分布をとる．）
返値：
dic: 顧客名と製品名のタプルをキーとし，分布の情報を値とした辞書； 分布の情報は以下のタプル
- fig: 誤差が最小の分布の密度関数と度数分布表の図
- fig2: 近似的な表形式の確率関数と度数分布表の図
- best_dist: 誤差最小の分布の密度関数
- hist_dist: 度風分布表のデータから構成された確率関数
- best_fit_name: 誤差最小の分布名
- best_fit_params: 誤差を最小にする分布のパラメータ
demand_df = pd.read_csv(folder+"demand.csv")
#demand_df = pd.read_csv(folder+"demand_all.csv")
cust_selected = [ demand_df["cust"].iloc[11]]
prod_selected = [ demand_df["prod"].iloc[1] ]
#cust_selected = ["前橋市"]
#prod_selected = ["F"]
#print(cust_selected, prod_selected)
dic = fit_demand(demand_df, cust_selected, prod_selected, agg_period ="1d", nbins=50, L=3)
if dic[ cust_selected[0], prod_selected[0] ] is not None:
fig, fig2, best_dist, hist_dist, best_fit_name, best_fit_params = dic[ cust_selected[0], prod_selected[0] ]
print(best_fit_name, best_fit_params)
plotly.offline.plot(fig);
数量割引の考慮
増分割引（incremental discount）: ある一定数量を超えると， 超えた分だけ単位あたりの購入費用が割り引きされる．
総量割引（all-units discount）: ある一定数量を超えると， 一度に購入したすべての量に対して割引が適用される．
増分割引
いま，発注量 $Q$ が分岐点 $\theta$ を超えると， その購入単価が $c_0$ から $c_1$ に（$\theta$ を超えた分だけ）割引されるものとする． 発注の固定費用を $K$ とし，購入費用の合計を，発注量 $Q$ の関数として $c(Q)$ とすると，$c(Q)$ は
$$ c(Q)= \left\{ \begin{array}{l l} 0 & Q=0 \\ K + c_0 Q & 0 <Q \leq \theta \\ K + c_0 \theta + c_1 (Q-\theta) & \theta <Q \end{array} \right. $$と書くことができる．
総費用は，購入費用と在庫費用の和である． 単位時間あたりの購入費用は，一度に $Q$ 単位を購入するときにかかる費用 $c(Q)$ をサイクル時間 $T$ で割った値であるので， $$ \frac{c(Q)}{T}=\frac{c(Q) d}{Q} $$ となる． 在庫費用は，購入費用に対する利子とその他のハンドリング費用や陳腐化費用を分けて考える． 購入費用に対する利子を除いた単位時間・単位数量あたりの在庫保管費用を $h$ とする． 購入費用のうち，発注固定費用を除いた部分は，$c(Q)-K$ であるので，単位数量あたりの購入費用は， $(c(Q)-K)/Q$ となる．利子率を $r$ とすると，単位在庫あたりの金利負担は $r (c(Q)-K)/Q$ となる． したがって在庫費用は，平均在庫量が $Q/2$ であることから， $$ \frac{1}{2} ( h + r (c(Q)-K) ) Q $$ となる． したがって，単位時間あたりの費用関数 $f(Q)$ は， $$ f(Q)= \frac{c(Q) d}{Q} + \frac{1}{2} ( h + r (c(Q)-K) ) Q $$ となる．
$Q>0$ としたとき，$c(Q)$ を以下のように，2本の直線の小さい方をとった関数に変形することができる． $$ c(Q) =\min \left\{ \begin{array}{l} K + c_0 Q \\ K_1 + c_1 Q \end{array} \right. $$ ここで，$K_1=K+(c_0-c_1) \theta$ である．
各直線ごとの問題は，通常の経済発注量モデルと同じ形になるので，最適解を容易に求めることができる．
たとえば，$c(Q)=K_1+c_1 Q$ の場合には， $$ f(Q)=d c_1 +\frac{d K_1}{Q} + \frac{1}{2} (h +r c_1) Q $$ となるので，$Q$ で微分して $0$ と置くことによって，最適発注量 $Q^*_1$ を得ることができる． $$ Q_1^* =\sqrt{ \frac{2 d K_1}{ h+r c_1 } } $$ また，最適費用は， $$ f(Q_1^*) = d c_1 + \sqrt{ 2 K_1 d (h+r c_1) } $$ である．
同様に，$c(Q)=K+c_0 Q$ の場合の最適発注量を求め，それを $Q_0^*$ とする． $$ Q_0^* =\sqrt{ \frac{2 d K}{ h+r c_0 } } $$
$c_0 >c_1$ であり，かつ $K <K_1$ であるので，$Q_0^* <Q_1^*$ となる．
$c(Q)=K_1+c_1 Q$ と $c(Q)=K+c_0 Q$ の場合の費用関数は，$Q=\theta$ で繋がっていることに注意すると， $f(Q_0^*)$ と $f(Q_1^*)$ をくらべて小さい方を選択すれば良いことが言える．
一般に，区分 $j$ に対する購入単価を $c_j (j=0,1,\ldots,m)$ とすると， $$ Q_j^* = \sqrt{ \frac{2 d K_j}{ h+r c_j } } $$ $$ K_j = K+(c_0-c_j) \theta_j $$ を区分ごとに計算し，最小費用の区分を求めれば良い．
総量割引
総量割引の場合には，発注量 $Q$ が分岐点 $\theta$ を超えると， その購入単価が $c_0$ から $c_1$ にすべて割引されるので，購入費用 $c(Q)$ は， $$ c(Q)= \left\{ \begin{array}{l l} 0 & Q=0 \\ K + c_0 Q & 0 <Q \leq \theta \\ K + c_1 Q & \theta <Q \end{array} \right. $$ と書くことができる．
$c_0 >c_1$ なので， $c(Q)= K + c_1 Q$ のときの最適発注量 $$ Q_1^* =\sqrt{ \frac{2 d K}{ h+r c_1 } } $$ が $\theta$ 以上のときには，$Q_1^*$ が最適になる．
$Q_1^* < \theta$ のときには，$c(Q)= K + c_0 Q$ のときの最適発注量 $$ Q_0^* =\sqrt{ \frac{2 d K}{ h+r c_0 } } $$ と発注量が $\theta$ のときの総費用をくらべて，小さい方を選択すれば良い．
一般に，区分 $j$ に対する購入単価を $c_j (j=0,1,\ldots,m)$ とすると， $$ Q_j^* =\sqrt{ \frac{2 d K}{ h+r c_j } } $$ を区分ごとに計算し，それが区分内にある場合には，$Q_j^*$ を候補とし， 区分外の場合には，分岐点 $\theta_j$ を候補とすれば良い．
経済発注量モデルを解くための関数 eoq
引数：
- K: 発注固定費用
- d: 需要量
- h: 在庫費用（割引モデルの場合には，金利分以外の在庫費用）
- b: 品切れ費用； Noneの場合には，品切れを許さないモデルになる．
- r: 利子率； 割引モデルの場合に使用する．
- c: 割引モデルの場合の区間ごとの割引費用；
- theta: 割引モデルの場合の区間の切れ目（分岐点）； 最初の区分は必ず0と仮定
- discount:割引の種類を表す文字列を入れる． "incremental"の場合は増量割引， "all"の場合は総量割引， None（既定値）の場合には割引がないモデルになる．
返値：
- Q: 最適発注量
- z: 最適値
K = 300.
d = 10.
h = 10.
b = None #40.
r = 0.01
c = [350., 200.]
theta = [0., 30.] #最初の区分は必ず0と仮定
eoq(K, d, h, b, r, c, theta, "incremental")
#eoq(K, d, h, b, r, c, theta)
(Q,R)方策と(s,S)方策のシミュレーションを行う関数 simulate_inventory
$(Q,R)$方策に基づくを在庫シミュレーションをNumPyを用いて行う。手持ち在庫量に発注済みの量を加えたものを在庫ポジションとよぶ． $(Q,R)$方策とは，正味在庫ポジション（在庫ポジション $-$ バックオーダー（品切れ）量）が発注点 $R$ を下回ったら、固定数量 $Q$ だけ発注する方策である．
一方， $(s,S)$方策とは，正味在庫ポジションが発注点 $s$ を下回ったら， 正味在庫ポジションが $S$ になるように発注する方策である．
保証リード時間 $GLT$ と補充リード時間 $RLT$ に対する $(s,S)$ 方策のシミュレーションは，以下の式に基づく．
- 在庫フロー： $I_t = I_{t-1} + d_{t-GLT} + Q_{t-RLT}$
- 正味在庫ポジション： $NI_t = NI_{t-1} -d_{t}$
- 発注量 $Q$: $NI_t < s$ のとき $Q_t = [S-NI_t]^+$
- 正味在庫ポジション （発注後）： $NI_t = NI_{t} + Q_t$
引数：
- n_samples: シミュレーションを行うサンプル数
- n_periods: 期間数
- mu: 需要の平均（需要は正規分布を仮定する．）
- sigma: 需要の標準偏差
- LT: リード時間
- Q: 発注ロットサイズ
- R: 発注点
- z: 安全在庫係数；安全在庫係数は基在庫レベル $S$ を求めるときに使われる。
- b: バックオーダー費用
- h: 在庫費用
- fc: 発注固定費用
- S: $(s,S)$方策のパラメータ； $s$ は $R$ と同じ． $S$ が None（既定値）の場合には， $(Q,R)$ 方策を適用する．
返値：
- cost: 総費用のサンプルごとの期待値を表すNumPy配列で，形状は（n_samples,)
- I: 在庫の時間的な変化を表すNumPy配列で、形状は (n_samples, n_periods+1)；在庫の可視化に用いる。
n_samples =1
n_periods =100
alpha = 3.
fc = 0.
mu = 100.
sigma = 6.
b = 100.
h = 2.
LT = 0
from scipy.stats import norm
omega = b/(b+h)
print("臨界率=", omega)
print("安全在庫係数=", norm.ppf(omega) )
z = norm.ppf(omega)
if LT ==0:
s = mu
S = mu
else:
s, S = approximate_ss(mu, sigma, LT-1, b, h, fc) #LT＝０には対応していない！ LT=1のときが近似公式の0に対応
print(s,S)
# min_Q = Q = np.sqrt(2*fc*mu/h/omega) #EOQ formula
# R = LT*mu + z*sigma*np.sqrt(LT) #z is derived by b and h
# print(Q,R)
cost, I = simulate_inventory(n_samples =n_samples, n_periods =n_periods, mu=mu, sigma=sigma, LT=LT, Q=None, R=s, b=b, h=h, fc=fc, S=S)
#cost.mean(), cost.std()
cost.mean()+alpha*cost.std()
plt.plot(I[0]); #在庫の推移の可視化
n_samples =1
n_periods =100
n_stages = 3
fc = 1.
mu = 100.
sigma = 10.
b = 100.
h = np.array([10., 5., 2.])
LT = np.array([1, 1, 1])
ELT = np.zeros(n_stages,dtype=int)
ELT[0] = LT[1]
for i in range(1,n_stages):
ELT[i] = ELT[i-1] + 1 + LT[i]
print("Echelon LT=", ELT)
s = np.zeros(n_stages)
S = np.zeros(n_stages)
for i in range(n_stages):
s[i], S[i] = approximate_ss(mu, sigma, ELT[i], b, h[i], fc)
print(s, S)
for i in range(1):
cost, I, T = multi_stage_simulate_inventory(n_samples, n_periods, mu, sigma, LT,
s, S, b =100., h = np.array([10., 5., 2.]), fc=fc)
print(cost.mean())
plt.plot(I[0,2]);
(Q,R)方策のパラメータ最適化関数 optimize_qr
$Q$ が与えられていないときには， $Q^*$ は経済発注量モデルを用いて初期化する．
$Q^*$ を固定し， 発注点 $R$ は，リード時間内の安全在庫量を初期値として1つずつ増やし，費用が増加した地点を最適な発注点 $R^*$ とする．
$R^*$ を固定し， 経済発注量の前後を探索し，費用が増加に転じた点を最適な発注量 $Q^*$ とする．
期待値だけでなく費用の標準偏差も考慮し， 「期待値 $+ \alpha$ 標準偏差」 を評価尺度とする．
シミュレーションベースなので，任意の分布，任意の確率過程， 任意の制約に対応できる．
引数：
- n_samples: シミュレーションを行うサンプル数
- n_periods: 期間数
- mu: 需要の平均（需要は正規分布を仮定する．）
- sigma: 需要の標準偏差
- LT: リード時間
- Q: 発注ロットサイズ； Qを固定したいときには整数を入れ，最適化したいときには None を入れる．
- R: 発注点； Rを固定したいときには整数を入れ，最適化したいときには None を入れる．
- z: 安全在庫係数； None（既定値）の場合には，在庫費用とバックオーダー費用から計算される．
- b: バックオーダー費用
- h: 在庫費用
- fc: 発注固定費用
- alpha: リスクの考慮の度合い；「期待値 $+ \alpha$ 標準偏差」 を評価尺度とする．
返値：
- R: 最適化された発注点
- Q: 最適化された発注量
n_samples =100
n_periods =100
alpha = 1.
fc = 10000.
mu = 100.
sigma = 10.
b = 100.
h = 1.
LT = 1
np.random.seed(123)
R,Q = optimize_qr(n_samples =n_samples, n_periods =n_periods, mu=mu, sigma=sigma, LT=LT, Q=None, R=None, z=z, b=b, h=h, fc=fc, alpha=alpha)
#print(Q, R)
#シミュレーションと可視化
cost, I = simulate_inventory(n_samples =n_samples, n_periods =n_periods, mu=mu, sigma=sigma, LT=LT, Q=Q, R=R, b=b, h=h)
plt.plot(I[0]); #在庫の推移の可視化
plt.hist(cost);
(s,S)方策のパラメータ最適化関数 optimize_ss
$Q$ が与えられていないときには， $Q^*$ は経済発注量モデルを用いて初期化する．
$Q^*$ を固定し， 発注点 $R=s$ は，リード時間内の安全在庫量を初期値として1つずつ増やし，費用が増加した地点を最適な発注点 $R^*$ とする．
$R^*$ を固定し， 「$R^* +$経済発注量」の前後を探索し，費用が増加に転じた点を最適な $S^*$ とする．
期待値だけでなく費用の標準偏差も考慮し， 「期待値 $+ \alpha$ 標準偏差」 を評価尺度とする．
シミュレーションベースなので，任意の分布，任意の確率過程， 任意の制約に対応できる．
引数：
- n_samples: シミュレーションを行うサンプル数
- n_periods: 期間数
- mu: 需要の平均（需要は正規分布を仮定する．）
- sigma: 需要の標準偏差
- LT: リード時間
- Q: 発注ロットサイズ； Qを固定したいときには整数を入れ，最適化したいときには None を入れる．
- R: 発注点； Rを固定したいときには整数を入れ，最適化したいときには None を入れる．
- z: 安全在庫係数； None（既定値）の場合には，在庫費用とバックオーダー費用から計算される．
- b: バックオーダー費用
- h: 在庫費用
- fc: 発注固定費用
- alpha: リスクの考慮の度合い；「期待値 $+ \alpha$ 標準偏差」 を評価尺度とする．
- S: 基在庫レベル； 最適化したいときには None（既定値）を入れる．
返値：
- s: 最適化された発注点
- S: 最適化された基在庫レベル
n_samples =100
n_periods =100
alpha = 1.
fc = 1.
mu = 100.
sigma = 10.
b = 100.
h = 10.
LT = 1
np.random.seed(123)
R,S = optimize_ss(n_samples =n_samples, n_periods =n_periods, mu=mu, sigma=sigma, LT=LT, Q=None, R=None, z=z, b=b, h=h, fc=fc, alpha=alpha, S=None)
s, S = approximate_ss(mu, sigma, LT, b, h, fc)
print(s, S)
#シミュレーションと可視化
for i in range(10):
cost, I = simulate_inventory(n_samples =n_samples, n_periods =n_periods, mu=mu, sigma=sigma, LT=LT, Q=None, R=R, b=b, h=h, fc=fc, S=S-i)
print(cost.mean())
plt.plot(I[0]); #在庫の推移の可視化
b,h
plt.hist(cost);
(s,S)方策の近似最適化関数 approximate_ss
回帰を用いてパラメータを探索した研究に基づく，近似 $s,S$ の算出
MANAGEMENT SCIENCE Vol. 30, No. 5, May 1984 Printed in U.S.A.
A REVISION OF THE POWER APPROXIMATION FOR COMPUTING (s, S) POLICIES
RICHARDEHRHARDT AND CHARLESMOSIER
we sets = min{sp,S0}a nd S = min{sp + Dp, S0}, where S0 = μL + Φ−1(p/(p + h)) a
引数：
- mu: 需要の平均（需要は正規分布を仮定する．）
- sigma: 需要の標準偏差
- LT: リード時間（0のときに次の期に到着する；シミュレーションのLT=1に相当）
- b: バックオーダー費用
- h: 在庫費用
- fc: 発注固定費用
返値：
- s: 最適化された発注点
- S: 最適化された基在庫レベル
深層学習を用いた拡張が課題となる．
n_samples =1
n_periods = 100
fc = 0.
mu = 100
sigma = 0
b = 100.
h = 1.
LT = 1
s, S = approximate_ss(mu, sigma, LT-1, b, h, fc)
print(s, S)
#cost, I = simulate_inventory(n_samples =n_samples, n_periods =n_periods, mu=mu, sigma=sigma, LT=LT, Q=None, R=s, b=b, h=h, S=S)
#plt.plot(I[0]); #在庫の推移の可視化
#print(cost.mean())
n_samples =100
n_periods =100
alpha = 0.
fc = 10000.
mu = 100.
sigma = 10.
b = 100.
h = 10.
LT = 3
omega = b/(b+h)
z = norm.ppf(omega)
print("R=", (LT)*mu+ z*sigma*np.sqrt(LT) )
R,Q = optimize_qr(n_samples =n_samples, n_periods =n_periods, mu=mu, sigma=sigma, LT=LT, Q=None, R=None, z=z, b=b, h=h, fc=fc, alpha=alpha)
np.random.seed(123)
cost, I = simulate_inventory(n_samples =n_samples, n_periods =n_periods, mu=mu, sigma=sigma, LT=LT, Q=Q, R=R, b=b, h=h)
print("R=",R, "Q=", Q)
print("RQ-policy cost=", cost.mean())
s,S = optimize_ss(n_samples =n_samples, n_periods =n_periods, mu=mu, sigma=sigma, LT=LT, Q=None, R=None, z=z, b=b, h=h, fc=fc, alpha=alpha, S=None)
np.random.seed(123)
cost, I = simulate_inventory(n_samples =n_samples, n_periods =n_periods, mu=mu, sigma=sigma, LT=LT, Q=None, R=s, b=b, h=h, S=S)
print("s=",s, "S=", S)
print("s,S policy cost=", cost.mean())
#plt.plot(I[0]); #在庫の推移の可視化
s, S = approximate_ss(mu, sigma, LT, b, h, fc)
np.random.seed(123)
cost, I = simulate_inventory(n_samples =n_samples, n_periods =n_periods, mu=mu, sigma=sigma, LT=LT, Q=None, R=s, b=b, h=h, S=S)
print("s=",s, "S=", S)
print("Approximate s,S policy cost=", cost.mean())
plt.plot(I[0]); #在庫の推移の可視化
多段階への拡張
$Q^*$ が決まっている場合： サイクル時間 $CT = Q^*/\mu$ を計算し， $L= \max \{ CT, LT \}$ をリード時間 $L$ として最適な基在庫レベル $S^*$ を定期発注方策最適化で求める． $L$分の平均在庫量 $L \mu$ を減じて安全在庫量を求め， さらにリード時間分の平均需要量 $\mu LT$ を加えたものが発注点 $R^*$ となる．
$Q^*$ を最適化した場合： 定期発注方策最適化で求めた $R^*$ を固定し， 最適な発注量 $Q^*$ を地点ごとに（上の関数optimize_qrを用いて）再最適化し，これを繰り返す．
多段階の $(Q,R)$方策と $(s,S)$ 方策のシミュレータを作り， パラメータの調整を行う．
fc = 10000.
mu = 100.
sigma = 10.
b = 1000.
h = 10.
LT = 3
s, S = approximate_ss(mu, sigma, LT, b, h, fc)
print(s, S)
hlist, slist, Slist = [], [], []
Rlist, Qlist =[],[]
for h in range(1,1000, 1):
hlist.append(h)
s, S = approximate_ss(mu, sigma, LT, b, h, fc)
#s,S = optimize_ss(n_samples =n_samples, n_periods =n_periods, mu=mu, sigma=sigma, LT=LT, Q=None, R=None, z=z, b=b, h=h, fc=fc, alpha=alpha, S=None)
slist.append(s)
Slist.append(S)
omega = b/(b+h)
z = norm.ppf(omega)
R = (LT+1)*mu+ z*sigma*np.sqrt(LT+1)
Rlist.append(R)
Qlist.append(s+np.sqrt(2*fc*mu/h/omega))
plt.plot(hlist,slist);
plt.plot(hlist,Slist);
#plt.plot(hlist,Rlist);
#plt.plot(hlist,Qlist);
# n_periods = 1000
# capacity = 10.
# convergence = 1e-5
# b = 100
# h = 1
# LT = 1
# critical_ratio = b/(b+h)
# S = best_dist.ppf(critical_ratio) + 5
# print(S, best_dist.ppf(0.01))
# demand = best_dist.rvs((n_samples,n_periods)) + 5.
# for iter_ in range(100):
# dC, cost, I = base_stock_simulation_using_dist(n_samples, n_periods, demand, capacity, LT, b, h, S)
# S = S - 0.1*dC
# print(f"{iter_}: {S:.2f} {dC:.3f} {cost:.2f}")
# if dC**2<=convergence:
# break
n_samples =100
n_periods =100
alpha = 1.
fc = 10000.
mu = 100.
sigma = 10.
b = 100.
h = 10.
LT = 3
np.random.seed(123)
from scipy.stats import norm
omega = b/(b+h)
print("臨界率=", omega)
print("安全在庫係数=", norm.ppf(omega) )
z = norm.ppf(omega)
capacity = np.inf
S = mu*LT + z*sigma*np.sqrt(LT)
convergence = 1e-5
print("Initial S=", S)
demand = np.random.normal(mu, sigma, (n_samples, n_periods))
print("t: S dS Cost")
for iter_ in range(100):
dC, cost, I = base_stock_simulation(n_samples, n_periods, demand, capacity, LT, b, h, S)
S = S - 1.*dC
print(f"{iter_}: {S:.2f} {dC:.3f} {cost:.2f}")
if dC**2<=convergence:
break
n_samples = 100
n_periods = 10000
mu = 100
sigma = 10
LT = 3
z = 2.43
b, h = 100, 1
capacity = 1000.
S = mu*LT + z*sigma*np.sqrt(LT)
convergence = 1e-5
print("Initial S=", S)
demand = np.random.normal(mu, sigma, (n_samples, n_periods))
print("t: S dS Cost")
for iter_ in range(100):
dC, cost, I = base_stock_simulation(n_samples, n_periods, demand, capacity, LT, b, h, S)
S = S - 1.*dC
print(f"{iter_}: {S:.2f} {dC:.3f} {cost:.2f}")
if dC**2<=convergence:
break
plt.plot(I[0,:1000]);
多段階在庫システムの定期発注方策に対するシミュレーションと微分値の計算関数 multi_stage_base_stock_simulation
多段階在庫システム対して実装したものが、multi_stage_base_stock_simulationである。
引数：
- n_samples: シミュレーションを行うサンプル数
- n_periods: 期間数
- demand: 需要地点ごとの需要量を保持した辞書。キーは顧客の番号。値は、需要のサンプルを表すNumPyの配列で、形状は (n_samples, n_periods)。
- capacity: 容量（生産量上限）を表すNumPy配列
- LT: 地点毎のリード時間を表すNumPy配列
- b: バックオーダー費用
- h: 在庫費用
- S: 基在庫レベル
返値：
- dC： 各地点における総費用の基準在庫レベルによる微分の期待値
- total_cost: 総費用の期待値
- I: 在庫の時間的な変化を表すNumPy配列で、形状はグラフの点の数をn_stagesとしたとき (n_samples, n_stages, n_periods+1)。
np.random.seed(123)
n_samples = 10
n_stages = 3
n_periods = 1000
mu = 100
sigma = 10
LT = np.array([1, 1, 1])
# 初期エシェロン在庫量は、最終需要地点以外には、定期発注方策の1日分の時間を加えておく。
# ELT = LT.cumsum()
# ELT = ELT + 1
# ELT[0] = ELT[0]-1
# ELT[2]+=1
ELT = np.zeros(n_stages)
ELT[0] = LT[1]
for i in range(1,n_stages):
ELT[i] = ELT[i-1] + 1 + LT[i]
print("Echelon LT=", ELT)
z = 1.33
b = 100
h = np.array([10., 5., 2.])
S = mu*(ELT) + z*sigma*np.sqrt(ELT)
print("S=", S)
capacity = np.array([120., 120., 130.])
demand=np.random.normal(mu, sigma, (n_samples, n_periods))
print("t: Cost |dC| S ")
epsilon=0.01
for iter_ in range(100):
dC, cost, I=multi_stage_base_stock_simulation(n_samples, n_periods, demand, capacity, LT, b, h, S)
sum_=0.
for j in range(n_stages):
sum_ += dC[j]**2
S[j]=S[j] - 1.*dC[j]
print(f"{iter_}: {cost:.2f}, {sum_:.3f}", S)
if sum_ <= epsilon:
break
from plotly.subplots import make_subplots
fig = make_subplots(rows=3, cols=1, subplot_titles= ["Stage " + str(i) for i in range(3)] )
for j in range(3):
trace = go.Scatter(
x= list(range(n_periods)),
y= I[0,j,:],
mode = 'markers + lines',
marker= dict(size= 5,
line= dict(width=1),
opacity= 0.3
),
name ="Stage "+str(j)
)
fig.add_trace(
trace,
row=j+1, col=1
)
fig.update_layout(height=1000,title_text=f"Inventories for all stages", showlegend=False)
plotly.offline.plot(fig)
ネットワーク型定期発注方策に対するシミュレーションと微分値の計算関数 network_base_stock_simulation
一般のネットワークに対して実装したものが、以下のnetwork_base_stock_simulationである。
引数：
- G: 有向グラフ（SCMGraphのインスタンス）
- n_samples: シミュレーションを行うサンプル数
- n_periods: 期間数
- demand: 需要地点ごとの需要量を保持した辞書。キーは顧客の番号。値は、需要のサンプルを表すNumPyの配列で、形状は (n_samples, n_periods)。
- capacity: 容量（生産量上限）を表すNumPy配列
- LT: 地点毎のリード時間を表すNumPy配列
- ELT: 地点毎のエシェロンリード時間を表すNumPy配列
- b: バックオーダー費用
- h: 在庫費用
- S: 基在庫レベル
- phi: 部品展開表における親製品を製造するために必要な子製品の量
- alpha: 上流の在庫地点が下流の在庫地点に在庫を配分する率
返値：
- dC： 各地点における総費用の基準在庫レベルによる微分の期待値
- total_cost: 総費用の期待値
- I: 在庫の時間的な変化を表すNumPy配列で、形状はグラフの点の数をn_stagesとしたとき (n_samples, n_stages, n_periods+1)。
G = SCMGraph()
z = 1.65
b = np.array([0.,100.,100.])
h = np.array([1., 5., 2.])
mu = np.array([300., 200., 100.])
sigma = np.array([12., 10., 15.])
capacity = np.array([1320., 1120., 1230.])
G.add_edges_from([(0, 1), (0, 2)])
LT = np.array([3, 1, 2])
ELT = np.zeros(len(G))
phi = {(0, 1): 1, (0, 2): 1}
alpha = {(0, 1): 0.5, (0, 2): 0.5}
n_samples = 100
n_periods = 1000
demand ={} #需要量
for i in G:
if G.out_degree(i) == 0:
demand[i] = np.random.normal(mu[i], sigma[i], (n_samples, n_periods))
# エシェロンリード時間と初期基在庫レベルの設定
ELT, S = initial_base_stock_level(G, LT, mu, z, sigma)
print("ELT=", ELT)
print("S=", S)
print("t: Cost |dC| S ")
epsilon = 0.01
for iter_ in range(100):
dC, cost, I = network_base_stock_simulation(G, n_samples, n_periods, demand, capacity, LT, ELT, b, h, S, phi, alpha)
sum_ = 0.
#print("dC=",dC)
for j in G:
sum_ += dC[j]**2
S[j] = S[j] - 1.*dC[j]
print(f"{iter_}: {cost:.2f}, {sum_:.5f}", S)
if sum_ <= epsilon:
break
from plotly.subplots import make_subplots
fig = make_subplots(rows=3, cols=1, subplot_titles= ["Stage " + str(i) for i in range(3)] )
for j in range(3):
trace = go.Scatter(
x= list(range(n_periods)),
y= I[0,j,:],
mode = 'markers + lines',
marker= dict(size= 5,
line= dict(width=1),
opacity= 0.3
),
name ="Stage "+str(j)
)
fig.add_trace(
trace,
row=j+1, col=1
)
fig.update_layout(height=1000,title_text=f"Inventories for all stages", showlegend=False)
plotly.offline.plot(fig);
G = SCMGraph()
G.add_edges_from([(0, 1), (0, 2)])
z = 1.65
h = np.array([1., 5., 2.])
mu = np.array([300., 200., 100.])
sigma = np.array([12., 10., 15.])
ProcTime = np.array([5, 5 ,5], int)
LTLB = np.array([0,0,0], int)
LTUB = np.array([10,1,2], int)
Lmax, MaxDemand, SafetyCost = max_demand_compute(G, ProcTime, LTLB, LTUB, z, mu, sigma, h)
print("Lmax=",Lmax)
print("MaxDemand=",MaxDemand)
print("SafetyCost=",SafetyCost)
G =SCMGraph()
G.add_edges_from([(0, 1), (0, 2)])
z = 1.65
h = np.array([1., 5., 2.])
mu = np.array([300., 200., 100.])
sigma = np.array([12., 10., 15.])
ProcTime = np.array([5, 5 ,5], int)
LTLB = np.array([0,0,0], int)
LTUB = np.array([10,1,2], int)
total_cost, Lstar, NRT = dynamic_programming_for_SSA(G, ProcTime, LTLB, LTUB, z, mu, sigma, h)
print("TotalCost=",total_cost)
print("Lstar",Lstar)
print("NRT",NRT)
G = SCMGraph()
n = 10
edges =[(i,i-1) for i in range(n-1,0,-1)]
G.add_edges_from( edges )
z = 1.65
h = np.arange(n,0,-1)
mu = np.array([100.]*n)
sigma = np.array([10.]*n)
ProcTime = np.array([5]*n, int)
LTLB = np.zeros(n, int)
LTUB = np.array([0]+[0]*(n-1), int)
total_cost, Lstar, NRT = dynamic_programming_for_SSA(G, ProcTime, LTLB, LTUB, z, mu, sigma, h)
print("TotalCost=",total_cost)
print("Lstar",Lstar)
print("NRT",NRT)
用語と記号：
- 保証リード時間 (guaranteer lead time, service time): $L$
- 生産時間 (production time): $T$
- 入庫リード時間 (incoming lead time): $LI$
- 補充リード時間 (replenishment lead time)
- 正味補充時間 (net replenishment time): $NRT$
- 最大補充リード時間： $MaxLI$
- 最小保証リード時間： $MinLT$
安全在庫配置モデルに対するタブー探索アルゴリズム tabu_search_for_SSA
引数：
- G: 閉路を含まない有向グラフ
- ProcTime: 処理時間（NumPyの配列）；以下の配列はグラフGの点を列挙したときの順序（グラフに点を追加した順序）の順に保管してあるものとする．
- LTUB: 保証リード時間上限（NumPyの配列）
- z: 安全在庫係数
- mu: 需要の平均（NumPyの配列）
- sigma: 需要の標準偏差（NumPyの配列）
- h: 在庫費用（NumPyの配列）
- max_iter: 最大反復数
- TLLB : タブー長の下限の初期値
- TLUB : タブー長の上限の初期値
- seed : 乱数の種
返値：
- best_cost：最良値
- best_sol: 最良解（在庫を保持するなら1、否なら0の配列）
- best_NRT: 最良解における正味補充時間（在庫日数）
- best_MaxLI: 最良解における最大補充リード時間
- best_MinLT: 最良解における最小保証リード時間
G = SCMGraph()
G.add_edges_from([(0, 1), (0, 2)])
z = 1.65
h = np.array([1., 1., 1.])
mu = np.array([300., 200., 100.])
sigma = np.array([12., 10., 15.])
LTUB = np.zeros(3)
ProcTime = np.array([5, 4 ,3], int) # 動的最適化の例題と合わせるため、生産時間から保証リード時間上限を減じてある
best_cost, best_sol, best_NRT, best_MaxLI, best_MinLT = tabu_search_for_SSA(G, ProcTime, LTUB, z, mu, sigma, h, max_iter = 10, TLLB =1, TLUB =3, seed = 1)
print("最良値", best_cost)
print("最良解", best_sol)
print("正味補充時間", best_NRT)
print("最大補充リード時間", best_MaxLI)
print("最小保証リード時間", best_MinLT)
G = SCMGraph()
n = 4
# 点の順序を0,1,2... とするため
for i in range(n):
G.add_node(i)
edges =[(i,i-1) for i in range(n-1,0,-1)]
G.add_edges_from( edges )
z = np.ones(n)
h = np.array([40,30,20,10])
mu = np.array([100.]*n)
sigma = np.array([10.]*n)
ProcTime = np.array([1,1,2,3], int)
LTUB = np.array([0,0,0,0],int)
best_cost, best_sol, best_NRT, best_MaxLI, best_MinLT = tabu_search_for_SSA(G, ProcTime, LTUB, z, mu, sigma, h, max_iter = 10, TLLB =1, TLUB =3, seed = 2)
print("最良値", best_cost)
print("最良解", best_sol)
print("正味補充時間", best_NRT)
print("最大補充リード時間", best_MaxLI)
print("最小保証リード時間", best_MinLT)
pos ={i:(n-i,0) for i in range(n) }
stage_df, bom_df = make_df_for_SSA(G, pos, ProcTime, LTUB, z, mu, sigma, h, best_NRT, best_MaxLI, best_MinLT)
stage_df.to_csv(folder_bom + "stage_ex0.csv")
bom_df.to_csv(folder_bom + "bom_e01.csv")
fig, G, pos = draw_graph_for_SSA_from_df(stage_df, bom_df)
#fig.show()
nx.draw(G, pos=pos)
# n_periods = 100, seed = 1, lr_find=False, max_lr =0.1, moms=(0.85,0.95))
#plotly.offline.plot(fig);
G = SCMGraph()
n =8
# 点の順序を0,1,2... とするため
for i in range(n):
G.add_node(i)
edges =[(i,i-1) for i in range(n-1,0,-1)]
G.add_edges_from( edges )
z = 1.65
#h = np.arange(n,0,-1)
h = np.array([20,20,10,10,10,5,5,1])
mu = np.array([100.]*n)
sigma = np.array([10.]*n)
ProcTime = np.array([5]*n, int)
LTUB = np.zeros(n)
best_cost, best_sol, best_NRT, best_MaxLI, best_MinLT = tabu_search_for_SSA(G, ProcTime, LTUB, z, mu, sigma, h, max_iter = 10, TLLB =1, TLUB =3, seed = 1)
print("最良値", best_cost)
print("最良解", best_sol)
print("正味補充時間", best_NRT)
print("最大補充リード時間", best_MaxLI)
print("最小保証リード時間", best_MinLT)
n=7
G = SCMGraph()
for i in range(n):
G.add_node(i)
G.add_edges_from([(0, 2), (1, 2), (2,4), (3,4), (4,5), (4,6)])
z = np.full(len(G),1.65)
h = np.array([1,1,3,1,5,6,6])
mu = np.array([200,200,200,200,200,100,100])
sigma = np.array([14.1,14.1,14.1,14.1,14.1,10,10])
LTUB = np.array([0,0,0,0,0,3,1], int)
ProcTime = np.array([6,2,3,3,3,3,3], int) # 動的最適化の例題と合わせるため、生産時間から保証リード時間上限を減じてある
best_cost, best_sol, best_NRT, best_MaxLI, best_MinLT = tabu_search_for_SSA(G, ProcTime, LTUB, z, mu, sigma, h, max_iter = 10, TLLB =1, TLUB =3, seed = 1)
print("最良値", best_cost)
print("最良解", best_sol)
print("正味補充時間", best_NRT)
print("最大補充リード時間", best_MaxLI)
print("最小保証リード時間", best_MinLT)
# pos = G.layout()
# stage_df, bom_df = make_df_for_SSA(G, pos, ProcTime, LTUB, z, mu, sigma, h, best_NRT, best_MaxLI, best_MinLT)
# stage_df.to_csv(folder_bom + "stage_ex1.csv")
# bom_df.to_csv(folder_bom + "bom_ex1.csv")
# fig, G, pos = draw_graph_for_SSA_from_df(stage_df, bom_df)
# #fig.show()
# nx.draw(G, pos=pos)
stage_df = pd.read_csv(folder_bom + "stage_ex1.csv")
bom_df = pd.read_csv(folder_bom + "bom_ex1.csv")
best_cost, stage_df, bom_df, fig = solve_SSA(stage_df, bom_df)
stage_df
cost, stage_df, I, cost_list, phi_list = periodic_inv_opt_fit_one_cycle(stage_df, bom_df, max_iter = 1000, n_samples = 10,
n_periods = 100, seed = 1, lr_find=False, max_lr =6.0, moms=(0.85,0.95))
fig = plot_simulation(stage_df, I, n_periods =100, samples=5)
plotly.offline.plot(fig);
fig = plot_inv_opt_lr_find(cost_list, phi_list)
plotly.offline.plot(fig);
Willemsのベンチマーク問題例の読み込み
Willemsのデータ（XMLファイル）を読み込みグラフのデータを組み立てる。
引数：
- file_name: Willemsのデータ番号；"01" から "38" までの「文字列」を入れる。
返値：
- G: SCMGraphのインスタンス；閉路を含まない有向グラフ
- pos: 点の座標；点の名前をキー、(x,y)座標を値としてもつ辞書
グラフは、点の属性として以下の値を保管してある。
- stageName: 点の名称
- stageCost: 点で付加される費用（浮動小数点数の文字列）
- relDepth: 点の深さ（下流（需要側）から0,1,2,...)
- stageClassification: 点の種類
- avgDemand: 平均需要量（整数の文字列）
- stDevDemand: 需要の標準偏差（整数の文字列）
- maxServiceTime: 最大サービス時間（保証リード時間）（整数の文字列）
- serviceLevel: サービス率（浮動小数点数の文字列）
- stageTime: 生産時間（整数の文字列）
- xPos: $x$ 座標（整数の文字列）
- yPos: $y$ 座標（整数の文字列）
G, pos = read_willems(file_name="../data/bom/03.xml")
nx.draw(G, pos=pos, with_labels=True, node_color="Yellow")
net = pyvis.network.Network(height='500px', width='500px', directed = True,
notebook = True, bgcolor='#ffffff', font_color=False, layout=1, heading='')
net.from_nx(G)
#net.show("example.html")
点の属性は、以下のように点に付随する属性の辞書、もしくはget_node_attributesメソッドを用いてアクセスできる。
for i in G: print( G.nodes[i]["stageTime"], end=" " )
print( nx.get_node_attributes(G, "maxServiceTime") )
在庫費用、需要量、需要の標準偏差の計算
ベンチマークデータからタブー探索で必要なデータを抽出し、必要なデータを準備する。 各地点には付加された費用と生産時間、需要地点には需要と標準偏差と保証リー時間上限（サービス時間）が入っている。 すべての点に対して、在庫費用と需要・標準偏差など、安全在庫配置モデルの計算に必要なデータを計算する関数を準備しておく。
引数：
- G: ベンチマーク問題例のグラフ
返値：
- ProcTime: 処理時間 （NumPy配列）
- LTUB: 保証リード時間上限 （NumPy配列）
- z: 安全在庫係数 （NumPy配列）
- mu: 需要の平均 （NumPy配列）
- sigma: 需要の標準偏差 （NumPy配列）
- h: 在庫費用 （NumPy配列）
mu, sigma, h, z, ProcTime, LTUB = extract_data_for_SSA(G)
best_cost, best_sol, best_NRT, best_MaxLI, best_MinLT = \
tabu_search_for_SSA(G, ProcTime, LTUB, z, mu, sigma, h, max_iter = len(G)*2,
TLLB =int(np.sqrt(len(G))/2), TLUB =int(np.sqrt(len(G))), seed = 1)
print("最良値", best_cost)
print("最良解", best_sol)
print("正味補充時間", best_NRT)
print("最大補充リード時間", best_MaxLI)
print("最小保証リード時間", best_MinLT)
nx.draw(G, pos=pos, with_labels=True, node_size=best_NRT*50,node_color=best_sol)
print("best obj.=", best_cost)
print("best sol.=", best_sol)
print("best NRT=", best_NRT)
fig = draw_graph_for_SSA(G, pos, best_NRT, best_MaxLI, best_MinLT)
#plotly.offline.plot(fig); #htmlを吐きたいときには、ここを生かす。
#fig.show()
安全在庫配置問題の結果をデータフレームに変換する関数 make_df_for_SSA
点(stage)の情報を保管したデータフレームと枝（部品展開表；bom）の情報を保管したデータフレームを構築する。
引数：
- G: 閉路を含まない有向グラフ
- pos: 点の座標を保管した辞書
- ProcTime: 処理時間 （NumPy配列）
- LTUB: 保証リード時間上限 （NumPy配列）
- z: 安全在庫係数 （NumPy配列）
- mu: 需要の平均 （NumPy配列）
- sigma: 需要の標準偏差 （NumPy配列）
- h: 在庫費用 （NumPy配列）
- best_NRT: 最良解における正味補充時間（在庫日数）
- best_MaxLI: 最良解における最大補充リード時間
- best_MinLT: 最良解における最小保証リード時間
返値：
点(stage)の情報を保管したデータフレーム
- name: 点名
- net_replenishment_time: 正味補充時間
- max_guaranteed_LT: 保証リード時間上限
- processing_time: 処理時間
- replenishment_LT: 補充リード時間
- guaranteed_LT: 保証リード時間
- z: 安全在庫係数
- average_demand: 需要の平均
- sigma: 需要の標準偏差
- h: 在庫費用
- b: 品切れ費用（在庫費用hとサービス率pをもとにp=b/(b+h)になるように計算）
- capacity: 生産容量（需要量の平均の10倍と設定しておく．）
- x: 点のx座標
- y: 点のy座標
枝（部品展開表；bom）の情報を保管したデータフレーム
- child: 子製品名
- parent: 親製品名
- units: 親製品を製造するのに必要な子製品の数（すべて1に設定）
- allocation: 在庫を複数の下流の点に配分するときの割合（在庫シミュレーション最適化で用いる．）
stage_df, bom_df = make_df_for_SSA(G, pos, ProcTime, LTUB, z, mu, sigma, h, best_NRT, best_MaxLI, best_MinLT)
#stage_df.to_csv(folder_bom + "ssa02.csv")
stage_df
bom_df.head()
fig, G, pos = draw_graph_for_SSA_from_df(stage_df, bom_df)
#plotly.offline.plot(fig)
nx.draw(G, pos=pos)
for i in range(1, 39):
fn = str(i)
if len(fn) == 1:
fn = "0"+fn
#print(i)
G, pos = read_willems(file_name=fn)
mu, sigma, h, z, ProcTime, LTUB = extract_data_for_SSA(G)
n = len(G)
best_NRT = np.zeros(n)
best_MaxLI = np.zeros(n)
best_MinLT =np.zeros(n)
stage_df, bom_df = make_df_for_SSA(
G, pos, ProcTime, LTUB, z, mu, sigma, h, best_NRT, best_MaxLI, best_MinLT)
stage_df.to_csv("../data/bom/ssa"+fn+".csv")
bom_df.to_csv("../data/bom/ssa_bom"+fn+".csv")
stage_df = pd.read_csv(folder_bom + "ssa03.csv")
bom_df = pd.read_csv(folder_bom + "ssa_bom03.csv")
best_cost, stage_df, bom_df, fig = solve_SSA(stage_df, bom_df)
stage_df
データフレームをもとに定期発注方策最適化を行う関数 periodic_inv_opt
定期発注方策の仮定の下では，保証リード時間は $0$ とみなす必要位がある． また，各点のリード時間は，安全在庫配置モデルの補充リード時間 (replenishment leadtime: best_MaxLI) に （定期発注方策なので） $1$ を加えたもの相当する．
ただし出荷を待っている（すでに下流からの注文がある；保証リード時間に相当する）在庫も含めた在庫レベルを考えるものとする．
保証リード時間中の需要は確定値なので考慮する必要はない． 保証リード時間 $+1$ 日をリード時間と考える．
最適な $(s,S)$ を近似的に求め，シミュレーションを行う．
安全在庫配置の結果がすでに入ったデータフレームをもとに，定期発注方策の最適化を行う．
需要が負になることが頻繁に発生する場合には，基在庫レベルを多少大きくする必要があり，どの程度大きくすれば良いかは実験を行う必要がある．
stage_df = pd.read_csv(folder_bom + "ssa02.csv", index_col=0)
bom_df = pd.read_csv(folder_bom + "ssa_bom02.csv", index_col=0)
best_cost, stage_df, bom_df, fig = solve_SSA(stage_df, bom_df)
cost, stage_df, I = periodic_inv_opt(stage_df, bom_df, max_iter = 1000, n_samples = 10, n_periods = 100, seed = 1)
stage_df = pd.read_csv(folder_bom + "ssa01.csv", index_col=0)
bom_df = pd.read_csv(folder_bom + "ssa_bom01.csv", index_col=0)
best_cost, stage_df, bom_df, fig = solve_SSA(stage_df, bom_df)
cost, stage_df, I = periodic_inv_opt(stage_df, bom_df, max_iter = 1000, n_samples = 10, n_periods = 100, seed = 1)
cost, stage_df, I, cost_list, phi_list = periodic_inv_opt_fit_one_cycle(stage_df, bom_df, max_iter = 1000, n_samples = 10,
n_periods = 100, seed = 1, lr_find=False, max_lr =3., moms=(0.85,0.95))
fig = plot_simulation(stage_df, I, n_periods =100, samples=5)
plotly.offline.plot(fig);
fig = plot_simulation(stage_df, I, n_periods =100, samples=5)
plotly.offline.plot(fig);
