在庫最適化と安全在庫配置システム MESSA (MEta Safety Stock Allocation system)

Streamlitで作成した簡易システムの紹介ビデオ

安全在庫配置システム: MESSA (MEta Safety Stock Allocation system)

はじめに

需要の不確実性に対処するための在庫が安全在庫である. 1拠点の安全在庫は古典であり,簡単な「公式」もしくはシミュレーションによって最適化できる. ここでは,サプライ・チェイン全体を考えて,各在庫地点の安全在庫を決定するための方策のパラメータの最適化を考える.

多くの研究者の長年の努力にも関わらず,多拠点での一般ネットワーク型の在庫最適化モデルはいまだ未解決の部分も多い. ここでは,適切に付加した仮定の下で,幾つかの最適化モデルを構築する.

基礎理論

サービスレベルと安全在庫係数

ここでは,小売店に顧客がやってくることによって需要が発生する場合を考える. 顧客は,毎日小売店にやってきて商品を購入するが,その量は気分によって変わるものとする. 一般には,商品が品切れしてしまうことは避けなければならない. 商品を余分に在庫しておく費用が,品切れしてしまう費用より小さい場合には, 商品を余分に在庫して気まぐれな顧客のために準備しておくことが常套手段である. このように,需要の不確実性に対処するためにもつ在庫を安全在庫(safety inventory, safety stock)とよぶ.

いま,簡単のため,毎日の需要量は,定常な分布にしたがっているものとする. ここで「定常」とは,日によって需要の分布が変化しないことを指す. 分布によっては,負の需要量になる可能性があるが,その場合には需要量を $0$ と定義する. たとえば, 需要が正規分布と仮定した場合には,負の部分を除いた正規分布(切断正規分布)を指すものとする.

Excelインターフェイス

在庫最適化用のExcelテンプレートを生成する関数 make_excel_messa

返値:

  • wb: 在庫最適化用のExcel Workbook

make_excel_messa[source]

make_excel_messa()

make_excel_messa関数の使用例

wb = make_excel_messa()
wb.save("messa-templete.xlsx")

Excel Workbookからデータフレームを生成する関数 prepare_df_for_messa

引数:

  • wb: Excel Workbook

返値:

  • stsge_df : 点の情報を保管したデータフレーム
  • bom_df : 枝の情報を保管したデータフレーム

prepare_df_for_messa[source]

prepare_df_for_messa(wb)

Excel Workbookから最適化用のデータを準備する関数 prepare_opt_for_messa

引数:

  • wb: Excel Workbook

返値:

  • G: 閉路を含まない有向グラフ
  • ProcTime: 処理時間(NumPyの配列);以下の配列はグラフGの点を列挙したときの順序(グラフに点を追加した順序)の順に保管してあるものとする.
  • LTUB: 保証リード時間上限(NumPyの配列)
  • z: 安全在庫係数
  • mu: 需要の平均(NumPyの配列)
  • sigma: 需要の標準偏差(NumPyの配列)
  • h: 在庫費用(NumPyの配列)

prepare_opt_for_messa[source]

prepare_opt_for_messa(wb)

prepare_opt_for_messaの使用例

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)

安全在庫配置の最適化し,Excelシートを生成する関数 messa_for_excel

引数:

  • wb: Excel Workbook
  • best_NRT: 最良解における正味補充時間(在庫日数)
  • best_MaxLI: 最良解における最大補充リード時間
  • best_MinLT: 最良解における最小保証リード時間

返値:

  • wb: Excel Workbook

messa_for_excel[source]

messa_for_excel(wb, best_NRT, best_MaxLI, best_MinLT)

messa_for_excel関数の使用例

