# ケーススタディ2


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

## 料金表

2021(2024)年から施行された標準料金を読み込み，pickleに保管しておく．

``` python
df.tail()
```

<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">km</th>
<th data-quarto-table-cell-role="th">area</th>
<th data-quarto-table-cell-role="th">2</th>
<th data-quarto-table-cell-role="th">4</th>
<th data-quarto-table-cell-role="th">10</th>
<th data-quarto-table-cell-role="th">20</th>
</tr>
</thead>
<tbody>
<tr>
<td data-quarto-table-cell-role="th">215</td>
<td>180</td>
<td>沖縄</td>
<td>38480</td>
<td>44990</td>
<td>59890</td>
<td>77070</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">216</td>
<td>190</td>
<td>沖縄</td>
<td>40010</td>
<td>46770</td>
<td>62270</td>
<td>80170</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">217</td>
<td>200</td>
<td>沖縄</td>
<td>41540</td>
<td>48540</td>
<td>64660</td>
<td>83280</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">218</td>
<td>210</td>
<td>沖縄</td>
<td>3050</td>
<td>3530</td>
<td>4700</td>
<td>6110</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">219</td>
<td>220</td>
<td>沖縄</td>
<td>7610</td>
<td>8810</td>
<td>11740</td>
<td>15270</td>
</tr>
</tbody>
</table>

</div>

``` python
def simple_cycles2(G, limit):
    subG = type(G)(G.edges())
    sccs = list(nx.strongly_connected_components(subG))
    while sccs:
        scc = sccs.pop()
        startnode = scc.pop()
        path = [startnode]
        blocked = set()
        blocked.add(startnode)
        stack = [(startnode, list(subG[startnode]))]

        while stack:
            thisnode, nbrs = stack[-1]

            if nbrs and len(path) < limit:
                nextnode = nbrs.pop()
                if nextnode == startnode:
                    yield path[:]
                elif nextnode not in blocked:
                    path.append(nextnode)
                    stack.append((nextnode, list(subG[nextnode])))
                    blocked.add(nextnode)
                    continue
            if not nbrs or len(path) >= limit:
                blocked.remove(thisnode)
                stack.pop()
                path.pop()
        subG.remove_node(startnode)
        H = subG.subgraph(scc)
        sccs.extend(list(nx.strongly_connected_components(H)))
```

``` python
dem = od_df.SHASHU_SAIDAISEKISAIRYO_CD_1.astype(int, errors="ignore")
dem.hist()
```

![](16case2_files/figure-commonmark/cell-16-output-1.png)

``` python
#中継地点の描画
#地点名も出すようにする！
data =[go.Scattermapbox(
        lat=lat_head,
        lon=lon_head,
        marker = {'size': num},
        hoverinfo='text',
        hovertext= head, 
        text= head,
        mode="markers",
        name="車庫"
    )]
data.append(    
    go.Scattermapbox(
        lat=lat_switch,
        lon=lon_switch,
        marker = {'size': 5},
        hoverinfo='text',
        text= switch,
        mode="markers",
        name="中継地点"
    )
)
fig = go.Figure(data=data, layout=layout)
plotly.offline.plot(fig);
```

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

ロジスティックネットワーク設計問題を例として，データの前処理やクラスタリングの方法による違いを調べる．

まずはデータを読み込む．列名はSCMOPTにあわせて変えておく．

``` python
cust_df = pd.read_csv(folder+ "Cust.csv")
cust_df.columns= ["name","lat","lon"]
cust_df["lat"] = cust_df["lat"].astype(float)
cust_df["lon"] = cust_df["lon"].astype(float)
#cust_df.head()
```

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

``` python
total_demand_df = pd.read_csv(folder+ "Cust-Prod.csv")
total_demand_df.columns=["cust","prod","demand"]
#total_demand_df.head()
```

``` python
dc_df = pd.read_csv(folder+ "DC.csv")
dc_df.columns=["name","fc","vc","lb","ub","lat","lon"]
dc_df.lb = 0.
#dc_df.head()
```

``` python
plnt_df = pd.read_csv(folder+ "Plnt.csv")
plnt_df.columns=["name","lat","lon"]
#plnt_df.head()
```

``` python
prod_df = pd.read_csv(folder+ "Prod.csv")
prod_df.columns=["name","weight"]
#prod_df
```

``` python
plnt_prod_df = pd.read_csv(folder+"Plnt-Prod.csv")
plnt_prod_df.columns=["plnt","prod","ub"]
#plnt_prod_df.head()
```

