# ケーススタディ


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

## ABC分析

注意：このケーススタディは，幾つかの企業のコンサルティングから抽出したもので，特定の企業を想定したものではない．データや登場人物は，架空のものである．

------------------------------------------------------------------------

あなたは新たに結成されたマーケティングチームに配属された新人社員だ．新しい上司から，需要と製品の膨大なデータを渡されて，分析するようにと命令された．
さて，何から始めれば良いだろうか？

まずは，データの概要を調べることから始める．
これをデータ探索と呼ぶ．本来はその後でデータの前処理を行うが，ここではSCMOPTのランダムデータ生成関数を用いる．

そのため，データの概要は以下のように既知である．

- 1000個の顧客と100個の製品があり，1年分の需要が与えられている．
  顧客と製品はパレート分布にしたがい，製品の方が偏った分布になっている．
- 全部で36500000の需要データがあり， 年次と週次の周期がある．
- 欠損値はない．

資料とデータへのリンク

https://www.dropbox.com/sh/pnkdixxsf8i3o8t/AADU1udy_D7IZf4PQBL0dg1Wa?dl=0

``` python
#ランダムで大規模な顧客需要の生成
random_seed = 1
cust_df = generate_cust(num_locations = 1000, random_seed = random_seed, prefecture = None, no_island=True)
prod_df = generate_many_prod(100, weight_bound=(1, 5), cust_value_bound =(5,10),
                              dc_value_bound=(2, 2), plnt_value_bound=(1, 1),
                              fc_bound=(100000,200000), random_seed=random_seed)
weekly_ratio = [0.0, 1.0, 1.2, 1.3, 0.9, 1.5, 0.2]  # 0 means Sunday
yearly_ratio = [1.0 + np.sin(i)*0.5 for i in range(13)]  # 0 means January
demand_df = generate_demand(cust_df, prod_df, cust_shape=1., prod_shape=0.9, weekly_ratio=weekly_ratio, yearly_ratio=yearly_ratio,
                            start="2019/01/01", periods=365, epsilon=1., random_seed=random_seed)
```

``` python
#可視化で確認
fig = plot_cust(cust_df)
plotly.offline.plot(fig);
```

顧客数1000，製品数100というのはそれほど大規模な問題例ではない．普通のメーカーのサブブランド程度の規模感だ．
それでも（たった）1年分の需要データで， 36500000 件にもなる．
csvファイルに保存しても，Excelでは開くことさえ出来ない．

サプライ・チェイン分析は，最初にこのような大きな問題例を，本質を損なわず，かつ取り扱いが容易になるようにすることから始めるのが良い．

通常の（実際の）需要データは，パレートの法則（全体の数値の大部分は、全体を構成するうちの一部の要素が生み出しているという理論。別名、80:20の法則、もう1つの別名、ばらつきの法則）に従う．
これを利用して取り扱う製品数を減らすための方法が，ABC分析である．

SCMOPTに含まれる基本情報分析システムSCBASにはABC分析をするための関数が準備されているので，それを使うことにする．
製品についてはAカテゴリーに8割の需要をもつものを入れたいので，閾値に
0.8, 0.1, 0.1 と入力する．

``` python
fig_prod, fig_cust, agg_df_prod, agg_df_cust, new_df = generate_figures_for_abc_analysis(
    demand_df, cumsum = False, cust_thres="0.7, 0.2, 0.1", prod_thres="0.8, 0.1, 0.1")
```

カテゴリーAに含まれる25個の製品だけで全体の80%以上をカバーしていることが分かる．
Aカテゴリーの製品だけを抜き出して保管しておく．
以降では，これらの製品に絞って解析を行う．

``` python
try:
    new_df.reset_index(inplace=True)
except:
    pass
demand_df = new_df[new_df.prod_ABC=="A"].copy() #copyをしないと後でwarningが出る
```

``` python
demand_df.head()
```

製品データフレームにABCなどの情報を入れて，Aカテゴリーの製品だけを抜き出しておく．

``` python
agg_df_prod.reset_index(inplace=True)
new_prod_df = add_abc(prod_df, agg_df_prod)
new_prod_df.head()
```

``` python
prod_df = new_prod_df[ new_prod_df.abc =="A" ].copy()
prod_df.head()
```

## 顧客の集約

25個の製品で全体の80%以上をカバーしている．
顧客数1000は解析するには多すぎるので，地点の近さによって集約を行う．

需要の集約を行い，顧客の需要量の合計を地点の重みとして，重み付きのクラスタリングを行う．

``` python
total_demand_df, demand_cust_df, demand_prod_df  = make_total_demand(demand_df, start="2019/01/01", finish="2050/12/31", num_of_days=365)
```

``` python
weight = demand_cust_df.demand.values
```

顧客間の移動時間と距離をOSRMを用いて計算しておく．

``` python
durations,  distances = compute_durations(cust_df)
```

集約する点の数を色々変えてkmean法でクラスタリングを行う．評価値（各地点から集約地点への重み付き距離の和）が変化しなくなる集約数を求める．
それには，エルボー（肘）法を用いる．

後で集約した地点を候補地点として実際の倉庫を決めるので，曲げた肘に相当する場所（施設数が30あたり）ではなく，少し大きめの50くらいが良い．
ここでは，集約数を50としてクラスタリングを行う．

集約のための方法としては，kmeans法を使うべきではない．kmeans法は高速であるためエルボー法と相性が良いが，クラスタリングは重み付き重心を求めているにすぎない．重心は距離の自乗なので，本来の重み付き距離の和を評価値にした場合とは異なった地点に収束する．重み付き直線（大圏）距離の和を最小にするための方法として反復Weiszfeld法と呼ばれる解法がある．これは不動点アルゴリズムを利用したもので，kmeans法より遅いが，より正確な地点を算出する．もう一つの方法として，階層的クラスタリングを用いたものがある．この方法は直線（大圏）距離を用いるのではなく，実際の道路距離を用いることができる．ここでは，地図から計算した本当の移動時間を尺度としてクラスタリングを行うことにする．

その結果を地図上に描画すると以下のようになる．

``` python
X, Y, partition, cost = hierarchical_clusterning(cust_df, weight, durations, num_of_facilities = 50, linkage="complete")
print(cost)
```

ピンクの点がオリジナルの顧客で，赤い点が集約された新しい顧客群（クラスター）である．
点の大きさが需要量の大きさを表している．点の大きさが需要量の大きさを表している．
これをもとに顧客と需要データの集約を行う．

``` python
aggregated_cust_df, aggregated_total_demand_df =  make_aggregated_df(cust_df, demand_df, total_demand_df, X, Y, partition, weight)
#aggregated_cust_df.head()
```

## ロジスティックネットワーク設計

ロジスティックネットワーク設計問題を解くための準備として，顧客と同じ位置に倉庫を生成しておこう．

``` python
aggregated_dc_df = generate_dc(aggregated_cust_df, aggregated_total_demand_df, prod_df, num_dc=10, lb_ratio=0.0, ub_ratio=1.3,vc_bounds=(0.0,0.1),fc_bounds=(1e+6,1e+7))
```

工場と生産データを生成する．

``` python
plnt_df, plnt_prod_df = generate_plnt(prod_df, aggregated_total_demand_df, lead_time_bound= (1,2))
```

ネットワークを生成して図示してみる．
直線距離ではなく，道路距離を用いる． 距離と移動時間の計算には，OSRM
(Open Streetmap Routing Machine) を使う．
あまり遠い点同士は使わないので，工場と倉庫間は1万km，倉庫と顧客間は300km以下の経路だけを生成する．

また，工場・倉庫間は20トントラック，倉庫・顧客間は4トントラックを使用していると仮定し，1kmあたりの費用をそれぞれ20円，10円，１時間あたりの費用を両方とも8000円と仮定して，1kgあたりの輸送費用を計算する．

最適化を行うためのネットワークは以下のようになる．オレンジ色の点が工場で，赤色が顧客である．簡単のため，ここでは倉庫を立地するための候補地点は顧客の位置と同じと仮定している．

``` python
durations, distances, node_df = compute_durations(aggregated_cust_df, plnt_df)
```

``` python
trans_df, graph, position  = make_network_using_road(aggregated_cust_df, aggregated_dc_df, plnt_df, durations, distances, plnt_dc_threshold = 10000., dc_cust_threshold = 300. )
```

需要量の桁数（小数点以下も）確認しておく．これがあまりにも大きいときには，スケーリングや数値の丸めを行う必要がある．
この場合は，需要の数値は最大6桁であり，小数点以下は10桁以上と長い．
最適化が困難な場合には，数値丸めを行う必要がある．

ロジスティックネットワーク設計問題を解く．単一ソース制約の有無で，得られるネットワーク（解）が変わる．
単一ソース制約がつくと，顧客に割り当てられる倉庫が1つになる．これによって，複数の倉庫から配送されることによるオペレーションの煩雑さが避けられる．
単一ソース制約がないと，複数の倉庫から配送されるが，それによって全体の費用は安くなる．
会社の方針によって決める必要があるが，ここでは両方で解いて比較してみよう．

