# 例題集：在庫最適化


<!-- WARNING: THIS FILE WAS AUTOGENERATED! DO NOT EDIT! -->

``` python
mean = 100.
sigma = 10.
normal = norm(mean,sigma) #需要量
price = 100.   #販売価格
purchase = 70. #仕入れ値
demand = normal.rvs(100)
lostsales = 0.
salvage = 10. #バイトに半値で売るとき 50 :-) 
loyality = 0.5
Q = 100 #発注量
QLB = 50
QUB = 130

retail_profit, seven_profit = [], [] 
for Q in range(QLB,QUB):
    
    profit_retail, profit_seven, lostsales_cost, salvage_value = 0,0,0,0
    for d in list(demand):
        purchase_cost = purchase*Q
        profit = price*min(Q,d)
        net_margin = profit-purchase_cost  #粗利益（通常会計 -> 本社利益と小売店利益は同じ）
        profit_retail +=  (net_margin + salvage*max(Q-d,0) -lostsales*max(d-Q,0))*(1.-loyality) #loyality = phi (=supplierへの利益分配率）)
        profit_seven  += (net_margin+ salvage*max(Q-d,0)- lostsales*max(d-Q,0) )*loyality
        lostsales_cost += lostsales*max(d-Q,0) #機会損出費用
        salvage_value += salvage*max(Q-d,0)    #売れ残り価値

    retail_profit.append( profit_retail)
    seven_profit.append( profit_seven)
    
#print(profit_retail, profit_seven, lostsales_cost, salvage_cost)

#最適な発注量を求める => プロット
# min_cost = np.max
retail_profit2, seven_profit2 = [], [] 
for Q in range(QLB,QUB):
    
    profit_retail, profit_seven, lostsales_cost, salvage_value = 0,0,0,0
    for d in list(demand):
        purchase_cost = purchase*min(Q,d) #売れた分だけ仕入れ値を払う（コンビニ会計）
        profit = price*min(Q,d)
        net_margin = profit-purchase_cost #粗利益（コンビニ会計だとnet_marginが大きくなる）
        profit_retail += net_margin*(1.-loyality) -purchase*max(Q-d,0)+salvage*max(Q-d,0) -lostsales*max(d-Q,0)*(1.-loyality) #売れ残りに対する原価を店舗利益から引く（コンビニ会計）
        profit_seven  += net_margin*loyality- lostsales*max(d-Q,0)*loyality
        lostsales_cost += lostsales*max(d-Q,0) #機会損出費用
        salvage_value += salvage*max(Q-d,0)    #売れ残り価値
    
    retail_profit2.append( profit_retail)
    seven_profit2.append( profit_seven)
    
#print(profit_retail, profit_seven, lostsales_cost, salvage_cost)
```

``` python
cost_df = pd.DataFrame({"発注量":list(range(QLB,QUB)),"小売利益":retail_profit, "本社利益":seven_profit,
                       "小売利益（コンビニ会計）":retail_profit2, "本社利益（コンビニ会計）":seven_profit2})
cost_df.head()
```

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }
&#10;    .dataframe tbody tr th {
        vertical-align: top;
    }
&#10;    .dataframe thead th {
        text-align: right;
    }
</style>

<table class="dataframe" data-quarto-postprocess="true" data-border="1">
<thead>
<tr style="text-align: right;">
<th data-quarto-table-cell-role="th"></th>
<th data-quarto-table-cell-role="th">発注量</th>
<th data-quarto-table-cell-role="th">小売利益</th>
<th data-quarto-table-cell-role="th">本社利益</th>
<th data-quarto-table-cell-role="th">小売利益（コンビニ会計）</th>
<th data-quarto-table-cell-role="th">本社利益（コンビニ会計）</th>
</tr>
</thead>
<tbody>
<tr>
<td data-quarto-table-cell-role="th">0</td>
<td>50</td>
<td>50000.0</td>
<td>50000.0</td>
<td>50000.0</td>
<td>50000.0</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">1</td>
<td>51</td>
<td>51000.0</td>
<td>51000.0</td>
<td>51000.0</td>
<td>51000.0</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">2</td>
<td>52</td>
<td>52000.0</td>
<td>52000.0</td>
<td>52000.0</td>
<td>52000.0</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">3</td>
<td>53</td>
<td>53000.0</td>
<td>53000.0</td>
<td>53000.0</td>
<td>53000.0</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">4</td>
<td>54</td>
<td>54000.0</td>
<td>54000.0</td>
<td>54000.0</td>
<td>54000.0</td>
</tr>
</tbody>
</table>

</div>

``` python
fig = px.scatter(cost_df, x="発注量",y=["小売利益", "本社利益","小売利益（コンビニ会計）", "本社利益（コンビニ会計）" ])
plotly.offline.plot(fig);
```

\#hide \## 生産計画モデル

最終工程の充填工程と前工程は独立に解ける．（中間の保管場所が制約にならないため．）

- 26の製品（最小と最大ロットサイズ）；製造名称が親製品に対応？

- 25の充填ライン

データに需要がない． 人数制約がない．

集合

- 製品の集合 *P*
- 期の集合 *T*
- ラインの集合 *L*

パラメータ

- 需要量 *d*<sub>*p**t*</sub>
- 段取り費用 *F*<sub>*p**l*</sub>
- 生産費用 *c*<sub>*p**l*</sub>
- 在庫費用 *h*<sub>*p*</sub>
- 最小，最大ロットサイズ
  *M**i**n**L*<sub>*p*</sub>, *M**a**x**L*<sub>*p*</sub>
- 最大人数 *M**a**x**W*<sub>*t*</sub>
- 最大生産時間 *U**B*<sub>*l**t*</sub>
- 必要人数 *a*<sub>*p**l*</sub>
- 生産時間 *w*<sub>*p**l*</sub>
- 段取り時間 *S**U**T*<sub>*p**l*</sub>

変数 - *t* 期の製造量 *x*<sub>*p**l**t*</sub> （整数） - *t*
期に段取りするとき 1 *y*<sub>*p**l**t*</sub> - *t* 期の在庫量
*I*<sub>*p**t*</sub>

*m**i**n**i**m**i**z**e*∑*h*<sub>*p*</sub>*I*<sub>*p**t*</sub> + ∑*F*<sub>*p**l*</sub>*y*<sub>*p**l**t*</sub> + ∑*c*<sub>*p**l*</sub>*x*<sub>*p**l**t*</sub>

*I*<sub>*p*, *t* − 1</sub> + ∑<sub>*l*</sub>*x*<sub>*p**l**t*</sub> = *d*<sub>*p**t*</sub> + *I*<sub>*p**t*</sub>   ∀*p*, *t*

*M**i**n**L*<sub>*p*</sub>*y*<sub>*p**l**t*</sub> ≤ *x*<sub>*p**l**t*</sub> ≤ *M**a**x**L*<sub>*p*</sub>*y*<sub>*p**l**t*</sub>   ∀*p*, *l*, *t*

∑<sub>*p*</sub>*w*<sub>*p**l*</sub>*x*<sub>*p**l**t*</sub> + *S**U**T*<sub>*p**l*</sub>*y*<sub>*p**l**t*</sub> ≤ *U**B*<sub>*l**t*</sub>   ∀*l*, *t*

∑<sub>*p*, *l*</sub>*a*<sub>*p**l*</sub>*x*<sub>*p**l**t*</sub> ≤ *M**a**x**W*<sub>*t*</sub>   ∀*t*

## 予測と在庫管理の融合

従来の需要予測の研究と在庫管理の研究は別々に行われてきた．
予測では，非定常な需要を仮定し，在庫管理では定常な分布を仮定する場合が多い．
これらの2つの（おそらくサプライ・チェインに対して最も重要かつ最も多くの研究が行われてきた）モデルを構築することは，永年の研究者たちの夢であったが，いまだに実務的で納得するものは出ていない．

従来の研究の多くは，以下の２つに分類される．

1.  綺麗な解析ができるような（実務的にはかなり特殊な）仮定の元での理論研究
2.  実際の問題をなんとか解決するための（理論的な背景はほとんど考慮されていない）ヒューリスティクス

以下では，理論的にもある程度きちんとした枠組みで，かつ実際問題を解決可能な方法について考える．方針は以下の通り．

1.  予測を行い誤差の分布を適当な分布で近似する．

2.  リード時間分の予測需要の合計と分布に基づく安全在庫の和が変動型の基在庫レベルになる．

3.  需要が負にならないようにスライドして，基在庫レベル最適化を行い，安全在庫量の最適化を行う．

4.  安全在庫と需要予測に基づく確定在庫を合わせることによって，動的な基在庫レベルを導く．

### Integrating Forecasting and Inventory Management

Traditional research on demand forecasting and inventory management have
been conducted separately. Forecasting often assumes a non-stationary
demand, while inventory management assumes a stationary distribution.
Building these two (probably the most important and most studied) models
for the supply chain has been a dream of researchers for a long time,
but nothing practical and convincing has been produced yet.

Most conventional research can be divided into two categories.

- Theoretical research based on assumptions that allow for clean
  analysis (which are quite specific from a practical standpoint)
- Heuristics to somehow solve real problems (with little or no
  consideration of theoretical background).

In the following, we will consider methods that are both theoretically
sound and capable of solving practical problems. The policy is as
follows.

1.  Make a forecast and approximate the distribution of the error with a
    suitable distribution.

2.  The sum of the total demand forecast for the lead time and the
    safety stock based on the distribution is the base stock level for
    the variable type.

3.  The base inventory level is optimized by sliding the demand so that
    it does not become negative, and the amount of safety stock is
    optimized.

4.  The dynamic base stock level is derived by combining the safety
    stock and the fixed stock based on the demand forecast.

### サンプル需要データの読み込み

