目次
- 目次
- はじめに
- シンプルな電力供給最適化問題
- cvxpyによる電力供給最適化問題解法サンプルコード
- scipyによる電力供給最適化問題解法サンプルコード
- 参考資料
- MyEnigma Supporters
はじめに
これまで、
様々な代表的な最適化問題を解いてきましたが、
今回は、
最小コストネットワークフロー問題という問題を解くことで、
シンプルな電力供給最適化問題を解いてみたいと思います。
シンプルな電力供給最適化問題
今回、実際にPythonコードで解いてみようと思う
シンプルな電力供給最適化問題は、
下記のようなグラフで表すことができるとします。
まず初めに、上記の各ノードは
それぞれ電気を必要とする街であるとします。
(都市でも、電気が必要な機械でもなんでもOKです。)
それぞれのノードは方向をもつエッジでつながっており、
それらのエッジは電力を送電するコストを表します。
この場合は各都市の距離と思っても良いと思います。
(遠くに電気を送る方がコストがかかる)
加えて、図には書いていませんが、
下記エッジには流すことができる最大電力を設定します。
今回は、3-5のノードのみ1しか電力を流せないとし、
その他のノードは無限(総電力量と比べて大きな値)としました。
ちなみに今回の例では、
すべてのノードの組み合わせが単方向リストですが、
両方リストにすることも簡単にできます。
最後に、各ノードの上と下に書かれている数値は、
現時点での電力保有量です。
マイナスの値は、電力が不足していることを表し、
プラスは電力を共有することができることを意味します。
上記の例では、ノード2が発電所を持っており、
そこから、各ノードに必要な電力を送電することを考えます。
今回は、ノード1,2,4,5が必要な電力量がわかっており、
それらを供給する分をすべてノード2が発電するという形にしました。
つまり、すべてのノードの電力量を足すと0になるようにできるとします。
上記の図もすべてのノードの電力量を足すと0になるようになっています。
では、上記の状態で、すべてのノードの電力量が0になるように送電し、
かつ、最も総コストが小さくなるように送電するには、
各エッジにどのように電気を送電すれば良いでしょうか?
上記の問題は、一般的にオペレーションズリサーチの分野では、
最小コストフロー問題と呼ばれます。
線形計画法の一種です。
上記の問題を線形計画問題として、
定式化すると下記のようになります。
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()
今回は最適化モデリングツールとして、
Python製のcvxpyを
ソルバーとしてはデフォルトのECOSを使っています。
下記が出力結果です。
上記の結果を図示すると、下記のようになります。
赤文字が最終的な電力送電量です。
すべてのノードに電力が供給されており、
かつ無駄な送電は発生していないように見えます。
このような最適化の計算を逐次リアルタイムに実施することにより、
電力の供給がされているのでしょうか?(実際の所は知らない)
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の結果と
同じになっていることがわかります。
参考資料
MyEnigma Supporters
もしこの記事が参考になり、
ブログをサポートしたいと思われた方は、
こちらからよろしくお願いします。