``` python
aggregated_flow_df, dc_df, cost_df, violation_df, status = solve_lnd(prod_df, aggregated_cust_df, aggregated_dc_df, plnt_df, plnt_prod_df, aggregated_total_demand_df, 
                           trans_df, dc_num= (0,13), single_sourcing=True, max_cpu = 10)
print("total cost=", cost_df.value.sum())
cost_df
```

``` python
aggregated_flow_df, aggregated_dc_df, cost_df, violation_df, status = solve_lnd(prod_df, aggregated_cust_df, aggregated_dc_df, plnt_df, plnt_prod_df, aggregated_total_demand_df, 
                           trans_df, dc_num= (0,12), single_sourcing=False, max_cpu = 60)
print("total cost=", cost_df.value.sum())
cost_df
```

この場合は，現場のオペレーションの容易さから，単一ソース制約ありを採用することにする．
得られたネットワークは以下のようになる．オレンジ色の点が工場で，赤色が開設した倉庫で，青色が顧客である．

### 続きは．．．他の最適化モデルとの連携

最適化をしただけで終わりではない．まずは，色々なケースを考えて最適化を行い，意思決定の補助とすることが重要だ．これはWhat
If分析と呼ばれる．たとえば，今後，人経費の高騰により倉庫での積替え費用が上昇することが見込まれるときに，どのような結果になるかをシミュレーションすることができる．また，自動運転車の導入により，遠距離のトラック輸送が容易になるシナリオも考えられる．そのような場合に，どのように対処すればよいのか，もしくは今からどのような準備をしておけば良いのかを最適化を部品として検討できるのである．

他のモデルとの連携も重要な課題である．

ロジスティックネットワーク設計モデルは，配送費用を単純な輸送費用として近似的に扱っている．実際には，トラックは複数の顧客を巡回輸送する．これには，別の最適化モデルである配送計画問題を解く必要がある．SCMOPTには，METRO
IVという高性能の配送最適化ソルバーが含まれているので，それとデータのやり取りをして，より精度の高い最適化モデルを構築できる．

また，上のロジスティックネットワーク設計モデルには，在庫の概念が含まれていなかった．（倉庫の容量制約としては考慮されていたが，費用としては輸送費用や積替え費用に含めた「線形な」費用として扱っていた．）実際には，在庫費用は非線形な費用であり，それを直接，ネットワーク設計モデルに組み込むことは難しい．幸い，SCMOPTには幾つかの在庫最適化モジュールが含まれている．それらとデータを介してやり取りをすることによって，在庫費用も含めた全体最適化が可能になるのである．

配送計画 -
デポごとに最適化をして，その費用を配分したものを本当の配送費用と考えて，ロジスティックネットワーク設計問題を再び解く．

在庫計画 -
輸送中も在庫地点として考えて，そこには在庫をおけないと固定；在庫を入れたロジスティックネットワーク設計問題を再び解く．

需要予測 -
今までの実績需要ではなく，先の需要を予測し，それにあったネットワーク設計を行う．

## 需要予測

あなたは，社内全体にプレゼンをしたところ，過去の需要に対する計画は役に立たないとの指摘を受けた．
将来の需要予測をしようと考えたが，
今，与えられているのは，過去1年分のデータだけなので，まずはそれを使ってこの先1年分の需要予測を行うことにした．

SCMOPTには，ベイズ推論に基づく B-Forecastと深層（機械）学習に基づく
D-Forecast
があるので，まずは特定の製品と顧客に対して，それらを利用して予測を行うことにする．

``` python
random_seed = 4
cust_df = generate_cust(num_locations = 1, random_seed = random_seed, prefecture = None, no_island=True)
prod_df = generate_many_prod(1, weight_bound=(1, 5), cust_value_bound =(5,10),
                              dc_value_bound=(2, 2), plnt_value_bound=(1, 1),
                              fc_bound=(100000,200000), random_seed=random_seed)
weekly_ratio = [0.0, 1.0, 1.2, 1.3, 0.9, 1.5, 0.2]  # 0 means Sunday
yearly_ratio = [1.0 + np.sin(i)*0.5 for i in range(13)]  # 0 means January
demand_df = generate_demand(cust_df, prod_df, cust_shape=0.9, prod_shape=0.9, weekly_ratio=weekly_ratio, yearly_ratio=yearly_ratio,
                            start="2019/01/01", periods=365, epsilon=1., random_seed=random_seed)
demand_df.tail()
```

### ベイズ推論で予測

ベイズ推論を使うと，傾向変動や季節変動を自分でコントロールすることができる．
まずは，今回使うデータは，日ベースのデータなので，週次と年次の2つの季節変動を加え，傾向変動を（区分的）線形関数と仮定して予測をしてみると，以下のような結果が得られる．

``` python
cust_selected = [ demand_df["cust"].iloc[0] ]
prod_selected = [ demand_df["prod"].iloc[0] ]
print(cust_selected, prod_selected)
dic = forecast_demand(demand_df, cust_selected, prod_selected, agg_period ="1d", forecast_periods =365, cumulative=False, lb=0., ub=None, growth ="linear",
                     daily_seasonality=False, weekly_seasonality=True, yearly_seasonality=True)
print(cust_selected[0], prod_selected[0])
if dic[ cust_selected[0], prod_selected[0] ] is not None:
    model, forecast = dic[ cust_selected[0], prod_selected[0] ]
    model.plot(forecast);
```

``` python
from fbprophet.plot import plot_plotly, plot_components_plotly
import plotly.offline as py

fig = plot_plotly(model, forecast)
plotly.offline.plot(fig);

fig = plot_components_plotly(model, forecast) 
plotly.offline.plot(fig);
```

どうやら，年度末の下降傾向が引き続き起こると考えて予測しているようだ．確認のため，予測要因ごとに図示してみると，週の波動や年次の波動はきちんと予測されているが，傾向変動だけが直線的に下降している．
需要が負になることはないので，今度は，下限を設定してから傾向をロジスティック曲線と仮定して予測してみる．

``` python
cust_selected = [ demand_df["cust"].iloc[0] ]
prod_selected = [ demand_df["prod"].iloc[0] ]
print(cust_selected, prod_selected)
dic = forecast_demand(demand_df, cust_selected, prod_selected, agg_period ="1d", forecast_periods =365, cumulative=False, lb=0., ub=None, growth ="logistic",
                     daily_seasonality=False, weekly_seasonality=True, yearly_seasonality=True)
print(cust_selected[0], prod_selected[0])
if dic[ cust_selected[0], prod_selected[0] ] is not None:
    model, forecast = dic[ cust_selected[0], prod_selected[0] ]
    model.plot(forecast);
```

予測量が負になることは避けられたが，やはり下降傾向と予測されている．
営業担当部長に聞くと，この顧客は定番製品に対しては，毎年同じように発注をしてくるらしい．
そこで，傾向変動なしで予測してみる．

``` python
cust_selected = [ demand_df["cust"].iloc[0] ]
prod_selected = [ demand_df["prod"].iloc[0] ]
print(cust_selected, prod_selected)
dic = forecast_demand(demand_df, cust_selected, prod_selected, agg_period ="1d", forecast_periods =365, cumulative=False, lb=0., ub=None, growth ="flat",
                     daily_seasonality=False, weekly_seasonality=True, yearly_seasonality=True)
print(cust_selected[0], prod_selected[0])
if dic[ cust_selected[0], prod_selected[0] ] is not None:
    model, forecast = dic[ cust_selected[0], prod_selected[0] ]
    model.plot(forecast);
```

``` python
fig = plot_plotly(model, forecast)
plotly.offline.plot(fig);
```

### １年分のデータで深層学習と機械学習

``` python
agg_period ="1d"
df, future = make_forecast_df(demand_df, customer=demand_df["cust"].iloc[0], product=demand_df["prod"].iloc[0], 
                              promo_df= None, agg_period="1d", forecast_periods =360)
```

``` python
learn = make_learn(df, horizon="2019-10-01", classification=False)
lr_min, lr_steep = learn.lr_find() 
print(lr_min, lr_steep)
```

``` python
train(learn, n_epoch = 100, max_lr=0.01, wd=2.)
```

``` python
fig, yhat, error = predict_using_dl(learn, df, future, classification=False)
print("error=", error)
plotly.offline.plot(fig);
```

``` python
feature = None
forest, importance = random_forest(df, horizon="2019-10-01", feature=feature, classification=False)
imp_df, fig = show_importance(importance, feature=feature, df=df)
plotly.offline.plot(fig);
fig, yhat_rf, error_rf =  predict_using_rf(df, future, forest, horizon="2019-10-01", feature=feature)
plotly.offline.plot(fig);
print("error", error_rf)
```