### 同一座標の顧客

同じ緯度・経度をもつ地点が多数存在する．これは，K-d木を用いると容易に確認できる．
ここでは，そのままのデータで解くが，k-メディアン法は，同じ地点があると2つともクラスターの中心（メディアン）として選択してしまうので，注意を要する．

``` python
df = enumerate_near_points(cust_df, 0.0001)
#df
```

### 地点間の道路距離と移動時間を計算

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

``` python
#ネットワーク生成
trans_df, graph, position  = make_network_using_road(cust_df, dc_df, plnt_df, durations, distances, plnt_dc_threshold = 10000., dc_cust_threshold = 200. )
# 大圏距離で代用する場合には，以下を生かす
#trans_df, graph, position  = make_network(cust_df, dc_df, plnt_df, plnt_dc_threshold = 10000., dc_cust_threshold = 500. )
```

### 需要の分布を確認

需要の分布を確認すると，ほとんどが2000以下なので，除外する．
0の需要だけを除外した場合は，数値のばらつきがあまりにも大きいために，Gurobiでさえ不安定になり，解くことができない（もしくは非常に時間がかかる）．

``` python
fig = px.histogram(total_demand_df,"demand")
plotly.offline.plot(fig);
```

``` python
total_demand_df = total_demand_df[ total_demand_df.demand >= 2000. ]
```

### 需要がない顧客を除外

これで773点の問題例が，約300点になった．

``` python
print(len(cust_df))
cust_df = remove_zero_total_demand(total_demand_df, cust_df)
print(len(cust_df))
```

    773
    299

``` python
flow_df, dc_df, cost_df, violation_df, status = solve_lnd(prod_df, cust_df, dc_df, plnt_df, plnt_prod_df, total_demand_df, 
                               trans_df, dc_num= (5,13), single_sourcing=True, max_cpu = 100)
```

``` python
cost_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">cost</th>
<th data-quarto-table-cell-role="th">value</th>
</tr>
</thead>
<tbody>
<tr>
<td data-quarto-table-cell-role="th">0</td>
<td>transportation (plant to dc)</td>
<td>9.736955e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">1</td>
<td>delivey (dc to customer)</td>
<td>3.606422e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">2</td>
<td>dc fixed</td>
<td>-0.000000e+00</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">3</td>
<td>dc variable</td>
<td>6.524372e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">4</td>
<td>infeasible penalty</td>
<td>0.000000e+00</td>
</tr>
</tbody>
</table>

</div>

``` python
#費用の可視化
fig = px.bar(cost_df, y="cost", x="value", orientation='h')
plotly.offline.plot(fig);
```

``` python
cost_df.value.sum()
```

    1986774926.5981956

``` python
fig = show_optimized_network(cust_df, dc_df, plnt_df, prod_df, flow_df, position)
plotly.offline.plot(fig);
```

### 階層的クラスタリング法で集約

データを読み直して，地点間の移動時間を計算した後で，クラスタリングを行いデータを集約する．

最適化後に，非集約し，本当の配送費用を計算する．

まずは，階層的クラスタリング法を適用する．これは，道路距離と時間を用いる．

``` python
#clusterning 
weight_of_prod ={p:w for (p,w) in zip(prod_df.name, prod_df.weight) } #prepare weight dic.
weight_of_cust = defaultdict(int)
for row in total_demand_df.itertuples():
    weight_of_cust[row.cust] += weight_of_prod[row.prod]*row.demand 
weight = [weight_of_cust[i] for i in weight_of_cust] 
durations, distances, node_df = compute_durations(cust_df, None)
```

``` python
fig = elbow_method(cust_df, weight, n_lb = 3, n_ub = 60, repetitions=1, method="hierarchical", durations=durations)
plotly.offline.plot(fig);
```

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

    382277301752.6439

``` python
fig = show_optimized_continuous_network(cust_df, X, Y, partition, weight= weight)
plotly.offline.plot(fig);
```

``` python
aggregated_cust_df, aggregated_total_demand_df =  make_aggregated_df(cust_df, None, total_demand_df, X, Y, partition, weight)
aggregated_dc_df = generate_dc(aggregated_cust_df, aggregated_total_demand_df, prod_df, num_dc=10, lb_ratio=0.0, ub_ratio=10.3,vc_bounds=(19,21),fc_bounds=(0,1))
durations, distances, node_df = compute_durations(aggregated_cust_df, plnt_df)
aggregated_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 = 500. )
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, 
                               aggregated_trans_df, dc_num= (5,13), single_sourcing=False, max_cpu = 100)
