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

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

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

YouTubeVideo("NBmMfXzccr8")

(Q,R)方策

  • LT: リード時間
  • LS: ロットサイズ(Qに相当する)

安全在庫 ($R=S$) を計算し、在庫シミュレーションをNumPyを用いて行う。 正味在庫レベル(実在庫と輸送中在庫の和)が基在庫レベル $S$ を下回ったら、固定数量 $Q$ だけ生産する。

引数:

  • n_samples: シミュレーションを行うサンプル数
  • n_periods: 期間数
  • mu: 需要の平均
  • sigma: 需要の標準偏差
  • LT: リード時間
  • LS: ロットサイズ ($=Q$)
  • z: 安全在庫係数;安全在庫は基在庫レベル $S$ を求めるときに使われる。
  • b: バックオーダー費用
  • h: 在庫費用

返値:

  • total_cost: 総費用の期待値
  • I: 在庫の時間的な変化を表すNumPy配列で、形状は (n_samples, n_periods+1);在庫の可視化に用いる。

simulate_inventory[source]

simulate_inventory(n_samples=1, n_periods=10, mu=100, sigma=10, LT=3, LS=300, z=1.65, b=100, h=1)

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

simulate_inventoryの使用例

cost, I = simulate_inventory(n_samples =10, n_periods =100, z=2.33) 
print("cost=",cost)
pd.DataFrame(I[0]).plot() #在庫の推移の可視化
S= 340.35678381635483
cost= 247.7860866103851
<AxesSubplot:>

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

ここでは,発注を行うことができるタイミングが限定されている定期発注方策を考える. 定期発注方策では,時間を離散値として取り扱う.実際には,日単位などで在庫を管理することが多いため, 定期発注方策は,多くの実際問題で適用可能なモデルである. 時間は,期とよばれる単位に分割されており,期を $1,2,\ldots,t,\ldots,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 \left\{ c, [ s- (I_{t}-D_t + T_{t})]^+ \right\} $$

ここで,$[ \cdot ]^+$ は,$\max \{0, \cdot\}$ を表す記号である.

  • $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)' \mathbf{1}[I_t<0] + h(I_t)' \mathbf{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)。

base_stock_simulation[source]

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

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

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={ } }

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

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

n_samples = 100
n_periods = 10000
mu = 100
sigma = 10
LT = 3
z = 2.43
b, h = 100, 1

capacity = 105
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 - 10*dC
    print(f"{iter_}: {S:.2f} {dC:.3f} {cost:.2f}")
    if dC**2<=convergence:
        break
Initial S= 342.08883462392373
t:   S      dS     Cost
0: 368.04 -2.596 73.98
1: 360.93 0.711 65.62
2: 356.71 0.422 61.47
3: 355.53 0.118 60.30
4: 355.39 0.014 60.22
5: 355.38 0.001 60.22

定期発注方策(多段階在庫モデル)

ここでは,単一段階モデルを直列多段階の場合に拡張する. この結果は,Glasserman―Tayur \cite{Glasserman1995} に基づく. 多くの実際問題は,直列多段階に分解して考えることができるので,このモデルは実務的にも有意義であるが, 以下で解説する一般のネットワークに対する基礎を与えるという意味でも重要である.

点の順番は下流(需要地点側)から $0,1,\ldots,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^i_{t+1} =I^i_{t} -q^{i-1}_{t} +q^i_{t-L_i} $$$$ T^i_{t+1} =T^i_{t} + q^i_{t} - q^i_{t-L_i} $$

となる. ただし,$i=0$ のときには,上式における, 上流の地点からの発注量 $q^{i-1}_{t}$ は需要量 $D_t$ に置き換える.

初期条件は,在庫量については $I^0_0=s^0$, $I^i_0=s^i-s^{i-1}, i=1,2,\ldots,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,\ldots,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^{i-1}_t/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^0_0)'=\mathbf{1}[ i^*=0], $$ $$(I_0^i)'= \mathbf{1}[ i^*=i] - \mathbf{1}[ i^*=i-1 ], i=1,\cdots,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である。

引数:

  • n_samples: シミュレーションを行うサンプル数
  • n_periods: 期間数
  • demand: 需要地点ごとの需要量を保持した辞書。キーは顧客の番号。値は、需要のサンプルを表すNumPyの配列で、形状は (n_samples, n_periods)。
  • capacity: 容量(生産量上限)を表すNumPy配列
  • LT: 地点毎のリード時間を表すNumPy配列
  • b: バックオーダー費用
  • h: 在庫費用
  • S: 基在庫レベル

返値:

  • dC: 各地点における総費用の基準在庫レベルによる微分の期待値
  • total_cost: 総費用の期待値
  • I: 在庫の時間的な変化を表すNumPy配列で、形状はグラフの点の数をn_stagesとしたとき (n_samples, n_stages, n_periods+1)。

multi_stage_base_stock_simulation[source]

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

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

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

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

# 多段階 基在庫レベル方策 在庫シミュレーション
np.random.seed(123)

n_samples = 10
n_stages = 3
n_periods = 1000
mu = 100
sigma = 10
LT = np.array([1, 2, 3])
# 初期エシェロン在庫量は、最終需要地点以外には、定期発注方策の1日分の時間を加えておく。
ELT = LT.cumsum()
ELT = ELT + 1
ELT[0] = ELT[0]-1
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 4 7]
S= [113.3        426.6        735.18849244]
t: Cost    |dC|      S 
0: 9845.73, 9570.201 [113.3499     426.7328     833.01579244]
1: 3754.00, 157.188 [112.7618     429.4649     845.23779244]
2: 3693.16, 48.284 [113.5357     431.512      851.83279244]
3: 3678.30, 20.909 [113.731      433.1613     856.09319244]
4: 3675.36, 13.231 [113.9966     434.1602     859.58069244]
5: 3675.81, 8.473 [114.1214     434.8663     862.40179244]
6: 3677.91, 5.566 [114.0982     435.7484     864.58989244]
7: 3680.46, 3.454 [114.2384     436.46       866.30109244]
8: 3682.41, 2.112 [114.3462     436.9272     867.67309244]
9: 3684.31, 1.474 [114.3998     437.2886     868.83109244]
10: 3686.13, 0.899 [114.487      437.7312     869.66529244]
11: 3687.37, 0.529 [114.5533     438.2183     870.20089244]
12: 3688.27, 0.286 [114.6411     438.5316     870.62479244]
13: 3688.90, 0.218 [114.6325     438.8289     870.98509244]
14: 3689.73, 0.141 [114.7732     438.9792     871.29909244]
15: 3689.96, 0.110 [114.7583     439.1067     871.60449244]
16: 3690.62, 0.086 [114.8388     439.185      871.87469244]
17: 3690.92, 0.071 [114.8078     439.1939     872.13879244]
18: 3691.48, 0.058 [114.8318     439.273      872.36569244]
19: 3691.87, 0.053 [114.8211     439.2968     872.59459244]
20: 3692.32, 0.048 [114.8299     439.3523     872.80529244]
21: 3692.71, 0.046 [114.866      439.3464     873.01709244]
22: 3692.99, 0.074 [114.8513     439.3027     873.28449244]
23: 3693.48, 0.103 [114.8279     439.2164     873.59219244]
24: 3694.04, 0.065 [114.8272     439.1844     873.84489244]
25: 3694.47, 0.042 [114.8069     439.2386     874.04199244]
26: 3694.92, 0.034 [114.8172     439.3403     874.19399244]
27: 3695.24, 0.028 [114.8322     439.4037     874.34659244]
28: 3695.52, 0.024 [114.8353     439.411      874.50119244]
29: 3695.80, 0.024 [114.8274     439.4167     874.65739244]
30: 3696.11, 0.018 [114.84       439.4374     874.78909244]
31: 3696.33, 0.020 [114.8293     439.4488     874.93139244]
32: 3696.63, 0.011 [114.8411     439.4969     875.02549244]
33: 3696.80, 0.011 [114.8364     439.5389     875.12019244]
34: 3697.02, 0.010 [114.8328     439.5895     875.20519244]

