(57)【要約】【課題】観測不可能な外乱による影響が生じた場合であっても、外乱の事前情報を表す標準偏差などの統計値を利用して、予測対象の出力の予測精度が低下することを抑制可能な予測装置を提供する。【解決手段】予測対象２の将来における出力を予測する予測装置３は、プロセッサ３１と、前記プロセッサ３１に接続され、前記予測対象２の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置３０と、を備える。前記プロセッサ３１は、過去に記憶された複数の前記入力計測値、及び複数の前記出力計測値から、移動平均フィルタを用いて前記入力に対する第１係数を同定する同定処理と、前記入力計測値、前記出力計測値、及び前記第１係数からなる予測モデルに基づいて、前記予測対象２の将来における出力を予測する予測処理と、を実行する。【選択図】図１Ａ
2020194535.tif 000002
【特許請求の範囲】
【請求項１】
予測対象の将来における出力を予測する予測装置であって、
プロセッサと、
前記プロセッサに接続され、前記予測対象の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置と、
を備え、
前記プロセッサは、
過去に記憶された複数の前記入力計測値、及び複数の前記出力計測値から、移動平均フィルタを用いて前記入力に対する第１係数を同定する同定処理と、
前記入力計測値、前記出力計測値、及び前記第１係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測する予測処理と、
を実行し、
前記同定処理において、観測雑音の共分散行列と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列とを用いる、
予測装置。
【請求項２】
予測対象の将来における出力を予測する予測装置であって、
プロセッサと、
前記プロセッサに接続され、前記予測対象の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置と、
を備え、
前記プロセッサは、
過去に記憶された複数の前記入力計測値、及び複数の前記出力計測値から、自己回帰移動平均フィルタを用いて前記入力に対する第１係数と、前記出力に対する第２係数とを同定する同定処理と、
前記入力計測値、前記出力計測値、前記第１係数、及び前記第２係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測する予測処理と、
を実行し、
前記同定処理において、観測雑音の共分散行列と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、前記第２係数に係るカーネル行列で重み付けした観測雑音の共分散行列とを用いる、
予測装置。
【請求項３】
予測対象の将来における出力を予測する予測装置であって、
プロセッサと、
前記プロセッサに接続され、前記予測対象の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置と、
を備え、
前記プロセッサは、過去に記憶された複数の前記入力計測値、及び複数の前記出力計測値から、無限インパルス応答フィルタを用いて前記入力に対する第１係数と、前記出力に対する第２係数とを同定する同定処理と、
前記入力計測値、前記出力計測値、前記第１係数、及び前記第２係数からなる無限インパルス応答フィルタの形式による予測モデルに基づいて、前記予測対象の将来における出力を予測する予測処理と、
を実行し、
前記同定処理において、観測雑音の共分散行列、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列を用いる、
予測装置。
【請求項４】
前記プロセッサは、前記予測処理において、現在から所定時間後の前記予測対象の出力を予測する、
請求項１から３の何れか一項に記載の予測装置。
【請求項５】
前記記録装置は、前記予測対象の複数種類の入力計測値、及び複数種類の出力計測値を記憶し、
前記プロセッサは、前記同定処理において、複数種類の前記入力それぞれに対する複数の前記第１係数と、複数種類の前記出力それぞれに対する複数の前記第２係数と、を同定する、
請求項２に記載の予測装置。
【請求項６】
請求項１から５の何れか一項に記載の予測装置と、
前記予測装置と通信可能に接続され、当該予測装置から受信した前記予測対象の出力の予測値に基づいて、当該予測対象の入力を調整する制御装置と、
を備える予測システム。
【請求項７】
前記予測対象は、電力系統に電力の供給を行う電源であり、
前記制御装置は、前記予測値に基づいて、前記電源が有するタービン装置の加減弁の開度を調整する、
請求項６に記載の予測システム。
【請求項８】
予測対象の将来における出力を予測する予測方法であって、
過去の複数の入力の計測値である入力計測値、及び複数の出力の計測値である出力計測値から、移動平均フィルタを用いて前記入力に対する第１係数を同定するステップと、
前記入力計測値、前記出力計測値、及び前記第１係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測するステップと、
を有し、
前記第１係数を同定するステップにおいて、観測雑音の共分散行列と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列とを用いる、
予測方法。
【請求項９】
予測対象の将来における出力を予測する予測方法であって、
過去の複数の入力の計測値である入力計測値、及び複数の出力の計測値である出力計測値から、自己回帰移動平均フィルタを用いて前記入力に対する第１係数と、前記出力に対する第２係数とを同定するステップと、
前記入力計測値、前記出力計測値、前記第１係数、及び前記第２係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測するステップと、
を有し、
前記第１係数及び前記第２係数を同定するステップにおいて、観測雑音の共分散行列と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、前記第２係数に係るカーネル行列で重み付けした観測雑音の共分散行列とを用いる、
予測方法。
【請求項１０】
予測対象の将来における出力を予測する予測方法であって、
過去の複数の入力の計測値である入力計測値、及び複数の出力の計測値である出力計測値から、無限インパルス応答フィルタを用いて前記入力に対する第１係数と、前記出力に対する第２係数とを同定するステップと、
前記入力計測値、前記出力計測値、前記第１係数、及び前記第２係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測するステップと、
を有し、
前記第１係数及び前記第２係数を同定するステップにおいて、観測雑音の共分散行列と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、を用いる、
予測方法。
【請求項１１】
プロセッサと、前記プロセッサに接続され、予測対象の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置と、を備える予測装置のコンピュータを機能させるプログラムであって、前記コンピュータに、
過去に記憶された複数の前記入力計測値、及び複数の前記出力計測値から、移動平均フィルタを用いて前記入力に対する第１係数を同定するステップと、
前記入力計測値、前記出力計測値、及び前記第１係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測するステップと、
を実行させ、
前記第１係数を同定するステップにおいて、観測雑音の共分散行列と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列とを用いる、
プログラム。
【請求項１２】
プロセッサと、前記プロセッサに接続され、予測対象の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置と、を備える予測装置のコンピュータを機能させるプログラムであって、前記コンピュータに、
過去に記憶された複数の前記入力計測値、及び複数の前記出力計測値から、自己回帰移動平均フィルタを用いて前記入力に対する第１係数と、前記出力に対する第２係数とを同定するステップと、
前記入力計測値、前記出力計測値、前記第１係数、及び前記第２係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測するステップと、
を実行させ、
前記第１係数及び前記第２係数を同定するステップにおいて、観測雑音の共分散行列と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、前記第２係数に係るカーネル行列で重み付けした観測雑音の共分散行列とを用いる、
プログラム。
【請求項１３】
プロセッサと、前記プロセッサに接続され、予測対象の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置と、を備える予測装置のコンピュータを機能させるプログラムであって、前記コンピュータに、
過去の複数の入力の計測値である入力計測値、及び複数の出力の計測値である出力計測値から、無限インパルス応答フィルタを用いて前記入力に対する第１係数と、前記出力に対する第２係数とを同定するステップと、
前記入力計測値、前記出力計測値、前記第１係数、及び前記第２係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測するステップと、
を実行させ、
前記第１係数及び前記第２係数を同定するステップにおいて、観測雑音の共分散行列と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、を用いる、
プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、予測装置、予測システム、予測方法、及びプログラムに関する。
【背景技術】
【０００２】
電力事業、ガス事業等のエネルギー事業分野、通信事業分野、タクシー、配送業投の運送事業分野等の様々な分野において、消費者の需要に応じて設備を稼働させるため、及び資源を適切に配分するために、将来の需要量を予測することが行われている。将来の需要量の予測には、予測対象の数値モデル（予測モデル）が使われる。
【０００３】
しかしながら、予測対象の特性は一定ではなく、経時的に変動する。したがって、ある時点で定めた数値モデルを使って算出した予測値と、実際に観測された値との間には乖離が生じる可能性がある。このため、予測精度を向上するために、予測モデルを時間の経過とともに更新することが考えられている（たとえば、特許文献１を参照）。
【０００４】
また、予測モデルの更新は、過去の実際の観測値に基づいて行われる。予測モデルの更新に使われるデータ（たとえば、入力及び出力の観測値）は、実用上はより少ないことが望まれる。少ないデータ（短期間の入力及び出力を示すデータ）で予測モデルを更新することが出来れば、予測対象の特性が短時間のうちに変動したとしても、予測モデルで変動をとらえることができるからである。たとえば、予測対象の特性が１０分間かけて変動する場合、過去１分間分の観測値からモデル更新できるならば、変動に合わせて予測モデルを更新することができる。一方、モデル更新に過去６０分間分の観測値が必要である場合、予測モデルが更新されるのは、予測対象が変動してから数十分後となるので、予測が現実と違ってしまう可能性がある。
【０００５】
近年では、少ないデータで予測を可能とする技術として、カーネル型システム同定法が考えられている（たとえば、非特許文献１及び非特許文献２を参照）。カーネル型システム同定法では、予測対象のインパルス応答をベイズ推定により同定して予測モデルを更新する。ここで、予測対象の出力は移動平均フィルタ（ＭＡフィルタ）の形式で表現される。具体的には、予測対象の入力をｕｉ（ｉ＝ｔ，ｔ−１，ｔ−２，…）、出力をｙｉ（ｉ＝ｔ，ｔ−１，ｔ−２，…）としたとき、出力ｙを、入力ｕの荷重付きの移動平均と、観測雑音との和で表現する。荷重付きの移動平均の荷重係数｛ａ１，ａ２，…，ａｎ｝は、予測対象の有限インパルス応答に相当する。観測雑音は平均がゼロ、標準偏差がσｗの正規分布で表されると仮定する。なお、標準偏差σは予め与えられているとする。移動平均フィルタによる予測モデルでは、荷重係数｛ａ１，ａ２，…，ａｎ｝の値を実際の入力及び出力の観測値ｙ、ｕに基づいて定める。カーネル型システム同定法では、荷重係数｛ａ１，ａ２，…，ａｎ｝に対し平均がゼロ、共分散行列Ｋ（ｎ×ｎ）の事前分布と、観測雑音に対し平均がゼロ、共分散行列σｗ２Ｉ（ｎ×ｎ）の事前分布と、を導入して、｛ａ１，ａ２，…，ａｎ｝の値を定める。このＫはカーネル行列と呼ばれ、予測対象に関する事前情報が組み込まれている。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２０１８−１６３５１５号公報
【非特許文献】
【０００７】
【非特許文献１】藤本悠介，杉江俊治：カーネル型システム同定のための入力設計に関する一考察，第59回自動制御連合講演会，北九州，2016.11.10-12, pp. 448--449 (2016.11.10)
【非特許文献２】G. Pillonetto, A. Chiuso and G. De Nicolao, “Prediction error identification of linear systems: A nonparametric Gaussian regression approach,” Automatica, 47(2), 291-305, 2011.
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、一般に、予測対象は複数の入力の影響を受ける。上述のような移動平均フィルタを使う予測モデルでは、特定の一つの入力（入力ｕ）に対して、予測対象の応答（出力ｙ）を予測する。言い換えれば、予測対象の応答は、特定の一つの入力にのみ影響されるとみなしている。このため、計測できない未知の入力（外乱）が加わると、応答は予測モデルで表すことができない。このため、従来の技術では、観測不可能な外乱による影響が生じると、予測精度が低下する可能性があった。
【０００９】
本発明は、このような課題に鑑みてなされたものであって、観測不可能な外乱による影響が生じた場合であっても、外乱の事前情報を表す標準偏差などの統計値を利用して、予測対象の出力の予測精度が低下することを抑制可能な予測装置、予測システム、予測方法、及びプログラムを提供する。
【課題を解決するための手段】
【００１０】
上記課題を解決するため、本発明は以下の手段を採用している。
本発明の第１の態様によれば、予測対象の将来における出力を予測する予測装置は、プロセッサと、前記プロセッサに接続され、前記予測対象の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置と、を備える。前記プロセッサは、過去に記憶された複数の前記入力計測値、及び複数の前記出力計測値から、移動平均フィルタを用いて前記入力に対する第１係数を同定する同定処理と、前記入力計測値、前記出力計測値、及び前記第１係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測する予測処理と、を実行する。前記同定処理については、出力ｙの予測値の事前分布の評価において、観測雑音の共分散行列σｗ２Ｉ（ｎ×ｎ）と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、を用いることを特徴とする。
【００１１】
本発明の第２の態様によれば、予測対象の将来における出力を予測する予測装置は、プロセッサと、前記プロセッサに接続され、前記予測対象の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置と、を備える。前記プロセッサは、過去に記憶された複数の前記入力計測値、及び複数の前記出力計測値から、自己回帰移動平均フィルタを用いて前記入力に対する第１係数と、前記出力に対する第２係数とを同定する同定処理と、前記入力計測値、前記出力計測値、前記第１係数、及び前記第２係数からなる自己回帰移動平均フィルタの形式による予測モデルに基づいて、前記予測対象の将来における出力を予測する予測処理と、を実行する。前記同定処理については、出力ｙの予測値の事前分布の評価において、観測雑音の共分散行列σｗ２Ｉ（ｎ×ｎ）に加え、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、さらに前記第２係数に係るカーネル行列で重み付けした観測雑音の共分散行列と、を用いることを特徴とする。
【００１２】
本発明の第３の態様によれば、予測対象の将来における出力を予測する予測装置は、プロセッサと、前記プロセッサに接続され、前記予測対象の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置と、を備える。前記プロセッサは、過去に記憶された複数の前記入力計測値、及び複数の前記出力計測値から、無限インパルス応答フィルタを用いて前記入力に対する第１係数と、前記出力に対する第２係数とを同定する同定処理と、前記入力計測値、前記出力計測値、前記第１係数、及び前記第２係数からなる無限インパルス応答フィルタの形式による予測モデルに基づいて、前記予測対象の将来における出力を予測する予測処理と、を実行する。前記同定処理については、出力ｙの予測値の事前分布の評価において、観測雑音の共分散行列σｗ２Ｉ（ｎ×ｎ）に加え、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列を用いることを特徴とする。
【００１３】
本発明の第４の態様によれば、前記プロセッサは、前記予測処理において、現在から所定時間後の前記予測対象の出力を予測する。
【００１４】
本発明の第５の態様によれば、前記記録装置は、前記予測対象の複数種類の入力計測値、及び複数種類の出力計測値を記憶する。前記プロセッサは、前記同定処理において、複数種類の前記入力それぞれに対する複数の前記第１係数と、複数種類の前記出力それぞれに対する複数の前記第２係数と、を同定する。
【００１５】
本発明の第６の態様によれば、予測システムは、第１から第５の何れか一の態様に記載の予測装置と、前記予測装置と通信可能に接続され、当該予測装置から受信した前記予測対象の出力の予測値に基づいて、当該予測対象の入力を調整する制御装置と、を備える。
【００１６】
本発明の第７の態様によれば、前記予測対象は、電力系統に電力の供給を行う電源であり、前記制御装置は、前記予測値に基づいて、前記電源が有するタービン装置の加減弁の開度を調整する。
【００１７】
本発明の第８の態様によれば、予測対象の将来における出力を予測する予測方法は、過去の複数の入力の計測値である入力計測値、及び複数の出力の計測値である出力計測値から、移動平均フィルタを用いて前記入力に対する第１係数を同定するステップと、前記入力計測値、前記出力計測値、及び前記第１係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測するステップと、を有する。前記第１係数を同定するステップについては、出力ｙの予測値の事前分布の評価において、観測雑音の共分散行列σｗ２Ｉ（ｎ×ｎ）と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、を用いることを特徴とする。
【００１８】
本発明の第９の態様によれば、予測対象の将来における出力を予測する予測方法は、過去の複数の入力の計測値である入力計測値、及び複数の出力の計測値である出力計測値から、自己回帰移動平均フィルタを用いて前記入力に対する第１係数と、前記出力に対す第２係数とを同定するステップと、前記入力計測値、前記出力計測値、前記第１係数、及び前記第２係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測するステップと、を有する。前記第１係数及び前記第２係数を同定するステップについては、出力ｙの予測値の事前分布の評価において、観測雑音の共分散行列σｗ２Ｉ（ｎ×ｎ）に加え、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、さらに前記第２係数に係るカーネル行列で重み付けした観測雑音の共分散行列と、を用いることを特徴とする。
【００１９】
本発明の第１０の態様によれば、予測対象の将来における出力を予測する予測方法は、過去の複数の入力の計測値である入力計測値、及び複数の出力の計測値である出力計測値から、無限インパルス応答フィルタを用いて前記入力に対する第１係数と、前記出力に対する第２係数とを同定するステップと、前記入力計測値、前記出力計測値、前記第１係数、及び前記第２係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測するステップと、を有する。前記第１係数及び前記第２係数を同定するステップにおいて、観測雑音の共分散行列と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、を用いることを特徴とする。
【００２０】
本発明の第１１の態様によれば、プロセッサと、前記プロセッサに接続され、予測対象の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置と、を備える予測装置のコンピュータを機能させるプログラムは、前記コンピュータに、過去に記憶された複数の前記入力計測値、及び複数の前記出力計測値から、移動平均フィルタを用いて前記入力に対する第１係数を同定するステップと、前記入力計測値、前記出力計測値、及び前記第１係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測するステップと、を実行させる。前記第１係数を同定するステップについては、観測雑音の共分散行列σｗ２Ｉ（ｎ×ｎ）と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、を用いることを特徴とする。
【００２１】
本発明の第１２の態様によれば、プロセッサと、前記プロセッサに接続され、予測対象の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置と、を備える予測装置のコンピュータを機能させるプログラムは、前記コンピュータに、過去に記憶された複数の前記入力計測値、及び複数の前記出力計測値から、自己回帰移動平均フィルタを用いて前記入力に対する第１係数と、前記出力に対す第２係数とを同定するステップと、前記入力計測値、前記出力計測値、前記第１係数、及び前記第２係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測するステップと、を実行させる。前記第１係数及び前記第２係数を同定するステップについては、出力ｙの予測値の事前分布の評価において、観測雑音の共分散行列σｗ２Ｉ（ｎ×ｎ）に加え、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、さらに前記第２係数に係るカーネル行列で重み付けした観測雑音の共分散行列と、を用いることを特徴とする。
【００２２】
本発明の第１３の態様によれば、プロセッサと、前記プロセッサに接続され、予測対象の入力の計測値である入力計測値、及び出力の計測値である出力計測値を記憶する記録装置と、を備える予測装置のコンピュータを機能させるプログラムは、前記コンピュータに、過去の複数の入力の計測値である入力計測値、及び複数の出力の計測値である出力計測値から、無限インパルス応答フィルタを用いて前記入力に対する第１係数と、前記出力に対する第２係数とを同定するステップと、前記入力計測値、前記出力計測値、前記第１係数、及び前記第２係数からなる予測モデルに基づいて、前記予測対象の将来における出力を予測するステップと、を実行させる。前記第１係数及び前記第２係数を同定するステップにおいて、観測雑音の共分散行列と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、を用いることを特徴とする。
【発明の効果】
【００２３】
上述の少なくとも一の態様によれば、観測不可能な外乱による影響が生じた場合であっても、外乱の標準偏差についての事前情報を利用して、予測対象の出力の予測精度が低下することを抑制可能である。
【図面の簡単な説明】
【００２４】
【図１Ａ】本発明の第１の実施形態に係る予測システムの機能構成を示す図である。
【図１Ｂ】本発明の第２の実施形態に係る予測システムの機能構成を示す図である。
【図２】本発明の第２の実施形態に係る予測システムの実用例を示す図である。
【図３】本発明の第３の実施形態に係る予測システムの機能構成を示す図である。
【図４】本発明の第４の実施形態に係る予測システムの機能構成を示す図である。
【図５】本発明の第５の実施形態に係る予測システムの機能構成を示す図である。
【図６】本発明の少なくとも一実施形態に係る予測装置、及び制御装置のハードウェア構成の一例を示す図である。
【発明を実施するための形態】
【００２５】
＜第１の実施形態＞
以下、本発明の第１の実施形態に係る予測装置、及び当該予測装置を具備する予測システムについて、図１Ａを参照しながら説明する。
【００２６】
（機能構成）
図１Ａは、本発明の第１の実施形態に係る予測システムの機能構成を示す図である。
図１Ａに示すように、予測システム１は、予測対象２の将来における出力を予測する予測装置３を備えている。
【００２７】
本実施形態に係る予測装置３は、カーネル型システム同定法の技術を利用して、将来の出力を予測する。また、予測装置３は、サーバ、パーソナルコンピュータ等のコンピュータを用いて構成され、記録装置３０と、プロセッサ３１とを備えている。
【００２８】
記録装置３０は、プロセッサ３１に接続され、予測対象２から所定周期ごとに受信する入力の計測値（以下、「入力計測値ｕ」と記載する）、及び出力の計測値（以下、「出力計測値ｙ」と記載する）を記憶する。
【００２９】
プロセッサ３１は、予測装置３の動作全体を司る。プロセッサ３１は、所定のプログラムに従って動作することにより、同定部３１０、予測部３１１としての機能を発揮する。
【００３０】
同定部３１０は、過去に記憶された複数の入力計測値ｕ、及び複数の出力計測値ｙから、移動平均フィルタ（ＭＡフィルタ）を用いて予測対象２の入力に対する第１係数ａを同定する同定処理Ｓ１を実行する。
【００３１】
予測部３１１は、入力計測値ｕ、出力計測値ｙ、第１係数ａからなる予測モデルに基づいて、予測対象２の将来における出力を予測する予測処理Ｓ２を実行する。
【００３２】
（カーネル型システム同定法について）
ここで、従来のカーネル型システム同定法における同定処理について説明する。従来のカーネル型システム同定法において、予測対象の時刻毎の入力計測値をｕｉ（ｉ＝ｔ，ｔ−１，ｔ−２，…）、出力計測値をｙｉ（ｉ＝ｔ，ｔ−１，ｔ−２，…）とすると、予測対象の１ステップ先（時刻ｔ＋１）における出力ｙｔ＋１は、以下の式（１）のように、移動平均フィルタ（ＭＡフィルタ）の形式で表現される。なお、明細書中における「＾ｙ」との表記は、以下に示す図、式中において「ｙ」にハット記号「＾」が付された表記に対応する。同様に、明細書中における「＾Ｋ」、「＾ａ」との表記は、以下に示す図、式中において「Ｋ」、「ａ」にハット記号「＾」が付された表記に対応する。
【００３３】
【数１】
2020194535.tif 000003
【００３４】
上式（１）は、予測対象の出力ｙを、入力計測値ｕの荷重付きの移動平均と、観測雑音ｗの和とで表現している。荷重付きの移動平均の荷重係数は｛ａ１，ａ２，…，ａｎ｝である。この荷重係数は予測対象のインパルス応答に相当する。観測雑音ｗは平均が「０」で、標準偏差が「σｗ」の正規分布、すなわちｗ〜Ｎ（０，σｗ２）と表されるものとする。「＾ｙｔ＋１」は「ｙｔ＋１」の予測値である。この「＾ｙｔ＋１」を求める式が予測モデルを表している。
【００３５】
カーネル型システム同定法は、入力の荷重係数｛ａ１，ａ２，…，ａｎ｝を多変量正規分布に従う確率変数としてベイズ推定するところに特徴がある。表記の簡単化のために、上式（１）の入力の荷重係数｛ａ１，ａ２，…，ａｎ｝をひとまとめにしたものを、列ベクトルａで表すことにする。列ベクトルａが従う前記多変量正規分布とは、より具体的には、平均が「０（ｎ×１）」、共分散行列が「Ｋ（ｎ×ｎ）」である。すなわち、ａ〜Ｎ（０，Ｋ）である。カーネル型システム同定法においては、共分散行列Ｋはカーネル行列と呼ばれており、予測対象に関する事前情報を表している。たとえば、一次安定スプラインカーネル（Tuned and Correlatedカーネル）に従うと、予測対象のインパルス応答が、時間に対して指数的に収束する時定数の事前情報からカーネル行列を定めることができる。具体的には、一次安定スプラインカーネルに従うと、カーネル行列Ｋのｉ行ｊ列要素は、以下の式（２）で与えられる。λは入力ｕに対する出力ｙのインパルス応答の大きさを表すパラメータであり、βはインパルス応答の収束の速さを表すパラメータである。これらの値は、予測対象の運転経験や運転の記録から、おおよその値を知ることは容易である。
【００３６】
【数２】
2020194535.tif 000004
【００３７】
また、ベクトル表現を導入すると、入力計測値ｕ及び出力計測値ｙは以下の式（３）で表される。なお、Ｎは計算に使用される観測点の数（ステップ数）であり、ｎはインパルス応答の次数（荷重係数ａの要素数）である。
【００３８】
【数３】
2020194535.tif 000005
【００３９】
これらのベクトルＵ，Ｙが分かっており、さらに、カーネル行列Ｋと、入力外乱の標準偏差σｖと、観測雑音の標準偏差σｗと、が事前情報として既知であるならば、ベイズ推定の意味において、荷重係数ａと出力計測値Ｙの予測値の事前分布は式（４）のように多変量正規分布として表される。式（４）において、記号Ｎの括弧の内側の第１項は期待値であり、第２項が共分散行列である。
【００４０】
【数４】
2020194535.tif 000006
【００４１】
このとき、前述の事前情報に加え、さらに出力計測値ｙを得たならば、ベイズ推論の意味において、荷重係数ａの最適推定値は事後分布の平均値として式（５）のように表される。式（５）の「＾ａ」を、予測モデルの定数（荷重係数ａ）に用いる。
【００４２】
【数５】
2020194535.tif 000007
【００４３】
式（４）において、左辺のａは式（１）の荷重係数｛ａ１，ａ２，…，ａｎ｝を列ベクトルとして表記したものである。式（４）の多変量正規分布において、共分散行列の（２，２）要素は出力計測値Ｙの予測値の共分散行列Ｖａｒ（Ｙ）を表している。本願の特徴はＶａｒ（Ｙ）の評価方法にあるので、その導出過程を説明しておく。式（４）のとおり、Ｖａｒ（Ｙ）＝ＵＫＵＴ＋σｗ２Ｉである。これは、観測雑音の時系列{ｗ１，ｗ２，…，ｗＮ}の列ベクトルＷ（Ｎ×１）で表すと、Ｖａｒ（Ｙ）＝Ｖａｒ（Ｕａ＋Ｗ）＝Ｅ（ＵａａＴＵＴ）＋Ｅ（ＵａＷＴ）＋Ｅ（ＷａＴＵＴ）＋Ｅ（ＷＷＴ）である。カーネル行列の定義からＥ（ａａＴ）＝Ｋ、観測雑音と入力は無相間だからＥ（ＵＷＴ）＝０、そして観測雑音の定義からＥ（ＷＷＴ）＝σｗ２Ｉである。したがって、Ｖａｒ（Ｙ）＝ＵＫＵＴ＋０＋０＋σｗ２Ｉとなる。
【００４４】
従来より、カーネル型システム同定法には、少ないデータ点数Ｎで予測モデルが得られるという利点がある。一般的な手法では、予測モデルの定数の数ｎが５から１０のときに、データ点数は２０００程度が必要となる場合がある。これに対し、カーネル型システム同定法では、データ点数Ｎを２００程度で同定することが可能となる。
【００４５】
しかしながら、従来のカーネル型システム同定法は、観測雑音を考慮するのみであり、入力に重畳する外乱があると、予測精度が低下する可能性があった。入力に重畳する外乱とは、入力を指令する信号の伝送に関わる雑音の他に、バルブなどの操作量の誤差などを含むものであり、以下では入力外乱と表記する。入力外乱に対処するため、本実施形態に係る予測装置３は、以下に説明するように、従来とは異なる同定処理Ｓ１及び予測処理Ｓ２を実行する。
【００４６】
（移動平均フィルタによる予測装置の処理）
本実施形態に係る予測装置３の同定部３１０は、移動平均フィルタ（ＭＡフィルタ）を用いるところは従来と同様である。しかし、本実施形態の移動平均フィルタでは、予測対象２の１ステップ先（時刻ｔ＋１）における出力「ｙｔ＋１」を、入力外乱ｖを考慮して以下の式（６Ａ）で表して、同定処理Ｓ１を行う。ここに、入力外乱ｖは、平均が「０」、標準偏差がσｖのガウス分布に従うものとする。これを、ｖ〜Ｎ（０，σｖ２）と表記する。なお、観測雑音ｗは従来の手法と同様である。「＾ｙｔ＋１」は「ｙｔ＋１」の予測値である。この「＾ｙｔ＋１」を求める式は、本実施形態における予測モデルを表している。
【００４７】
【数６Ａ】
2020194535.tif 000008
【００４８】
本実施形態では、入力外乱が予測値に及ぼす影響を考慮する点において、従来の手法と異なっている。式（６Ａ）の移動平均フィルタの形式のモデルについて、カーネル行列Ｋと、入力外乱の標準偏差σｖと、観測雑音の標準偏差σｗと、入力計測値ｕとが事前情報として既知であるならば、ベイズ推定の意味において、荷重係数ａ（第１係数）と出力計測値Ｙの予測値の事前分布は式（６Ｂ）のように多変量正規分布で表される。式（６Ｂ）において、記号Ｎの括弧の内側の第１項は期待値であり、第２項が共分散行列である。
【００４９】
式（６Ｂ）の多変量正規分布において、共分散行列の（２，２）要素は出力計測値Ｙの予測値の共分散行列Ｖａｒ（Ｙ）を表している。本願の特徴はＶａｒ（Ｙ）の評価方法にあるので、その導出過程を説明しておく。式（６Ａ）のとおり、Ｖａｒ（Ｙ）＝ＵＫＵＴ＋ＤＫ＋σｗ２Ｉである。これは、観測雑音の時系列{ｗ１，ｗ２，…，ｗＮ}を列ベクトルＷ（Ｎ×１）で表し、同様に、入力外乱の時系列{ｖ１，ｖ２，…，ｖＮ}を列ベクトルＶ（Ｎ×１）で表すと、Ｖａｒ（Ｙ）＝Ｖａｒ（（Ｕ＋Ｖ）ａ＋Ｗ）＝Ｅ（ＵａａＴＵＴ）＋Ｅ（ＶａａＴＶＴ）＋Ｅ（（Ｕ＋Ｖ）ａＷＴ）＋Ｅ（ＷａＴ（ＵＴ＋ＶＴ））＋Ｅ（ＷＷＴ）である。カーネル行列の定義からＥ（ａａＴ）＝Ｋ、入力は観測雑音および入力外乱と無相間だからＥ（（Ｕ＋Ｖ）ＷＴ）＝０、そして観測雑音の定義からＥ（ＷＷＴ）＝σｗ２Ｉである。したがって、Ｖａｒ（Ｙ）＝ＵＫＵＴ＋Ｅ（ＶＫＶＴ）＋０＋０＋σｗ２Ｉとなる。Ｅ（ＶＫＶＴ）は前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列であり、式（６Ｂ）では、ＤＫと記している。
【００５０】
【数６Ｂ】
2020194535.tif 000009
【００５１】
共分散行列を構成する行列ＤＫは本実施形態が特徴とする、第１係数ａに係るカーネル行列Ｋで重み付けした入力外乱の共分散行列である。行列ＤＫの構造は式（６Ｃ）で表され、その要素は式（６Ｄ）で計算される。式（６Ｄ）において、記号Ｔｒは行列のトレースまたは跡、すなわち行列の対角要素の総和を表している。また、ｅｐとｅｑはＮ次元線形空間のｐ番目とｑ番目の基底を表す行ベクトルであり、「ｐ＝ｑ」ならば「ｅｐｅｑＴ＝０」、「ｐ≠ｑ」ならば「ｅｐｅｑＴ＝１」である。
【００５２】
【数６Ｃ】
2020194535.tif 000010
【００５３】
【数６Ｄ】
2020194535.tif 000011
【００５４】
このとき、前述の事前情報に加え、さらに出力計測値ｙを得たならば、ベイズ推論の意味において、第１係数ａの最適推定値は式（６Ｅ）のように事後分布の平均値として表される。従来の技術による第１係数ａの最適値には、式（５）に示すとおり行列ＤＫがないので、入力外乱は考慮されておらず、したがって第１係数ａを正しく推定できなかった。両者の違いは、Ｙの事前分布の評価において、Ｖａｒ（Ｙ）として、第１係数ａに係るカーネル行列Ｋで重み付けした入力外乱の共分散行列ＤＫを考慮するか否かに起因する。本実施形態に係る技術では、標準偏差σｖと、第１係数ａに係るカーネル行列Ｋと、から、入力外乱が出力Ｙの事前分布に及ぼす影響を考慮するので、特に入力外乱があるときの第1係数ａの推定精度を向上することができる。
【００５５】
【数６Ｅ】
2020194535.tif 000012
【００５６】
上式（６Ｅ）の「＾ａ」の値を、予測モデルの定数（第１係数ａ）に用いる。
【００５７】
同定部３１０は、所定時間（たとえば、１０分）ごとに上述の同定処理Ｓ１を実行して予測モデルを更新する。これにより、同定部３１０は、予測対象２の特性の変化に応じた予測モデルを常に提供することができる。なお、所定時間は、予測対象２の特性の変動周期等に応じて任意に設定される。
【００５８】
また、予測部３１１は、同定部３１０により同定された定数（第１係数ａ）と、入力計測値及ｕ及び出力計測値ｙとからなる予測モデルに基づいて、予測対象２の１ステップ先（時刻ｔ＋１）の出力を予測する予測処理Ｓ２を実行する。
【００５９】
予測部３１１が予測した予測値「＾ｙｔ＋１」は、予測対象２の制御装置２１０に送信される。そうすると、制御装置２１０は、この予測値「＾ｙｔ＋１」に基づいて、予測対象２の出力が適正な値となるように制御、調整する。
【００６０】
（作用効果）
以上のように、本実施形態に係る予測システム１の予測装置３は、入力に対する第１係数ａの同定処理Ｓ１については、観測雑音の共分散行列σｗ２Ｉ（ｎ×ｎ）と、前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、を用いることを特徴とする。これにより、観測不可能な外乱による影響が生じた場合であっても、外乱の標準偏差についての事前情報を利用して、予測対象の出力の予測精度が低下することを抑制することができる。
【００６１】
＜第２の実施形態＞
以下、本発明の第２の実施形態に係る予測装置、及び当該予測装置を具備する予測システムについて、図１Ｂを参照しながら説明する。前述した第１の実施形態との差異は、第１の実施形態が移動平均フィルタを利用するのに対し、本実施形態では自己回帰移動平均フィルタを利用するところにある。これは、予測装置３の同定部３１０にて行う同定処理Ｓ１と予測部３１１にて行う予測処理Ｓ２の演算にのみ影響を及ぼすので、それのみを説明する。それ以外の機能と構成は第１の実施形態と同一である。
【００６２】
（自己回帰移動平均フィルタによる予測装置の処理）
本実施形態に係る予測装置３の同定部３１０は、前述の移動平均フィルタに代えて、自己回帰移動平均フィルタ（ＡＲＭＡフィルタ）を用いて同定処理Ｓ１を行う。自己回帰移動平均フィルタでは、予測対象２の１ステップ先（時刻ｔ＋１）における出力「ｙｔ＋１」を以下の式（６Ｆ）で表す。なお、観測雑音ｗは従来の手法と同様である。「＾ｙｔ＋１」は「ｙｔ＋１」の予測値である。この「＾ｙｔ＋１」を求める式は、本実施形態における予測モデルを表している。
【００６３】
【数６Ｆ】
2020194535.tif 000013
【００６４】
本実施形態では、過去の出力計測値ｙを予測に利用する点において、前実施形態の移動平均フィルタを用いる手法と異なっている。これに伴い、同定部３１０は、入力に対する係数（以下、「第１係数」とも記載する）のみならず、出力に対する係数（以下、「第２係数」とも記載する）も同定する。これを実現するために、入力に対する係数ａの共分散行列をカーネル行列Ｋａにより、出力に対する係数ｂの共分散行列をカーネル行列Ｋｂにより、個別に表記する。
【００６５】
カーネル行列Ｋａ及びＫｂのｉ行ｊ列要素は、それぞれ以下の式（８）及び（９）で与えられる。λａ、λｂ、βａ、βｂについては、値について事前情報が有り、おおよその値は分かっているものとする。なお、λａ及びλｂの値は、同一であってもよいし、異なっていてもよい。同様に、βａ及びβｂの値も、同一であってもよいし、異なっていてもよい。さらに、式（６Ｆ）では第１係数｛ａ１，ａ２，…｝と第２係数｛ｂ１，ｂ２，…｝は、ともに同じ長さｎであると記している。しかし、長さが異なっていても良い。例えば、第１係数の長さは０であっても良い。
【００６６】
【数８】
2020194535.tif 000014
【００６７】
【数９】
2020194535.tif 000015
【００６８】
また、次の式（１０）で表されるベクトル表現を導入すると、予測モデルの定数［ａＴ，ｂＴ］Ｔの事前分布は、式（１１）で表される。
【００６９】
式（１１）の多変量正規分布において、共分散行列の（２，２）要素は出力計測値Ｙの予測値の共分散行列Ｖａｒ（Ｙ）を表している。式（１１）では、Ｖａｒ（Ｙ）＝ＵａＫａＵａＴ＋ＵｂＫｂＵｂＴ＋ＤＫａ＋ＤＫｂ＋σｗ２Ｉと記している。これは以下のようにして導出される。観測雑音の時系列{ｗ１，ｗ２，…，ｗＮ}を列ベクトルＷ（Ｎ×１）で表し、同様に、入力外乱の時系列{ｖ１，ｖ２，…，ｖＮ}を列ベクトルＶ（Ｎ×１）で表すと、Ｖａｒ（Ｙ）＝Ｖａｒ（（Ｕａ＋Ｖ）ａ＋（Ｕｂ＋Ｗ）ｂ＋Ｗ）＝Ｅ（ＵａＫａＵａＴ）＋Ｅ（ＵｂＫｂＵｂＴ）＋Ｅ（ＶＫａＶＴ）＋Ｅ（ＷＫｂＷＴ）＋Ｅ（（Ｕａ＋Ｖ）ａｂＴ（Ｕｂ＋Ｗ）Ｔ）＋Ｅ（（Ｕｂ＋Ｗ）ｂａＴ（Ｕａ＋Ｖ）Ｔ）＋Ｅ（ＷＷＴ）である。したがって、Ｖａｒ（Ｙ）＝ＵａＫａＵａＴ＋ＵｂＫｂＵｂＴ＋ＤＫａ＋ＤＫｂ＋０＋０＋σｗ２Ｉとなる。Ｅ（ＶＫａＶＴ）は前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列であり、ＤＫａと記している。Ｅ（ＷＫｂＷＴ）は前記第２係数に係るカーネル行列で重み付けした観測雑音の共分散行列であり、ＤＫｂと記している。
【００７０】
【数１０】
2020194535.tif 000016
【００７１】
【数１１】
2020194535.tif 000017
【００７２】
【数１２】
2020194535.tif 000018
【００７３】
この事前分布に加えて、出力計測値ｙを得たならば、ベイズ推論の意味において、第１係数ａと第２係数ｂの最適推定値は式（１２）のように表される。式（１２）において、ＤＫａは前記第１係数に係るカーネル行列Ｋａで重み付けした入力外乱の共分散行列であり、第１の実施形態で述べたとおりである。ＤＫｂは前記第２係数に係るカーネル行列Ｋｂで重み付けした観測雑音ｗの共分散行列であり、本実施例で特徴とするものである。ＤＫａの計算は、第１の実施例におけるＤＫの計算において、ＫをＫａに代えて実施する。ＤＫｂの計算は、第１の実施例におけるＤＫの計算において、ＫをＫｂに代え、さらにσｖをσｗに代えて実施する。本願の技術の新規性として、例えば、非特許文献２を参照する。非特許文献２の式（３６）および式（３７）が本実施形態の式（１２）に対応しているが、非特許文献２の式にはＤＫａとＤＫｂに相当する項は記されていない。ＤＫａとＤＫｂは、それぞれ入力外乱と観測雑音が出力計測値Yの予測値の事前分布に及ぼす影響を表すものであるから、これらを欠く技術では、前記事前分布は実際に対して誤差を持つ。その結果、従来の技術では、予測モデルを正しく同定することができない。
【００７４】
なお、上式（１０）において、Ｎは計算に使用される観測点の数（ステップ数）であり、ｎは移動平均または自己回帰の要素数である。たとえば、「Ｎ＝１００」、「ｎ＝１０」等の値が任意に設定される。このとき、同定部３１０は、記録装置３０からＮステップ分の計測値からなる入力計測値ベクトルｕ及び出力計測値ベクトルｙを抽出し、これらから以下の式（１０）に表されるベクトルを得る処理Ｓ１１を実行しているとする。また、上式（１１）及び（１２）において、「σｗ」は観測雑音ｗの標準偏差であり、値について事前情報を有しており、おおよその値が分かっているとする。同様に「σｖ」は入力外乱ｖの標準偏差であり、値について事前情報を有しており、おおよその値が分かっているものとする。
【００７５】
上式（１２）の「＾ａ」及び「＾ｂ」の値を、予測モデルの定数（第１係数ａ及び第２係数ｂ）に用いる。
【００７６】
同定部３１０は、所定時間（たとえば、１０分）ごとに上述の同定処理Ｓ１を実行して予測モデルを更新する。これにより、同定部３１０は、予測対象２の特性の変化に応じた予測モデルを常に提供することができる。なお、所定時間は、予測対象２の特性の変動周期等に応じて任意に設定される。
【００７７】
なお、本実施形態では、第１係数ａ及び第２係数ｂの長さはともにｎであり等しい値である態様を例として説明したが、これに限られることはない。他の実施形態では、第１係数ａ及び第２係数ｂの長さは異なっていてもよい。
【００７８】
また、予測部３１１は、同定部３１０により同定された定数（第１係数ａ及び第２係数ｂ）と、入力計測値及ｕ及び出力計測値ｙとからなる予測モデルに基づいて、予測対象２の１ステップ先（時刻ｔ＋１）の出力を予測する予測処理Ｓ２を実行する。
【００７９】
予測部３１１が予測した予測値「＾ｙｔ＋１」は、予測対象２の制御装置２１０に送信される。そうすると、制御装置２１０は、この予測値「＾ｙｔ＋１」に基づいて、予測対象２の出力が適正な値となるように制御、調整する。
【００８０】
（実用例）
図２は、本発明の第２の実施形態に係る予測システムの実用例を示す図である。
以下、図２を参照しながら、本実施形態に係る予測システム１の実用方法の一例について説明する。なお、以下に説明する実用例は、第１の実施形態に適用してもよい。
【００８１】
図２に示すように、本実施形態に係る予測システム１は、電力の需給調整を適切に実施するために、発電所Ｇにより発電される電力の出力変動を予測する態様であるとする。電力の需給調整とは、電力系統Ｌ１の周波数を一定に維持することを示す。電力系統Ｌ１には、電力を消費する不図示の需要家（工場、一般家庭等）が接続され、これら需要家はそれぞれ必要に応じて電力を消費する。電力系統Ｌ１の一般的な性質として、電力の消費（需要）が発電（供給）を超過すると電力系統Ｌ１の周波数が低下し、逆に供給が需要を超過すると周波数が増加する。発電所Ｇは、周波数の変動が一定範囲（たとえば、基準値±０．２Ｈｚ）を逸脱しないように、時々刻々と変化する需要に応じて発電量を調整している。しかしながら、発電所Ｇの調整能力には限度があるため、需給調整を行う場合、将来の周波数変動を予測することは有効である。周波数が増加又は低下することが予め分かっていれば、先回りして発電量を低下又は増加させることができるからである。
【００８２】
このような知見をふまえ、本実施形態に係る予測システム１は、発電所Ｇにおける将来の周波数変動を予測するものである。
【００８３】
図２に示すように、予測システム１は、予測装置３と、計測器５０とを備えている。
予測装置３は、将来において発電所Ｇから電力系統Ｌ１へ供給される電力の周波数（出力）を予測する。
計測器５０は、発電所Ｇと電力系統Ｌ１との接続点に設置され、当該接続点において発電所Ｇから電力系統Ｌ１に供給される有効電力と、当該接続点における電力系統Ｌ１の周波数とを計測可能とする。
【００８４】
また、発電所には、発電した電力を電力系統Ｌ１に供給する複数の電源２１、２２、２３、…が設けられている。電源２１、２２、２３、…はいずれも同一の構成を有しているとする。このため、ここでは、複数の電源のうち電源２１の構成を例に説明する。電源２１は、制御装置２１０と、タービン装置２１１（たとえば、ガスタービン、蒸気タービン等）と、発電機２１２と、加減弁２１３とを有している。
【００８５】
制御装置２１０は、コンピュータにより構成され、タービン装置２１１及び発電機２１２の運転制御を行う。特に、制御装置２１０は、発電機２１２の回転速度（出力の周波数に対応）を常時モニタリングして、当該回転速度が一定に保たれるように、タービン装置２１１への燃料、又は蒸気の供給量を自動調整する。また、制御装置２１０は、予測装置３から受信した周波数の予測値「＾ｙｔ＋１」に基づいて、タービン装置２１１への燃料又は蒸気の供給量を自動調整する。
【００８６】
加減弁２１３は、制御装置２１０から受信した制御信号に基づいて動作することにより、タービン装置２１１に供給される燃料、又は蒸気の供給量を変更する。
【００８７】
電源２１は、電力系統Ｌ１に接続されている。電源２１と電力系統Ｌ１との接続点には、計測器５０が設置されている。電源２１と電力系統Ｌ１との接続点に設置された計測器５０は、当該電源２１から電力系統Ｌ１へと出力される有効電力の計測値、及び、周波数の計測値を取得する。計測器５０は、所定の通信網（インターネット回線等）を介して、有効電力及び周波数の計測値を予測装置３に送信する。予測装置３は発電所Ｇの内部に配置してもよい。さらに、予測装置３は制御装置２１０の内部に配置してもよい。
【００８８】
同様に、他の電源２２、２３、…と電力系統Ｌ１との接続点に設置された計測器５０は、当該電源２２、２３、…の各々から電力系統Ｌ１へと出力される有効電力の計測値、及び周波数の計測値を取得し、予測装置３に送信する。ここで、有効電力の計測値及び周波数の計測値は、それぞれ図１の入力計測値ｕ、出力計測値ｙに対応する。
【００８９】
また、制御装置２１０が予測装置３から受信した周波数の予測値「＾ｙｔ＋１」に基づいて、タービン装置２１１への燃料又は蒸気の供給量を自動調整する処理の詳細について説明する。たとえば、制御装置２１０は、周波数の予測値「＾ｙｔ＋１」を加減弁２１３の開度に利用する。周波数の基準値に対する偏差Δｆに応じて発電所Ｇが追加的に発生する出力ΔＰ１は、一般に調停率とよばれる定数δにより、以下の式（１３）で関連付けられる。ここで、ｆｎは電力系統Ｌ１の基準周波数、Ｐｎは電源２１，２２，２３，…の定格出力である。
【００９０】
【数１３】
2020194535.tif 000019
【００９１】
同様に、制御装置２１０は、回転速度の基準値からの偏差に基づいて、以下の式（１４）で算出したΔＰ２にしたがって、燃料又は蒸気の加減弁２１３の開度を調整している。ここで、Ｒｎは発電機２１２の定格回転数である。
【００９２】
【数１４】
2020194535.tif 000020
【００９３】
制御装置２１０は、上述の各式で求めたΔＰ１、及びΔＰ２の加重平均に基づいて、加減弁２１３の開度指令値を計算する。制御装置２１０は、このように求めた開度指令値を制御信号として加減弁２１３に送信することにより、タービン装置２１１への燃料又は蒸気の供給量を調整する。これにより、制御装置２１０は、電力系統Ｌ１の周波数（出力）を適切に維持できるように、電源２１，２２，２３，…の有効電力（入力）を精度よく調整することができる。すなわち、制御装置２１０は、発電所Ｇにおける調整能力を向上させることができる。
【００９４】
（作用効果）
以上のように、本実施形態に係る予測システム１の予測装置３は、プロセッサ３１と、プロセッサ３１に接続され、予測対象２の入力の計測値である入力計測値ｕ、及び出力の計測値である出力計測値ｙを記憶する記録装置３０と、を備える。プロセッサ３１は、過去に記憶された複数の入力計測値ｕ、及び複数の出力計測値ｙから、自己回帰移動平均フィルタを用いて入力に対する第１係数ａと、出力に対する第２係数ｂとを同定する同定処理Ｓ１と、入力計測値ｕ、出力計測値ｙ、第１係数ａ、及び第２係数ｂからなる予測モデルに基づいて、予測対象２の将来における出力を予測する予測処理Ｓ２と、を実行する。  第１実施例で述べた移動平均フィルタによる予測は、入力計測値ｕのみから出力を予測するので、観測不可能な外乱の影響によって出力計測値ｙが変動した場合であっても、これを予測に使うことが出来ない。したがって、外乱が大きい場合には予測精度が低下する。しかし、一般に自己回帰移動平均フィルタによる予測では、過去の出力計測値ｙも予測につかうので、観測不可能な外乱の影響によって出力計測値ｙが変動した場合であっても、その変動を予測に反映できる点は有利である。さらに、第２の実施形態で述べた自己回帰移動平均フィルタによる予測は、事前分布の評価において、第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列と、第２係数に係るカーネル行列で重み付けした観測雑音の共分散行列と、を用いるところに特徴がある。この特徴により、予測対象２の入力計測値ｕと出力計測値ｙとを関係付ける第１係数ａ、及び第２係数ｂを、観測不可能な入力外乱および出力外乱の影響を考慮に入れつつ決定することができる。これにより、第２の実施例による予測装置３は外乱による予測精度の低下が抑制される。
【００９５】
また、予測システム１は、予測装置３と、予測装置３と通信可能に接続され、当該予測装置３から受信した予測対象２の出力の予測値に基づいて、当該予測対象２の入力を調整する制御装置２１０と、を備える。
このようにすることで、予測システム１は、予測対象２の出力が適切な値となるように、入力を精度よく調整することができる。
【００９６】
また、予測対象２は、たとえば、電力系統Ｌ１に電力の供給を行う電源２１，２２，２３，…である。制御装置２１０は、電力系統Ｌ１の周波数（出力）の予測値に基づいて、電源２１，２２，２３，…が有するタービン装置２１１の加減弁２１３の開度を調整する。
このようにすることで、予測システム１の制御装置２１０は、電力系統Ｌ１の周波数（出力）を適切に維持できるように、電源２１，２２，２３，…の有効電力（入力）を精度よく調整することができる。すなわち、制御装置２１０は、発電所Ｇにおける調整能力を向上させることができる。
【００９７】
＜第３の実施形態＞
以下、本発明の第３の実施形態に係る予測装置、及び当該予測装置を具備する予測システムについて、図１Ｃを参照しながら説明する。前述した第２の実施形態との差異は、第２の実施形態が自己回帰移動平均フィルタ（ＡＲＭＡフィルタ）を利用するのに対し、本実施形態では無限インパルス応答フィルタ（ＩＩＲフィルタ）を利用するところにある。これは、予測装置３の同定部３１０にて行う同定処理Ｓ１と予測部３１１にて行う予測処理Ｓ２の演算にのみ影響を及ぼすので、それのみを説明する。それ以外の機能と構成は第１の実施形態と同一である。
【００９８】
（無限インパルス応答フィルタ（ＩＩＲフィルタ）による予測装置の処理）
本実施形態に係る予測装置３の同定部３１０は、前述の移動平均フィルタ、自己回帰移動平均フィルタ（ＡＲＭＡフィルタ）に代えて、無限インパルス応答フィルタを用いて同定処理Ｓ１を行う。無限インパルス応答フィルタでは、予測対象２の１ステップ先（時刻ｔ＋１）における出力「ｙｔ＋１」を以下の式（１５Ａ）で表す。「＾ｙｔ＋１」は「ｙｔ＋１」の予測値である。この「＾ｙｔ＋１」を求める式は、本実施形態における予測モデルを表している。
【００９９】
【数１５Ａ】
2020194535.tif 000021
【０１００】
無限インパルス応答フィルタと自己回帰移動平均フィルタとの違いは、内部状態量ｘの有無にある。自己回帰移動平均フィルタは、入力観測値ｕと出力の観測値ｙとに基づいて出力を予測した。これに対し、無限インパルス応答フィルタでは、出力の観測値に代えて、内部状態ｘを利用する。これにより、無限インパルス応答フィルタは、移動平均フィルタと同様に、入力観測値ｕのみから出力を予測することが可能となり実装が容易である。さらに、移動平均フィルタとの比較では、無限インパルス応答フィルタは、一般に、移動平均フィルタに比べてフィルタの次数（荷重係数の数ｎ）が小さくて済むので、数値演算の点で有利である。
【０１０１】
カーネル行列Ｋａ及びＫｂは、第２の実施形態と同様である。予測モデルの定数［ａＴ，ｂＴ］Ｔの事前分布は、多変量正規分布として式（１５Ｂ）のように表される。
【０１０２】
【数１５Ｂ】
2020194535.tif 000022
【０１０３】
式（１５Ｂ）の多変量正規分布において、共分散行列の（２，２）要素は出力計測値Ｙの予測値の共分散行列Ｖａｒ（Ｙ）を表している。式（１５Ａ）では、Ｖａｒ（Ｙ）＝ＵａＫａＵａＴ＋ＵｂＫｂＵｂＴ＋ＤＫａ＋σｗ２Ｉと記している。これは以下のように導出したものである。観測雑音の時系列{ｗ１，ｗ２，…，ｗＮ}を列ベクトルＷ（Ｎ×１）で表し、同様に、入力外乱の時系列{ｖ１，ｖ２，…，ｖＮ}を列ベクトルＶ（Ｎ×１）で表すと、Ｖａｒ（Ｙ）＝Ｖａｒ（（Ｕａ＋Ｖ）ａ＋Ｕｂｂ＋Ｗ）＝Ｅ（ＵａａａＴＵａＴ）＋Ｅ（ＵｂｂｂＴＵｂＴ）＋Ｅ（ＶａａＴＶＴ）＋Ｅ（（Ｕａ＋Ｖ）ａＷＴ）＋Ｅ（ＷａＴ（Ｕａ＋Ｖ）Ｔ）＋Ｅ（ＷＷＴ）である。したがって、Ｖａｒ（Ｙ）＝ＵａＫａＵａＴ＋ＵｂＫｂＵｂＴ＋ＤＫａ＋０＋０＋σｗ２Ｉとなる。Ｅ（ＶａａＴＶＴ）は前記第１係数に係るカーネル行列で重み付けした入力外乱の共分散行列であり、Ｅ（ＶＫａＶＴ）またはＤＫａと記している。この事前分布に加えて、出力計測値ｙを得たならば、ベイズ推論の意味において、第１係数ａと第２係数ｂの最適推定値は事後分布の平均値として式（１５Ｃ）のように表される。
【０１０４】
【数１５Ｃ】
2020194535.tif 000023
【０１０５】
（作用効果）
以上のように、本実施形態に係る予測システム１の予測装置３は、プロセッサ３１と、プロセッサ３１に接続され、予測対象２の入力の計測値である入力計測値ｕ、及び出力の計測値である出力計測値ｙを記憶する記録装置３０と、を備える。プロセッサ３１は、過去に記憶された複数の入力計測値ｕ、及び複数の出力計測値ｙから、無限インパルス応答フィルタの形式のモデルを用いて入力に対する第１係数ａと、出力に対する第２係数ｂとを同定する同定処理Ｓ１と、入力計測値ｕ、出力計測値ｙ、第１係数ａ、及び第２係数ｂからなる予測モデルに基づいて、予測対象２の将来における出力を予測する予測処理Ｓ２と、を実行する。無限インパルス応答フィルタを利用することにより、移動平均フィルタと同様に、入力観測値ｕのみから出力を予測することが可能であるから実装が容易である。さらに、移動平均フィルタとの比較では、一般に、無限インパルス応答フィルタは、移動平均フィルタに比べて、フィルタの次数（荷重係数の数ｎ）が小さくて済むから数値演算の点で有利である、などの利点が得られる。
【０１０６】
＜第４の実施形態＞
次に、本発明の第４の実施形態に係る予測システム１について、図３を参照しながら説明する。
第１及び第２の実施形態と共通の構成要素には同一の符号を付して詳細説明を省略する。
本実施形態に係る予測装置３は、予測対象２のｍステップ先（時刻ｔ＋ｍ）における出力の予測を行う点において、第１〜第３の実施形態と異なっている。
【０１０７】
（予測装置の処理）
図３は、本発明の第４の実施形態に係る予測システムの機能構成を示す図である。
図３に示すように、本実施形態では、予測対象２のｍステップ先（時刻ｔ＋ｍ）における出力の予測値「＾ｙｍ＋１」を以下の式（１６）で表す。この「＾ｙｔ＋１」を求める式は、本実施形態における予測モデルを表している。
【０１０８】
【数１６】
2020194535.tif 000024
【０１０９】
上式（１６）において、右辺第１項はｍステップ先の予測を表している。また、｛ｕｔ＋１，ｕｔ＋２，…，ｕｔ＋ｍ−１｝は将来の入力である。｛＾ｙｔ＋１，＾ｙｔ＋２，…，＾ｙｔ＋ｍ−１｝は無限インパルス応答フィルタによる将来の出力の予測値、または自己回帰移動平均フィルタによる将来の出力の予測値であってもよい。上式（１６）において、右辺第２項は自己回帰移動平均フィルタによる過去の出力の推定値を表している。また、｛ｕｔ―ｎ＋ｍ，ｕｔ―ｎ＋ｍ＋１，…，ｕｔ｝は過去の入力値、｛ｙｔ―ｎ＋ｍ，ｙｔ―ｎ＋ｍ＋１，…，ｙｔ｝は過去の出力の計測値である。
【０１１０】
予測装置３の予測部３１１は、予測処理Ｓ２において、上式（１６）で表される予測モデルに基づいて、予測対象２のｍステップ先（時刻ｔ＋ｍ）の出力を予測する。なお、ｍの値は、たとえば予測システム１の運用者が予測対象２に応じて任意の値を設定してもよい。
【０１１１】
予測部３１１が予測した予測値「＾ｙｔ＋ｍ」は、予測対象２の制御装置２１０に送信される。そうすると、制御装置２１０は、この予測値「＾ｙｔ＋ｍ」に基づいて、予測対象２の出力が適正な値となるように制御（調整）する。
【０１１２】
（作用効果）
以上のように、本実施形態に係る予測装置３のプロセッサ３１（予測部３１１）は、予測処理Ｓ２において、現在から所定時間後（ｍステップ先）の前記予測対象の出力を予測する。
このようにすることで、予測装置３は、第１の実施形態よりも長時間先における予測対象２の出力を予測することができる。
【０１１３】
また、予測システム１の制御装置２１０は、ｍステップ先の出力の予測値「＾ｙｔ＋ｍ」に基づいて、早い段階から予測対象２の出力を調整することができる。
【０１１４】
なお、第２の実施形態に係る実用例は、本実施形態にも適用可能である。この場合、予測システム１は、早い段階から電力系統Ｌ１の周波数（出力）を適切に維持できるように、電源２１，２２，２３，…の有効電力（入力）を調整することができる。
【０１１５】
＜第５の実施形態＞
次に、本発明の第５の実施形態に係る予測システム１について、図４を参照しながら説明する。
第１〜第４の実施形態と共通の構成要素には同一の符号を付して詳細説明を省略する。
本実施形態に係る予測装置３は、複数種類の入力、及び複数種類の出力を持つシステム等を予測対象２として扱う点において、第１〜第４の実施形態と異なっている。
【０１１６】
図４は、本発明の第５の実施形態に係る予測システムの機能構成を示す図である。
図４に示すように、本実施形態に係る予測装置３の記録装置３０は、予測対象２から複数種類の入力計測値（入力ベクトル）、及び複数種類の出力計測値（出力ベクトル）を記憶する。
【０１１７】
また、入力の次数をｎｕ、入力ベクトルを［ｕ１，…，ｕｎｕ］、出力の次数をｎｙ、出力ベクトルを［ｙ１，…，ｙｎｙ］、入力外乱を［ｖ１，…，ｖｎｕ］、その標準偏差を［σｖ１，…，σｖｎｕ］、観測雑音を［ｙ１，…，ｙｎｙ］、その共分散行列を［σｗ１，…，σｗｎｙ］と記すと、出力の予測式（予測モデル）は係数ベクトルＡｉ，Ｂｉ（ｉ＝１，２，…，ｎ）により、以下の式（１７）のように表される。
【０１１８】
【数１７】
2020194535.tif 000025
【０１１９】
また、本実施形態では、一次安定スプラインカーネルに従いカーネル行列ＫＡ、ＫＢを以下の式（１８）のように定める。
【０１２０】
【数１８】
2020194535.tif 000026
【０１２１】
さらに、カーネル行列ＫＡ及びＫＢのｉ行ｊ列要素は、それぞれ以下の式（１９）及び（２０）で与えられる。
【０１２２】
【数１９】
2020194535.tif 000027
【０１２３】
【数２０】
2020194535.tif 000028
【０１２４】
また、次の式（２１）で表されるベクトル表現を導入すると、予測モデルの定数［ＡＴ，ＢＴ］Ｔと出力Ｙの事前分布は、式（２２）で表される。この事前分布に加えて、出力計測値Ｙを得たならば、ベイズ推論の意味において、予測モデルの定数［ＡＴ，ＢＴ］Ｔの最適推定値は式（２３）のように表される。式（２２）において、ＤＫＡはカーネル行列ＫＡで重み付けした入力外乱の共分散行列であり、ＤＫｂはカーネル行列ＫＢで重み付けした入力外乱の共分散行列であり、本実施例で特徴とするものである。ＤＫＡとＤＫＢの計算は、第１の実施形態と同様である。なお、明細書中における「＾Ａ」、「＾Ｂ」との表記は、本実施形態において説明する図、式中において「Ａ」、「Ｂ」にハット記号「＾」が付された表記に対応する。
【０１２５】
【数２１】
2020194535.tif 000029
【０１２６】
【数２２】
2020194535.tif 000030
【０１２７】
【数２３】
2020194535.tif 000031
【０１２８】
上式（２３）において、「＾Ａ」はサイズｎ×ｎｕの行列であり、「＾Ｂ」はサイズｎ×ｎｙの行列である。これら「＾Ａ」及び「＾Ｂ」の値を、予測モデルの定数（第１係数Ａ及び第２係数Ｂ）に用いる。また、予測モデルのＡｉ，Ｂｊ（ｉ＝１，２，…，ｎ）は、それぞれ、Ａ、Ｂのｉ行ベクトルである。
【０１２９】
同定部３１０は、所定時間ごとに上述の同定処理Ｓ１を実行して予測モデルを更新する。これにより、同定部３１０は、予測対象２の特性の変化に応じた予測モデルを常に提供することができる。なお、所定時間は、予測対象２の特性の変動周期等に応じて任意に設定される。
【０１３０】
なお、本実施形態では、第１係数Ａ及び第２係数Ｂの行数はともにｎであり等しい値である例が示されているが、これに限られることはない。他の実施形態では、第１係数Ａ及び第２係数Ｂの行数は異なっていてもよい。
【０１３１】
また、予測部３１１は、同定部３１０により同定された定数（第１係数Ａ及び第２係数Ｂ）と、入力計測値及ｕ（入力ベクトル）及び出力計測値ｙ（出力ベクトル）とからなる予測モデルに基づいて、予測対象２の１ステップ先（時刻ｔ＋１）の複数種類（ｎｙ種類）の出力それぞれの値を予測する予測処理Ｓ２を実行する。
【０１３２】
予測部３１１が予測した予測値［＾ｙ；１，…，＾ｙ；ｎｙ］ｔ＋１は、予測対象２の制御装置２１０に送信される。そうすると、制御装置２１０は、この予測値［＾ｙ；１，…，＾ｙ；ｎｙ］ｔ＋１に基づいて、予測対象２の各出力値が適正な値となるように制御、調整する。
【０１３３】
なお、本実施形態では、予測部３１１が１ステップ先（時刻ｔ＋１）の出力を予測する態様を例として説明したが、これに限られることはない。予測部３１１は、第２の実施形態と同様に、ｍステップ先（時刻ｔ＋ｍ）の出力を予測するようにしてもよい。
【０１３４】
（作用効果）
以上のように、本実施形態に係る予測装置３の記録装置３０は、予測対象２の複数種類の入力計測値ｕ、及び複数種類の出力計測値ｙを記憶し、プロセッサ３１（同定部３１０）は、同定処理Ｓ１において、複数種類の入力それぞれに対する複数の第１係数Ａと、複数種類の出力それぞれに対する複数の第２係数Ｂと、を同定する。
このようにすることで、予測装置３は、予測対象２が複数の入力及び複数種類の出力を持つシステム等であった場合でも、この予測対象２の複数種類の出力それぞれを予測することができる。
【０１３５】
なお、第２の実施形態に係る実用例は、本実施形態にも適用可能である。この場合、予測装置３は、入力計測値として、計測器５０から受信した有効電力の計測値に加え、他の計測値（たとえば、発電機２１２の発電量等）を利用してもよい。さらに、他の発電所から受信した有効電力の計測値や、太陽光や風力などの自然エネルギーによる発電を支配する日射や風力などを利用しても良い。同様に、予測装置３は、出力計測値として、計測器５０から受信した周波数の計測値に加え、他の計測値（たとえば、発電機２１２の回転速度等）を利用してもよい。さらに、他の発電所から受信した周波数や発電機の回転速度を利用しても良い。
【０１３６】
＜ハードウェア構成＞
図５は、本発明の少なくとも一実施形態に係る予測装置、及び制御装置のハードウェア構成の一例を示す図である。
図５に示すように、コンピュータ９００は、プロセッサ９０１、メインメモリ９０２、ストレージ９０３、インタフェース９０４を備える。
【０１３７】
上述の予測装置３及び制御装置２１０は、それぞれコンピュータ９００に実装される。そして、上述した各処理部の動作は、プログラムの形式でストレージ９０３に記憶されている。プロセッサ９０１は、プログラムをストレージ９０３から読み出してメインメモリ９０２に展開し、当該プログラムに従って上記処理を実行する。また、プロセッサ９０１は、プログラムに従って、上述した各記憶部に対応する記憶領域をメインメモリ９０２に確保する。
【０１３８】
プログラムは、コンピュータ９００に発揮させる機能の一部を実現するためのものであってもよい。たとえば、プログラムは、ストレージ９０３に既に記憶されている他のプログラムとの組み合わせ、または他の装置に実装された他のプログラムとの組み合わせによって機能を発揮させるものであってもよい。なお、他の実施形態においては、コンピュータ９００は、上記構成に加えて、または上記構成に代えてＰＬＤ（Programmable Logic Device）などのカスタムＬＳＩ（Large Scale Integrated Circuit）を備えてもよい。ＰＬＤの例としては、ＰＡＬ(Programmable Array Logic)、ＧＡＬ(Generic Array Logic)、ＣＰＬＤ(Complex Programmable Logic Device)、ＦＰＧＡ（Field Programmable Gate Array）が挙げられる。この場合、プロセッサ９０１によって実現される機能の一部または全部が当該集積回路によって実現されてよい。
【０１３９】
ストレージ９０３の例としては、磁気ディスク、光磁気ディスク、光ディスク、半導体メモリ等が挙げられる。ストレージ９０３は、コンピュータ９００のバスに直接接続された内部メディアであってもよいし、インタフェース９０４または通信回線を介してコンピュータ９００に接続される外部メディア９１０であってもよい。また、このプログラムが通信回線によってコンピュータ９００に配信される場合、配信を受けたコンピュータ９００が当該プログラムをメインメモリ９０２に展開し、上記処理を実行してもよい。少なくとも１つの実施形態において、ストレージ９０３は、一時的でない有形の記憶媒体である。
【０１４０】
また、当該プログラムは、前述した機能の一部を実現するためのものであってもよい。さらに、当該プログラムは、前述した機能をストレージ９０３に既に記憶されている他のプログラムとの組み合わせで実現するもの、いわゆる差分ファイル（差分プログラム）であってもよい。
【０１４１】
以上、本発明の実施形態について詳細に説明したが、本発明の技術的思想を逸脱しない限り、これらに限定されることはなく、多少の設計変更等も可能である。
【符号の説明】
【０１４２】
１  予測システム
２  予測対象
２１，２２，２３  電源
２１０  制御装置
２１１  タービン装置
２１２  発電機
２１３  加減弁
３  予測装置
３０  記録装置
３１  プロセッサ
３１０  同定部
３１１  予測部
５０  計測器
【図１Ａ】
2020194535.tif 000032
【図１Ｂ】
2020194535.tif 000033
【図２】
2020194535.tif 000034
【図３】
2020194535.tif 000035
【図４】
2020194535.tif 000036
【図５】
2020194535.tif 000037
【図６】
2020194535.tif 000038