```

``` python
fig = show_optimized_network(aggregated_cust_df, aggregated_dc_df, plnt_df, prod_df, flow_df, position)
plotly.offline.plot(fig);
```

``` python
cost_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">cost</th>
<th data-quarto-table-cell-role="th">value</th>
</tr>
</thead>
<tbody>
<tr>
<td data-quarto-table-cell-role="th">0</td>
<td>transportation (plant to dc)</td>
<td>9.761203e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">1</td>
<td>delivey (dc to customer)</td>
<td>2.418022e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">2</td>
<td>dc fixed</td>
<td>-0.000000e+00</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">3</td>
<td>dc variable</td>
<td>6.505497e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">4</td>
<td>infeasible penalty</td>
<td>0.000000e+00</td>
</tr>
</tbody>
</table>

</div>

上の費用のうち，顧客への配送費用は，集約した点に対するものなので，本来の費用ではない．本来の費用を再計算する．

``` python
# Single Source を仮定
from_dc = {} #集約した顧客に入る倉庫の名称を入れた辞書（集約した顧客名をキーとする）
for row in flow_df.itertuples():
    if row.to_node[:5]=="Cust_":
        from_dc[row.to_node[5:]] = row.from_node[3:]
#開設した集約倉庫に最も近い顧客を探索
from scipy.spatial import KDTree 
points = [(i,j)  for (i,j)  in zip(cust_df.lat, cust_df.lon)]
tree = KDTree(points)
dc_cust_dic={} #集約した倉庫の名前を入れると，それに最も近い，もとの顧客indexを返す辞書
for row in dc_df[ dc_df.open_close==1 ].itertuples():
    dis, near = tree.query( (row.lat, row.lon) )
    dc_cust_dic[row.name] = cust_df.loc[int(near),"index"]
durations, distances, node_df = compute_durations(cust_df, None)
#dc_set = set(from_dc.values())    #from_dc: 各集約顧客が，どのDCから来ているかを表す辞書
#agg_cust_set = set(from_dc.keys())
dc_per_dis = 10./4000
dc_per_time = 8000./4000
total_del = 0.
for row in aggregated_cust_df.itertuples():
    dc_name = from_dc[row.name]
    #print(dc_name, row.name)
    i = dc_cust_dic[dc_name]
    for c in row.customers:
        j = cust_df.loc[c,"index"]
        if distances[i][j] < 9999999:
            dis = distances[i][j]/1000.
            time = durations[i][j]/3600.
            cost_per_kg = dis*dc_per_dis + time*dc_per_time
            total_del += weight[j]*cost_per_kg 
total_del
```

    436088454.2303338

配送費用は，クラスタリングした場合より大きくなり，総費用は以下のように増加する．
これは，需要が2000未満のものを削除して最適化した場合の1.4%増であり，良い近似になっていると考えられる．

``` python
total = cost_df.value[0] + cost_df.value[3] + total_del
print(total, total/1986774926)
```

    2020658919.3316379 1.0170547719765404

### k-メディアン法で集約

k-median法で集約する
k-メディアン法も道路距離と移動時間をもとに計算を行う．同じ座標の地点があると，何も割り当てられないクラスターができることがあるので，多少多めに（この場合は１つ多く）施設数を入れておく．

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

``` python
X, Y, partition, best_ub, lb_list, ub_list, phi_list = solve_k_median(cust_df, weight, durations, num_of_facilities=51, 
                                 max_iter=10000, max_lr=0.01, moms=(0.85,0.95), convergence=1e-5,  lr_find = True, adam=False)
```

``` python
fig = plot_k_median_lr_find(phi_list, lb_list)
plotly.offline.plot(fig);
#lr_finder => 0.01 , adam => 50000
```

``` python
X, Y, partition, best_ub, lb_list, ub_list, phi_list = solve_k_median(cust_df, weight, durations, num_of_facilities=51, 
                                 max_iter=2000, max_lr=0.05, moms=(0.85,0.95), convergence=1e-5,  lr_find = False, adam=False)
```

``` python
fig = plot_k_median(lb_list, ub_list)
plotly.offline.plot(fig);
```

``` python
fig = show_optimized_continuous_network(cust_df, X, Y, partition, weight= weight)
plotly.offline.plot(fig);
```

