ケーススタディ2

SCMOPTのケーススタディ 2

料金表

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

df.tail()

	km	area	2	4	10	20
215	180	沖縄	38480	44990	59890	77070
216	190	沖縄	40010	46770	62270	80170
217	200	沖縄	41540	48540	64660	83280
218	210	沖縄	3050	3530	4700	6110
219	220	沖縄	7610	8810	11740	15270

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

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

#中継地点の描画
#地点名も出すようにする！
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にあわせて変えておく．

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

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

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

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

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

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

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つともクラスターの中心（メディアン）として選択してしまうので，注意を要する．

df = enumerate_near_points(cust_df, 0.0001)
#df

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

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

#ネットワーク生成
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でさえ不安定になり，解くことができない（もしくは非常に時間がかかる）．

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

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

需要がない顧客を除外

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

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

773
299

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)

cost_df

	cost	value
0	transportation (plant to dc)	9.736955e+08
1	delivey (dc to customer)	3.606422e+08
2	dc fixed	-0.000000e+00
3	dc variable	6.524372e+08
4	infeasible penalty	0.000000e+00

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

cost_df.value.sum()

1986774926.5981956

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

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

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

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

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

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

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

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

382277301752.6439

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

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)

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

cost_df

	cost	value
0	transportation (plant to dc)	9.761203e+08
1	delivey (dc to customer)	2.418022e+08
2	dc fixed	-0.000000e+00
3	dc variable	6.505497e+08
4	infeasible penalty	0.000000e+00

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

# 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%増であり，良い近似になっていると考えられる．

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

2020658919.3316379 1.0170547719765404

k-メディアン法で集約

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

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

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)

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

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)

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

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

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)

cost_df

	cost	value
0	transportation (plant to dc)	9.647053e+08
1	delivey (dc to customer)	2.503247e+08
2	dc fixed	-0.000000e+00
3	dc variable	6.269581e+08
4	infeasible penalty	-0.000000e+00

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

aggregated_cust_df.head()

	name	lat	lon	customers	demand
0	cust0	34.577843	135.438274	[51, 112, 125, 151, 152, 206, 256, 265, 266, 267, 270, 285]	8.997956e+06
1	cust1	41.826084	140.725339	[1]	1.824079e+06
2	cust2	36.031805	139.630817	[2, 15, 64, 81, 85, 129, 238, 293]	4.023094e+06
3	cust3	34.674853	135.634251	[8, 35, 82, 86, 103, 104, 120, 150, 259, 264]	8.146061e+06
4	cust4	35.982864	139.860838	[11, 68, 77, 154, 174, 182, 260, 275, 287]	6.499158e+06

# 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

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

2021081255.6539354

2005167825/1986774926

1.0092576661600168

k-means法で集約

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

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)

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

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

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

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)

cost_df

	cost	value
0	transportation (plant to dc)	9.935685e+08
1	delivey (dc to customer)	2.343569e+08
2	dc fixed	-0.000000e+00
3	dc variable	6.559424e+08
4	infeasible penalty	-0.000000e+00

# 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

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

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

2022563380.9275873

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法は（重み付き）直線距離の和を最小化する．

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

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

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

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)

cost_df

	cost	value
0	transportation (plant to dc)	9.902503e+08
1	delivey (dc to customer)	2.799122e+08
2	dc fixed	-0.000000e+00
3	dc variable	6.400399e+08
4	infeasible penalty	-0.000000e+00

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

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