Adamによる最適化

深層学習で定番のAdamで最適化をしてみたが、収束は悪い。

# Adam  (slower than a simple code)
np.random.seed(123)

n_samples = 10
n_stages = 3
n_periods = 1000
mu = 100
sigma = 10
LT = np.array([1, 2, 3])
# 初期エシェロン在庫量は、最終需要地点以外には、定期発注方策の1日分の時間を加えておく。
ELT = LT.cumsum()
ELT = ELT + 1
ELT[0] = ELT[0]-1
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

convergence = 1e-5
alpha = 1.
beta_1 = 0.9
beta_2 = 0.999
epsilon = 1e-8
m_t = np.zeros(n_stages)
v_t = np.zeros(n_stages)

for t in range(100):
    dC, cost, I = multi_stage_base_stock_simulation(
        n_samples, n_periods, demand, capacity, LT, b, h, S)
    g_t = dC  # 勾配
    # 移動平均の更新
    m_t = beta_1*m_t + (1-beta_1)*g_t
    # 勾配の二乗の移動平均の更新
    v_t = beta_2*v_t + (1-beta_2)*(g_t**2)
    m_cap = m_t/(1-beta_1**(t+1))
    v_cap = v_t/(1-beta_2**(t+1))
    S=S - (alpha*m_cap)/(np.sqrt(v_cap)+epsilon)
    print(f"{t}: {cost:.2f}", S, dC)
    if np.dot(g_t, g_t) <= convergence:
        break