``` python
aggregated_cust_df, aggregated_total_demand_df =  make_aggregated_df(cust_df, None, total_demand_df, X, Y, partition, weight)
aggregated_dc_df = generate_dc(aggregated_cust_df, aggregated_total_demand_df, prod_df, num_dc=10, lb_ratio=0.0, ub_ratio=10.3,vc_bounds=(19,21),fc_bounds=(0,1))
durations, distances, node_df = compute_durations(aggregated_cust_df, plnt_df)
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 = 500. )
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= (5,13), single_sourcing=True, max_cpu = 100)
```

``` python
cost_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">cost</th>
<th data-quarto-table-cell-role="th">value</th>
</tr>
</thead>
<tbody>
<tr>
<td data-quarto-table-cell-role="th">0</td>
<td>transportation (plant to dc)</td>
<td>9.647053e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">1</td>
<td>delivey (dc to customer)</td>
<td>2.503247e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">2</td>
<td>dc fixed</td>
<td>-0.000000e+00</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">3</td>
<td>dc variable</td>
<td>6.269581e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">4</td>
<td>infeasible penalty</td>
<td>-0.000000e+00</td>
</tr>
</tbody>
</table>

</div>

``` python
fig = show_optimized_network(aggregated_cust_df, aggregated_dc_df, plnt_df, prod_df, flow_df, position)
plotly.offline.plot(fig);
```

``` python
aggregated_cust_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">name</th>
<th data-quarto-table-cell-role="th">lat</th>
<th data-quarto-table-cell-role="th">lon</th>
<th data-quarto-table-cell-role="th">customers</th>
<th data-quarto-table-cell-role="th">demand</th>
</tr>
</thead>
<tbody>
<tr>
<td data-quarto-table-cell-role="th">0</td>
<td>cust0</td>
<td>34.577843</td>
<td>135.438274</td>
<td>[51, 112, 125, 151, 152, 206, 256, 265, 266, 267, 270, 285]</td>
<td>8.997956e+06</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">1</td>
<td>cust1</td>
<td>41.826084</td>
<td>140.725339</td>
<td>[1]</td>
<td>1.824079e+06</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">2</td>
<td>cust2</td>
<td>36.031805</td>
<td>139.630817</td>
<td>[2, 15, 64, 81, 85, 129, 238, 293]</td>
<td>4.023094e+06</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">3</td>
<td>cust3</td>
<td>34.674853</td>
<td>135.634251</td>
<td>[8, 35, 82, 86, 103, 104, 120, 150, 259, 264]</td>
<td>8.146061e+06</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">4</td>
<td>cust4</td>
<td>35.982864</td>
<td>139.860838</td>
<td>[11, 68, 77, 154, 174, 182, 260, 275, 287]</td>
<td>6.499158e+06</td>
</tr>
</tbody>
</table>

</div>

``` python
# Single Source を仮定
try:
    cust_df.reset_index(inplace=True)
except:
    pass
from_dc = {} #集約した顧客に入る倉庫の名称を入れた辞書（集約した顧客名をキーとする）
for row in flow_df.itertuples():
    if row.to_node[:5]=="Cust_":
        from_dc[row.to_node[5:]] = row.from_node[3:]
#開設した集約倉庫に最も近い顧客を探索
from scipy.spatial import KDTree 
points = [(i,j)  for (i,j)  in zip(cust_df.lat, cust_df.lon)]
tree = KDTree(points)
dc_cust_dic={} #集約した倉庫の名前を入れると，それに最も近い，もとの顧客indexを返す辞書
for row in dc_df[ dc_df.open_close==1 ].itertuples():
    dis, near = tree.query( (row.lat, row.lon) )
    dc_cust_dic[row.name] = cust_df.loc[int(near),"index"]
durations, distances, node_df = compute_durations(cust_df, None)
#dc_set = set(from_dc.values())    #from_dc: 各集約顧客が，どのDCから来ているかを表す辞書
#agg_cust_set = set(from_dc.keys())
dc_per_dis = 10./4000
dc_per_time = 8000./4000
total_del = 0.
for row in aggregated_cust_df.itertuples():
    try:
        dc_name = from_dc[row.name]
    except:
        continue
    i = dc_cust_dic[dc_name]
    for c in row.customers:
        j = cust_df.loc[c,"index"]
        if distances[i][j] < 9999999:
            dis = distances[i][j]/1000.
            time = durations[i][j]/3600.
            cost_per_kg = dis*dc_per_dis + time*dc_per_time
            total_del += weight[j]*cost_per_kg 
total_del
```

    429417810.7833409

