# 在庫最適化モジュール optinv と安全在庫配置システム MESSA ## Streamlitで作成した簡易システムの紹介ビデオ 安全在庫配置システム: MESSA (MEta Safety Stock Allocation system) ![](03inventory_files/figure-commonmark/cell-3-output-1.jpeg) ## はじめに 需要の不確実性に対処するための在庫が安全在庫である． １拠点の安全在庫は古典であり，簡単な「公式」もしくはシミュレーションによって最適化できる． ここでは，サプライ・チェイン全体を考えて，各在庫地点の安全在庫を決定するための方策のパラメータの最適化を考える． 多くの研究者の長年の努力にも関わらず，多拠点での一般ネットワーク型の在庫最適化モデルはいまだ未解決の部分も多い． ここでは，適切に付加した仮定の下で，幾つかの最適化モデルを構築する． ## 基礎理論 ### サービスレベルと安全在庫係数 ここでは，小売店に顧客がやってくることによって需要が発生する場合を考える． 顧客は，毎日小売店にやってきて商品を購入するが，その量は気分によって変わるものとする． 一般には，商品が品切れしてしまうことは避けなければならない． 商品を余分に在庫しておく費用が，品切れしてしまう費用より小さい場合には， 商品を余分に在庫して気まぐれな顧客のために準備しておくことが常套手段である． このように，需要の不確実性に対処するためにもつ在庫を**安全在庫**（safety inventory, safety stock）とよぶ． いま，簡単のため，毎日の需要量は，定常な分布にしたがっているものとする． ここで「定常」とは，日によって需要の分布が変化しないことを指す． 分布によっては，負の需要量になる可能性があるが，その場合には需要量を 0 と定義する． たとえば， 需要が正規分布と仮定した場合には，負の部分を除いた正規分布（切断正規分布）を指すものとする． ## Excelインターフェイス ### 在庫最適化用のExcelテンプレートを生成する関数 make_excel_messa 返値： - wb: 在庫最適化用のExcel Workbook ------------------------------------------------------------------------ source ### make_excel_messa ``` python def make_excel_messa( ): ``` *Call self as a function.* #### make_excel_messa関数の使用例 ``` python wb = make_excel_messa() wb.save("messa-templete.xlsx") ``` ### Excel Workbookからデータフレームを生成する関数 prepare_df_for_messa 引数： - wb: Excel Workbook 返値： - stsge_df : 点の情報を保管したデータフレーム - bom_df : 枝の情報を保管したデータフレーム ------------------------------------------------------------------------ source ### prepare_df_for_messa ``` python def prepare_df_for_messa( wb ): ``` *Call self as a function.* ### Excel Workbookから最適化用のデータを準備する関数 prepare_opt_for_messa 引数： - wb: Excel Workbook 返値： - G: 閉路を含まない有向グラフ - ProcTime: 処理時間（NumPyの配列）；以下の配列はグラフGの点を列挙したときの順序（グラフに点を追加した順序）の順に保管してあるものとする． - LTUB: 保証リード時間上限（NumPyの配列） - z: 安全在庫係数 - mu: 需要の平均（NumPyの配列） - sigma: 需要の標準偏差（NumPyの配列） - h: 在庫費用（NumPyの配列） ------------------------------------------------------------------------ source ### prepare_opt_for_messa ``` python def prepare_opt_for_messa( wb ): ``` *Call self as a function.* #### prepare_opt_for_messaの使用例 ``` python 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 ------------------------------------------------------------------------ source ### messa_for_excel ``` python def messa_for_excel( wb, best_NRT, best_MaxLI, best_MinLT ): ``` *Call self as a function.* #### messa_for_excel関数の使用例 ``` python 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") ``` ``` python 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 ``` python 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の図オブジェクト ------------------------------------------------------------------------ source ### generate_pyvis_for_messa ``` python def generate_pyvis_for_messa( wb, G, best_sol:NoneType=None, best_MaxLI:NoneType=None, best_MinLT:NoneType=None, notebook:bool=True ): ``` *Call self as a function.* #### generate_pyvis_for_messa関数の使用例 ``` python net = generate_pyvis_for_messa(wb, G, best_sol, best_MaxLI, best_MinLT) #net.show("example.html") ``` ## 需要予測と合わせた動的在庫最適化 ### 発注固定費用を考えない場合 \#hide 需要予測 *ŷ**t* と需要量 *y**t* の差（誤差）を *e**t* とし， 定常性を仮定する． この分布を，平均 0，標準偏差 *σ* の正規分布とすると，基在庫レベル $S = z \sigma \sqrt{NRT}$ が定まる． ただし *z* は安全在庫係数である． ここで *N**R**T* は正味補充時間であり，補充リード時間 *R**L**T* から保証リード時間 *G**L**T* を減じたものと定義される． 安全在庫のための発注量 *Q**t* と予測に基づく動的な需要のための発注量を分けて考える． 安全在庫のための在庫方策は，予測誤差 *e**t* を需要とした基在庫方策になる． - *I**P**t*: 期首在庫ポジション； 期末在庫ポジション（常に *S*）から需要量 *e**t* を減じた値 *I**P**t* = *S* − *e**t* - *Q**t*: 安全在庫分の発注量； *S* − *I**P**t* と設定する． 負の値もとる． *Q**t* = *S* − *I**P**t* = *S* − (*S* − *e**t*) = *e**t* つまり，誤差の分だけ発注すれば良い． - *O**t*: 動的な需要を考慮した発注量； リード時間後の需要の予測値 *d̂**t* + *R**L**T* − *G**L**T* 分だけ発注するものとする． 負の場合も考えられるが，その場合には 0 とし， 次の期の発注量から減じる． *O**t* = \[*d̂**t* + *R**L**T* − *G**L**T* + *Q**t*\]+ *d̂**t* + *R**L**T* − *G**L**T* + *Q**t* \< 0 のとき： *e**t* + 1 := *e**t* + 1 + *d̂**t* + *R**L**T* − *G**L**T* + *Q**t* 在庫のフローは以下の式で得られる． *I**t* − 1 + *O**t* − *R**L**T* − *d**t* − *G**L**T* = *I**t* 初期在庫量は， 過去の発注量がすべて 0 と仮定したとき， 安全在庫量 *S* に$\sum\_{t=-GLT}^{RLT} \hat{d}\_{t}$ を加えたものとする． $$ I_0 = \sum\_{t=-GLT}^{RLT} \hat{d}\_{t} + S $$ ### 発注固定費用を考慮した場合 \#hide 発注費用を考慮した場合では，動的な需要に対して動的ロットサイズ決定問題を解いて，発注量を決める． サイクル時間 *C**T* は，需要を定常と仮定したときの発注間隔であり，経済発注量の公式で簡単に求めることができる． 安全在庫分は，サイクル時間 *C**T* と正味補充時間 *N**R**T* の大きい方を *L* としたとき，安全在庫係数 *z* と検証期間における予測誤差の標準偏差 *σ* を用いて，以下のように決める． $$ z \sigma \sqrt{T} $$ - *O**t*: 動的な需要を考慮した発注量； まず，需要の予測値 *d̂* に対して動的ロットサイズ決定問題を解くことによって求める． それを正味補充時間 *N**R**T* 分だけ前にずらすことによって *O**t* を決める． - *O*′*t*: 誤差補正した発注量； 需要予測 *ŷ**t* と需要量 *y**t* の差（誤差）を *e**t* を保存し，その累積値 *E**t* を計算する． 動的発注量 *O**t* が正なら， \[*O**t* + *E**t* − 1\]+ を *O*′*t* とし， 累積量を 0 に戻す． ただし，誤差補正発注量が負になる場合には，その分だけ累積誤差に加えておく． 式で表すと以下のようになる． $$ O'\_t = \left\\ \begin{array}{ll} \[O_t+E_t\]^+ & O_t \> 0 \\ 0 & O_t = 0 \end{array} \right. $$ $$ E_t = \left\\ \begin{array}{ll} e_t & O_t \> 0, O_t+E\_{t-1} \geq 0 \\ O_t+E\_{t-1}+e_t & O_t \> 0, O_t+E\_{t-1} \< 0 \\ E\_{t-1} + e_t & O_t = 0 \end{array} \right. $$ 在庫量は，補充リード時間を *R**L**T*，保証リード時間を *G**L**T* としたとき， 以下の式で得られる． *I**t* − 1 + *O*′*t* − *R**L**T* − *d**t* − *G**L**T* = *I**t* ### Wagner-Whitin法による動的ロットサイズ決定問題の求解関数 ww 動的ロットサイズ決定モデルの最適化（モデルについては動的ロットサイズ決定モデルlotsize参照） 需要予測に基づく不確実性を考慮した在庫最適化（forecast）で利用 引数: - demand: 多期間需要を表すリスト - fc: 固定費用（定数かリスト）; 関数内で需要と同じ長さのNumPy配列に変換される． - vc: 変動費用（定数かリスト）; 関数内で需要と同じ長さのNumPy配列に変換される． - h: 在庫費用（定数かリスト）; 関数内で需要と同じ長さのNumPy配列に変換される． 返値： - cost: 最適値 - order: 各期の発注量を入れたリスト ------------------------------------------------------------------------ source ### ww ``` python def ww( demand, fc:float=100.0, vc:float=0.0, h:float=5.0 ): ``` *Call self as a function.* #### ww関数の使用例 ``` python 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: 最良誤差のパラメータ ------------------------------------------------------------------------ source ### best_distribution ``` python def best_distribution( data ): ``` *Call self as a function.* ### best_distribution関数の適用例 ``` python 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) ![](03inventory_files/figure-commonmark/cell-22-output-1.png) ### 分布に適合した表形式の確率分布を求める関数 best_histogram 定常分布をもつ需要データから表形式（度数分布）による確率分布を求める． データ数が少ない（50未満）の場合には，平均値を合わせるため左に0.5だけ*x*軸の区分をシフトさせる． 引数： - data: データの配列 - nbins: ビン数（既定値は50とデータ数の小さい方） 返値： - fig: データの度数分布表（ヒストグラム）と分布の密度関数の図 - hist_dist: 度数分布（パラメータ固定） ------------------------------------------------------------------------ source ### best_histogram ``` python def best_histogram( data, nbins:int=50 ): ``` *Call self as a function.* ### best_histogram関数の適用例 ``` python data = np.random.normal(100,10,1000) fig, hist_dist = best_histogram(data, nbins = 20) plotly.offline.plot(fig); ``` ![](03inventory_files/figure-commonmark/cell-25-output-1.png) 得られた分布に基いてランダムな需要を発生させる． ``` python best_dist.mean(),hist_dist.mean() ``` (100.48413228614606, 100.20054738268533) 得られた分布に対して，在庫管理に必要な統計量を計算する． ``` python 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: 集約を行う期。例えば１週間を１期とするなら”1w”とする。既定値は1日（“1d”） - nbins: ヒストグラムのビン数（既定値は50） - L: リード時間（在庫管理で用いる．L日前からL-1前の需要の合計の分布をとる．） 返値： - dic: 顧客名と製品名のタプルをキーとし，分布の情報を値とした辞書； 分布の情報は以下のタプル - fig: 誤差が最小の分布の密度関数と度数分布表の図 - fig2: 近似的な表形式の確率関数と度数分布表の図 - best_dist: 誤差最小の分布の密度関数 - hist_dist: 度風分布表のデータから構成された確率関数 - best_fit_name: 誤差最小の分布名 - best_fit_params: 誤差を最小にする分布のパラメータ ------------------------------------------------------------------------ source ### fit_demand ``` python def fit_demand( demand_df, cust_selected, prod_selected, agg_period:str='1d', nbins:int=50, L:int=1 ): ``` *需要分布の推定（顧客と製品に対して）:リード時間Lを入力とし，その間の需要の分布を計算する必要がある！* ### fit_demand関数の適用例 ``` python 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:4138: RuntimeWarning: divide by zero encountered in log ``` python 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) ![](03inventory_files/figure-commonmark/cell-31-output-1.png) ## 経済発注量モデル ### **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* が分岐点 *θ* を超えると， その購入単価が *c*0 から *c*1 に（*θ* を超えた分だけ）割引されるものとする． 発注の固定費用を *K* とし，購入費用の合計を，発注量 *Q* の関数として *c*(*Q*) とすると，*c*(*Q*) は $$ c(Q)= \left\\ \begin{array}{l l} 0 & Q=0 \\ K + c_0 Q & 0 \ 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)*θ* である． 各直線ごとの問題は，通常の経済発注量モデルと同じ形になるので，最適解を容易に求めることができる． たとえば，*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* = *θ* で繋がっていることに注意すると， *f*(*Q*0\*) と *f*(*Q*1\*) をくらべて小さい方を選択すれば良いことが言える． 一般に，区分 *j* に対する購入単価を *c**j*(*j* = 0, 1, …, *m*) とすると， $$ Q_j^\* = \sqrt{ \frac{2 d K_j}{ h+r c_j } } $$ *K**j* = *K* + (*c*0 − *c**j*)*θ**j* を区分ごとに計算し，最小費用の区分を求めれば良い． #### 総量割引 総量割引の場合には，発注量 *Q* が分岐点 *θ* を超えると， その購入単価が *c*0 から *c*1 にすべて割引されるので，購入費用 *c*(*Q*) は， $$ c(Q)= \left\\ \begin{array}{l l} 0 & Q=0 \\ K + c_0 Q & 0 \0 \> *c*1 なので， *c*(*Q*) = *K* + *c*1*Q* のときの最適発注量 $$ Q_1^\* =\sqrt{ \frac{2 d K}{ h+r c_1 } } $$ が *θ* 以上のときには，*Q*1\* が最適になる． *Q*1\* \< *θ* のときには，*c*(*Q*) = *K* + *c*0*Q* のときの最適発注量 $$ Q_0^\* =\sqrt{ \frac{2 d K}{ h+r c_0 } } $$ と発注量が *θ* のときの総費用をくらべて，小さい方を選択すれば良い． 一般に，区分 *j* に対する購入単価を *c**j*(*j* = 0, 1, …, *m*) とすると， $$ Q_j^\* =\sqrt{ \frac{2 d K}{ h+r c_j } } $$ を区分ごとに計算し，それが区分内にある場合には，*Q**j*\* を候補とし， 区分外の場合には，分岐点 *θ**j* を候補とすれば良い． ## 経済発注量モデルを解くための関数 eoq 引数： - K: 発注固定費用 - d: 需要量 - h: 在庫費用（割引モデルの場合には，金利分以外の在庫費用） - b: 品切れ費用； Noneの場合には，品切れを許さないモデルになる． - r: 利子率； 割引モデルの場合に使用する． - c: 割引モデルの場合の区間ごとの割引費用； - theta: 割引モデルの場合の区間の切れ目（分岐点）； 最初の区分は必ず0と仮定 - discount:割引の種類を表す文字列を入れる． “incremental”の場合は増量割引， “all”の場合は総量割引， None（既定値）の場合には割引がないモデルになる． 返値： - Q: 最適発注量 - z: 最適値 ------------------------------------------------------------------------ source ### eoq ``` python def eoq( K, d, h, b, r, c, theta, discount:NoneType=None ): ``` *Call self as a function.* ### eoq関数の使用例 ``` python 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* になるように発注する方策である． 保証リード時間 *G**L**T* と補充リード時間 *R**L**T* に対する (*s*, *S*) 方策のシミュレーションは，以下の式に基づく． - 在庫フロー： *I**t* = *I**t* − 1 + *d**t* − *G**L**T* + *Q**t* − *R**L**T* - 正味在庫ポジション： *N**I**t* = *N**I**t* − 1 − *d**t* - 発注量 *Q*: *N**I**t* \< *s* のとき *Q**t* = \[*S* − *N**I**t*\]+ - 正味在庫ポジション （発注後）： *N**I**t* = *N**I**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)；在庫の可視化に用いる。 ------------------------------------------------------------------------ source ### simulate_inventory ``` python def simulate_inventory( n_samples:int=1, n_periods:int=10, mu:float=100.0, sigma:float=10.0, LT:int=3, Q:float=300.0, R:float=100.0, b:float=100.0, h:float=1.0, fc:float=10000.0, S:NoneType=None ): ``` *(Q,R)方策と(s,S)方策による在庫シミュレーション* ### simulate_inventory関数の使用例 ``` python 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() ``` 臨界率= 0.9803921568627451 安全在庫係数= 2.0619165008094615 100.0 100.0 0.0 ``` python plt.plot(I[0]); #在庫の推移の可視化 ``` ![](03inventory_files/figure-commonmark/cell-37-output-1.png) ### 多段階の場合のシミュレーション ------------------------------------------------------------------------ source ### multi_stage_simulate_inventory ``` python def multi_stage_simulate_inventory( n_samples:int=1, n_periods:int=10, mu:float=100.0, sigma:float=10.0, LT:list=[1, 1, 1], s:NoneType=None, S:NoneType=None, b:float=100.0, h:list=[10.0, 5.0, 2.0], fc:float=10000.0 ): ``` *多段階ネットワークに対する(s,S)方策による在庫シミュレーション* ``` python 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]); ``` Echelon LT= [1 3 5] [218.88226596 433.36782555 650.50643525] [218.88226596 433.36782555 650.50643525] 2800.004477798352 ![](03inventory_files/figure-commonmark/cell-39-output-2.png) ## (Q,R)方策のパラメータ最適化関数 optimize_qr *Q* が与えられていないときには， *Q*\* は経済発注量モデルを用いて初期化する． *Q*\* を固定し， 発注点 *R* は，リード時間内の安全在庫量を初期値として1つずつ増やし，費用が増加した地点を最適な発注点 *R*\* とする． *R*\* を固定し， 経済発注量の前後を探索し，費用が増加に転じた点を最適な発注量 *Q*\* とする． 期待値だけでなく費用の標準偏差も考慮し， 「期待値 +*α* 標準偏差」 を評価尺度とする． シミュレーションベースなので，任意の分布，任意の確率過程， 任意の制約に対応できる． 引数： - n_samples: シミュレーションを行うサンプル数 - n_periods: 期間数 - mu: 需要の平均（需要は正規分布を仮定する．） - sigma: 需要の標準偏差 - LT: リード時間 - Q: 発注ロットサイズ； Qを固定したいときには整数を入れ，最適化したいときには None を入れる． - R: 発注点； Rを固定したいときには整数を入れ，最適化したいときには None を入れる． - z: 安全在庫係数； None（既定値）の場合には，在庫費用とバックオーダー費用から計算される． - b: バックオーダー費用 - h: 在庫費用 - fc: 発注固定費用 - alpha: リスクの考慮の度合い；「期待値 +*α* 標準偏差」 を評価尺度とする． 返値： - R: 最適化された発注点 - Q: 最適化された発注量 ------------------------------------------------------------------------ source ### optimize_qr ``` python def optimize_qr( n_samples:int=1, n_periods:int=10, mu:float=100.0, sigma:float=10.0, LT:int=3, Q:NoneType=None, R:NoneType=None, z:NoneType=None, b:float=100.0, h:float=1.0, fc:float=10000.0, alpha:float=1.0 ): ``` *(Q,R)方策のパラメータ最適化* ### optimize_qr関数の使用例 ``` python 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]); #在庫の推移の可視化 ``` ![](03inventory_files/figure-commonmark/cell-41-output-1.png) ``` python #費用のヒストグラム plt.hist(cost); ``` ![](03inventory_files/figure-commonmark/cell-42-output-1.png) ## (s,S)方策のパラメータ最適化関数 optimize_ss *Q* が与えられていないときには， *Q*\* は経済発注量モデルを用いて初期化する． *Q*\* を固定し， 発注点 *R* = *s* は，リード時間内の安全在庫量を初期値として1つずつ増やし，費用が増加した地点を最適な発注点 *R*\* とする． *R*\* を固定し， 「*R*\*+経済発注量」の前後を探索し，費用が増加に転じた点を最適な *S*\* とする． 期待値だけでなく費用の標準偏差も考慮し， 「期待値 +*α* 標準偏差」 を評価尺度とする． シミュレーションベースなので，任意の分布，任意の確率過程， 任意の制約に対応できる． 引数： - n_samples: シミュレーションを行うサンプル数 - n_periods: 期間数 - mu: 需要の平均（需要は正規分布を仮定する．） - sigma: 需要の標準偏差 - LT: リード時間 - Q: 発注ロットサイズ； Qを固定したいときには整数を入れ，最適化したいときには None を入れる． - R: 発注点； Rを固定したいときには整数を入れ，最適化したいときには None を入れる． - z: 安全在庫係数； None（既定値）の場合には，在庫費用とバックオーダー費用から計算される． - b: バックオーダー費用 - h: 在庫費用 - fc: 発注固定費用 - alpha: リスクの考慮の度合い；「期待値 +*α* 標準偏差」 を評価尺度とする． - S: 基在庫レベル； 最適化したいときには None（既定値）を入れる． 返値： - s: 最適化された発注点 - S: 最適化された基在庫レベル ------------------------------------------------------------------------ source ### optimize_ss ``` python def optimize_ss( n_samples:int=1, n_periods:int=10, mu:float=100.0, sigma:float=10.0, LT:int=3, Q:NoneType=None, R:NoneType=None, z:NoneType=None, b:float=100.0, h:float=1.0, fc:float=10000.0, alpha:float=1.0, S:NoneType=None ): ``` *(s,S)方策のパラメータ最適化* ### optimize_ss関数の使用例 ``` python 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.8822659611578 218.8822659611578 763.7424638470262 764.4567798777746 764.5266594870563 765.2828889754412 768.528185341823 770.2273226521506 778.5564569609345 784.1213191908391 788.7634416969585 796.8654287662399 ![](03inventory_files/figure-commonmark/cell-44-output-2.png) ``` python b,h ``` (100.0, 10.0) ``` python #費用のヒストグラム plt.hist(cost); ``` ![](03inventory_files/figure-commonmark/cell-46-output-1.png) ### (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: 最適化された基在庫レベル 深層学習を用いた拡張が課題となる． ------------------------------------------------------------------------ source ### approximate_ss ``` python def approximate_ss( mu:float=100.0, sigma:float=10.0, LT:int=0, b:float=100.0, h:float=1.0, fc:float=10000.0 ): ``` *Call self as a function.* ### approximate_ss関数の使用例 ``` python 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 ``` python #費用のシミュレーション 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]); #在庫の推移の可視化 ``` R= 323.12595676092786 R= 323 Q= 461 RQ-policy cost= 4714.398125108141 s= 313 S= 772 s,S policy cost= 4686.765178752777 s= 349.55651923599106 S= 768.2721946765712 Approximate s,S policy cost= 4855.83862088315 ![](03inventory_files/figure-commonmark/cell-49-output-2.png) ### 多段階への拡張 *Q*\* が決まっている場合： サイクル時間 *C**T* = *Q*\*/*μ* を計算し， *L* = max {*C**T*, *L**T*} をリード時間 *L* として最適な基在庫レベル *S*\* を定期発注方策最適化で求める． *L*分の平均在庫量 *L**μ* を減じて安全在庫量を求め， さらにリード時間分の平均需要量 *μ**L**T* を加えたものが発注点 *R*\* となる． *Q*\* を最適化した場合： 定期発注方策最適化で求めた *R*\* を固定し， 最適な発注量 *Q*\* を地点ごとに（上の関数optimize_qrを用いて）再最適化し，これを繰り返す． 多段階の (*Q*, *R*)方策と (*s*, *S*) 方策のシミュレータを作り， パラメータの調整を行う． ### (s,S)方策の最適パラメータの可視化 ``` python 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 ``` python 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); ``` ![](03inventory_files/figure-commonmark/cell-51-output-1.png) ## 定期発注方策（単一段階モデル） \#hide ここでは，発注を行うことができるタイミングが限定されている定期発注方策を考える． 定期発注方策では，時間を離散値として取り扱う．実際には，日単位などで在庫を管理することが多いため， 定期発注方策は，多くの実際問題で適用可能なモデルである． 時間は，期とよばれる単位に分割されており，期を 1, 2, …, *t*, …, *t*max と記す． はじめに，1つの小売店に顧客が需要を引き起こす場合を考える． モデルに必要な記号を導入する． - *D**t*: *t* 期における需要量．任意の分布もしくは確率過程． ここでは，シミュレーションを行いながら，基在庫レベルを最適化するので，シミュレーション可能な任意の 需要を仮定して良い． - *L*: リード時間（lead time）；*t* 期の期末に発注された商品は，*t* + *L* + 1 期の期首に 到着するものと仮定する．この発注から到着の間の時間差をリード時間とよぶ． - *s*: 基在庫レベル． 各期の期末に，期末在庫ポジションが基在庫レベル *s* になるように，上流の在庫地点に対して，発注するものとする． 期首から期末にかけて需要 *D**t* が発生し，その後で発注するものと仮定する． - *c*: 発注量の上限．容量ともよばれる． - *q**t*: *t* 期における発注量． *t* 期の期末に，期末の在庫ポジション（在庫量に輸送中在庫量を加えたもの）*I**t* − *D**t* + *T**t* が基在庫レベル *s* になるように発注するものとする．ただし，発注量の上限を超えた場合には， 上限いっぱいの量を発注するものとする． これを式で表すと，以下のようになる． *q**t* = min {*c*, \[*s* − (*I**t* − *D**t* + *T**t*)\]+} ここで，\[⋅\]+ は，max {0, ⋅} を表す記号である． - *I**t*: *t* 期の期首における正味在庫量；期首に以前注文した量を受け取った後の在庫量で測定するものとする． なお，正味在庫量は，「実在庫量 − バックオーダー（品切れ）量」を表す， したがって，*t* + 1 期の期首在庫量 *I**t* + 1 は，*t* 期の期首在庫量 *I**t* から需要量 *D**t* を減じた後， *t* + 1期の期首に *t* − *L* 期に発注された量 *q**t* − *L* を受け取るので， 在庫量と需要量，発注量の関係式（在庫フロー保存式）は，以下のようになる． *I**t* + 1 = *I**t* − *D**t* + *q**t* − *L* - *T**t*: *t* 期の期首における輸送中在庫量； 在庫量と同様に，以前注文した量を受け取った後で測定するものとする． *t* + 1 期の期首の輸送中在庫量は，前の期の輸送中在庫量 *T**t* に，新たに発注した量 *q**t* を加え， *t* + 1 期の期首に到着する量 *q**t* − *L* を減じるので，以下のように計算できる． *T**t* + 1 = *T**t* + *q**t* − *q**t* − *L* - *b*: 単位期間，単位バックオーダー量あたりの品切れ費用． - *h*: 単位期間，単位在庫量あたりの在庫保管費用． *t* 期の費用を表す確率変数 *C**t* は，期首における（過去の注文を受け取った後の）正味在庫量で測定するものと仮定すると， *C**t* = *b*\[*I**t*\]− + *h*\[*I**t*\]+ と書くことができる． 目的関数は，*t*max 期の期待費用 $$ \frac{1}{ t\_{\max} } \sum\_{t=1}^{t\_{\max}} \mathbf{E}\[ C_t \] $$ であり，これを最小化する． ここで，モデルの仮定に関する注意を述べておく． 期中に需要が発生し，期首に（注文を受け取った後で）在庫を評価すると仮定していたので， 期首の在庫が 0 でも品切れが生じないことになる（期中に予約をして，翌日の朝に受け取ると解釈もできるが， 多少無理がある）． また，発注を行うタイミングも，実際には期末でなく，期中に行われることが多い． たとえば，あるスーパーでは， 午前中の暇な時間帯に（POSシステムがあるにも関わらず）実在庫を調べ，発注を行っている． ここで行ったモデルの仮定は，記号の簡潔さと従来の研究との整合性のためであり， 実際にはシミュレーションを行いながら，パラメータの最適化を行うので，ほとんどの現実問題に対処可能である． 変数は基在庫レベル *s* である． 目的関数の *s* による微分値 (*I**t*)′ は、単一段階、モデルにおいてはすべて 1 になるので、 *t* 期の費用関数 *C**t* の *s* による微分は， (*C**t*)′ = −*b*(*I**t*)′**1**\[*I**t* \< 0\] + *h*(*I**t*)′**1**\[*I**t* \> 0\] となる． また，目的関数（*C**t* の期待値）の *s* による微分は，確率 1 で微分値の期待値 $$ \mathbf{E}\left\[ \frac{1}{t\_{\max}} \sum\_{t=1}^{t\_{\max} } (C_t)' \right\] $$ に収束することが示される． したがって，(*C**t*)′ の期待値をとり，これが負ならば基在庫レベル *s* を増やし，正ならば減らせば良いことが分かる． 実際には，非線形計画による探索を行う必要があるが，通常の（微分値のみを利用した比較的単純な）非線形最適化の手法が適用できる． 上のアルゴリズムを実装したものが、以下のbase_stock_simulationである。 引数： - n_samples: シミュレーションを行うサンプル数 - n_periods: 期間数 - demand: 需要のサンプルを表すNumPyの配列で、形状は (n_samples, n_periods)。 - capacity: 容量（生産量上限） - LT: リード時間 - b: バックオーダー費用 - h: 在庫費用 - S: 基在庫レベル 返値： - dC： 各地点における総費用の基準在庫レベルによる微分の期待値 - total_cost: 総費用の期待値 - I: 在庫の時間的な変化を表すNumPy配列で、形状は (n_samples, n_periods+1)。 \#hide ““” Solve the base-stock level optimization problem using the simulation based approach (infinitestimal perterbation analysis: IPA) The algorithm is based on: @article{Glasserman1995, author={P. Glasserman and S. Tayur}, title={Sensitivity Analysis for Base Stock Levels in Multi-Echelon Production-Inventory Systems}, journal={Management Science}, year={1995}, volume={41}, pages={263--281}, annote={ } } ““” ------------------------------------------------------------------------ source ### base_stock_simulation ``` python def base_stock_simulation( n_samples, n_periods, demand, capacity, LT, b, h, S ): ``` *単一段階在庫システムの定期発注方策に対するシミュレーションと微分値の計算* ------------------------------------------------------------------------ source ### base_stock_simulation_using_dist ``` python def base_stock_simulation_using_dist( n_samples, n_periods, demand, capacity, LT, b, h, S ): ``` *単一段階在庫システムの定期発注方策に対するシミュレーションと微分値の計算* ``` python # n_samples = 10 # 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. ``` ``` python # print("t: S dS Cost") # 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を下回ったら終了する。 ``` python #(Q,R)方策と同じ仮定で実験 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 ``` python 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 ``` python plt.plot(I[0,:1000]); ``` ![](03inventory_files/figure-commonmark/cell-65-output-1.png) ## 定期発注方策（多段階在庫モデル） \#hide ここでは，単一段階モデルを直列多段階の場合に拡張する． この結果は，Glasserman―Tayur に基づく． 多くの実際問題は，直列多段階に分解して考えることができるので，このモデルは実務的にも有意義であるが， 以下で解説する一般のネットワークに対する基礎を与えるという意味でも重要である． 点の順番は下流（需要地点側）から 0, 1, …, *n* − 1 と並んでいるものとし， *n* 番目には，点 *n* − 1 に品目を供給する外部の点があるものと仮定する． この *n* 番目の点には十分な在庫があり，品切れはしないものとする． 多段階の問題では，基在庫レベルのかわりに，エシェロン基在庫レベルを用いる． 各段階におけるエシェロン在庫ポジションが，エシェロン基在庫レベルより小さいときには， エシェロン基在庫レベルになるように，上流の地点に発注を行う． この方策を，エシェロン基在庫方策とよぶ． 直列に並んだ在庫地点を段階とよび，下流から数えて *i* 番目の在庫地点を，第 *i* 段階とよぶ． 単一段階のときの記号に， 上添字 *i* をつけたものを，多段階用の記号とする． たとえば，第 *i* 段階の *t* 期における正味在庫量は， *I**t**i* と記す． この記号を用いると，エシェロン基在庫方策は，エシェロン基在庫レベル *s**i* を与えたとき， 期末のエシェロン在庫ポジション $$ \sum\_{j=0}^i (I_t^j +T_t^j) -D_t $$ を *s**i* にするように発注を行う方策と定義できる． 実際には，発注量の上限（容量）*c**i* と上流の在庫地点の在庫量 *I**t**i* + 1 によっても制限を受けるので， 発注量 *q**t**i* は， $$ q^i_t = \min \left\\ c_i, \left\[ s^i + D_t- \sum\_{j=1}^i (I^j\_{t} + T^j\_{t}) \right\]^+, \[ I_t^{i+1}\]^+ \right\\ $$ となる．ただし，*i* = *n* − 1 のときには，*I**t**i* + 1 を除く． 正味在庫と輸送中在庫の再帰式は，それぞれ *I**t* + 1*i* = *I**t**i* − *q**t**i* − 1 + *q**t* − *L**i**i* *T**t* + 1*i* = *T**t**i* + *q**t**i* − *q**t* − *L**i**i* となる． ただし，*i* = 0 のときには，上式における， 上流の地点からの発注量 *q**t**i* − 1 は需要量 *D**t* に置き換える． 初期条件は，在庫量については *I*00 = *s*0， *I*0*i* = *s**i* − *s**i* − 1, *i* = 1, 2, …, *n* − 1 と設定し， その他の変数については，すべて 0 と設定する． *t*期の費用を表す確率変数 *C**t* は， $$ C_t = b \[I^0_t\]^- + h^0 \[I^0_t\]^+ + h^0 T_t^0 +\sum\_{j=1}^{n-1} h^j (I^j_t +T^{j}\_t) $$ と書くことができる． ここで， *b* は品切れ費用であり，最終需要地点（点0）でのみ課せられるものとする． *h**i* は *i*段階の在庫費用である（エシェロン在庫費用でないことに注意されたい）． また，ここでは，輸送中在庫量に対する在庫費用は，下流の在庫地点の在庫費用で計算するものと仮定している． 目的関数は，*t*max 期の期待費用 $$ \frac{1}{t\_{\max}} \sum\_{t=1}^{t\_{\max}} \mathbf{E}\[ C_t \] $$ であり，これを最小化する． 変数は基在庫レベル *s**i*, *i* = 0, 1, …, *n* − 1 である．シミュレーションをしながら， 目的関数の各 *i*( = *i*\*) に対する *s**i*( = *s*\*) による微分値の期待値を計算する． 生産量計算式の微分は，以下のような場合分けを必要とする． $$ \frac{d q^i_t}{d s^\*} = \left\\ \begin{array}{l l } 0 & 容量制約が制約のとき \\ 0 & 発注量が 0 のとき \\ (I_t^{i+1})' & i+1 番目の地点の在庫が制約のとき \\ \mathbf{1}\[ i=i^\*\] - \sum\_{j=0}^i \left( \frac{d I^j\_{t}}{d s^\*} +\frac{d T^j\_{t}}{d s^\*} \right) & それ以外 \end{array} \right. $$ 在庫フロー保存式を *s*\* で微分すると， $$ \frac{d I^i\_{t+1}}{d s^\*} = \frac{d I^i\_{t}}{d s^\*} - \frac{d q^{i-1}\_t}{d s^\*} + \frac{d q^i\_{t-L_i}}{d s^\*} $$ となる．ただし，*i* = 0 のときには，右辺の第2項目 *d**q**t**i* − 1/*d**s*\* は除くものとする． 同様に，輸送中在庫の計算式を *s*\* で微分すると， $$ \frac{d T^i\_{t+1}}{d s^\*} = \frac{d T^i\_{t}}{d s^\*} + \frac{d q^i\_{t}}{d s^\*} - \frac{d q^i\_{t-L_i}}{d s^\*} $$ となる． 上の再帰式を使うと，微分値をシミュレーションをしながら順次更新していくことができる． 在庫量の微分値の初期値は，在庫量の初期値の設定から、 (*I*00)′ = **1**\[*i*\* = 0\], (*I*0*i*)′ = **1**\[*i*\* = *i*\] − **1**\[*i*\* = *i* − 1\], *i* = 1, ⋯, *n* − 1 となる． 他の変数の微分値は，すべて 0 である． *t* 期の費用関数 *C**t* の *s*\* による微分は， $$ (C_t)' = -b (I^0_t)' \mathbf{1}\[ I^0_t\<0\] + h_0 (I^0_t)' \mathbf{1}\[ I^0_t\>0 \] + h_0 (T_t^0)' +\sum\_{i=1}^{n-1} h_i \left\\ (T^i_t)' +(I^{i-1}\_t)' \right \\ $$ となる． 目的関数（*C**t* の期待値）の *s*\* による微分は，確率 1 で微分値の期待値 $$ \mathbf{E}\left\[ \frac{1}{t\_{\max}} \sum\_{t=1}^{t\_{\max}} (C_t)' \right\] $$ に収束するので，(*C**t*)′ の期待値をもとに，*s*\* の更新を行う． 上のアルゴリズムを実装したものが、以下のmulti_stage_base_stock_simulationである。 ### 多段階在庫システムの定期発注方策に対するシミュレーションと微分値の計算関数 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)。 ------------------------------------------------------------------------ source ### multi_stage_base_stock_simulation ``` python def 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は、微分値を返すので、最急降下法を適用して最適値を探索する。 微分値の２乗和がepsilonを下回ったら終了する。 ``` python # 多段階 基在庫レベル方策 在庫シミュレーション 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 ] \#hide \### Adamによる最適化 深層学習で定番のAdamで最適化をしてみたが、収束は悪い。 ``` python #可視化（最初のサンプルを表示） 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' ![](03inventory_files/figure-commonmark/cell-70-output-1.png) ## 定期発注方策（一般ネットワークモデル） \#hide 実際には，組立工程があったり，複数の地点から発注を受けたるする場合がある． 以下では，一般のネットワークの場合に，上の手法を拡張することを考える． ある地点 *i* が複数の地点から発注を受け，さらにその地点自身でも需要が発生する場合を考える． すべての需要（発注）に応えるだけの在庫量がない場合に， 複数の需要（発注）に対して，どのように在庫を割り振るかを決める必要がある． これを**分配ルール**（allocation rule）とよぶ． 代表的な分配ルールには，以下のものがある． - 優先順位法: 直接需要ならびに（複数の）下流の在庫地点に対して，優先順位を付けておき，その順に需要（発注）量を満たす． - 発注量比例配分法: 需要（発注）量に応じて比例配分をする． この方法だと，供給配分に起因する鞭効果が発生し，安全在庫量が増えてしまう危険性がある． - 実績比例配分法: 過去の実績に応じて，在庫の比例配分を行う．この方法だと，ある顧客の過去の実績が多く，かつ今回の発注分が少ない場合に，過剰な配分をしてしまい，（この顧客用に確保された）在庫が残っているにもかかわらず，他の顧客の需要を満たすことができないという理不尽な現象が起きる危険性がある． - 最適配分法: 線形計画問題を解くことによって，過去の実績から計算された優先度と配分量の積の合計を最大化するように配分する． 以下では，実績比例配分法を例として解説する． ネットワークを表現するための記号を導入しておく． - *P**A**R**E**N**T**i*: ネットワークにおける点 *i* の直後の点の集合であり，点 *i* が直接供給を行う地点の集合を表す． - *A**N**C**E**S**T**O**R**i*: 点 *i* からネットワークの有向枝を辿って到達可能な点に対応する点集合． 点 *i* 自身も含めた点集合と定義する． - *C**H**I**L**D**i*: ネットワークにおける点 *i* の直前の点の集合であり，点 *i* が発注を行う点の集合を表す。 *C**H**I**L**D**i* に発注したすべての部品が揃わないと生産できないものと仮定する。 - *ϕ**i**j*: 点 *j* に対応する品目を 1単位製造するために必要な点 *i* の品目の量． これは，ネットワークが部品展開表を表すときに用いられるパラメータであり，親品目 *j* を製造するために必要 な子品目 *i* の量を，枝 (*i*, *j*) に付加したものである． 品目の量を表す単位として，*i*-unit，*j*-unitを導入しておくと分かりやすい． ちなみに，*ϕ**i**j* の単位（次元）は，\[*i*-units / *j*-unit \] である． 上流の地点が存在しない，言い換えれば *C**H**I**L**D**i* が空の点に対しては，その上流の地点に十分な 在庫があるものとし，品切れはしないものとする． *t* 期の期首在庫量が与えられている状態で，需要が発生する．需要はすべての点で発生する可能性があるものとする． （つまり，部品に対しても需要があるものと仮定する．）点 *i* における *t*期の需要を *D**t**i* と記す． 期末に，割当ルールによって，下流の地点に対する在庫の配分を行い， 次いで発注量を決定する．これによって，次期の期首在庫 *I**t* + 1*i* が定まる． 在庫地点 *i* の在庫地点の *t* 期における正味在庫量を *I**t**i* と記す． この記号を用いると，エシェロン基在庫方策は，エシェロン基在庫レベル *s**i* を与えたとき， *t* 期末のエシェロン在庫ポジション *E**i**t* = ∑*j* ∈ *A**N**C**E**S**T**O**R**i**ρ**i**j*(*I**t**j* + *T**t**j* − *D**t**j*) を *s**i* にするように発注を行う方策と定義できる． 複数の下流の在庫地点からの注文に対して， 現在の在庫で対応できない場合には，過去の実績に応じて，比例配分を行うものとする． 在庫地点 *i* の在庫を下流の在庫地点 *j* の配分するときの比率（配分比）を *α**i**j* とする． また，在庫地点 *i* が直接の需要 *D**t**i* を持っている場合には， その需要に対する配分比を *α**i*0 とする． 配分比は， *α**i*0 + ∑*j* ∈ *P**A**R**E**N**T**i**α**i**j* = 1 を満たすように設定する． 発注量の上限（容量）*c**i* と上流の在庫地点の在庫量によっても制限を受けるので， 発注量 *q**t**i* は， $$ q_t^i =\min \left\\ \begin{array}{l } c_i \\ \left\[ s_i - E_t^i \right\]^+ \\ \min\_{k \in CHILD_i} \alpha\_{ki} \[ I_t^k \]^+ /\phi\_{ki} \end{array} \right. $$ と決められる． 正味在庫量の計算は， *I**t* + 1*i* = *I**t**i* − ∑*j* ∈ *P**A**R**E**N**T**i**ϕ**i**j**q**t**j* − *D**t**i* + *q**t* − *L**i**i* と行う． 最後に輸送中在庫量は， *T**t* + 1*i* = *T**t**i* + *q**t**i* − *q**t* − *L**i**i* と更新する． 変数の初期値は，在庫量については， *I*0*i* = *s**i* − ∑*j* ∈ *P**A**R**E**N**T**i**ϕ**i**j**s**j* と設定し，その他の変数は，すべて 0 に設定する． *t* 期の費用を表す確率変数 *C**t* は， $$ C_t = \sum\_{i=1}^n b^i \[ I^i\_{t} \]^- + \sum\_{i=1}^n h^i \left\\ \[I^i_t\]^+ + T^{i}\_t \right\\ $$ と書くことができる． ここで，*b**i* は，点 *i* に対する直接需要の品切れ費用であり， *h**i* は *i* 段階の在庫費用である． また，輸送中在庫量に対する在庫費用は，下流の在庫地点の在庫費用で計算するものと仮定している． 微分値の計算は，直列多段階の場合と同様に，以下のようになる． 実績比例配分法の場合には，発注量の微分値の更新式は， $$ \frac{d q^i_t}{d s^\*} = \left\\ \begin{array}{l l } 0 & 容量制約が制約のとき \\ 0 & 発注量が 0 のとき \\ \alpha\_{ki} (I_t^{k})' & 上流の点 k の在庫が制約のとき \\ \mathbf{1}\[ i=i^\*\] - \sum\_{j \in ANCESTOR_j } \rho\_{ij} \left( \frac{d I^j\_{t}}{d s^\*} +\frac{d T^j\_{t}}{d s^\*} \right) & それ以外 \end{array} \right. $$ となる． 在庫の微分値の更新式は， $$ \frac{d I^i\_{t+1}}{d s^\*} = \frac{d I^i\_{t}}{d s^\*} - \sum\_{j \in PARENT_i} \phi\_{ij} \frac{d q^{j}\_t}{d s^\*} + \frac{d q^i\_{t-L_i}}{d s^\*} $$ となる． 最後に，輸送中在庫量の更新式は，直列の場合と同様に， $$ \frac{d T^i\_{t+1}}{d s^\*} = \frac{d T^i\_{t}}{d s^\*} + \frac{d q^i\_{t}}{d s^\*} - \frac{d q^i\_{t-L_i}}{d s^\*} $$ となる． 各点 *i* の在庫量の，点 *i*\* のエシェロン基在庫レベル *s*\* による微分値の初期値は， $$ \frac{d I_0^i}{d s^\*} = \mathbf{1} \[i=i^\*\] -\sum\_{j \in PARENT_i} \phi\_{ij} \mathbf{1}\[ j=i^\*\] $$ と設定し，その他の変数の微分値は，すべて 0 に設定する． *t* 期の費用関数 *C**t* の *s*\* による微分は， $$ (C_t)' = - \sum\_{i=1}^n b^i (I^i_t)' \mathbf{1}\[I^i_t \<0\] + \sum\_{i=1}^n h^i \left\\ (I^i_t)' \mathbf{1}\[ I^i_t\>0 \] + (T^i_t)' \right \\ $$ となる． 前節までと同様に， 目的関数（*C**t* の期待値）の *s*\* による微分は，確率 1 で微分値の期待値 に収束するので，(*C**t*)′ の期待値をもとに，*s*\* の更新をする非線形最適化を行う． また，収束には時間がかかるので，良い初期解から出発することが肝要である． ここでは，各点上での生産時間を 1 期（日）としたときの点 *i* におけるエシェロンリード時間 *E**L**T**i*， 累積需要の標準偏差（の推定値） *σ**i*，安全在庫係数の推定値 *z**i* を用いて，以下のように決めるものとする． $$ (ELT_i-1) \left( \sum\_{j \in ANCESTOR_i } \rho\_{ij} D^j \right) + \sigma_i z_i \sqrt{ELT_i-1} $$ ここで，エシェロンリード時間から 1 期（日）を減じるのは， 需要地点以外の下流の地点をもつ在庫地点では，在庫量が 0 だと，下流の地点の生産ができないので， 必ず一定量の在庫をもつ必要があるが，最終需要地点における在庫は 0 に近づけるようにするためである． ### ネットワーク型定期発注方策に対するシミュレーションと微分値の計算関数 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)。 ------------------------------------------------------------------------ source ### network_base_stock_simulation ``` python def network_base_stock_simulation( G, n_samples, n_periods, demand, capacity, LT, ELT, b, h, S, phi, alpha ): ``` *ネットワーク型定期発注方策に対するシミュレーションと微分値の計算* ### 初期基在庫レベルとエシェロンリード時間を計算する関数 initial_base_stock_level 定期発注方策を仮定したときの、エシェロンリード時間を計算し、それをもとに基在庫レベルを計算して返す関数。 エシェロンリード時間は、最終需要地点以外の点に対しては、1日を加えて計算する。 ------------------------------------------------------------------------ source ### initial_base_stock_level ``` python def initial_base_stock_level( G, LT, mu, z, sigma ): ``` *初期基在庫レベルとエシェロンリード時間の計算* ### network_base_stock_simulation関数を利用した在庫最適化 上で作成した関数を用いて、最適化を行う。network_base_stock_simulationは、微分値を返すので、最急降下法を適用して最適値を探索する。 微分値の２乗和がepsilonを下回ったら終了する。 ``` python 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 ``` ELT= [6. 1. 2.] S= [1848.49989691 216.5 235.00178567] t: Cost |dC| S 0: 3019.95, 17.08926 [1847.62530563 217.49740973 238.91707905] 1: 3006.19, 5.82522 [1846.8067274 217.51097292 241.18753291] 2: 3000.89, 2.75480 [1846.02908948 217.50571214 242.65383682] 3: 2998.09, 1.62961 [1845.29737183 217.49139168 243.69977714] 4: 2996.28, 1.06417 [1844.60260002 217.47909228 244.46221195] 5: 2995.01, 0.72428 [1843.94828582 217.46319584 245.00617509] 6: 2994.07, 0.54494 [1843.32732419 217.45502006 245.40527113] 7: 2993.30, 0.43767 [1842.730011 217.45163114 245.68966556] 8: 2992.67, 0.35866 [1842.1640293 217.45166499 245.88543662] 9: 2992.12, 0.31603 [1841.61648417 217.44484704 246.01264073] 10: 2991.63, 0.28328 [1841.09019154 217.44626505 246.09199254] 11: 2991.18, 0.25641 [1840.5857877 217.44053849 246.13615565] 12: 2990.78, 0.23289 [1840.1034789 217.43167377 246.15003918] 13: 2990.40, 0.21733 [1839.63730977 217.43089205 246.1457369 ] 14: 2990.04, 0.19542 [1839.19629673 217.42285795 246.11627852] 15: 2989.72, 0.18277 [1838.77053194 217.42063879 246.07771032] 16: 2989.41, 0.16754 [1838.36472953 217.41741412 246.02429507] 17: 2989.13, 0.15640 [1837.97370051 217.41310713 245.96531407] 18: 2988.86, 0.14614 [1837.59639766 217.4098182 245.90386064] 19: 2988.61, 0.13544 [1837.23391051 217.40119201 245.84089142] 20: 2988.37, 0.12520 [1836.8848687 217.40228999 245.78289211] 21: 2988.14, 0.11705 [1836.54740286 217.40023727 245.72661672] 22: 2987.93, 0.11064 [1836.21885629 217.39874373 245.67472198] 23: 2987.71, 0.10247 [1835.90308804 217.39500345 245.62232573] 24: 2987.52, 0.09693 [1835.59618264 217.40030404 245.57027426] 25: 2987.32, 0.09327 [1835.29460321 217.40133181 245.52210444] 26: 2987.13, 0.08379 [1835.00956215 217.39951229 245.47171922] 27: 2986.96, 0.08055 [1834.72877887 217.39764763 245.43037613] 28: 2986.79, 0.07318 [1834.46321624 217.39425339 245.37892772] 29: 2986.64, 0.06854 [1834.20645251 217.39259184 245.32788804] 30: 2986.49, 0.06316 [1833.95949445 217.3895499 245.28136566] 31: 2986.35, 0.05834 [1833.72243734 217.38608544 245.23517302] 32: 2986.21, 0.05304 [1833.49542662 217.38626567 245.19632617] 33: 2986.08, 0.04822 [1833.27944085 217.3870057 245.15668505] 34: 2985.96, 0.04516 [1833.0698292 217.38347793 245.12192358] 35: 2985.85, 0.04205 [1832.86827297 217.38137318 245.08422204] 36: 2985.74, 0.03913 [1832.67413149 217.38133212 245.04626724] 37: 2985.63, 0.03586 [1832.4898897 217.3808698 245.00250501] 38: 2985.54, 0.03389 [1832.31055889 217.37865597 244.96099329] 39: 2985.45, 0.03082 [1832.14016615 217.37908883 244.91872806] 40: 2985.36, 0.02814 [1831.97646732 217.3801999 244.88204575] 41: 2985.28, 0.02513 [1831.82350194 217.37479015 244.8407684 ] 42: 2985.21, 0.02265 [1831.67971848 217.37178324 244.79644754] 43: 2985.14, 0.02048 [1831.54170044 217.37143618 244.7586381 ] 44: 2985.08, 0.02011 [1831.40387548 217.377215 244.72575092] 45: 2985.01, 0.01788 [1831.27363076 217.37646674 244.6955142 ] 46: 2984.95, 0.01682 [1831.14750214 217.37581502 244.66525365] 47: 2984.89, 0.01627 [1831.02248436 217.37510563 244.64002069] 48: 2984.83, 0.01352 [1830.91097271 217.37036853 244.60741999] 49: 2984.78, 0.01248 [1830.80391694 217.37170715 244.57547878] 50: 2984.73, 0.01153 [1830.70035915 217.36954778 244.54720106] 51: 2984.69, 0.01034 [1830.60259375 217.37002131 244.51926937] 52: 2984.64, 0.00965 [1830.50919519 217.36718115 244.48901697] ``` python #可視化（最初のサンプルを表示） 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); ``` ![](03inventory_files/figure-commonmark/cell-75-output-1.png) \#hide \## 顧客の需要予測と安全在庫 上では、定常なランダム需要を仮定して、安全在庫量の最適化を行なったが、実際問題においては需要は定常とは限らない。 過去の需要系列demand_dfをもとに、Prophetによる予測を行う関数 forecastを準備する。 比較として、（需要の定常性を仮定して）平均と標準偏差も算出する。 引数： - demand_df: 需要の系列を表すデータフレーム。顧客Cust、製品Prod,日時Date、需要量demandの列をもつものと仮定する。 - agg_period: 集約を行う期。例えば１週間を１期とするなら”1w”とする。既定値は1日（“1d”） - forecast_period: 予測を行う未来の期数。既定値は1 返値： 以下の情報を列情報にもつデータフレーム - mean: 需要の平均 - std : 需要の標準偏差 - 0..forecast_period-1: Prophetによる予測から計算された未来のforecast_periods 期分の最大在庫量 \#hide \### forecast関数の使用例 １週間を1期として、30期を予測。 非常に時間がかかるので、注意。 \#hide \## Prophetによる予測と需要の平均・標準偏差を用いた安全在庫量の比較と可視化 Prophetによる可視化 \#hide \## 非定常な需要予測をもとにした基在庫レベルの最適化 以下では、実際問題を想定して、非定常な需要をもとに、安全在庫量を計算することを考える。 過去の需要系列demand_dfをもとに、Prophetによる予測サンプルを生成し、それに基づいたシミュレーションベースの基在庫レベルの 最適化を行う関数base_stock_prophetを準備する。 引数： - demand_df: 需要の系列を表すデータフレーム。顧客Cust、製品Prod,日時Date、需要量demandの列をもつものと仮定する。 - agg_period: 集約を行う期。例えば１週間を１期とするなら”1w”とする。既定値は1日（“1d”） - forecast_period: 予測を行う未来の期数。既定値は1 TODO: 今は、リード時間等は決め打ちだが、入力値とする。 返値： 以下の情報を列情報にもつデータフレーム - mean: 需要の平均 - std : 需要の標準偏差 - S: 基在庫レベル - cost: 1期間の平均費用 \#hide \### base_stock_prophet関数の使用例 ## 安全在庫配置モデル ### 動的最適化 \#hide ここでは，ネットワークが木状になっている場合について考え，動的最適化による厳密解法を示す． 木ネットワークモデルは，閉路を含まない有向グラフ *G* = (*N*, *A*) 上で義される． ここでは，*G* を無向グラフにしたものが木であると仮定する． 点集合 *N* は在庫地点を表し，枝 (*i*, *j*) ∈ *A* が存在するとき，在庫地点 *i* は在庫地点 *j* に補充 を行うことを表す．複数の在庫地点から補充を受ける点においては，補充された各々の品目を用いて 別の品目を生産すると考える． 点 *i* における品目の生産時間は，*T**i* 日である． *T**i* には各段階における生産時間の他に，待ち時間および輸送時間も含めて考える． このとき点 *j* は，複数の在庫地点から補充を受けるので，点 *j* が品目を発注してから，すべての品目が揃うまで生産を開始できない． 点上で生産される品目は，点によって唯一に定まるものと仮定する． このとき，点と品目を表す添え字は同一であると考えられるので， 有向グラフ *G* = (*N*, *A*) は部品展開表を意味することになる． (*i*, *j*) ∈ *A* のとき，点 *j* 上の品目 *j* は，点 *i* から補充される品目 *i* をもとに生産される． 品目 *j* を生産するのに必要な品目 *i* の量を *ϕ**i**j* と記す． 需要は，後続する点をもたない点 *j* 上で発生するものとし， その 1 日あたりの需要量は，期待値 *μ**j* の定常な分布をもつものとする． 点 *j* における *t* 日間における需要の最大値を *D**j*(*t*) 記す． 直接の需要をもたない点（後続点をもつ点） *i* に対する需要量の期待値 *μ**i* は， *μ**i* = ∑(*i*, *j*) ∈ *A**ϕ**i**j**μ**j* と計算できる． 点 *i* に対する*t* 日間における需要の最大値 *D**i*(*t*) は， *D**i*(*t*) = *t**μ**i* + (∑(*i*, *j*) ∈ *A**ϕ**i**j*(*D**j*(*t*) − *t**μ**j*)*p*)1/*p* と計算されるものと仮定する． ここで，*p*( ≥ 1) は定数である． *p* が大きいほど，需要の逆相関が強いことを表す． 点 *j* 上での需要が常に同期していると仮定した場合には *p* = 1 が， 点 *j* 上での需要が独立な正規分布にしたがうと仮定した場合には *p* = 2 が推奨される． 点 *i* の保証リード時間を *L**i* と記す． ここでは，各点 *i* に対して，保証リード時間の下限は 0 とし，上限は *L̄**i* であるとする． 点 *j* が品目を発注してから，すべての品目が 揃うまでの時間（日数）を入庫リード時間とよび， *L**I**j* と記す． 点 *j* における入庫リード時間 *L**I**j* は，以下の式を満たす． *L**i* ≤ *L**I**j*   ∀(*i*, *j*) ∈ *A* 入庫リード時間 *L**I**i* に生産時間 *T**i* を加えたものが，補充の指示を行ってから生産を完了するまでの時間となる． これを，補充リード時間とよぶ． 補充リード時間から、保証リード時間 *L**i* を減じた時間内の最大需要に相当する在庫を保持していれば，在庫切れの心配がないことになる． 補充リード時間から *L**i* を減じた時間（*L**i* + 1 + *T**i* − *L**i*）を正味補充時間とよぶ． 点 *i* における安全在庫量 *I**i* は，正味補充時間内における最大需要量から平均需要量を減じた量であるので， *I**i* = *D*(*L**I**i* + *T**i* − *L**i*) − (*L**I**i* + *T**i* − *L**i*)*μ**i* となる． 記述を簡単にするために，点 *i* において，保証リード時間が *L*， 入庫リード時間が *L**I* のときの安全在庫費用を表す，以下の関数 *H**C**i*(*L*, *L**I*) を導入しておく． *H**C**i*(*L*, *L**I*) = *h**i*{*D*(*L**I* + *T**i* − *L*) − (*L**I* + *T**i* − *L*)*μ**i*} 上の記号を用いると，木ネットワークモデルにおける安全在庫配置問題は，以下のように定式化できる． $$ \begin{array}{l l l } minimize & \sum\_{i \in N} HC_i (L_i,LI_i) & \\ s.t. & L_i \leq LI\_{i}+T_i & \forall i \in N \\ & L_i \leq LI_j & \forall (i,j) \in A \\ & 0 \leq L_i \leq \bar{L}\_i & \forall i \in N \end{array} $$ ここで，最初の制約式は，正味補充時間が非負であることを表す． 木ネットワークモデルにおける安全在庫配置問題に対する動的最適化アルゴリズムは，以下のように構成できる． 与えられた有向グラフ *G* = (*N*, *A*) を無向グラフにしたものを *Ḡ* = (*N*, *E*) と記す． 動的計画法は，グラフ *Ḡ* の葉側から順次再帰方程式を計算していく． 以下では，点集合 *N* は上のアルゴリズムによって得られた順序で 1, 2, ⋯, |*N*| と番号がついているものと仮定する． 各点 *k* ∈ *N* に対して，点集合 {1, 2, ⋯, *k*} の中で，点 *k* から *Ḡ* 上で到達可能な点の集合 を *N**k* と記す． 点集合 *N**k* が下流（需要地点側）の点に対して保証リード時間 *L* で補充するときの 最小費用を *f**k*(*L*)， 点集合 *N**k* が上流（供給地点側）の点から入庫リード時間 *L**I* で補充を受けるときの 最小費用を *g**k*(*L**I*) と書く． ここで，*f**k*(*L*) は，保証リード時間 *L* に対して非増加関数になっているとは限らない． これを非増加関数にするために，保証リード時間 *L* 「以下」で補充するときの最小費用 *f**k*−(*L*) を導入する． 関数 *f**k*−(*L*) は，以下のように計算される． *f**k*−(*L*) = min*T**k* − *L* ≤ *x* ≤ *L**f**k*(*x*) ここで，下限 *T**k* − *L* は，正味補充時間 *x* + *T**k* − *L* が非負であることから導かれたものである． 同様に，関数 *g**k*(*L**I*) を，入庫リード時間 *L**I* に対して非減少関数としたものを *g**k*−(*L**I*) とする． 関数 *g**k*−(*L**I*) は，以下のように計算される． *g**k*−(*L**I*) = min*L**I* ≤ *y* ≤ *L**I* + *T**k**g**k*(*y*) ここで，上限 *L**I* + *T**k* は，正味補充時間 *L**I* + *T**k* − *y* が非負であることから導かれたものである． 点 *k* における保証リード時間が *L*，入庫リード時間が *L**I* のときの 点集合 *N**k* 全体での最小費用を *c**k*(*L*, *L**I*) と書く． このとき，*c*, *f*, *g*, *f*−, *g*− 間の再帰方程式は，以下のようになる． *c**k*(*L*, *L**I*) = *H**C**k*(*L*, *L**I*) + ∑(*i*, *k*) ∈ *A*, *i* \< *k**f**i*−(*L**I*) + ∑(*k*, *j*) ∈ *A*, *j* \< *k**g**j*−(*L*) *f**k*(*L*) = min*L**I**c**k*(*L*, *L**I*) *g**k*(*L**I*) = min*L**c**k*(*L*, *L**I*) *f**k*−(*L*) = min*T**k* − *L* ≤ *x* ≤ *L**f**k*(*x*) *g**k*−(*L**I*) = min*L**I* ≤ *y* ≤ *L**I* + *T**k**g**k*(*y*) 上の動的最適化アルゴリズムの計算量は，補充リード時間の最大値を *L*max としたとき，*O*(|*N*|*L*max2) となる． このモデルの拡張として，枝 (*i*, *j*) 上での輸送時間 *τ**i**j* を考慮することが考えられる． これは，上の再帰方程式を *c**k*(*L*, *L**I*) = *H**C**k*(*L*, *L**I*) + ∑(*i*, *k*) ∈ *A*, *i* \< *k**f**i*−(*L**I* − *τ**i**k*) + ∑(*k*, *j*) ∈ *A*, *j* \< *k**g**j*−(*L* + *τ**k**j*) と変更するだけで良い．もしくは、ダミーの点を点間に挿入することによっても処理できる。 上のアルゴリズムは，入力サイズを |*N*| としたときには，*L*max = ∑*i* ∈ *N**T**i* であるので，擬多項式時間である． ### 最大需要の計算と安全在庫費用の計算 引数： - G: 有向グラフ - ProcTime: 処理時間 - LTLB: 保証リード時間の下限 - LTUB: 保証リード時間の上限 - z: 安全在庫係数 - mu: 需要の平均 - sigma: 需要の標準偏差 - h: 在庫費用 返値： - Lmax: 正味補充時間（在庫日数）の上限 - MaxDemand: t (=0..Lmax) 日間の最大需要量 - SafetyCost: t (=0..Lmax) 日間の安全在庫費用 ------------------------------------------------------------------------ source ### max_demand_compute ``` python def max_demand_compute( G, ProcTime, LTLB, LTUB, z, mu, sigma, h ): ``` *computing the max demand for t days and the safety stock cost when the net reprenishment time is t (t=0,…,MaxNRT(i))* This function is used in the safety stock allocation problem. Current verstion assumes that the demand has a normal distribution. retuen Lmax (maximum net reprenishment time +1), MaxDemand, SafetyCost ### 最大需要計算関数 max_demand_compute の適用例 ``` python 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) ``` Lmax= 11 MaxDemand= [[ 0. 319.8 628.00142853 934.29460599 1239.6 1544.27414595 1848.49989691 2152.38587596 2456.00285707 2759.4 3062.61309767] [ 0. 216.5 423.33452378 628.57883832 833. 1036.89512163 1240.41658076 1443.65489663 1646.66904756 1849.5 2052.17758139] [ 0. 124.75 235.00178567 342.86825749 449.5 555.34268244 660.62487113 765.48234495 870.00357134 974.25 1078.26637209]] SafetyCost= [[ 0. 19.8 28.00142853 34.29460599 39.6 44.27414595 48.49989691 52.38587596 56.00285707 59.4 62.61309767] [ 0. 82.5 116.6726189 142.89419162 165. 184.47560814 202.08290378 218.27448316 233.34523779 247.5 260.88790696] [ 0. 49.5 70.00357134 85.73651497 99. 110.68536489 121.24974227 130.9646899 140.00714267 148.5 156.53274418]] ### 安全在庫配置モデルに対する動的最適化アルゴリズム 引数： - G: 閉路を含まない有向グラフ（ただし無向グラフにしたものが木である必要がある。） - ProcTime: 処理時間（NumPyの配列） - LTLB: 保証リード時間の下限（NumPyの配列） - LTUB: 保証リード時間の上限（NumPyの配列） - z: 安全在庫係数 - mu: 需要の平均（NumPyの配列） - sigma: 需要の標準偏差（NumPyの配列） - h: 在庫費用（NumPyの配列） 返値： - total_cost: 最適費用 - Lstar: 最適保証リード時間 - NRT: 最適正味補充時間（在庫日数） ------------------------------------------------------------------------ source ### dynamic_programming_for_SSA ``` python def dynamic_programming_for_SSA( G, ProcTime, LTLB, LTUB, z, mu, sigma, h ): ``` *Solve the safety stock allocation problem by the dynamic programming algorithm.* The network should be a tree. The algorithm is based on: @incollection{Graves2003, author={S. C. Graves and S. Willems}, title={Supply Chain Design: Safety Stock Placement and Supply Chain Configuratation}, year={2003}, volume={11}, pages={95–132}, editor={A. G. {de} Kok and S.C. Graves}, publisher={Elsevier}, series={Handbook in Operations Research and Management Science}, chapter={3}, booktitle={Supply Chain Management: {D}esign, Coordination and Operation}, annote={ } } ### dynamic_programming_for_SSA 関数の適用例 - ３点の例題 - 10点の直列多段階モデル ``` python 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) ``` TotalCost= 295.0106609291552 Lstar defaultdict(.>, {0: 0, 1: 1, 2: 2}) NRT defaultdict(.>, {0: 5, 1: 4, 2: 3}) ``` python 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) ``` TotalCost= 2029.231689581059 Lstar defaultdict(.>, {9: 0, 8: 0, 7: 0, 6: 0, 5: 0, 4: 0, 3: 0, 2: 0, 1: 0, 0: 0}) NRT defaultdict(.>, {9: 5, 8: 5, 7: 5, 6: 5, 5: 5, 4: 5, 3: 5, 2: 5, 1: 5, 0: 5}) ### 安全在庫配置モデルに対するメタ解法 \#hide より実践的な方法として、一般ネットワークに対応したメタ解法が考えられる。以下では、 @book{book, author = {Minner, Stefan}, year = {2000}, month = {01}, pages = {}, title = {Strategic Safety Stocks in Supply Chains}, volume = {490}, doi = {10.1007/978-3-642-58296-7}, series = {Lecture Notes in Economics and Mathematical Systems}, publisher ={Springer} } に基づいたタブー探索によるメタ解法を示す。NumPyで近傍を高速に計算している点が特徴である。 準備として、安全在庫配置問題の定式化を，正味補充時間を *N**T**R**i* (net replenishment time)を変数とした定式化に変換することを考える。 *L**i* = max*k* : (*k*, *i*) ∈ *A*{*L**k*} + *T**i* − *N**R**T**i* を元の(保証リード時間と入庫リード時間を変数とした)定式化に代入する。 全ての供給地点（入次数が 0 の点）から点 *i* に至るパスの集合を *P**a**t**h**T**o**i*、全ての供給地点から需要地点（出次数が 0 の点）に至る パスの集合を *P**A**T**H* とする。パス *p* に含まれる点の集合を *N**p* とすると、以下の定式化を得る。 $$ \begin{array}{l l l } minimize & \sum\_{i \in N} \sum h_i \sigma_i \sqrt{NRT_i} & \\ s.t. & \max\_{p \in PathTo_i} \sum\_{i \in N_p} (T_i - NRT_i ) \geq 0 & \forall i \in N \\ & \sum\_{i \in N_p} (NRT_i - T_i ) \geq 0 & \forall p \in PATH \\ & NRT_i \geq 0 & \forall i \in N \end{array} $$ 点 *i* への最大入庫リード時間は、 max*p* ∈ *P**a**t**h**T**o**i*∑*i* ∈ *N**p* \\ {*i*}(*T**i* − *N**R**T**i*) であるので、それに生産時間 *T**i* を加えたものが、最大補充リード時間 *M**a**x**L**I**i* になる。 また、点 *i* から需要地点に至るパスの集合を *P**a**t**h**F**r**o**m**i* と記すと、点 *i* の最小保証リード時間 *M**i**n**L**T**i* は、 min*p* ∈ *P**a**t**h**F**r**o**m**i*∑*i* ∈ *N**p* \\ {*i*}(*N**R**T**i* − *T**i*) と書くことができる。 点 *i* の正味補充時間 *N**T**R**i* は、補充リード時間から保証リード時間を減じたもの（と0の大きい方）であるので、 上流の地点からの最大補充リード時間 *M**a**x**L**I**i* と下流の地点への最小保証リード時間 *M**i**n**L**T**i* を用いて、 以下のように計算できる。 *N**R**T**i* = \[*M**a**x**L**I**i* − *M**i**n**L**T**i*\]+ 点ごとに、正の在庫を配置するか否かを表すベクトル *b* を与えたとき、最大補充リード時間 *M**a**x**L**I**i* は下流の地点から順に以下のように計算できる。 $$ MaxLI_i = \left\\ \begin{array}{ l l} T_i & \forall i \in 供給地点 \\ T_i + \max\_{ (k,i) \in A } (1-b_k) MaxLI_k & それ以外 \end{array} \right. $$ これは、上流の地点が在庫地点の場合には、入庫リード時間が 0 になり、在庫地点でない場合には、正味補充時間が 0 になることを利用している。 また、最小保証リード時間 *M**i**n**L**T**i* は上流の地点から順に以下のように計算できる。 $$ MinLT_i = \left\\ \begin{array}{ l l} LTUB_i & \forall i \in 需要地点 \\ \min\_{ (i,j) \in A } NRT_j + MinLT_j - T_j & それ以外 \end{array} \right. $$ ここで， *L**T**U**B**i* は地点 *i* に対する保証リード時間の上限である． 同時に、正味補充時間 *N**R**T**i* も以下のように計算する。 *N**R**T**i* = \[*M**a**x**L**I**i* − *M**i**n**L**T**i*\]+ 全ての地点の正味補充時間が計算されたら、総費用は、安全在庫係数 *z*、需要の標準偏差 *σ**i*、在庫費用 *h**i* を用いて、以下のように評価できる。 $$ \sum h_i z \sigma_i \sqrt{ NRT_i } $$ 近傍や探索方法については、種々のメタ解法が考えられるが、以下ではNumPyを用いて一度に評価値を計算するタブー探索を構築する。 また、複数製品に対して同時に評価値を計算することも可能である。 決定変数 *b* は、需要地点に対しては除外することができる。 *b* が使われるのは、 *M**a**x**L**I* の計算の部分だけであり、 需要地点には後続する点がないため、需要地点に対する *b* は解の評価に影響を与えないためである。 上の方法では、供給地点、需要地点のリード時間は 0 と仮定しているが、ダミーの点を追加したり、生産時間を変更することによって、保証リード時間の上下限付きの問題も解ける。 用語と記号： - 保証リード時間　(guaranteer lead time, service time): *L* - 生産時間 (production time): *T* - 入庫リード時間 (incoming lead time): *L**I* - 補充リード時間 (replenishment lead time) - 正味補充時間 (net replenishment time): *N**R**T* - 最大補充リード時間： *M**a**x**L**I* - 最小保証リード時間： *M**i**n**L**T* ### 安全在庫配置モデルに対するタブー探索アルゴリズム 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: 最良解における最小保証リード時間 ------------------------------------------------------------------------ source ### tabu_search_for_SSA ``` python def tabu_search_for_SSA( G, ProcTime, LTUB, z, mu, sigma, h, max_iter:int=100, TLLB:int=1, TLUB:int=10, seed:int=1 ): ``` *一般のネットワークの安全在庫配置モデルに対するタブー探索（mu未使用；NRT日の最大在庫量をシミュレーションもしくは畳み込みによって事前計算？もしくは正規分布で近似）* ### tabu_search_for_SSA関数の使用例 - ３点の例題 - 10点の直列多段階モデル - ベンチマーク問題例 ``` python #3点の例題 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) ``` 最良値 119.50357133746822 最良解 [0 1 1] 正味補充時間 [0. 9. 8.] 最大補充リード時間 [5. 9. 8.] 最小保証リード時間 [5. 0. 0.] ``` python #4点の直列多段階の例題 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) ``` 最良値 973.2050807568877 最良解 [1 0 0 1] 正味補充時間 [4. 0. 0. 3.] 最大補充リード時間 [4. 3. 2. 3.] 最小保証リード時間 [0. 3. 2. 0.] ![](03inventory_files/figure-commonmark/cell-93-output-2.png) ``` python #cost, stage_df, I, cost_list, phi_list = periodic_inv_opt_fit_one_cycle(stage_df, bom_df, max_iter = 100, n_samples = 10, # n_periods = 100, seed = 1, lr_find=False, max_lr =0.1, moms=(0.85,0.95)) ``` ``` python #fig = plot_simulation(stage_df, I, n_periods =100, samples=5) #plotly.offline.plot(fig); ``` ``` python #10点の直列多段階モデル 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) ``` 最良値 1905.4467494843116 最良解 [1 0 1 0 0 0 0 1] 正味補充時間 [10. 0. 25. 0. 0. 0. 0. 5.] 最大補充リード時間 [10. 5. 25. 20. 15. 10. 5. 5.] 最小保証リード時間 [ 0. 5. 0. 20. 15. 10. 5. 0.] ``` python #7点の一般型ネットワークの例題 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) ``` 最良値 514.8330943986502 最良解 [1 1 0 0 1 0 1] 正味補充時間 [6. 2. 0. 0. 6. 0. 2.] 最大補充リード時間 [6. 2. 3. 3. 6. 3. 3.] 最小保証リード時間 [0. 0. 3. 3. 0. 3. 1.] ``` python 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 ```
Unnamed: 0 name net_replenishment_time max_guaranteed_LT processing_time replenishment_LT guaranteed_LT z average_demand sigma h b capacity x y
0 0 0 6.0 0 6 6.0 0.0 1.65 200 14.1 1 8.8 2000.0 0 0
1 1 1 2.0 0 2 2.0 0.0 1.65 200 14.1 1 8.8 2000.0 0 2
2 2 2 0.0 0 3 3.0 3.0 1.65 200 14.1 3 26.3 2000.0 1 0
3 3 3 0.0 0 3 3.0 3.0 1.65 200 14.1 1 8.8 2000.0 0 1
4 4 4 6.0 0 3 6.0 0.0 1.65 200 14.1 5 43.9 2000.0 2 0
5 5 5 0.0 3 3 3.0 3.0 1.65 100 10.0 6 52.7 1000.0 3 0
6 6 6 2.0 1 3 3.0 1.0 1.65 100 10.0 6 52.7 1000.0 3 1
``` python 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)) ``` ELT= [21. 17. 13. 13. 11. 1. 3.] S= [4306.61362354 3495.92405238 2683.88315042 2683.88315042 2277.16127575 116.5 328.57883832] t: Cost |dC| S 0: 16082.26, 2384121.79078 [4306.61362354 3495.92405238 2683.88315042 2683.88315042 2277.16127575 116.5 328.57883832] 50: 16064.14, 55008642.72020 [4304.02992277 3498.69251131 2679.18513171 2686.70283345 2274.05919829 119.66406743 344.12554439] 100: 16065.15, 5026167.93997 [4305.81174737 3498.76846583 2678.40095109 2686.61752507 2273.55957572 120.84683667 340.51778989] 150: 16061.73, 109999436.60082 [4303.16094465 3502.71988907 2672.66577621 2685.7481414 2279.55090957 120.56158991 341.11925034] 200: 16065.26, 5222011.53967 [4312.76590192 3492.32576853 2681.96668542 2684.31474986 2272.03751252 120.58780939 340.94018907] 250: 16062.70, 28978730.42492 [4312.2184163 3492.26969664 2681.9034483 2682.24835951 2272.76515243 120.68572406 341.01899329] 300: 16083.57, 11181822.94933 [4311.15839432 3493.40225594 2680.87064495 2703.05431135 2273.45772924 120.83500766 341.05489796] 350: 16079.60, 5321.32095 [4307.69066274 3497.77039769 2676.08815218 2700.18398616 2277.70718962 120.48910017 340.87407204] 400: 16075.63, 63551.03463 [4307.64834193 3497.78911077 2676.32680439 2696.23356042 2277.16193522 120.79860651 341.06201143] 450: 16071.07, 11328.58722 [4308.05463524 3497.72018746 2676.22034128 2691.43264714 2276.99734823 120.79266494 340.98458816] 500: 16065.01, 29646.92335 [4307.82339303 3497.72585609 2676.17683586 2685.7119275 2277.17212027 120.96624478 341.06602026] 550: 16058.02, 105022.66152 [4307.27597593 3497.85808953 2676.04251307 2679.3697765 2277.57701891 120.59324681 341.0362919 ] 600: 16126.80, 78916.21887 [4307.84910082 3497.47376626 2676.02237295 2747.63728074 2276.2455183 120.83557065 340.82623601] 650: 16128.17, 75923.78927 [4307.12100365 3497.49710856 2675.84794371 2749.80816501 2277.04754176 120.91189266 340.91285357] 700: 16126.23, 105024.49925 [4306.96412636 3497.52463691 2675.65866046 2748.08164634 2277.39564944 119.95046856 341.00884094] 750: 16125.26, 25019.64204 [4307.6174081 3497.33181103 2675.81217298 2746.63589639 2276.7266567 120.89359202 340.92325567] 800: 16123.20, 107188.28867 [4306.65730431 3497.56396823 2675.5592795 2745.52972634 2277.63620685 120.55485679 340.81960078] 850: 16122.84, 105773.51829 [4307.15776745 3497.487481 2675.42919782 2744.77972635 2277.49292265 120.53375714 341.01864899] 900: 16122.64, 7777.36243 [4307.44336725 3497.45996464 2675.33777062 2744.35799326 2277.41409729 120.54642398 340.86674134] 950: 16122.86, 14175.35550 [4307.9300071 3497.38395884 2675.23658045 2744.191592 2277.26791974 120.55581911 340.85063693] Best: 16023.43, 0.00965 [4307.50553981 3497.77958792 2676.18044936 2656.18886109 2277.33810851 121.21002819 341.04645203] ``` python fig = plot_simulation(stage_df, I, n_periods =100, samples=5) plotly.offline.plot(fig); ``` ``` python 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* 座標（整数の文字列） ------------------------------------------------------------------------ source ### read_willems ``` python def read_willems( file_name ): ``` *Willemsのベンチマーク問題例の読み込み関数* ### read_wellems関数の使用例 ``` python G, pos = read_willems(file_name="../data/bom/03.xml") nx.draw(G, pos=pos, with_labels=True, node_color="Yellow") ``` ![](03inventory_files/figure-commonmark/cell-103-output-1.png) ``` python #Pyvisの利用 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メソッドを用いてアクセスできる。 ``` python for i in G: print( G.nodes[i]["stageTime"], end=" " ) print( nx.get_node_attributes(G, "maxServiceTime") ) ``` 10 10 28 15 10 0 0 0 {'Retail_0001': '0', 'Retail_0002': '0', 'Retail_0003': '0'} ### 在庫費用、需要量、需要の標準偏差の計算 ベンチマークデータからタブー探索で必要なデータを抽出し、必要なデータを準備する。 各地点には付加された費用と生産時間、需要地点には需要と標準偏差と保証リー時間上限（サービス時間）が入っている。 すべての点に対して、在庫費用と需要・標準偏差など、安全在庫配置モデルの計算に必要なデータを計算する関数を準備しておく。 引数： - G: ベンチマーク問題例のグラフ 返値： - ProcTime: 処理時間 （NumPy配列） - LTUB: 保証リード時間上限 （NumPy配列） - z: 安全在庫係数 （NumPy配列） - mu: 需要の平均 （NumPy配列） - sigma: 需要の標準偏差 （NumPy配列） - h: 在庫費用 （NumPy配列） ------------------------------------------------------------------------ source ### extract_data_for_SSA ``` python def extract_data_for_SSA( G ): ``` *ベンチマークデータからタブー探索で必要なデータを抽出する関数* ### extract_data_for_SSA関数とタブー探索の実行例 ``` python 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) ``` cost= 17431180.517088737 [1 1 1 1 1 1 1 1 1 0 0 1 0 1 1 0 0] 最良値 14635043.250984032 最良解 [1 1 1 1 0 1 0 0 0 0 1 0 1 1 0 0 0] 正味補充時間 [ 3. 13.2 16.3 16.2 0. 47.5 0. 0. 0. 0. 16. 0. 31. 57. 0. 0. 0. ] 最大補充リード時間 [ 3. 13.2 16.3 16.2 57. 47.5 12. 12. 45. 37.5 53.5 26. 31. 59. 2. 2. 2. ] 最小保証リード時間 [ 0. 0. 0. 0. 57. 0. 12. 12. 45. 37.5 37.5 37.5 0. 2. 2. 2. 2. ] ``` python #在庫の有無を色で区別し、正味補充時間（在庫量）をノードの大きさで表して描画 nx.draw(G, pos=pos, with_labels=True, node_size=best_NRT*50,node_color=best_sol) ``` ![](03inventory_files/figure-commonmark/cell-108-output-1.png) ``` python print("best obj.=", best_cost) print("best sol.=", best_sol) print("best NRT=", best_NRT) ``` best obj.= 14635043.250984032 best sol.= [1 1 1 1 0 1 0 0 0 0 1 0 1 1 0 0 0] best NRT= [ 3. 13.2 16.3 16.2 0. 47.5 0. 0. 0. 0. 16. 0. 31. 57. 0. 0. 0. ] ### 安全在庫配置問題の結果を描画する関数 draw_graph_for_SSA 点の大きさが在庫量（正味補充時間）を表し、点の色（色の濃いものほど小さい）が保証リード時間を表す。 引数： - G: 閉路を含まない有向グラフ - pos: 点の座標を保管した辞書 - best_NRT: 最良解における正味補充時間（在庫日数） - best_MaxLI: 最良解における最大補充リード時間 - best_MinLT: 最良解における最小保証リード時間 返値： - Plotlyの図オブジェクト ------------------------------------------------------------------------ source ### draw_graph_for_SSA ``` python def draw_graph_for_SSA( G, pos, best_NRT, best_MaxLI, best_MinLT ): ``` *安全在庫配置問題の結果をPlotlyで描画する関数* ### draw_graph_for_SSAの使用例 ``` python 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: 在庫を複数の下流の点に配分するときの割合（在庫シミュレーション最適化で用いる．） ------------------------------------------------------------------------ source ### make_df_for_SSA ``` python def make_df_for_SSA( G, pos, ProcTime, LTUB, z, mu, sigma, h, best_NRT, best_MaxLI, best_MinLT ): ``` *点(stage)と枝(bom)を保管したデータフレームを返す関数* ### make_df_for_SSAの使用例 ``` python 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 ```
name net_replenishment_time max_guaranteed_LT processing_time replenishment_LT guaranteed_LT z average_demand sigma h b capacity x y
0 Manuf_0001 10.0 0.0 10.0 10.0 0.0 1.644854 298.0 36.633651 65.0 565.2 2980.0 176 32
1 Manuf_0002 10.0 0.0 10.0 10.0 0.0 1.644854 120.0 2.236068 62.0 539.2 1200.0 176 96
2 Part_0001 28.0 0.0 28.0 28.0 0.0 1.644854 418.0 36.701831 12.0 104.4 4180.0 32 32
3 Part_0002 15.0 0.0 15.0 15.0 0.0 1.644854 418.0 36.701831 5.0 43.5 4180.0 32 96
4 Part_0003 10.0 0.0 10.0 10.0 0.0 1.644854 418.0 36.701831 9.0 78.3 4180.0 32 160
5 Retail_0001 0.0 0.0 0.0 0.0 0.0 1.644854 253.0 36.620000 65.0 565.2 2530.0 304 32
6 Retail_0002 0.0 0.0 0.0 0.0 0.0 1.644854 45.0 1.000000 127.0 1104.4 450.0 304 96
7 Retail_0003 0.0 0.0 0.0 0.0 0.0 1.644854 75.0 2.000000 62.0 539.2 750.0 304 160
``` python #bom_df.to_csv(folder_bom + "ssa_bom02.csv") bom_df.head() ```
child parent units allocation
0 Manuf_0001 Retail_0001 1 0.5
1 Manuf_0001 Retail_0002 1 0.5
2 Manuf_0002 Retail_0002 1 0.5
3 Manuf_0002 Retail_0003 1 0.5
4 Part_0001 Manuf_0001 1 0.5
### データフレームからデータの図を描画する関数 draw_graph_for_SSA_from_df 引数： - 点(stage)の情報を保管したデータフレーム - 枝（部品展開表；bom）の情報を保管したデータフレーム 返値： - fig : Plotlyの図オブジェクト - G : グラフ - pos : 点の座標を表す辞書 ------------------------------------------------------------------------ source ### draw_graph_for_SSA_from_df ``` python def draw_graph_for_SSA_from_df( stage_df, bom_df ): ``` *データフレームを入れると、安全在庫配置問題のデータをPlotlyで描画する関数* ### draw_graph_for_SSA_from_df関数の使用例 ``` python fig, G, pos = draw_graph_for_SSA_from_df(stage_df, bom_df) #plotly.offline.plot(fig) nx.draw(G, pos=pos) ``` ![](03inventory_files/figure-commonmark/cell-116-output-1.png) ### 全てのベンチマークデータをcsvファイルに変換 ``` python 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") ``` ### データフレームを入力すると、一般のネットワークの安全在庫配置モデルに対するタブー探索を行う関数　solve_SSA 引数： - stsge_df : 点の情報を保管したデータフレーム - bom_df : 枝の情報を保管したデータフレーム 返値： - stage_df : （最適化された後の）点の情報を保管したデータフレーム - bom_df : （最適化された後の）枝の情報を保管したデータフレーム - fig : Plotlyの図オブジェクト ------------------------------------------------------------------------ source ### solve_SSA ``` python def solve_SSA( stage_df, bom_df ): ``` *データフレームを入力とした一般のネットワークの安全在庫配置モデルに対するタブー探索* ### solve_SSA関数の使用例 ``` python 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 ```
Unnamed: 0 name net_replenishment_time max_guaranteed_LT processing_time replenishment_LT guaranteed_LT z average_demand sigma h b capacity x y
0 0 Dist_0001 3.0 0.0 1.0 3.0 0.0 1.644854 70.0 35.000000 2750.0 23913.9 700.0 384 272
1 1 Dist_0002 13.2 0.0 1.2 13.2 0.0 1.644854 126.0 132.300000 4103.0 35679.6 1260.0 480 64
2 2 Dist_0003 16.3 0.0 4.3 16.3 0.0 1.644854 57.0 51.300000 4162.0 36192.6 570.0 480 144
3 3 Dist_0004 16.2 0.0 4.2 16.2 0.0 1.644854 46.0 45.000000 4247.4 36935.3 460.0 480 240
4 4 Manuf_0001 0.0 0.0 12.0 57.0 57.0 1.644854 126.0 132.300000 2400.0 20870.3 1260.0 192 32
5 5 Manuf_0002 47.5 0.0 10.0 47.5 0.0 1.644854 173.0 76.692177 2350.0 20435.5 1730.0 192 272
6 6 Manuf_0003 0.0 0.0 10.0 12.0 12.0 1.644854 126.0 132.300000 3953.0 34375.2 1260.0 384 64
7 7 Manuf_0004 0.0 0.0 10.0 12.0 12.0 1.644854 103.0 68.239944 3912.0 34018.6 1030.0 384 192
8 8 Part_0001 0.0 0.0 45.0 45.0 45.0 1.644854 126.0 132.300000 1100.0 9565.6 1260.0 96 32
9 9 Part_0002 0.0 0.0 37.5 37.5 37.5 1.644854 299.0 152.921483 600.0 5217.6 2990.0 96 112
10 10 Part_0003 16.0 0.0 53.5 53.5 37.5 1.644854 299.0 152.921483 400.0 3478.4 2990.0 96 192
11 11 Part_0004 0.0 0.0 26.0 26.0 37.5 1.644854 173.0 76.692177 1100.0 9565.6 1730.0 96 272
12 12 Part_0005 31.0 0.0 31.0 31.0 0.0 1.644854 229.0 148.862285 1500.0 13044.0 2290.0 192 144
13 13 Trans_0001 57.0 0.0 2.0 59.0 2.0 1.644854 126.0 132.300000 2400.0 20870.3 1260.0 288 32
14 14 Trans_0002 0.0 0.0 2.0 2.0 2.0 1.644854 126.0 132.300000 1503.0 13070.0 1260.0 288 112
15 15 Trans_0003 0.0 0.0 2.0 2.0 2.0 1.644854 103.0 68.239944 1502.0 13061.3 1030.0 288 192
16 16 Trans_0004 0.0 0.0 2.0 2.0 2.0 1.644854 173.0 76.692177 2350.0 20435.5 1730.0 288 272
![](03inventory_files/figure-commonmark/cell-120-output-1.png) ## データフレームをもとに定期発注方策最適化を行う関数　periodic_inv_opt 定期発注方策の仮定の下では，保証リード時間は 0 とみなす必要位がある． また，各点のリード時間は，安全在庫配置モデルの補充リード時間 (replenishment leadtime: best_MaxLI) に （定期発注方策なので） 1 を加えたもの相当する． ただし出荷を待っている（すでに下流からの注文がある；保証リード時間に相当する）在庫も含めた在庫レベルを考えるものとする． 保証リード時間中の需要は確定値なので考慮する必要はない． 保証リード時間 +1 日をリード時間と考える． 最適な (*s*, *S*) を近似的に求め，シミュレーションを行う． 安全在庫配置の結果がすでに入ったデータフレームをもとに，定期発注方策の最適化を行う． 需要が負になることが頻繁に発生する場合には，基在庫レベルを多少大きくする必要があり，どの程度大きくすれば良いかは実験を行う必要がある． ------------------------------------------------------------------------ source ### periodic_inv_opt ``` python def periodic_inv_opt( stage_df, bom_df, max_iter:int=1, n_samples:int=10, n_periods:int=100, seed:int=1 ): ``` *データフレームを入力とした定期発注方策最適化* ## periodic_inv_opt関数の使用例 ``` python 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) ``` ELT= [51. 3. 3. 3. 3. 1. 1. 1. 1. 19. 19. 19. 19.] S= [1348319.56810962 34742.07917634 16612.97028229 51290.49501161 53718.70272284 17130.14553439 8127.40124938 23756.17464127 3150.96049975 133232.97749839 64708.72283258 220679.57299807 234805.54429147] t: Cost |dC| S 0: 220675718.04, 341808445420091.93750 [1348320.56810962 34743.07917634 16613.97028229 51289.49501161 53717.70272284 17131.14553439 8128.40124938 23757.17464127 3149.96049975 133231.97749839 64707.72283258 220680.57299807 234806.54429147] 50: 220518415.34, 745192019079272.87500 [1348369.74940509 34792.85221809 16664.8282914 51239.63768253 53667.72090384 17181.14553428 8178.40124932 23807.17464124 3099.97095921 133182.24568996 64657.06116374 220730.44691115 234856.53155207] 100: 220354712.04, 545116807024090.43750 [1348426.76713248 34846.91481195 16716.94689668 51187.45789964 53616.50062463 17231.14553417 8228.40124925 23857.17464121 3050.044785 133127.36879553 64603.34998017 220781.00090685 234905.83987621] 150: 220191360.04, 816167149409793.25000 [1348486.82370908 34903.48503105 16766.10557845 51133.51485836 53565.50242318 17277.19121652 8278.40124919 23907.17464117 3000.22634602 133069.85025368 64550.7454711 220832.68490125 234953.49107123] 200: 220031267.65, 524126465993651.37500 [1348542.95615425 34955.9460657 16815.80172057 51081.56513545 53514.01504129 17313.78270815 8327.85942526 23957.17398411 2950.52739226 133016.37778189 64498.76830765 220882.89683515 235002.88942521] 250: 219863018.77, 851656363689291.25000 [1348612.21089618 35006.21452572 16862.70239911 51027.50135256 53460.33050642 17351.19918629 8369.84554242 24007.15874206 2900.91843682 132960.78488283 64444.2078238 220930.75367586 235049.94351966] 300: 219700343.56, 491721644568249.31250 [1348676.41473432 35054.97028877 16909.60862722 50974.98279207 53408.28295656 17373.56656534 8411.73101577 24056.33644305 2851.52246168 132907.26589133 64390.73173766 220977.72557878 235095.97078381] 350: 219542527.16, 484599717516724.68750 [1348736.35266099 35102.34503693 16956.66288597 50924.01894217 53357.83028046 17396.07668408 8445.93033377 24105.04489955 2802.22025908 132855.793889 64338.46719159 221024.0100242 235141.25257648] 400: 219385962.74, 941456836707623.00000 [1348794.9972844 35150.88708824 17003.98943326 50872.34184194 53306.37285049 17420.11896177 8480.91118128 24152.88404331 2753.03918282 132803.81058383 64286.38134399 221071.74828293 235188.42507887] 450: 219231581.50, 1292884402927184.75000 [1348853.50496323 35195.70244227 17051.45155343 50820.34418959 53254.74315809 17445.57594705 8513.21949266 24200.89941181 2703.97199643 132754.33492422 64234.18980277 221119.93761815 235235.8310807 ] 500: 219077632.32, 518212719275614.12500 [1348911.93773799 35240.93227256 17097.03219103 50768.58804249 53203.15333844 17472.34115796 8538.87860671 24249.15353218 2654.8591634 132704.45854483 64183.09232032 221167.71096212 235283.11500071] 550: 218927945.80, 485300108160533.00000 [1348966.0166179 35285.31383805 17144.33879554 50718.35296326 53153.23104983 17500.32266857 8565.58203945 24297.48613848 2605.84288308 132656.52583678 64132.59678114 221215.37673921 235330.12264657] 600: 218785108.88, 564728180060374.50000 [1349018.55399569 35322.15998852 17187.18826913 50674.21697735 53110.22096048 17529.43982656 8593.41530842 24345.721259 2556.96370265 132614.06722426 64085.44475664 221254.63665692 235367.40000979] 650: 218646284.66, 231028727702383.43750 [1349070.90351672 35352.84581936 17226.14782567 50634.78730121 53072.46039344 17559.62092361 8622.31320316 24394.12631278 2508.00635225 132575.53572489 64040.66882231 221286.60962038 235396.3918021 ] 700: 218506805.66, 261632579367210.09375 [1349123.51945598 35385.14764584 17265.96791165 50594.24685425 53033.32750908 17590.80151025 8652.21613975 24442.68290959 2459.09382769 132535.61299023 63995.17601281 221319.80658156 235426.92666981] 750: 218369502.68, 256688191996597.75000 [1349175.01467117 35416.09715995 17305.43544488 50554.70182601 52995.33893496 17622.92314284 8683.06987149 24491.37060831 2410.33081751 132496.76185324 63950.27649986 221351.79321731 235455.93472848] 800: 218238726.74, 76741322279162.53125 [1349222.7074572 35442.901349 17340.43378657 50519.26270029 52961.25523777 17655.93243016 8714.82462737 24540.17268668 2361.73199689 132461.77924135 63909.55459394 221379.28820056 235480.69272982] 850: 218123687.56, 252671648888372.81250 [1349258.70671942 35462.56084521 17362.58570206 50493.34870657 52935.24328265 17689.78029324 8747.43444326 24589.07531308 2313.27083062 132435.63810722 63881.27524148 221398.5689691 235499.65316048] 900: 218009137.59, 68610489926442.88281 [1349294.39416997 35482.04231764 17384.15586794 50467.40216142 52909.2094108 17724.42138067 8780.85663644 24638.06691778 2264.76814395 132409.64951475 63853.40192051 221417.95847274 235518.71131709] 950: 217894318.98, 195265332282450.40625 [1349331.44805033 35501.08900369 17405.30090711 50440.96572054 52882.59493555 17759.81359948 8815.05138677 24687.13772632 2216.40572571 132383.31962837 63825.11154321 221437.14500597 235537.67686744] Best: 217782429.64, 149495319804164.03125 [1349367.47399352 35517.50819437 17426.82285353 50416.49339415 52858.10594517 17794.46036713 8848.57052732 24732.04221619 2170.10604192 132358.93389869 63796.96827159 221453.67157211 235553.69765306] ``` python 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) ``` ELT= [13. 13. 43. 30. 25. 1. 1. 1.] S= [ 4091.25968171 1573.2612339 18369.86692342 12870.65539681 10751.84569932 313.23453982 46.64485363 78.28970725] t: Cost |dC| S 0: 738253.04, 1246180.05146 [ 4092.25968171 1574.2612339 18370.86692342 12871.65539681 10752.84569932 314.23453982 45.64485363 77.28970725] 50: 713359.98, 66398.11736 [ 4126.42972125 1582.6009784 18420.72971972 12921.39842509 10802.46119396 360.61652664 44.05167729 86.44553083] 100: 697639.86, 53830.07941 [ 4091.67521411 1584.22271517 18470.1237238 12968.54931585 10851.5716812 397.82352518 44.38684851 96.95976312] 150: 685425.19, 48819.45366 [ 4079.83481905 1586.63059694 18512.88671884 13015.14479448 10903.8863555 446.18484549 44.49041005 101.36175361] 200: 671502.46, 35828.72853 [ 4105.09263065 1588.59693178 18557.05478751 13057.59583735 10948.05285451 492.24195436 44.68381728 91.43024867] 250: 659084.98, 26829.23061 [ 4139.92605168 1589.98668442 18599.73444809 13096.62095046 10985.33379269 523.59614578 44.91317741 91.86145222] 300: 650148.26, 19159.74856 [ 4165.5160091 1591.75562153 18637.47297218 13130.09797031 11022.13602295 546.41056921 45.03762899 93.43725329] 350: 643304.81, 15687.78590 [ 4189.55033689 1592.70953476 18674.42229426 13160.04183158 11051.55120041 566.08940204 45.18650575 92.97595812] 400: 638532.00, 10794.45977 [ 4204.45349226 1593.43260336 18709.22350723 13191.90346841 11074.28284291 577.87592218 45.32749954 94.76898651] 450: 635303.54, 7460.30899 [ 4219.73677263 1594.94109142 18739.05818407 13218.34505493 11090.14206376 585.36319824 45.50282587 96.32977209] 500: 634971.64, 6489.75384 [ 4205.15995699 1594.66460962 18764.41560121 13242.9950394 11104.76389076 595.0605535 45.52592229 94.76459883] 550: 632963.55, 3295.00148 [ 4217.08463465 1596.33267547 18784.72278404 13262.46129873 11130.39204887 600.38535139 45.66693841 97.64674733] 600: 631952.77, 2317.41725 [ 4224.15853318 1596.71048501 18802.49007556 13277.91652818 11149.06447296 603.45915015 45.78075323 97.78152212] 650: 631286.34, 1530.61933 [ 4229.94756846 1597.52954626 18815.89877426 13292.54637619 11166.39287595 605.57515493 45.90637947 96.65280796] 700: 631066.34, 1517.05126 [ 4230.92904393 1597.63591293 18827.5361712 13307.72669112 11177.41762992 606.05643722 45.90453356 95.79483507] 750: 631333.14, 992.71457 [ 4226.21435871 1597.88033477 18840.39992139 13322.96803148 11187.32570124 606.11049996 45.94290883 95.9827774 ] 800: 630968.29, 666.99091 [ 4233.53446195 1598.39470045 18851.41302688 13334.88998105 11197.63796615 606.12822126 46.00327799 96.49321949] 850: 630994.24, 535.97391 [ 4235.19800183 1598.66085629 18862.17833695 13345.50332417 11205.21368173 606.3251161 46.08802725 96.52798901] 900: 631246.78, 650.50730 [ 4233.23452117 1598.47753639 18872.58430379 13354.65422521 11213.35212038 606.41638563 46.08940061 96.72492291] 950: 631469.82, 426.99272 [ 4232.59142704 1598.93457249 18881.19858988 13363.216574 11227.65568071 606.53896245 46.0892974 96.96500643] Best: 631419.44, 363.69285 [ 4235.8734129 1599.09495873 18886.92290645 13373.24500118 11236.88720748 606.63677805 46.14271796 96.99789479] ------------------------------------------------------------------------ source ### periodic_inv_opt_fit_one_cycle ``` python def periodic_inv_opt_fit_one_cycle( stage_df, bom_df, max_iter:int=1, n_samples:int=10, n_periods:int=100, seed:int=1, lr_find:bool=False, max_lr:float=1.0, moms:tuple=(0.85, 0.95) ): ``` *データフレームを入力とした定期発注方策最適化にfit_one_cycle法を適用* ### periodic_inv_opt_fit_one_cycle関数の使用例 ``` python 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)) ``` ELT= [13. 13. 43. 30. 25. 1. 1. 1.] S= [ 4091.25968171 1573.2612339 18369.86692342 12870.65539681 10751.84569932 313.23453982 46.64485363 78.28970725] t: Cost |dC| S 0: 738253.04, 1246180.05146 [ 4091.25968171 1573.2612339 18369.86692342 12870.65539681 10751.84569932 313.23453982 46.64485363 78.28970725] 50: 729962.63, 89221.59767 [ 4103.85176968 1579.8573079 18382.68942081 12883.48247112 10764.70974338 325.54292727 43.20516334 84.22070204] 100: 717473.79, 64478.15708 [ 4122.22831209 1582.4721447 18409.8506545 12910.6163079 10791.94966657 350.22769019 44.20632408 86.84664873] 150: 701812.48, 58172.65323 [ 4120.19269689 1584.16946547 18451.2273357 12951.53437628 10832.40539252 370.21700749 44.31117571 87.94913286] 200: 690953.13, 60815.35292 [ 4062.70653571 1585.69807974 18501.9607133 13000.27053152 10894.66454178 423.63561966 44.24871301 102.28241792] 250: 669454.27, 35903.22928 [ 4110.0131565 1589.35990371 18561.27325685 13063.44568125 10952.60156984 480.86736111 44.84719815 90.04694133] 300: 651959.38, 22369.98247 [ 4160.03111767 1590.95128037 18628.021709 13124.54412233 11010.7422393 527.18429571 44.96473601 93.08156545] 350: 640868.74, 14428.94592 [ 4188.07897284 1593.45999619 18697.3150612 13179.58213487 11068.59723696 553.86982831 45.29400135 93.93846431] 400: 634387.40, 6546.35564 [ 4219.32210446 1595.22121831 18758.49281169 13234.61171226 11097.53514323 568.46997814 45.53821725 95.20422733] 450: 631948.45, 2488.80843 [ 4225.03604559 1596.63585148 18802.02811649 13275.93793071 11146.58359853 577.4433366 45.77948835 98.02081483] 500: 631081.45, 1775.35928 [ 4230.42993272 1597.67066657 18831.48844513 13309.58154627 11180.65425739 574.25880735 45.90536989 95.88547543] 550: 630970.64, 617.93082 [ 4234.89376666 1598.45718056 18858.40419835 13338.65504692 11202.45543356 571.04654927 46.01808733 97.27355649] 600: 631544.47, 534.24797 [ 4230.62877671 1598.58264321 18881.2790734 13357.60957547 11221.70974986 570.3504795 46.07435167 98.29520824] 650: 631406.22, 337.12891 [ 4238.26664247 1599.30027922 18893.87019139 13375.21210064 11239.7887459 570.65236523 46.18509758 98.43290482] 700: 631465.29, 191.68175 [ 4241.63093738 1599.38044819 18902.08059305 13387.88206348 11249.01060311 571.435331 46.18127829 98.20658432] 750: 631534.94, 127.62068 [ 4243.77805685 1599.62032251 18908.61889896 13396.99837585 11255.72444415 572.09966484 46.25620134 97.39047516] 800: 631600.65, 106.00507 [ 4244.95997097 1599.74597316 18912.93543413 13402.38345209 11260.70871733 572.66442053 46.28662202 97.54121488] 850: 631596.31, 97.50065 [ 4246.6751567 1599.83671964 18915.27856784 13406.08338141 11263.49953519 573.02046203 46.30120998 97.53524439] 900: 631628.51, 99.24852 [ 4246.81121286 1599.8218352 18916.23618929 13408.21981863 11265.13857997 573.24120237 46.2968931 97.53178233] 950: 631648.49, 96.24211 [ 4246.69363253 1599.82581378 18916.60874107 13409.06371862 11265.84744558 573.32375665 46.29486373 97.53590301] Best: 630910.88, 0.00965 [ 4234.58753372 1597.78313105 18826.26660408 13302.34402449 11175.91912062 574.98567419 45.92049737 96.15625735] ``` python fig = plot_simulation(stage_df, I, n_periods =100, samples=5) plotly.offline.plot(fig); ``` ### 学習率探索の結果を可視化する関数 plot_inv_opt_lr_find 引数： - cost_list: 費用の変化を入れたリスト - ss_list: 学習率（ステップサイズ）設定のパラメータの変化を入れたリスト 返値： - Plotlyの図オブジェクト ------------------------------------------------------------------------ source ### plot_inv_opt_lr_find ``` python def plot_inv_opt_lr_find( cost_list, ss_list ): ``` *学習率探索の結果を可視化する関数* ![](03inventory_files/figure-commonmark/cell-128-output-1.png) ### 探索を可視化する関数 plot_inv_opt 引数： - cost_list: 費用の変化を入れたリスト 返値： - fig: 費用の変化を可視化した図オブジェクト ------------------------------------------------------------------------ source ### plot_inv_opt ``` python def plot_inv_opt( cost_list ): ``` *Call self as a function.* ![](03inventory_files/figure-commonmark/cell-130-output-1.png) ### シミュレーション最適化の図表示関数 plot_simulation ------------------------------------------------------------------------ source ### plot_simulation ``` python def plot_simulation( stage_df, I, n_periods:int=1, samples:int=5, stage_id_list:NoneType=None ): ``` *Call self as a function.* ``` python fig = plot_simulation(stage_df, I, n_periods =100, samples=5) plotly.offline.plot(fig); ``` ![](03inventory_files/figure-commonmark/cell-134-output-1.png)