MyEnigma

とある自律移動システムエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

最大流問題(Max Flow Problem)とPython, Juliaによる解法サンプルプログラム

目次

はじめに

これまで、最適化問題の代表的な問題である

ポートフォリオ最適化や、

myenigma.hatenablog.com

ナップサック問題問題などを

myenigma.hatenablog.com

最適化ライブラリで解く方法を紹介してきましたが、

今回は最大流問題(Max Flow Problem)という

こちらも有名な最適化問題の概要と、

それをプログラムで解く方法について紹介したいと思います。

 

最大流問題(Max Flow Problem)

最大流問題(Max Flow Problem)は、

代表的な最適化問題の一つです。

dic.nicovideo.jp

- 最大フロー問題 - Wikipedia

 

各エッジの最大流量が指定された有効グラフの中で

スタートからゴールまで最流量が大きくなるような

各エッジの流量を求める問題です。

例えば、下記のような有効グラフにおいて、

image141.png (347×209)

sからtまで最も流量が大きくなる流量を計算する問題です。

 

この最大流量問題は、

線形計画問題であるため、

myenigma.hatenablog.com

大規模なノードの問題でも高速に最適値を計算することができます。

 

この最大流量最適化は、

ネットワークの最適化に使用されたり、

リソースの最適配分などに利用されます。

 

Juliaによるサンプルプログラム

今回は、以前紹介したJuliaという言語と、

myenigma.hatenablog.com

最適化ライブラリJuMPを使って、

myenigma.hatenablog.com

最大流問題を解いてみることにします。

 

今回は下記のような、

有効グラフにおいて、ノード0からノード5までの

流量が最も大きくなる各エッジの流量を求めたいと思います。

ford_fulkerson11.png (394×183)

 

下記がサンプルプログラムです。

基本的には、

各エッジの流量を変数とし、

正の値であることと、

各エッジの最大流量の制約を追加します。

そして、スタートノード(0)とゴールノード(5)以外の

ノード(1-4)の流入と流出は収支が合う必要があるため、

流入と流出の等号制約を追加します。

あとはそれを解くだけです。

前述の通り、最大流量問題は線形計画問題なので、

Clpソルバーで解くことができます。

COIN-OR Linear Programming Solver

 

#
# max flow problem optimization
#
# author: Atsushi Sakai
#

using JuMP
using Clp
solver = ClpSolver()


function main()
    println(PROGRAM_FILE," start!!")

    model = Model(solver=solver)
    @variable(model, 0<=x01<=16)
    @variable(model, 0<=x02<=13)
    @variable(model, 0<=x12<=10)
    @variable(model, 0<=x21<=4)
    @variable(model, 0<=x13<=12)
    @variable(model, 0<=x32<=9)
    @variable(model, 0<=x24<=14)
    @variable(model, 0<=x43<=7)
    @variable(model, 0<=x35<=20)
    @variable(model, 0<=x45<=4)


    @constraint(model, x01 + x21 == x12 + x13 )
    @constraint(model, x02 + x32 + x12 == x21 + x24)
    @constraint(model, x13 + x43 == x32 + x35)
    @constraint(model, x24 == x43 + x45)

    @objective(model, Max, sum(x35+x45))
    status = solve(model)

    println("x01:",getvalue(x01))
    println("x02:",getvalue(x02))
    println("x12:",getvalue(x12))
    println("x21:",getvalue(x21))
    println("x13:",getvalue(x13))
    println("x32:",getvalue(x32))
    println("x24:",getvalue(x24))
    println("x43:",getvalue(x43))
    println("x35:",getvalue(x35))
    println("x45:",getvalue(x45))
    println("final flow:",getvalue(x45)+getvalue(x35))

end


if contains(@__FILE__, PROGRAM_FILE)
    @time main()
end

上記のコードは下記のリポジトリでも公開しています。

github.com

 

上記のコードを実行すると、

下記の結果が表示されます。

f:id:meison_amsl:20170917143311p:plain

 

上記の結果の通り、

各エッジの流量と、最終的な最大流量が計算できました。

 

scipyによる最大流問題の解法

続いて、同じ問題をscipy.optimize.linprogで解いてみたいと思います。

docs.scipy.org

 

JuMPと異なり、最適化のモデリングはしてくれないので、

手計算で線形計画法の標準形に計算しなおす必要がありました。

myenigma.hatenablog.com

import numpy as np
from scipy.optimize import linprog

c = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0])
A_eq = np.array([
    [1.0, 0.0, -1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 1.0, 1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 1.0, -1.0, 0.0, 1.0, -1.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, 0.0, -1.0],
])
b_eq = np.zeros((4,))
lb = np.zeros((10,))
ub = np.array([16.0, 13.0, 10.0, 4.0, 12.0, 9.0, 14.0, 7.0, 20.0, 4.0])
bounds = [(ilb, iub) for ilb, iub in zip(lb, ub)]

result = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
print(result)

下記が実行結果です。

f:id:meison_amsl:20200405220512p:plain

若干、収束判定のしきい値などにより、誤差はありますが、

ほぼJuMPと同じ結果が得られていることがわかります。

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

 

MyEnigma Supporters

もしこの記事が参考になり、

ブログをサポートしたいと思われた方は、

こちらからよろしくお願いします。

myenigma.hatenablog.com