『エネルギー・リスクマネジメントの数理モデル』の紹介。

『エネルギー・リスクマネジメントの数理モデル』が面白かったので紹介。

エネルギー・リスクマネジメントの数理モデル (確率工学シリーズ)

これは数理最適化 Advent Calendar 2022の7日目の記事です。

概要

タイトルにエネルギーと入ってるので、エネルギーに特化した話しか書かれてないのかと思いきや、実は不確実性に関する数理最適化の話がかなり充実している。 なので、エネルギー分野に限らず、不確実性を扱わないといけない場合にはかなり参考になると思う。

章立ては以下：

まず「1章 確率計画法の基礎」「2章 2段階確率計画問題の解法」で確率計画法の基本をおさえたうえで、「3章 リスクマネジメント」ではCVaRの概念を導入してリスクマネジメントに関する確率計画法を考えられるようにする。 ここまでは一般的な確率計画法の話（数値例として模擬店の話が書かれてるけど、全然エネルギーの話ではない）。

さらに「4章 ロバスト最適化」では確率計画法とは違ったアプローチであるロバスト最適化について説明されていて、「5章 リアルオプション」では金融の世界でよく出てくるオプションの話とそれに対してロバスト最適化を応用した話が書かれている。 これもエネルギー関係なく一般的に使える話。

そして「6章 小売電気事業者の電力調達」「7章 電源投資の経済性評価」「8章 エネルギーサプライチェーンマネジメント」と、申し訳程度に電力関係の応用事例が紹介されている。 ぶっちゃけ7章、8章は自分の関心のあるところではなかったので読んでないんだけどね（＾＾；

この本に関して、個人的にはCVaRに関して理解がちゃんとできたのが大きな収穫だった。 椎名先生の『確率計画法』でもCVaRの説明は書かれてるんだけど、全然理解できなかったので。

あと、ロバスト最適化も名前と概要は聞くけどちゃんと数式では理解してなかったので、それを理解できたのは大きかった。 とくにロバスト最適化というと二次錐計画問題の応用例として語られることが多い印象があったんだけど、不確実性集合として凸多面体を考えることもでき、この場合は線形計画問題の双対理論をうまく使うことで線形計画問題に帰着できるというのがとても面白かった。

競馬の話をロバスト化してみる

この不確実性集合が凸多面体になる場合のロバスト最適化が個人的にツボだったので、具体的な例でちょっと紹介してみたい。

ということで、とりあげるのは先日書いた競馬の話。

ここではオッズが固定値で与えられてたけど、オッズは投票によって変動するので、幅があると考えてみる。

すなわち、1着の馬を当てるとして、1〜3番のどれかが1着になると予想し、オッズも次のようになると予想したとする：

• 1番のオッズは現在17倍で、最終的に16〜18倍だろう
• 2番のオッズは現在5倍で、最終的に2〜8倍だろう
• 3番のオッズは現在8倍で、最終的に6〜10倍だろう
• 1〜3番のオッズの合計は28〜32倍だろう

そして、賭けは1口100円で、当てたときには最低500円はプラスになってほしいとして、購入金額の合計を最小化したい。

ここでもし現在のオッズにしたがって購入すると、最適解は(1, 2, 2)となる。 実際、現在のオッズが維持されるのであれば、1番が勝てば1,200円、2番が勝てば500円、3番が勝てば1,100円儲かる計算になっている。

けど、ここでもし2番のオッズが最終的には2倍まで下がってしまったとすると、400円しか戻ってこないので100円の損失になってしまう。 こういうのを避けたい。

ちなみに、この問題に限れば、勝つ馬はどれか1頭だけなので勝った馬以外のオッズはリターンに影響がなく、したがって実際にはオッズが一番低くなった場合を想定して問題を解けば十分ではあったりする。 つまり、オッズが(16, 2, 6)と想定して最適化問題を解けばいい（オッズの合計が24と制約を満たしてないようだけど、ここでは関係ない）。

ただ、せっかくなのでこの問題をロバスト最適化で解いてみたい。

ロバスト最適化による定式化

ここでは馬の集合を、オッズをとする。 すると、オッズのとりうる値の集合は次のようになる：

これは行列とベクトルを

とすると、

と書け、これは凸多面体になっている。