wb = load_workbook("messa-ex0.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-ex0-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) 
381.5586541289441 795.0168907487716
print(cost)
[ 919.89862591  948.05684698  920.90360513  901.53117422  978.40281571
  985.80314351  982.80524091 1004.41837426  908.45463009  910.88346913
 1013.30517203 1002.63451923  944.20356419  950.92319381  991.4487325
 1005.74315039  989.71830694  894.95722159  932.32456478 1034.34868868
  989.93721087  887.14958143 1005.24264336  930.87790909  985.53006173
  964.73001228  867.6148398   889.51345778  927.08272682  989.39805795
  947.09357108  944.11876153  851.58416805  912.99713672  918.22448857
  920.1760374   938.63267923 1004.94235121  934.5637913  1028.30223003
  899.44690617  987.70898023  940.43204155  967.66536805  973.9129936
  955.72520684  992.98751856  936.96843298  934.65462183  945.06470443
  981.24517606  900.81826348  988.70784619  921.04706172  936.89958353
  947.90388996  989.12432527  872.83433521  932.77813591  862.17094196
  985.72413895  936.89237018 1023.71661428  964.87732328  953.76558634
 1011.30512442  990.32441593  955.50739932  994.94277939  932.80300147
  877.27644321 1025.84460488  923.69863305  936.88988765  891.83696978
  949.58100957 1004.1200027   900.381117   1025.9255643   962.63342321
  967.47201364  947.00494303  937.47668896  915.74201878  919.43516852
  957.93705096 1076.91097663  970.81206875  995.17904569  974.26974145
  996.66401229  885.12782484  926.88501749  928.31426422  919.09466732
  916.24530175  973.52307689  943.04661681  901.73595883  868.15223943]

Pyvizの図を生成する関数 generate_pyvis_for_messa

引数:

  • wb: Excel Workbook
  • G: 閉路を含まない有向グラフ
  • best_cost:最良値
  • best_sol: 最良解(在庫を保持するなら1、否なら0の配列)
  • best_MaxLI: 最良解における最大補充リード時間
  • best_MinLT: 最良解における最小保証リード時間
  • notebook: Jupyter用の図を生成するとき True (既定値 True)

返値:

  • net: pyvisの図オブジェクト

generate_pyvis_for_messa[source]

generate_pyvis_for_messa(wb, G, best_sol=None, best_MaxLI=None, best_MinLT=None, notebook=True)

generate_pyvis_for_messa関数の使用例

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: 各期の発注量を入れたリスト

ww[source]

ww(demand, fc=100.0, vc=0.0, h=5.0)

ww関数の使用例

fixed= 3
variable= [1,1,3,3,3] 
h = 1. 
demand=[5,7,3,6,4]
ww(demand, fc =fixed, vc=variable, h=h)
(57.0, array([ 5., 16.,  0.,  0.,  4.]))

分布に最もよく適合した連続確率分布を求める関数 best_distribution

定常分布をもつ需要データに最もよく合う連続確率分布を求める.

引数:

  • data: データの配列

返値:

  • fig: データの度数分布表(ヒストグラム)と最良誤差分布の密度関数の図(適合する分布がない場合には,度数分布表のみ)
  • frozon_dist: 最良誤差分布(パラメータ固定)
  • best_fit_name: 最良誤差分布の名前
  • best_fit_params: 最良誤差のパラメータ

best_distribution[source]

best_distribution(data)

best_distribution関数の適用例

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)
powerlognorm (11.799333308942117, 0.20671017783771561, 16.468598592930654, 116.65707079233306)

分布に適合した表形式の確率分布を求める関数 best_histogram

定常分布をもつ需要データから表形式(度数分布)による確率分布を求める. データ数が少ない(50未満)の場合には,平均値を合わせるため左に0.5だけ$x$軸の区分をシフトさせる.

引数:

  • data: データの配列
  • nbins: ビン数(既定値は50とデータ数の小さい方)

返値:

  • fig: データの度数分布表(ヒストグラム)と分布の密度関数の図
  • hist_dist: 度数分布(パラメータ固定)

best_histogram[source]

best_histogram(data, nbins=50)

best_histogram関数の適用例

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()
(100.48413228614606, 100.20054738268533)

得られた分布に対して,在庫管理に必要な統計量を計算する.

print( best_dist.ppf(0.9), hist_dist.ppf(0.9) ) #90%品切れしないための在庫量
print( best_dist.sf(200), hist_dist.sf(200) ) #在庫が200のときの品切れ率
112.80769882762985 112.77535790405301
0.0 0.0

需要分布に対して最適な適合をもつ分布を返す関数 fit_demand

引数:

  • demand_df: 需要データフレーム
  • cust_selected: 需要の分布を求めたい顧客名のリスト
  • prod_selected: 需要の分布を求めたい製品名のリスト
  • agg_period: 集約を行う期。例えば1週間を1期とするなら"1w"とする。既定値は1日("1d")
  • nbins: ヒストグラムのビン数(既定値は50)
  • L: リード時間(在庫管理で用いる.L日前からL-1前の需要の合計の分布をとる.)

返値:

  • dic: 顧客名と製品名のタプルをキーとし,分布の情報を値とした辞書; 分布の情報は以下のタプル

    • fig: 誤差が最小の分布の密度関数と度数分布表の図
    • fig2: 近似的な表形式の確率関数と度数分布表の図
    • best_dist: 誤差最小の分布の密度関数
    • hist_dist: 度風分布表のデータから構成された確率関数
    • best_fit_name: 誤差最小の分布名
    • best_fit_params: 誤差を最小にする分布のパラメータ

