在庫最適化と安全在庫配置システム 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= 260.0801609745146
<matplotlib.axes._subplots.AxesSubplot at 0x7ff7686d5e20>

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

ここでは,発注を行うことができるタイミングが限定されている定期発注方策を考える. 定期発注方策では,時間を離散値として取り扱う.実際には,日単位などで在庫を管理することが多いため, 定期発注方策は,多くの実際問題で適用可能なモデルである. 時間は,期とよばれる単位に分割されており,期を $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: 367.78 -2.569 74.84
1: 360.93 0.685 65.85
2: 356.88 0.405 62.05
3: 355.65 0.124 60.94
4: 355.51 0.014 60.85
5: 355.53 -0.002 60.85

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

ここでは,単一段階モデルを直列多段階の場合に拡張する. この結果は,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_{0=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を下回ったら終了する。

# 多段階 基在庫レベル方策 在庫シミュレーション
import numpy as np
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"Inventorys for all stages", showlegend=False)

#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日を加えて計算する。

def initial_base_stock_evel(G, LT):
    ''' 
    初期基在庫レベルとエシェロンリード時間の計算
    '''
    for i in G.up_order():
        if G.out_degree(i) == 0:
            ELT[i] = LT[i]
        else:
            max_succ_LT = 0
            for j in G.successors(i):
                max_succ_LT = max(ELT[j], max_succ_LT)
            ELT[i] = max_succ_LT+LT[i]+1  # add one day cycle time
    S = ELT*mu + z*sigma * np.sqrt(ELT)
    return ELT, S

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 = 10
n_periods = 100

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_evel(G, LT)
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: 3064.50, 25.24981 [1848.01014691  217.92475     239.79553567]
1: 3044.39, 5.14276 [1847.64264691  218.17725     242.01903567]
2: 3039.40, 1.87026 [1847.27264691  218.32725     243.32703567]
3: 3037.15, 0.53619 [1847.09377191  218.286125    244.03591067]
4: 3036.21, 0.41027 [1846.90939691  218.2505      244.64828567]
5: 3035.45, 0.20768 [1846.72002191  218.219875    245.06166067]
6: 3034.96, 0.14048 [1846.52364691  218.19625     245.38003567]
7: 3034.57, 0.04757 [1846.44539691  218.1595      245.58028567]
8: 3034.39, 0.04757 [1846.36714691  218.12275     245.78053567]
9: 3034.21, 0.04757 [1846.28889691  218.086       245.98078567]
10: 3034.04, 0.02947 [1846.31264691  217.94725     246.07903567]
11: 3034.08, 0.00931 [1846.31314691  218.04175     246.09853567]

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

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

過去の需要系列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: 処理時間
  • LTLB: 保証リード時間の下限
  • LTUB: 保証リード時間の上限
  • z: 安全在庫係数
  • mu: 需要の平均
  • sigma: 需要の標準偏差
  • h: 在庫費用

返値:

  • 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 0x7ff6e8ed71f0>, {0: 0, 1: 1, 2: 2})
NRT defaultdict(<function dynamic_programming_for_SSA.<locals>.<lambda> at 0x7ff6e8ed73a0>, {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 0x7ff6e8ed7550>, {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 0x7ff6e8ed7160>, {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} 0 & \forall i \in 需要地点 \\ \min_{ (i,j) \in A } NRT_j + MinLT_j - T_j & それ以外 \end{array} \right. $$

同時に、正味補充時間 $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: 処理時間
  • LTUB: 保証リード時間上限
  • z: 安全在庫係数
  • mu: 需要の平均
  • sigma: 需要の標準偏差
  • h: 在庫費用
  • 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., 5., 2.])
mu = np.array([300., 200., 100.])
sigma = np.array([12., 10., 15.])
LTUB = np.zeros(3)

ProcTime = np.array([5, 4 ,3], int) # 動的最適化の例題と合わせるため、生産時間から保証リード時間上限を減じてある

tabu_search_for_SSA(G, ProcTime,  LTUB, z, mu, sigma, h, max_iter = 10, TLLB =1, TLUB =3, seed = 1)
iter cost TLLB TLUB LTM
0 387.507 1 3 0.0
1 295.011 1 3 0.0
Tabu List Clear!
3 387.507 1 2 0.0
Tabu List Clear!
5 295.011 1 2 0.0
Tabu List Clear!
7 387.507 1 2 0.0
Tabu List Clear!
9 295.011 1 2 0.0
(295.0106609291552,
 array([1, 1, 1]),
 array([5., 4., 3.]),
 array([5., 4., 3.]),
 array([0., 0., 0.]))
#10点の直列多段階モデル
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)
LTUB = np.zeros(n)