``` python
cost_df.value[0] + cost_df.value[3] + total_del
```

    2021081255.6539354

``` python
2005167825/1986774926
```

    1.0092576661600168

### k-means法で集約

k-means法は重心法ともよばれ，座標をもとに平面上の適当な地点に施設を配置する．

``` python
weight_of_prod ={p:w for (p,w) in zip(prod_df.name, prod_df.weight) } #prepare weight dic.
weight_of_cust = defaultdict(int)
for row in total_demand_df.itertuples():
    weight_of_cust[row.cust] += weight_of_prod[row.prod]*row.demand 
weight = [weight_of_cust[i] for i in weight_of_cust] 
durations, distances, node_df = compute_durations(cust_df, None)
```

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

``` python
X, Y, partition, cost = kmeans(cust_df, weight, num_of_facilities =50)
```

``` python
fig = show_optimized_continuous_network(cust_df, X, Y, partition, weight= weight)
plotly.offline.plot(fig);
```

``` python
aggregated_cust_df, aggregated_total_demand_df =  make_aggregated_df(cust_df, None, total_demand_df, X, Y, partition, weight)
aggregated_dc_df = generate_dc(aggregated_cust_df, aggregated_total_demand_df, prod_df, num_dc=10, lb_ratio=0.0, ub_ratio=10.3,vc_bounds=(19,21),fc_bounds=(0,1))
durations, distances, node_df = compute_durations(aggregated_cust_df, plnt_df)
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 = 500. )
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= (5,13), single_sourcing=True, max_cpu = 100)
```

``` python
cost_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">cost</th>
<th data-quarto-table-cell-role="th">value</th>
</tr>
</thead>
<tbody>
<tr>
<td data-quarto-table-cell-role="th">0</td>
<td>transportation (plant to dc)</td>
<td>9.935685e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">1</td>
<td>delivey (dc to customer)</td>
<td>2.343569e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">2</td>
<td>dc fixed</td>
<td>-0.000000e+00</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">3</td>
<td>dc variable</td>
<td>6.559424e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">4</td>
<td>infeasible penalty</td>
<td>-0.000000e+00</td>
</tr>
</tbody>
</table>

</div>

``` python
# Single Source を仮定
#cust_df.reset_index(inplace=True)
from_dc = {} #集約した顧客に入る倉庫の名称を入れた辞書（集約した顧客名をキーとする）
for row in flow_df.itertuples():
    if row.to_node[:5]=="Cust_":
        from_dc[row.to_node[5:]] = row.from_node[3:]
#開設した集約倉庫に最も近い顧客を探索
from scipy.spatial import KDTree 
points = [(i,j)  for (i,j)  in zip(cust_df.lat, cust_df.lon)]
tree = KDTree(points)
dc_cust_dic={} #集約した倉庫の名前を入れると，それに最も近い，もとの顧客indexを返す辞書
for row in dc_df[ dc_df.open_close==1 ].itertuples():
    dis, near = tree.query( (row.lat, row.lon) )
    dc_cust_dic[row.name] = cust_df.loc[int(near),"index"]
durations, distances, node_df = compute_durations(cust_df, None)
#dc_set = set(from_dc.values())    #from_dc: 各集約顧客が，どのDCから来ているかを表す辞書
#agg_cust_set = set(from_dc.keys())
dc_per_dis = 10./4000
dc_per_time = 8000./4000
total_del = 0.
for row in aggregated_cust_df.itertuples():
    #try:
    dc_name = from_dc[row.name]
    #except:
    #    print("error")
    #    continue
    i = dc_cust_dic[dc_name]
    for c in row.customers:
        j = cust_df.loc[c,"index"]
        if distances[i][j] < 9999999:
            dis = distances[i][j]/1000.
            time = durations[i][j]/3600.
            cost_per_kg = dis*dc_per_dis + time*dc_per_time
            total_del += weight[j]*cost_per_kg 
        else:
            print(i,j,"error")
