容量制約付きの最小木に対する定式化とソルバーによる求解
import networkx as nx
import matplotlib.pyplot as plt
from gurobipy import Model, quicksum, GRB
#from mypulp import Model, quicksum, GRB
from collections import OrderedDict, defaultdict
容量制約付き有向最小木問題
最小木問題において、枝に向きがある場合には、最小有向木問題 (minimum arborescence problem) になる。
この問題は多項式時間で解くことができるが、これに枝の容量制約を付加したものが、容量制約付き有向最小木問題である。 ベンチマーク問題例は、
http://people.brunel.ac.uk/~mastjjb/jeb/orlib/capmstinfo.html
から入手できる。
有向グラフ $D =(N,A)$ と特定の点（根と呼ぶ； $0$ と仮定）が与えられている。他の点は1単位の需要をもち、それは根から選択された枝を経由して運ばなければならない。枝を通過可能なフロー量の上限（容量） $Q$ が与えられているとき、枝上の総費用を最小にする有向木を求めよ。
定式化
- $D =(N,A)$: 有向グラフ
- $n$ : 点の数
- $0$ : グラフの根を表す点
- $c_{ij}$: 枝 $(i,j)$ の費用
- $Q$: 枝の容量
- $x_{ij}$: 枝を使用するとき1、それ以外のとき0
- $y_{ij}$: 枝上を流れるフロー量
with open("../data/capmst/capmst1.txt") as f:
lines = f.readlines()
n_line = 0
for iter_ in range(20):
prob = lines[n_line]
n_line += 1
n = int(lines[n_line]) # number of nodes
G = nx.DiGraph()
n_line += 1
for i in range(n + 1):
array = []
while 1:
line = lines[n_line].split()
array.extend(line)
if len(array) >= n + 1:
break
n_line += 1
for j, cost in enumerate(array):
if i != j:
G.add_edge(i, j, weight=int(cost))
n_line += 1
n_line += 1
nx.write_weighted_edgelist(G, "../data/capmst/" + prob[:-6].strip() + ".txt")
例として，tc40-1.txtの問題を読み込んで解いてみる．容量Qは別途設定するが，以下のコードのコメントにあるように，点数nによって変える必要がある．
fname = "../data/capmst/tc40-1.txt"
G = nx.read_weighted_edgelist(fname, create_using=nx.DiGraph, nodetype=int)
n = len(G)
# Q=3, 5, or 10 for n=40
# Q=5, 10, or 20 for n=80
Q = 10
print("n=", n, "Q=", Q)
model = Model()
x, y = {}, {}
for (i, j) in G.edges():
if i != j and j != n - 1: # ルートであるn-1に入る枝はない
x[i, j] = model.addVar(vtype="B", name=f"x[{i},{j}]")
y[i, j] = model.addVar(vtype="C", name=f"y[{i},{j}]")
model.update()
model.addConstr(quicksum(x[i, j] for (i, j) in x) == n - 1)
for j in range(0, n - 1): # ルート以外の点は入次数が1
model.addConstr(
quicksum(x[i, j] for i in range(n) if i != j) == 1, name=f"in_degree[{j}]"
)
for j in range(0, n - 1): # ルート以外の点に対するフロー保存式
model.addConstr(
quicksum(y[i, j] for i in range(n) if i != j)
- quicksum(y[j, i] for i in range(0, n - 1) if i != j)
== 1,
name=f"flow_conservarion[{j}]",
)
for i in range(n):
for j in range(0, n - 1):
if i != j:
model.addConstr(x[i, j] <= y[i, j], name=f"lower_bound[{i},{j}]")
if i == n - 1: # ルートに入るフロー量はQ以下
model.addConstr(y[i, j] <= Q * x[i, j], name=f"upper_bound[{i},{j}]")
else: # ルート以外の点に入るフロー量はQ-1以下
model.addConstr(
y[i, j] <= (Q - 1) * x[i, j], name=f"upper_bound[{i},{j}]"
)
model.setObjective(quicksum(G[i][j]["weight"] * x[i, j] for (i, j) in x), GRB.MINIMIZE)
model.optimize()
SolGraph = nx.DiGraph()
for (i, j) in x:
if x[i, j].X > 0.1:
SolGraph.add_edge(i, j)
nx.draw(SolGraph)