さて、購入口数をとしたとき、この最適化問題は次のように書ける：

ここで厄介なのがという部分で、は凸多面体なので要素は無数に存在し、素直に考えると制約が無限に必要となってしまう。 けど、これをうまく解決するのが面白いところ。

とりあえず1番の馬についてこの制約を書き直してみると、これはリターンが一番悪かった場合について満たされるべき制約という意味なので、次のように書ける：

これをちょっと変形すると、次のようになる：

ところで、というのは

という線形計画問題の最適値なので、これを主問題と考えたときの双対問題

の最適値で上から抑えられる（弱双対定理）。 （なお、は双対変数）

そこで、この制約は次のように書き換えることが可能：

こうするとという変数がうまいこと消え去ってるのが分かるかと思う。 双対理論をうまく使い、主問題の変数をなくして双対問題の変数に置き換え、という変数同士の積を問題から消して一次の制約に落とし込んでいる。 いやー、お見事って感じ。

これを2, 3番の馬についても適用すると、元の問題は次のように書き換えられる：

ここまで来ればこれはもうただの線形計画問題なので、あとはソルバーで解けばいい。

PuLPを使った求解

コードにして解いてみたのが以下：

import pandas as pd
import pulp

# データ ----------
unit_cost = 100  # 1口あたりの金額
min_profit = 500  # 勝った場合に最低得たい金額

horses = ["1", "2", "3"]

odds_table = pd.DataFrame({
    "min": [16., 2., 6.],   # total: 24
    "max": [18., 8., 10.],  # total: 36
}, index=horses)
odds_total_min = 28
odds_total_max = 32

D = pd.DataFrame({
    "1": [-1, 1, 0, 0, 0, 0, -1, 1],
    "2": [0, 0, -1, 1, 0, 0, -1, 1],
    "3": [0, 0, 0, 0, -1, 1, -1, 1],
}, dtype=float)
b = pd.Series([-16, 18, -2, 8, -6, 10, -28, 32], dtype=float)
aux_index = list(D.index)

# 問題 ----------
problem = pulp.LpProblem(sense=pulp.LpMinimize)

# 変数 ----------
# 購入口数
buy = pulp.LpVariable.dicts("buy", horses, lowBound=0, cat=pulp.LpInteger)
# 補助変数
aux = {}
for horse in horses:
    aux[horse] = pulp.LpVariable.dicts(f"aux_{horse}", aux_index, lowBound=0)

# 目的変数 ----------
total_buy = sum(unit_cost * buy[horse] for horse in horses)
problem += total_buy

# 制約 ----------
# どの馬が勝ったとしても条件を満たす
for horse in horses:
    problem += (
        unit_cost * sum(b[i] * aux[horse][i] for i in aux_index)
        <=
        - total_buy - min_profit
    )
    for other in horses:
        target = - buy[horse] if horse == other else 0
        problem += (
            sum(D[other][i] * aux[horse][i] for i in aux_index)
            ==
            target
        )

# 求解 ----------
solver = pulp.PULP_CBC_CMD(msg=False)
status = problem.solve(solver)

# 出力 ----------
print(pulp.LpStatus[status])
print(f"total: {total_buy.value()}")
for horse in horses:
    print(f"{horse}: {buy[horse].value()}")

これを実行すると次のような出力が得られる：

Optimal
total: 1700.0
1: 2.0
2: 11.0
3: 4.0

つまり、1,700円使って1番を2口、2番を11口、3番を4口買うとよくて、実際、1番が勝ってオッズが16倍になってた場合は1,500円、2番が勝ってオッズが2倍になってた場合は500円、3番が勝ってオッズが6倍になってた場合は700円の儲けがあり、リターンがマイナスになってしまうことが防げている。

---

ちなみにロバスト最適化を業務で試してみたこともあるんだけど、この不確実性集合が小さくできないとなかなかいい解にならない感じはあった。 実際に使っていくときにはそのあたりが肝になってきそう。

あとせっかくなのでちょっと宣伝。

今回の記事では線形計画問題や確率計画法についてはとくに説明しなかったけど、それらについては以下の同人誌に書いたのでよかったら読んでほしい。

今日はここまで！

エネルギー・リスクマネジメントの数理モデル (確率工学シリーズ)