#total_del
cost_df.value[0] + cost_df.value[3] + total_del
```

    304 413 error
    304 420 error
    304 429 error

    2022563380.9275873

道路距離が非常に大きい地点間で配送を行っている．
これは，離島（対馬）への道路距離が計算できないためである．

``` python
cost_df.value[0] + cost_df.value[3] + total_del
```

    2022563380.9275873

``` python
fig = show_optimized_network(aggregated_cust_df, aggregated_dc_df, plnt_df, prod_df, flow_df, position)
plotly.offline.plot(fig);
```

### 反復 Weiszfeld法で集約

k-means法が施設と顧客の距離の（重み付き）自乗和の最小化であるのに対して，
Weiszfeld法は（重み付き）直線距離の和を最小化する．

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

``` python
X, Y, partition, cost = repeated_weiszfeld(cust_df, weight, 
            num_of_facilities = 50, epsilon=0.0001, max_iter = 30, numpy=True, seed=2)
print(cost)
```

    2462506476.093891
    CPU times: user 2min 38s, sys: 83.3 ms, total: 2min 38s
    Wall time: 2min 38s

``` python
fig = show_optimized_continuous_network(cust_df, X, Y, partition, weight= weight)
plotly.offline.plot(fig);
```

``` python
aggregated_cust_df, aggregated_total_demand_df =  make_aggregated_df(cust_df, None, total_demand_df, X, Y, partition, weight)
aggregated_dc_df = generate_dc(aggregated_cust_df, aggregated_total_demand_df, prod_df, num_dc=10, lb_ratio=0.0, ub_ratio=10.3,vc_bounds=(19,21),fc_bounds=(0,1))
durations, distances, node_df = compute_durations(aggregated_cust_df, plnt_df)
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 = 500. )
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= (5,13), single_sourcing=True, max_cpu = 100)
```

``` python
cost_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">cost</th>
<th data-quarto-table-cell-role="th">value</th>
</tr>
</thead>
<tbody>
<tr>
<td data-quarto-table-cell-role="th">0</td>
<td>transportation (plant to dc)</td>
<td>9.902503e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">1</td>
<td>delivey (dc to customer)</td>
<td>2.799122e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">2</td>
<td>dc fixed</td>
<td>-0.000000e+00</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">3</td>
<td>dc variable</td>
<td>6.400399e+08</td>
</tr>
<tr>
<td data-quarto-table-cell-role="th">4</td>
<td>infeasible penalty</td>
<td>-0.000000e+00</td>
</tr>
</tbody>
</table>

</div>

``` python
fig = show_optimized_network(aggregated_cust_df, aggregated_dc_df, plnt_df, prod_df, flow_df, position)
plotly.offline.plot(fig);
```

``` python
# Single Source を仮定
cust_df.reset_index(inplace=True)
from_dc = {} #集約した顧客に入る倉庫の名称を入れた辞書（集約した顧客名をキーとする）
for row in flow_df.itertuples():
    if row.to_node[:5]=="Cust_":
        from_dc[row.to_node[5:]] = row.from_node[3:]
#開設した集約倉庫に最も近い顧客を探索
from scipy.spatial import KDTree 
points = [(i,j)  for (i,j)  in zip(cust_df.lat, cust_df.lon)]
tree = KDTree(points)
dc_cust_dic={} #集約した倉庫の名前を入れると，それに最も近い，もとの顧客indexを返す辞書
for row in dc_df[ dc_df.open_close==1 ].itertuples():
    dis, near = tree.query( (row.lat, row.lon) )
    dc_cust_dic[row.name] = cust_df.loc[int(near),"index"]
durations, distances, node_df = compute_durations(cust_df, None)
#dc_set = set(from_dc.values())    #from_dc: 各集約顧客が，どのDCから来ているかを表す辞書
#agg_cust_set = set(from_dc.keys())
dc_per_dis = 10./4000
dc_per_time = 8000./4000
total_del = 0.
for row in aggregated_cust_df.itertuples():
    try:
        dc_name = from_dc[row.name]
    except:
        print("error")
        continue
    i = dc_cust_dic[dc_name]
    for c in row.customers:
        j = cust_df.loc[c,"index"]
        if distances[i][j] < 9999999:
            dis = distances[i][j]/1000.
            time = durations[i][j]/3600.
            cost_per_kg = dis*dc_per_dis + time*dc_per_time
            total_del += weight[j]*cost_per_kg 
        else:
            print(i,j,"error")
#total_del
cost_df.value[0] + cost_df.value[3] + total_del
```

    122 413 error
    122 420 error
    122 429 error

    1992036715.0329504

k-means法と同様に，道路距離が非常に大きい地点間で配送を行っている．
これは，離島（対馬）への道路距離が計算できないためである．