MyEnigma

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

Pythonを使ってシンプルな電力供給最適化問題を解いてみる

目次

はじめに

これまで、

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

myenigma.hatenablog.com

myenigma.hatenablog.com

今回は、

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

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

 

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

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

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

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

f:id:meison_amsl:20171022054108p:plain

 

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

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

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

 

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

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

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

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

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

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

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

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

ちなみに今回の例では、

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

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

 

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

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

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

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

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

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

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

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

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

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

 

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

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

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

 

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

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

線形計画法の一種です。

myenigma.hatenablog.com

 

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

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

f:id:meison_amsl:20171022060601p:plain

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

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

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

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

制約条件としては、

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

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

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

 

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

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

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]

    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()

github.com

 

今回は最適化モデリングツールとして、

Python製のcvxpyを

myenigma.hatenablog.com

ソルバーとしてはデフォルトのECOSを使っています。

 

下記が出力結果です。

f:id:meison_amsl:20171022071028p:plain

上記の結果を図示すると、下記のようになります。

赤文字が最終的な電力送電量です。

f:id:meison_amsl:20171022070635p:plain

 

すべてのノードに電力が供給されており、

かつ無駄な送電は発生していないように見えます。

 

このような最適化の計算を逐次リアルタイムに実施することにより、

電力の供給がされているのでしょうか?(実際の所は知らない)

 

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

上記の問題を

手計算で線形計画問題の標準形に変換し、

scipyのoptimize.linprogで解いたのが、

下記のサンプルコードです。

import numpy as np
from scipy.optimize import linprog

c = np.array([4.0, 2.0, 1.0, 1.0, 2.0, 5.0, 2.0, 4.0])
A_eq = np.array([
    [-1.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0],
    [1.0, 0.0, -1.0, 0.0, 0.0, 0.0,  0.0, 1.0],
    [0.0, 1.0, 1.0, -1.0, -1.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 1.0, 0.0, -1.0, -1.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, -1.0],
])
b_eq = np.array([5.0, -11.0, 1.0, 2.0, 3.0])
lb = np.zeros((8,))
ub = np.array([None, None, None, None, 1.0, None, None, None])
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)
print(np.round(result.x, 2))

出力は下記の通りでcvxpyの結果と

同じになっていることがわかります。

f:id:meison_amsl:20200406202841p:plain

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

 

MyEnigma Supporters

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

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

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

myenigma.hatenablog.com