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

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

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

はじめに

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

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

基礎理論

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

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

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

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

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

引数:

  • data: データの配列

返値:

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

best_distribution[source]

best_distribution(data)

best_distribution関数の適用例

data = np.random.normal(100,10,1000)
fig, best_dist, best_fit_name, best_fit_params  = best_distribution(data)
plotly.offline.plot(fig);
print(best_fit_name, best_fit_params)
chi (94.21006051510034, -38.185731338435374, 14.332790810006365)

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

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

引数:

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

返値:

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

best_histogram[source]

best_histogram(data, nbins=50)

best_histogram関数の適用例

data = np.random.normal(100,10,1000)
fig, hist_dist = best_histogram(data, nbins = 20)
plotly.offline.plot(fig);

得られた分布に基いてランダムな需要を発生させる.

best_dist.mean(),hist_dist.mean()
(100.6427245224179, 100.00992105985364)

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

print( best_dist.ppf(0.9), hist_dist.ppf(0.9) ) #90%品切れしないための在庫量
print( best_dist.sf(200), hist_dist.sf(200) ) #在庫が200のときの品切れ率
113.2742615567341 112.910431775983
2.968721686659403e-16 0.0

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

引数:

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

返値:

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

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

fit_demand[source]

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

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

fit_demand関数の適用例

demand_df = pd.read_csv(folder+"demand_.csv") 
#demand_df = pd.read_csv(folder+"demand_all.csv")
cust_selected = [ demand_df["cust"].iloc[11]]
prod_selected = [ demand_df["prod"].iloc[1] ]
#cust_selected = ["前橋市"]
#prod_selected = ["F"]
#print(cust_selected, prod_selected)
dic  = fit_demand(demand_df, cust_selected, prod_selected, agg_period ="1d", nbins=50, L=3)
if dic[ cust_selected[0], prod_selected[0] ] is not None:
    fig, fig2, best_dist, hist_dist, best_fit_name, best_fit_params  = dic[ cust_selected[0], prod_selected[0] ]
print(best_fit_name, best_fit_params)
plotly.offline.plot(fig);
johnsonsb (0.08752303969230585, 1.4375734324106038, 1835.074189743212, 2195.4575332581912)

(Q,R)方策のシミュレーションを行う関数 simulate_inventory

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

引数:

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

返値:

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

simulate_inventory[source]

