サービスネットワーク設計 SENDO (SErvice Network Design Optimizer)

簡易システムの紹介ビデオ

サービスネットワーク設計最適化システム SENDO

YouTubeVideo("GY0L8frFNko")

モデルの定式化

サービスネットワーク設計問題の定式化を示す. ターミナル間の需要量の予測値をもとに,各ターミナルに対して,他のターミナル行きの荷物を,どのターミナルに輸送すれば良いのかを指示する. 以下に,定式化を行うために必要な記号と集合を説明した後に,定式化を示す.

集合:

  • $G=(N,A)$: 有向グラフ:;点集合はターミナルの集合,枝集合は移動可能な点対の集合.

  • $K$: 発ターミナル・着ターミナルの組(これを品種と呼ぶ)の集合.

  • $P_k$: 品種 $k \in K$ に対するパスの集合. すべてのパスの集合を $P$ とする.

  • $A_{pk}$: 品種 $k$ に対するパス $p$ (これをパス $p,k$ と呼ぶ)が使用する枝の集合.

パラメータ:

  • $o_k, d_k$: 品種 $k$ の発ターミナルと着ターミナル
  • $D_{k}$: 品種 $k$ の需要量(発ターミナルから着ターミナルまでの荷量)
  • $c_{ij}$: 枝 $(i,j)$ の固定費用(トラック1台あたりの移動費用). 簡単のため最も多い10トントラックのみと仮定する.

  • $C_{pk}$: 品種 $k$ に対するパス $p$ の変動費用(需要量 $1$ 単位あたりの費用)

  • $M_{ij}$: トラック1台あたりの容量
  • $Succ_{pk}$: 品種 $k$ に対するパス $p$ の最初の枝を除去したパス(後続パス)

変数:

  • $x_{pk}$: 品種 $k$ に対してパス $p$ が使われるとき $1$,それ以外のとき $0$
  • $y_{ij}$: 枝 $(i,j)$ 上を移動するトラックの台数(整数変数)

定式化:

$$ \begin{array}{l l l } minimize & \sum_{(i,j) \in A} c_{ij} y_{ij} + \sum_{k \in K} \sum_{p \in P_k} C_{pk} D_{k} x_{pk} & \\ s.t. & \sum_{p \in P_k} x_{pk} = 1 & \forall k \in K \\ & \sum_{p,k:(i,j) \in A_{pk}} D_{k} x_{pk} \leq M_{ij} y_{ij} & \forall (i,j) \in A\\ & x_{pk} \leq x_{p',k'} & \forall (p',k') = Succ_{pk}\\ &\lceil D_k/M_{ij} \rceil x_{pk} \leq y_{ij} & \forall (i,j) \in A_{pk}, p \in P_k, k \in K\\ & \sum_{(i,j) \in (S,N \setminus S)} y_{ij} \geq \lceil \sum_{k: o_k \in S, d_k \in N \setminus S } D_k/M \rceil & \forall S \subset N, S \neq \emptyset \end{array} $$

制約の意味は以下の通り。

  1. パス選択制約:各品種はいずれかのパスによって流さなければならない.

  2. 容量制約: 各枝上を流れる品種の量の合計は,その枝上を移動したトラックの容量の合計を超えてはいけない.

  3. 入木制約: パス $p,k$ $=(n_1,n_2,\ldots, n_m)$ が選ばれたときには,その後続パス $p', k'$ $=(n_2,\ldots, n_m)$も選択されなければならない. つまり、同じ着地をもつ品種が一度合流したら,同じパス上を移動しなければならない.

  4. 強化制約: 品種 $k \in K$ に対するパスの集合は,ちょうど1つ選択する必要があるが,それらのパスが通過する枝 $(i,j)$ に対して需要 $D_k$ を満たす分のトラックが移動しなければならない.

  5. カットセット制約: 点の部分集合 $S$ とそれ以外の点集合 $N \setminus S$ へ移動する需要量を満たすだけのトラックが,$S$ から $N \setminus S$ へ移動しなければならない.これは,すべてのトラック容量 $M_{ij}$ が同一量 $M$ であるときに使われ,カットセット制約とよばれる.

データの読み込み

  • cust_df : 顧客データフレーム
  • od_df : OD需要データフレーム
od_df = pd.read_csv(folder +"od.csv", index_col=0)
Demand = od_df.values
#Demand
cust_df = pd.read_csv(folder+"Cust.csv", index_col=0)
n = len(cust_df)
cust_df.head()
name lat lon
id
1 札幌市 43.06417 141.34694
2 青森市 40.82444 140.74000
3 盛岡市 39.70361 141.15250
4 仙台市 38.26889 140.87194
5 秋田市 39.71861 140.10250
# OD需要のヒストグラムの描画
trace1 = go.Histogram(
    x= [ Demand[i,j] for i in range(n) for j in range(n) ], 
    opacity=0.5,
    name="OD間需要"
)

data = [trace1]
layout = go.Layout(
    title='OD間需要量',
    xaxis=dict(
        title="需要量"
    ),
    yaxis=dict(
        title='件数'
    )
)
fig = go.Figure(data=data, layout=layout)
#fig.show()