プロモーションデータを付加した需要データを読み込む．顧客と製品を選択肢，1日単位の需要系列を生成する．
2019年の年初から2020年の年末までのデータであり，2020年10月以降を検証用として用いる．

予測手法は自動機械学習を用いる．この部分は深層学習やベイズ推論に置き換えても良い．

### Loading sample demand data

Load the demand data with promotion data added. Choose customers and
products, and generate a daily demand series. The data covers the period
from the beginning of 2019 to the end of 2020, and the data after
October 2020 is used for validation.

The forecasting method is based on automatic machine learning. This part
can be replaced by deep learning or Bayesian inference.

``` python
promo_df = pd.read_csv(folder+"promo.csv", index_col=0)
demand_df = pd.read_csv(folder+"demand_with_promo_all.csv")
c = "仙台市" 
p = "C"
print(c,p)
agg_period ="1d"
df, future = make_forecast_df(demand_df, customer=c, product=p,promo_df= promo_df, agg_period="1d", forecast_periods = 10)
horizon="2020/10/01"
best, result_df = automl(df, horizon, 5)
```

<style  type="text/css" >
    #T_05252_ th {
          text-align: left;
    }#T_05252_row0_col0,#T_05252_row0_col6,#T_05252_row1_col0,#T_05252_row1_col1,#T_05252_row1_col2,#T_05252_row1_col3,#T_05252_row1_col4,#T_05252_row1_col5,#T_05252_row1_col6,#T_05252_row2_col0,#T_05252_row2_col1,#T_05252_row2_col2,#T_05252_row2_col3,#T_05252_row2_col4,#T_05252_row2_col5,#T_05252_row3_col0,#T_05252_row3_col1,#T_05252_row3_col2,#T_05252_row3_col3,#T_05252_row3_col4,#T_05252_row3_col5,#T_05252_row3_col6,#T_05252_row4_col0,#T_05252_row4_col1,#T_05252_row4_col2,#T_05252_row4_col3,#T_05252_row4_col4,#T_05252_row4_col5,#T_05252_row4_col6,#T_05252_row5_col0,#T_05252_row5_col1,#T_05252_row5_col2,#T_05252_row5_col3,#T_05252_row5_col4,#T_05252_row5_col5,#T_05252_row5_col6,#T_05252_row6_col0,#T_05252_row6_col1,#T_05252_row6_col2,#T_05252_row6_col3,#T_05252_row6_col4,#T_05252_row6_col5,#T_05252_row6_col6,#T_05252_row7_col0,#T_05252_row7_col1,#T_05252_row7_col2,#T_05252_row7_col3,#T_05252_row7_col4,#T_05252_row7_col5,#T_05252_row7_col6,#T_05252_row8_col0,#T_05252_row8_col1,#T_05252_row8_col2,#T_05252_row8_col3,#T_05252_row8_col4,#T_05252_row8_col5,#T_05252_row8_col6,#T_05252_row9_col0,#T_05252_row9_col1,#T_05252_row9_col2,#T_05252_row9_col3,#T_05252_row9_col4,#T_05252_row9_col5,#T_05252_row9_col6,#T_05252_row10_col0,#T_05252_row10_col1,#T_05252_row10_col2,#T_05252_row10_col3,#T_05252_row10_col4,#T_05252_row10_col5,#T_05252_row10_col6,#T_05252_row11_col0,#T_05252_row11_col1,#T_05252_row11_col2,#T_05252_row11_col3,#T_05252_row11_col4,#T_05252_row11_col5,#T_05252_row11_col6,#T_05252_row12_col0,#T_05252_row12_col1,#T_05252_row12_col2,#T_05252_row12_col3,#T_05252_row12_col4,#T_05252_row12_col5,#T_05252_row12_col6,#T_05252_row13_col0,#T_05252_row13_col1,#T_05252_row13_col2,#T_05252_row13_col3,#T_05252_row13_col4,#T_05252_row13_col5,#T_05252_row13_col6,#T_05252_row14_col0,#T_05252_row14_col1,#T_05252_row14_col2,#T_05252_row14_col3,#T_05252_row14_col4,#T_05252_row14_col5,#T_05252_row14_col6,#T_05252_row15_col0,#T_05252_row15_col1,#T_05252_row15_col2,#T_05252_row15_col3,#T_05252_row15_col4,#T_05252_row15_col5,#T_05252_row15_col6,#T_05252_row16_col0,#T_05252_row16_col1,#T_05252_row16_col2,#T_05252_row16_col3,#T_05252_row16_col4,#T_05252_row16_col5,#T_05252_row16_col6,#T_05252_row17_col0,#T_05252_row17_col1,#T_05252_row17_col2,#T_05252_row17_col3,#T_05252_row17_col4,#T_05252_row17_col5,#T_05252_row17_col6,#T_05252_row18_col0,#T_05252_row18_col1,#T_05252_row18_col2,#T_05252_row18_col3,#T_05252_row18_col4,#T_05252_row18_col5,#T_05252_row18_col6{
            text-align:  left;
            text-align:  left;
        }#T_05252_row0_col1,#T_05252_row0_col2,#T_05252_row0_col3,#T_05252_row0_col4,#T_05252_row0_col5,#T_05252_row2_col6{
            text-align:  left;
            text-align:  left;
            background-color:  yellow;
        }#T_05252_row0_col7,#T_05252_row1_col7,#T_05252_row2_col7,#T_05252_row3_col7,#T_05252_row4_col7,#T_05252_row5_col7,#T_05252_row6_col7,#T_05252_row7_col7,#T_05252_row9_col7,#T_05252_row10_col7,#T_05252_row11_col7{
            text-align:  left;
            text-align:  left;
            background-color:  lightgrey;
        }#T_05252_row8_col7,#T_05252_row12_col7,#T_05252_row13_col7,#T_05252_row14_col7,#T_05252_row15_col7,#T_05252_row16_col7,#T_05252_row17_col7,#T_05252_row18_col7{
            text-align:  left;
            text-align:  left;
            background-color:  yellow;
            background-color:  lightgrey;
        }</style>

