PuLPで解く最大流問題と最小費用流問題

「[PuLPで解く線形計画法と連立方程式とナップサック問題](https://qiita.com/maskot1977/items/41f5e72309b5274f54a0)」と「[PuLPで解く最適割当問題と最適輸送問題](https://qiita.com/maskot1977/items/995588c4eedf2a7b3a01)」の続編として、今度は最大流問題と最小費用流問題を解いてみます。

# インストール


```python
!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


# 必要なライブラリのインポート


```python
import pulp
import numpy as np
import matplotlib.pyplot as plt
```

# 問題設定

座標 `(0, 0)` から `(7, 7)` の 8 x 8 のグリッド状の頂点（ノード）を考えます。これらの頂点間に、次のように縦・横・斜めの辺（エッジ）が存在するものとします。


```python
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])
```

![PuLPで解く最大流問題と最小費用流問題_5_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/e4279eb1-0523-5e85-f99a-bc4c5001e90e.png)

上の図で、辺に書かれている数字は、その辺の容量（その辺に沿って流すことができる「品物」の最大量）とします。



# 最大流問題

最大流問題（Maximum flow problem）とは、供給点（ソース) から需要点（シンク）まで流すことができる「品物」の最大量を求める問題です。

ここでは、原点 `(0, 0)` を供給点（ソース) とし、点 `(7, 7)` を需要点（シンク）とします。

まず PuLP を使って、問題設定するためのインスタンスを作成しましょう。

```python
import pulp

prob = pulp.LpProblem('max_flow_problem', sense = pulp.LpMaximize)
```

辺ごとの重み（品物の輸送量）の最小値がゼロ、最大値がその辺の容量となるようにします。

```python
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)` から出てくる品物の総和」を最大化することと同義です。

```python
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)` に品物が出ていくものとし、逆流はしないものとします。

```python
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` に向かう辺の流量を同じであるとしておきます。

```python
for k1, v1 in val.items():
    for k2, v2 in v1.items():
        prob += val[k1][k2] == val[k2][k1]
```


これで目的関数の入力と制約条件の入力が完了したので、ソルバーを実行します。

```python
result = prob.solve()
pulp.LpStatus[result]
```




    'Optimal'



結果は次のようになりました。

```python
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])
```

![PuLPで解く最大流問題と最小費用流問題_13_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/912a543a-afac-0222-be40-d643302b3f77.png)

    
品物の流れがゼロになった辺は省略し、たとえば容量が `33` の辺に品物が `2` だけ流れた時は `2 / 33` のように表示しています。流れの多い辺は太く表示しています。
    


# 最小費用流問題

最小費用流問題 (Minimum cost flow problem) は、最大流問題が最大の流量を求めるのに対して、流量を固定したときにそのコストの総量を最小化する問題です。

ここでは、上の最大流問題で求まった最大流量で固定することにしましょう。

```python
total_flow = prob.objective.value()
```

今回新たに加えるのは、各辺において `1` の量の品物を流すためのコストです。

```python
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
```

次のように、最大流問題のときとほぼ同じコードで解けます。違う点は、目的関数が「辺のコスト」と「流量」の積和であるという点です。


```python
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'



結果は次のようになりました。

```python
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])
```

![PuLPで解く最大流問題と最小費用流問題_18_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/9ca7ea32-d0b4-390d-f2f4-af0373d0b06d.png)

    
最大流問題で得られた結果と、最小費用流問題で得られた結果は、流量は同じですが微妙に異なりますね。どちらも、解が複数ある場合はその中の１つしか出力してくれないと思います。

最小費用流問題の総流量を別の数字に変更すると当然違う結果が得られます。試してみてください。