結果を営業部長に見せたところ，この顧客のこの製品に対しては，こんな感じで十分だろうとの評価をもらったが，そもそも毎年同じ需要なら，予測を行う必要はないだろう．
深層学習と機械学習も誤差が大きく，年次の季節変動は重要でないと判断されている．さらには，新製品などは，1年分もデータがとれないものがあるとの指摘を受けた．
1年分もデータがない場合には，年次の季節変動は役に立たない．
さて，どうすれば良いだろうか？

ひとつの方法は，できるだけ過去に遡ったデータを集めることだ．少なくとも2年分のデータがなければ，年次の季節変動は予測できない．新製品に対しては，過去のデータがないので，類似の製品の過去のデータを使うことが考えられる．
そのような場合に使える手法として，深層学習がある．深層学習では，時系列データをそのまま扱うのではなく，曜日，月末，年末などから特徴を自動的に抽出することによって予測を行う．

深層（機械）学習でも，顧客と製品を1つに絞って予測することもできる．2年分のデータがあると仮定して，ベイズ推論と比較してみよう．

### ベイズ推論で予測

結果は１年分とほぼ同じだが，線形傾向でも予想は負にならない．2020年10月以降のデータで交差検証を行うと，2乗平均誤差は163程度になった．

``` python
random_seed = 4
cust_df = generate_cust(num_locations = 1, random_seed = random_seed, prefecture = None, no_island=True)
prod_df = generate_many_prod(1, weight_bound=(1, 5), cust_value_bound =(5,10),
                              dc_value_bound=(2, 2), plnt_value_bound=(1, 1),
                              fc_bound=(100000,200000), random_seed=random_seed)
weekly_ratio = [0.0, 1.0, 1.2, 1.3, 0.9, 1.5, 0.2]  # 0 means Sunday
yearly_ratio = [1.0 + np.sin(i)*0.5 for i in range(13)]  # 0 means January
demand_df = generate_demand(cust_df, prod_df, cust_shape=0.9, prod_shape=0.9, weekly_ratio=weekly_ratio, yearly_ratio=yearly_ratio,
                            start="2019/01/01", periods=365*2, epsilon=1., random_seed=random_seed)
demand_df.tail()

cust_selected = [ demand_df["cust"].iloc[0] ]
prod_selected = [ demand_df["prod"].iloc[0] ]
print(cust_selected, prod_selected)
dic = forecast_demand(demand_df, cust_selected, prod_selected, agg_period ="1d", forecast_periods =365, cumulative=False, lb=0., ub=None, growth ="linear",
                     daily_seasonality=False, weekly_seasonality=True, yearly_seasonality=True)
print(cust_selected[0], prod_selected[0])
if dic[ cust_selected[0], prod_selected[0] ] is not None:
    model, forecast = dic[ cust_selected[0], prod_selected[0] ]
    model.plot(forecast);
```

交差検証を行い，誤差を評価しておこう．

``` python
from fbprophet.diagnostics import cross_validation
from fbprophet.plot import plot_cross_validation_metric
from fbprophet.diagnostics import performance_metrics
df_cv = cross_validation(model, initial='612 days', period='90 days', horizon = '90 days')
#df_cv.head()
```

``` python
df_p = performance_metrics(df_cv, rolling_window=0.1)
df_p.mse.mean()
```

### １つの顧客と製品に対して深層学習で予測

埋め込み層の後に，200層,
100層の全結合層を入れたニューラルネットで学習させた．2020年10月以前が訓練データ，それ以降が検証データである．最大学習率0.01，重み減衰2のone
cycle法で100エポック訓練したところ，2乗平均誤差は25程度になった．

``` python
agg_period ="1d"
df, future = make_forecast_df(demand_df, customer=demand_df["cust"].iloc[0], product=demand_df["prod"].iloc[0], 
                              promo_df= None, agg_period="1d", forecast_periods =360)
```

``` python
learn = make_learn(df, horizon="2020-10-01", classification=False)
lr_min, lr_steep = learn.lr_find() 
print(lr_min, lr_steep)
```

``` python
train(learn, n_epoch = 100, max_lr=0.01, wd=2.)
```

``` python
fig, yhat, error = predict_using_dl(learn, df, future, classification=False)
print("error=", error)
plotly.offline.plot(fig);
```

### 機械学習による予測

同様のデータに対して，ランダム森による学習を行ったところ，2乗平均誤差は50程度になった．特徴の重要性を可視化すると以下のようになり，特徴は順に，曜日，年初からの日にち，月，週が効いていることが分かる．

``` python
feature = None
forest, importance = random_forest(df, horizon="2020-10-01", feature=feature, classification=False)
imp_df, fig = show_importance(importance, feature=feature, df=df)
plotly.offline.plot(fig);
fig, yhat_rf, error_rf =  predict_using_rf(df, future, forest, horizon="2020-10-01", feature=feature)
#plotly.offline.plot(fig);
print("error", error_rf)
```

### 全ての製品を一度に予測

深層学習では，顧客名や製品名もカテゴリーデータとして用いることができる．これらの特徴を「埋め込み」を用いて適当な大きさのベクトルとして表現して，予測精度の向上を図るのである．ある製品の需要だけが1年分のデータしかないものとしよう．他の製品と同時に深層学習を行うことによって，どれだけ予測の精度が向上するか確認してみる．

上の深層学習と同じ程度の訓練を行うと，100製品に対する合計の2乗平均誤差が362となり，以下のような予測が得られた．この手法を用いると，販売実績の少ない新製品に対する予測も可能になる．

``` python
random_seed = 4
cust_df = generate_cust(num_locations = 1, random_seed = random_seed, prefecture = None, no_island=True)
prod_df = generate_many_prod(100, weight_bound=(1, 5), cust_value_bound =(5,10),
                              dc_value_bound=(2, 2), plnt_value_bound=(1, 1),
                              fc_bound=(100000,200000), random_seed=random_seed)
weekly_ratio = [0.0, 1.0, 1.2, 1.3, 0.9, 1.5, 0.2]  # 0 means Sunday
yearly_ratio = [1.0 + np.sin(i)*0.5 for i in range(13)]  # 0 means January
demand_df = generate_demand(cust_df, prod_df, cust_shape=0.9, prod_shape=0.9, weekly_ratio=weekly_ratio, yearly_ratio=yearly_ratio,
                            start="2019/01/01", periods=365*2, epsilon=1., random_seed=random_seed)
customer=demand_df["cust"].iloc[0]
product=demand_df["prod"].iloc[0]
demand_df["date"] = pd.to_datetime(demand_df["date"])
```

``` python
demand_df["delete"] = (demand_df["prod"]==product) & (demand_df.date<=pd.to_datetime("2019-12-31"))
```

``` python
demand_df = demand_df[  demand_df.delete==False ]
demand_df.drop("delete", axis=1, inplace=True)
```

``` python
all_df, future = make_forecast_df_all(demand_df, promo_df=None, agg_period="1d",forecast_periods=365)
```

``` python
future.head()
```

``` python
learn = make_learn(all_df, horizon="2020-10-01", classification=False)
lr_min, lr_steep = learn.lr_find() 
print(lr_min, lr_steep)
```

``` python
train(learn, n_epoch = 50, max_lr=0.001, wd=2.)
```

``` python
fig, predict, error = predict_using_all_dl(learn, df, future, customer=customer, product = product, classification=False)
plotly.offline.plot(fig);
```

## 生産計画

ロジスティックネットワーク設計の結果を，主力の工場長に知らせたところ，繁忙期にはこんなに生産できないとの返事が返ってきた．
ロジスティックネットワーク設計モデルでは，年間が需要量が工場の生産能力以下であるという条件で解いていたが，
我社の製品は売れる月とそうでない月があるようだ．
さて，どのようなアプローチをとれば良いのであろうか？

まずは，今まで年ベースで考えていたデータを，工場の生産計画単位である月ベースにすることから始めることにする．

``` python
# データの読み込み
#demand_df = pd.read_csv(folder+"case-demand.csv", compression='gzip')
prod_df = pd.read_csv(folder+"case-prod.csv")
cust_df = pd.read_csv(folder+"case-cust.csv")
agg_df_prod = pd.read_csv(folder+"case-agg-df-prod.csv")
flow_df = pd.read_csv(folder+"case-flow.csv", index_col=0)
```

``` python
aggregated_cust_df = pd.read_csv(folder+"case-agg-cust.csv")
trans_df  = pd.read_csv(folder+"case-trans.csv")
plnt_prod_df  = pd.read_csv(folder+"case-plnt-prod.csv")
```

ロジスティックネットワーク設計モデルは，顧客の需要を集約して（クラスタリングして）扱っていた．元データの需要は，各顧客・製品に対する日単位での需要量を保持していたが，これをクラスター単位に集約する．
さらに，ロジスティックネットワーク設計で求めた最適なフロー（工場・倉庫，倉庫・顧客間の製品の移動）から，工場（ならびに倉庫）単位での製品の需要量を計算する．

``` python
aggregated_cluster_df = aggregate_demand_by_cluster(cust_df, aggregated_cust_df, demand_df)
```