<table id="T_05252_" data-quarto-postprocess="true">
<thead>
<tr>
<th class="blank level0" data-quarto-table-cell-role="th"></th>
<th class="col_heading level0 col0"
data-quarto-table-cell-role="th">Model</th>
<th class="col_heading level0 col1"
data-quarto-table-cell-role="th">MAE</th>
<th class="col_heading level0 col2"
data-quarto-table-cell-role="th">MSE</th>
<th class="col_heading level0 col3"
data-quarto-table-cell-role="th">RMSE</th>
<th class="col_heading level0 col4"
data-quarto-table-cell-role="th">R2</th>
<th class="col_heading level0 col5"
data-quarto-table-cell-role="th">RMSLE</th>
<th class="col_heading level0 col6"
data-quarto-table-cell-role="th">MAPE</th>
<th class="col_heading level0 col7" data-quarto-table-cell-role="th">TT
(Sec)</th>
</tr>
</thead>
<tbody>
<tr>
<td id="T_05252_level0_row0" class="row_heading level0 row0"
data-quarto-table-cell-role="th">et</td>
<td id="T_05252_row0_col0" class="data row0 col0">Extra Trees
Regressor</td>
<td id="T_05252_row0_col1" class="data row0 col1">0.8184</td>
<td id="T_05252_row0_col2" class="data row0 col2">1.4329</td>
<td id="T_05252_row0_col3" class="data row0 col3">1.1970</td>
<td id="T_05252_row0_col4" class="data row0 col4">0.9259</td>
<td id="T_05252_row0_col5" class="data row0 col5">0.3454</td>
<td id="T_05252_row0_col6" class="data row0 col6">0.3228</td>
<td id="T_05252_row0_col7" class="data row0 col7">0.1400</td>
</tr>
<tr>
<td id="T_05252_level0_row1" class="row_heading level0 row1"
data-quarto-table-cell-role="th">catboost</td>
<td id="T_05252_row1_col0" class="data row1 col0">CatBoost
Regressor</td>
<td id="T_05252_row1_col1" class="data row1 col1">0.9641</td>
<td id="T_05252_row1_col2" class="data row1 col2">1.6140</td>
<td id="T_05252_row1_col3" class="data row1 col3">1.2704</td>
<td id="T_05252_row1_col4" class="data row1 col4">0.9166</td>
<td id="T_05252_row1_col5" class="data row1 col5">0.4705</td>
<td id="T_05252_row1_col6" class="data row1 col6">0.3642</td>
<td id="T_05252_row1_col7" class="data row1 col7">1.2500</td>
</tr>
<tr>
<td id="T_05252_level0_row2" class="row_heading level0 row2"
data-quarto-table-cell-role="th">rf</td>
<td id="T_05252_row2_col0" class="data row2 col0">Random Forest
Regressor</td>
<td id="T_05252_row2_col1" class="data row2 col1">0.9471</td>
<td id="T_05252_row2_col2" class="data row2 col2">1.8677</td>
<td id="T_05252_row2_col3" class="data row2 col3">1.3666</td>
<td id="T_05252_row2_col4" class="data row2 col4">0.9035</td>
<td id="T_05252_row2_col5" class="data row2 col5">0.3581</td>
<td id="T_05252_row2_col6" class="data row2 col6">0.2946</td>
<td id="T_05252_row2_col7" class="data row2 col7">0.1700</td>
</tr>
<tr>
<td id="T_05252_level0_row3" class="row_heading level0 row3"
data-quarto-table-cell-role="th">lightgbm</td>
<td id="T_05252_row3_col0" class="data row3 col0">Light Gradient
Boosting Machine</td>
<td id="T_05252_row3_col1" class="data row3 col1">1.0635</td>
<td id="T_05252_row3_col2" class="data row3 col2">2.0255</td>
<td id="T_05252_row3_col3" class="data row3 col3">1.4232</td>
<td id="T_05252_row3_col4" class="data row3 col4">0.8953</td>
<td id="T_05252_row3_col5" class="data row3 col5">0.4909</td>
<td id="T_05252_row3_col6" class="data row3 col6">0.3828</td>
<td id="T_05252_row3_col7" class="data row3 col7">0.1600</td>
</tr>
<tr>
<td id="T_05252_level0_row4" class="row_heading level0 row4"
data-quarto-table-cell-role="th">xgboost</td>
<td id="T_05252_row4_col0" class="data row4 col0">Extreme Gradient
Boosting</td>
<td id="T_05252_row4_col1" class="data row4 col1">1.0059</td>
<td id="T_05252_row4_col2" class="data row4 col2">2.0651</td>
<td id="T_05252_row4_col3" class="data row4 col3">1.4370</td>
<td id="T_05252_row4_col4" class="data row4 col4">0.8933</td>
<td id="T_05252_row4_col5" class="data row4 col5">0.4459</td>
<td id="T_05252_row4_col6" class="data row4 col6">0.3132</td>
<td id="T_05252_row4_col7" class="data row4 col7">0.3900</td>
</tr>
<tr>
<td id="T_05252_level0_row5" class="row_heading level0 row5"
data-quarto-table-cell-role="th">gbr</td>
<td id="T_05252_row5_col0" class="data row5 col0">Gradient Boosting
Regressor</td>
<td id="T_05252_row5_col1" class="data row5 col1">1.1372</td>
<td id="T_05252_row5_col2" class="data row5 col2">2.7493</td>
<td id="T_05252_row5_col3" class="data row5 col3">1.6581</td>
<td id="T_05252_row5_col4" class="data row5 col4">0.8579</td>
<td id="T_05252_row5_col5" class="data row5 col5">0.4964</td>
<td id="T_05252_row5_col6" class="data row5 col6">0.4625</td>
<td id="T_05252_row5_col7" class="data row5 col7">0.0800</td>
</tr>
<tr>
<td id="T_05252_level0_row6" class="row_heading level0 row6"
data-quarto-table-cell-role="th">dt</td>
<td id="T_05252_row6_col0" class="data row6 col0">Decision Tree
Regressor</td>
<td id="T_05252_row6_col1" class="data row6 col1">1.3407</td>
<td id="T_05252_row6_col2" class="data row6 col2">5.5824</td>
<td id="T_05252_row6_col3" class="data row6 col3">2.3627</td>
<td id="T_05252_row6_col4" class="data row6 col4">0.7114</td>
<td id="T_05252_row6_col5" class="data row6 col5">0.7269</td>
<td id="T_05252_row6_col6" class="data row6 col6">0.5754</td>
<td id="T_05252_row6_col7" class="data row6 col7">0.0100</td>
</tr>
<tr>
<td id="T_05252_level0_row7" class="row_heading level0 row7"
data-quarto-table-cell-role="th">ada</td>
<td id="T_05252_row7_col0" class="data row7 col0">AdaBoost
Regressor</td>
<td id="T_05252_row7_col1" class="data row7 col1">2.3201</td>
<td id="T_05252_row7_col2" class="data row7 col2">7.0186</td>
<td id="T_05252_row7_col3" class="data row7 col3">2.6493</td>
<td id="T_05252_row7_col4" class="data row7 col4">0.6372</td>
<td id="T_05252_row7_col5" class="data row7 col5">0.9424</td>
<td id="T_05252_row7_col6" class="data row7 col6">0.5857</td>
<td id="T_05252_row7_col7" class="data row7 col7">0.0700</td>
</tr>
<tr>
<td id="T_05252_level0_row8" class="row_heading level0 row8"
data-quarto-table-cell-role="th">ridge</td>
<td id="T_05252_row8_col0" class="data row8 col0">Ridge Regression</td>
<td id="T_05252_row8_col1" class="data row8 col1">2.1316</td>
<td id="T_05252_row8_col2" class="data row8 col2">7.4200</td>
<td id="T_05252_row8_col3" class="data row8 col3">2.7240</td>
<td id="T_05252_row8_col4" class="data row8 col4">0.6164</td>
<td id="T_05252_row8_col5" class="data row8 col5">0.9256</td>
<td id="T_05252_row8_col6" class="data row8 col6">0.9212</td>
<td id="T_05252_row8_col7" class="data row8 col7">0.0000</td>
</tr>
<tr>
<td id="T_05252_level0_row9" class="row_heading level0 row9"
data-quarto-table-cell-role="th">br</td>
<td id="T_05252_row9_col0" class="data row9 col0">Bayesian Ridge</td>
<td id="T_05252_row9_col1" class="data row9 col1">2.1351</td>
<td id="T_05252_row9_col2" class="data row9 col2">7.4732</td>
<td id="T_05252_row9_col3" class="data row9 col3">2.7337</td>
<td id="T_05252_row9_col4" class="data row9 col4">0.6137</td>
<td id="T_05252_row9_col5" class="data row9 col5">0.9270</td>
<td id="T_05252_row9_col6" class="data row9 col6">0.9264</td>
<td id="T_05252_row9_col7" class="data row9 col7">0.0100</td>
</tr>
<tr>
<td id="T_05252_level0_row10" class="row_heading level0 row10"
data-quarto-table-cell-role="th">lar</td>
<td id="T_05252_row10_col0" class="data row10 col0">Least Angle
Regression</td>
<td id="T_05252_row10_col1" class="data row10 col1">3.7772</td>
<td id="T_05252_row10_col2" class="data row10 col2">18.6936</td>
<td id="T_05252_row10_col3" class="data row10 col3">4.3236</td>
<td id="T_05252_row10_col4" class="data row10 col4">0.0337</td>
<td id="T_05252_row10_col5" class="data row10 col5">1.1358</td>
<td id="T_05252_row10_col6" class="data row10 col6">1.2991</td>
<td id="T_05252_row10_col7" class="data row10 col7">0.0100</td>
</tr>
<tr>
<td id="T_05252_level0_row11" class="row_heading level0 row11"
data-quarto-table-cell-role="th">huber</td>
<td id="T_05252_row11_col0" class="data row11 col0">Huber Regressor</td>
<td id="T_05252_row11_col1" class="data row11 col1">3.0573</td>
<td id="T_05252_row11_col2" class="data row11 col2">19.5931</td>
<td id="T_05252_row11_col3" class="data row11 col3">4.4264</td>
<td id="T_05252_row11_col4" class="data row11 col4">-0.0128</td>
<td id="T_05252_row11_col5" class="data row11 col5">0.8909</td>
<td id="T_05252_row11_col6" class="data row11 col6">0.6111</td>
<td id="T_05252_row11_col7" class="data row11 col7">0.0200</td>
</tr>
<tr>
<td id="T_05252_level0_row12" class="row_heading level0 row12"
data-quarto-table-cell-role="th">en</td>
<td id="T_05252_row12_col0" class="data row12 col0">Elastic Net</td>
<td id="T_05252_row12_col1" class="data row12 col1">4.0496</td>
<td id="T_05252_row12_col2" class="data row12 col2">19.6190</td>
<td id="T_05252_row12_col3" class="data row12 col3">4.4293</td>
<td id="T_05252_row12_col4" class="data row12 col4">-0.0142</td>
<td id="T_05252_row12_col5" class="data row12 col5">1.2366</td>
<td id="T_05252_row12_col6" class="data row12 col6">1.1734</td>
<td id="T_05252_row12_col7" class="data row12 col7">0.0000</td>
</tr>
<tr>
<td id="T_05252_level0_row13" class="row_heading level0 row13"
data-quarto-table-cell-role="th">par</td>
<td id="T_05252_row13_col0" class="data row13 col0">Passive Aggressive
Regressor</td>
<td id="T_05252_row13_col1" class="data row13 col1">3.8849</td>
<td id="T_05252_row13_col2" class="data row13 col2">19.8872</td>
<td id="T_05252_row13_col3" class="data row13 col3">4.4595</td>
<td id="T_05252_row13_col4" class="data row13 col4">-0.0280</td>
<td id="T_05252_row13_col5" class="data row13 col5">1.1862</td>
<td id="T_05252_row13_col6" class="data row13 col6">1.0219</td>
<td id="T_05252_row13_col7" class="data row13 col7">0.0000</td>
</tr>
<tr>
<td id="T_05252_level0_row14" class="row_heading level0 row14"
data-quarto-table-cell-role="th">lr</td>
<td id="T_05252_row14_col0" class="data row14 col0">Linear
Regression</td>
<td id="T_05252_row14_col1" class="data row14 col1">4.1256</td>
<td id="T_05252_row14_col2" class="data row14 col2">20.4303</td>
<td id="T_05252_row14_col3" class="data row14 col3">4.5200</td>
<td id="T_05252_row14_col4" class="data row14 col4">-0.0561</td>
<td id="T_05252_row14_col5" class="data row14 col5">1.2505</td>
<td id="T_05252_row14_col6" class="data row14 col6">1.2217</td>
<td id="T_05252_row14_col7" class="data row14 col7">0.0000</td>
</tr>
<tr>
<td id="T_05252_level0_row15" class="row_heading level0 row15"
data-quarto-table-cell-role="th">lasso</td>
<td id="T_05252_row15_col0" class="data row15 col0">Lasso
Regression</td>
<td id="T_05252_row15_col1" class="data row15 col1">4.1243</td>
<td id="T_05252_row15_col2" class="data row15 col2">20.4572</td>
<td id="T_05252_row15_col3" class="data row15 col3">4.5230</td>
<td id="T_05252_row15_col4" class="data row15 col4">-0.0575</td>
<td id="T_05252_row15_col5" class="data row15 col5">1.2509</td>
<td id="T_05252_row15_col6" class="data row15 col6">1.2182</td>
<td id="T_05252_row15_col7" class="data row15 col7">0.0000</td>
</tr>
<tr>
<td id="T_05252_level0_row16" class="row_heading level0 row16"
data-quarto-table-cell-role="th">llar</td>
<td id="T_05252_row16_col0" class="data row16 col0">Lasso Least Angle
Regression</td>
<td id="T_05252_row16_col1" class="data row16 col1">4.7760</td>
<td id="T_05252_row16_col2" class="data row16 col2">26.0867</td>
<td id="T_05252_row16_col3" class="data row16 col3">5.1075</td>
<td id="T_05252_row16_col4" class="data row16 col4">-0.3485</td>
<td id="T_05252_row16_col5" class="data row16 col5">1.3798</td>
<td id="T_05252_row16_col6" class="data row16 col6">1.5266</td>
<td id="T_05252_row16_col7" class="data row16 col7">0.0000</td>
</tr>
<tr>
<td id="T_05252_level0_row17" class="row_heading level0 row17"
data-quarto-table-cell-role="th">omp</td>
<td id="T_05252_row17_col0" class="data row17 col0">Orthogonal Matching
Pursuit</td>
<td id="T_05252_row17_col1" class="data row17 col1">4.1709</td>
<td id="T_05252_row17_col2" class="data row17 col2">28.0720</td>
<td id="T_05252_row17_col3" class="data row17 col3">5.2983</td>
<td id="T_05252_row17_col4" class="data row17 col4">-0.4511</td>
<td id="T_05252_row17_col5" class="data row17 col5">1.2779</td>
<td id="T_05252_row17_col6" class="data row17 col6">1.0241</td>
<td id="T_05252_row17_col7" class="data row17 col7">0.0000</td>
</tr>
<tr>
<td id="T_05252_level0_row18" class="row_heading level0 row18"
data-quarto-table-cell-role="th">knn</td>
<td id="T_05252_row18_col0" class="data row18 col0">K Neighbors
Regressor</td>
<td id="T_05252_row18_col1" class="data row18 col1">7.1319</td>
<td id="T_05252_row18_col2" class="data row18 col2">62.3846</td>
<td id="T_05252_row18_col3" class="data row18 col3">7.8984</td>
<td id="T_05252_row18_col4" class="data row18 col4">-2.2248</td>
<td id="T_05252_row18_col5" class="data row18 col5">1.7197</td>
<td id="T_05252_row18_col6" class="data row18 col6">2.7013</td>
<td id="T_05252_row18_col7" class="data row18 col7">0.0000</td>
</tr>
</tbody>
</table>

