「PuLPで解く線形計画法と連立方程式とナップサック問題」と「PuLPで解く最適割当問題と最適輸送問題」の続編として、今度は最大流問題と最小費用流問題を解いてみます。
インストール
!pip install pulp
Collecting pulp
[?25l Downloading https://files.pythonhosted.org/packages/14/c4/0eec14a0123209c261de6ff154ef3be5cad3fd557c084f468356662e0585/PuLP-2.4-py3-none-any.whl (40.6MB)
[K |████████████████████████████████| 40.6MB 76kB/s
[?25hCollecting amply>=0.1.2
Downloading https://files.pythonhosted.org/packages/f3/c5/dfa09dd2595a2ab2ab4e6fa7bebef9565812722e1980d04b0edce5032066/amply-0.1.4-py3-none-any.whl
Requirement already satisfied: pyparsing in /usr/local/lib/python3.7/dist-packages (from amply>=0.1.2->pulp) (2.4.7)
Requirement already satisfied: docutils>=0.3 in /usr/local/lib/python3.7/dist-packages (from amply>=0.1.2->pulp) (0.17.1)
Installing collected packages: amply, pulp
Successfully installed amply-0.1.4 pulp-2.4
必要なライブラリのインポート
import pulp
import numpy as np
import matplotlib.pyplot as plt
問題設定
座標 (0, 0) から (7, 7) の 8 x 8 のグリッド状の頂点（ノード）を考えます。これらの頂点間に、次のように縦・横・斜めの辺（エッジ）が存在するものとします。
grid_size = 8
plt.figure(figsize=(12,12))
checkpoints = []
for x1 in range(grid_size):
for y1 in range(grid_size):
if x1 == y1:
pass
checkpoints.append((x1, y1))
plt.scatter([x1], [y1])
edges = {}
for p1 in checkpoints:
for p2 in checkpoints:
if p1[0] == p2[0] and abs(p1[1] - p2[1]) == 1:
plt.plot([p1[0], p2[0]], [p1[1], p2[1]])
elif p1[1] == p2[1] and abs(p1[0] - p2[0]) == 1:
plt.plot([p1[0], p2[0]], [p1[1], p2[1]])
elif p1[1] - p2[1] == 1 and p1[0] - p2[0] == 1:
plt.plot([p1[0], p2[0]], [p1[1], p2[1]])
else:
continue
if p1 not in edges.keys():
edges[p1] = {}
if p2 not in edges.keys():
edges[p2] = {}
weight = int(np.random.rand() * 50)
if p2 not in edges[p1]:
edges[p1][p2] = weight
if p1 not in edges[p2]:
edges[p2][p1] = weight
plt.text((p1[0] + p2[0]) / 2, (p1[1] + p2[1]) / 2, edges[p1][p2])
上の図で、辺に書かれている数字は、その辺の容量（その辺に沿って流すことができる「品物」の最大量）とします。
最大流問題
最大流問題（Maximum flow problem）とは、供給点（ソース) から需要点（シンク）まで流すことができる「品物」の最大量を求める問題です。
ここでは、原点 (0, 0) を供給点（ソース) とし、点 (7, 7) を需要点（シンク）とします。
まず PuLP を使って、問題設定するためのインスタンスを作成しましょう。
import pulp
prob = pulp.LpProblem('max_flow_problem', sense = pulp.LpMaximize)
辺ごとの重み（品物の輸送量）の最小値がゼロ、最大値がその辺の容量となるようにします。
val = {}
for k1, v1 in edges.items():
val[k1] = {}
for k2, v2 in v1.items():
val[k1][k2] = pulp.LpVariable("({},{})-({},{})".format(
k1[0], k1[1], k2[0], k2[1]
), lowBound = 0, upBound = v2)
最大化すべき値は、供給点（ソース) から需要点（シンク）まで流すことができる「品物」の最大量ですが、これは「供給点 (0, 0) から出てくる品物の総和」を最大化することと同義です。
prob += pulp.lpSum([v for v in val[(0, 0)].values()])
「供給点（ソース) と需要点（シンク）以外の頂点（ノード）において、入ってくる品物の量の総和と出ていく品物の量の総和は同じである」という制約条件を加えます。
ここで、ノード (i, j) は (i - 1, j), (i, j - 1), (i - 1, j - 1) から品物が入ってきて、(i + 1, j), (i, j + 1), (i + 1, j + 1) に品物が出ていくものとし、逆流はしないものとします。
for x, y in checkpoints:
if x == grid_size - 1 and y == grid_size - 1:
continue
elif x == 0 and y == 0:
continue
sahen = 0
if (x, y - 1) in val[(x, y)]:
sahen += val[(x, y)][(x, y - 1)]
if (x - 1, y) in val[(x, y)]:
sahen += val[(x, y)][(x - 1, y)]
if (x - 1, y - 1) in val[(x, y)]:
sahen += val[(x, y)][(x - 1, y - 1)]
uhen = 0
if (x, y + 1) in val[(x, y)]:
uhen += val[(x, y)][(x, y + 1)]
if (x + 1, y) in val[(x, y)]:
uhen += val[(x, y)][(x + 1, y)]
if (x + 1, y + 1) in val[(x, y)]:
uhen += val[(x, y)][(x + 1, y + 1)]
prob += sahen == uhen
今回の問題設定では、実装を簡単にするため、頂点 k1 から k2 に向かう辺と、 k2 から k1 に向かう辺の流量を同じであるとしておきます。
for k1, v1 in val.items():
for k2, v2 in v1.items():
prob += val[k1][k2] == val[k2][k1]
これで目的関数の入力と制約条件の入力が完了したので、ソルバーを実行します。
result = prob.solve()
pulp.LpStatus[result]
'Optimal'
結果は次のようになりました。
plt.figure(figsize=(12,12))
plt.title("Optimal value = {}".format(prob.objective.value()))
for k1, v1 in val.items():
for k2, v2 in v1.items():
if v2.value() > 0:
plt.plot([k1[0], k2[0]], [k1[1], k2[1]], linewidth=int(v2.value()), alpha=0.5)
#print(k1, k2, v2.value())
plt.text(
(k1[0] + k2[0]) / 2,
(k1[1] + k2[1]) / 2,
"{} / {}".format(int(v2.value()), edges[k1][k2])
)
for x1 in range(grid_size):
for y1 in range(grid_size):
plt.scatter([x1], [y1])
品物の流れがゼロになった辺は省略し、たとえば容量が 33 の辺に品物が 2 だけ流れた時は 2 / 33 のように表示しています。流れの多い辺は太く表示しています。
最小費用流問題
最小費用流問題 (Minimum cost flow problem) は、最大流問題が最大の流量を求めるのに対して、流量を固定したときにそのコストの総量を最小化する問題です。
ここでは、上の最大流問題で求まった最大流量で固定することにしましょう。
total_flow = prob.objective.value()
今回新たに加えるのは、各辺において 1 の量の品物を流すためのコストです。
costs = {}
for k1, v1 in edges.items():
for k2, v2 in v1.items():
cost = int(np.random.rand() * 100)
if k1 not in costs.keys():
costs[k1] = {}
if k2 not in costs.keys():
costs[k2] = {}
costs[k1][k2] = cost
costs[k2][k1] = cost
次のように、最大流問題のときとほぼ同じコードで解けます。違う点は、目的関数が「辺のコスト」と「流量」の積和であるという点です。
import pulp
prob = pulp.LpProblem('minimum_flow_problem', sense = pulp.LpMinimize)
sum_cost = 0
for k1, v1 in val.items():
for k2, v2 in v1.items():
sum_cost = costs[k1][k2] * v2
prob += sum_cost
prob += pulp.lpSum([v for v in val[(0, 0)].values()]) == total_flow
for x, y in checkpoints:
if x == grid_size - 1 and y == grid_size - 1:
continue
elif x == 0 and y == 0:
continue
sahen = 0
if (x, y - 1) in val[(x, y)]:
sahen += val[(x, y)][(x, y - 1)]
if (x - 1, y) in val[(x, y)]:
sahen += val[(x, y)][(x - 1, y)]
if (x - 1, y - 1) in val[(x, y)]:
sahen += val[(x, y)][(x - 1, y - 1)]
uhen = 0
if (x, y + 1) in val[(x, y)]:
uhen += val[(x, y)][(x, y + 1)]
if (x + 1, y) in val[(x, y)]:
uhen += val[(x, y)][(x + 1, y)]
if (x + 1, y + 1) in val[(x, y)]:
uhen += val[(x, y)][(x + 1, y + 1)]
prob += sahen == uhen
for k1, v1 in val.items():
for k2, v2 in v1.items():
prob += val[k1][k2] == val[k2][k1]
result = prob.solve()
pulp.LpStatus[result]
'Optimal'
結果は次のようになりました。
plt.figure(figsize=(12,12))
plt.title("Optimal value = {}".format(prob.objective.value()))
for k1, v1 in val.items():
for k2, v2 in v1.items():
if v2.value() > 0:
plt.plot([k1[0], k2[0]], [k1[1], k2[1]], linewidth=int(v2.value()), alpha=0.5)
#print(k1, k2, v2.value())
plt.text((k1[0] + k2[0]) / 2, (k1[1] + k2[1]) / 2, "{} / {}".format(int(v2.value()), edges[k1][k2]))
for x1 in range(grid_size):
for y1 in range(grid_size):
plt.scatter([x1], [y1])
最大流問題で得られた結果と、最小費用流問題で得られた結果は、流量は同じですが微妙に異なりますね。どちらも、解が複数ある場合はその中の１つしか出力してくれないと思います。
最小費用流問題の総流量を別の数字に変更すると当然違う結果が得られます。試してみてください。