Echelon LT= [1 4 7]
S= [113.3        426.6        735.18849244]
t: Cost    |dC|      S 
0: 9845.73 [114.2999998  427.59999992 736.18849244] [-4.99000e-02 -1.32800e-01 -9.78273e+01]
1: 9747.78 [115.28600435 428.60128156 737.18846888] [-3.89000e-02 -1.43800e-01 -9.77393e+01]
2: 9649.89 [116.18730905 429.60046128 738.1884141 ] [-1.69000e-02 -1.32800e-01 -9.76733e+01]
3: 9552.14 [117.09778289 430.59830083 739.18828071] [-2.90000e-02 -1.31700e-01 -9.75193e+01]
4: 9454.56 [117.91682207 431.59540818 740.18802839] [-7.00000e-03 -1.31700e-01 -9.73543e+01]
5: 9357.14 [118.43254903 432.59207671 741.18765806] [ 2.2700e-02 -1.3170e-01 -9.7263e+01]
6: 9259.86 [118.63673136 433.58846972 742.18713757] [ 3.370e-02 -1.317e-01 -9.712e+01]
7: 9162.79 [118.66511526 434.58468711 743.18641602] [ 2.50000e-02 -1.31700e-01 -9.69353e+01]
8: 9065.87 [118.52615842 435.58483468 744.18547473] [ 3.04000e-02 -1.42700e-01 -9.67977e+01]
9: 8969.12 [118.07813924 436.58797442 745.18429019] [  0.106   -0.1427 -96.6533]
10: 8872.60 [117.45315845 437.58913069 746.18286932] [  0.2017  -0.1317 -96.562 ]
11: 8776.29 [116.72487951 438.5917337  747.18118585] [  0.3006  -0.1391 -96.4225]
12: 8680.30 [115.91945341 439.59964287 748.1790772 ] [  0.2039  -0.1501 -96.0508]
13: 8584.64 [115.07224771 440.61181186 749.17642679] [  0.1548  -0.1501 -95.7157]
14: 8489.36 [114.22254632 441.62306159 750.17316092] [  0.0967  -0.1391 -95.4266]
15: 8394.49 [113.36035341 442.64148549 751.16919098] [  0.1125  -0.1596 -95.1029]
16: 8300.08 [112.46348613 443.66158281 752.16442796] [  0.1596  -0.1486 -94.765 ]
17: 8206.13 [111.50964095 444.68730551 753.15886995] [  0.294   -0.1596 -94.5474]
18: 8112.74 [110.67136278 445.71328672 754.15222687] [-3.22000e-02 -1.48600e-01 -9.39022e+01]
19: 8019.81 [109.8378124  446.72691189 755.14461355] [  0.0957  -0.1213 -93.8044]
20: 7927.39 [108.93587152 447.72994886 756.13611592] [  0.2702  -0.1215 -93.6927]
21: 7835.60 [107.96709202 448.71804241 757.12661763] [  0.361   -0.1105 -93.3105]
22: 7744.63 [106.96848476 449.69383949 758.11597722] [  0.2359  -0.1118 -92.8761]
23: 7654.47 [106.00508944 450.65834073 759.10407744] [ 8.91000e-02 -1.10700e-01 -9.24554e+01]
24: 7565.15 [105.29724786 451.61302701 760.09075484] [ -0.2128  -0.1107 -91.9555]
25: 7476.54 [104.77429624 452.5527995  761.07598012] [ -0.1674  -0.0997 -91.6159]
26: 7388.61 [104.70716386 453.46573979 762.05941071] [-5.98700e-01 -7.77000e-02 -9.08546e+01]
27: 7300.81 [104.97346519 454.3619933  763.04082171] [-7.12800e-01 -8.87000e-02 -9.02015e+01]
28: 7213.19 [105.40506626 455.24823321 764.02021017] [ -0.4923  -0.096  -89.8097]
29: 7126.09 [106.035865   456.14676495 764.99703547] [ -0.8141  -0.1308 -88.7271]
30: 7039.48 [106.80644659 457.05624498 765.97114402] [ -0.7811  -0.1308 -88.0671]
31: 6953.51 [107.69795433 457.97548575 766.94227028] [ -0.9163  -0.1308 -87.2499]
32: 6868.36 [108.66315274 458.90961422 767.91026665] [ -0.7266  -0.1418 -86.5486]
33: 6784.23 [109.64150025 459.85067936 768.87504911] [ -0.4186  -0.1308 -85.9106]
34: 6701.30 [110.59259025 460.7780012  769.83645531] [ -0.2341  -0.0978 -85.1711]
35: 6619.56 [111.46197023 461.68603662 770.79445954] [-1.16000e-02 -8.68000e-02 -8.45796e+01]
36: 6538.89 [112.19850608 462.56948875 771.74904012] [ 1.83500e-01 -7.58000e-02 -8.39937e+01]
37: 6459.34 [112.76806036 463.43856262 772.7000379 ] [  0.3282  -0.0868 -83.2474]
38: 6381.11 [113.15471804 464.27145051 773.64720509] [ 4.31900e-01 -5.38000e-02 -8.23831e+01]
39: 6304.10 [113.35158443 465.05497263 774.59043925] [ 5.24100e-01 -3.18000e-02 -8.16613e+01]
40: 6228.28 [113.46804896 465.76701656 775.52922993] [ 2.21100e-01  1.20000e-03 -8.04673e+01]
41: 6153.74 [113.51201643 466.46627735 776.46313595] [ 2.18500e-01 -6.41000e-02 -7.92994e+01]
42: 6080.78 [113.5012389  467.16264797 777.39163127] [ 1.78500e-01 -7.51000e-02 -7.79944e+01]
43: 6009.36 [113.50411313 467.87313714 778.31431924] [-4.42000e-02 -9.71000e-02 -7.67817e+01]
44: 5939.22 [113.47840584 468.6176261  779.2312854 ] [  0.0978  -0.1266 -76.0582]
45: 5870.40 [113.34697508 469.36207105 780.14263114] [  0.3752  -0.0842 -75.366 ]
46: 5803.02 [113.13427513 470.10030627 781.0481063 ] [  0.3265  -0.0757 -74.3138]
47: 5737.03 [112.80901965 470.92291128 781.94755244] [  0.4741  -0.208  -73.3391]
48: 5672.46 [112.50529879 471.84122795 782.84052614] [ 2.8000e-02 -2.5300e-01 -7.2056e+01]
49: 5609.31 [112.30890347 472.82863995 783.72669007] [ -0.2687  -0.231  -70.8463]
50: 5547.67 [112.08922982 473.86557253 784.60612774] [  0.142   -0.2131 -70.0429]
51: 5487.68 [111.84617934 474.94474454 785.47865639] [  0.1488  -0.2131 -68.9827]
52: 5429.13 [111.89121364 476.05441939 786.34337885] [ -0.9042  -0.2021 -67.1927]
53: 5371.63 [112.22649195 477.16677563 787.19985018] [ -1.1614  -0.1581 -65.7805]
54: 5315.35 [112.72334311 478.30942206 788.04805181] [ -0.822   -0.2115 -64.7575]
55: 5260.46 [113.24055692 479.47782044 788.88817258] [ -0.2719  -0.2115 -63.9436]
56: 5207.16 [113.76921929 480.65274826 789.72005358] [ -0.2412  -0.1811 -62.8167]
57: 5155.30 [114.25213467 481.82125284 790.54365994] [-1.17000e-02 -1.59100e-01 -6.18032e+01]
58: 5104.96 [114.69111578 482.977677   791.35889375] [-2.40000e-03 -1.48100e-01 -6.07345e+01]
59: 5055.96 [115.12460575 484.13580966 792.16554989] [ -0.1378  -0.1709 -59.5643]
60: 5008.37 [115.55532376 485.2949408  792.96319406] [ -0.1457  -0.1705 -58.1708]
61: 4962.38 [116.0519975  486.46058028 793.75136596] [ -0.424   -0.1815 -56.7265]
62: 4917.75 [116.60449295 487.6265179  794.5299519 ] [ -0.41   -0.172 -55.573]
63: 4874.45 [117.19701123 488.79297175 795.29901527] [ -0.3669  -0.173  -54.5811]
64: 4832.51 [117.81873049 489.94741451 796.05850447] [ -0.3376  -0.151  -53.4994]
65: 4791.84 [118.51023648 491.06507472 796.80822818] [ -0.5274  -0.107  -52.2976]
66: 4752.63 [119.19865975 492.08435908 797.54835964] [-2.44400e-01 -8.00000e-03 -5.14146e+01]
67: 4714.76 [119.88319218 493.01437393 798.27879453] [-2.38900e-01 -8.00000e-03 -5.03091e+01]
68: 4678.27 [120.60735719 493.8707797  798.9991389 ] [-4.20200e-01 -1.84000e-02 -4.89514e+01]
69: 4643.11 [121.3518048  494.64464233 799.70939459] [-3.54200e-01  3.60000e-03 -4.79174e+01]
70: 4609.21 [122.11325431 495.34358068 800.40950786] [-3.48400e-01  3.60000e-03 -4.68452e+01]
71: 4576.63 [122.84633572 495.96638916 801.0995204 ] [-1.67400e-01  1.46000e-02 -4.58602e+01]
72: 4545.24 [123.50232056 496.52006529 801.7797531 ] [ 3.21000e-02  1.46000e-02 -4.51247e+01]
73: 4514.85 [124.0974967  497.00268956 802.45024978] [-2.40000e-03  2.56000e-02 -4.41882e+01]
74: 4485.53 [124.65646468 497.40448443 803.11081072] [ -0.0739   0.047  -43.0601]
75: 4457.26 [125.25084517 497.73284834 803.76116859] [ -0.3393   0.047  -41.8707]
76: 4430.10 [125.8594928  497.99449385 804.40144916] [ -0.2683   0.047  -40.9957]
77: 4404.04 [126.5128657  498.18703668 805.03125828] [ -0.3958   0.058  -39.7132]
78: 4379.18 [127.22612142 498.31680392 805.65045502] [ -0.4786   0.058  -38.6184]
79: 4355.37 [127.90026052 498.38875824 806.25960122] [ -0.1056   0.059  -38.0904]
80: 4332.50 [128.48895127 498.40817377 806.85893747] [  0.0806   0.059  -37.3526]
81: 4310.48 [128.99794201 498.3712925  807.44843493] [  0.0884   0.07   -36.4254]
82: 4289.27 [129.52313709 498.2832776  808.0277009 ] [ -0.24    0.07  -35.217]
83: 4268.98 [130.11910959 498.14880539 808.59640512] [ -0.4599   0.07   -34.0401]
84: 4249.71 [130.74542346 497.96356179 809.15474128] [ -0.3281   0.081  -33.2589]
85: 4231.37 [131.36624853 497.72371696 809.70265704] [ -0.2003   0.092  -32.3087]
86: 4213.88 [131.96147538 497.4429035  810.24038609] [ -0.1219   0.081  -31.5841]
87: 4197.13 [132.5488356  497.12580813 810.76801062] [ -0.1784   0.0798 -30.7674]
88: 4181.12 [133.13274835 496.77845725 811.28559111] [ -0.1914   0.0763 -29.9479]
89: 4165.83 [133.73283778 496.40356914 811.79309881] [ -0.2632   0.0763 -29.0731]
90: 4151.20 [134.33619605 496.0036041  812.29085484] [ -0.2199   0.0763 -28.4674]
91: 4137.23 [134.99271695 495.5894259  812.77879859] [ -0.4112   0.0653 -27.6051]
92: 4123.90 [135.68946308 495.16220825 813.25726473] [ -0.3822   0.0653 -27.0401]
93: 4111.16 [136.4740279  494.74769329 813.72618204] [ -0.5931   0.0344 -26.2043]
94: 4099.06 [137.31210291 494.33567211 814.18580631] [ -0.4931   0.0454 -25.6113]
95: 4087.60 [138.05439255 493.92575039 814.63704018] [  0.0581   0.0454 -25.5135]
96: 4076.59 [138.61605377 493.49095391 815.08022256] [  0.3845   0.0784 -25.0699]
97: 4065.87 [138.99458031 493.03287357 815.51574508] [  0.4565   0.0792 -24.6917]
98: 4055.40 [139.16308319 492.55359233 815.94394652] [  0.6206   0.0792 -24.3058]
99: 4045.15 [139.13783592 492.05531565 816.36487984] [  0.6515   0.0788 -23.7423]
#可視化(最初のサンプルを表示)
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)
#fig.show()
Image("../figure/multi_stage_base_stock.png")

定期発注方策(一般ネットワークモデル)

実際には,組立工程があったり,複数の地点から発注を受けたるする場合がある. 以下では,一般のネットワークの場合に,上の手法を拡張することを考える.