``` python
fig, error, result = predict_using_automl(df, future, best[0])
#plotly.offline.plot(fig);
```

<style  type="text/css" >
</style>

<table id="T_4b60d_" data-quarto-postprocess="true">
<thead>
<tr>
<th class="blank level0" data-quarto-table-cell-role="th"></th>
<th class="col_heading level0 col0"
data-quarto-table-cell-role="th">Model</th>
<th class="col_heading level0 col1"
data-quarto-table-cell-role="th">MAE</th>
<th class="col_heading level0 col2"
data-quarto-table-cell-role="th">MSE</th>
<th class="col_heading level0 col3"
data-quarto-table-cell-role="th">RMSE</th>
<th class="col_heading level0 col4"
data-quarto-table-cell-role="th">R2</th>
<th class="col_heading level0 col5"
data-quarto-table-cell-role="th">RMSLE</th>
<th class="col_heading level0 col6"
data-quarto-table-cell-role="th">MAPE</th>
</tr>
</thead>
<tbody>
<tr>
<td id="T_4b60d_level0_row0" class="row_heading level0 row0"
data-quarto-table-cell-role="th">0</td>
<td id="T_4b60d_row0_col0" class="data row0 col0">Extra Trees
Regressor</td>
<td id="T_4b60d_row0_col1" class="data row0 col1">0.8184</td>
<td id="T_4b60d_row0_col2" class="data row0 col2">1.4329</td>
<td id="T_4b60d_row0_col3" class="data row0 col3">1.1970</td>
<td id="T_4b60d_row0_col4" class="data row0 col4">0.9259</td>
<td id="T_4b60d_row0_col5" class="data row0 col5">0.3454</td>
<td id="T_4b60d_row0_col6" class="data row0 col6">0.3228</td>
</tr>
</tbody>
</table>

### 検証用データに対する誤差分布の解析

検証用データ期間に対して1日単位の誤差を求め，最も適合する連続分布を求める．

リード時間内の誤差の和の分布を求め，同様に最も適合する連続分布を求める．

``` python
train_idx, valid_idx = prepare_index(df, horizon)
print(len(valid_idx))
```

    90

``` python
error = result[:len(df)].Label - df.demand
error[valid_idx].hist(density=True);
```

![](91ex_inv_files/figure-commonmark/cell-29-output-1.png)

``` python
error[valid_idx]
```

    640   -0.88
    641   -1.10
    642   -0.18
    643    0.47
    644   -0.52
           ... 
    725    0.02
    726    0.01
    727   -0.85
    728    0.05
    729    0.06
    Length: 90, dtype: float64

``` python
#1日分の需要量（誤差）分布を求める
fig, best_dist, best_fit_name, best_fit_params  = best_distribution(error[valid_idx])
#plotly.offline.plot(fig);
print(best_fit_name, best_fit_params)
```

    dgamma (0.5146404593477012, 0.019999999999999997, 1.4358365014276395)

![](91ex_inv_files/figure-commonmark/cell-32-output-1.png)

``` python
d = best_dist.rvs((5,10000))
pd.Series(d.sum(axis=0)).hist(density=True);
```

![](91ex_inv_files/figure-commonmark/cell-33-output-1.png)

``` python
n = norm(loc=best_dist.mean()*5,scale=best_dist.std()*np.sqrt(5))
d2 = best_dist.rvs(10000)
pd.Series(d2).hist(density=True);
```

![](91ex_inv_files/figure-commonmark/cell-34-output-1.png)

``` python
#正味補充（リード）時間内の最大在庫量を計算
Lmax = 50
MaxDemand = [0.] 
for L in range(1,min(Lmax,len(data))): #リード時間
    df_error = pd.Series(error[valid_idx])
    data = df_error.rolling(window=L).sum()[L-1:].values #L日前から直前までの需要の合計
    fig, hist  = best_histogram(data)
    #plotly.offline.plot(fig2);
    MaxDemand.append( hist.ppf(0.95) )
```

``` python
#非減少にする
MaxDemand2 = []
max_ = MaxDemand[0]
for d in MaxDemand:
    max_ = max(max_, d)
    if d < max_: 
        MaxDemand2.append(max_)
    else:
        MaxDemand2.append(d)
pd.Series(MaxDemand2).plot();
```

![](91ex_inv_files/figure-commonmark/cell-36-output-1.png)

``` python
L = 5#リード時間
df_error = pd.Series(error[valid_idx])
data = df_error.rolling(window=L).sum()[L-1:].values #L日前から直前までの需要の合計
pd.Series(data).hist();
```

![](91ex_inv_files/figure-commonmark/cell-37-output-1.png)

``` python
fig2, best_dist2, best_fit_name2, best_fit_params2  = best_distribution(data)
plotly.offline.plot(fig2);
print(best_fit_name2, best_fit_params2)
```

    dgamma (0.7658916979282773, 0.10999999999999988, 1.766101465922829)

![](91ex_inv_files/figure-commonmark/cell-39-output-1.png)

### 基在庫レベルの設定

誤差分布から定常な需要分布を生成する．

リード時間内の誤差分布から，初期基在庫レベルを以下の式で生成する．

品切れ費用 *b*，在庫費用から臨界率を計算する．

*ω* = *b*/(*b* + *h*)

リード時間内の需要分布に対して，品切れ率が臨界率に一致するように安全在庫量を求める．

また，需要が負にならないように，それ以下になる確率が1%の値を 0
と設定することによって，需要を定数だけ大きくする．
毎日，この定数分だけの需要が発生すると仮定して，需要データを生成する．

この需要系列を利用して，定常状態の最適な基在庫レベルをシミュレーションに基づく最適化で行う．

``` python
n_samples = 10
n_periods = 1000
capacity = 10.
convergence = 1e-5
b = 100
h = 1 
LT = L
critical_ratio = b/(b+h)
S = best_dist2.ppf(critical_ratio) 
lb = best_dist.ppf(0.01)
print(S, lb )
demand = best_dist.rvs((n_samples,n_periods)) - lb #需要を分布が0以上になるようにシフトさせる．
S = S - lb*LT #基在庫レベルはリード時間内の定常需要だけ大きくする．
print(S)
```

    6.106113612640433 -3.923975750210528
    17.87804086327202

