(57)【特許請求の範囲】
【請求項１】
記憶装置、演算装置及び出力装置を備えたコンピュータによって階層的ベイズ推定によりパラメータに分布を持たせたガウシアンプロセスを用いて電力負荷の予測分布の平均と分散を予測する電力負荷予測方法であって、
前記記憶装置に記憶された事前分布に基づいて、ＭＡＰ推定またはハイブリッドモンテカルロ法により前記演算装置が構築した前記ガウシアンプロセスにより、前記演算装置が前記予測分布の平均μＮ＋１と分散σ２Ｎ＋１を下記の式のように算出する第１のステップと、
0004694984.tif 000035
但し、
ｔＮは目的変数
ＣＮは時刻Ｎの事前分布Ｐ（ｔＮ）の分散共分散行列
η，ｒは、ＣＮ＋１とＣＮとの間の補正項のベクトルとスカラ
ＣＮ＋１は時刻Ｎ＋１の事後分布Ｐ（ｔＮ＋１）の分散共分散行列
算出した前記予測分布の平均と分散とを前記出力装置を介して出力する第２のステップと
を備えたことを特徴とする電力負荷予測方法。
【請求項２】
前記演算装置は、前記第１のステップにおいて、ＭＡＰ推定によりＭＡＰ推定値を求めることで、ガウシアンプロセスを構築する
請求項１記載の電力負荷予測方法。
【請求項３】
前記演算装置は、前記第１のステップにおいて、ガウシアンプロセスによる予測分布の平均と分散の計算に必要な変数を、ハイブリッドモンテカルロ法により算出する
請求項１記載の電力負荷予測方法。
【請求項４】
記憶装置、演算装置及び出力装置を備えたコンピュータに階層的ベイズ推定によりパラメータに分布を持たせたガウシアンプロセスを用いて電力負荷の予測分布の平均と分散を予測させる電力負荷予測処理プログラムであって、
前記コンピュータに、
前記記憶装置に記憶された事前分布に基づいて、ＭＡＰ推定またはハイブリッドモンテカルロ法により前記演算装置が構築した前記ガウシアンプロセスにより、前記演算装置が前記予測分布の平均μＮ＋１と分散σ２Ｎ＋１を下記の式のように算出する第１のステップと、
0004694984.tif 000036
但し、
ｔＮは目的変数
ＣＮは時刻Ｎの事前分布Ｐ（ｔＮ）の分散共分散行列
η，ｒは、ＣＮ＋１とＣＮとの間の補正項のベクトルとスカラ
ＣＮ＋１は時刻Ｎ＋１の事後分布Ｐ（ｔＮ＋１）の分散共分散行列
算出した前記予測分布の平均と分散とを前記出力装置を介して出力する第２のステップと
を実行させるための電力負荷予測処理プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、階層的ベイズ推定を用いたガウシアンプロセスにより電力負荷予測の上下限値と平均値とを求める電力負荷予測方法、及び電力負荷予測処理プログラムに関する。
【背景技術】
【０００２】
電力システムにおいて、需要と供給のバランスをとるため、発電計画は、電力システムの安定度と経済性が重要である。従来から、未来の事象を予測するための手法として、過去の実データである事前プロセスをモデル化し、そのモデルの内部変数を推定するといった手法が知られている（特許文献１）。
【０００３】
【特許文献１】特開２００５−７８５１９号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
近年、電力自由化や気象条件の変動により、電力システムにおける不確定性の要因が増大している。電力自由化の新しい環境では、余分な発電量を回避し、利益最大化・リスク最小化を目指す傾向にある。その結果、発電計画の元となる電力負荷予測においても、不確定性を考慮した新しい電力負荷予測手法が切望されている。しかしながら、従来の手法では、電力負荷予測の上下限は系統的に求めることができないため、電力負荷を効率よく推定することができないという課題があった。
【０００５】
本発明は、このような点に鑑みなされたもので、電力負荷予測の上下限を求めることができ、電力負荷を効率よく推定することが可能な電力負荷予測方法を提供することを目的とする。
【課題を解決するための手段】
【０００６】
本発明に係る電力負荷予測方法は、記憶装置、演算装置及び出力装置を備えたコンピュータを用い電力負荷を予測する電力負荷予測方法であって、前記記憶装置に記憶された事前分布に基づいて、ＭＡＰ推定またはハイブリッドモンテカルロ法により前記演算装置が構築したガウシアンプロセスにより、前記演算装置が予測分布の平均と分散を算出するステップと、算出した予測分布の平均と分散とを前記出力装置を介して出力するステップとを備えたことを特徴とする。
【０００７】
前記演算装置は、ＭＡＰ推定によりＭＡＰ推定値を求めることで、ガウシアンプロセスを構築する。
【０００８】
前記演算装置は、ハイブリッドモンテカルロ法により、サンプリングによって予測分布の平均と分散を算出する。
【０００９】
本発明に係る電力負荷予測処理プログラムは、記憶装置、演算装置及び出力装置を備えたコンピュータに用い電力負荷を予測させる電力負荷予測処理プログラムであって、前記コンピュータに、前記記憶装置に記憶された事前分布に基づいて、ＭＡＰ推定またはハイブリッドモンテカルロ法により前記演算装置が構築したガウシアンプロセスにより、前記演算装置が予測分布の平均と分散を算出するステップと、算出した予測分布の平均と分散とを前記出力装置を介して出力するステップとを実行させるためのものである。
【発明の効果】
【００１０】
本発明によれば、電力負荷予測の上下限を求めることができ、電力負荷を精度よく推定することが可能な電力負荷予測方法を提供することができる。
【発明を実施するための最良の形態】
【００１１】
以下、図面を参照しながら本発明の実施形態を詳細に説明する。
本実施形態の電力負荷予測方法では、階層的ベイズ推定を用いたガウシアンプロセスを用いて電力負荷予測の上下限値と平均値とを求める。この際、本実施形態では、後述するＭＡＰ推定あるいはハイブリッドモンテカルロ法により、ガウシアンプロセスの平均値と分散を求める。
【００１２】
図１は、本実施形態に係る電力負荷予測方法を実現するためのコンピュータシステムを示す図である。このシステムは、入力装置（又は受信装置）１、演算装置２、記憶装置３及び出力装置（又は送信装置）４を備えて構成されている。入力装置１は、事前分布などの電力負荷予測に用いられる各種の情報を入力する際に用いられる。演算装置２は、ＭＡＰ推定またはハイブリッドモンテカルロ法によりガウシアンプロセスを構築し、予測分布の平均と分散を算出する処理などを行う。出力装置４は、予測分布の平均と分散を表示あるいは送信する際に用いられる。
【００１３】
ここで、各用語について簡単に説明する。
「ベイズモデル」とは、パラメータとデータを確率変数とみなして予測値を点で推定することではなく、分布として推定するためのモデルである。また、「ガウシアンプロセス（以下、「ＧＰ」という。）」とは、ベイズモデルにおけるパラメータの事前分布として正規分布を仮定するモデルである。階層的ベイズ推定モデルとはパラメータについても事前分布を仮定して、即ち、パラメータの分布のパラメータ（ハイパーパラメ−タ）を考慮したベイズモデルである。なお、階層的ベイズ推定は、特定分布の影響を緩和する働きをする。「ＭＡＰ推定」は、勾配情報を用いて事後確率が最大化するようにハイパーパラメータを決定する推定手法である。そして、「ハイブリッドモンテカルロ法」は、マルコフ連鎖モンテカルロ法の一種であり、サンプリングにより事後分布の平均と分散を決定する手法である。
【００１４】
次に、ＧＰについて詳しく説明する。ここでは、ベイズモデルにおける線形回帰と，近年研究が活発なカーネルマシンとの関連性からのＧＰの定義について説明する。
【００１５】
まず、ニューラルネットワークと同様に、Ｈ個の基底関数ΦＨ＝｛Φｈ（・）｝,（ｈ=１,…,Ｈ）を用い、その重み付き総和で回帰を行うことを考える。但し、Ｈは非常に大きい値であるとする。さらに、Ｎ個の説明変数ｘＮ＝｛ｘｎ｝,（ｎ＝１,…,Ｎ）が与えられているものとする。これらの基底関数と入力変数より、Φｈ（ｘｎ）を要素に持つ行列Ｒを定義する。この行列Ｒの要素Ｒｎｈは、次式のようになる。
【００１６】
【数１】
0004694984.tif 000002
【００１７】
よって、基底関数の重み付き総和をｙ（ｘ）＝｛ｙｎ｝,（ｎ＝１,…, Ｎ） とすると、次式（式２）のように書くことができる。なお、式２において、ｕＨ＝｛ｕｈ｝, （ｈ=１,…,Ｈ）は、基底関数に対応する重みである。
【００１８】
【数２】
0004694984.tif 000003
【００１９】
さらに、｛ｕｈ｝に事前分布を設定する。この事前分布は、次式（式３）に示すような、それぞれの重みに独立に平均「０」であり、かつ分散「σｕ２」の正規分布であると仮定する。なお、式３において、「Ｉ」は単位行列を意味する。
【００２０】
【数３】
0004694984.tif 000004
【００２１】
すると、式２より「ｙ」も同様に平均ゼロの正規分布に従うことが分かる。「ｙ」の分散共分散行列を「Ｑ」とすると、次式（式４）のようになる。
【００２２】
【数４】
0004694984.tif 000005
【００２３】
「ｙ」の分散共分散行列「Ｑ」が式４のようになるので、「ｙ」の事前分布Ｐ（ｙ）は、次式（式５）のように表すことができる。
【００２４】
【数５】
0004694984.tif 000006
【００２５】
次に、目的変数ｔＮ＝｛ｔｎ｝,(ｎ＝１,…,Ｎ)について考える。「ｔｎ」が対応する「ｙｎ」に平均０、分散「σν２」の正規雑音が加わったものとすると、やはり「ｔ」の事前分布も、次式（式６）に示すように正規分布となる。
【００２６】
【数６】
0004694984.tif 000007
【００２７】
ここで、式６の分散共分散行列を「Ｃ」と表すことにすると、この分散共分散行列「Ｃ」は、次式（式７）のようになる。
【００２８】
【数７】
0004694984.tif 000008
【００２９】
もし、基底関数の数Ｈがデータ数Ｎよりも小さければ、ＲＲＴは必ず正則ではなくなる。しかし、ノイズ項「σν２」はフルランクであることから、「Ｃ」は必ず正則となり、数値的な安定性が得られる。このような通常の設定に何らかの条件を加えて不良設定を解消しようという試みのことは「正則化」と呼ばれ、ラジアル基底関数ネットワークの導出の背景となった概念でもある。
【００３０】
さて、分散共分散行列である「Q」について考える。「Q」の（ｉ，ｊ）番目の要素「Qｉｊ」は、次式（式８）のようになる。
【００３１】
【数８】
0004694984.tif 000009
【００３２】
この式８から、実際に「Q」を計算するために必要なことはΦｈ（ｘ） 同士の内積計算であることが分かる。ここで、この内積計算が、次式（式９）のような、共分散関数に相当するカーネル関数で表すことができると仮定する。
【００３３】
【数９】
0004694984.tif 000010
【００３４】
もちろん、式９は任意の関数について成り立つわけではない。式９が成り立つための必要十分条件は、カーネル関数が対称半正定値を満たすことである。このような条件のことを「Ｍｅｒｃｅｒ条件」といい、「Ｍｅｒｃｅｒ条件」を満たすカーネルを特に「Ｍｅｒｃｅｒカーネル」という。また、このようにΦｈ（ｘ）の内積をカーネルで置き換え、Φｈ（ｘ）のような計算を回避することを「カーネルトリック」という。
【００３５】
対称半正定値性は通常の分散共分散行列についても当てはまるため、ＧＰは通常の線形回帰をカーネルにより拡張したモデルと見なすことができる。「Ｍｅｒｃｅｒカーネル」としては、例えば次式（式１０）に示す「ガウシアンカーネル」がある。但し、式１０において、「ω」はガウシアンカーネルの幅である。
【００３６】
【数１０】
0004694984.tif 000011
【００３７】
以上のことから、式７の要素は次式（式１１）のように表すことができる。
【００３８】
【数１１】
0004694984.tif 000012
【００３９】
次に、得られた「ｔＮ」から新たな目標値「ｔＮ＋１」を予測する問題を考える。「ｔＮ＋１」の条件付き分布Ｐ（ｔＮ＋１|ｔＮ）は、ベイズの定理から次式（式１２）に示すように表される。
【００４０】
【数１２】
0004694984.tif 000013
【００４１】
式１２における右辺の分子は仮定より正規分布である。また、「ｔＮ」はすでに得られているので、分母は定数となる。よって、分布Ｐ（ｔＮ＋１|ｔＮ）は、次式（式１３）のように表すことができる。なお、式１３において、「ＣＮ＋１」は、「ｔ１」から「ｔＮ＋１」までの分散共分散行列である。
【００４２】
【数１３】
0004694984.tif 000014
【００４３】
ここで、「ＣＮ＋１」と「ＣＮ＋１−１」を次式（式１４）のように分割する。
【００４４】
【数１４】
0004694984.tif 000015
【００４５】
分割逆行列の一般公式により、「ＣＮ＋１−１」の各要素は次式（式１５）により「ＣＮ−１」で表すことができ、「ＣＮ＋１−１」を「ＣＮ−１」を用いて表すことができる。
【００４６】
【数１５】
0004694984.tif 000016
【００４７】
式１４および式１５から、式１３の分布における平均「μＮ＋１」と分散「σ２Ｎ＋１」は、次式（式１６）のようになる。
【００４８】
【数１６】
0004694984.tif 000017
【００４９】
上記のように、ベイズモデルにおける線形回帰と、近年研究が活発なカーネルマシンとの関連性からＧＰが定義される。
【００５０】
上述したように、式１０に示す共分散関数、すなわちカーネル関数の形式が決定すれば、ＧＰにおける予測分布は行列・ベクトル演算により簡単に求めることができる。よって、ＧＰにおける予測には共分散関数が重要な役割を果たす。しかし、このままでは共分散関数におけるパラメータ（例えば，式１０における「ω」）を決定することができない。ＧＰでは共分散関数のパラメータもさらに確率変数とみなし、事前分布を設定してベイズの定理より事後分布を求めることでパラメータを決定する。
【００５１】
上述したモデルは全て共分散関数のパラメータが決定した下での議論であったため、これらは共分散関数のパラメータの確率分布に依存することになる。このように、通常のパラメータよりもさらに上位に設定されたパラメータを「ハイパーパラメータ」といい、これらをベイズの枠組みでとらえることを「階層ベイズモデル」という。
【００５２】
次に、ハイパーパラメータを決定し、予測を行う方法について述べる。
ここでは、ハイパーパラメータθ=｛θｋ｝, （ｋ＝１,…,Ｋ）に事前分布Ｐ（θ）を設定した場合の予測分布の算出法について述べる。この場合、予測分布は、ベイズの定理から次式（式１７）の様になる。なお、式１７において、「Ｄ」は与えられたデータである。
【００５３】
【数１７】
0004694984.tif 000018
【００５４】
さらに、Ｐ（θ｜Ｄ）については、ベイズの定理より次式（式１８）が成り立つ。
【００５５】
【数１８】
0004694984.tif 000019
【００５６】
ここで、Ｐ（θ｜Ｄ）は、ハイパーパラメータでデータをどれだけ尤もらしく表現できるかを表す量（ハイパーパラメータ周辺尤度）であり、「ＡＢＩＣ」や「ｅｖｉｄｅｎｃｅ」と呼ばれることもある。この項の対数は次式（式１９）で表される。なお、式１９において、「det」は、行列の行列式である。
【００５７】
【数１９】
0004694984.tif 000020
【００５８】
式１９から分かるように、式１７の積分は非常に複雑となり、解析的に行うことは不可能である。そのため、近似計算が必要となる。以下、近似を行うための２つの手法について説明する。
【００５９】
「ＭＡＰ(Maximum A Posteriori)推定」
先ず、ＭＡＰ推定について説明する。ＭＡＰ推定とは、事後確率が最大となる値を点推定する手法である。すなわち、ＭＡＰ推定値を「θMAP」とすると、Ｐ(ｔＮ＋１｜ｘＮ＋１,Ｄ)をＰ(ｔＮ＋１｜ｘＮ＋１,Ｄ,θMAP)とする近似を行う。「θ」を固定するので、予測分布は正規分布となる。これを行うためには式１９とハイパーパラメータ事前分布の対数を用いて式１８の最大化を行うようにすればよい。式１９の偏導関数は次式（式２０）で与えられる。なお、式２０において、Ｔｒは行列のトレースである。
【００６０】
【数２０】
0004694984.tif 000021
【００６１】
よって、ｌｎＰ（θ）の偏導関数を求めることができれば、共役勾配法などの最適化法を用いて式１８の最大化を行うことができるようになる。
【００６２】
「ハイブリッドモンテカルロによるサンプリング」
次に、ハイブリッドモンテカルロによるサンプリングについて説明する。式１７の予測分布を近似するもう１つの方法として、マルコフ連鎖モンテカルロ法を用いたサンプリングにより期待値などを計算する方法がある。高い確率を持つ領域から効率よくサンプルをとることができれば、サンプルが少数であっても十分に精度の高い近似値が得られる。
【００６３】
しかし、この領域は高次元領域におけるごく小さな領域であり、単純なランダムウォークではこの領域からサンプルをとること自体が難しくなる。本例では、式１９のように導関数を計算できる。この場合は、マルコフ連鎖モンテカルロ法の１種であるハイブリッドモンテカルロ（以下、「ＨＭＣ」という。）を用いるのが効果的である。「ＨＭＣ」は，Metropolis基準とハミルトニアンダイナミカルシステムのシミュレートを組み合わせたハイブリッド手法である。一般に、モンテカルロ法では事後分布からのサンプルを得て、それから次式（式２１）を用いて式１７の近似を行う。
【００６４】
【数２１】
0004694984.tif 000022
【００６５】
なお、式２１において、「Ｓ」は予測に用いるサンプル数であり、「θ（Ｓ）」はｓ番目のθのサンプルである。
【００６６】
「ＨＭＣ」では勾配情報を高い確率をもつ領域を見つけ出すために進むべき方向の決定に用いる。一方、モーメント変数ｐ＝｛ｐｋ｝，（ｋ＝１，…，Ｋ）を導入し、探索する運動の様子を決定する。モーメント変数より、運動エネルギーｋｉ（ｐ）が次式（式２２）で定まる。但し、式２２において、「ｍ＝｛ｍｋ｝，（ｋ＝１，…，Ｋ）」は、仮想的質量を意味する。
【００６７】
【数２２】
0004694984.tif 000023
【００６８】
そして、この運動エネルギーの確率分布Ｐ（ｐ）は次式（式２３）に示すように、ギブス分布でモデル化できる。なお、式２３において、「Ｚｋ」は規格化定数である。
【００６９】
【数２３】
0004694984.tif 000024
【００７０】
また、式２２における「ｍ」は便宜上全て「１」と置いてよく、「１」と置いた場合には式２３の確率分布は標準正規分布となる。一方、「θ」を位置変数に対応させ、位置エネルギーをＥ（θ）＝−ｌｎＰ（θ｜Ｄ）とし、その確率分布をＰ（θ｜Ｄ）に比例するｅｘｐ（−Ｅ（θ））とする。すると、位置エネルギーと運動エネルギーから、次式（式２４）のように、総エネルギーであるハミルトニアンが決定できる。
【００７１】
【数２４】
0004694984.tif 000025
【００７２】
このハミルトニアンを元に仮想時間「τ」を用いたダイナミカルシステムのシミュレーションを行う。このシステムは次式（式２５）の微分方程式により支配される。
【００７３】
【数２５】
0004694984.tif 000026
【００７４】
しかし、この微分方程式は一般に複雑な方程式となり、直接これを解くことはできない。通常は、次式（式２６）のleapfrogと呼ばれる離散近似法がとられている。なお、式２６において、「ε」はステップサイズである。
【００７５】
【数２６】
0004694984.tif 000027
【００７６】
これをＬ回繰り返し、現在の状態（θ，ｐ）から時間ετ後の提案状態（θ＊，ｐ＊）を作り出す。これをMetropolis基準に従ってこの状態の受容・棄却を決定する。Metropolis基準は次式（式２７）の確率により提案状態を受容する方法である。
【００７７】
【数２７】
0004694984.tif 000028
【００７８】
さらに、今回は目的の分布により早く収束させるため、Horowitsによるpersistenceパラメータα（０≦α≦１）を導入する。通常のＨＭＣではＬ回の反復を行うごとに新しいモーメント変数のサンプリングを行うが、パラメータαを導入する場合は次式（式２８）によりモーメントを更新する。なお、式２８において、「ν」はＮ（０，ｍ）から得られる確率変数である。
【００７９】
【数２８】
0004694984.tif 000029
【００８０】
このパラメータαの導入により、θ空間における次の状態への探索する方向が類似することになり、従ってランダムウォークを避けることに貢献している。得られたサンプルから、予測平均の推定値の事後分布の平均と分散を次式（式２９，式３０）により推定する。
【００８１】
【数２９】
0004694984.tif 000030
【００８２】
【数３０】
0004694984.tif 000031
【００８３】
ｓ番目の予測分散のサンプルアルゴリズムの初期は目的の分布に収束していないと考えられるため、通常は初期サンプルの何％かは予測に用いずに捨てる。これをburn-inという。
【００８４】
以下、本例の電力負荷予測のための処理方法について説明する。
本例では、ＧＰを用いて短期電力負荷予測を行い、負荷の不確定性をエラーバーにより表現可能な電力負荷予測方法を提供する。電力負荷予測問題では季節変動などの不確定要素が含まれるため、ある程度の誤差は必ず発生する。その誤差を抑えることも確かに重要だが、観測値が取り得るであろう範囲を予測して不確定性を表現し、観測値の傾向をつかむことも重要であると思われる。また、これらの情報は運用者にとっても有用な情報となりえると考えられる。ＧＰのモデル構築には、ＭＡＰ推定とＨＭＣを用いる。それぞれの手法の流れを以下に示す。
【００８５】
「ＭＡＰ推定の手順」
図２は、ＭＡＰ推定による電力負荷予測処理の例を示すフローチャートである。電力負荷予測処理において、演算装置２は、先ず、ハイパーパラメータの事前分布の形式とそのパラメータ、及び共分散関数の形式を決定する（ステップＳ１１）。次いで、演算装置２は、適当な初期値から出発して式１８を最大化し、θＭＡＰを得る（ステップＳ１２）。さらに、演算装置２は、θＭＡＰを用いて予測分布の平均と分散を式１６により求める（ステップＳ１３）。
【００８６】
「ＨＭＣの手順」
図３は、ＨＭＣによる電力負荷予測処理の例を示すフローチャートである。電力負荷予測処理において、演算装置２は、先ず、ハイパーパラメータの事前分布の形式とそのパラメータ、及び共分散関数の形式を決定する（ステップＳ２１ａ）。また、演算装置２は、ＨＭＣの反復回数、leapfrogの反復回数、ステップサイズ、persistence，burn-inの割合を決定する（ステップＳ２１ｂ）。なお、θの初期値については、乱数などの適当な方法で決定する。また，演算装置２は、モーメント変数ｐの初期値を標準正規乱数により発生させる（ステップＳ２２）。そして、演算装置２は、ステップＳ２１ｂで決定したＨＭＣの反復回数を満たすまで、以下（ステップＳ２３〜Ｓ２８）を実行する。
【００８７】
ステップＳ２３〜Ｓ２８では、演算装置２は、先ず、現在の状態を（θ，ｐ）とし，ステップサイズ「ε」のleapfrogをＬ回行い、提案状態（θ＊，ｐ＊）を作成する（ステップＳ２３）。次いで、演算装置２は、−（Ｈ（θ＊，ｐ＊）−Ｈ（θ，ｐ））＜０であれば提案状態を受容し、（θ，ｐ）←（θ＊，ｐ＊）とする（ステップＳ２４）。
【００８８】
また、演算装置２は、ｒｎｄ（）＜ｅｓｐ［−（Ｈ（θ＊，ｐ＊）−Ｈ（θ，ｐ））］であれば提案状態を受容し、（θ，ｐ）←（θ＊，ｐ＊）とする（ステップＳ２５）。そうでなければ、演算装置２は、提案状態を棄却し、（θ，ｐ）←（θ，−ｐ）とする（ステップＳ２６）。なお、ｒｎｄ（）は、区間［０，１］の一様乱数とする。
【００８９】
次いで、演算装置２は、式２８によりモーメントを更新する（ステップＳ２７）。そして、演算装置２は、burn-inを行わないのであれば、θをサンプルとして採用する（ステップＳ２８）。
【００９０】
そして、演算装置２は、ステップＳ２１ｂで決定したＨＭＣの反復回数を満たしていれば（ステップＳ２９）、式２９と式３０を用いて事後分布の予測平均と予測分散を評価する（ステップＳ３０）。
【００９１】
以下、本実施形態によるＧＰを用いて短期電力負荷予測を行い負荷の不確定性をエラーバーにより表現する電力価格予測のシミュレーション結果について説明する。
【００９２】
「シミュレーション条件」
先ず、シミュレーション条件について説明する。本シミュレーションでは、サンプルデータとして、平成１１年から平成１４年の夏季（６月〜９月）の平日データの最大電力負荷を用い、そのうち平成１４年のデータをテストデータとして用いた。ただし、お盆などの特異日のデータはデータから除いてある。図４は、使用する入力変数を示す図である。なお、曜日については、図５に示すようにビット化した。よって、入力次元は１１次元である。
【００９３】
また、ＧＰに用いる共分散関数として、次式（式３１）の関数を用いた。なお、式３１において、「Ｇ」は入力次元であり、ｌｎｖ，ｌｎｗ，ｌｎａはハイパーパラメータであり、「δｉｊ」はクロネッカーのデルタである。
【００９４】
【数３１】
0004694984.tif 000032
【００９５】
式３１の共分散関数は「ｖ」，「ｗ」，「ａ」が正の場合に対称半正定値となる。よって、ｌｎｖ，ｌｎｗ，ｌｎａをハイパーパラメータとすることで「ｖ」，「ｗ」，「ａ」の正値性が自然に満たされるよう設計してある。
【００９６】
また、各入力次元に独立に重み「ｗ」を設けてある。こうすることで，予測に高く寄与した変数の「ｗ」は大きな値を、一方で予測にあまり寄与しない変数の「ｗ」は小さな値をとることになる。すなわち、予測に柔軟性を持たせることができる。また，曜日を前述の様にビット化したため、どの曜日が予測に高く寄与したかの知見を得られることが期待できる。
【００９７】
上記の「ｗ」の事前分布はそれぞれ独立なガンマ分布とし、「ｖ」と「ａ」の事前分布はそれぞれ独立な正規分布とした。ガンマ分布は式３２で表され、正規分布は式３３で表される。なお、「α」は形状母数、「β」は尺度母数、「μ」は平均、「σ」は標準偏差、「Γ（）」はガンマ関数を示す。
【００９８】
【数３２】
0004694984.tif 000033
【００９９】
【数３３】
0004694984.tif 000034
【０１００】
図６は、ＧＰのハイパーパラメータの事前分布に関する母数を示す説明図である。また、図７は、ＨＭＣを用いたＧＰ（以下、ＧＰ（ＨＭＣ）という。）におけるパラメータを示す説明図である。なお、ＭＡＰ推定を用いたＧＰ（以下、ＧＰ（ＭＡＰ）という。）の最適化には、共役勾配法を用いた。また、誤差を計算する際には、ＧＰ（ＭＡＰ）・ＧＰ（ＨＭＣ）における誤差の計算には得られた分布の平均値を用いた。
【０１０１】
一方、比較手法としては中間層が１層の多層パーセプトロン（以下、ＭＬＰという。）とラジアル基底関数ネットワーク（以下、ＲＢＦＮという。）を用いた。ＭＬＰの学習アルゴリズムは誤差逆伝播法を用い、ＲＢＦＮの学習アルゴリズムは、基底関数の中心と幅の初期値をまずK-meansで与え、その更新と結合加重の計算を交互に行った。中心と幅の更新は共役勾配法を、結合加重の計算には最小二乗法を用いた。また、ＭＬＰ・ＲＢＦＮ共に結合加重の二乗和をペナルティとして評価関数に加える正則化を施した（ペナルティ係数：λ）。これら比較手法のパラメータを図８に示す。ただし、入出力変数は、全ての手法において平均０、分散１に標準化して用いた。
【０１０２】
「シミュレーション結果」
上記のシミュレーション条件によるシミュレーション結果について説明する。先ず、図９に各手法の平均誤差と最大誤差をそれぞれ示す。図９により、ＧＰ（ＭＡＰ）は比較手法であるＭＬＰとＲＢＦＮと比べて、平均誤差において0.4063％と0.3512％、最大誤差において1.2796％と0.9535％を上回る結果を残した。
【０１０３】
同様に、ＧＰ（ＨＭＣ）は平均誤差において0.4086％と0.3535％、最大誤差において1.3038％と1.2437％上回る結果を示した。よって、ＧＰがＭＬＰとＲＢＦＮよりも短期電力負荷予測において高い予測能力を持つことが分かる。また、ＧＰに注目すると、平均誤差・最大誤差共に若干だがＧＰ（ＨＭＣ）がＧＰ（ＭＡＰ）を上回る結果となった。
【０１０４】
また、ＧＰ（ＭＡＰ）とＧＰ（ＨＭＣ）における平成１４年度全日の予測結果と不確定性の指標であるエラーバー（３σ）を図１０と図１１にそれぞれ示す。これらの図に示されている様に、エラーバーを用いて不確定性を表現可能な点が本例の特徴である。
【０１０５】
誤差の報告で述べたことと同様に、１日ごとの予測結果も両手法の差異は若干であった。エラーバーに注目すると、全日において両手法とも観測値が±３σ内に入っていることが見て取れる。実際、観測値が±ｎσに入っている割合は、±σ内が73.42％、±2σが93.67％、±３σが100.00％であった。これらの結果より、ＧＰは翌日最大電力の予測結果だけではなく、その不確定性の算出に関しても高い能力を示すことがわかる。
【０１０６】
最後に、 ＧＰ（ＭＡＰ）におけるハイパーパラメータ「w」の平均値を図１２に示す。同図より、ＧＰ（ＭＡＰ）においては最高気温・平均気温・不快指数・当日最大電力が予測に高く寄与していることが読み取れる。一方で，ビット化した曜日に注目すると，どの曜日もハイパーパラメータ「w」の値は小さいものの月曜日と金曜日が他の曜日よりも高い値を示した。よって、曜日の特色として週明けと週末は他の曜日よりも重要な情報をもっているといえる。
【０１０７】
以上より、ＧＰは従来法と比べて予測能力が高いだけでなく、予測値の上下限値であるエラーバーに関し付加的な情報が得られる点でも優れていることが示された。
【０１０８】
上述したように、本実施形態では、ＧＰを適用して短期電力負荷予測を行うようにしたので、以下のような効果が得られる。すなわち、ＧＰにより共分散関数において各次元に重みを決定し，それをベイズの定理から決定したので、従来よりも柔軟な予測を行うことができる。また、ＧＰによる短期電力負荷予測により電力負荷予測を行うこととしたので、実データと比較した結果、平均誤差及び最大誤差を共に削減することができ、精度の高い予測を行うことができる。さらに、予測を点推定ではなく予測分布を用いて行うこととしたので、エラーバーの算出が可能となり、短期電力負荷予測における不確定性を表現することができるようになった。
【図面の簡単な説明】
【０１０９】
【図１】本発明の一実施形態に係る電力負荷予測システムのブロック図である。
【図２】ＭＡＰ推定による電力負荷予測処理の例を示すフローチャートである。
【図３】ＨＭＣによる電力負荷予測処理の例を示すフローチャートである。
【図４】入力変数の例を示す説明図である。
【図５】曜日についてビット化した情報の例を示す説明図である。
【図６】ＧＰのハイパーパラメータの事前分布に関する母数を示す説明図である。
【図７】ＧＰ（ＨＭＣ）におけるパラメータを示す説明図である。
【図８】比較手法のパラメータの例を示す説明図である。
【図９】各手法の平均誤差と最大誤差を示す説明図である。
【図１０】ＧＰ（ＭＡＰ）における予測結果とエラーバーを示す説明図である。
【図１１】ＧＰ（ＨＭＣ）における予測結果とエラーバーを示す説明図である。
【図１２】ＧＰ（ＭＡＰ）におけるハイパーパラメータの平均値を示す説明図である。
【符号の説明】
【０１１０】
１…入力装置
２…演算装置
３…記憶装置
４…出力装置
【図１】
0004694984.tif 000037
【図２】
0004694984.tif 000038
【図３】
0004694984.tif 000039
【図４】
0004694984.tif 000040
【図５】
0004694984.tif 000041
【図６】
0004694984.tif 000042
【図７】
0004694984.tif 000043
【図８】
0004694984.tif 000044
【図９】
0004694984.tif 000045
【図１０】
0004694984.tif 000046
【図１１】
0004694984.tif 000047
【図１２】
0004694984.tif 000048