ある地点 $i$ が複数の地点から発注を受け,さらにその地点自身でも需要が発生する場合を考える. すべての需要(発注)に応えるだけの在庫量がない場合に, 複数の需要(発注)に対して,どのように在庫を割り振るかを決める必要がある. これを分配ルール(allocation rule)とよぶ. 代表的な分配ルールには,以下のものがある.

  • 優先順位法: 直接需要ならびに(複数の)下流の在庫地点に対して,優先順位を付けておき,その順に需要(発注)量を満たす.
  • 発注量比例配分法: 需要(発注)量に応じて比例配分をする. この方法だと,供給配分に起因する鞭効果が発生し,安全在庫量が増えてしまう危険性がある.
  • 実績比例配分法: 過去の実績に応じて,在庫の比例配分を行う.この方法だと,ある顧客の過去の実績が多く,かつ今回の発注分が少ない場合に,過剰な配分をしてしまい,(この顧客用に確保された)在庫が残っているにもかかわらず,他の顧客の需要を満たすことができないという理不尽な現象が起きる危険性がある.
  • 最適配分法: 線形計画問題を解くことによって,過去の実績から計算された優先度と発注量の積の合計を最大化するように配分する.

以下では,実績比例配分法を例として解説する.

ネットワークを表現するための記号を導入しておく.

  • $PARENT_i$: ネットワークにおける点 $i$ の直後の点の集合であり,点 $i$ が直接供給を行う地点の集合を表す.
  • $ANCESTOR_{i}$: 点 $i$ からネットワークの有向枝を辿って到達可能な点に対応する点集合. 点 $i$ 自身も含めた点集合と定義する.
  • $CHILD_i$: ネットワークにおける点 $i$ の直前の点の集合であり,点 $i$ が発注を行う点の集合を表す。 $CHILD_i$ に発注したすべての部品が揃わないと生産できないものと仮定する。

  • $\phi_{ij}$: 点 $j$ に対応する品目を $1$単位製造するために必要な点 $i$ の品目の量. これは,ネットワークが部品展開表を表すときに用いられるパラメータであり,親品目 $j$ を製造するために必要 な子品目 $i$ の量を,枝 $(i,j)$ に付加したものである. 品目の量を表す単位として,$i$-unit,$j$-unitを導入しておくと分かりやすい. ちなみに,$\phi_{ij}$ の単位(次元)は,$[ i$-units $/$ $j$-unit $]$ である.

上流の地点が存在しない,言い換えれば $CHILD_i$ が空の点に対しては,その上流の地点に十分な 在庫があるものとし,品切れはしないものとする.

$t$ 期の期首在庫量が与えられている状態で,需要が発生する.需要はすべての点で発生する可能性があるものとする. (つまり,部品に対しても需要があるものと仮定する.)点 $i$ における $t$期の需要を $D_t^i$ と記す. 期末に,割当ルールによって,下流の地点に対する在庫の配分を行い, 次いで発注量を決定する.これによって,次期の期首在庫 $I_{t+1}^i$ が定まる.

在庫地点 $i$ の在庫地点の $t$ 期における正味在庫量を $I_t^i$ と記す. この記号を用いると,エシェロン基在庫方策は,エシェロン基在庫レベル $s^i$ を与えたとき, $t$ 期末のエシェロン在庫ポジション $$ E_i^t = \sum_{j \in ANCESTOR_i } \rho_{ij} ( I_t^j +T_t^j -D_t^j ) $$ を $s_i$ にするように発注を行う方策と定義できる.

複数の下流の在庫地点からの注文に対して, 現在の在庫で対応できない場合には,過去の実績に応じて,比例配分を行うものとする.

在庫地点 $i$ の在庫を下流の在庫地点 $j$ の配分するときの比率(配分比)を $\alpha_{ij}$ とする. また,在庫地点 $i$ が直接の需要 $D_t^i$ を持っている場合には, その需要に対する配分比を $\alpha_{i0}$ とする.

配分比は, $$ \alpha_{i0}+ \sum_{j \in PARENT_i} \alpha_{ij} =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 - \sum_{j \in PARENT_i} \phi_{ij} 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 -\sum_{j \in PARENT_i} \phi_{ij} 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$ におけるエシェロンリード時間 $ELT_i$, 累積需要の標準偏差(の推定値) $\sigma_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である。

引数:

  • 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)。

network_base_stock_simulation[source]

network_base_stock_simulation(G, n_samples, n_periods, demand, capacity, LT, ELT, b, h, S, phi, alpha)

ネットワーク型定期発注方策に対するシミュレーションと微分値の計算

初期基在庫レベルとエシェロンリード時間を計算する関数

定期発注方策を仮定したときの、エシェロンリード時間を計算し、それをもとに基在庫レベルを計算して返す関数。 エシェロンリード時間は、最終需要地点以外の点に対しては、1日を加えて計算する。

initial_base_stock_level[source]

initial_base_stock_level(G, LT, mu, z, sigma)

初期基在庫レベルとエシェロンリード時間の計算

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

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

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.18, 17.49920 [1847.62848441  217.6321625   238.93345817]
1: 3005.20, 5.65867 [1846.81417691  217.57742     241.16786567]
2: 2999.98, 2.76635 [1846.04694566  217.54595125  242.64323692]
3: 2997.20, 1.49623 [1845.32856066  217.53493625  243.63320192]
4: 2995.52, 0.94987 [1844.65297441  217.5231225   244.33556817]
5: 2994.34, 0.69698 [1844.00798316  217.51431375  244.86555942]
6: 2993.42, 0.51765 [1843.39660566  217.50129125  245.24463692]
7: 2992.68, 0.42572 [1842.80979878  217.49309813  245.52978379]
8: 2992.04, 0.36018 [1842.24721566  217.48693125  245.73868692]
9: 2991.48, 0.30672 [1841.71369003  217.47690688  245.88691254]
10: 2990.99, 0.27152 [1841.20141691  217.47608     245.98226567]
11: 2990.54, 0.24467 [1840.70978378  217.47141313  246.03657879]
12: 2990.13, 0.22343 [1840.23792941  217.4658675   246.06397317]
13: 2989.76, 0.20697 [1839.78303441  217.4622625   246.06930817]
14: 2989.41, 0.19432 [1839.34245441  217.4558925   246.05624817]
15: 2989.08, 0.17402 [1838.92677566  217.44352125  246.02338692]
16: 2988.78, 0.16365 [1838.52434441  217.4399525   245.98237817]
17: 2988.49, 0.15234 [1838.13728566  217.44096125  245.93211692]
18: 2988.22, 0.14172 [1837.76518003  217.44066688  245.87506254]
19: 2987.97, 0.13010 [1837.40935316  217.43984375  245.81600942]
20: 2987.73, 0.12253 [1837.06290972  217.44013719  245.76593286]
21: 2987.50, 0.11261 [1836.73264753  217.43369938  245.70681504]
22: 2987.29, 0.10671 [1836.41113784  217.43005906  245.64914473]
23: 2987.09, 0.10097 [1836.09910816  217.43268875  245.58913442]
24: 2986.89, 0.09478 [1835.79712034  217.43052656  245.52928223]
25: 2986.70, 0.08938 [1835.50444503  217.42955188  245.46827754]
26: 2986.53, 0.08300 [1835.22155347  217.43139344  245.41380911]
27: 2986.36, 0.07914 [1834.94601816  217.42902875  245.35708442]
28: 2986.19, 0.07502 [1834.67817534  217.42842156  245.29980723]
29: 2986.03, 0.06897 [1834.42177191  217.422675    245.24333067]
30: 2985.89, 0.06463 [1834.17314659  217.42490031  245.19029598]
31: 2985.75, 0.05896 [1833.93761503  217.42558188  245.13130754]
32: 2985.61, 0.05554 [1833.70886847  217.42262844  245.07471411]
33: 2985.49, 0.05119 [1833.48938566  217.41671125  245.02007692]
34: 2985.37, 0.04757 [1833.27696784  217.41527906  244.97061473]
35: 2985.26, 0.04456 [1833.07092003  217.41062688  244.92498254]
36: 2985.15, 0.04212 [1832.87020222  217.41009469  244.88218036]
37: 2985.04, 0.03920 [1832.67669753  217.40759938  244.84032504]
38: 2984.93, 0.03415 [1832.49888878  217.39780813  244.79093379]
39: 2984.85, 0.03053 [1832.32974691  217.3951      244.74715567]
40: 2984.76, 0.02796 [1832.16752753  217.39281938  244.70665504]
41: 2984.68, 0.02601 [1832.01042941  217.3906675   244.67021317]
42: 2984.60, 0.02378 [1831.85920316  217.39209375  244.64013942]
43: 2984.53, 0.02228 [1831.71226066  217.39448625  244.61394192]
44: 2984.45, 0.02147 [1831.56797316  217.40052375  244.58916942]
45: 2984.37, 0.02022 [1831.42782566  217.40452125  244.56535692]
46: 2984.30, 0.01794 [1831.29711941  217.4064275   244.53618317]
47: 2984.23, 0.01676 [1831.17113816  217.40465875  244.50636442]
48: 2984.17, 0.01507 [1831.05197691  217.40657     244.47686567]
49: 2984.12, 0.01395 [1830.93749222  217.40695469  244.44779036]
50: 2984.07, 0.01243 [1830.82989378  217.40255313  244.41896879]
51: 2984.02, 0.01130 [1830.72906034  217.40083656  244.38542223]
52: 2983.98, 0.01038 [1830.63123878  217.39925813  244.35702379]
53: 2983.94, 0.00925 [1830.53806722  217.39722969  244.33315536]
#可視化(最初のサンプルを表示)
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)
#fig.show()
Image("../figure/multi_stage_base_stock_network.png")

