posted with カエレバ

矢部博数理工学社 2006-04

Amazonで検索

楽天市場で検索

Yahooショッピングで検索

はじめに

これまで、

様々な代表的な最適化問題を解いてきましたが、

myenigma.hatenablog.com

今回は、

最小コストネットワークフロー問題という問題を解くことで、

シンプルな電力供給最適化問題を解いてみたいと思います。

シンプルな電力供給最適化問題

今回、実際にPythonコードで解いてみようと思う

シンプルな電力供給最適化問題は、

下記のようなグラフで表すことができるとします。

f:id:meison_amsl:20171022054108p:plain

まず初めに、上記の各ノードは

それぞれ電気を必要とする街であるとします。

(都市でも、電気が必要な機械でもなんでもOKです。)

それぞれのノードは方向をもつエッジでつながっており、

それらのエッジは電力を送電するコストを表します。

この場合は各都市の距離と思っても良いと思います。

(遠くに電気を送る方がコストがかかる)

加えて、図には書いていませんが、

下記エッジには流すことができる最大電力を設定します。

今回は、3-5のノードのみ1しか電力を流せないとし、

その他のノードは無限(総電力量と比べて大きな値)としました。

ちなみに今回の例では、

すべてのノードの組み合わせが単方向リストですが、

両方リストにすることも簡単にできます。

最後に、各ノードの上と下に書かれている数値は、

現時点での電力保有量です。

マイナスの値は、電力が不足していることを表し、

プラスは電力を共有することができることを意味します。

上記の例では、ノード2が発電所を持っており、

そこから、各ノードに必要な電力を送電することを考えます。

今回は、ノード1,2,4,5が必要な電力量がわかっており、

それらを供給する分をすべてノード2が発電するという形にしました。

つまり、すべてのノードの電力量を足すと0になるようにできるとします。

上記の図もすべてのノードの電力量を足すと0になるようになっています。

では、上記の状態で、すべてのノードの電力量が0になるように送電し、

かつ、最も総コストが小さくなるように送電するには、

各エッジにどのように電気を送電すれば良いでしょうか？

上記の問題は、一般的にオペレーションズリサーチの分野では、

最小コストフロー問題と呼ばれます。

Minimum-cost flow problem - Wikipedia

線形計画法の一種です。

myenigma.hatenablog.com

上記の問題を線形計画問題として、

定式化すると下記のようになります。

f:id:meison_amsl:20171022060601p:plain

xは各エッジの送電量で最適化変数です。

cは各エッジのコスト、bは現在の各ノードの電力量、

uは各エッジの最大送電量です。

コスト関数としては、各エッジのコストと送電料の和を最小化し、

制約条件としては、

各ノードに流入する量と流出する量と、

現在の電力量が辻褄があうようにする条件と、

各エッジの送電料の最大量と最小値(0)を設定しています。

Pythonによるシンプルな電力供給最適化問題解法サンプルコード

前述のシンプルな電力供給最適化問題を解く

Pythonコードが下記です。

"""

Minimal Cost Network Flow

author: Atsushi Sakai

"""

import cvxpy
import numpy as np


def main():
    print(__file__ + " start!!")

    ni = [1, 1, 2, 3, 3, 4, 4, 5]
    nj = [2, 3, 3, 4, 5, 1, 5, 2]
    cij = np.array([4, 2, 1, 1, 2, 5, 2, 4])
    inf = 100000.0
    uij = [inf, inf, inf, inf, 1, inf, inf, inf]
    bi = [-5, 11, -1, -2, -3]
    ni = [1, 1, 2, 3, 3, 4, 4, 5]

    print("ni", ni)
    print("nj", nj)
    print("cij", cij)
    print("uij", uij)
    print("bi", bi)

    x = cvxpy.Variable(len(ni))
    objective = cvxpy.Minimize(cij.T * x)
    constraints = [0 <= x, x <= uij]
    for iin in range(1, len(bi) + 1):
        inp = sum([x[i - 1] for i in range(1, len(nj)) if ni[i - 1] == iin])
        out = sum([x[i - 1] for i in range(1, len(ni)) if nj[i - 1] == iin])
        print(iin, bi[iin - 1])
        constraints += [inp - out == bi[iin - 1]]

    prob = cvxpy.Problem(objective, constraints)

    result = prob.solve(solver=cvxpy.ECOS)
    print("Opt result:", result)
    print("optimal parameter:\n", [float(np.round(i)) for i in x.value])
    print("status:" + prob.status)


if __name__ == '__main__':
    main()