``` python
#クラスターの需要データから，倉庫ならびに工場での需要量を計算する．
aggregated_demand_df, G, Prod = aggregate_demand(aggregated_cluster_df, flow_df)
```

``` python
#aggregated_demand_df.set_index("date")
aggregated_demand_df[["cust","prod","demand"]].head()
```

我社の工場は３つあり，
1つは主力工場ですべての製品を作ることができ，他の2つは，約半分の製品だけを作ることが可能である．
そのため，製品ごとのサプライ・チェインを描画すると，以下のような図になる．
左から工場(2箇所），倉庫(13箇所），顧客(50箇所）であり，2つの工場で役割分担をしていることが確認できる．

``` python
p ="prod44"
pos = G[p].layout()
nx.draw(G[p], pos=pos,  node_color="Blue", node_size=10)
```

``` python
#安全在庫量の簡易計算
safety_df = compute_safety_stock(aggregated_cluster_df, aggregated_demand_df, trans_df, plnt_prod_df, Prod, G)
```

``` python
safety_df.head()
```

さて，主力の工場に割り当てられている倉庫（顧客）の需要量を計算してみる．製品は全部で24種類あり，すべてABCのAカテゴリーなので比較的大きな量になっている．
月別に集約してから製品別に可視化してみよう．

``` python
#生産計画（工場ごと）: 小田原工場のみ
demand_odawara = aggregated_demand_df[ aggregated_demand_df.cust=="Plnt_Odawara" ].copy()
pd.pivot_table(demand_odawara, index="cust", columns="prod", values="demand", aggfunc=sum)
```

``` python
try:
    demand_odawara.reset_index(inplace=True)
except:
    pass
demand_odawara["date"] = pd.to_datetime(demand_odawara["date"])
demand_odawara.set_index("date",inplace=True)
```

``` python
monthly_demand = demand_odawara.groupby(["prod"]).resample("MS").sum()
monthly_demand.head()
```

``` python
try:
    monthly_demand.reset_index(inplace=True)
except:
    pass
fig = px.line(monthly_demand , x="date", y="demand", color="prod")
#plotly.offline.plot(fig);
Image("../figure/case-monthly-demand.png")
```

``` python
# #製品ごと（本来ならば工場ごと）に集約
demand_ = pd.pivot_table(monthly_demand, index= "prod", columns ="date",  values="demand", aggfunc=sum)
demand = demand_.values
_, T = demand.shape 
print("T=",T)
```

``` python
demand_.head()
```

``` python
#生産工程の容量を計算
average_demand = demand.sum()/T
average_demand
```

``` python
prod_df.head()
```

``` python
cycle_time= prod_df["lot_size"].mean()/prod_df["average_demand"].mean()
cycle_time #CT = 10 weeks
```

``` python
(0.1*average_demand/3600 + 24/2.5 )/(24*30.5)
```

``` python
capacity = 24*30.5
prod_time_bound =(1/36000, 1/36000)
setup_time_bound = (1,1)
prod_cost_bound = (1/2000., 1/1000.)
print("capacity=",capacity)
#資源データフレーム
#期ごとに容量が変化するモデル
#print(capacity)
name_, period_, capacity_ = [], [], []
s = 1
for t in range(T):
    name_.append( f"Res{s}")
    period_.append(t)
    capacity_.append( capacity)
resource_df = pd.DataFrame(data={"name": name_, "period": period_, "capacity": capacity_})

# #生産情報データフレーム生成
items = list(prod_df.name)
num_item = len(items)
production_df = pd.DataFrame(data={"name": items,
                             "ProdTime": [random.uniform(prod_time_bound[0], prod_time_bound[1]) for i in range(num_item)],
                             "SetupTime": [random.randint(setup_time_bound[0], setup_time_bound[1]) for i in range(num_item)],
                             "ProdCost": [random.uniform(prod_cost_bound[0], prod_cost_bound[1]) for i in range(num_item)],
                             "SetupCost": list( prod_df.fixed_cost ) 
                             })
production_df.set_index("name", inplace=True)
production_df.head()
prod_df.to_csv(folder+"case-prod.csv")
production_df.to_csv(folder+"case-production.csv")
# bom_df.to_csv(folder+"bom.csv")
resource_df.to_csv(folder+"case-resource.csv")
```

``` python
production_df.head()
```

生産時間，段取り時間，生産（変動）費用，生産固定（段取り）費用を入れ，ロットサイズ最適化システム
OptLotで求解してみる．
結果から，在庫をうまく使えば生産容量を超えることなしに，すべての需要を満たすことができることが分かる．

``` python
model, T = lotsizing(prod_df, production_df, bom_df=None, demand =  demand, resource_df=resource_df, max_cpu= 10)
```

``` python
print("model status",model.Status)
x, I, y, slack, surplus = model.__data
violated, production, inventory, fig_inv, fig_capacity = show_result_for_lotsizing(model, T, prod_df, production_df, bom_df=None, resource_df=resource_df)
#plotly.offline.plot(fig_capacity);
#plotly.offline.plot(fig_inv);
```

### ローリングホライズン方式

上の結果を工場長にプレゼンしたところ，年末にかけて生産量が減少し，在庫も減っていることが指摘された．これは，年度末の在庫量を，工場で指定された安全在庫量に設定したからであるが，実際には少なすぎるようだ．さて，どのようにこの問題点を乗り切れば良いだろうか？

これは，無限期間の動的な問題を，12ヶ月の静的な問題に帰着させたことによる．実際の会社は，12月で閉める訳ではない．解決法としては色々考えられるが，一番簡単な方法は，より簡単には，年末の在庫量を安全在庫量ではなく，目標在庫量に変更する方法である．しかし，この方法だと目標在庫量の設定に曖昧さが残る．より実用的な方法としてローリングホライズン方式がある．この方法は，12ヶ月の静的な問題を，毎月時間をずらしながら解いていくものである．毎月１２ヶ月分の問題を解き直すので，動的な環境の変化にも追従できる．

## スケジューリング最適化

工場長の承認も得れたので現場の方に生産計画を提出したところ，生産資源をきちんと考えないと，生産可能かどうか判断できないとの指摘を受けた．
さて，どうすれば良いだろうか？

月次の生産最適化モデルは，タクティカルレベルの意思決定であり，より細かい生産計画はオペレーショナルレベルのスケジューリング最適化で行う必要がある．
ここでは，SCMOPTに含まれるスケジューリング最適化システム OptSeq
を使って問題を解決することにしよう．

OptSeqは様々な付加条件に対応できるが，ここでは簡単のため単一の資源制約のもとで，スケジューリング最適化をしてみる．結果は，以下のような詳細なガントチャートが得られる．もちろん，上述のローリングホライズン方式を用いて，適当な期間ごとに再最適化を行い，スケジュールを見直す必要がある．

``` python
activity =[ ]
prod = list(production.index)
T = len(production.columns)
duration = {}
#print(T)
for t in range(T):
    for i, vol in enumerate(production[t]):
        if vol > 0.001:
            duration[i,t] = round(vol * 60)
```

``` python
model = optseq.Model()
act, mode, resource = {}, {}, {} 
resource["machine"] = model.addResource("machine",capacity = 1)
for (i,t) in duration:
    act[i,t] = model.addActivity(f"act[{i},{t}]", duedate = 30*24*(t+1)*60)
    mode[i,t] = optseq.Mode(f"mode{i},{t}", duration= duration[i,t])
    mode[i,t].addResource(resource["machine"], requirement=1)
    act[i,t].addModes(mode[i,t])
```

``` python
model.Params.TimeLimit = 10
model.optimize()
```

``` python
fig = optseq.make_gannt(model, start = "2019/01/01", period ="minutes")
fig.update_layout(
        autosize=False,
        width=6000,
        height=2000,
        paper_bgcolor="LightSteelBlue",
    )
plotly.offline.plot(fig);
```

## 安全在庫配置

社内で新しい倉庫の立地について検討したところ，営業部と倉庫の担当者から在庫が考慮されていないのではないかと問われた．
さて，あなたは次に何をすれば良いのだろうか？

SCMOPTの基本分析システム
SCBASには，基本的な在庫解析が含まれているので，まずはそれを使ってみよう．

リスク共同管理分析は，在庫をサプライ・チェインの上流（工場側）に置くべきか，下流（顧客側）に置くべきかの簡易判断をしてくれる．この会社の場合は（需要をランダムに同じ分布で生成したので当然だが），需要の相関が強いため，上流に置いてもあまり在庫の削減はできないことが分かる．

需要が無相関だと，顧客数を n としたとき $1/\sqrt{n}$
になることが知られているが，この場合はほとんど削減されない．一般に，同じ時期に同じようなキャンペーンを行う場合には，相関が高いと考えられる．

``` python
inv_reduction_df = risk_pooling_analysis(demand_df, agg_period="1w")
```

``` python
inv_reduction_df.sort_values("reduction_ratio", ascending=False).head()
inv_reduction_df.to_csv(folder+"case-inv-reduction.csv")
```

``` python
fig = show_inventory_reduction(inv_reduction_df)
plotly.offline.plot(fig);
```

次にSCBASの在庫解析関数で，（工場が1つと仮定した場合の）工場における生産ロットサイズや安全在庫量を計算してみよう．これは，古典的な経済発注量モデルと新聞売り子モデルを用いたものである．

以下のように，製品別の需要量，標準偏差，在庫費用，ロットサイズ，安全在庫量などが計算される．

``` python
prod_df = inventory_analysis(prod_df, demand_df, inv_reduction_df, z = 1.65, LT = 1, r = 0.3, num_days=7)
```

``` python
prod_df.head()
```

``` python
prod_df.to_csv(folder+"case-prod-all.csv")
```

``` python
# 分布型を確認
d = aggregated_cluster_df[ (aggregated_cluster_df["prod"]=="prod1") & (aggregated_cluster_df.cust=="cust0") ]
d.hist(bins=100);
```

``` python
# 分布型を確認（工場）
d = aggregated_demand_df[ (aggregated_demand_df["prod"]=="prod1") & (aggregated_demand_df.cust=="Plnt_Odawara") ]
d.hist(bins=100);
```

次に，安全在庫配置システム MESSA (MEta Safety Stock Allocation
system)を利用して，ネットワーク全体での在庫の配置を最適化しよう．

準備として，ロジスティックネットワーク設計モデルで得られた最適フローをもとに，集約された顧客（クラスター），倉庫，工場ごとの（1日あたりの）製品の需要の平均と標準偏差を計算する．
これが，安全在庫配置を最適化するための基本データとなる．
このモデルで決めるのは，在庫保管地点ごとの保証リード時間（サービス時間）である．

このモデルでは，製品ごとに在庫補充の流れを表すネットワークを構築する．
今回の例では，各製品の在庫は必ず1つの地点から補充されているので，地点における作業時間にリード時間を加えておく．（この仮定は必ずしも成立するとは限らない．ロジスティックネットワーク設計モデルで単一ソース制約を付加しても，工場の容量制約があれば，倉庫は2つの工場から補充を受ける可能性がある．
その場合には，工場と倉庫間に輸送を表す点を追加し，その作業時間を工場・倉庫間のリード時間と設定する．）

他にも，在庫費用や品切れ費用を準備する必要がある．在庫費用は評価しにくいデータであるが，ここでは製品の価値と在庫保管率を元に計算することにする．
製品の価値は地点別に計算されていると仮定する．ここでは，工場，倉庫，顧客における価値が徐々に増加するように設定しておく．

品切れ（バックオーダー）費用は，製品の価値と同じと考える．つまり，品切れしている場合は製品の価値と同等額が失われたものと考える．
在庫費用は，製品の価値に在庫保管比率を乗じたものとする．
ここで保管費率は以下の量の和である．

- 利子率:投資額利率
- 保険料率: 製品の種類および企業の方針によって異なる．
- 消耗費率および陳腐化率: 製品の腐敗，破損，目減りなどを考慮して計算
- 税率: 在庫に課せられる法的な税率（日本では0)

サービス率は品切れを起こさない確率であり，古典的な新聞売り子モデルから，品切れ費用と在庫費用から（つまり保管比率から）計算できる．
安全在庫係数は，需要の分布が決まればサービス率から計算できる．ここでは．需要は正規分布と仮定して計算する．
（顧客の需要は正規分布とは限らないが，クラスターに集約された顧客や倉庫，工場に集約された需要は正規分布に近いと考えられる．）

``` python
#顧客（クラスター）に対する需要の平均と標準偏差を計算
p = Prod[0]
cluster_dem = aggregated_cluster_df[ aggregated_cluster_df["prod"] == p ]
cluster_sum_std_df = pd.pivot_table(cluster_dem, index= "cust", values="demand", aggfunc=[np.average, np.std])
cluster_sum_std_df.head()
```

``` python
cluster_sum_std_df.reset_index(level=0,inplace=True)
cluster_sum_std_df["cust"] = "Cust_" + cluster_sum_std_df["cust"]
cluster_sum_std_df.head()
```

``` python
# 工場・倉庫，製品ごとの合計需要量と標準偏差を計算
p = Prod[0]
#print(p)
agg_dem = aggregated_demand_df[ aggregated_demand_df["prod"] == p ]
deamnd_sum_std_df = pd.pivot_table( agg_dem , index= "cust", values="demand", aggfunc=[np.average, np.std])
deamnd_sum_std_df.reset_index(level=0,inplace=True)
```

``` python
deamnd_sum_std_df.head()
```

``` python
# 合体
stage_df = pd.concat( [ cluster_sum_std_df, deamnd_sum_std_df] )
stage_df.columns =["name", "average_demand", "sigma"]
stage_df["cv"] = stage_df.sigma/ stage_df.average_demand #変動係数を追加
#stage_df
```

``` python
#点の処理時間を保持する辞書を作成
#単一ソースを仮定するので，リード時間も作業時間に加える
stage_time_dic ={}
piv = pd.pivot_table(trans_df, index=["to_node", "kind"], values=["lead_time","stage_time"])
piv.reset_index(inplace=True)
for row in  piv.itertuples():
    if row.kind == "plnt-dc":
        stage_time_dic["DC_"+ str(row.to_node)] = int(row.stage_time) + int(row.lead_time)
    else:
         stage_time_dic["Cust_"+str(row.to_node) ] = int(row.stage_time) + int(row.lead_time)
```

``` python
#工場での生産リード時間を計算，製品の情報をprod_infoに保管
try:
    prod_df.set_index("name", inplace=True)
except:
    pass
prod_info = prod_df.loc[p]
plant_LT = round(prod_info.lot_size/prod_info.average_demand)
print(prod_info)
```

``` python
#各地点の作業時間，在庫費用，品切れ費用，安全在庫係数を，保管比率と価値から計算
from scipy.stats import norm

try:
    stage_df.set_index("name", inplace=True)
except:
    pass

r = 0.1 #保管比率は与える
stage_time_list, z_list, h_list, b_list = [], [] ,[], []
for row in stage_df.itertuples():
    if row.Index[:4]=="Plnt":
        stage_time_list.append(plant_LT)
    else:
        stage_time_list.append( stage_time_dic[row.Index] )
        
    if row.Index[:4] == "Plnt":
        b = prod_info.plnt_value
    elif row.Index[:2] == "DC":
        b = prod_info.dc_value #=1
    else:
        b = prod_info.cust_value #価値 = 品切れ費用
    h = r*b     #在庫費用
    p_ = b/(b+h) #サービス率
    z = norm.ppf(p_) #安全在庫係数
    z_list.append(z)
    b_list.append(b)
    h_list.append(h)

stage_df["processing_time"] = stage_time_list
stage_df["z"] = z_list
stage_df["h"] = h_list
stage_df["b"] = b_list
stage_df.reset_index(inplace=True)
```

``` python
# bom_dfを構築
bom_df = flow_df[ flow_df["prod"] ==p ].copy()
bom_df.columns =["child", "parent", "prod", "flow"]
bom_df["units"] = 1
total_flow_from = defaultdict(float) #工場もしくは倉庫から出るフローの合計を保管する辞書
for row in bom_df.itertuples():
    total_flow_from[ row.child ] += row.flow 
#配分率を計算
allocation_list =[]
for row in bom_df.itertuples():
    if total_flow_from[ row.child ] > 0.001:
        allocation_list.append( row.flow/total_flow_from[ row.child ] )
    else:
        allocation_list.append(1.)
bom_df["allocation"] = allocation_list
```

在庫保管地点のデータと，地点間の情報を保管した部品展開表（BOM）データを元に，安全在庫配置最適化を行う．

``` python
# 他のデータも準備して安全在庫配置問題を求解
stage_df['max_guaranteed_LT'] = 0
pos = G[p].layout()
x_list, y_list =[],[]
for row in stage_df.itertuples():
    x_list.append( pos[row.name][0] )
    y_list.append( pos[row.name][1] )
stage_df["x"] = x_list
stage_df["y"] = y_list
best_cost, stage_df, bom_df, fig = solve_SSA(stage_df, bom_df)
```

工場と倉庫での製品の価値が同等と仮定した場合，以下のような結果が得られる．
点の大きさが正味補充時間（安全在庫を保持する日数）を，点の色の濃さが保証リード時間を表している
この場合は，工場には在庫を置かず倉庫にまとめるのが最適であることが分かる．工場の保証リード時間は長く，倉庫が発注をしてから生産を開始するため在庫はなしで運用する．
そのかわり，倉庫には20日分程度の在庫を保持する．

``` python
stage_df
```

倉庫における製品の価値が工場での価値の2倍と仮定した場合は，以下のような結果が得られる．

この場合は，すべての地点の保証リード時間が0となり，すべての地点で自分の作業時間分の安全在庫をもつのが最適になる．

安全在庫配置システムMESSAには，シミュレーションをしながら基在庫レベルを最適化するための方法も入っている．
この手法は，各地点での1日あたりの処理量上限の制約も付加できるが，ここでは十分な処理能力があると仮定して，最適化を行う．

この手法では在庫がない場合に，下流の地点の注文をどのように配分するかを定義する必要がある．ここでは，ロジスティックネットワーク設計での製品のフロー量をもとに（実績量をもとに）配分を行うものとする．

``` python
# エシェロン基在庫レベルを計算
stage_df["capacity"] = 999999999999.
cost, stage_df, I = periodic_inv_opt(stage_df, bom_df, max_iter = 10, n_samples = 10, n_periods = 100, seed = 1)
```

各在庫保管地点での在庫の変化の図と最適な基在庫レベルが得られる．1つの顧客と倉庫ならびに工場でのシミュレーション結果の一部を表示する．

基在庫レベルを用いた在庫の運用は簡単であり，各地点の在庫ポジション（在庫量にすでに発注済みで未到着の量を加え，バックオーダー量を減じた量）が基在庫レベルになるように発注するだけである．

## 配送最適化

ロジスティックネットワーク設計の結果を上司にプレゼンしたところ，最適化された倉庫から実際に配送可能なのか質問された．さらに取締役へ報告したところ，配送費用が正しく評価されているか調べるように命令された．さて，あなたは次に何をすれば良いのだろうか？

SCMOPTには配送最適化ソルバー METRO
IVが入っている．ロジスティックネットワーク設計モデルで得られた倉庫をデポと考え，配送最適化モデルを解くことによって，実際に配送が可能かを検証し，実際の配送費用を計算してみる．

まず，ネットワーク設計開設すると決められた倉庫と，それらの倉庫に割り当てられた顧客の情報をもとにして，顧客のクラスタリングを行う．
各倉庫をデポとした独立な配送最適化問題を作成して，配送最適化ソルバー
METRO IVで求解する．

顧客の需要量は，年間稼働日数を200日と仮定して日割りにし，
各デポ（倉庫）には，20トン車(容量は20000kg）のトラックを，必要台数分だけ配置する．1つの顧客の総需要量がトラックの容量を超過した場合には，需要を分割して複数の顧客を設定する．

1日の稼働時間は10時間とし，0から36000秒の時間枠を設定する．各顧客上での作業時間は0とし，荷物は配達のみとする．
地点間の移動時間と距離は，OSRMを使って計算する．また，ルートが見つからない地点間に対しては，大圏距離の5倍の距離を設定し，そこを時速50kmで移動すると仮定して移動時間を計算する．

目的関数は稼働時間の合計を最小化するものとし，配送しきれない顧客がいた場合にはペナルティを与えるものとする．

配送計画のデータを生成する．

``` python
#list(map( len, customer_df_dic.values())) #number of customers
```

配送最適化モデルをデポごとに解く．

``` python
shipment_df = break_df = pd.DataFrame()

total_cost = total_distance = total_service = 0
capacity = 20000
vehicle_dic, job_dic = {}, {}
for instance, depot_name in enumerate(depots):
    print(i,depot_name)
    #depot_name =depots[10]
    job_df = customer_df_dic[ depot_name ].copy()
    job_df["location"] = "[" + job_df.lon.astype(str)  + ","+ job_df.lat.astype(str) +"]"
    job_df["location_index"] = None 
    job_df["daily_demand"] = job_df.demand/200*3 #年間稼働日は200日と仮定；製品の平均重量は 3と仮定
    job_df["daily_demand"] = job_df["daily_demand"].astype(int)
    job_df["pickup"] = "[0]"
    job_df["delivery"] = "[" + job_df.daily_demand.astype(str) +"]"
    job_df["time_windows"] = "[[0,36000]]" #10 hours
    job_df["service"] = 0  #10 minutes
    job_df["skills"] = "[0]"
    job_df["priority"] = 100

    #容量オーバーの需要を分割
    over_jobs = job_df[ job_df.daily_demand > capacity ]
    job_df.drop( over_jobs.index, inplace=True)
    for i in range(len(over_jobs)):
        row = over_jobs.iloc[i]
        dem = row.daily_demand
        while dem > capacity:
            new_row = row.copy()
            new_row.daily_demand = capacity
            dem -= capacity
            job_df = job_df.append(new_row)
        new_row = row.copy()
        new_row.daily_demand = dem
        job_df = job_df.append(new_row)
    job_df["daily_demand"] = job_df["daily_demand"].astype(int)
    job_df["delivery"] = "[" + job_df.daily_demand.astype(str) +"]"

    try:
        dc_df.set_index("name", inplace=True)
    except:
        pass
    depot_info = dc_df.loc[depot_name[3:],:]
    num_vehicles = round(job_df.daily_demand.sum()/capacity)
    name_list = []
    for i in range(num_vehicles):
        name_list.append( f"vehicle{i}" )
    vehicle_df = pd.DataFrame({"name": name_list})
    vehicle_df["start"] = f"[{depot_info.lon},{depot_info.lat}]"
    vehicle_df["end"] = f"[{depot_info.lon},{depot_info.lat}]"
    vehicle_df["start_index"] = None
    vehicle_df["end_index"] = None
    vehicle_df["skills"] = "[0]"
    vehicle_df["breaks"] = "[]"
    vehicle_df["capacity"] = "[20000]"
    vehicle_df["time_window"] = "[0,36000]" 

    #辞書に保管
    job_dic[instance] = job_df 
    vehicle_dic[instance] = vehicle_df

    matrix = False #時間行列を用いる場合 True
    map_style = 0
    #shipment_df = break_df = pd.DataFrame()
    model = build_model_for_vrp(job_df, shipment_df=None, vehicle_df=vehicle_df, break_df=None)

    #最適化
    input_dic, output_dic, error = optimize_vrp(model,matrix=matrix)

    if error != "":
        print(error) #Streamlitでエラーを表示
        loc_error = "b'[Error] Unfound route(s) from location"
        if str(error)[:40] == loc_error:
            #location error
            #prepare kd-tree
            from scipy.spatial import KDTree  
            xy=[]
            for i in job_df.location: #本来ならnode_dfで保持する．
                lon, lat =ast.literal_eval(i)
                xy.append( (lon,lat) )
            tree = KDTree(xy)

            lon, lat = str(error)[41:].split(";")
            point =float(lon[1:]), float(lat[:-4])
            dis, idx = tree.query(point)
            print( f"Undound route at point {idx}" ) # + str(job_df.iloc[idx]) )
            job_df.drop(job_df.index[idx], inplace=True)  # delete the point
    else:
        if matrix:
            fig = make_fig_for_vrp(input_dic, output_dic, show_mode="route", map_style=map_style)
        else:
            fig = make_fig_for_vrp(input_dic, output_dic, show_mode="road", map_style=map_style)
        #plotly.offline.plot(fig);
        #fig.write_html(f"map{instance}.html", include_plotlyjs=True)
        #route_df_dic, unassigned_job_df, unassigned_shipment_df = make_route_for_vrp(job_df, shipment_df, vehicle_df, break_df, output_dic)

        #fig = gannt_for_vrp(route_df_dic, start= "2020/01/01 8:00" )
        #plotly.offline.plot(fig);
        #fig.write_html(f"gannt{instance}.html")

    #     fig = make_tw_fig_for_vrp(job_df, shipment_df, start= "2020/01/01 8:00")
    #     plotly.offline.plot(fig);
        cost = output_dic["summary"]["cost"]
        service = output_dic["summary"]["service"]
        distance = output_dic["summary"]["distance"]
        print("cost=", cost)
        total_cost += int(cost)
        total_service += int(service)
        total_distance += distance 
    #break 
print("total_cost, service, distance =", total_cost, total_service, total_distance)
```

### 全体最適化

すべてのデポの問題を一度に最適化することもできる．計算時間は10分程度かかり，結果は総移動時間は7017639秒，総移動距離は139976582kmとなり，ほんの少しだが悪化する．
これは，ジョブ数が1385の大規模な問題例であるため，複数デポから配送可能な柔軟性よりも，求解が不十分であったことによる解の悪化が上回ったと考えられる．運用上も日々割り当てデポを変えることは煩雑なので，デポごとに最適化することが推奨される．

結果の図から，北海道に2箇所のデポがあり，統合可能であると社長から指摘された．
これは，ロジスティックネットワーク設計において倉庫・顧客間の最大道路距離を300kmと設定していたためであり，再計算の必要があるだろう．沖縄デポは数件しか配達していないので，
自社で倉庫をもたずに，3PLに委託することも検討すべきだ．

実際には，工場から倉庫までの輸送費用も正確に計算し，ロジスティックネットワーク設計を再び解く必要があるだろう．また，倉庫の候補地点も，顧客上に配置するのではなく，より適切な候補地点を探す必要がある．このように，複数の最適化モデルを交互に利用しながら，徐々に詳細な意思決定を行っていくのである．

``` python
((6797198)/3600 * 8000 +  135201515/1000*20 )*200 /3.379672e+09
```

``` python
((7017639)/3600 * 8000 +  139976582/1000*20 )*200 /3.379672e+09
```

``` python
((114652)/3600 * 8000 +  22109115/1000*20 )*200 /3.238943e+09
```

全部のデポに対して一度に解く．

``` python
matrix = False #時間行列を用いる場合 True
map_style = 1
#shipment_df = break_df = pd.DataFrame()
model = build_model_for_vrp(job_df, shipment_df=None, vehicle_df=vehicle_df, break_df=None)

#最適化
input_dic, output_dic, error = optimize_vrp(model,matrix=matrix,explore=5)
```

``` python
if matrix:
    fig = make_fig_for_vrp(input_dic, output_dic, show_mode="route", map_style=map_style)
else:
    fig = make_fig_for_vrp(input_dic, output_dic, show_mode="route", map_style=map_style)
plotly.offline.plot(fig);
fig.write_html(f"map_all.html", include_plotlyjs=True)

route_df_dic, unassigned_job_df, unassigned_shipment_df = make_route_for_vrp(job_df, shipment_df, vehicle_df, break_df, output_dic)

fig = gannt_for_vrp(route_df_dic, start= "2020/01/01 8:00" )
#plotly.offline.plot(fig);
fig.write_html(f"gannt{instance}.html")

fig = make_tw_fig_for_vrp(job_df, shipment_df, start= "2020/01/01 8:00")
plotly.offline.plot(fig);
cost = output_dic["summary"]["cost"]
service = output_dic["summary"]["service"]
distance = output_dic["summary"]["distance"]
print(cost, service, distance)
```

### アナリティクス

実際の配送最適化には，様々な付加条件が付く．それらに対処するためには，単に最適化を行うだけではなく，データの分析（アナリティクス）が重要になる．

たとえば，デポの位置を決めるためには，顧客の需要量の大きさや，デポからの移動時間が重要になる．通常は，単に土地が安いだけでなく，利便性を考えてインターの近所を物色することになる．その際には，デポ（の候補地点）から顧客への移動時間のヒストグラムを見ると良い．

デポに配置する運搬車の数も重要な意思決定変数である．どのような種類のトラックを何代準備しておくかは，中期の意思決定モデルであり，配送最適化は短期の意思決定モデルである．そのため，需要の不確実性だけでなく，傭車の費用なども考慮して，最適化を行う必要がある．また，時間枠も重要な要因である．(積込みと積み降ろしを考慮した）需要量と時間枠は相互に影響するので，時間帯ごとの運搬車の必要量を概算しておくことが重要になる．そのためには，以下のような時間帯別のトラックの必要台数の概算を見ると良い．

他にも実際には，時間枠のための待ち時間や休憩制約（何時から何時までの間に何分の休憩をとる必要があるという条件）を考慮する必要がある．もちろん，時間枠や休憩は複数考える必要がある．これを可視化するためには，ガントチャートを見ると良い．

## サービスネットワーク設計

注意：このケーススタディは，幾つかの企業のコンサルティングから抽出したもので，特定の企業を想定したものではない．データや登場人物は，架空のものである．

------------------------------------------------------------------------

前職でロジスティクス改革に成功したあなたは，ロジスティクス専門カンパニーにヘッドハンティングされた．この会社は，広域の3PLを行う老舗であるが，長年の感と経験に基づく意思決定に限界を感じており，最近データ収集のプロジェクトを終えたようだ．
新任早々，あなたは，収集された全拠点間データを渡され，分析するようにと命令された．さて，何から始めれば良いだろうか？

例によって，データの概要を調べることから始めよう．
与えられたデータは，以前のメーカーの例とは異なり，今度のデータは拠点間のデータのようだ．

拠点は全部で10箇所あって，主に北日本中心に営業をしている．この表で与えられているubは拠点の積替え容量，vcは積替えの変動費用を表している．また，lat,lonは緯度・経度であり，これをもとに地図上で移動時間と距離を計算する．

拠点を地図上に描画しておく．

拠点間に集約されたOD需要量のデータも与えられている．ここでODとは，Origin-DEstinationの略で，需要の発地点と着地点を表している．
このように，広域拠点間輸送最適化は，工場から顧客への一方方向のモノの流れだけではなく，多対多のモノの流れを考える必要がある．

ここでは，OD間の需要は，都道府県の人口をもとに重力法と呼ばれる手法で計算したものである．
重心法では，人口を質量と考え，質量の重い者同士は引き合う力が強く，かつ距離に反比例する形で需要が定まると仮定している．
OD需要量をヒートマップで図示すると以下のようになる．
この図から，大きな都市間の需要が大きいことが確認できる．

トラックの種類は１種類とし，容量は1000と推定されている．
地点間の移動距離と移動時間は，道路距離を地図から計算することによって算出する．
移動費用は，トラックの燃料費を除く変動費は距離に比例し，人件費と燃料費は移動時間に比例するものと仮定して，1kmあたりの費用は20円,
1時間あたりの費用は8000円と設定して計算する．

解くべき問題は，サービスネットワーク問題と呼ばれる問題であり，拠点で積替えを行い荷物を集約してトラックの積載率を向上させることによって総費用を最小化することを目的としている．
費用は，トラックの費用の他には，積替え費用がある．

この問題は数理最適化モデルとして定式化でき，市販のソルバーで解くことがでいると期待された．しかし，このタイプの問題は非常に難しく，たった10拠点の問題例でも，許容時間内に解を算出することができなかった．

そこで，拠点間のパスで使う可能性のあるものだけを生成し，近似的に求解する方法をとることにした．拠点間ごとに費用の安い方から10本のパスを生成したとき，積替え費用98853円，
トラック費用1656448円，総費用1755301円の解が得られた．

さらに，必要なパスを随時追加する勾配スケーリング法と呼ばれる専用解法を試したところ，積替え費用102847円，
トラック費用1633083円，総費用1735930円の解が得られた．
地図上に描画するとトラックの移動は複雑で，他の拠点から前橋市に向かう荷物の移動経路だけを描画すると，途中で積替えを行って荷物を集約していることが分かる．

``` python
try:
# just MIP
    #path_df, vehicle_df, cost_df, fig = solve_sndp(dc_df, od_df, cost_per_dis=20, cost_per_time=8000, capacity=1000., max_cpu = 100, scaling=False, k=10, alpha=0.9, max_iter = 10, osrm = True)
# slope scaling + column generation
    path_df, vehicle_df, cost_df, fig = solve_sndp(dc_df, od_df, cost_per_dis=20, cost_per_time=8000, capacity=1000., max_cpu = 100, scaling=True, k=5, alpha=0.5, max_iter = 10, osrm = True)
except SolverError as e:
    print(e)
```

### 配送最適化再び

サービスネットワーク設計の結果を上司にプレゼンしたところ，荷物を途中で積み合わせることにいよって費用が大幅に削減できることは理解したが，トラックの輸送経路はどうするのかと質問された．さらに取締役へ報告したところ，配送費用が正しく評価されているか調べるように命令された．

あなたは，デジャブをみているような錯覚を覚えたが，当然，配送最適化を行わなければならない．
今回は，以前の会社で行ったミルクラン形式の配送計画ではなく，長距離輸送を伴う計画を建てなければならない．
長距離輸送では，輸送業者との契約によりデポに戻ってこなくても良い場合もある．もちろん，法令に伴った休憩も考慮する必要がある．

幸い，METRO
IVは，休憩や始点（終点）を指定しないパス型の輸送も考慮できるように設計されている．
輸送業者に聞いたところ，長距離の場合には，帰りに運ぶ荷物を自分たちで見つけたいので，出発地点に戻らなくても良いということが分かった．

``` python
capacity = 1000
amount, pickup_point, pickup_location, delivery_point, delivery_location =[],[],[],[],[]
#,skills,priority
#delivery_service, delivery_time_windows
#pickup_service, pickup_time_windows, 
for row in vehicle_df.itertuples():
    #full trucks
    for i in range(int(row.number)-1):
        amount.append(f"[{capacity}]")
        pickup_point.append( row[4] )
        delivery_point.append( row[5] )
        pickup_location.append( f"[{dc_df.iloc[row.from_id,6]},{dc_df.iloc[row.from_id,5]}]" ) #lon, lat
        delivery_location.append( f"[{dc_df.iloc[row.to_id,6]},{dc_df.iloc[row.to_id,5]}]" ) 
    #fractional delivery
    amount.append(f"[{round(row.flow % capacity)}]")
    pickup_point.append( row[4] )
    delivery_point.append( row[5] )
    pickup_location.append( f"[{dc_df.iloc[row.from_id,6]},{dc_df.iloc[row.from_id,5]}]" ) #lon, lat
    delivery_location.append( f"[{dc_df.iloc[row.to_id,6]},{dc_df.iloc[row.to_id,5]}]" ) 
shipment_df = pd.DataFrame({"amount":amount, "pickup_point":pickup_point, "pickup_location":pickup_location, 
              "delivery_point":delivery_point, "delivery_location":delivery_location})
```

``` python
shipment_df["delivery_time_windows"] = "[[0,72000]]" #20 hours
shipment_df["pickup_time_windows"] = "[[0,72000]]" #20 hours
shipment_df["delivery_service"] = 600 #10 minutes 
shipment_df["pickup_service"] = 600
shipment_df["skills"] = "[0]"
shipment_df["priority"] = 100
shipment_df["pickup_index"] = None
shipment_df["delivery_index"] = None
shipment_df.head()
```

``` python
#トラックの台数
num_trucks = np.zeros(n, int)
for row in vehicle_df.itertuples():
    num_trucks[row.from_id] += row.number
num_trucks
```

``` python
name, start = [], []

for i in range(n):
    for v in range(num_trucks[i]):
        name.append(dc_df.iloc[i,0]+f":vehicle{v}" )
        start.append( f"[{dc_df.iloc[i,6]},{dc_df.iloc[i,5]}]" ) 
vehicle_df = pd.DataFrame({"name":name, "start":start})
vehicle_df["end"] = "[]"
vehicle_df["start_index"] = None
vehicle_df["end_index"] = None
vehicle_df["skills"] = "[0]"
vehicle_df["breaks"] = "[0,1,2]"
vehicle_df["capacity"] = f"[{capacity}]"
vehicle_df["time_window"] = "[0,46800]" #13 hours
```

``` python
#最大13h 4hに30分の休憩
description, tw = [],[]
t = 0
for i in range(3):
    description.append(f"break{i}")
    t += 4*3600
    tw.append( f"[{t},{t+1800}]" ) 

break_df = pd.DataFrame( {"description": description, "time_window": tw} )
#break_df.to_csv(folder+"break.csv")
break_df["service"] = 1800
break_df
```

### 配送最適化実行

``` python
matrix = False
job_df = pd.DataFrame()
model = build_model_for_vrp(job_df, shipment_df, vehicle_df, break_df)

input_dic, output_dic, error = optimize_vrp(model,matrix=matrix)
```

``` python
#output_dic["summary"]["cost"], output_dic["summary"]["distance"]
```

``` python
606255*8000/3600+ 20*12344977/1000
```

``` python
(606255*8000/3600+ 20*12344977/1000)/1633083
```

``` python
fig = make_fig_for_vrp(input_dic, output_dic, show_mode="route")
plotly.offline.plot(fig);
fig.write_html(f"map-sendo.html", include_plotlyjs=True)
route_df_dic, unassigned_job_df, unassigned_shipment_df = make_route_for_vrp(job_df, shipment_df, vehicle_df, break_df, output_dic)

fig = gannt_for_vrp(route_df_dic, start= "2020/01/01 8:00" )
plotly.offline.plot(fig);
```

求解したところ，全部で22台のトラックで処理できることが分かった．

総移動時間は606255秒，移動距離は12344977mであり，サービスネットワーク設計ではｍ1kmあたりの費用は20円,
1時間あたりの費用は8000円と設定していたので，

    606255*8000/3600+ 20*12344977/1000 = 1594133 円

となる．サービスネットワーク設計でのトラック費用は1633083円であったので，98%（誤差2%）であり，良い近似であることが確認される．

ガントチャートと，ルートの一部を描画すると，以下のようになる．休憩を法令に準じて４時間おきに３０分入れているので，
無駄な待ち時間が挿入されているが，その後に仕事が無い場合には，切れば良い．

また，ルートの一部を地図上に描画すると，少ない荷物は積み合わせ輸送していることが確認できる．

## シフト最適化

サービスネットワーク設計で成功したあなたに，今回は拠点でのシフト最適化の依頼が来た．
拠点では，入庫された荷物を行き先別に仕分けし，出庫する作業が行われているが，仕分け作業に膨大な時間が割かれている．
現状では，非常に忙しいピーク時間帯があるため多くのスタッフを当てている状態だ．
ただし，ピーク時以外は，行う作業がなく手持ち無沙汰になっている作業場もあるようだ．

ここでは，午前9時から21時の12時間を考え，スタッフの適切配置を考えることにする．現状では，30日の希望シフトを出してもらって，ほぼ希望通りに来てもらい，当日来た人数をもとに配置を経験的に(Excelで）決定しているが，
これにSCMOPTに含まれるOptShiftを用いた最適化を適用することにする．以下では，OptShiftの用語を用いて解説するものとし，OptShiftの用語にならい，作業員はスタッフ，作業はジョブと呼ぶものとする．

計画期間は30日とし，１ヶ月分の計画を建てることにする．これは，今までのルール通りにしたいという現場からの要望に基づくものである．
作業量は曜日によって異なるので，休日，曜日などの情報をもった日データを準備する．ちなみに，今回考える5月のスケジュールは休日が多いので，日本のカレンダーデータを用いることにする．

次に，9時から21時の期（1時間毎に時間を区切ったデータ）を生成する．
全部で12期の問題となる．

この会社では，休憩時間は４時間を超えると1時間の休憩をとるという規則がある．これは，以下の表のようになっており，5時間働いたら１時間休憩，10時間働いたら２時間休憩となっている．
また，期間(period)が3から始まるのは，一度出社したら，最低でも3時間は働くことを表している．

ジョブは，仕分けが中心だが，冷凍の仕分けと方面別の仕分け（2つ）に別れている．これらのジョブは，忙しくなる時間帯が異なるので，できればスタッフを複数のジョブに割り当てたいのだが，現状は忙しいときにはフル回転で，仕事がないときには休憩という感じになっている．

スタッフは全員時給1200円であり，人によって休日や最大稼働時間，出勤可能時間帯が異なる．
また，希望する休日や，できるジョブの種類も設定できるようになっている．
ここでは，20人のスタッフが出勤可能とする．

最後に，各ジョブに対する必要人数を現場からの聞き取りで収集する．ここでは，平日，日曜，休日で作業量が大きく異なることから，3通りの時間帯別の必要人数データを準備する．

以上のデータを入力し，バイト代の総額を目的関数と設定する．さらに，制約逸脱のペナルティ費用も設定し，なるべく制約を満たしてくれる解を探索することにする．

ヒヤリングをすると，各制約を逸脱したときのペナルティは，以下のように推定される．もちろん，これらのペナルティを色々変更して，なるべく現場にあった解を探索する必要がある．

- スタッフの人数の下限：
  人数が足りないとスタッフに過剰の負担がかかることになる．これが続くとバイトが集まらなくなる可能性があるので，1人足りないと1万円のペナルティと設定する．
- スタッフの人数の上限： １人多いとその分の時給が会社の損になる．
  時給は1200円だが，余裕をもって作業することができると，離職率が減るので，必要人数より1人多いと1000円のペナルティと設定する．
- ジョブチェンジ:
  スタッフは途中でジョブを変わることができる．当然，作業場は異なるので，移動する必要があり，移動時間分の時給が会社の損出になる．ここでは，作業場が比較的近いので，ペナルティは（頻繁なジョブチェンジを避けるために）1円と設定する．
- 休憩の違反：
  稼働時間に対して決められた休憩時間は「必ず」守るものとする．
  また，休憩が最初と最後に入ることは望ましくはないが，現場ではそれほど問題にしていないので，8000円と見積もることにする．
- 月間稼働日数：
  今月働きたい日数は，税金の関係で重要だ．これを逸脱すると来月以降に出勤してくれなくなるかもしれない．シフト最適化は動的な問題であり，今は1ヶ月の「静的な」問題を考えて，最適化しているが，来月以降も考慮して解かないと，年度末に大変なことになる（バイトが誰も来てくれない！）．よって，スタッフ自身が，自分で今月はこのくらい働こうという希望日数を入力し，それからの逸脱をペナルティとする．ここではそれを，5000円と設定する．

以上のパラメータをもとに最適化すると，以下の結果が得られる．
総費用は，116万程度でこれは現状より10%くらい低いようだ．スタッフの必要人数を超過した人数の合計は，30日間で264時間となった．ジョブは3種類あるので，ひとつのジョブあたり１日３人くらいの余裕が必要ということになる．
他の制約逸脱はジョブチェンジの回数があり，30日で69回となった．1日あたりにすると2回程度であり，作業場が近いため，問題はないと判断される．

例として，4日目のシフトのガントチャートと，各ジョブに割り当てられた人数と必要人数を図示してみる．
この最適化された結果では，余剰なスタッフが割り当てられているが，そこにジョブが入る可能性もあるので，このまま使うことにした．

このように月次のシフト最適化は比較的簡単に解くことができ，得られる結果も，様々なオペレーショナルな条件を満たしたものとなる．
より大規模な問題の場合には，さらに様々な工夫を追加する必要があるが，制約最適化ソルバーSCOPを使えば，それも比較的容易である．

\#hide
求解する前に，おおよその必要人数とスタッフの希望勤務時間の関係を調べておくこともできる．

\#hide 本事例で用いたシステムは，シフト最適化システム OptShiftです．
OptShiftは，サプライ・チェイン最適化システムSCMOPT に含まれています．