顧客の需要予測と安全在庫

上では、定常なランダム需要を仮定して、安全在庫量の最適化を行なったが、実際問題においては需要は定常とは限らない。

過去の需要系列demand_dfをもとに、Prophetによる予測を行う関数 forecastを準備する。 比較として、(需要の定常性を仮定して)平均と標準偏差も算出する。

引数:

  • demand_df: 需要の系列を表すデータフレーム。顧客Cust、製品Prod,日時Date、需要量demandの列をもつものと仮定する。
  • agg_period: 集約を行う期。例えば1週間を1期とするなら"1w"とする。既定値は1日("1d")
  • forecast_period: 予測を行う未来の期数。既定値は1

返値:

以下の情報を列情報にもつデータフレーム

  • mean: 需要の平均
  • std : 需要の標準偏差
  • 0..forecast_period-1: Prophetによる予測から計算された未来のforecast_periods 期分の最大在庫量

forecast[source]

forecast(demand_df, agg_period='1d', forecast_periods=1)

forecast関数の使用例

1週間を1期として、30期を予測。

非常に時間がかかるので、注意。

# ProphetのWarningがたくさん出る場合には、以下を生かす。
if False: #実験時にはTrueにする。
    import warnings
    warnings.simplefilter('ignore')
    demand_df = pd.read_csv(folder+"demand.csv")
    forecast_df = forecast(demand_df, agg_period ="1w", forecast_periods =30)
    forecast_df.to_csv(folder+"forecast_yhat.csv")
    forecast_df.head()

Prophetによる予測と需要の平均・標準偏差を用いた安全在庫量の比較と可視化

Prophetによる可視化

forecast_df = pd.read_csv(folder+"forecast.csv",index_col=0)
forecast_df.head()
cust prod mean std 0 1 2 3 4 5 ... 20 21 22 23 24 25 26 27 28 29
0 那覇市 F 1.365385 1.188653 3.080633 3.069456 3.050913 3.045405 3.044113 3.219213 ... 3.565478 3.644279 3.682303 3.781509 3.845952 3.952716 3.931932 4.123730 4.380080 4.413510
1 那覇市 J 1.269231 1.254254 2.778413 2.882829 2.818556 2.776102 2.759031 2.879232 ... 4.970052 4.962857 5.160201 5.585258 5.624556 5.843050 6.193982 6.326028 6.624619 6.851527
2 那覇市 A 1.519231 1.128775 2.886841 2.910441 2.805579 2.835632 2.948129 2.905958 ... 4.588084 5.084178 5.221664 5.184114 5.669602 5.741231 6.039508 6.066876 6.382044 6.709195
3 那覇市 G 3.711538 2.538616 11.859874 12.478805 11.984558 12.825155 12.427594 12.216680 ... 11.792271 12.131251 12.408490 12.546387 12.383458 12.020243 12.712764 11.540315 12.457729 12.247556
4 那覇市 D 1.288462 1.226130 4.006041 3.950352 3.969163 3.909864 3.985566 3.904234 ... 5.521809 5.525722 5.905974 6.300353 6.622292 6.643114 7.051544 7.078838 7.130381 7.580375

5 rows × 34 columns

# 1期分の安全在庫量を安全在庫係数1.65(5%の品切れ率)で計算
forecast_df["Stationary"] = forecast_df["mean"] + 1.65*forecast_df["std"]
# 定常時の安全在庫量とProhetによる1期分の最大需要量を1つの列に変換
inv_pivot = pd.pivot_table(forecast_df, index="cust", values=["Stationary","0"], aggfunc=sum)
inv_stack = pd.DataFrame(inv_pivot.stack(), columns=["demand"])
inv_stack.reset_index(inplace=True)
#fig = px.bar(inv_stack, x='Cust', y="demand",color='level_1',barmode="group")
#fig.show()
#Prophetの最大需要量を顧客別に(製品は積み上げ棒グラフとして)可視化
fig = px.bar(forecast_df, x='cust', y="0",color='prod')
#fig.show()
#定常を仮定した場合の最大需要量を顧客別に(製品は積み上げ棒グラフとして)可視化
fig = px.bar(forecast_df, x='cust', y='Stationary',color='prod')
#fig.show()
Image("../figure/max_demand_prophet.png")
Image("../figure/max_demand_stationary.png")

非定常な需要予測をもとにした基在庫レベルの最適化

以下では、実際問題を想定して、非定常な需要をもとに、安全在庫量を計算することを考える。

過去の需要系列demand_dfをもとに、Prophetによる予測サンプルを生成し、それに基づいたシミュレーションベースの基在庫レベルの 最適化を行う関数base_stock_prophetを準備する。

引数:

  • demand_df: 需要の系列を表すデータフレーム。顧客Cust、製品Prod,日時Date、需要量demandの列をもつものと仮定する。
  • agg_period: 集約を行う期。例えば1週間を1期とするなら"1w"とする。既定値は1日("1d")
  • forecast_period: 予測を行う未来の期数。既定値は1

TODO: 今は、リード時間等は決め打ちだが、入力値とする。

返値:

以下の情報を列情報にもつデータフレーム

  • mean: 需要の平均
  • std : 需要の標準偏差
  • S: 基在庫レベル
  • cost: 1期間の平均費用

base_stock_prophet[source]

base_stock_prophet(demand_df, agg_period='1d', forecast_periods=1)

base_stock_prophet関数の使用例

# ProphetのWarningがたくさん出る場合には、以下を生かす。
import warnings
warnings.simplefilter('ignore')
if False: #実験時にはTrueにする。
    demand_df = pd.read_csv(folder+"demand.csv")
    base_stock_df = base_stock_prophet(demand_df, agg_period ="1w", forecast_periods =10)
    base_stock_df.head()