``` python
n_samples = 10
n_periods = 1000
capacity = 1000.
convergence = 1e-5
b = 100
h = 1 
LT = L
critical_ratio = b/(b+h)
S = best_dist2.ppf(critical_ratio) 
lb = best_dist.ppf(0.01)
print(S, lb )
demand = best_dist.rvs((n_samples,n_periods)) - lb #需要を分布が0以上になるようにシフトさせる．
S = S - lb*LT #基在庫レベルはリード時間内の定常需要だけ大きくする．
print(S)
print("t:   S      dS     Cost")
for iter_ in range(100):
    dC, cost, I = base_stock_simulation(n_samples, n_periods, demand, capacity, LT, b, h, S)
    S = S - .1*dC
    print(f"{iter_}: {S:.2f} {dC:.3f} {cost:.2f}")
    if dC**2<=convergence:
        break
```

    6.106113612640433 -3.923975750210528
    17.87804086327202
    t:   S      dS     Cost
    0: 17.87 0.091 7.54
    1: 17.86 0.081 7.54
    2: 17.85 0.071 7.54
    3: 17.85 0.071 7.54
    4: 17.84 0.071 7.54
    5: 17.83 0.061 7.54
    6: 17.83 0.061 7.54
    7: 17.82 0.061 7.54
    8: 17.82 0.061 7.54
    9: 17.81 0.051 7.54
    10: 17.81 0.051 7.54
    11: 17.80 0.051 7.54
    12: 17.80 0.030 7.54
    13: 17.79 0.030 7.54
    14: 17.79 0.010 7.54
    15: 17.79 0.010 7.54
    16: 17.79 0.010 7.54
    17: 17.79 0.010 7.54
    18: 17.79 0.010 7.54
    19: 17.79 0.010 7.54
    20: 17.79 0.010 7.54
    21: 17.79 0.010 7.54
    22: 17.78 0.010 7.54
    23: 17.78 0.010 7.54
    24: 17.78 0.010 7.53
    25: 17.78 0.010 7.53
    26: 17.78 0.010 7.53
    27: 17.78 0.000 7.53

``` python
pd.DataFrame(I[0]).plot(); #在庫の推移の可視化
```

![](91ex_inv_files/figure-commonmark/cell-42-output-1.png)

### 動的基在庫レベルの計算

定数分だけシフトした需要分を減じることによって，予測の誤差に対応するための安全在庫料量を計算する．

予測値は確定的な値として計算されるので，リード時間分の確定値の和だけ在庫を保持する．これに安全在庫分を加えたものが動的な基在庫レベルになる．

日々の運用は，在庫ポジションをこの値になるように発注（生産）を行う．

``` python
#安全在庫量
safety = S + lb*LT
print(safety)
```

    6.008873612640432

``` python
#未来の確定値の需要予測（リード時間分の累積値）
future = pd.Series(result[len(df):].Label)
data = future.rolling(window=L).sum()[L-1:].values #L日前から直前までの需要の合計
data
```

    array([13.55, 13.47, 10.12,  6.89, 14.99, 21.15, 20.65, 20.65])

``` python
#動的基在庫レベル 
data + safety
```

    array([19.55887361, 19.47887361, 16.12887361, 12.89887361, 20.99887361,
           27.15887361, 26.65887361, 26.65887361])

``` python
#検証データでシミュレーション
future = result.Label[valid_idx] 
data = future.rolling(window=L).sum()[L-1:].values #L日前から直前までの需要の合計
DS = data + safety #dynamic base stock level 
dem = df.demand[valid_idx]
```

``` python
I = np.zeros( len(dem)+1 )
I[0] = DS[0] #エシェロン在庫ポジション
for t, d in enumerate(dem):
    I[t+1] = I[t] - d
    order = max(DS[t+1] -I[t+1],0) 
    I[t+1] = I[t+1] + order
    #print(t,d,I[t+1])
    if t>=len(DS)-2:
        break
```

``` python
pd.Series(I[:-3]).plot();
```

![](91ex_inv_files/figure-commonmark/cell-48-output-1.png)

``` python
dem.plot();
```

![](91ex_inv_files/figure-commonmark/cell-49-output-1.png)

### 正規分布に近い需要関数の例

上の例ではzip分布に近い需要を仮定していた．以下では，正規分布に近い場合を考える．

``` python
demand_df = pd.read_csv(folder+"demand_normal.csv") 
c = "仙台市" 
p = "A"
print(c,p)
agg_period ="1d"
df, future = make_forecast_df(demand_df, customer=c, product=p,promo_df= None, agg_period="1d", forecast_periods = 0)
df.demand.max()
horizon="2019/10/01"
best, result_df = automl(df, horizon, 5)
```

<style  type="text/css" >
    #T_dc265_ th {
          text-align: left;
    }#T_dc265_row0_col0,#T_dc265_row1_col0,#T_dc265_row1_col1,#T_dc265_row1_col2,#T_dc265_row1_col3,#T_dc265_row1_col4,#T_dc265_row1_col5,#T_dc265_row1_col6,#T_dc265_row2_col0,#T_dc265_row2_col1,#T_dc265_row2_col2,#T_dc265_row2_col3,#T_dc265_row2_col4,#T_dc265_row2_col5,#T_dc265_row2_col6,#T_dc265_row3_col0,#T_dc265_row3_col1,#T_dc265_row3_col2,#T_dc265_row3_col3,#T_dc265_row3_col4,#T_dc265_row3_col5,#T_dc265_row3_col6,#T_dc265_row4_col0,#T_dc265_row4_col1,#T_dc265_row4_col2,#T_dc265_row4_col3,#T_dc265_row4_col4,#T_dc265_row4_col5,#T_dc265_row4_col6,#T_dc265_row5_col0,#T_dc265_row5_col1,#T_dc265_row5_col2,#T_dc265_row5_col3,#T_dc265_row5_col4,#T_dc265_row5_col5,#T_dc265_row5_col6,#T_dc265_row6_col0,#T_dc265_row6_col1,#T_dc265_row6_col2,#T_dc265_row6_col3,#T_dc265_row6_col4,#T_dc265_row6_col5,#T_dc265_row6_col6,#T_dc265_row7_col0,#T_dc265_row7_col1,#T_dc265_row7_col2,#T_dc265_row7_col3,#T_dc265_row7_col4,#T_dc265_row7_col5,#T_dc265_row7_col6,#T_dc265_row8_col0,#T_dc265_row8_col1,#T_dc265_row8_col2,#T_dc265_row8_col3,#T_dc265_row8_col4,#T_dc265_row8_col5,#T_dc265_row8_col6,#T_dc265_row9_col0,#T_dc265_row9_col1,#T_dc265_row9_col2,#T_dc265_row9_col3,#T_dc265_row9_col4,#T_dc265_row9_col5,#T_dc265_row9_col6,#T_dc265_row10_col0,#T_dc265_row10_col1,#T_dc265_row10_col2,#T_dc265_row10_col3,#T_dc265_row10_col4,#T_dc265_row10_col5,#T_dc265_row10_col6,#T_dc265_row11_col0,#T_dc265_row11_col1,#T_dc265_row11_col2,#T_dc265_row11_col3,#T_dc265_row11_col4,#T_dc265_row11_col5,#T_dc265_row11_col6,#T_dc265_row12_col0,#T_dc265_row12_col1,#T_dc265_row12_col2,#T_dc265_row12_col3,#T_dc265_row12_col4,#T_dc265_row12_col5,#T_dc265_row12_col6,#T_dc265_row13_col0,#T_dc265_row13_col1,#T_dc265_row13_col2,#T_dc265_row13_col3,#T_dc265_row13_col4,#T_dc265_row13_col5,#T_dc265_row13_col6,#T_dc265_row14_col0,#T_dc265_row14_col1,#T_dc265_row14_col2,#T_dc265_row14_col3,#T_dc265_row14_col4,#T_dc265_row14_col5,#T_dc265_row14_col6,#T_dc265_row15_col0,#T_dc265_row15_col1,#T_dc265_row15_col2,#T_dc265_row15_col3,#T_dc265_row15_col4,#T_dc265_row15_col5,#T_dc265_row15_col6,#T_dc265_row16_col0,#T_dc265_row16_col1,#T_dc265_row16_col2,#T_dc265_row16_col3,#T_dc265_row16_col4,#T_dc265_row16_col5,#T_dc265_row16_col6,#T_dc265_row17_col0,#T_dc265_row17_col1,#T_dc265_row17_col2,#T_dc265_row17_col3,#T_dc265_row17_col4,#T_dc265_row17_col5,#T_dc265_row17_col6,#T_dc265_row18_col0,#T_dc265_row18_col1,#T_dc265_row18_col2,#T_dc265_row18_col3,#T_dc265_row18_col4,#T_dc265_row18_col5,#T_dc265_row18_col6{
            text-align:  left;
            text-align:  left;
        }#T_dc265_row0_col1,#T_dc265_row0_col2,#T_dc265_row0_col3,#T_dc265_row0_col4,#T_dc265_row0_col5,#T_dc265_row0_col6{
            text-align:  left;
            text-align:  left;
            background-color:  yellow;
        }#T_dc265_row0_col7,#T_dc265_row3_col7,#T_dc265_row5_col7,#T_dc265_row6_col7,#T_dc265_row8_col7,#T_dc265_row9_col7,#T_dc265_row10_col7,#T_dc265_row11_col7,#T_dc265_row14_col7,#T_dc265_row16_col7{
            text-align:  left;
            text-align:  left;
            background-color:  lightgrey;
        }#T_dc265_row1_col7,#T_dc265_row2_col7,#T_dc265_row4_col7,#T_dc265_row7_col7,#T_dc265_row12_col7,#T_dc265_row13_col7,#T_dc265_row15_col7,#T_dc265_row17_col7,#T_dc265_row18_col7{
            text-align:  left;
            text-align:  left;
            background-color:  yellow;
            background-color:  lightgrey;
        }</style>