tabu_search_for_SSA(G, ProcTime, LTUB, z, mu, sigma, h, max_iter = 10, TLLB =1, TLUB =3, seed = 1)
iter cost TLLB TLUB LTM
0 1365.12 1 3 0.0
1 1048.427 1 3 0.0
2 820.822 1 3 0.0
3 611.832 1 3 0.0
4 423.751 1 3 0.0
5 258.266 1 3 0.0
6 116.673 1 3 0.0
7 454.365 1 3 0.0
8 531.712 1 3 0.0
9 454.365 1 3 0.0
(116.67261889578035,
 array([0, 0, 0, 0, 0, 0, 0, 0, 0, 1]),
 array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., 50.]),
 array([ 5., 10., 15., 20., 25., 30., 35., 40., 45., 50.]),
 array([ 5., 10., 15., 20., 25., 30., 35., 40., 45.,  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="02")
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") )
30 14 14 14 14 5 5 5 5 15 15 15 15 {'Retail_0001': '20', 'Retail_0002': '20', 'Retail_0003': '20', 'Retail_0004': '20'}

在庫費用、需要量、需要の標準偏差の計算

ベンチマークデータからタブー探索で必要なデータを抽出し、必要なデータを準備する。 各地点には付加された費用と生産時間、需要地点には需要と標準偏差と保証リー時間上限(サービス時間)が入っている。 すべての点に対して、在庫費用と需要・標準偏差など、安全在庫配置モデルの計算に必要なデータを計算する関数を準備しておく。

引数:

  • 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)
iter cost TLLB TLUB LTM
0 27894349.902 1 3 0.0
1 27304807.83 1 3 0.0
2 27029688.196 1 3 0.0
3 27029688.196 1 3 8464.5883
4 27029688.196 1 3 16929.1766
5 27029688.196 1 3 25393.7649
6 27029688.196 1 3 33858.3533
7 27029688.196 1 3 42322.9416
8 27029688.196 1 3 50787.5299
9 27029688.196 1 3 59252.1182
10 27029688.196 1 3 67716.7065
11 27029688.196 1 3 76181.2948
12 27029688.196 1 3 84645.8832
13 27183931.54 1 3 84645.8832
14 27183931.54 1 3 93110.4715
15 27029688.196 1 3 93110.4715
16 27029688.196 1 3 101575.0598
17 27029688.196 1 3 110039.6481
18 27029688.196 1 3 118504.2364
19 27029688.196 1 3 126968.8247
20 27029688.196 1 3 135433.413
21 27304807.83 1 3 135433.413
22 27304807.83 1 3 143898.0014
23 27304807.83 1 3 152362.5897
24 27029688.196 1 3 152362.5897
25 27183931.54 1 3 152362.5897
#在庫の有無を色で区別し、正味補充時間(在庫量)をノードの大きさで表して描画
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.= 27029688.195813652
best sol.= [1 1 0 0 1 1 1 1 1 1 1 1 1]
best NRT= [30.  0.  0.  0.  0.  0.  0.  0.  0. 14. 14. 14. 14.]

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

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

引数:

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

安全在庫配置問題の結果をデータフレームに変換する関数

点(stage)の情報を保管したデータフレームと枝(部品展開表;bom)の情報を保管したデータフレームを構築する。

