import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
from gurobipy import Model, quicksum, GRB, multidict, tuplelist
# from mypulp import Model, quicksum, GRB, multidict, tuplelist
from collections import defaultdict
G = nx.DiGraph()
G.add_node(0, demand=-10)
G.add_node(4, demand=10)
capacity = {(0, 1): 5, (0, 2): 8, (1, 4): 8, (2, 1): 2, (2, 3): 5, (3, 4): 6}
G.add_weighted_edges_from(
[(0, 1, 10), (0, 2, 5), (1, 4, 1), (2, 1, 3), (2, 3, 1), (3, 4, 6)]
)
for (i, j) in G.edges():
G[i][j]["capacity"] = capacity[i, j]
pos = {0: (0, 1), 1: (1, 2), 2: (1, 0), 3: (2, 0), 4: (2, 2)}
edge_labels = {}
for (i, j) in G.edges():
edge_labels[i, j] = f"{ G[i][j]['weight'] }({ G[i][j]['capacity']})"
plt.figure()
nx.draw(G, pos=pos, with_labels=True, node_size=1000)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
plt.show()
ここで考えるネットワーク上の最適化問題は最小費用流問題とよばれ， サプライ・チェインの本質とも言える「最小費用のフロー」を求める問題である． 最小費用流問題の目的は，最大流問題と同様に，ある基準を最適化する「フロー」 を求めることであるが，最短路問題と同様に，目的は費用の最小化である． つまり，最小費用流問題は，最大流問題と最短路問題の2つの特徴をあわせもった 問題と考えられる．
最短路問題のように始点と終点を決めるのではなく，ここではもう少し一般的に各点ごとに流出量を定義することにしよう． 流出量は非負とは限らず，負の流出量は流入量を表す． 流出量の総和は，ネットワーク内に「もの」が貯まらない（もしくは不足しない）ためには，$0$ になる必要がある．
最小費用流問題（minimum cost flow problem）
$n$ 個の点から構成される点集合 $V$ および $m$ 本の枝から構成される枝集合 $E$, $V$ および $E$ から成る有向グラフ $G=(V,E)$， 枝上に定義される費用関数 $c: E \rightarrow \mathbf{R}$， 枝上に定義される非負の容量関数 $u: E \rightarrow \mathbf{R}_+ \cup \{ \infty \} $， 点上に定義される流出量関数 $b: V \rightarrow \mathbf{R}$ が与えられたとき， 「実行可能フロー」で，費用の合計が最小になるものを求めよ． ただし，$\sum_{i \in V} b_i =0$ を満たすものとする．
上の問題の定義を完結させるためには，「実行可能フロー」を厳密に定義する必要がある． （最小費用流問題に対する）実行可能フロー （feasible flow）とは枝上に定義された 実数値関数 $x: E \rightarrow \mathbf{R}$ で， 以下の性質を満たすものを指す．
フロー整合条件: $$ \sum_{j: ji \in E} x_{ji} - \sum_{j: ij \in E} x_{ij} =b_i $$
容量制約と非負制約: $$ 0 \leq x_{e} \leq u_{e} \ \ \ \forall e \in E $$
networkXだと，最小費用流問題を解くためにはネットワーク単体法の関数（network_simplex）と容量スケーリング法（capacity_scaling）が準備されている． 点の属性として需要demand，枝の属性として容量capacityと費用weightを定義する． network_simplexは費用が浮動小数点数の場合には非常に遅くなるので，数値を整数に丸めてから実行する必要がある． capacity_scalingは，費用が浮動小数点数でも性能に差はないが，一般にはnetwork_simplexより遅い．
cost, flow = nx.algorithms.flow.network_simplex(G)
for i in G.nodes():
for j in flow[i]:
edge_labels[i, j] = f"{flow[i][j]}"
plt.figure()
nx.draw(G, pos=pos, with_labels=True, node_size=1000)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
plt.show()
flow = nx.algorithms.flow.max_flow_min_cost(G, s=0, t=4)
for i in G.nodes():
for j in flow[i]:
edge_labels[i, j] = f"{flow[i][j]}"
plt.figure()
nx.draw(G, pos=pos, with_labels=True, node_size=1000)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
plt.show()
輸送問題
最小費用流問題は，輸送問題(transportation problem)とよばれる古典的な問題を特殊形として含む．
いま，顧客数を $n$，工場数を $m$ とし， 顧客を $i=1,2,\ldots,n$，工場を $j=1,2,\ldots,m$ と番号で表すものとする． また，顧客の集合を $I=\{1,2,\ldots,n \}$，工場の集合を $J=\{ 1,2,\ldots,m \}$ とする． 顧客 $i$ の需要量を $d_i$，顧客 $i$ と施設 $j$ 間に $1$ 単位の需要が移動するときにかかる 輸送費用を $c_{ij}$，工場 $j$ の容量 $M_j$ とする． また， $x_{ij}$ を工場 $j$ から顧客 $i$ に輸送される量を表す連続変数する．
上の記号および変数を用いると，輸送問題は以下の線形最適化問題として定式化できる．
$$ \begin{array}{l l l} minimize & \displaystyle\sum_{i \in I} \displaystyle\sum_{j \in J} c_{ij} x_{ij} & \\ s.t. & \displaystyle\sum_{j \in J} x_{ij} =d_i & \forall i \in I \\ & \displaystyle\sum_{i \in I} x_{ij} \leq M_j & \forall j \in J \\ & x_{ij} \geq 0 & \forall i \in I, j \in J \end{array} $$目的関数は輸送費用の和の最小化であり，最初の制約は需要を満たす条件， 2番目の制約は工場の容量制約である．
networkXの最小費用流を使う場合には，需要の合計が $0$ になるように， ダミーの点を追加する．
以下では，数理最適化による定式化を示す．
I, d = multidict({1: 80, 2: 270, 3: 250, 4: 160, 5: 180})
J, M = multidict({1: 500, 2: 500, 3: 500})
c = {
(1, 1): 4,
(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,
}
model = Model("transportation")
x = {}
for i in I:
for j in J:
x[i, j] = model.addVar(vtype="C", name="x(%s,%s)" % (i, j))
model.update()
for i in I:
model.addConstr(quicksum(x[i, j] for j in J) == d[i], name="Demand(%s)" % i)
for j in J:
model.addConstr(quicksum(x[i, j] for i in I) <= M[j], name="Capacity(%s)" % j)
model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), GRB.MINIMIZE)
model.optimize()
print("Optimal value:", model.ObjVal)
EPS = 1.0e-6
for (i, j) in x:
if x[i, j].X > EPS:
print(
"sending quantity %10s from factory %3s to customer %3s" % (x[i, j].X, j, i)
)
G = nx.DiGraph()
sum_ = 0
for j in M:
sum_ -= M[j]
G.add_node(f"plant{j}", demand=-M[j])
for i in d:
sum_ += d[i]
G.add_node(f"customer{i}", demand=d[i])
# add dummy customer with demand sum_
G.add_node("dummy", demand=-sum_)
for (i, j) in c:
G.add_edge(f"plant{j}", f"customer{i}", weight=c[i, j])
for j in M:
G.add_edge(f"plant{j}", "dummy", weight=0)
cost, flow = nx.flow.network_simplex(G)
edge_labels = {}
for (i, j) in G.edges():
if flow[i][j] > 0.001:
edge_labels[i, j] = f"{flow[i][j]}"
pos = nx.bipartite_layout(G, nodes=[f"plant{j}" for j in M])
plt.figure()
nx.draw(G, pos=pos, with_labels=True, node_size=2000, node_color="y")
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
plt.show()
下限制約付き最小費用流問題
実務においてはしばしばフロー量の下限制約が付加されるときがある． ここでは，簡単なデータの変換によって下限なしの問題に帰着できることを示す．
いま，枝 $(i,j)$ 間のフロー量 $x_{ij}$ が $\ell_{ij} \leq x_{ij} \leq u_{ij}$ を満たさなければならないものとする．
このとき， $x_{ij} - \ell_{ij}$ を新たなフローを表す変数とし，以下のようにデータを変更する．
- 枝の容量 $u_{ij}$ を $u_{ij} - \ell_{ij}$ とする．
- 点 $i$ の需要量 $d_i$ に $\ell_{ij}$ を加える．
- 点 $j$ の需要量 $d_j$ から $\ell_{ij}$ を減じる（言い換えれば，$\ell_{ij}$ だけ供給する）．
すると， $(i,j)$ 上には下限 $\ell_{ij}$ の超過分だけが流れ，得られた結果に下限を加えたものが，もとの（下限付きの）問題の解になる．
例として，最初の例題の $(0,2)$ 間のフロー量の下限を $6$ としたものを考える．
G = nx.DiGraph()
G.add_node(0, demand=-10)
G.add_node(4, demand=10)
capacity = {(0, 1): 5, (0, 2): 8, (1, 4): 8, (2, 1): 2, (2, 3): 5, (3, 4): 6}
G.add_weighted_edges_from(
[(0, 1, 10), (0, 2, 5), (1, 4, 1), (2, 1, 3), (2, 3, 1), (3, 4, 6)]
)
for (i, j) in G.edges():
G[i][j]["capacity"] = capacity[i, j]
pos = {0: (0, 1), 1: (1, 2), 2: (1, 0), 3: (2, 0), 4: (2, 2)}
edge_labels = {}
for (i, j) in G.edges():
edge_labels[i, j] = f"{ G[i][j]['weight'] }({ G[i][j]['capacity']})"
plt.figure()
nx.draw(G, pos=pos, with_labels=True, node_size=1000)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
plt.show()
G.add_node(0, demand=-4) # -10+6
G.add_node(4, demand=10)
G.add_node(2, demand=-6)
cost, flow = nx.algorithms.flow.network_simplex(G)
for i in G.nodes():
for j in flow[i]:
if i == 0 and j == 2:
flow[i][j] += 6
edge_labels[i, j] = f"{flow[i][j]}"
plt.figure()
nx.draw(G, pos=pos, with_labels=True, node_size=1000)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
plt.show()
G = nx.DiGraph()
G.add_node(0, demand=-10)
G.add_node(4, demand=10)
capacity = {(0, 1): 5, (0, 2): 8, (1, 4): 8, (2, 1): 2, (2, 3): 5, (3, 4): 6}
G.add_weighted_edges_from(
[(0, 1, 10), (0, 2, 5), (1, 4, 1), (2, 1, 3), (2, 3, 1), (3, 4, 6)]
)
for (i, j) in G.edges():
G[i][j]["capacity"] = capacity[i, j]
pos = {0: (0, 1), 1: (1, 2), 2: (1, 0), 3: (2, 0), 4: (2, 2)}
edge_labels = {}
for (i, j) in G.edges():
edge_labels[i, j] = f"{ G[i][j]['weight'] }({ G[i][j]['capacity']})"
nx.draw(G, pos=pos, with_labels=True, node_size=1000)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
cost, flow = nx.algorithms.flow.network_simplex(G)
for i in G.nodes():
for j in flow[i]:
edge_labels[i, j] = f"{flow[i][j]}"
nx.draw(G, pos=pos, with_labels=True, node_size=1000)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
paths = []
flow_values = []
while 1:
inflow = defaultdict(int)
outflow = defaultdict(int)
excess = defaultdict(int)
for (i, j) in G.edges():
outflow[i] += flow[i][j]
inflow[j] += flow[i][j]
for i in range(5):
excess[i] = outflow[i] - inflow[i]
max_excess = 0
max_node = -1
min_excess = 0
min_node = -1
for i in G.nodes():
if max_excess < excess[i]:
max_node = i
max_excess = excess[i]
if min_excess > excess[i]:
min_node = i
min_excess = excess[i]
if max_excess == 0:
break
print("Find a path from ", max_node, " to ", min_node)
D = nx.Graph()
for (i, j) in G.edges():
if flow[i][j] > 0:
D.add_edge(i, j)
path = nx.shortest_path(D, source=max_node, target=min_node)
print("Shortest path= ", path)
paths.append(path)
min_flow = 9999999
i = path[0]
for j in path[1:]:
min_flow = min(min_flow, flow[i][j])
i = j
print("min_flow=", min_flow)
flow_values.append(min_flow)
i = path[0]
for j in path[1:]:
flow[i][j] -= min_flow
i = j
print("new flow=", flow)
print(paths, flow_values)