<table id="T_dc265_" data-quarto-postprocess="true">
<thead>
<tr>
<th class="blank level0" data-quarto-table-cell-role="th"></th>
<th class="col_heading level0 col0"
data-quarto-table-cell-role="th">Model</th>
<th class="col_heading level0 col1"
data-quarto-table-cell-role="th">MAE</th>
<th class="col_heading level0 col2"
data-quarto-table-cell-role="th">MSE</th>
<th class="col_heading level0 col3"
data-quarto-table-cell-role="th">RMSE</th>
<th class="col_heading level0 col4"
data-quarto-table-cell-role="th">R2</th>
<th class="col_heading level0 col5"
data-quarto-table-cell-role="th">RMSLE</th>
<th class="col_heading level0 col6"
data-quarto-table-cell-role="th">MAPE</th>
<th class="col_heading level0 col7" data-quarto-table-cell-role="th">TT
(Sec)</th>
</tr>
</thead>
<tbody>
<tr>
<td id="T_dc265_level0_row0" class="row_heading level0 row0"
data-quarto-table-cell-role="th">lasso</td>
<td id="T_dc265_row0_col0" class="data row0 col0">Lasso Regression</td>
<td id="T_dc265_row0_col1" class="data row0 col1">100.2628</td>
<td id="T_dc265_row0_col2" class="data row0 col2">15136.4229</td>
<td id="T_dc265_row0_col3" class="data row0 col3">123.0302</td>
<td id="T_dc265_row0_col4" class="data row0 col4">0.7174</td>
<td id="T_dc265_row0_col5" class="data row0 col5">0.0801</td>
<td id="T_dc265_row0_col6" class="data row0 col6">0.0656</td>
<td id="T_dc265_row0_col7" class="data row0 col7">0.0100</td>
</tr>
<tr>
<td id="T_dc265_level0_row1" class="row_heading level0 row1"
data-quarto-table-cell-role="th">ridge</td>
<td id="T_dc265_row1_col0" class="data row1 col0">Ridge Regression</td>
<td id="T_dc265_row1_col1" class="data row1 col1">115.1678</td>
<td id="T_dc265_row1_col2" class="data row1 col2">21126.2422</td>
<td id="T_dc265_row1_col3" class="data row1 col3">145.3487</td>
<td id="T_dc265_row1_col4" class="data row1 col4">0.6056</td>
<td id="T_dc265_row1_col5" class="data row1 col5">0.0939</td>
<td id="T_dc265_row1_col6" class="data row1 col6">0.0716</td>
<td id="T_dc265_row1_col7" class="data row1 col7">0.0000</td>
</tr>
<tr>
<td id="T_dc265_level0_row2" class="row_heading level0 row2"
data-quarto-table-cell-role="th">llar</td>
<td id="T_dc265_row2_col0" class="data row2 col0">Lasso Least Angle
Regression</td>
<td id="T_dc265_row2_col1" class="data row2 col1">124.3677</td>
<td id="T_dc265_row2_col2" class="data row2 col2">22642.7084</td>
<td id="T_dc265_row2_col3" class="data row2 col3">150.4749</td>
<td id="T_dc265_row2_col4" class="data row2 col4">0.5773</td>
<td id="T_dc265_row2_col5" class="data row2 col5">0.1000</td>
<td id="T_dc265_row2_col6" class="data row2 col6">0.0847</td>
<td id="T_dc265_row2_col7" class="data row2 col7">0.0000</td>
</tr>
<tr>
<td id="T_dc265_level0_row3" class="row_heading level0 row3"
data-quarto-table-cell-role="th">lightgbm</td>
<td id="T_dc265_row3_col0" class="data row3 col0">Light Gradient
Boosting Machine</td>
<td id="T_dc265_row3_col1" class="data row3 col1">135.1054</td>
<td id="T_dc265_row3_col2" class="data row3 col2">26361.0052</td>
<td id="T_dc265_row3_col3" class="data row3 col3">162.3607</td>
<td id="T_dc265_row3_col4" class="data row3 col4">0.5079</td>
<td id="T_dc265_row3_col5" class="data row3 col5">0.1054</td>
<td id="T_dc265_row3_col6" class="data row3 col6">0.0912</td>
<td id="T_dc265_row3_col7" class="data row3 col7">0.0400</td>
</tr>
<tr>
<td id="T_dc265_level0_row4" class="row_heading level0 row4"
data-quarto-table-cell-role="th">omp</td>
<td id="T_dc265_row4_col0" class="data row4 col0">Orthogonal Matching
Pursuit</td>
<td id="T_dc265_row4_col1" class="data row4 col1">130.6787</td>
<td id="T_dc265_row4_col2" class="data row4 col2">27112.9591</td>
<td id="T_dc265_row4_col3" class="data row4 col3">164.6601</td>
<td id="T_dc265_row4_col4" class="data row4 col4">0.4938</td>
<td id="T_dc265_row4_col5" class="data row4 col5">0.1072</td>
<td id="T_dc265_row4_col6" class="data row4 col6">0.0880</td>
<td id="T_dc265_row4_col7" class="data row4 col7">0.0000</td>
</tr>
<tr>
<td id="T_dc265_level0_row5" class="row_heading level0 row5"
data-quarto-table-cell-role="th">gbr</td>
<td id="T_dc265_row5_col0" class="data row5 col0">Gradient Boosting
Regressor</td>
<td id="T_dc265_row5_col1" class="data row5 col1">140.5485</td>
<td id="T_dc265_row5_col2" class="data row5 col2">28375.8903</td>
<td id="T_dc265_row5_col3" class="data row5 col3">168.4514</td>
<td id="T_dc265_row5_col4" class="data row5 col4">0.4703</td>
<td id="T_dc265_row5_col5" class="data row5 col5">0.1097</td>
<td id="T_dc265_row5_col6" class="data row5 col6">0.0951</td>
<td id="T_dc265_row5_col7" class="data row5 col7">0.0400</td>
</tr>
<tr>
<td id="T_dc265_level0_row6" class="row_heading level0 row6"
data-quarto-table-cell-role="th">catboost</td>
<td id="T_dc265_row6_col0" class="data row6 col0">CatBoost
Regressor</td>
<td id="T_dc265_row6_col1" class="data row6 col1">151.5904</td>
<td id="T_dc265_row6_col2" class="data row6 col2">32636.8823</td>
<td id="T_dc265_row6_col3" class="data row6 col3">180.6568</td>
<td id="T_dc265_row6_col4" class="data row6 col4">0.3907</td>
<td id="T_dc265_row6_col5" class="data row6 col5">0.1175</td>
<td id="T_dc265_row6_col6" class="data row6 col6">0.1030</td>
<td id="T_dc265_row6_col7" class="data row6 col7">1.2700</td>
</tr>
<tr>
<td id="T_dc265_level0_row7" class="row_heading level0 row7"
data-quarto-table-cell-role="th">en</td>
<td id="T_dc265_row7_col0" class="data row7 col0">Elastic Net</td>
<td id="T_dc265_row7_col1" class="data row7 col1">171.0394</td>
<td id="T_dc265_row7_col2" class="data row7 col2">42394.4023</td>
<td id="T_dc265_row7_col3" class="data row7 col3">205.8990</td>
<td id="T_dc265_row7_col4" class="data row7 col4">0.2086</td>
<td id="T_dc265_row7_col5" class="data row7 col5">0.1362</td>
<td id="T_dc265_row7_col6" class="data row7 col6">0.1163</td>
<td id="T_dc265_row7_col7" class="data row7 col7">0.0000</td>
</tr>
<tr>
<td id="T_dc265_level0_row8" class="row_heading level0 row8"
data-quarto-table-cell-role="th">et</td>
<td id="T_dc265_row8_col0" class="data row8 col0">Extra Trees
Regressor</td>
<td id="T_dc265_row8_col1" class="data row8 col1">177.6587</td>
<td id="T_dc265_row8_col2" class="data row8 col2">43153.4575</td>
<td id="T_dc265_row8_col3" class="data row8 col3">207.7341</td>
<td id="T_dc265_row8_col4" class="data row8 col4">0.1944</td>
<td id="T_dc265_row8_col5" class="data row8 col5">0.1327</td>
<td id="T_dc265_row8_col6" class="data row8 col6">0.1196</td>
<td id="T_dc265_row8_col7" class="data row8 col7">0.1400</td>
</tr>
<tr>
<td id="T_dc265_level0_row9" class="row_heading level0 row9"
data-quarto-table-cell-role="th">ada</td>
<td id="T_dc265_row9_col0" class="data row9 col0">AdaBoost
Regressor</td>
<td id="T_dc265_row9_col1" class="data row9 col1">185.1509</td>
<td id="T_dc265_row9_col2" class="data row9 col2">47749.3467</td>
<td id="T_dc265_row9_col3" class="data row9 col3">218.5162</td>
<td id="T_dc265_row9_col4" class="data row9 col4">0.1086</td>
<td id="T_dc265_row9_col5" class="data row9 col5">0.1402</td>
<td id="T_dc265_row9_col6" class="data row9 col6">0.1265</td>
<td id="T_dc265_row9_col7" class="data row9 col7">0.0500</td>
</tr>
<tr>
<td id="T_dc265_level0_row10" class="row_heading level0 row10"
data-quarto-table-cell-role="th">rf</td>
<td id="T_dc265_row10_col0" class="data row10 col0">Random Forest
Regressor</td>
<td id="T_dc265_row10_col1" class="data row10 col1">191.4697</td>
<td id="T_dc265_row10_col2" class="data row10 col2">51454.3014</td>
<td id="T_dc265_row10_col3" class="data row10 col3">226.8354</td>
<td id="T_dc265_row10_col4" class="data row10 col4">0.0394</td>
<td id="T_dc265_row10_col5" class="data row10 col5">0.1429</td>
<td id="T_dc265_row10_col6" class="data row10 col6">0.1296</td>
<td id="T_dc265_row10_col7" class="data row10 col7">0.1700</td>
</tr>
<tr>
<td id="T_dc265_level0_row11" class="row_heading level0 row11"
data-quarto-table-cell-role="th">xgboost</td>
<td id="T_dc265_row11_col0" class="data row11 col0">Extreme Gradient
Boosting</td>
<td id="T_dc265_row11_col1" class="data row11 col1">191.4150</td>
<td id="T_dc265_row11_col2" class="data row11 col2">53848.4531</td>
<td id="T_dc265_row11_col3" class="data row11 col3">232.0527</td>
<td id="T_dc265_row11_col4" class="data row11 col4">-0.0053</td>
<td id="T_dc265_row11_col5" class="data row11 col5">0.1472</td>
<td id="T_dc265_row11_col6" class="data row11 col6">0.1308</td>
<td id="T_dc265_row11_col7" class="data row11 col7">0.3300</td>
</tr>
<tr>
<td id="T_dc265_level0_row12" class="row_heading level0 row12"
data-quarto-table-cell-role="th">par</td>
<td id="T_dc265_row12_col0" class="data row12 col0">Passive Aggressive
Regressor</td>
<td id="T_dc265_row12_col1" class="data row12 col1">191.9421</td>
<td id="T_dc265_row12_col2" class="data row12 col2">54763.1891</td>
<td id="T_dc265_row12_col3" class="data row12 col3">234.0154</td>
<td id="T_dc265_row12_col4" class="data row12 col4">-0.0224</td>
<td id="T_dc265_row12_col5" class="data row12 col5">0.1515</td>
<td id="T_dc265_row12_col6" class="data row12 col6">0.1276</td>
<td id="T_dc265_row12_col7" class="data row12 col7">0.0000</td>
</tr>
<tr>
<td id="T_dc265_level0_row13" class="row_heading level0 row13"
data-quarto-table-cell-role="th">br</td>
<td id="T_dc265_row13_col0" class="data row13 col0">Bayesian Ridge</td>
<td id="T_dc265_row13_col1" class="data row13 col1">195.9908</td>
<td id="T_dc265_row13_col2" class="data row13 col2">57561.6403</td>
<td id="T_dc265_row13_col3" class="data row13 col3">239.9201</td>
<td id="T_dc265_row13_col4" class="data row13 col4">-0.0746</td>
<td id="T_dc265_row13_col5" class="data row13 col5">0.1564</td>
<td id="T_dc265_row13_col6" class="data row13 col6">0.1327</td>
<td id="T_dc265_row13_col7" class="data row13 col7">0.0000</td>
</tr>
<tr>
<td id="T_dc265_level0_row14" class="row_heading level0 row14"
data-quarto-table-cell-role="th">huber</td>
<td id="T_dc265_row14_col0" class="data row14 col0">Huber Regressor</td>
<td id="T_dc265_row14_col1" class="data row14 col1">196.6167</td>
<td id="T_dc265_row14_col2" class="data row14 col2">58014.9136</td>
<td id="T_dc265_row14_col3" class="data row14 col3">240.8629</td>
<td id="T_dc265_row14_col4" class="data row14 col4">-0.0831</td>
<td id="T_dc265_row14_col5" class="data row14 col5">0.1570</td>
<td id="T_dc265_row14_col6" class="data row14 col6">0.1332</td>
<td id="T_dc265_row14_col7" class="data row14 col7">0.0100</td>
</tr>
<tr>
<td id="T_dc265_level0_row15" class="row_heading level0 row15"
data-quarto-table-cell-role="th">lr</td>
<td id="T_dc265_row15_col0" class="data row15 col0">Linear
Regression</td>
<td id="T_dc265_row15_col1" class="data row15 col1">197.9349</td>
<td id="T_dc265_row15_col2" class="data row15 col2">58133.7969</td>
<td id="T_dc265_row15_col3" class="data row15 col3">241.1095</td>
<td id="T_dc265_row15_col4" class="data row15 col4">-0.0853</td>
<td id="T_dc265_row15_col5" class="data row15 col5">0.1571</td>
<td id="T_dc265_row15_col6" class="data row15 col6">0.1340</td>
<td id="T_dc265_row15_col7" class="data row15 col7">0.0000</td>
</tr>
<tr>
<td id="T_dc265_level0_row16" class="row_heading level0 row16"
data-quarto-table-cell-role="th">lar</td>
<td id="T_dc265_row16_col0" class="data row16 col0">Least Angle
Regression</td>
<td id="T_dc265_row16_col1" class="data row16 col1">212.4354</td>
<td id="T_dc265_row16_col2" class="data row16 col2">61619.4643</td>
<td id="T_dc265_row16_col3" class="data row16 col3">248.2327</td>
<td id="T_dc265_row16_col4" class="data row16 col4">-0.1503</td>
<td id="T_dc265_row16_col5" class="data row16 col5">0.1696</td>
<td id="T_dc265_row16_col6" class="data row16 col6">0.1315</td>
<td id="T_dc265_row16_col7" class="data row16 col7">0.0100</td>
</tr>
<tr>
<td id="T_dc265_level0_row17" class="row_heading level0 row17"
data-quarto-table-cell-role="th">dt</td>
<td id="T_dc265_row17_col0" class="data row17 col0">Decision Tree
Regressor</td>
<td id="T_dc265_row17_col1" class="data row17 col1">223.6154</td>
<td id="T_dc265_row17_col2" class="data row17 col2">71599.5055</td>
<td id="T_dc265_row17_col3" class="data row17 col3">267.5808</td>
<td id="T_dc265_row17_col4" class="data row17 col4">-0.3367</td>
<td id="T_dc265_row17_col5" class="data row17 col5">0.1597</td>
<td id="T_dc265_row17_col6" class="data row17 col6">0.1475</td>
<td id="T_dc265_row17_col7" class="data row17 col7">0.0000</td>
</tr>
<tr>
<td id="T_dc265_level0_row18" class="row_heading level0 row18"
data-quarto-table-cell-role="th">knn</td>
<td id="T_dc265_row18_col0" class="data row18 col0">K Neighbors
Regressor</td>
<td id="T_dc265_row18_col1" class="data row18 col1">235.7868</td>
<td id="T_dc265_row18_col2" class="data row18 col2">79601.6484</td>
<td id="T_dc265_row18_col3" class="data row18 col3">282.1376</td>
<td id="T_dc265_row18_col4" class="data row18 col4">-0.4861</td>
<td id="T_dc265_row18_col5" class="data row18 col5">0.1833</td>
<td id="T_dc265_row18_col6" class="data row18 col6">0.1647</td>
<td id="T_dc265_row18_col7" class="data row18 col7">0.0000</td>
</tr>
</tbody>
</table>