引数:

  • G: 閉路を含まない有向グラフ
  • pos: 点の座標を保管した辞書
  • ProcTime: 処理時間 (NumPy配列)
  • LTUB: 保証リード時間上限 (NumPy配列)
  • z: 安全在庫係数 (NumPy配列)
  • mu: 需要の平均 (NumPy配列)
  • sigma: 需要の標準偏差 (NumPy配列)
  • h: 在庫費用 (NumPy配列)
  • best_NRT: 最良解における正味補充時間(在庫日数)
  • best_MaxLI: 最良解における最大補充リード時間
  • best_MinLT: 最良解における最小保証リード時間

返値:

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

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.head()
name net_replenishment_time max_guaranteed_LT processing_time replenishment_LT guaranteed_LT z average_demand sigma h x y
0 Manuf_0001 30.0 0.0 30.0 30.0 0.0 1.750686 22700.0 15246.638974 80.00 64 48
1 Manuf_0002 0.0 0.0 14.0 14.0 15.0 1.750686 4000.0 7500.000000 92.55 416 48
2 Manuf_0003 0.0 0.0 14.0 14.0 15.0 1.750686 2000.0 3500.000000 92.50 416 112
3 Manuf_0004 0.0 0.0 14.0 14.0 15.0 1.750686 8000.0 9000.000000 92.35 416 176
4 Manuf_0005 0.0 0.0 14.0 14.0 15.0 1.750686 8700.0 9108.238029 92.35 416 240
#bom_df.to_csv(folder_bom + "ssa_bom02.csv")
bom_df.head()
child parent
0 Manuf_0001 Trans_0001
1 Manuf_0001 Trans_0002
2 Manuf_0001 Trans_0003
3 Manuf_0001 Trans_0004
4 Manuf_0002 Retail_0001

データフレームからタブー探索を行う関数

引数:

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

tabu_search_for_SSA_from_df[source]

tabu_search_for_SSA_from_df(stage_df, bom_df, max_iter=None, TLLB=None, TLUB=None, seed=1)

データフレームからタブー探索を行う関数

tabu_search_for_SSA_from_df関数の使用例

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, best_sol, best_NRT, best_MaxLI, best_MinLT = tabu_search_for_SSA_from_df(stage_df, bom_df)
print(best_cost)
iter cost TLLB TLUB LTM
0 33226181.195 1 3 0.0
1 29630717.895 1 3 0.0
2 26349600.316 1 3 0.0
3 23670297.557 1 3 0.0
4 23670297.557 1 3 11418.2899
5 19889447.976 1 3 11418.2899
6 19889447.976 1 3 22836.5797
7 19889447.976 1 3 34254.8696
8 21874527.988 1 3 34254.8696
9 21874527.988 1 3 45673.1594
10 21874527.988 1 3 57091.4493
11 21874527.988 1 3 68509.7391
12 21874527.988 1 3 79928.029
13 21874527.988 1 3 91346.3188
14 19889447.976 1 3 91346.3188
15 19889447.976 1 3 102764.6087
16 19889447.976 1 3 114182.8985
17 19889447.976 1 3 125601.1884
18 19889447.976 1 3 137019.4782
19 19889447.976 1 3 148437.7681
20 19889447.976 1 3 159856.0579
21 19889447.976 1 3 171274.3478
22 19889447.976 1 3 182692.6376
23 21874527.988 1 3 182692.6376
24 21874527.988 1 3 194110.9275
25 21874527.988 1 3 205529.2174
19889447.976331897

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

引数:

  • 点(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)
