import random
import math
from collections import defaultdict
from gurobipy import Model, quicksum, GRB
# from mypulp import Model, quicksum, GRB
import networkx as nx
import matplotlib.pyplot as plt
G = nx.grid_2d_graph(3, 4)
nx.draw_spectral(
G, with_labels=True, node_color="w", node_size=1000, edge_color="g", width=10
)
plt.show()
中国郵便配達人問題
枝巡回問題において，積載量制約がなく，枝が無向かつ連結な場合には，中国郵便配達人問題（Chinese postman problem）とよばれる． この問題は，中国人の数学者 Mei-Gu Kwan（管 梅谷）によってはじめて考えられた問題である．
枝上に費用をもつ一般の無向グラフを，Euler閉路をもつように最小費用で枝を追加することを考える． 無向グラフがEuler閉路もつためには，点の次数が偶数である必要があった． 次数が奇数の点集合に対して，点間の最短路の費用をもつ枝をはったグラフを作り，それに対して最小費用のマッチングを求める． マッチングに対応するパスを，元のグラフに追加することによって，点の次数はすべて偶数になるので，Euler閉路は簡単に求めることができる．
NetworkXには，上のアルゴリズムがある．Euler閉路をもつように枝が追加されたグラフは，eulerize関数で計算できる．
NewG = nx.eulerize(G)
print("Eulerian?", nx.is_eulerian(NewG))
D = nx.MultiDiGraph()
for e in nx.eulerian_circuit(NewG):
D.add_edge(e[0], e[1])
print(e[0],"=>" ,end="")
print(e[1])
nx.draw_spectral(D, with_labels=True, node_size=1000, edge_color="g", width=1)
空輸送最小化問題
いま，荷が発地から着地まで積み替えることなしに運ばれるものとする．．
直送される荷を，$1$ 台の運搬車によって運ぶとき，運搬車が余分な距離を走らないようにするためには， 荷物を運ばないで移動している，いわゆる空輸送を最小化すれば良い．
これは，Euler閉路を求めることに他ならない．有向グラフが（有向の）Euler閉路をもつためには， グラフの各点の入次数（点に入ってくる枝の本数）と出次数（点から出て行く枝の本数）が一致していれば良い． すなわち，空輸送の最小化は，なるべく少ない（総費用が小さい）枝を追加してグラフの入次数と出次数が一致するようにする問題に帰着される．
これは，出次数から入次数を減じた量を需要量（負の量は供給量）とした輸送問題を解くことに他ならない． 輸送問題は，最小費用流問題に帰着できるので，networkXで容易に求解できる．
例として，2つのデポから10の顧客への直送の問題をランダムに生成する．
random.seed(1)
def distance(x1, y1, x2, y2):
"""distance: euclidean distance between (x1,y1) and (x2,y2)"""
return int(math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2))
# ランダムに 顧客を配置
n = 10
x = dict([(i, 100 * random.random()) for i in range(n)])
y = dict([(i, 100 * random.random()) for i in range(n)])
# 両端の中央にデポを配置
n_depots = 2
x.update({n: 0, n + 1: 100})
y.update({n: 50, n + 1: 50})
N = n + n_depots
pos = {}
# 複数の輸送要求がある場合には MultiDiGraph を用いる
G = nx.DiGraph()
for i in range(N):
G.add_node(i)
pos[i] = (x[i], y[i])
c = {}
for i in range(N):
for j in range(N):
if j != i:
c[i, j] = distance(x[i], y[i], x[j], y[j])
# 遠いデポからのみ輸送を定義
for i in range(n):
max_ = -1
max_depot = -1
for j in range(n, N):
if max_ < c[i, j]:
max_ = c[i, j]
max_depot = j
G.add_edge(max_depot, i)
plt.figure()
nx.draw(G, pos=pos, node_size=300, with_labels=True, node_color="yellow")
plt.show()
次に，各点の出次数と入次数の差を計算し，networkXのネットワーク単体法 network_simplex で最小費用流問題を解き，得られた枝を追加する． これによって，各点の入次数と出次数が一致し，Euler閉路を求めることができる．
B = nx.DiGraph()
for i in range(N):
B.add_node(i, demand=G.out_degree(i) - G.in_degree(i))
for i in range(n):
for j in range(n, N):
B.add_edge(i, j, weight=c[i, j])
cost, flow = nx.network_simplex(B)
G2 = G.copy()
for i in flow:
for j in flow[i]:
if flow[i][j] == 1: # 複数の輸送がある場合には，並列枝を追加
G2.add_edge(i, j)
nx.draw(G2, pos=pos, node_size=300, with_labels=True, node_color="yellow")
for i in nx.eulerian_circuit(G2, source=n):
print(i)
実際には，1台の運搬車ですべての荷を運ぶことは時間的に難しい． 以下では，上で生成したグラフの単純閉路をすべて列挙し，すべての荷が運べるように最小費用の閉路を選ぶことによって，現実的な解を得る． これは，集合分割アプローチを適用したことに他ならない．
Cycles = defaultdict(list) # 枝をキー，枝を含む閉路のリストを値とした辞書
Cost = defaultdict(int) # 閉路の費用
for cycle in nx.simple_cycles(G2):
total = 0
for idx, i in enumerate(cycle[:-1]):
edge = i, cycle[idx + 1]
total += c[edge]
if edge in G.edges():
Cycles[edge].append(cycle)
edge = cycle[-1], cycle[0]
total += c[edge]
if edge in G.edges():
Cycles[edge].append(cycle)
Cost[tuple(cycle)] = total
Cost
model = Model()
x = {}
for cycle in Cost:
x[cycle] = model.addVar(vtype="B", name=f"Cost[{cycle}]")
model.update()
for edge in Cycles:
model.addConstr(quicksum(x[tuple(cycle)] for cycle in Cycles[edge]) == 1)
model.setObjective(quicksum(Cost[cycle] for cycle in Cost), GRB.MINIMIZE)
model.optimize()
for cycle in x:
if x[cycle].X > 0.1:
print(cycle, x[cycle].X, Cost[cycle])