simulate_inventory(dist, 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(hist_dist, n_samples =10, n_periods =100, z=2.33) 
print("cost=",cost)
pd.DataFrame(I[0]).plot(); #在庫の推移の可視化
S= 340.35678381635483
cost= 249.03296616909125

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

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

def base_stock_simulation_using_dist(n_samples, n_periods, demand, capacity, LT, b, h, S):
    """
    単一段階在庫システムの定期発注方策に対するシミュレーションと微分値の計算
    """
    I = np.zeros((n_samples, n_periods+1))
    T = np.zeros((n_samples, n_periods+1))
    I[:, 0] = S  # initial inventory
    production = np.zeros((n_samples, LT))
    sum_dC = 0.
    for t in range(n_periods):
        I[:, t+1] = I[:, t] - demand[:, t] + \
            production[:, (t-LT) % LT]  # 在庫量の更新
        prod = np.minimum(capacity, S+demand[:, t]-I[:, t]-T[:, t])  # 生産量の計算

        #print(t,demand[:,t],I[:,t],ITI[:,t],prod,production[:, (t-LT) % LT])

        T[:, t+1] = T[:, t] + prod - production[:, (t-LT) % LT]  # 輸送中在庫量の更新
        production[:, t % LT] = prod  # 生産量の更新

        dC = np.where(I[:, t] < 0, -b, h)
        sum_dC += dC.sum()

    total_cost = (-1*b*I[I < 0].sum() + h*I[I > 0].sum())/n_periods/n_samples
    return sum_dC/n_samples/n_periods, total_cost, I
n_samples = 10
n_periods = 1000
capacity = 10.
convergence = 1e-5
b = 100
h = 1 
LT = 1
critical_ratio = b/(b+h)
S = best_dist.ppf(critical_ratio) + 5
print(S, best_dist.ppf(0.01))
demand = best_dist.rvs((n_samples,n_periods)) + 5.
8.012476283763686 -3.778827559110546
print("t:   S      dS     Cost")
for iter_ in range(100):
    dC, cost, I = base_stock_simulation_using_dist(n_samples, n_periods, demand, capacity, LT, b, h, S)
    S = S - 0.1*dC
    print(f"{iter_}: {S:.2f} {dC:.3f} {cost:.2f}")
    if dC**2<=convergence:
        break
t:   S      dS     Cost
0: 8.18 0.535 3.88
1: 8.13 0.485 3.85
2: 8.09 0.404 3.83
3: 8.06 0.344 3.81
4: 8.03 0.242 3.80
5: 8.01 0.222 3.80
6: 7.99 0.172 3.79
7: 7.98 0.162 3.79
8: 7.96 0.152 3.79
9: 7.95 0.091 3.79
10: 7.95 0.081 3.79
11: 7.94 0.051 3.78
12: 7.94 0.051 3.78
13: 7.93 0.041 3.78
14: 7.93 0.010 3.78
15: 7.93 0.010 3.78
16: 7.93 0.010 3.78
17: 7.93 0.010 3.78
18: 7.93 0.000 3.78

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: 366.95 -2.486 73.01
1: 360.02 0.693 65.04
2: 356.09 0.393 61.18
3: 355.05 0.105 60.16
4: 355.03 0.002 60.10

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

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

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

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

一般のネットワークに対して実装したものが、以下のnetwork_base_stock_simulationである。

引数:

  • G: 有向グラフ(SCMGraphのインスタンス)
  • n_samples: シミュレーションを行うサンプル数
  • n_periods: 期間数
  • demand: 需要地点ごとの需要量を保持した辞書。キーは顧客の番号。値は、需要のサンプルを表すNumPyの配列で、形状は (n_samples, n_periods)。
  • capacity: 容量(生産量上限)を表すNumPy配列
  • LT: 地点毎のリード時間を表すNumPy配列
  • ELT: 地点毎のエシェロンリード時間を表すNumPy配列
  • b: バックオーダー費用
  • h: 在庫費用
  • S: 基在庫レベル
  • phi: 部品展開表における親製品を製造するために必要な子製品の量
  • alpha: 上流の在庫地点が下流の在庫地点に在庫を配分する率

返値:

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

network_base_stock_simulation[source]

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

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

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

定期発注方策を仮定したときの、エシェロンリード時間を計算し、それをもとに基在庫レベルを計算して返す関数。 エシェロンリード時間は、最終需要地点以外の点に対しては、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.95, 17.08926 [1847.62530563  217.49740973  238.91707905]
1: 3006.19, 5.82522 [1846.8067274   217.51097292  241.18753291]
2: 3000.89, 2.75480 [1846.02908948  217.50571214  242.65383682]
3: 2998.09, 1.62961 [1845.29737183  217.49139168  243.69977714]
4: 2996.28, 1.06417 [1844.60260002  217.47909228  244.46221195]
5: 2995.01, 0.72428 [1843.94828582  217.46319584  245.00617509]
6: 2994.07, 0.54494 [1843.32732419  217.45502006  245.40527113]
7: 2993.30, 0.43767 [1842.730011    217.45163114  245.68966556]
8: 2992.67, 0.35866 [1842.1640293   217.45166499  245.88543662]
9: 2992.12, 0.31603 [1841.61648417  217.44484704  246.01264073]
10: 2991.63, 0.28328 [1841.09019154  217.44626505  246.09199254]
11: 2991.18, 0.25641 [1840.5857877   217.44053849  246.13615565]
12: 2990.78, 0.23289 [1840.1034789   217.43167377  246.15003918]
13: 2990.40, 0.21733 [1839.63730977  217.43089205  246.1457369 ]
14: 2990.04, 0.19542 [1839.19629673  217.42285795  246.11627852]
15: 2989.72, 0.18277 [1838.77053194  217.42063879  246.07771032]
16: 2989.41, 0.16754 [1838.36472953  217.41741412  246.02429507]
17: 2989.13, 0.15640 [1837.97370051  217.41310713  245.96531407]
18: 2988.86, 0.14614 [1837.59639766  217.4098182   245.90386064]
19: 2988.61, 0.13544 [1837.23391051  217.40119201  245.84089142]
20: 2988.37, 0.12520 [1836.8848687   217.40228999  245.78289211]
21: 2988.14, 0.11705 [1836.54740286  217.40023727  245.72661672]
22: 2987.93, 0.11064 [1836.21885629  217.39874373  245.67472198]
23: 2987.71, 0.10247 [1835.90308804  217.39500345  245.62232573]
24: 2987.52, 0.09693 [1835.59618264  217.40030404  245.57027426]
25: 2987.32, 0.09327 [1835.29460321  217.40133181  245.52210444]
26: 2987.13, 0.08379 [1835.00956215  217.39951229  245.47171922]
27: 2986.96, 0.08055 [1834.72877887  217.39764763  245.43037613]
28: 2986.79, 0.07318 [1834.46321624  217.39425339  245.37892772]
29: 2986.64, 0.06854 [1834.20645251  217.39259184  245.32788804]
30: 2986.49, 0.06316 [1833.95949445  217.3895499   245.28136566]
31: 2986.35, 0.05834 [1833.72243734  217.38608544  245.23517302]
32: 2986.21, 0.05304 [1833.49542662  217.38626567  245.19632617]
33: 2986.08, 0.04822 [1833.27944085  217.3870057   245.15668505]
34: 2985.96, 0.04516 [1833.0698292   217.38347793  245.12192358]
35: 2985.85, 0.04205 [1832.86827297  217.38137318  245.08422204]
36: 2985.74, 0.03913 [1832.67413149  217.38133212  245.04626724]
37: 2985.63, 0.03586 [1832.4898897   217.3808698   245.00250501]
38: 2985.54, 0.03389 [1832.31055889  217.37865597  244.96099329]
39: 2985.45, 0.03082 [1832.14016615  217.37908883  244.91872806]
40: 2985.36, 0.02814 [1831.97646732  217.3801999   244.88204575]
41: 2985.28, 0.02513 [1831.82350194  217.37479015  244.8407684 ]
42: 2985.21, 0.02265 [1831.67971848  217.37178324  244.79644754]
43: 2985.14, 0.02048 [1831.54170044  217.37143618  244.7586381 ]
44: 2985.08, 0.02011 [1831.40387548  217.377215    244.72575092]
45: 2985.01, 0.01788 [1831.27363076  217.37646674  244.6955142 ]
46: 2984.95, 0.01682 [1831.14750214  217.37581502  244.66525365]
47: 2984.89, 0.01627 [1831.02248436  217.37510563  244.64002069]
48: 2984.83, 0.01352 [1830.91097271  217.37036853  244.60741999]
49: 2984.78, 0.01248 [1830.80391694  217.37170715  244.57547878]
50: 2984.73, 0.01153 [1830.70035915  217.36954778  244.54720106]
51: 2984.69, 0.01034 [1830.60259375  217.37002131  244.51926937]
52: 2984.64, 0.00965 [1830.50919519  217.36718115  244.48901697]
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);