fit_demand[source]

fit_demand(demand_df, cust_selected, prod_selected, agg_period='1d', nbins=50, L=1)

需要分布の推定(顧客と製品に対して):リード時間Lを入力とし,その間の需要の分布を計算する必要がある!

fit_demand関数の適用例

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)
/Users/mikiokubo/Library/Caches/pypoetry/virtualenvs/scmopt-KVX_UVCf-py3.8/lib/python3.8/site-packages/scipy/stats/_continuous_distns.py:4024: RuntimeWarning:

divide by zero encountered in log

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);
johnsonsu (3.3868462264434953, 0.24284658331908165, 97.00000053369794, 5.641728376788476e-07)

経済発注量モデル

EOQ公式(economic ordering quantity formula)

  • 需要量は1日あたり $d$ 単位である.
  • 発注の際には,発注量によらない固定的な費用(これを発注費用とよぶ) $K$ 円が課せられる.
  • 在庫保管費用は保管されている在庫量に比例してかかり, 品目 $1$ 個あたりの保管費用は $1$日で $h$ 円とする.

最適発注量: $$ Q^* = \sqrt{ \frac{2Kd}{h} } $$

最適値: $$ \sqrt{ 2 K h d} $$

バックオーダーを考慮

  • 品切れ費用はバックオーダーの量に比例してかかり, 品目 $1$ 個あたりの品切れ費用は $1$ 日あたり $b$ 円とする.

臨界率(critical ratio): $$ \omega = \frac{b}{b+h} $$

最適値: $$ \sqrt{ 2 K h d \omega} $$

最適発注量: $$ Q^* =\sqrt{ \frac{2Kd}{h \omega} } $$

数量割引の考慮

  • 増分割引(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: 最適値

eoq[source]

eoq(K, d, h, b, r, c, theta, discount=None)

eoq関数の使用例

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)
(89.44271909999159, 3073.3126291998988)

(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);在庫の可視化に用いる。

simulate_inventory[source]

simulate_inventory(n_samples=1, n_periods=10, mu=100.0, sigma=10.0, LT=3, Q=300.0, R=100.0, b=100.0, h=1.0, fc=10000.0, S=None)

(Q,R)方策と(s,S)方策による在庫シミュレーション

simulate_inventory関数の使用例

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=0には対応していない! 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()
臨界率= 0.9803921568627451
安全在庫係数= 2.0619165008094615
100.0 100.0
0.0
plt.plot(I[0]); #在庫の推移の可視化

多段階の場合のシミュレーション

multi_stage_simulate_inventory[source]

multi_stage_simulate_inventory()

多段階ネットワークに対する(s,S)方策による在庫シミュレーション

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: 最適化された発注量

optimize_qr[source]

optimize_qr(n_samples=1, n_periods=10, mu=100.0, sigma=10.0, LT=3, Q=None, R=None, z=None, b=100.0, h=1.0, fc=10000.0, alpha=1.0)

(Q,R)方策のパラメータ最適化

optimize_qr関数の使用例

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: 最適化された基在庫レベル

optimize_ss[source]

optimize_ss(n_samples=1, n_periods=10, mu=100.0, sigma=10.0, LT=3, Q=None, R=None, z=None, b=100.0, h=1.0, fc=10000.0, alpha=1.0, S=None)

(s,S)方策のパラメータ最適化

optimize_ss関数の使用例

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]); #在庫の推移の可視化
Use base stock policy!
218.88226462598004 218.88226462598004
268.83869384305
267.09802453439505
266.49604253145344
271.79213720384644
275.20204386181666
286.2810742804162
295.2959839051329
299.99477028405556
319.5278164677335
339.65342633672446
b,h
(100.0, 10.0)
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: 最適化された基在庫レベル

深層学習を用いた拡張が課題となる.

approximate_ss[source]

approximate_ss(mu=100.0, sigma=10.0, LT=0, b=100.0, h=1.0, fc=10000.0)

approximate_ss関数の使用例

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())
99.13000106299978 99.13000106299978
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, z=z, 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, z=z, 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, z=z, b=b, h=h, S=S) 
print("s=",s, "S=", S)
print("Approximate s,S policy cost=", cost.mean())
plt.plot(I[0]); #在庫の推移の可視化
R= 323.12595676092786
R= 324 Q= 464
RQ-policy cost= 4488.511312454572
s= 333 S= 792
s,S policy cost= 4313.729385300712
s= 349.5565241426376 S= 768.2721993964062
Approximate s,S policy cost= 4406.340519095063