# 1期分の安全在庫量を安全在庫係数1.65(5%の品切れ率)で計算
#base_stock_df.to_csv(folder+"base_stock.csv")
base_stock_df = pd.read_csv(folder + "base_stock.csv")
base_stock_df["stationary"] = base_stock_df["mean"] + 1.65*base_stock_df["std"]
# 定常時の安全在庫量とProhetによる1期分の最大需要量を1つの列に変換
inv_pivot = pd.pivot_table(base_stock_df, index="cust", values=["stationary","S"], aggfunc=sum)
inv_stack = pd.DataFrame(inv_pivot.stack(), columns=["demand"])
inv_stack.reset_index(inplace=True)
fig = px.bar(inv_stack, x='cust', y="demand",color='level_1',barmode="group")
#fig.show()
Image("../figure/base_stock_vs_prophet.png")
#Prophetの最大需要量を顧客別に(製品は積み上げ棒グラフとして)可視化
fig = px.bar(base_stock_df, x='cust', y="S",color='prod')
#fig.show()
#定常を仮定した場合の最大需要量を顧客別に(製品は積み上げ棒グラフとして)可視化
fig = px.bar(base_stock_df, x='cust', y='stationary',color='prod')
#fig.show()

安全在庫配置モデル

動的最適化

ここでは,ネットワークが木状になっている場合について考え,動的最適化による厳密解法を示す.

木ネットワークモデルは,閉路を含まない有向グラフ $G=(N,A)$ 上で義される. ここでは,$G$ を無向グラフにしたものが木であると仮定する. 点集合 $N$ は在庫地点を表し,枝 $(i,j) \in A$ が存在するとき,在庫地点 $i$ は在庫地点 $j$ に補充 を行うことを表す.複数の在庫地点から補充を受ける点においては,補充された各々の品目を用いて 別の品目を生産すると考える.

点 $i$ における品目の生産時間は,$T_i$ 日である. $T_i$ には各段階における生産時間の他に,待ち時間および輸送時間も含めて考える.

このとき点 $j$ は,複数の在庫地点から補充を受けるので,点 $j$ が品目を発注してから,すべての品目が揃うまで生産を開始できない.

点上で生産される品目は,点によって唯一に定まるものと仮定する. このとき,点と品目を表す添え字は同一であると考えられるので, 有向グラフ $G=(N,A)$ は部品展開表を意味することになる. $(i,j) \in A$ のとき,点 $j$ 上の品目 $j$ は,点 $i$ から補充される品目 $i$ をもとに生産される. 品目 $j$ を生産するのに必要な品目 $i$ の量を $\phi_{ij}$ と記す.

需要は,後続する点をもたない点 $j$ 上で発生するものとし, その $1$ 日あたりの需要量は,期待値 $\mu_j$ の定常な分布をもつものとする. 点 $j$ における $t$ 日間における需要の最大値を $D_j(t)$ 記す.

直接の需要をもたない点(後続点をもつ点) $i$ に対する需要量の期待値 $\mu_i$ は, $$ \mu_i = \sum_{(i,j) \in A} \phi_{ij} \mu_j $$ と計算できる.

点 $i$ に対する$t$ 日間における需要の最大値 $D_i(t)$ は, $$ D_i(t)= t \mu_i +\left( \sum_{(i,j) \in A} \phi_{ij} (D_j(t)-t \mu_j)^p \right)^{1/p} $$ と計算されるものと仮定する.

ここで,$p (\geq 1)$ は定数である. $p$ が大きいほど,需要の逆相関が強いことを表す. 点 $j$ 上での需要が常に同期していると仮定した場合には $p=1$ が, 点 $j$ 上での需要が独立な正規分布にしたがうと仮定した場合には $p=2$ が推奨される.

点 $i$ の保証リード時間を $L_i$ と記す. ここでは,各点 $i$ に対して,保証リード時間の下限は $0$ とし,上限は $\bar{L}_i$ であるとする.

点 $j$ が品目を発注してから,すべての品目が 揃うまでの時間(日数)を入庫リード時間とよび, $LI_j$ と記す.

点 $j$ における入庫リード時間 $LI_j$ は,以下の式を満たす. $$ L_i \leq LI_j \ \ \ \forall (i,j) \in A $$

入庫リード時間 $LI_i$ に生産時間 $T_i$ を加えたものが,補充の指示を行ってから生産を完了するまでの時間となる. これを,補充リード時間とよぶ. 補充リード時間から、保証リード時間 $L_i$ を減じた時間内の最大需要に相当する在庫を保持していれば,在庫切れの心配がないことになる. 補充リード時間から $L_i$ を減じた時間($L_{i+1}+T_i-L_i$)を正味補充時間とよぶ.

点 $i$ における安全在庫量 $I_i$ は,正味補充時間内における最大需要量から平均需要量を減じた量であるので, $$ I_i= D(LI_{i}+T_i-L_i) - (LI_{i}+T_i-L_i) \mu_i $$ となる.

記述を簡単にするために,点 $i$ において,保証リード時間が $L$, 入庫リード時間が $LI$ のときの安全在庫費用を表す,以下の関数 $HC_i(L,LI)$ を導入しておく.

$$ HC_i (L,LI)=h_i \left\{ D(LI+T_i-L) - (LI+T_i-L) \mu_i \right\} $$

上の記号を用いると,木ネットワークモデルにおける安全在庫配置問題は,以下のように定式化できる.

$$ \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)$ を無向グラフにしたものを $\bar{G}=(N,E)$ と記す. 動的計画法は,グラフ $\bar{G}$ の葉側から順次再帰方程式を計算していく.

以下では,点集合 $N$ は上のアルゴリズムによって得られた順序で $1,2,\cdots,|N|$ と番号がついているものと仮定する. 各点 $k \in N$ に対して,点集合 $\{1,2,\cdots,k\}$ の中で,点 $k$ から $\bar{G}$ 上で到達可能な点の集合 を $N_k$ と記す. 点集合 $N_k$ が下流(需要地点側)の点に対して保証リード時間 $L$ で補充するときの 最小費用を $f_k(L)$, 点集合 $N_k$ が上流(供給地点側)の点から入庫リード時間 $LI$ で補充を受けるときの 最小費用を $g_k(LI)$ と書く.

ここで,$f_k(L)$ は,保証リード時間 $L$ に対して非増加関数になっているとは限らない. これを非増加関数にするために,保証リード時間 $L$ 「以下」で補充するときの最小費用 $f^-_k(L)$ を導入する. 関数 $f^-_k(L)$ は,以下のように計算される.

$$ f^-_k(L)= \min_{T_k -L \leq x \leq L} f_k(x) $$

ここで,下限 $T_k -L$ は,正味補充時間 $x+T_k-L$ が非負であることから導かれたものである.

同様に,関数 $g_k(LI)$ を,入庫リード時間 $LI$ に対して非減少関数としたものを $g^-_k(LI)$ とする. 関数 $g^-_k(LI)$ は,以下のように計算される.

$$ g^-_k(LI)= \min_{LI \leq y \leq LI+T_k} g_k(y) $$

ここで,上限 $LI+T_k$ は,正味補充時間 $LI+T_k-y$ が非負であることから導かれたものである.

点 $k$ における保証リード時間が $L$,入庫リード時間が $LI$ のときの 点集合 $N_k$ 全体での最小費用を $c_k(L,LI)$ と書く. このとき,$c,f,g,f^-,g^-$ 間の再帰方程式は,以下のようになる.

