サービスレベルと安全在庫係数
ここでは，小売店に顧客がやってくることによって需要が発生する場合を考える． 顧客は，毎日小売店にやってきて商品を購入するが，その量は気分によって変わるものとする． 一般には，商品が品切れしてしまうことは避けなければならない． 商品を余分に在庫しておく費用が，品切れしてしまう費用より小さい場合には， 商品を余分に在庫して気まぐれな顧客のために準備しておくことが常套手段である． このように，需要の不確実性に対処するためにもつ在庫を安全在庫（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)