第k最短路を解くための関数

引数:

  • G : NetworkXのグラフ
  • source : 始点
  • sink : 終点
  • k : 第k最短路(単純パス)を求める。
  • weight : 移動距離を評価するための枝の属性名

返値:

  • cost_list : k本のパスの費用を入れたリスト
  • path_list : k本のパス(点列のタプル)を入れたリスト

k_th_sp[source]

k_th_sp(G, source, sink, k, weight='weight')

Find k-th shortest paths and returns path and its cost

サービスネットワーク設計モデル

サービスネットワーク設計問題に対するパス型のモデル

引数:

  • K : 発生点と到着点の対を入れたリスト
  • all_paths = set of all paths
  • paths[o,d] = set of all paths from o to d
  • path_od[p] = dictionary that maps a path (tuple) to its origin-destination pair (tuple)
  • path_id[p] = dictionary that maps a path to its ID
  • paths_through[v,w] = set of paths through edge (v,w)

sndp[source]

sndp(K, all_paths, paths, path_id, path_od, paths_through, Demand, Distance, transfer_cost, capacity, relax=False, cpu=10.0)

サービスネットワーク設計問題に対するパス型のモデル

add_shortest_paths[source]

add_shortest_paths(G, K, all_paths, paths, path_id, path_od, paths_through, k=1)

各地点間の第k最短路を計算し、追加する関数

メイン

solve_sndp[source]

solve_sndp(cust_df, Demand, Distance, transfer_cost=1.0, capacity=10000.0, cpu=10.0, Scaling=True, k=10, alpha=0.5, max_iter=100)

サービスネットワーク設計問題を解くためのメイン関数

使用例

cust_df = pd.read_csv(folder+"Cust.csv", index_col=0)
n = len(cust_df)
od_df = pd.read_csv(folder +"od.csv", index_col=0)

n= 10
cust_df = cust_df.iloc[:n,:]
od_df = od_df.iloc[:n,:n]

Demand = od_df.values

Distance = np.zeros( (n,n) )
for i, row1 in enumerate(cust_df.itertuples()):
    for j, row2 in enumerate(cust_df.itertuples()):
        Distance[i,j] = distance( (row1.lat, row1.lon), (row2.lat,row2.lon)).km

# just MIP
x, y, pos, paths, path_id = solve_sndp(cust_df, Demand, Distance, transfer_cost = 1.0, capacity = 10000., 
                                        cpu = 100., Scaling=False, k=5, alpha=0.5)

# slope scaling + column generation
#x, y, pos, paths, path_id = solve_sndp(cust_df, Demand, Distance, transfer_cost = 1.0, capacity = 10000., 
#                                        cpu = 10., Scaling=True, k=5, alpha=0.9, max_iter=10)
# of paths 450
obj fun. val. = 57831.26947215112

可視化関数

draw_snd[source]

draw_snd(cust_df, x, y, pos, paths, path_id, destination)

可視化関数

fig = draw_snd(cust_df, x, y, pos, paths, path_id, destination="札幌市")
#fig.show()
#Image("../figure/draw_snd.png")
# #強化制約(パス変数とトラック台数変数との関係式)
# force ={}
# for (b1,b2) in K:
#     #print("Variablle: now processing",b1.name,b2.name)
#     for (count,o,d) in ODpath[b1.name,b2.name]:
#         edges = edgeInPath[count, o, d ] 
#         if len(edges)<=1:
#             for (i,j) in edges:
#                 force[count,o,d,i,j] = m.addConstr( (int( (demand_dic[o,d]-0.0001)/800 ) +1 ) *x[count,o,d] <=  y[i,j] )
#                 force[count,o,d,i,j].Lazy = 3
                
# # 点ごとのカットセット制約(強化用)
# # 出る量に対する制約
# cutset1 ={} 
# for b1 in B0:
#     total_demand = 0.
#     for b2 in B0:
#         if b1 != b2:
#             total_demand += demand_dic[b1.name,b2.name]
#     cutset1[b1.name] = m.addConstr( (int( (total_demand -0.0001)/800 ) +1 )  <=  gurobi.quicksum( y[i,j] for (i,j) in y if i==b1.name ) )
#     #cutset1[b1.name].Lazy = 3 
# # 入ってくる量に対する制約
# cutset2 ={} 
# for b1 in B0:
#     total_demand = 0.
#     for b2 in B0:
#         if b1 != b2:
#             total_demand += demand_dic[b2.name,b1.name]
#     cutset2[b1.name] = m.addConstr( (int( (total_demand -0.0001)/800 ) +1 )  <=  gurobi.quicksum( y[j,i] for (j,i) in y if i==b1.name ) )
#     #cutset2[b1.name].Lazy = 3             
# #m.Params.Method = 2 # barrier method
# m.Params.Heuristics = 0.8
# m.Params.TIME_LIMIT = 3600.0*2
# m.optimize()
# print("obj fun. val. =", m.ObjVal)
# #for (i,o,d) in x:
# #    if x[i,o,d].X > 0:
# #        print(i,o,d,x[i,o,d].X)
# #for (i,j) in y:
# #    if y[i,j].X >0:
# #        print(i,j,y[i,j].X)