$$ c_k(L,LI) = HC_k(L,LI) +\sum_{(i,k) \in A, i<k} f^-_i(LI) +\sum_{(k,j) \in A, j<k} g^-_j(L) $$$$ f_k(L) =\min_{LI} c_k(L,LI) $$$$ g_k(LI) =\min_{L} c_k(L,LI) $$$$ f^-_k(L)= \min_{T_k -L \leq x \leq L} f_k(x) $$$$ g^-_k(LI)= \min_{LI \leq y \leq LI+T_k} g_k(y) $$

上の動的最適化アルゴリズムの計算量は,補充リード時間の最大値を $L_{\max}$ としたとき,$O(|N| L_{\max}^2)$ となる.

このモデルの拡張として,枝 $(i,j)$ 上での輸送時間 $\tau_{ij}$ を考慮することが考えられる. これは,上の再帰方程式を $$ c_k(L,LI) = HC_k (L,LI) +\sum_{(i,k) \in A, i<k} f^-_i(LI-\tau_{ik}) +\sum_{(k,j) \in A, j<k} g^-_j(L+\tau_{kj}) $$ と変更するだけで良い.もしくは、ダミーの点を点間に挿入することによっても処理できる。

上のアルゴリズムは,入力サイズを $|N|$ としたときには,$L_{\max}=\sum_{i \in N} T_i$ であるので,擬多項式時間である.

最大需要の計算と安全在庫費用の計算

引数:

  • G: 有向グラフ
  • ProcTime: 処理時間
  • LTLB: 保証リード時間の下限
  • LTUB: 保証リード時間の上限
  • z: 安全在庫係数
  • mu: 需要の平均
  • sigma: 需要の標準偏差
  • h: 在庫費用

返値:

  • Lmax: 正味補充時間(在庫日数)の上限
  • MaxDemand: t (=0..Lmax) 日間の最大需要量
  • SafetyCost: t (=0..Lmax) 日間の安全在庫費用

max_demand_compute[source]

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 の適用例

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: 最適正味補充時間(在庫日数)

dynamic_programming_for_SSA[source]

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 関数の適用例

  • 3点の例題
  • 10点の直列多段階モデル
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(<function dynamic_programming_for_SSA.<locals>.<lambda> at 0x7ffe205fb040>, {0: 0, 1: 1, 2: 2})
NRT defaultdict(<function dynamic_programming_for_SSA.<locals>.<lambda> at 0x7ffe205fb1f0>, {0: 5, 1: 4, 2: 3})
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(<function dynamic_programming_for_SSA.<locals>.<lambda> at 0x7ffe205fb430>, {9: 0, 8: 0, 7: 0, 6: 0, 5: 0, 4: 0, 3: 0, 2: 0, 1: 0, 0: 0})
NRT defaultdict(<function dynamic_programming_for_SSA.<locals>.<lambda> at 0x7ffe205fb280>, {9: 5, 8: 5, 7: 5, 6: 5, 5: 5, 4: 5, 3: 5, 2: 5, 1: 5, 0: 5})

安全在庫配置モデルに対するメタ解法

より実践的な方法として、一般ネットワークに対応したメタ解法が考えられる。以下では、

@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}
}

に基づいたタブー探索によるメタ解法を示す。NumPuyで近傍を高速に計算している点が特徴である。

準備として、安全在庫配置問題の定式化を,正味補充時間を $NTR_i$ (net replenishment time)を変数とした定式化に変換することを考える。

$$ L_i = \max_{k: (k,i) \in A} \{ L_k \} + T_i - NRT_i $$

を元の(保証リード時間と入庫リード時間を変数とした)定式化に代入する。

全ての供給地点(入次数が $0$ の点)から点 $i$ に至るパスの集合を $PathTo_i$、全ての供給地点から需要地点(出次数が $0$ の点)に至る パスの集合を $PATH$ とする。パス $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 \in PathTo_i} \sum_{i \in N_p \setminus \{ i \}} (T_i - NRT_i ) $$ であるので、それに生産時間 $T_i$ を加えたものが、最大補充リード時間 $MaxLI_i$ になる。

また、点 $i$ から需要地点に至るパスの集合を $PathFrom_i$ と記すと、点 $i$ の最小保証リード時間 $MinLT_i$ は、 $$ \min_{p \in PathFrom_i} \sum_{i \in N_p \setminus \{ i \}} (NRT_i - T_i ) $$ と書くことができる。

点 $i$ の正味補充時間 $NTR_i$ は、補充リード時間から保証リード時間を減じたもの(と$0$の大きい方)であるので、 上流の地点からの最大補充リード時間 $MaxLI_i$ と下流の地点への最小保証リード時間 $MinLT_i$ を用いて、 以下のように計算できる。

$$ NRT_i = \left[ MaxLI_i - MinLT_i \right]^+ $$

点ごとに、正の在庫を配置するか否かを表すベクトル $b$ を与えたとき、最大補充リード時間 $MaxLI_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$ になることを利用している。

また、最小保証リード時間 $MinLT_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. $$

ここで, $LTUB_i$ は地点 $i$ に対する保証リード時間の上限である.

同時に、正味補充時間 $NRT_i$ も以下のように計算する。

$$ NRT_i = \left[ MaxLI_i - MinLT_i \right]^+ $$

全ての地点の正味補充時間が計算されたら、総費用は、安全在庫係数 $z$、需要の標準偏差 $\sigma_i$、在庫費用 $h_i$ を用いて、以下のように評価できる。

$$ \sum h_i z \sigma_i \sqrt{ NRT_i } $$

近傍や探索方法については、種々のメタ解法が考えられるが、以下ではNumPyを用いて一度に評価値を計算するタブー探索を構築する。 また、複数製品に対して同時に評価値を計算することも可能である。

決定変数 $b$ は、需要地点に対しては除外することができる。 $b$ が使われるのは、 $MaxLI$ の計算の部分だけであり、 需要地点には後続する点がないため、需要地点に対する $b$ は解の評価に影響を与えないためである。

上の方法では、供給地点、需要地点のリード時間は $0$ と仮定しているが、ダミーの点を追加したり、生産時間を変更することによって、保証リード時間の上下限付きの問題も解ける。

用語と記号:

  • 保証リード時間 (guaranteer lead time, service time): $L$
  • 生産時間 (production time): $T$
  • 入庫リード時間 (incoming lead time): $LI$
  • 補充リード時間 (replenishment lead time)
  • 正味補充時間 (net replenishment time): $NRT$
  • 最大補充リード時間: $MaxLI$
  • 最小保証リード時間: $MinLT$

安全在庫配置モデルに対するタブー探索アルゴリズム

引数:

  • 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: 最良解における最小保証リード時間

tabu_search_for_SSA[source]

tabu_search_for_SSA(G, ProcTime, LTUB, z, mu, sigma, h, max_iter=100, TLLB=1, TLUB=10, seed=1)

一般のネットワークの安全在庫配置モデルに対するタブー探索

tabu_search_for_SSA関数の使用例

  • 3点の例題
  • 10点の直列多段階モデル
  • ベンチマーク問題例
#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)
iter cost TLLB TLUB LTM
0 119.504 1 3 0.0
1 120.142 1 3 0.0
Tabu List Clear!
3 119.504 1 2 0.0
Tabu List Clear!
5 120.142 1 2 0.0
Tabu List Clear!
7 119.504 1 2 0.0
Tabu List Clear!
9 120.142 1 2 0.0
最良値 119.50357133746822
最良解 [0 1 1]
正味補充時間 [0. 9. 8.]
最大補充リード時間 [5. 9. 8.]
最小保証リード時間 [5. 0. 0.]
#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)
iter cost TLLB TLUB LTM
0 2250.602 1 3 0.0
1 2055.846 1 3 0.0
2 1947.783 1 3 0.0
3 1969.521 1 3 0.0
4 1947.783 1 3 0.0
5 1980.377 1 3 0.0
6 1905.447 1 3 0.0
7 1947.294 1 3 0.0
8 2042.342 1 3 0.0
9 2002.825 1 3 0.0
最良値 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.]

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$ 座標(整数の文字列)

