Steiner木問題に対する定式化とソルバー
import networkx as nx
from gurobipy import Model, quicksum, GRB
# from mypulp import Model, quicksum, GRB
import matplotlib.pyplot as plt
import random
Steiner木問題に対する定式化
Steiner木問題は，最小木問題の拡張であり，$NP$-困難である． 最小木問題が，すべての点を通る全域木を求めるのに対して，Stiner木問題では，与えられた点（ターミナルと呼ばれる）の部分集合を繋ぐ木を求める． ターミナル以外の点は使っても良いし，使わなくても良い．これが問題を難しくする．定義は以下のようになる．
無向グラフ $G(V,E)$ と枝上の費用関数 $c: E \rightarrow \mathbf{R}$ が与えられたとき， 点の部分集合 $VT \subseteq V$ と枝の部分集合 $ET \subseteq E$ から構成される連結グラフで
$$ \sum_{e \in ET} c_e $$を最小にするものを求める．ただし，特定の点 $0$ は必ず $VT$ に含まれるものと仮定する．
最小木問題に対する多品種流定式化を用いると，Steiner木問題を混合整数最適化問題として定式化できる．
枝 $e \in E$ がSteiner木に含まれるとき $1$，それ以外のとき $0$ を表す $0$-$1$ 変数 $x_e$ を導入する． 点$0$ から点$k$ へ流れるフローが，枝$(i,j)$ 上を $i$ から $j$ の向きで通過する量を表す実数変数を $f_{ij}^k$ と したとき，多品種流定式化は，以下のようになる．
$$ \begin{array}{l l l } minimize & \sum_{e \in E} c_e x_e & \\ s.t. & \sum_{(j,i) \in \delta (\{i\})} f_{ji}^k -\sum_{(i,j) \in \delta (\{i\})} f_{ij}^k = \left\{ \begin{array}{l } -1 \\ 0 \\ 1 \\ \end{array} \right. & \begin{array}{l } i=0, \forall k \in VT \setminus \{0\} \\ \forall i \in V \setminus \{0,k\}, k \in V \setminus \{0\} \\ i=k, \forall k \in VT \setminus \{0\} \\ \end{array} \\ & f_{ij}^k +f_{ji}^{k'} \leq x_e & \forall k,k' \in V \setminus \{0\}, e=(i,j) \in E \\ & x_e \in \{0,1\} & \forall e \in E \\ & f_{ij}^k \geq 0, f_{ji}^k \geq 0 & \forall k \in V \setminus \{0\}, e=(i,j) \in E \end{array} $$ランダムに枝上の費用（重み）を設定した格子グラフに対して，四隅の点をターミナルとしたSteiner木問題を解いてみる．
m, n = 5, 5
lb, ub = 1, 20
random.seed(1)
G = nx.grid_2d_graph(m, n)
pos = {(i, j): (i, j) for (i, j) in G.nodes()}
for (i, j) in G.edges():
G[i][j]["weight"] = random.randint(lb, ub)
Terminal = [(0, 0), (m - 1, 0), (0, n - 1), (m - 1, n - 1)]
root = Terminal[0]
model = Model()
x, f = {}, {}
for (i, j) in G.edges():
x[i, j] = model.addVar(vtype="B", name=f"x[{i},{j}]")
for k in Terminal[1:]:
f[i, j, k] = model.addVar(vtype="C", name=f"f[{i},{j},{k}]")
f[j, i, k] = model.addVar(vtype="C", name=f"f[{j},{i},{k}]")
model.update()
for k in Terminal[1:]:
for i in G.nodes():
if i == root:
model.addConstr(
quicksum(f[i, j, k] for j in G[i]) - quicksum(f[j, i, k] for j in G[i])
== 1
)
elif i == k:
model.addConstr(
quicksum(f[i, j, k] for j in G[i]) - quicksum(f[j, i, k] for j in G[i])
== -1
)
else:
model.addConstr(
quicksum(f[i, j, k] for j in G[i]) - quicksum(f[j, i, k] for j in G[i])
== 0
)
for (i, j) in G.edges():
for k in Terminal[1:]:
model.addConstr(f[i, j, k] + f[j, i, k] <= x[i, j])
model.setObjective(
quicksum(G[i][j]["weight"] * x[i, j] for (i, j) in G.edges()), GRB.MINIMIZE
)
model.optimize()
edges = []
print(model.ObjVal)
for (i, j) in G.edges():
if x[i, j].X > 0.01:
edges.append((i, j))
plt.figure()
nx.draw(G, pos=pos, node_size=100)
edge_labels = {}
for (i, j) in G.edges():
edge_labels[i, j] = f"{ G[i][j]['weight'] }"
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
nx.draw(G, pos=pos, width=5, edgelist=edges, edge_color="orange")
plt.show()
from networkx.algorithms.approximation.steinertree import steiner_tree
steiner = steiner_tree(G, terminal_nodes=Terminal)
total = 0
for (i, j) in steiner.edges():
total += G[i][j]["weight"]
print(total)
plt.figure()
nx.draw(G, pos=pos, with_labels=False, node_size=100)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
nx.draw(G, pos=pos, width=5, edgelist=steiner.edges(), edge_color="orange", node_size=0)
plt.show()
賞金収集有向Steiner木問題
ここでは，応用（枝の向きによって費用が異なる場合もあるので）を意識して有向グラフ上で考える． 有向グラフ $D(N,A)$， 点上の利得（賞金）関数 $p: N \rightarrow \mathbf{R}$， 枝上の費用関数 $c: S \rightarrow \mathbf{R}$ が与えられたとき， 点の部分集合 $NT \subseteq N$ と枝の部分集合 $AT \subseteq A$ から構成される連結グラフで
$$ \sum_{i \in NT} p_i - \sum_{(i,j) \in AT} c_{ij} $$を最大にするものを求める．ただし，特定の点 $0$ は必ず $NT$ に含まれるものと仮定する．
最小木問題に対する多品種流定式化を用いると，賞金収集有向Steiner木問題を混合整数最適化問題として定式化できる．
有向グラフ $D=(N,A)$ 上で定式化する．有向枝 $(i,j) \in A$ がSteiner木に含まれるとき $1$，それ以外のとき $0$ を表す $0$-$1$ 変数 $x_{ij}$ を導入する． さらに，点 $i$ がSteiner木に含まれるとき $1$ の$0$-$1$ 変数 $y_{i}$ を準備する． 点$0$ から点$k$ へ流れるフローが，枝$(i,j)$ 上を通過する量を表す実数変数を $f_{ij}^k$ と したとき，賞金収集有向Steiner木問題に対する多品種流定式化は，以下のようになる．
$$ \begin{array}{l l l } maximize & \sum_{i \in N \setminus \{0\}} p_i y_i - \sum_{(i,j) \in A} c_{ij} x_{ij} & \\ s.t. & \sum_{j} x_{ij} = y_i & \forall i \in N \setminus \{0\} \\ & \sum_{j} f_{ji}^k -\sum_{j} f_{ij}^k = \left\{ \begin{array}{l } -y_k \\ 0 \\ y_k \\ \end{array} \right. & \begin{array}{l } i=0, \forall k \in N \setminus \{0\} \\ \forall i \in N \setminus \{0,k\}, k \in V \setminus \{0\} \\ i=k, \forall k \in N \setminus \{0\} \\ \end{array} \\ & f_{ij}^k \leq x_{ij} & \forall k \in N \setminus \{0\},(i,j) \in A \\ & x_{ij} \in \{0,1\} & \forall (i,j) \in A \\ & y_i \in \{0,1\} & \forall i \in N \\ & f_{ij}^k \geq 0 & \forall k \in V \setminus \{0\}, (i,j) \in A \\ \end{array} $$m, n = 5, 5
lb, ub = 1, 20
plb, pub = 1, 20
root = (0, 0)
G = nx.grid_2d_graph(m, n)
pos = {(i, j): (i, j) for (i, j) in G.nodes()}
for (i, j) in G.edges():
G[i][j]["weight"] = random.randint(lb, ub)
for i in G.nodes():
if i != root:
G.nodes[i]["prize"] = random.randint(plb, pub)
D = G.to_directed()
model = Model()
x, y, f = {}, {}, {}
for (i, j) in D.edges():
x[i, j] = model.addVar(vtype="B", name=f"x[{i},{j}]")
for k in D.nodes():
if k != root:
f[i, j, k] = model.addVar(vtype="C", name=f"f[{i},{j},{k}]")
for i in D.nodes():
if i != root:
y[i] = model.addVar(vtype="B", name=f"y[{i}]")
model.update()
for i in D.nodes():
if i != root:
model.addConstr(quicksum(x[j, i] for j in D[i]) == y[i])
for k in D.nodes():
if k == root:
continue
for i in D.nodes():
if i == root:
model.addConstr(
quicksum(f[i, j, k] for j in D[i]) - quicksum(f[j, i, k] for j in D[i])
== y[k]
)
elif i == k:
model.addConstr(
quicksum(f[i, j, k] for j in D[i]) - quicksum(f[j, i, k] for j in D[i])
== -y[k]
)
else:
model.addConstr(
quicksum(f[i, j, k] for j in D[i]) - quicksum(f[j, i, k] for j in D[i])
== 0
)
for (i, j) in D.edges():
for k in D.nodes:
if k != root:
model.addConstr(f[i, j, k] <= x[i, j])
model.setObjective(
quicksum(D.nodes[i]["prize"] * y[i] for i in D.nodes() if i != root)
- quicksum(D[i][j]["weight"] * x[i, j] for (i, j) in D.edges()),
GRB.MAXIMIZE,
)
model.optimize()
edges = []
print(model.ObjVal)
for (i, j) in D.edges():
if x[i, j].X > 0.01:
edges.append((i, j))
plt.figure()
nx.draw(G, pos=pos, with_labels=False, node_size=1000, node_color="yellow")
edge_labels = {}
for (i, j) in G.edges():
edge_labels[i, j] = f"{ G[i][j]['weight'] }"
node_labels = {}
G.nodes[root]["prize"] = 0
for i in G.nodes():
node_labels[i] = G.nodes[i]["prize"]
nx.draw_networkx_labels(G, pos, labels=node_labels)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
nx.draw(D, pos=pos, width=5, edgelist=edges, edge_color="red", node_size=0)
plt.show()