多段階への拡張

$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)$ 方策のシミュレータを作り, パラメータの調整を行う.

(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)
398.3997427604041 817.1154180141727
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);

定期発注方策(単一段階モデル)

base_stock_simulation[source]

base_stock_simulation(n_samples, n_periods, demand, capacity, LT, b, h, S)

単一段階在庫システムの定期発注方策に対するシミュレーションと微分値の計算

base_stock_simulation_using_dist[source]

base_stock_simulation_using_dist(n_samples, n_periods, demand, capacity, LT, b, h, S)

単一段階在庫システムの定期発注方策に対するシミュレーションと微分値の計算

# 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

base_stock_simulation関数を利用した在庫最適化

上で作成した関数を用いて、最適化を行う。base_stock_simulationは、微分値を返すので、最急降下法を適用して最適値を探索する。 微分値の2乗がepsilonを下回ったら終了する。

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
臨界率= 0.9090909090909091
安全在庫係数= 1.3351777361189363
Initial S= 323.12595676092786
t:   S      dS     Cost
0: 323.56 -0.439 377.45
1: 323.51 0.056 377.35
2: 323.52 -0.010 377.35
3: 323.53 -0.010 377.35
4: 323.52 0.012 377.35
5: 323.53 -0.010 377.35
6: 323.53 0.001 377.35
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
Initial S= 342.08883462392373
t:   S      dS     Cost
0: 341.86 0.228 46.55
1: 341.66 0.200 46.50
2: 341.49 0.175 46.47
3: 341.33 0.153 46.44
4: 341.20 0.133 46.42
5: 341.09 0.114 46.40
6: 340.99 0.098 46.39
7: 340.90 0.086 46.38
8: 340.83 0.073 46.37
9: 340.77 0.062 46.37
10: 340.71 0.053 46.36
11: 340.67 0.046 46.36
12: 340.63 0.039 46.36
13: 340.59 0.033 46.36
14: 340.57 0.029 46.36
15: 340.54 0.024 46.36
16: 340.52 0.021 46.35
17: 340.50 0.018 46.35
18: 340.49 0.015 46.35
19: 340.47 0.012 46.35
20: 340.46 0.011 46.35
21: 340.46 0.008 46.35
22: 340.45 0.007 46.35
23: 340.44 0.006 46.35
24: 340.44 0.005 46.35
25: 340.43 0.004 46.35
26: 340.43 0.004 46.35
27: 340.43 0.003 46.35
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)。

multi_stage_base_stock_simulation[source]

multi_stage_base_stock_simulation(n_samples, n_periods, demand, capacity, LT, b, h, S)

多段階在庫システムの定期発注方策に対するシミュレーションと微分値の計算

multi_stage_base_stock_simulation関数を利用した在庫最適化

上で作成した関数を用いて、最適化を行う。multi_stage_base_stock_simulationは、微分値を返すので、最急降下法を適用して最適値を探索する。 微分値の2乗和がepsilonを下回ったら終了する。

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
Echelon LT= [1. 3. 5.]
S= [113.3        323.03627574 529.7397041 ]
t: Cost    |dC|      S 
0: 2757.29, 92.626 [112.3865     326.23927574 538.7692041 ]
1: 2733.44, 27.943 [113.318      329.28117574 542.9908041 ]
2: 2724.89, 9.334 [113.6215     331.12797574 545.4055041 ]
3: 2724.59, 4.356 [114.0775     332.28437574 547.0821041 ]
4: 2724.66, 2.075 [114.235      333.08457574 548.2694041 ]
5: 2725.86, 1.013 [114.33       333.64117574 549.1028041 ]
6: 2726.98, 0.687 [114.3505     334.09137574 549.7981041 ]
7: 2728.21, 0.347 [114.434      334.44877574 550.2592041 ]
8: 2728.95, 0.186 [114.5095     334.67267574 550.6198041 ]
9: 2729.49, 0.102 [114.531      334.81887574 550.9031041 ]
10: 2730.01, 0.085 [114.5225     334.95277574 551.1627041 ]
11: 2730.59, 0.049 [114.592      335.01267574 551.3633041 ]
12: 2730.75, 0.041 [114.6175     335.05577574 551.5587041 ]
13: 2731.03, 0.034 [114.6415     335.06957574 551.7409041 ]
14: 2731.26, 0.029 [114.6415     335.07687574 551.9096041 ]
15: 2731.54, 0.021 [114.6275     335.10497574 552.0495041 ]
16: 2731.84, 0.009 [114.6635     335.13047574 552.1310041 ]
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)
'temp-plot.html'