``` python
fig, error, result = predict_using_automl(df, future, best[0])
plotly.offline.plot(fig);
```

<style  type="text/css" >
</style>

<table id="T_602fc_" data-quarto-postprocess="true">
<thead>
<tr>
<th class="blank level0" data-quarto-table-cell-role="th"></th>
<th class="col_heading level0 col0"
data-quarto-table-cell-role="th">Model</th>
<th class="col_heading level0 col1"
data-quarto-table-cell-role="th">MAE</th>
<th class="col_heading level0 col2"
data-quarto-table-cell-role="th">MSE</th>
<th class="col_heading level0 col3"
data-quarto-table-cell-role="th">RMSE</th>
<th class="col_heading level0 col4"
data-quarto-table-cell-role="th">R2</th>
<th class="col_heading level0 col5"
data-quarto-table-cell-role="th">RMSLE</th>
<th class="col_heading level0 col6"
data-quarto-table-cell-role="th">MAPE</th>
</tr>
</thead>
<tbody>
<tr>
<td id="T_602fc_level0_row0" class="row_heading level0 row0"
data-quarto-table-cell-role="th">0</td>
<td id="T_602fc_row0_col0" class="data row0 col0">Lasso Regression</td>
<td id="T_602fc_row0_col1" class="data row0 col1">100.2628</td>
<td id="T_602fc_row0_col2" class="data row0 col2">15136.4229</td>
<td id="T_602fc_row0_col3" class="data row0 col3">123.0302</td>
<td id="T_602fc_row0_col4" class="data row0 col4">0.7174</td>
<td id="T_602fc_row0_col5" class="data row0 col5">0.0801</td>
<td id="T_602fc_row0_col6" class="data row0 col6">0.0656</td>
</tr>
</tbody>
</table>

``` python
train_idx, valid_idx = prepare_index(df, horizon)
print(len(valid_idx))
error = result.Label - df.demand
error[valid_idx].hist(density=True);
```

    91

![](91ex_inv_files/figure-commonmark/cell-52-output-2.png)

``` python
fig, best_dist, best_fit_name, best_fit_params  = best_distribution(error[valid_idx])
plotly.offline.plot(fig);
print(best_fit_name, best_fit_params)
```

    mielke (2.6729912087197274, 11.291010366009026, -373.2172202339409, 503.0204118518145)

![](91ex_inv_files/figure-commonmark/cell-54-output-1.png)

``` python
L = 3 #リード時間
df_error = pd.Series(error[valid_idx])
data = df_error.rolling(window=L).sum()[L-1:].values #L日前から直前までの需要の合計
```

``` python
fig2, best_dist2, best_fit_name2, best_fit_params2  = best_distribution(data)
plotly.offline.plot(fig2);
print(best_fit_name2, best_fit_params2)
```

    gausshyper (0.5815996180621221, 3.799377992229532, -5.653814263450716, 2.624922751643501, -605.2475585937501, 1374.3047245376556)

![](91ex_inv_files/figure-commonmark/cell-57-output-1.png)

### 基在庫レベルの設定