stage_df["net_replenishment_time"] = best_NRT
stage_df["max_guaranteed_LT"] = best_MaxLI
stage_df["replenishment_LT"] = best_MinLT 
stage_df
name net_replenishment_time max_guaranteed_LT processing_time replenishment_LT guaranteed_LT z average_demand sigma h x y
0 Manuf_0001 0.0 30.0 30.0 30.0 0.0 1.750686 22700.0 15246.638974 80.00 64 48
1 Manuf_0002 0.0 44.0 14.0 44.0 0.0 1.750686 4000.0 7500.000000 92.55 416 48
2 Manuf_0003 0.0 44.0 14.0 44.0 0.0 1.750686 2000.0 3500.000000 92.50 416 112
3 Manuf_0004 0.0 44.0 14.0 44.0 0.0 1.750686 8000.0 9000.000000 92.35 416 176
4 Manuf_0005 0.0 44.0 14.0 44.0 0.0 1.750686 8700.0 9108.238029 92.35 416 240
5 Retail_0001 0.0 49.0 5.0 49.0 0.0 1.750686 4000.0 7500.000000 92.55 592 48
6 Retail_0002 34.0 54.0 5.0 20.0 0.0 1.750686 2000.0 3500.000000 92.50 592 112
7 Retail_0003 0.0 49.0 5.0 49.0 0.0 1.750686 8000.0 9000.000000 184.70 592 176
8 Retail_0004 34.0 54.0 5.0 20.0 0.0 1.750686 700.0 1400.000000 92.35 592 240
9 Trans_0001 0.0 59.0 15.0 59.0 0.0 1.750686 4000.0 7500.000000 80.55 240 48
10 Trans_0002 74.0 74.0 15.0 0.0 0.0 1.750686 2000.0 3500.000000 80.50 240 112
11 Trans_0003 0.0 59.0 15.0 59.0 0.0 1.750686 8000.0 9000.000000 80.35 240 176
12 Trans_0004 74.0 74.0 15.0 0.0 0.0 1.750686 8700.0 9108.238029 80.35 240 240

全てのベンチマークデータを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")

Dash GUI

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

引数:

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

返値:

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

solve_SSA[source]

solve_SSA(stage_df, bom_df)

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

solve_SSA関数の使用例

import pandas as pd
stage_df = pd.read_csv(folder_bom + "ssa02.csv")
bom_df = pd.read_csv(folder_bom + "ssa_bom02.csv")
stage_df, bom_df, fig = solve_SSA(stage_df, bom_df)
stage_df.head()
iter cost TLLB TLUB LTM
0 27894349.902 1 3 0.0
1 27304807.83 1 3 0.0
2 27029688.196 1 3 0.0
3 27029688.196 1 3 8464.5883
4 27029688.196 1 3 16929.1766
5 27029688.196 1 3 25393.7649
6 27029688.196 1 3 33858.3533
7 27029688.196 1 3 42322.9416
8 27029688.196 1 3 50787.5299
9 27029688.196 1 3 59252.1182
10 27029688.196 1 3 67716.7065
11 27029688.196 1 3 76181.2948
12 27029688.196 1 3 84645.8832
13 27183931.54 1 3 84645.8832
14 27183931.54 1 3 93110.4715
15 27029688.196 1 3 93110.4715
16 27029688.196 1 3 101575.0598
17 27029688.196 1 3 110039.6481
18 27029688.196 1 3 118504.2364
19 27029688.196 1 3 126968.8247
20 27029688.196 1 3 135433.413
21 27304807.83 1 3 135433.413
22 27304807.83 1 3 143898.0014
23 27304807.83 1 3 152362.5897
24 27029688.196 1 3 152362.5897
25 27183931.54 1 3 152362.5897
name net_replenishment_time max_guaranteed_LT processing_time replenishment_LT guaranteed_LT z average_demand sigma h x y
0 Manuf_0001 30.0 0.0 30.0 30.0 0.0 1.750686 22700.0 15246.638974 80.00 64 48
1 Manuf_0002 0.0 0.0 14.0 14.0 15.0 1.750686 4000.0 7500.000000 92.55 416 48
2 Manuf_0003 0.0 0.0 14.0 14.0 15.0 1.750686 2000.0 3500.000000 92.50 416 112
3 Manuf_0004 0.0 0.0 14.0 14.0 15.0 1.750686 8000.0 9000.000000 92.35 416 176
4 Manuf_0005 0.0 0.0 14.0 14.0 15.0 1.750686 8700.0 9108.238029 92.35 416 240
#fig.show()
Image("../figure/solve_SSA.png")