read_willems[source]

read_willems(file_name='01')

Willemsのベンチマーク問題例の読み込み関数

read_wellems関数の使用例

G, pos = read_willems(file_name="01")
nx.draw(G, pos=pos, with_labels=True, node_color="Yellow")

点の属性は、以下のように点に付随する属性の辞書、もしくはget_node_attributesメソッドを用いてアクセスできる。

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配列)

extract_data_for_SSA[source]

extract_data_for_SSA(G)

ベンチマークデータからタブー探索で必要なデータを抽出する関数

extract_data_for_SSA関数とタブー探索の実行例

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)
iter cost TLLB TLUB LTM
0 23335.673 1 2 0.0
1 19827.322 1 2 0.0
2 20411.781 1 2 0.0
3 19827.322 1 2 0.0
4 20483.295 1 2 0.0
5 19827.322 1 2 0.0
6 20411.781 1 2 0.0
7 20407.165 1 2 0.0
8 20411.781 1 2 0.0
9 19827.322 1 2 0.0
10 20483.295 1 2 0.0
11 19827.322 1 2 0.0
12 20411.781 1 2 0.0
13 19827.322 1 2 0.0
14 20483.295 1 2 0.0
15 19827.322 1 2 0.0
最良値 19827.322283475798
最良解 [1 1 1 1 1 0 0 0]
正味補充時間 [10. 10. 28. 15. 10.  0.  0.  0.]
最大補充リード時間 [10. 10. 28. 15. 10.  0.  0.  0.]
最小保証リード時間 [0. 0. 0. 0. 0. 0. 0. 0.]
#在庫の有無を色で区別し、正味補充時間(在庫量)をノードの大きさで表して描画
nx.draw(G, pos=pos, with_labels=True, node_size=best_NRT*50,node_color=best_sol)
print("best obj.=", best_cost)
print("best sol.=", best_sol)
print("best NRT=", best_NRT)
best obj.= 19827.322283475798
best sol.= [1 1 1 1 1 0 0 0]
best NRT= [10. 10. 28. 15. 10.  0.  0.  0.]

安全在庫配置問題の結果を描画する関数 draw_graph_for_SSA

点の大きさが在庫量(正味補充時間)を表し、点の色(色の濃いものほど小さい)が保証リード時間を表す。

引数:

  • G: 閉路を含まない有向グラフ
  • pos: 点の座標を保管した辞書
  • best_NRT: 最良解における正味補充時間(在庫日数)
  • best_MaxLI: 最良解における最大補充リード時間
  • best_MinLT: 最良解における最小保証リード時間

返値:

  • Plotlyの図オブジェクト

draw_graph_for_SSA[source]

draw_graph_for_SSA(G, pos, best_NRT, best_MaxLI, best_MinLT)

安全在庫配置問題の結果をPlotlyで描画する関数

draw_graph_for_SSAの使用例

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: 在庫を複数の下流の点に配分するときの割合(在庫シミュレーション最適化で用いる.)

make_df_for_SSA[source]

make_df_for_SSA(G, pos, ProcTime, LTUB, z, mu, sigma, h, best_NRT, best_MaxLI, best_MinLT)

点(stage)と枝(bom)を保管したデータフレームを返す関数

make_df_for_SSAの使用例

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
#bom_df.to_csv(folder_bom + "ssa_bom02.csv")
bom_df.head()
child parent units allocation
0 Manuf_0001 Trans_0001 1 1.0
1 Manuf_0002 Trans_0004 1 1.0
2 Manuf_0003 Dist_0002 1 1.0
3 Manuf_0004 Dist_0003 1 0.5
4 Manuf_0004 Dist_0004 1 0.5

データフレームからデータの図を描画する関数

引数:

  • 点(stage)の情報を保管したデータフレーム
  • 枝(部品展開表;bom)の情報を保管したデータフレーム

返値:

  • fig : Plotlyの図オブジェクト
  • G : グラフ
  • pos : 点の座標を表す辞書

draw_graph_for_SSA_from_df[source]

draw_graph_for_SSA_from_df(stage_df, bom_df)

データフレームを入れると、安全在庫配置問題のデータをPlotlyで描画する関数

draw_graph_for_SSA_from_df関数の使用例

fig, G, pos = draw_graph_for_SSA_from_df(stage_df, bom_df)
#fig.show()
nx.draw(G, pos=pos)

全てのベンチマークデータをcsvファイルに変換する。

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")

データフレームを入力すると、一般のネットワークの安全在庫配置モデルに対するタブー探索を行う関数

引数:

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

返値:

  • stage_df : (最適化された後の)点の情報を保管したデータフレーム
  • bom_df : (最適化された後の)枝の情報を保管したデータフレーム
  • fig : Plotlyの図オブジェクト

solve_SSA[source]

solve_SSA(stage_df, bom_df)

データフレームを入力とした一般のネットワークの安全在庫配置モデルに対するタブー探索

solve_SSA関数の使用例

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
iter cost TLLB TLUB LTM
0 16575084.763 2 4 0.0
1 15759981.899 2 4 0.0
2 15045932.1 2 4 0.0
3 14835552.33 2 4 0.0
4 14835552.33 2 4 3015.7752
5 14835552.33 2 4 6031.5504
6 14635043.251 2 4 6031.5504
7 14635043.251 2 4 9047.3255
8 14635043.251 2 4 12063.1007
9 14635043.251 2 4 15078.8759
10 14635043.251 2 4 18094.6511
11 14635043.251 2 4 21110.4262
12 14635043.251 2 4 24126.2014
13 14635043.251 2 4 27141.9766
14 14835552.33 2 4 27141.9766
15 14835552.33 2 4 30157.7518
16 14835552.33 2 4 33173.5269
17 14835552.33 2 4 36189.3021
18 14835552.33 2 4 39205.0773
19 14835552.33 2 4 42220.8525
20 14835552.33 2 4 45236.6276
21 14835552.33 2 4 48252.4028
22 14835552.33 2 4 51268.178
23 14635043.251 2 4 51268.178
24 14635043.251 2 4 54283.9532
25 14635043.251 2 4 57299.7283
26 14635043.251 2 4 60315.5035
27 14635043.251 2 4 63331.2787
28 14635043.251 2 4 66347.0539
29 14635043.251 2 4 69362.829
30 14635043.251 2 4 72378.6042
31 14845423.021 2 4 72378.6042
32 14845423.021 2 4 75394.3794
33 14845423.021 2 4 78410.1546
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
#fig.show()
Image("../figure/solve_SSA.png")

データフレームをもとに定期発注方策最適化を行う関数 periodic_inv_opt

各点のリード時間は,安全在庫配置モデルの補充リード時間 (replenishment leadtime: best_MaxLI) に (定期発注方策なので) $1$ を加えたもの相当する.

ただし出荷を待っている(すでに下流からの注文がある)在庫も含めた在庫レベルを考えるものとする.

安全在庫配置の結果がすでに入ったデータフレームをもとに,定期発注方策の最適化を行う.

需要が負になることが頻繁に発生する場合には,基在庫レベルを多少大きくする必要があり,どの程度大きくすれば良いかは実験を行う必要がある.

periodic_inv_opt[source]

periodic_inv_opt(stage_df, bom_df, max_iter=1, n_samples=10, n_periods=100, seed=1