はじめに
ここでは，ロジスティクス・ネットワーク設計問題に対する包括的モデルを示す． モデルの目的は，単位期間（通常は年）ベースのデータをもとに， ロジスティクス・ネットワークの形状を決めることにある． モデルを求解することによって得られるのは，倉庫の設置の是非， 地点間別の各製品の単位期間内の総輸送量， 工場別の各製品の単位期間内の総生産量である．
ロジスティクス・ネットワーク設計モデルでは， 以下の意思決定項目に対する最適化を同時に行う．
- 各製品をどの工場でどれだけ生産するか？
- 各製品をどの倉庫（配送センターや中継拠点の総称）で保管するか？
- 各顧客群の各製品の需要を，どの地点から運ぶか？
- 倉庫をどこに新設するか？（または移転，閉鎖するか？）
- 複数の倉庫の候補地点からどれを選択するか？
このような意思決定は，単に現状のロジスティクス・ネットワーク の見直しを行う場合だけでなく，以下の状況においても有効に用いられる．
- 吸収・合併の後のロジスティクス・ネットワークの再編成
- 新製品投入時の意思決定（どこで製造して，どのような流通チャネルで流すか）
- ロジスティクスにおける戦略的提携（ロジスティクス・パートナーシップ）
顧客データ
品が最終的に消費される場所を顧客（customer）とよぶ． 顧客は，対象となるシステムに依存して，小売店であったり個々の家庭であったりする． モデルとして扱う場合には，同じ施設からサービスされると考えられる顧客の集まりを 1つの顧客に集約することができる． ここでは，個々の顧客を同じサプライ・チェインでサービス を行っても差し支えない程度に集約したものを顧客群（customer group）とよぶこともあるが，以下では単に顧客と呼ぶ．
通常は，需要の大きい顧客や特殊なサービスを要求する顧客は1つの顧客として扱い， 残りの顧客を地区ごとに集約して扱うことが多い． これは，同じ地区に属する顧客への配送は，（大口や特殊な顧客以外は） 同一のデポから行われる場合が多いことに起因する． 個々の顧客の数は数万から数十万におよぶこともあるが， 集約後の顧客群の数は，多くても数百になるようにする．
顧客群に集約するための簡便法として郵便番号を用いる方法もあるが， できれば1つの顧客群に含まれる顧客数および需要量が，なるべくバランスするように することが望ましい． ここでは，幾つかのクラスタリング手法 (Weizfeld法， $k$-means法, 階層的クラスタリング法, $k$-メディアン法） を用いる．
例題で用いる顧客データの列は， 名前 (name) と緯度・経度 (lat, lon; latitude, longitudeの略) である．
cust_df = pd.read_csv(folder + "Cust.csv", index_col=0)
cust_df.head()
製品データ
ロジスティクス・ネットワーク内を流れるものを製品（product）とよぶ． 製品の種類は，対象とするロジスティクス・システムにもよるが，通常は膨大なものになるが， ストラテジックレベルの意思決定モデルを扱う場合には，通常同じ特性をもつ製品を集約して扱う． 製品を集約したものを，個々の製品と区別するために製品（product group）とよぶこともある． 通常，製品の数が数十程度になるように集約する． また，ABC分析をして，影響の少ない製品を除いて最適化を行うこともある．
製品データは，以下の列をもつが，ロジスティクス・ネットワーク設計モデルで使うのは，名称 (name)と重量 (weight)だけである．
- name: 製品の名称
- weight: 製品の重量 (単位はkg)
- volume: 製品の容量 (単位は 𝑚3 ）
- cust_value: 顧客上での製品の価値
- dc_value: 倉庫上での製品の価値
- plnt_value: 工場における製品の価値
- fixed_cost: 製品を生産する際の段取り（固定）費用
prod_df = pd.read_csv(folder + "Prod.csv", index_col=0)
prod_df.head()
倉庫データ
最終需要地点および供給地点以外の場所で，かつ製品の製造・加工を行わない地点を総称して 倉庫（warehouse）とよぶことにする． 実際には，倉庫は流通センター（DC: Distribution Centerの略），デポ，クロスドッキング地点，港，空港，ハブなど様々な形態をとる．
倉庫の役割と機能には，以下のものが考えられる．
調達（購入）の規模の経済性のため． 大量購入に伴う割引が行われるとき，大量に購入した製品を一時的に保管しておく際に，倉庫の保管機能が用いられる．
混載による規模の経済性のため．複数の小売店（スーパーなど）で消費される製品を， 直接製造業者（メーカー）から輸送すると，小ロットの輸送のため費用が増大する． メーカーから車立て輸送を行うためには，複数の小売店の需要をあわせる必要があるが， そのための中継地点として倉庫の仕分け機能が用いられる． メーカーから運ばれた複数の製品は，倉庫で方面別に仕分けされて，小売店に運ばれる． この際，複数の製品を混載して輸送するので，やはり小ロット輸送による費用の増大を避けることができる． この仕分け機能に特化した倉庫をクロスドッキング地点（cross docking point）とよぶ． クロスドッキング地点は，保管機能をもたず，メーカーから小売店が発注した量と同じ量を受け取り， 即座に仕分けして小売店に輸送する．
リスク共同管理効果を得るため． 複数の小売店の在庫を，集約して倉庫に保管することによって，安全在庫量を減らすことができる．これをリスク共同管理（risk pooling）効果とよぶ． この場合には，倉庫は，リスク共同管理によってサプライ・チェイン全体での在庫を減らすために用いられる．
顧客へのサービスレベル向上のため． たとえば，コンビニエンスストアの配送センターは，弁当が約$3$時間以内で配達可能なように設置される． これによって，$1$日に数回の弁当の配送が可能になり，顧客に出来たての弁当を届けることが可能になる． このように，倉庫は小売店（コンビニエンスストア）へのリード時間を減少させるために用いられる．
単一ソース条件． 顧客（特に小売店）に対して，複数の倉庫からの配送を行うことを許さないという条件が付加されることがある． これは，荷捌き場の制限や，荷物の受け入れを複数回行う手間を省くためであるが，そのためには，複数のメーカーの製品を倉庫で積み替えて 運ぶことが必要になる．
製品に付加価値をつけるため． 最終製品への分化を，なるべくサプライ・チェインの下流で行うことによって， 顧客サービスを向上させ，安全在庫を減少させることができる． これを遅延差別化（delayed differentiation, postponement）とよぶ． 遅延差別化の方法として， 工場で最終製品まで製造してまうのではなく，顧客に近い倉庫で，製造の最終工程を行うことがあげられる． たとえば，値札やラベルを付けたり，簡単な加工を行うことが代表例である．
倉庫データは，倉庫の配置可能地点を表し，以下の列をもつ．
- name: 倉庫名称
- lb: 容量下限
- ub: 容量上限
- fc: 固定費用
- vc: 変動費用
- lat: 緯度
- lon: 経度
dc_df = pd.read_csv(folder + "DC.csv", index_col=0)
dc_df.head()
plnt_df = pd.read_csv(folder + "Plnt.csv", index_col=0)
plnt_df.head()
demand_df = pd.read_csv(folder + "demand.csv", index_col=0)
demand_df.head()
plnt_prod_df = pd.read_csv(folder + "Plnt-Prod.csv", index_col=0)
plnt_prod_df.head()
連続施設配置問題
平面上に施設が配置可能と仮定した施設問題に対する近似解法を考える。
施設数が1つの場合には、Weiszfeld法と呼ばれる反復法で近似解が得られることが知られている。 複数施設に対しては、k-means法と似た反復解法が考えられる。ここでは、重み付き複数施設配置問題に対するWeiszfeld法を構築する。
以下では、1つの施設を選択するためのWeiszfeld法について解説する。
各顧客 $i$ は平面上に分布しているものとし，その座標を $(X_i,Y_i)$ とする． 顧客は重み $w_i$ をもち，目的関数は施設と顧客の間の距離に重みを乗じたものの和とする． 顧客と施設間の距離は $\ell_p$ ノルム（$1 \leq p$）を用いて計算されるものとする． 東京都市圏の道路ネットワークでは $p \approx 1.8901$ と推定されており， 実務的には Euclidノルム（$p=2$）に迂回係数（$1.3$ 程度） を乗じたものを使う場合が多い．
顧客の集合 $I$ に対して，単一の倉庫の配置地点 $(X,Y)$ を決定する問題は， $$ f(X,Y) = \sum_{i \in I} w_i \left\{ (X-X_i)^p +(Y-Y_i)^p \right\}^{1/p} $$ を最小にする $(X,Y) \in \mathbf{R}^2$ を求める問題になる． これは，以下の Weiszfeld法を拡張した解法を用いることによって容易に求解できる．
上の関数は凸関数であるので， $(X,Y)$ が最適解である必要十分条件は，$(X,Y)$ が $$ \frac{\partial f(X,Y)}{\partial X} =0 $$ および $$ \frac{\partial f(X,Y)}{\partial Y} =0 $$ を満たすことである． $$ \frac{\partial f(X,Y)}{\partial X} = \sum_{i \in I} w_i (X-X_i) \frac{ |X-X_i|^{p-2}}{\left\{ (X-X_i)^p +(Y-Y_i)^p \right\}^{(p-1)/p} } =0 $$ は陽的に解くことができないが， 以下の $X$ に対する繰り返し式を示唆している． $$ X^{(q+1)}= \frac{ \sum_{i \in I} w_i |X^{(q)}-X_i|^{p-2} X_i/ \left\{ (X^{(q)}-X_i)^p +(Y^{(q)}-Y_i)^p \right\}^{(p-1)/p} }{ \sum_{i \in I} w_i |X^{(q)}-X_i|^{p-2}/ \left\{ (X^{(q)}-X_i)^p +(Y^{(q)}-Y_i)^p \right\}^{(p-1)/p} } $$
$Y$ についても同様の式を導くことができる．
これらの式を利用した解法は，一般に不動点アルゴリズムとよばれ， $1 \leq p \leq 2$ のとき大域的最適解に収束する．
Weiszfeld法を実行する関数 weiszfeld_numpy と weiszfeld
施設数が1つの場合には、Weiszfeld法と呼ばれる反復法で近似解が得られることが知られている。 複数施設に対しては、$k$-means法と似た反復解法が考えられる。ここでは、重み付き複数施設配置問題に対するWeiszfeld法を構築する。
移動距離は大圏距離distanceを用いる。これはgeopyから読み込んであると仮定する。 例題においては，顧客の重みは、年間の製品需要量に製品の重量を乗じたものとする。
引数：
- cust_df : 顧客データフレーム；緯度経度情報を利用する。
- weight : 顧客の重み；重み付きの大圏距離の和を最小化するものとする。
- num_of_facilities : 平面上に配置する施設の数を指定する。
- epsilon : 誤差上限
- max_iter : 反復回数の上限
- seed: 乱数の種
- X0 : 施設の緯度のリストの初期値
- Y0 : 施設の経度のリストの初期値
返値：
- X : 施設の緯度のリスト
- Y : 施設の経度のリスト
- partition : 顧客の割り当てられた施設の番号を格納したリスト
- cost: 費用（重み付きの距離の総和）
weiszfeld_numpyはNumPyを用いた実装であり， weiszfeldは通常の実装である．
cust_df = pd.read_csv(folder + "Cust.csv", index_col=0) #read customer data for using Lat/Lng
cust_df.head()
prod_df = pd.read_csv(folder + "Prod.csv", index_col=0) #read product data (for using weight)
prod_df.head()
total_demand = pd.read_csv(folder + "total_demand.csv")
total_demand.head()
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.itertuples():
weight_of_cust[row.cust] += weight_of_prod[row.prod]*row.demand
#weight = [weight_of_cust[row.name] for row in cust_df.itertuples()]
weight = [weight_of_cust[i] for i in weight_of_cust]
# 関数呼び出し
X, Y, partition, cost = weiszfeld(cust_df, weight, num_of_facilities = 10, epsilon=0.0001, max_iter = 10, seed=1)
print(cost)
# X, Y, partition, cost = weiszfeld(cust_df, weight, num_of_facilities = 10, epsilon=0.0001, max_iter = 10, seed=1, X0=X, Y0=Y)
# print(cost)
#X, Y, partition, cost = weiszfeld_numpy(cust_df, weight, num_of_facilities = 10, epsilon=0.0001, max_iter = 10, seed=1)
#print(cost)
# X, Y, partition, cost = weiszfeld_numpy(cust_df, weight, num_of_facilities = 10, epsilon=0.0001, max_iter = 10, seed=123, X0=X, Y0=Y)
# print(cost)
反復Weiszfeld法を実行する関数 repeated_weiszfeld
Weiszfeld法を一定回数初期解を変えながら実行し，最も良い結果を返す関数
引数：
- cust_df : 顧客データフレーム；緯度経度情報を利用する。
- weight : 顧客の重み；重み付きの大圏距離の和を最小化するものとする。
- num_of_facilities : 平面上に配置する施設の数を指定する。
- epsilon : 誤差上限
- max_iter : 反復回数の上限
- numpy: NumPyを用いる場合にTrue
- seed: 乱数の種
返値：
- X : 施設の緯度のリスト
- Y : 施設の経度のリスト
- partition : 顧客の割り当てられた施設の番号を格納したリスト
- cost: 費用（重み付きの距離の総和）
%time
X, Y, partition, cost = repeated_weiszfeld(cust_df, weight,
num_of_facilities = 10, epsilon=0.0001, max_iter = 30, numpy=False, seed=2)
print(cost)
demand_df = pd.read_csv(folder + "demand.csv")
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)
weight = demand_cust_df.demand.values
#total_demand_df = make_total_demand(demand_df)
total_demand_df.head()
%%time
X, Y, partition, cost = repeated_weiszfeld(cust_df, weight, num_of_facilities = 8, epsilon=0.0001, max_iter = 3, numpy=True, seed=1)
print(cost)
durations, distances, node_df = compute_durations(cust_df, host=host)
X, Y, partition, cost = hierarchical_clusterning(cust_df, weight, durations, num_of_facilities = 10, linkage="complete")
print(cost)
def transportation(C, capacity):
M = capacity
n,m = C.shape
C = np.ceil(C) # for network simplex
G =nx.DiGraph()
sum_ = 0
for j in range(m):
sum_ -= M
G.add_node(f"plant{j}", demand=-M)
for i in range(n):
sum_ += 1
G.add_node(i, demand= 1)
#add dummy customer with demand sum_
G.add_node("dummy", demand= -sum_)
for i in range(n):
for j in range(m):
G.add_edge(f"plant{j}", i, weight=C[i,j])
for j in range(m):
G.add_edge(f"plant{j}", "dummy", weight = 0)
cost, flow = nx.flow.network_simplex(G)
return cost, flow
c = {(1,1):4.5, (1,2):6, (1,3):9,
(2,1):5, (2,2):4, (2,3):7,
(3,1):6, (3,2):3, (3,3):4,
(4,1):8, (4,2):5, (4,3):3,
(5,1):10, (5,2):8, (5,3):4,
}
n = 5
m = 3
C = np.zeros( (n,m) )
for (i,j) in c:
C[i-1,j-1] = c[i,j]
cost, flow = transportation(C, capacity=3)
def find_median(C, flow, n, m):
total_cost = 0
for j,f in enumerate(flow):
if j>=m:
break
cluster =[]
for i in range(n):
if flow[f][i]==1:
cluster.append(i)
if len(cluster)>1:
#find the best location for the cluster
nodes = np.array(cluster)
subC = C[np.ix_(nodes,nodes)]
best_location = nodes[subC.sum(axis=0).argmin()]
cost = subC.sum(axis=0).min()
total_cost+=cost
return total_cost
c = np.array(durations)
C = c*weight.reshape((n,1))
n = len(C)
k = 10
m = k
M = int(len(cust_df)/k*1.2) #容量制約
idx = np.array([0,1,2])
idx = np.array([ 0,72,71,70,69,68,67,66,65,64])
cost, flow = transportation(C[:,idx], capacity=M)
find_median(C, flow, n, m)
k-メディアン問題を解くための関数 solve_k_median
引数：
- cust_df : 顧客データフレーム；緯度経度情報を利用する。
- weight : 顧客の需要量
- cost: 移動費用（地図で求めたdurationsを入力する．）
- num_of_facilities : 平面上に配置する施設の数を指定する。
- max_iter =100: 最大反復数
- max_lr =0.01: 最大学習率（ステップサイズ）
- moms=(0.85,0.95): 慣性項の下限と上限
- convergence=1e-5: 収束判定のためのパラメータ
- lr_find = False: 学習率探索を行う場合 True
- adam = False: Adam法によるパラメータの調整を行う場合 True，上界と下界の差をもとに決める場合 False
返値：
- X : 施設の緯度のリスト
- Y : 施設の経度のリスト
- partition : 顧客の割り当てられた施設の番号を格納したリスト
- best_ub: 費み付きの費用の総和
- lb_list: 下界値の変化を入れたリスト
- ub_list: 上界値の変化を入れたリスト
- phi_list: 学習率探索の場合の学習率の変化を入れたリスト
cust_df = generate_cust(num_locations = 100, random_seed = 3, prefecture = None , no_island=True)
weight=np.random.random(len(cust_df))*10
durations, distances, node_df = compute_durations(cust_df, host=host)
%%time
X, Y, partition, best_ub, lb_list, ub_list, phi_list = solve_k_median(cust_df, weight, durations, num_of_facilities=10,
max_iter=10000, max_lr=20., moms=(0.85,0.95), convergence=1e-5, lr_find = False, adam=True, capacity = None)
fig = plot_k_median_lr_find(phi_list, lb_list)
plotly.offline.plot(fig);
fig = plot_k_median(lb_list, ub_list)
plotly.offline.plot(fig);