``` python
n_samples = 10
n_periods = 1000
capacity = 1000.
convergence = 1e-5
b = 100
h = 1 
LT = L
critical_ratio = b/(b+h)
S = best_dist2.ppf(critical_ratio) 
lb = best_dist.ppf(0.01)
print(S, lb )
demand = best_dist.rvs((n_samples,n_periods)) - lb #需要を分布が0以上になるようにシフトさせる．
S = S - lb*LT #基在庫レベルはリード時間内の定常需要だけ大きくする．
print(S)
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
```

    585.0727636548719 -283.39988265249485
    1435.2724116123563
    t:   S      dS     Cost
    0: 1428.10 0.717 580.84
    1: 1421.03 0.707 575.75
    2: 1413.96 0.707 570.74
    3: 1407.49 0.646 566.01
    4: 1401.53 0.596 562.02
    5: 1395.67 0.586 558.49
    6: 1390.22 0.545 555.13
    7: 1385.27 0.495 552.26
    8: 1381.03 0.424 549.95
    9: 1377.49 0.354 548.29
    10: 1374.26 0.323 547.06
    11: 1371.13 0.313 546.03
    12: 1368.40 0.273 545.12
    13: 1365.77 0.263 544.37
    14: 1363.24 0.253 543.70
    15: 1360.92 0.232 543.10
    16: 1358.90 0.202 542.59
    17: 1357.18 0.172 542.21
    18: 1355.56 0.162 541.92
    19: 1354.05 0.152 541.67
    20: 1352.73 0.131 541.46
    21: 1351.42 0.131 541.29
    22: 1350.41 0.101 541.13
    23: 1349.70 0.071 541.05
    24: 1348.99 0.071 541.00
    25: 1348.29 0.071 540.95
    26: 1347.88 0.041 540.91
    27: 1347.48 0.041 540.89
    28: 1347.17 0.030 540.88
    29: 1346.87 0.030 540.87
    30: 1346.66 0.020 540.86
    31: 1346.46 0.020 540.85
    32: 1346.26 0.020 540.85
    33: 1346.16 0.010 540.85
    34: 1346.05 0.010 540.85
    35: 1345.95 0.010 540.84
    36: 1345.85 0.010 540.84
    37: 1345.75 0.010 540.84
    38: 1345.65 0.010 540.84
    39: 1345.54 0.010 540.84
    40: 1345.44 0.010 540.84
    41: 1345.34 0.010 540.84
    42: 1345.24 0.010 540.84
    43: 1345.14 0.010 540.84
    44: 1345.14 0.000 540.83

``` python
pd.DataFrame(I[0]).plot(); #在庫の推移の可視化
```

![](91ex_inv_files/figure-commonmark/cell-59-output-1.png)

### 動的基在庫レベルの設定と検証データでシミュレーション

``` python
#安全在庫量
safety = S + lb*LT
print(safety)
#検証データでシミュレーション
future = result.Label[valid_idx] 
data = future.rolling(window=L).sum()[L-1:].values #L日前から直前までの需要の合計
DS = data + safety #dynamic base stock level 
dem = df.demand[valid_idx]
I = np.zeros( len(dem)+1 )
I[0] = DS[0] #エシェロン在庫ポジション
for t, d in enumerate(dem):
    I[t+1] = I[t] - d
    order = max(DS[t+1] -I[t+1],0) 
    I[t+1] = I[t+1] + order
    #print(t,d,I[t+1])
    if t>=len(DS)-2:
        break
pd.Series(I[:-3]).plot();
```

    494.93576365487047

![](91ex_inv_files/figure-commonmark/cell-60-output-2.png)

## ネットワーク型の例

``` python
#7点の一般型ネットワークの例題
n = 7
G = SCMGraph()
for i in range(n):
    G.add_node(i)
G.add_edges_from([(0, 2), (1, 2), (2,4), (3,4), (4,5), (4,6)])

z = np.full(len(G),1.65)
h = np.array([1,1,3,1,5,6,6])
mu = np.array([200,200,200,200,200,100,100])
sigma = np.array([14.1,14.1,14.1,14.1,14.1,10,10])
LTUB = np.array([0,0,0,0,0,3,1], int)

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

best_cost, best_sol, best_NRT, best_MaxLI, best_MinLT  = tabu_search_for_SSA(G, ProcTime,  LTUB, z, mu, sigma, h, max_iter = 10, TLLB =1, TLUB =3, seed = 1)
print("最良値", best_cost)
print("最良解", best_sol)
print("正味補充時間", best_NRT)
print("最大補充リード時間", best_MaxLI)
print("最小保証リード時間", best_MinLT)

pos = G.layout()
stage_df, bom_df = make_df_for_SSA(G, pos, ProcTime, LTUB, z, mu, sigma, h, best_NRT, best_MaxLI, best_MinLT)
stage_df.to_csv(folder_bom + "stage_ex1.csv")
bom_df.to_csv(folder_bom + "bom_ex1.csv")

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

    最良値 514.8330943986502
    最良解 [1 1 0 0 1 0 1]
    正味補充時間 [6. 2. 0. 0. 6. 0. 2.]
    最大補充リード時間 [6. 2. 3. 3. 6. 3. 3.]
    最小保証リード時間 [0. 0. 3. 3. 0. 3. 1.]

![](91ex_inv_files/figure-commonmark/cell-61-output-2.png)

``` python
stage_df = pd.read_csv(folder_bom + "stage_ex1.csv")
bom_df = pd.read_csv(folder_bom + "bom_ex1.csv")
best_cost, stage_df, bom_df, fig = solve_SSA(stage_df, bom_df)
stage_df
```

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }
&#10;    .dataframe tbody tr th {
        vertical-align: top;
    }
&#10;    .dataframe thead th {
        text-align: right;
    }
</style>

<table class="dataframe" data-quarto-postprocess="true" data-border="1">
<thead>
<tr style="text-align: right;">
<th data-quarto-table-cell-role="th"></th>
<th data-quarto-table-cell-role="th">Unnamed: 0</th>
<th data-quarto-table-cell-role="th">name</th>
<th data-quarto-table-cell-role="th">net_replenishment_time</th>
<th data-quarto-table-cell-role="th">max_guaranteed_LT</th>
<th data-quarto-table-cell-role="th">processing_time</th>
<th data-quarto-table-cell-role="th">replenishment_LT</th>
<th data-quarto-table-cell-role="th">guaranteed_LT</th>
<th data-quarto-table-cell-role="th">z</th>
<th data-quarto-table-cell-role="th">average_demand</th>
<th data-quarto-table-cell-role="th">sigma</th>
<th data-quarto-table-cell-role="th">h</th>
<th data-quarto-table-cell-role="th">b</th>
<th data-quarto-table-cell-role="th">capacity</th>
<th data-quarto-table-cell-role="th">x</th>
<th data-quarto-table-cell-role="th">y</th>
</tr>
</thead>
<tbody>
<tr>
<td data-quarto-table-cell-role="th">0</td>
<td>0</td>
<td>0</td>
<td>6.0</td>
<td>0</td>
<td>6</td>
<td>6.0</td>
<td>0.0</td>
<td>1.65</td>
<td>200</td>
<td>14.1</td>
<td>1</td>
<td>8.8</td>
<td>2000.0</td>
<td>0</td>
<td>0</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">1</td>
<td>1</td>
<td>1</td>
<td>2.0</td>
<td>0</td>
<td>2</td>
<td>2.0</td>
<td>0.0</td>
<td>1.65</td>
<td>200</td>
<td>14.1</td>
<td>1</td>
<td>8.8</td>
<td>2000.0</td>
<td>0</td>
<td>2</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">2</td>
<td>2</td>
<td>2</td>
<td>0.0</td>
<td>0</td>
<td>3</td>
<td>3.0</td>
<td>3.0</td>
<td>1.65</td>
<td>200</td>
<td>14.1</td>
<td>3</td>
<td>26.3</td>
<td>2000.0</td>
<td>1</td>
<td>0</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">3</td>
<td>3</td>
<td>3</td>
<td>0.0</td>
<td>0</td>
<td>3</td>
<td>3.0</td>
<td>3.0</td>
<td>1.65</td>
<td>200</td>
<td>14.1</td>
<td>1</td>
<td>8.8</td>
<td>2000.0</td>
<td>0</td>
<td>1</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">4</td>
<td>4</td>
<td>4</td>
<td>6.0</td>
<td>0</td>
<td>3</td>
<td>6.0</td>
<td>0.0</td>
<td>1.65</td>
<td>200</td>
<td>14.1</td>
<td>5</td>
<td>43.9</td>
<td>2000.0</td>
<td>2</td>
<td>0</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">5</td>
<td>5</td>
<td>5</td>
<td>0.0</td>
<td>3</td>
<td>3</td>
<td>3.0</td>
<td>3.0</td>
<td>1.65</td>
<td>100</td>
<td>10.0</td>
<td>6</td>
<td>52.7</td>
<td>1000.0</td>
<td>3</td>
<td>0</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">6</td>
<td>6</td>
<td>6</td>
<td>2.0</td>
<td>1</td>
<td>3</td>
<td>3.0</td>
<td>1.0</td>
<td>1.65</td>
<td>100</td>
<td>10.0</td>
<td>6</td>
<td>52.7</td>
<td>1000.0</td>
<td>3</td>
<td>1</td>
</tr>
</tbody>
</table>

</div>

``` python
cost, stage_df, I, cost_list, phi_list = periodic_inv_opt_fit_one_cycle(stage_df, bom_df, max_iter = 1000, n_samples = 10, 
                                                                        n_periods = 100, seed = 1, lr_find=False, max_lr =6.0, moms=(0.85,0.95))
```

    ELT= [21. 17. 13. 13. 11.  1.  3.]
    S= [4306.61362354 3495.92405238 2683.88315042 2683.88315042 2277.16127575
      116.5         328.57883832]
    t: Cost    |dC|      S 
    0: 16082.26, 2384121.79078 [4306.61362354 3495.92405238 2683.88315042 2683.88315042 2277.16127575
      116.5         328.57883832]
    Best: 16054.62 [4294.07257219 3509.66262327 2667.975591   2687.85892252 2286.37514486
      124.01255267  341.92807231]

![](91ex_inv_files/figure-commonmark/cell-64-output-1.png)