第9章：エネルギー貯蔵 II
章の概要
本章は、第 8 章のモデルを拡張し、風力発電所とグリッドという2つの電源からのエネルギーを組み合わせて、蓄電デバイスの助けを借りながら時間依存の負荷に応える、より複雑なエネルギー貯蔵問題から出発する。この問題を特徴づけているのは、1時間ごとに更新される24時間先の風力予測が与えられている点である。これらの予測は時間とともにかなり変化し得るため、単に予測が不完全であるというだけでなく、予測そのものが変化するという新たな不確実性の源が生じる。
まず、風力予測の不確実性をモデル化するための2つの方法を検討する。1つ目はガウス過程回帰と呼ばれる手法であり、今後24時間に発電されると見込まれる風力エネルギー量のような連続的な過程をモデル化するのに有用である。
2つ目の方法は、「隠れ状態モデル」と呼ばれるものに基づく、強力でありながら驚くほど単純な手法を用いる。これにより、交差時間（crossing times）として知られる予測の性質を再現できる。交差時間とは、予測値が実際値を上回っている、あるいは下回っている時間の長さを指す。これは貯蔵問題をモデル化する上で重要な挙動である。
方策を設計するにあたっては、第 6 章で最初に紹介した手法を応用する。そこでは決定論的先読み方策から出発し、時間とともにうまく機能するようにパラメータを導入した。今回のエネルギー問題では、風力から得られるエネルギーの予測値の最良推定値に、何時間先を予測しているかに応じた係数を掛け合わせて計画を立てる。これにより、調整が必要な24個のパラメータを持つ決定論的最適化モデルが得られる。
ナラティブ
ここからは、図 9.1 に示されるやや複雑なエネルギー貯蔵問題を解いていく。これまでのようにグリッドとの間でエネルギーを売買するだけの貯蔵システムとは対照的に、今回は風力発電所とグリッドからのエネルギーを用いて建物の時間依存の負荷を満たすという問題に直面し、これらの異なるプロセスを平準化するための単一の蓄電デバイスが用いられる。
この問題はもう一つの特徴的な性質も示す。すなわち、すべての外生プロセス（風力、価格、負荷、気温）が、それぞれ異なる予測可能性を持つ動的過程から生じるという点である。
- 負荷 – 負荷（エネルギー需要）は、時刻（人々が出社し始める午前8時までに建物を特定の温度にする必要がある、など）や気温に依存する、比較的予測可能なパターンに従う。
- 気温 – 気温は、時刻や季節に依存する、ある程度予測可能な過程であるが、ある程度の精度で予測できる局地的な気象条件も反映する。
- 風力 – 風力予測サービスを提供するベンダーが存在するが、その予測はあまり正確ではなく、1日のうちでもかなり速く変化する（図 9.2 参照）。
- 価格 – グリッド上の電力価格は需給を反映しており、電力供給者は需要に迅速に対応するよう設計されている。しかし、短期的な不足によって価格が平均の10倍から100倍にまで急騰するスパイクが生じることがある。また、負荷が発電機の出力を下げる速度よりも速く落ち込む期間もあり、余剰電力が非常に低い、時にはマイナスの価格で売却されることもある。
これらの各プロセスは、それぞれ異なる精度で予測することができる。風力発電の予測は最も精度が低く、夜間に風が強くなる傾向はあるものの、高値も低値も昼夜を問わずいつでも発生し得る。負荷は、主に人間の活動によるものだが気温にも起因して、時刻と強い相関を持つ。なお、夏の午後の暑さはエアコン使用によって昼の真ん中にピークを生み出す一方、冬季には（電気暖房でまかなわれる場合の）暖房負荷をかえって減少させることがある点に注意されたい。気温もまた、太陽の出入りにより強い時刻依存性を持つが、気象前線の通過によって変動が生じることもある。
我々の問題は、各時点でグリッドからどれだけ購入するか（あるいはグリッドに売り戻すか）、そしてどれだけ貯蔵するか（一部の系統運用者ではこれらの決定を5分刻みで行うこともある）を決めることである。エネルギー需要を満たす必要があるが、それ以外の点では、グリッドまたは風力発電所からエネルギーを購入するコストを差し引いた、エネルギー販売による収益を最大化したい。
問題の枠組み化
3つの枠組み化の質問に対する答えは以下の通りである。
- メトリクス: 需要を満たすことで得られる収益から、グリッドからのエネルギー購入コストを差し引いた、総期待利益を最大化したい。
- 決定: 6つの決定変数がある。グリッドから蓄電への電力の流れ、グリッドから負荷への電力の流れ、風力発電所から蓄電への電力の流れ、風力発電所から負荷への電力の流れ、蓄電からグリッドへの電力の流れ、そして蓄電から負荷への電力の流れである。
- 不確実性: 以下の動的情報過程がある。負荷（エネルギー需要）、気温（負荷に影響する）、風力発電所からのエネルギー、負荷を満たすことで受け取る価格、そしてグリッドからエネルギーを購入するコストである。
基本モデル
状態変数
まず、時刻 $t$ におけるシステムのスナップショットをモデル化する。これには、時刻 $t$ にバッテリーに貯蔵されているエネルギー量（MWh単位）である $R_t$、時刻 $t$ におけるエネルギー需要（負荷）（MW単位）である $L_t$、時刻 $t$ における気温である $\tau_t$、時刻 $t$ における風力からのエネルギー（MW単位）である $w_t$、時刻 $t$ に建物への負荷を満たすためにMWhあたりで支払われる金額である $p^{load}_t$、そしてグリッドから電力を購入するコスト（これはグリッドへ売り戻した場合に支払われる価格でもある）である $c^{grid}_t$、が含まれる。
対象とする問題は（日次サイクルのため）非常に時間依存性が強いので、問題のダイナミクスをモデル化するためにも、また将来起こり得ることを見越して決定を下すためにも、予測を利用する必要がある。図 9.2 に示した風力の例のように、ローリング形式の予測が一連与えられていると仮定する。負荷（$L$）、気温（$\tau$）、風力（$w$）、市場価格（$p$）、グリッド価格（$G$）の予測を以下のようにモデル化する。$f^L_{tt’}$ は、時刻 $t$ 時点で分かっている情報のもとでの、時刻 $t’ > t$ における負荷 $L_t$（MW単位）の予測；$f^\tau_{tt’}$ は、時刻 $t$ 時点で分かっている情報のもとでの、時刻 $t’ > t$ における気温 $\tau_t$ の予測；$f^w_{tt’}$ は、時刻 $t$ 時点で分かっている情報のもとでの、時刻 $t’ > t$ における風力発電量 $w_t$（MW単位）の予測；$f^p_{tt’}$ は、時刻 $t$ 時点で分かっている情報のもとでの、$t’ > t$ における市場価格 $p^{load}_t$（＄/MWh単位）の予測；そして $f^G_{tt’}$ は、時刻 $t$ 時点で分かっている情報のもとでの、$t’ > t$ におけるグリッド価格 $c^{grid}_t$（＄/MWh単位）の予測である。
すべての予測は、計画期間 $t, t+1, \ldots, t+H$ 全体にわたるベクトルであり、$H$ は指定された計画期間（例：24時間）である。$X \in \Xcal = \lbrace L, T, W, P, G\rbrace $ に対する予測のベクトルを $f^X_t$ とする。
このとき、状態変数は次のようになる。
\[S_t = (\underbrace{R_t}_{R_t}, \underbrace{(L_t, \tau_t, w_t, p^{load}_t, c^{grid}_t)}_{I_t}, \underbrace{(f^L_t, f^T_t, f^w_t, f^P_t, f^G_t )}_{B_t}).\]ここでは、制御可能な資源 $R_t$（我々の物理状態変数）、負荷・気温・風力発電量・価格のスナップショット（情報変数 $I_t$としてまとめることができる）、そして将来に関する一種の信念 $B_t$を表す予測、をグループ化している。
我々は、比較的高次元の状態変数を扱っていることにすぐ気づく。5分刻みで計画を立てる場合、24時間のローリング予測は288個の要素を持つことになる。これは、状態 $S_t$ にいることの価値 $V_t(S_t)$ を推定しようとする人が直面する課題を示唆している。
決定変数
システムの決定変数は次の通りである。時刻 $t$ に風力発電所からバッテリーへ移動する電力量である $x^{wr}_t$；時刻 $t$ に風力発電所から負荷（建物）へ移動する電力量である $x^{w\ell}_t$；時刻 $t$ にグリッドからバッテリーへ移動する電力量である $x^{gr}_t$；時刻 $t$ にバッテリーからグリッドへ移動する電力量である $x^{rg}_t$；時刻 $t$ にグリッドから負荷へ移動する電力量である $x^{g\ell}_t$；時刻 $t$ にバッテリーから負荷へ移動する電力量である $x^{r\ell}_t$；そして、カバーされなかった負荷（「負荷遮断（load shedding）」として知られる）である $x^{loss}_t$。
これらの変数は、以下の制約条件のもとで決定されなければならない。
\[\begin{align} x^{w\ell}_t + x^{g\ell}_t + \frac{1}{\eta} x^{r\ell}_t + x^{loss}_t &= L_t, \label{eq:energysystem1}\\ x^{r\ell}_t + x^{rg}_t &\leq \eta R_t, \label{eq:energysystem2} \end{align}\] \[\begin{align} x^{wr}_t + x^{gr}_t &\leq \frac{1}{\eta} (R^{max} - R_t), \label{eq:energysystem3}\\ x^{rg}_t &\leq \eta R_t, \label{eq:energysystem3a}\\ x^{w\ell}_t + x^{wr}_t &\leq w_t, \label{eq:energysystem4}\\ x^{wr}_t + x^{gr}_t &\leq \frac{1}{\eta} u^{charge}, \label{eq:energysystem5}\\ x^{r\ell}_t + x^{rg}_t &\leq \eta u^{discharge}, \label{eq:energysystem6}\\ x^{wr}_t, x^{w\ell}_t, x^{gr}_t, x^{rg}_t, x^{g\ell}_t, x^{r\ell}_t &\geq 0. \label{eq:energysystem7} \end{align}\]式 $\eqref{eq:energysystem1}$ は、負荷（建物）に供給する電力を負荷量以下に制限する（建物を過負荷にすることはできない）。変数 $x^{loss}$ は、負荷をカバーできなかった量を表す。この制約は、バッテリーから取り出したエネルギーについての変換損失を捉えている。式 $\eqref{eq:energysystem2}$ は、変換損失を調整した上で、バッテリー貯蔵から実際に貯蔵されている量以上の電力を取り出すことはできないことを表す。式 $\eqref{eq:energysystem3}$ は、同様に変換損失を調整した上で、バッテリーに移動できる量を利用可能な容量に制限する。式 $\eqref{eq:energysystem3a}$ は、貯蔵からグリッドへ戻せる量を制限する。式 $\eqref{eq:energysystem4}$ は、風力発電所からの量を、その時点で風力発電所が発電している量に制限する。式 $\eqref{eq:energysystem5}$–$\eqref{eq:energysystem6}$ は、バッテリーへの出入りの流れを充放電レートに制限する。式 $\eqref{eq:energysystem7}$ は、各変数に非負制約を課す。
外生情報
外生情報の1つ目の源は、任意の過程「$X$」について、実測値と予測値の差である。
\[X = (L, \tau, w, p^{load}, c^{grid}).\]$X_t$ をその過程、$\varepsilon^X_{t+1}$ を実測値と予測値の差とする。このとき、次のようになる。
\[\varepsilon^X_{t+1} = X_{t+1} - f^X_{t,t+1}.\]$\varepsilon^X_{t+1}$ は、過去データから抽出したサンプルを用いてモデル化することもできるし、何らかの仮定された分布に従うと仮定してモデル化することもできる。
外生情報の2つ目の源は、時間を1ステップ進めるごとに予測が変化することである。ここでも、各情報過程 $X$ に対する予測のベクトルを $f^X_t$ とする。ここで $f^X_{tt}$ は時刻 $t$ における実測値である。$t$ と $t+1$ の間で、時刻 $t’$ に対する予測の変化を $\fhat^X_{t+1,t’}$ とすると、次のようになる。
\[\fhat^X_{t+1,t'} = f^X_{t+1,t'} - f^X_{tt'},~ t'=t, t+1, \ldots, t+H.\]外生的な変化 $\fhat^X_{t+1,t’}$ は、時間期間 $t’$ をまたいで相関を持つ。もしそうでなければ、計画期間 $t’=t, \ldots, t+H$ にわたってプロットした予測は、図 9.2 の風力予測に見られるような滑らかさを示さなくなるだろう。予測の不確実性のモデル化という問題については、後ほど改めて取り上げる。
このことから、外生情報を次のように書くことができる。
\[W_{t+1,X} = (\varepsilon^X_{t+1}, \fhat^X_{t+1,t'}), t' > t,\]ここで $X$ は、各変数（負荷、気温、風力、市場価格、グリッド価格）に対応する。
遷移関数
資源状態変数の推移は次のように与えられる。
\[\begin{align} R_{t+1} = R_t + \eta (x^{wr}_t + x^{gr}_t) - \frac{1}{\eta} (x^{rg}_t + x^{r\ell}_t).\label{eq:energytransitionII1} \end{align}\]$L_t$、$\tau_t$、$w_t$、$p^{load}_t$、$c^{grid}_t$ の各変数は、予測を用いて推移する。例えば、負荷 $L_t$ の推移は次のように書くことができる。
\[\begin{align} L_{t+1} = f^L_{t,t+1} + \varepsilon^L_{t+1}, \label{eq:energytransitionII2} \end{align}\]$\tau_t$、$w_t$、$p^{load}_t$、$c^{grid}_t$ についても同様の式を作成できる。
予測の推移は次のように書く。
\[\begin{align} f^X_{t+1,t'} = f^X_{tt'} + \fhat^X_{t+1,t'}, ~X\in\Xcal, ~t'=t+1, \ldots, t+1+H, \label{eq:energytransitionII3} \end{align}\]ここで $X=L, \tau, w, p^{load}$ かつ $c^{grid}$ である。式 $\eqref{eq:energytransitionII3}$ は、文献において「予測推移のマルチンゲールモデル」として知られている。「マルチンゲール」という用語は、無作為な偏差 $\fhat^X_{t+1,t’}$ が平均としてゼロであると仮定していることから、$f^X_{tt’}$ が $f^X_{t+1,t’}$ の不偏推定値であるという我々の仮定を指している。
式 $\eqref{eq:energytransitionII1}$、$\eqref{eq:energytransitionII2}$、$\eqref{eq:energytransitionII3}$（すべての予測 $X$ について）が、遷移関数
\[S_{t+1} = S^M(S_t,x_t,W_{t+1}).\]を構成する。
目的関数
時刻 $t$ における我々の利益関数は次で与えられる。
\[C(S_t,x_t) = (x^{w\ell}_t + x^{g\ell}_t + \eta x^{r\ell}_t) p^{load}_t - (x^{g\ell}_t + x^{gr}_t)c^{grid}_t,\]市場価格$p^{load}_t$とグリッド価格$c^{grid}_t$は状態変数$S_t$に含まれる。我々の目的関数は依然として次式で与えられる正準問題である。
\[\max_\pi \E \sum_{t=0}^T C(S_t,X^\pi(S_t))\]これまでと同様に、$S_{t+1} = S^M(S_t,X^\pi(S_t),W_{t+1})$は式$\eqref{eq:energytransitionII1}$、$\eqref{eq:energytransitionII2}$、$\eqref{eq:energytransitionII3}$によって与えられる。
不確実性のモデル化
以下では、時間にわたる不確実性をモデル化する2つのスタイルについて説明する。まず、時に ガウス過程回帰 と呼ばれる手法を用いて、時間を通じた予測誤差の相関をモデル化する方法を説明する。この方法により、時間が進むにつれて予測ベクトルが自然な形で発展することが保証される。
続いて、我々が非常に現実的なサンプルパスを確率過程に与えることを見出した、隠れ状態マルコフモデルについて説明する。このモデルは誤差分布を精緻に再現するだけでなく、交差時間（crossing times）、すなわち実際の過程(例えば風速)が予測などのベンチマークを上回るか下回るかを維持する時間を捉える点でも非常に優れている。予測が実際値を上回るか下回るかの時間を適切に捉えることができれば、それは時間にわたる相関を捉えていることを意味する。
この節では、不確実性のモデル化のためのツールボックスに備えておくべき、いくつかの強力な確率的モデリング手法を示す。確率的モデリングは技術的に高度なものになりうるが、本節はそれを反映している。読者への注意として、本節は不確実性のモデル化に関する他の節よりもはるかに込み入った内容であることを申し添えておく。
予測誤差に対するガウス過程回帰
ガウス過程回帰(略してGPR)は、相関を持つ正規分布確率変数の系列を生成するための単純な手法である。GPRは、連続的な曲面を推定しようとする際に特に有用であり、曲面上のある点が期待より高ければ、その近傍の点も期待より高くなる、といった性質を持つ。
$X_{t’}$を時刻$t’$における任意の外生過程(価格、負荷、気温、風力エネルギー)の実際の結果とし、$f^X_{tt’}$を時刻$t < t’$に行われた$X_{t’}$の予測とする。$X_{t’}$と予測$f^X_{tt’}$との差を表すある誤差$\varepsilon_{t’-t}$を仮定するのが一般的である。そして$\varepsilon^X_{t’-t}$について、例えば次のようなモデルを仮定することになる。
\[\varepsilon^X_{t'-t} \sim N(0, (t'-t) \sigma^2_X).\]我々は、これとはやや異なるアプローチを取り、予測の変化$\fhat^X_{t+1,t’}$の分布が次式で記述されると仮定する。
\[\fhat^X_{t+1,t'} \sim N(0, \sigma^2_X).\]次に、予測の変化$\fhat^X_{t+1,t’}$が時刻$t’$を通じて次の共分散関数で相関しているものと仮定する。
\[\begin{align} Cov(\fhat^X_{t+1,t'},\fhat^X_{t+1,t''}) = \sigma^2_X e^{-\beta\vert t''-t'\vert }. \label{eq:forecastcovariancefunction} \end{align}\]式$\eqref{eq:forecastcovariancefunction}$における共分散関数$Cov(\fhat^X_{t+1,t’},\fhat^X_{t+1,t’’})$は、時間にわたる共分散が、2つの時点の差が大きくなるにつれて減少する相関を示すという性質を捉えている。この単純なモデルは、データから、あるいは判断によって推定しなければならない調整可能なパラメータ$\beta$を導入する。例えば、$\beta$のさまざまな値に対する共分散の値をプロットし、妥当と思われる値を選ぶことも可能であろう。
この共分散関数を用いて、要素$\Sigma^X_{t’t’’} = \sigma^2_X e^{-\beta\vert t’‘-t’\vert }$を持つ共分散行列$\Sigma^X$を作成することができる。コレスキー分解と呼ばれる手法を用いることで、予測の変化について相関を持つサンプルを簡単に作成することができる。まず、共分散行列$\Sigma^X$の、いわば「平方根」に相当するものを作成し、これを下三角行列$L$に格納する。pythonでは、NumPyパッケージを用いて、次のpythonコマンドを使用する。
L = scipy.linalg.cholesky(Sigma_X, lower=True)
$\Sigma^X = L^T L$が成り立つことに注意されたい。これが$L$を$\Sigma^X$の平方根とみなす理由である。
次に、平均0、分散1の正規分布に従う独立な確率変数の系列$Z_{\tau}$($\tau = 1, \ldots, H$について)を生成する。ここで$Z=(Z_{t+1}, Z_{t+2}, \ldots, Z_{t+H})^T$を、これらの独立に分布する標準正規確率変数からなる列ベクトルとする。予測の変化について相関を持つサンプルは、次式を用いて作成することができる。
\[\begin{pmatrix} \fhat^X_{t+1,t+1} \\ \fhat^X_{t+1,t+2} \\ \vdots \\ \fhat^X_{t+1,t+H} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{pmatrix} + L Z.\]この式により、式$\eqref{eq:forecastcovariancefunction}$の指数減衰関数に従って相関する、予測の変化のサンプル集合$\fhat^X_{t+1,t+1}, \ldots, \fhat^X_{t+1,t+H}$が得られる。その結果として、予測誤差の分散が次式に従って時間とともに線形に増加していく、発展する予測の集合が得られることになる。
\[Var(\varepsilon^X_{t'-t}) = (t'-t) \sigma^2_X.\]発展する予測$f^X_{t,t’}, f^X_{t+1,t’}, \ldots$は、図9.2で見た発展する風力予測において確認した挙動を示すことになる。
隠れ状態マルコフモデル
エネルギー分野において確率モデルを構築する際の課題の一つに、交差時間(crossing time)として知られる性質を捉えることが挙げられる。これは、実際の過程(例えば価格や風速)が予測などのあるベンチマークを上回るか下回るかの時間である。図9.3は、風力過程に関する上方交差時間(up-crossing time)を示している。
標準的な時系列モデリングを用いて交差時間を再現しようとしたが、うまくいかなかった。うまくいったのは、隠れ状態変数$S^C_t$を持つマルコフモデルの開発であり、これは過程がベンチマークを上回るか下回るかの動態を捉えるように較正される。この手順は以下のステップからなる。
ステップ1 – 実際の過程をベンチマークと比較し、実際の過程がベンチマークを上回るあるいは下回るように移り変わる時点を見つけ、その過程が上回っていた(A)か下回っていた(B)か、およびその継続時間を捉えたデータセットを出力する。これらの期間を短・中・長(S/M/L)の3つのバケットに集約し、各セグメントにA/BとS/M/Lのラベルを付け、6つの状態を作成する。これらが「隠れ状態」と呼ばれるのは、時刻$t$において実際の過程がベンチマークを上回っているか下回っているかは分かる一方で、その長さが短・中・長のいずれであるかは、過程がベンチマークを交差した後にならないと分からないからである。
ステップ2 – $S^C_t$の過去の系列を用いて、1段階遷移行列$P^C[S^C_{t+1}\vert S^C_t]$、すなわち交差過程が現在状態$S^C_t$にあるという条件のもとで値$S^C_{t+1}$を取る確率を計算する。
ステップ3 – 実際の過程(例えば風速)を、経験的累積分布に基づき、例えば5つのバケットに集約する。$W^g_t$を、集約された風速(1から5までの数)とする。
ステップ4 – 過去のデータから、$W^g_t$と$S^C_t$が与えられたときの風速の条件付き分布、すなわち$F^W[W_{t+1}\vert W^g_t, S^C_t]$、$W^g_t$と$S^C_t$が与えられたときの風速$W_{t+1}$の経験的累積分布を計算する。
1段階遷移行列$P^C[S^C_{t+1}\vert S^C_t]$と条件付き累積分布$F^W[W_{t+1}\vert W^g_t, S^C_t]$を用いることで、我々はまず$S^C_t$が与えられたときの隠れ状態変数$S^C_{t+1}$をシミュレートすることによって(これらは30通りしかないことに注意)、我々の確率過程をシミュレートすることができる。次に、風速$W_t$から、集約された風速$W^g_t$を求め、条件付き累積分布$F^W[W_{t+1}\vert W^g_t, S^C_t]$から実際の風速$W_{t+1}$をサンプリングすることができる。
このロジックは、風力やグリッド価格をモデル化するさまざまなデータセットにわたって、誤差分布(実際値対ベンチマーク)、および上方交差・下方交差の両方の分布を正確に再現することが確認されている。図9.4は、ある特定のデータセットにおけるこれらの分布を示している。
方策の設計
この問題における最大の複雑さは、予測が状態変数の一部であり、これによって予測のローリングな発展を明示的にモデル化できるようになっている点である。難点は、これによって状態変数が高次元になってしまうことである。
予測を扱うための最も一般的なアプローチでは、予測を固定することによって将来を近似する先読みモデルを用いる。最も一般的な2つの戦略は次の通りである。
- 先読みの定式化の中に予測を組み込んだ決定論的先読み。
- 潜在変数を伴う確率的先読み – 予測を用いて確率的先読みモデルを構築し、それを古典的な動的計画法を用いて解くことができる。先読みモデルにおいては、予測の時間にわたる発展をモデル化するのではなく、予測をモデルに固定する。モデル(先読みモデルを含む)の中である変数を無視する場合、それは 潜在変数(latent variable) と呼ばれる。
これら2つの手法はいずれも、先読みモデル内で予測の発展をモデル化しないという点で、予測を潜在変数として用いている。確率的先読みモデルの難点は、それを解くのがより困難であるという点である。我々のエネルギー貯蔵問題を短い時間刻み(5分、あるいはそれ以下かもしれない)で最適化しているのであれば、これは厳密な動的計画法や近似的な動的計画法といった手法にとって問題を引き起こしうる。
このため、予測を伴う時間依存の問題を扱うための最も一般的な戦略は、我々の動的最短経路問題で行ったのとちょうど同じように、決定論的先読みモデルを解くことである。以下でそのようなモデルを説明し、その後、決定論的モデルが不確実性をより良く扱えるようにするパラメータ化されたバージョンを導入する。
ここで、ローリング予測を用いた計画立案がオペレーションズマネジメントにおいて非常に一般的であることに触れておく。奇妙なことに、教科書はほぼ一様にローリング予測の適切なモデル化を無視している。例えば、在庫計画に関するほとんどすべての書籍は、状態変数を在庫と同一視している。もし予測が各時間期間ごとに更新されているのであれば、予測を更新するために必要な情報を状態変数の中で表現しなければならないことを認識しているのは、ごく一部だけである。この情報を無視するならば、我々は実質的に予測を一定に保つ先読みモデルを作っていることになる。
決定論的先読み
我々は、第6章で最初に導入したのと同じ表記スタイルを用いる。そこでは、我々が解こうとしている問題である 基本モデル(base model) と、基本モデルを解くための方策の一形態として解く 先読みモデル(lookahead model) とを区別した。
我々の基本モデルの正準的な定式化が次式であったことを思い出されたい。
\[\max_\pi \E \left\{\sum_{t=0}^T C(S_t,X^\pi(S_t))\vert S_0\right\},\]ここで$S_{t+1} = S^M(S_t,X^\pi(S_t),W_{t+1})$であり、外生情報過程$(S_0, W_1, W_2, \ldots, W_T)$を持つ。変数$W_t$は状態$S_t$および/または決定$x_t$に依存する場合があることに注意されたい。もしそうであれば、変数$W_{t+1}$は$S_t$と$x_t$が分かった後に、その場で生成されなければならない。
我々は、決定論的先読みモデルを定式化することによって方策を作成する。ここではすべての変数にチルダを付して表記し、我々が決定を下している時刻$t$と、先読みモデル内部の時間変数である時刻$t’$の両方によって添字付けする。したがって、$\xtilde_{tt’}$、すなわち時刻$t$に生成される先読みモデルにおける時刻$t’$での決定; $\ctilde_{tt’}$、すなわち$\xtilde_{tt’}$に対するコスト係数; そして$\Rtilde_{tt’}$、すなわち時刻$t$に生成される先読みモデルにおける時刻$t’$でのバッテリー内のエネルギー、を定義することになる。
$x_t = \xtilde_{tt}$、$c_t = \ctilde_{tt}$、以下同様であることに注意されたい。
我々は、次の線形計画問題として決定論的先読み方策$X^{DLA}_t(S_t)$を作成する。
\[\begin{align} X^{DLA}_t(S_t) = \argmax_{x_t, (\xtilde_{tt'},t'=t+1, \ldots, t+H)} \left(C(S_t,x_t) + \sum_{t'=t+1}^{t+H} C(\Stilde_{tt'},\xtilde_{tt'})\right), \label{eq:energydetlookahead0} \end{align}\]ここで
\[C(S_t,x_t) = (x^{w\ell}_t + x^{g\ell}_t + \eta x^{r\ell}_t) p^{load}_t - (x^{g\ell}_t + x^{gr}_t)c^{grid}_t,\] \[C(\Stilde_{tt'},\xtilde_{tt'}) = (\xtilde^{w\ell}_{tt'} + \xtilde^{g\ell}_{tt'} + \eta \xtilde^{r\ell}_{tt'}) \ptilde^{load}_{tt'} - (\xtilde^{g\ell}_{tt'} + \xtilde^{gr}_{tt'})\ctilde^{grid}_{tt'}.\]この問題は、$x_t$に対する制約$\eqref{eq:energysystem1}$–$\eqref{eq:energysystem7}$、および全ての$t’ = t+1, \ldots, t+H$に対する$\xtilde_{tt’}$に対する以下の制約に従って解かれなければならない。
\[\begin{align} \Rtilde_{t,t'+1} -\Big(R_{tt'}+\eta (x^{wr}_{tt'} + x^{gr}_{tt'}) - \frac{1}{\eta} (x^{r\ell}_{tt'} + x^{rg}_{tt'})\Big) &= 0, \label{eq:energydetlookahead1}\\ \xtilde^{w\ell}_{tt'} + \xtilde^{g\ell}_{tt'} + \frac{1}{\eta} \xtilde^{r\ell}_{tt'} + \xtilde^{loss}_{tt'} &\leq f^L_{tt'}, \label{eq:energydetlookahead2} \\ \xtilde^{r\ell}_{tt'} + \xtilde^{rg}_{tt'} &\leq \eta \Rtilde_{tt'}, \label{eq:energydetlookahead3}\\ \xtilde^{wr}_{tt'} + \xtilde^{gr}_{tt'} &\leq \frac{1}{\eta} (R^{max} - \Rtilde_{tt'}), \label{eq:energydetlookahead4} \end{align}\] \[\begin{align} \xtilde^{rg}_t &\leq \eta \Rtilde_{tt'} \label{eq:energydetlookahead4a}\\ \xtilde^{w\ell}_{tt'} + \xtilde^{wr}_{tt'} &\leq f^W_{tt'}, \label{eq:energydetlookahead5}\\ \xtilde^{wr}_{tt'} + \xtilde^{gr}_{tt'} &\leq \frac{1}{\eta} u^{charge}, \label{eq:energydetlookahead6}\\ \xtilde^{r\ell}_{tt'} + \xtilde^{rg}_{tt'} &\leq \eta u^{discharge}, \label{eq:energydetlookahead7}\\ \xtilde_{tt'} &\geq 0. \label{eq:energydetlookahead8} \end{align}\]これらの式は、基本制約$\eqref{eq:energysystem1}$–$\eqref{eq:energysystem7}$の式を反映したものであり、唯一の変更点は、実際の風力$W_t$の代わりに、$\xtilde_{tt’}$、$\Rtilde_{tt’}$のような先読み変数や$f^W_{tt’}$のような予測を用いていることである。
式$\eqref{eq:energydetlookahead0}$–$\eqref{eq:energydetlookahead8}$によって記述されるモデルは比較的単純な線形計画問題であり、これを解くためのパッケージが、MatlabやPythonといった言語で現在では利用可能になっている。
$X^{DLA}_t(S_t)$のような先読み方策は、このような動的で時間変化する問題において広く用いられている。これらは、決定論的最短経路問題で最初に示したように、ローリング方式で解く必要がある。このため、これらは時に「ローリングホライズン手法」あるいは「レシーディングホライズン手法」と呼ばれる。これらの先読み方策に基づく「モデル予測制御」という分野全体が存在する。
このエネルギー貯蔵問題のような応用では、決定論的先読みモデルの使用によって、不確実性を考慮できていないという懸念が生じる。例えば、風力の急激な低下や電力系統における価格の高騰から身を守るために、バッテリーに余分なエネルギーを貯蔵しておきたい場合があるだろう。次節では、不確実性を扱うために決定論的先読みモデルをどのように利用できるかを説明する。
パラメータ化先読み
決定論的先読みが不確実性を扱わないという問題に対処する非常に簡単な方法がある。必要なのは、不確実性のためにモデル（あるいは解）をどのように修正すべきかを考えることである。例えば、予期しない変動に対処するために、将来のために余分な貯蔵の代償を支払いたいと思うかもしれない。もちろん、必要になるかもしれない今この瞬間にエネルギーを貯蔵しておくよう、モデルに強制することはできない。また、あまり正確ではないかもしれない予測を割り引いて評価したい場合もあるだろう。
これらの変更は、制約式$\eqref{eq:energydetlookahead1}$－$\eqref{eq:energydetlookahead8}$を以下のもので置き換えることによって導入できる。
\[\begin{align} \Rtilde_{t,t'+1} -\Big(R_{tt'}+\eta (x^{wr}_{tt'} + x^{gr}_{tt'}) - \frac{1}{\eta} (x^{r\ell}_{tt'} + x^{rg}_{tt'})\Big) &= 0, \label{eq:energydetlookaheadmod1}\\ \xtilde^{w\ell}_{tt'} + \xtilde^{g\ell}_{tt'} + \frac{1}{\eta} \xtilde^{r\ell}_{tt'} + \xtilde^{loss}_{tt'} &= \theta^L_{t'-t} f^L_{tt'}, \label{eq:energydetlookaheadmod2}\\ \xtilde^{r\ell}_{tt'} + \xtilde^{rg}_{tt'} &\leq \eta \Rtilde_{tt'}, \label{eq:energydetlookaheadmod3}\\ \xtilde^{wr}_{tt'} + \xtilde^{gr}_{tt'} &\leq \frac{1}{\eta} (R^{max} - \Rtilde_{tt'}), \label{eq:energydetlookaheadmod4}\\ \xtilde^{rg}_t &\leq \eta \Rtilde_{tt'} \label{eq:energydetlookaheadmod4a}\\ \xtilde^{w\ell}_{tt'} + \xtilde^{w\ell}_{tt'}&\leq \theta^W_{t'-t} f^W_{tt'}, \label{eq:energydetlookaheadmod5}\\ \xtilde^{wr}_{tt'} + \xtilde^{gr}_{tt'} &\leq \frac{1}{\eta}u^{charge}, \label{eq:energydetlookaheadmod6}\\ \xtilde^{r\ell}_{tt'} + \xtilde^{rg}_{tt'} &\leq \eta u^{discharge}, \label{eq:energydetlookaheadmod7}\\ \xtilde_{tt'} &\geq 0. \label{eq:energydetlookaheadmod8} \end{align}\]制約式$\eqref{eq:energydetlookaheadmod2}$および$\eqref{eq:energydetlookaheadmod5}$の右辺を修正するためにパラメータを導入したことに注意されたい。ここでは、負荷と風力の予測を修正するための係数$\theta^L_{t’-t}$および$\theta^W_{t’-t}$を導入しており、この係数は何時間先の予測かによってインデックス付けされている。次に、予備を維持するために貯蔵中のエネルギーをすべて使い切る能力を制限したいという考えのもとで、制約式$\eqref{eq:energydetlookaheadmod3}$を修正した。
パラメータ化された制約式$\eqref{eq:energydetlookaheadmod1}$－$\eqref{eq:energydetlookaheadmod8}$のもとで解かれる先読み方策を$X^{DLA-P}(S_t\vert \theta)$と表すことにする。これらのパラメータ化をどのように導入するかを決定した後（これはあらゆるパラメトリックモデルの背後にある技術である）、$\theta$の最良値を見つけるという問題が残る。これは、7章で扱ったパラメータ探索の問題である。
我々は、ベクトル$\theta=(\theta^L, \theta^W)$の最良値を探索する最適化されたパラメータ化決定論的先読みを用いた場合の相対的な改善を、これらのパラメータを1.0に設定する基本方策と比較して計算した。実験では、$\theta^L_{t’-t} = 1$と設定し、風力予測の係数$\theta^W_{t’-t}$のみを最適化した。
その結果を図9.5に示す。これによると、平均で約30パーセントの性能改善が見られる。重要な点は、この改善が現場で意思決定を行う際の追加の複雑さを伴わないということである。ここでは説明していない唯一のステップは、パラメータベクトル$\theta$を調整しなければならないという点である。残念ながら、$\theta$の最適化のプロセスは容易ではない。
本章で学んだこと
- 本章では、ローリング予測を伴う、はるかに複雑なエネルギー貯蔵問題を導入した。
- 「予測進化のマルチンゲールモデル」を導入した。これは、将来の予測が時間とともに進化し、予測の期待変化がゼロである（ただし実際の変化は正または負である）と仮定するものである。
- 次に、隠れセミマルコフモデルを説明し、予測が実際の値を上回っている、あるいは下回っている時間である「交差時間」を再現するのに役立てた。
- 決定論的先読み方策を導入し、続いてパラメータ化決定論的先読みを導入した。ここでは各予測に対する係数を導入した（これはDLA/CFAのハイブリッドのもう一つの例である）。
- 調整済みのDLA/CFA方策が、純粋な（未調整の）決定論的先読みを約30パーセント上回ることを示した。
演習問題
復習問題
- 「予測進化のマルチンゲールモデル」とは何を意味するか。
- 時刻$t$における負荷$L_t$などの各量の予測全体、$t'=t, \ldots, t+H$に対する$f^L_{tt'}$は、状態変数に含まれる。それはなぜか。［ヒント：予測とそれが予測している変数についての遷移方程式を見よ。］
- 変数$x_t$、$t=0, \ldots, T$と、$t' = t, \ldots, t+H$に対する変数$\xtilde_{tt'}$との違いは何か。
- 「交差時間」とは何か。
- 「コレスキー分解」とは何か、また何に用いられるか。
- 風力エネルギーの隠れ状態マルコフモデルにおいて、隠れている状態は何か。それがなぜ隠れているのかを説明せよ。
- 上記のパラメータ化先読み方策は、どのクラスの方策に属するか。最良の調整パラメータの組を見つけるために、どのような目的関数が用いられるか。
問題解決演習
- 本章の問題の条件のもとで意思決定を行うためのパラメータ化方策を設計してみよ。ルールやパラメータ化された関数の形であれば何を用いてもよい。唯一の制限は、何かに対して最適化を行うことが許されない（すなわち、方策の中で$\argmax_x$を使用できない）という点である。
- 我々のパラメータ化先読みは、予測の前に係数を導入することに限定されていた。エネルギー貯蔵装置がその容量に近づきすぎないようにする（これにより予測を超える風力の急増を貯蔵できるようになる）、あるいはゼロに近づきすぎないようにする（風力の落ち込みに備えて）といった、加算的な調整を導入することもできる。代替のパラメータ化を提案し、あなたの構造がなぜ価値を付加しうるかについて論じよ。