JP2022173863A - モジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラム - Google Patents
モジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラム Download PDFInfo
- Publication number
- JP2022173863A JP2022173863A JP2021079844A JP2021079844A JP2022173863A JP 2022173863 A JP2022173863 A JP 2022173863A JP 2021079844 A JP2021079844 A JP 2021079844A JP 2021079844 A JP2021079844 A JP 2021079844A JP 2022173863 A JP2022173863 A JP 2022173863A
- Authority
- JP
- Japan
- Prior art keywords
- prediction
- pairwise
- output
- model
- input
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 183
- 230000008569 process Effects 0.000 claims abstract description 59
- 230000002194 synthesizing effect Effects 0.000 claims abstract description 14
- 238000011156 evaluation Methods 0.000 claims description 34
- 238000001308 synthesis method Methods 0.000 claims description 31
- 238000012937 correction Methods 0.000 claims description 18
- 238000013480 data collection Methods 0.000 claims description 12
- 238000013500 data storage Methods 0.000 claims description 11
- 230000004044 response Effects 0.000 claims description 4
- 238000012546 transfer Methods 0.000 claims description 4
- 238000013277 forecasting method Methods 0.000 claims description 2
- 230000009466 transformation Effects 0.000 claims description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 37
- 238000013459 approach Methods 0.000 description 28
- 238000013075 data extraction Methods 0.000 description 21
- 238000010586 diagram Methods 0.000 description 21
- 238000012544 monitoring process Methods 0.000 description 19
- 230000000694 effects Effects 0.000 description 17
- 230000008901 benefit Effects 0.000 description 15
- 238000012935 Averaging Methods 0.000 description 10
- 230000006870 function Effects 0.000 description 10
- 239000010865 sewage Substances 0.000 description 9
- 239000000284 extract Substances 0.000 description 8
- 230000010354 integration Effects 0.000 description 8
- 238000012545 processing Methods 0.000 description 8
- 235000020681 well water Nutrition 0.000 description 8
- 239000002349 well water Substances 0.000 description 8
- FFBHFFJDDLITSX-UHFFFAOYSA-N benzyl N-[2-hydroxy-4-(3-oxomorpholin-4-yl)phenyl]carbamate Chemical compound OC1=C(NC(=O)OCC2=CC=CC=C2)C=CC(=C1)N1CCOCC1=O FFBHFFJDDLITSX-UHFFFAOYSA-N 0.000 description 7
- 230000015572 biosynthetic process Effects 0.000 description 7
- 238000004422 calculation algorithm Methods 0.000 description 7
- 230000015556 catabolic process Effects 0.000 description 7
- 238000012417 linear regression Methods 0.000 description 7
- 238000003786 synthesis reaction Methods 0.000 description 7
- 238000004458 analytical method Methods 0.000 description 6
- 230000036961 partial effect Effects 0.000 description 6
- 239000000126 substance Substances 0.000 description 6
- 230000002159 abnormal effect Effects 0.000 description 5
- 230000002411 adverse Effects 0.000 description 5
- 235000000332 black box Nutrition 0.000 description 5
- 230000008859 change Effects 0.000 description 5
- 238000005259 measurement Methods 0.000 description 5
- 239000013598 vector Substances 0.000 description 5
- 230000009467 reduction Effects 0.000 description 4
- 229910000831 Steel Inorganic materials 0.000 description 3
- 230000009471 action Effects 0.000 description 3
- 238000010276 construction Methods 0.000 description 3
- 238000010219 correlation analysis Methods 0.000 description 3
- 238000007405 data analysis Methods 0.000 description 3
- 238000012217 deletion Methods 0.000 description 3
- 230000037430 deletion Effects 0.000 description 3
- 238000009826 distribution Methods 0.000 description 3
- 230000005611 electricity Effects 0.000 description 3
- 230000005284 excitation Effects 0.000 description 3
- 238000010801 machine learning Methods 0.000 description 3
- 238000004519 manufacturing process Methods 0.000 description 3
- 238000005086 pumping Methods 0.000 description 3
- 239000010959 steel Substances 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 230000001364 causal effect Effects 0.000 description 2
- 239000002131 composite material Substances 0.000 description 2
- 230000007423 decrease Effects 0.000 description 2
- 238000013135 deep learning Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 230000007774 longterm Effects 0.000 description 2
- 238000012423 maintenance Methods 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000008450 motivation Effects 0.000 description 2
- 230000002250 progressing effect Effects 0.000 description 2
- 239000004065 semiconductor Substances 0.000 description 2
- 238000003860 storage Methods 0.000 description 2
- 238000012706 support-vector machine Methods 0.000 description 2
- 241000408659 Darpa Species 0.000 description 1
- 238000004378 air conditioning Methods 0.000 description 1
- 238000007630 basic procedure Methods 0.000 description 1
- 230000006399 behavior Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 230000003111 delayed effect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000001687 destabilization Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000006866 deterioration Effects 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 235000013305 food Nutrition 0.000 description 1
- 230000001939 inductive effect Effects 0.000 description 1
- 230000000670 limiting effect Effects 0.000 description 1
- 230000007257 malfunction Effects 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000013488 ordinary least square regression Methods 0.000 description 1
- 230000002085 persistent effect Effects 0.000 description 1
- 238000013439 planning Methods 0.000 description 1
- 238000010248 power generation Methods 0.000 description 1
- 239000000047 product Substances 0.000 description 1
- 238000000746 purification Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000007637 random forest analysis Methods 0.000 description 1
- 230000002829 reductive effect Effects 0.000 description 1
- 238000000611 regression analysis Methods 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 230000007480 spreading Effects 0.000 description 1
- 238000003892 spreading Methods 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
- 230000002459 sustained effect Effects 0.000 description 1
- 208000024891 symptom Diseases 0.000 description 1
- 238000012731 temporal analysis Methods 0.000 description 1
- 238000000700 time series analysis Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B23/00—Testing or monitoring of control systems or parts thereof
- G05B23/02—Electric testing or monitoring
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Automation & Control Theory (AREA)
- Testing And Monitoring For Control Systems (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
【課題】 予測結果を合理的に説明するとともに、容易に予測結果を調整可能とするモジュラー型時系列データ予測装置を提供する。【解決手段】 実施形態による装置は、複数のプロセス変数の時系列データを所定の周期で収集保存し、複数の前記プロセス変数の中から予測対象となる出力変数を選択し、複数のプロセス変数の中から入力変数の候補を選択する機能2と、時系列データから抽出された出力変数と複数の入力変数との同定用データを用いて、1入力1出力のペアワイズ予測モデルを定義する機能4と、複数のペアワイズ予測モデルの予測値の合成法を定義する機能5と、複数のペアワイズ予測モデルに複数の入力変数の予測用データを入力して、複数の入力変数のそれぞれに対応する予測値を演算する機能6と、定義された合成法により複数の予測値を合成することにより出力変数の予測値を演算する機能7と、を有する。【選択図】図2
Description
本発明の実施形態は、モジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラムに関する。
上下水道システム、雨水排水システム、電力システム、交通システムなどのインフラシステム、あるいは、鉄鋼プロセス、石油化学プロセス、半導体製造プロセス、などのプロセス系の産業プラントなどでは、通常、複数のプロセス状態を測定する複数のオンラインセンサが設置されている。プロセス監視制御システム(SCADA: Supervisory Control And Data Acquisition)と呼ばれるシステムは、上記のインフラシステムや産業プラントに設置されたセンサ群の計測により得られるプロセスデータ(流量、温度、水質、操作量など)を取得し、時系列データとしてサーバ上に保持している。プロセス監視制御システムは、通常、これらの時系列データをトレンドグラフとして監視員に提供することが多い。
監視員は、上記のような直接的な時系列データのトレンド監視に加えて、予測情報を利用するプラント監視を行うことがある。予測情報を利用するプラント監視システムによれば、監視中のリアルタイムの実測値(計測値)だけでなく、監視対象の将来の値を予測して、予測結果を運転管理者に提供することにより、適切な運転計画の実現を支援することができる。予測情報を利用するプラント監視システムは、プラントの異常兆候検出および異常診断と並んで、アドバンスト監視の典型的な手段である。
予測情報を利用するプラント監視が用いられるインフラシステムの例として、例えば、浄水システムにおける水需要の予測、下水処理場への汚水流入予測、雨天時の下水処理場および雨水排水ポンプ場への雨水流入予測、雨天時の河川水位やダム水位の予測、あるいは、気象レーダなどによる降雨の短期予測(ナウキャスト)、交通システムにおける交通量予測や渋滞予測、電力システムにおける電力需要の予測、あるいは、ビルなどの施設内の食堂や店舗などの顧客数予測、などのシステムがあり、インフラ領域の様々な分野で広く使われている。
予測情報を利用するプラント監視システムは、予測情報そのものを付加価値として提供するだけでなく、プラントの運転計画やプラントの制御に予測情報を利用することで、より効率の良い運用やより安全性の高い運用を実現するためにも用いられ得る。例えば、水需要や電力需要等の需要予測は、水運用計画や発電および蓄電計画(EMS:エネルギーマネージメントシステム)などの最適化に利用することで、効率の良い運用を実現するために用いることができる。また、例えば降雨予測や流入予測は、雨水排水ポンプの運転制御のための入力情報として予測値を用いることで、浸水の回避や抑制などのリスク低減に貢献する。
インフラシステムだけでなくプロセス系のプラントである鉄鋼プロセスや石油化学プロセスなどにおいても、製品品質や歩留まりの予測などに時系列データ解析を用いることも多く、その情報は生産効率(歩留まり)向上のために利用されることも多い。
時系列データの予測は、インフラシステムやプロセス系の産業プラントで極めて重要な役割を果たし、プラントの監視データを利用したデータドリブンの(ブラックボックス的アプローチの)予測方法として、従来から様々な方法が適用されている。
従来からプラント監視データなどの時系列データを用いた予測技術は広く用いられているが、近年の統計的機械学習分野を中心とするAI関連の理論、技術、手法の急速な進歩に伴い、時系列データ解析に対するAI手法も同時並行的に進展し、分野および領域を問わず、様々なアドバンストなAI手法を時系列データの予測にも適用する動きが加速している。特に、NNの発展形である深層学習ネットワーク(DNN:ディープラーニングネットワーク)などは、昨今注目を集めている代表的な方法であり、DNNはAIの代名詞としても使われる場合もある。
AIの代表的な手法であるDNNの他にも、重回帰分析で問題になる多重共線性による不安定化(悪条件問題)を回避する方法として(先に述べたPLSとは異なる観点で)導入された正則化に基づく各種の方法、特に、古くからあるL2正則化(リッジ回帰)に加え、多数の説明変数候補の中から必要な説明変数を自動選択することのできるL1正則化手法であるLassoなどの方法、Lassoとは異なる視点でベイズ推論的な手法で説明変数の自動選択を行うRVM(適合ベクトルマシン)に基づく回帰手法(予測手法)であるRVR(適合ベクトル回帰)などの方法、最適化を用いて識別問題(分類問題)をロバスト化したSVM(サポートベクトルマシン)を回帰問題(予測問題)に転用したSVR(サポートベクトル回帰)などの方法、従来の重回帰などのような線形回帰を非線形に拡張し、さらに確率分布を考えてノンパラメトリック回帰問題に帰着させ、点予測ではなく予測分布を出力するガウス過程回帰などの方法、あるいは、様々な予測サブシステムを組み合わせて予測性能を向上させるバギングやブースティングと呼ばれる各種の方法(ランダムフォレスト、アダブースト、XGブーストなど)など、多様なアドバンストな時系列データの予測手法も用いられる様になってきている。
このように、アドバンストなAIに基づく予測手法が次々と開発され、従来から用いられる方法と比較して各段に予測精度が向上するという報告も数多くなされる様になり、AI関連技術に関する期待が高まっている。
一方で、AI関連技術が実際の現実的なシステム(インフラシステムやプロセス系のプラント)において、必ずしも容易に展開できない理由の一つとして、統計的およびAI的な解析手法が持つ本質的な性質であるブラックボックス的なアプローチの欠点として指摘されている、「結果に対する説明性の欠如」の問題がある。そして、この「結果に対する説明性の欠如」の問題に対する意識が急速に高まっており、この問題を解決すべく、最近では「説明可能AI(XAI:eXplainable AI)の概念が米国DARPAにより提唱され、注目が集まっている。
インフラ系のシステムでは、物理的、化学的な知見に立脚したホワイトボックス的なモデルを用いる事も多いが、これは、「Accountability(説明責任)」が求められることに関連していると考えられる。実際、ホワイトボックス的なモデルによる予測精度がブラックボックス的なモデルによる予測精度よりも必ずしも高いわけではなく、また、ホワイトボックス的なモデルが必ずしも、実際の物理現象を正確かつ忠実に再現しているとも限らない。それにも関わらず、ホワイトボックスモデルが時として好まれ、利用される理由は、ホワイトボックスモデルに含まれる「パラメータ」には物理的および化学的な意味があるため「合理的な説明がしやすく、他者に納得してもらいやすい」ことに加え、このような物理的および化学的な意味を持つことにより、精度が十分でない場合に、どのように調整すればよいかという方針や指針をたてやすい、すなわち、調整しやすい、という利点があるためであると考えられる。実際、物理的および化学的に意味を持つパラメータは、パラメータの値自体に物理的および化学的意味があるため、理解がしやすく、そのパラメータの値の取りうる範囲を想定することや、値の増減が結果に与える影響を物理的および化学的知見に基づいて考えることができる。
上記のように、「合理的な説明のしやすさと調整のしやすさ」がホワイトボックスモデルを用いることの大きな動機付けになっている可能性がある。もしAI的な手法によるブラックボックスモデルにおいても「合理的な説明のしやすさと調整のしやすさ」という要素を付加することができれば、その応用範囲は大きく広がると考えられる。
本発明の実施形態は、上記事情を鑑みて成されたものであって、予測結果を合理的に説明するとともに、容易に予測結果を調整可能とするモジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラムを提供することを目的とする。
実施形態によるモジュラー型時系列データ予測装置は、複数のプロセス変数を所定の周期で計測する複数のプロセスセンサを有するシステム又はプロセスに適用される装置であって、複数の前記プロセス変数の時系列データを所定の周期で収集し保存するとともに、複数の前記プロセス変数の中から予測対象となる少なくとも一つの出力変数を選択する出力変数データ選択部と、複数の前記プロセス変数の中から複数の入力変数の候補を選択する入力変数データ選択部と、を含むデータ収集保存部と、前記時系列データから抽出された前記出力変数と複数の前記入力変数との同定用データを用いて、1入力1出力のペア毎にペアワイズ予測モデルのパラメータを同定して複数の前記ペアワイズ予測モデルを定義するペアワイズ予測モデル同定部と、複数の前記ペアワイズ予測モデルから出力されるペアワイズ予測値の合成法を定義する予測モデル合成法定義部と、時間の進行方向に所定の周期又はリアルタイムで前記時系列データから抽出された複数の前記入力変数の予測用データを、複数の前記ペアワイズ予測モデルに入力して、複数の前記入力変数のそれぞれに対応する前記ペアワイズ予測値を演算するペアワイズ出力変数予測部と、前記合成法により複数の前記ペアワイズ予測値を合成して前記出力変数の予測値を演算する合成出力変数予測部と、を有する。
本実施形態では、調整可能な時系列デー予測を行うための基本的手段として、モジュラー型の時系列データ予測装置について説明する。本実施形態のモジュラー型の時系列予測装置は、予測モデル(システム)を部分モデル(サブシステム)に分解可能にして、各部分モデルを容易に着脱可能にするものである。
モジュラー型の時系列データ予測装置によれば、例えば、複数の監視変数中のいずれか一つの変数のデータが欠測していたり、複数の監視変数にアウトライア(異常データ)が含まれていたり、監視変数がノイズ等を多量に含む不安定な計測データである場合、その影響が全体に波及することを回避できる。また、モジュラー型の時系列データ予測装置によれば、例えば、説明変数が必要であるか否(不要)かを適宜判断して、部分モデルを容易に脱着したり、予測モデルを部分的に調整したりすることができる。
このような利点は、モジュラー型の構成によりモデルが「分解可能」であるという事によって得られる利点である。予測モデルにおいてモジュラー型の構成をとる利点は、単に分解可能ということによって得られる利点だけでなく、これに加えて、このモジュラー型の構成をとることで、物理法則に矛盾せずに合理的で納得できるモデルを構築することが可能になり、さらに、このことにより調整を容易にすることができる点である。
予測モデルにおいてモジュラー型の構成をとる利点について以下に説明する。
一般に、ブラックボックス的な手法が説明可能(Explainable)でない理由は、入力および出力という外部情報だけしかわからず、内部構造が見えない不可視である事が原因であると考えられている。すなわち、ブラックボックス的な手法は、単に内部構造が見えないだけでなく、そもそも、外部の入出力関係からだけでは、内部構造を一意に決めることができず、同じ入出力関係を与える複数(無数)の内部構造がありうる。このため、内部構造の解釈が難しくなり、場合によっては物理法則に矛盾する不合理なモデルになってしまうことがある。これは、モデルの一意性の問題として、制御理論分野の中に含まれるであるシステム同定分野等で知られた事項である。このようなモデルの一意性の概念は、例えば、「ある入力データと出力データのペアに対して、その入出力関係を表すモデルを唯一に決定することができるか」という問題であり、モデルの可同定問題と呼ばれている。
一般に、ブラックボックス的な手法が説明可能(Explainable)でない理由は、入力および出力という外部情報だけしかわからず、内部構造が見えない不可視である事が原因であると考えられている。すなわち、ブラックボックス的な手法は、単に内部構造が見えないだけでなく、そもそも、外部の入出力関係からだけでは、内部構造を一意に決めることができず、同じ入出力関係を与える複数(無数)の内部構造がありうる。このため、内部構造の解釈が難しくなり、場合によっては物理法則に矛盾する不合理なモデルになってしまうことがある。これは、モデルの一意性の問題として、制御理論分野の中に含まれるであるシステム同定分野等で知られた事項である。このようなモデルの一意性の概念は、例えば、「ある入力データと出力データのペアに対して、その入出力関係を表すモデルを唯一に決定することができるか」という問題であり、モデルの可同定問題と呼ばれている。
可同定問題は、モデルの構造(モデルを表現する数式の形)を予め限定した場合、例えば、モデルを線形伝達関数モデルや時系列モデルに限定した場合、そのモデルに含まれるパラメータ(例えば伝達関数に含まれる係数)を一意に決めることができるかという問題となる。この場合、可同定問題の中でも特にパラメータ可同定性の問題と呼ばれる。このパラメータ可同定性の問題は、現実のデータを扱う問題でも頻繁に現れる問題であり、以下の3つのケースは特に遭遇することが多い問題である。
一つ目は、先に述べた多重共線性の問題であり、これは、入力変数(説明変数)間に強い相関を持つ場合の問題である。これは、例えば、雨水ポンプ場への雨水流入量を予測する場合に、入力変数としてレーダ雨量情報を用いて、近接する異なる2つのメッシュの降雨量を入力変数とするような場合などに典型的に現れる問題であり、現実のデータでは、説明変数間に強い相関を持つことは極めて多い。
二つ目は、制御理論の分野でよく知られている、持続的励振条件(Persistent Exciting Condition、PE条件)と呼ばれる条件が成立しない場合の問題であり、これは、モデル化の対象プロセスに対して、十分な周波数成分を含まない入力データが適用されている場合に生じる問題である。極端ではあるが、現実にも認められる例として、入力データが一定値である場合、通常出力データも一定値になるが、このような一定値の入出力関係を表す伝達関数モデルは無数に存在するため、係数を一意に決めることができない。例えば、y(t)=b1×u(t-1)+b2×u(t-2)と表される簡単な伝達関数モデル(時系列モデル)を考える場合、u(t-1)=u(t-2)=k(一定値)であるとb1とb2を一意に決めることはできない。この概念を一般化したものが持続的励振条件というパラメータ可同定条件であるが、このような条件が成立しない場合も多い。
三つ目は、モデルの構造の複雑さと比較して、データの得られる量が少なくデータの質が良くない場合である。典型的な例として、ある対象が非線形的な挙動を示すことが物理的な観点から予めわかっているが、得られるデータが少なかったり、ある特定の運用条件の付近のデータしか得られなかったりすることにより、非線形の係数を同定できない様な場合である。例えば、y(t)=c1×u(t)+c2×u(t)^2と書かれることがわかっているが、データとして得られるu(t)の値が小さい値の場合、u(t)^2<<u(t)となるため、c2を同定することが困難になる様な場合である。
上記三つの例は、パラメータの可同定性が欠落する典型的なケースであるが、現実の問題では、このような状況に遭遇することが極めて多い。モデル化の対象が、比較的小さな機械系システムおよび電気系システムの場合は、可同定性を向上させるための実験を行い、可同定性を確保した上で、モデルを構築することが常道である。しかし、上下水道、雨水排水、電力システム、などの巨大な社会インフラシステムの場合、モデル化の対象となるデータが自然現象から得られるデータ(降雨、日射量、風速など)であったり、対象プラントが稼働している状態のデータであったりすることが多く、パラメータの可同定性を向上させるための実験などを行うことができず、モデル化に利用できるデータは与えられたものであることが多いため、可同定性の欠落を回避することが本質的に難しい場合が多い。この結果、社会インフラシステムの場合には、可同定性が欠落した状態でモデルを構築せざるを得ない。
以下、可同定性の欠落による問題の一例として「多重共線性」について、従来の技術と本実施形態のモジュラー型時系列データ予測装置とを対比する。
「多重共線性」とは、複数の説明変数(入力変数)間に相関関係(一次従属関係)があることにより、回帰モデルのパラメータを一意に決めることができなくなるという問題である。これを最も簡単な例を用いて説明する。
「多重共線性」とは、複数の説明変数(入力変数)間に相関関係(一次従属関係)があることにより、回帰モデルのパラメータを一意に決めることができなくなるという問題である。これを最も簡単な例を用いて説明する。
一つの出力yが二つの入力u1、u2の重回帰で表される以下の重回帰モデルを考える。
y=a1×u1+a2×u2 (1)
(1)式は最も簡単な重回帰モデルであり、a1とa2とは回帰係数と呼ばれるパラメータである。ここで、物理的な具体的イメージを持つために、例えば、出力yを雨水ポンプ場への雨水流入量、入力u1および入力u2を各々異なる箇所の降雨量であることを想定して考えてみる。この場合、降雨量が増加すれば流入量が増加することは、物理的に考えて自明であるから、a1とa2とは正の値を持たなければならないことは容易にわかる。ここで、降雨量u1と降雨量u2とに従属関係のある多重共線性の問題がある場合を想定する。
y=a1×u1+a2×u2 (1)
(1)式は最も簡単な重回帰モデルであり、a1とa2とは回帰係数と呼ばれるパラメータである。ここで、物理的な具体的イメージを持つために、例えば、出力yを雨水ポンプ場への雨水流入量、入力u1および入力u2を各々異なる箇所の降雨量であることを想定して考えてみる。この場合、降雨量が増加すれば流入量が増加することは、物理的に考えて自明であるから、a1とa2とは正の値を持たなければならないことは容易にわかる。ここで、降雨量u1と降雨量u2とに従属関係のある多重共線性の問題がある場合を想定する。
この多重共線性問題は、実際の現象においても頻繁に見られることであり、この場合、例えば、入力u1と入力u2とをレーダ雨量データの隣接メッシュの降雨量であることを想定してみれば、この問題が特殊な問題ではなく、むしろ自然な問題であることが容易にわかる。例えば、国土交通省が配信するXRAINと呼ばれる気象レーダでは、1メッシュのサイズが250m×250mであるため、隣接メッシュの降雨量とは互いに数百メートル離れた箇所の降雨量を意味しており、その降雨量がほぼ同じであることは物理的に考えて極めて自然な事である。2つの降雨量データがレーダ雨量の隣接メッシュの降雨量ではなくて、2か所の地上雨量データであったとしても、通常雨水排水区の大きさが数キロ四方程度以内であることを考えると比較的近い値を持つことは容易に想像できる。従って、このような場合u1≒u2となっていることは、自然な事であり、これは多重共線性問題そのものである。
そこで、ここでは簡単のため、u1=u2という完全な多重共線性の状態にあった場合を想定する。そして、例えば、20mm/hの強度の雨量に対して、10m3/sの雨水流入があるとする。今、u1=u2であるので、u:=u1=u2と定義すると、(1)式は、以下の(2)式の様に書き換えられる。
y=(a1+a2)×u (2)
上記(2)式は単回帰の形であり、u:=u1=u2であるから、a1+a2=0.5となる。しかし、a1+a2=0.5という制約条件を満足してさえいれば、a1とa2とがいかなる値であっても、(2)式は全く同じ出力結果を出すため、a1とa2の各値はa1+a2=0.5を満たす範囲で任意の値を持つことができる。そのため、例えば、a1=3.5、a2=-3でも、a1=100、a2=-99.5となっても問題がなく、全く同じ入出力関係を与えるため、各々のパラメータ値の相違を入出力データから識別することができなくなる。
上記(2)式は単回帰の形であり、u:=u1=u2であるから、a1+a2=0.5となる。しかし、a1+a2=0.5という制約条件を満足してさえいれば、a1とa2とがいかなる値であっても、(2)式は全く同じ出力結果を出すため、a1とa2の各値はa1+a2=0.5を満たす範囲で任意の値を持つことができる。そのため、例えば、a1=3.5、a2=-3でも、a1=100、a2=-99.5となっても問題がなく、全く同じ入出力関係を与えるため、各々のパラメータ値の相違を入出力データから識別することができなくなる。
一方、a1とa2とは、各雨量に対する係数であるから、本来、すくなくとも正の値を持つべきであり、2地点の降雨量が同じであるなら、その値はほぼ等しくなるべきであることも物理的に考えれば妥当である。しかし、このように、入力変数間の相関関係(一次従属関係)により、パラメータ値を一意に決めることができず、本来のパラメータ値を正しく同定できなくなる事が多重共線性の問題の本質である。
このような多重共線性などのパラメータ可同定性の欠如により、入出力関係から一意にパラメータを決めることができなくなる問題は、システム同定分野や統計分野(統計的機械学習分野≒AI分野)では、数学的な「悪条件問題」として扱われ、この悪条件を回避するための各種のアルゴリズムが開発されている。先に述べた正則化やPLSなどの方法は、このような悪条件回避のための代表的な技術的手段(テクニック)であり、アドバンストなAI手法の多くには、正則化や(PLSで用いられる)次元圧縮などの手法(テクニック)が採用されている。
しかし、例えば、正則化やPLSなどのテクニックを駆使したアドバンストなアルゴリズムを(1)式の問題に適用したとしても、上記の問題に関して、a1>0、a2>0、かつa1≒a2となるような値を推定するとは限らず、むしろ、そのような値にならない可能性の方が高い。これは、正則化や次元圧縮などは、多重共線性(やこれを一般化した可同定性の欠如)の問題を、数学的な悪条件の問題として捉え、これを回避するためのある種の“テクニック”であり、悪条件を回避してパラメータを推定(同定)できるようにしているだけであるからである。すなわち、多重共線性を持つ重回帰問題に対して古典的な最小2乗法を適用すると、「答えを求めることができない=パラメータ値を求めることができない」という問題になるため、悪条件を回避するというテクニック(正則化や次元圧縮)をアルゴリズムに導入することで、「答えを求めることができる=パラメータ値を同定できる」問題にしているだけであり、同定した値が物理的に合理的に納得できるものであるか否かを直接考慮しているわけではないからである。
実際、例えば、正則化のテクニックでは、悪条件を回避するために、パラメータを同定する際に、予測誤差を評価するだけでなく、「パラメータの値自身が大きくなりすぎない」というパラメータに関するノルム(大きさを図る指標)を小さくする、という評価指標を加えることで悪条件を回避しており、これはパラメータの物理的な意味づけや解釈を行うこととは直接関係していない。
一方、このような問題に対して、合理的に納得のできるパラメータ値を得る方法について以下に説明する。はじめに、(1)式のモデルではなく、(3)式と(4)式との単回帰モデルを考える。
y=a1´×u1 (3)
y=a2´×u2 (4)
y=a1´×u1 (3)
y=a2´×u2 (4)
この時、u1=u2であり、u1=u2=20mm/hに対して、y=10m3/sを出力するモデルは、単位系をそのままで考えると、a1´=a2´=0.5であることは明らかである。従って、このような関係を持つ場合に、(3)式と(4)式と単回帰を適用すると、a1´=a2´=0.5に近い値がパラメータとして同定(学習)されることになる。この0.5という値は、(回帰モデルなので物理的な意味は明確ではないが)合理的で納得できる値である。
次に、(3)式と(4)式とを用いて、(1)式の重回帰モデルを構築することを考える。
今回の例では、入力変数が2か所の降雨量となり、その和をとっているため、入力となる降雨量が倍になるので、各降雨量に対する係数のa1とa2とは、一つの降雨量を入力とした単回帰モデルの係数の半分となり、a1=a2=0.25となると考える事が物理的に妥当な納得できるパラメータ値であると考えられる。
今回の例では、入力変数が2か所の降雨量となり、その和をとっているため、入力となる降雨量が倍になるので、各降雨量に対する係数のa1とa2とは、一つの降雨量を入力とした単回帰モデルの係数の半分となり、a1=a2=0.25となると考える事が物理的に妥当な納得できるパラメータ値であると考えられる。
このようにして、(3)式と(4)式とでu1とyとの関係、および、u2とyとの関係を単回帰で各々求めて、2つの予測モデルを構築する。そして、次に2つの予測モデルの出力yの平均を計算する。このように、構成的(Constructive)な手順を踏むと、a1とa2との値は、各々0.25と、物理的に妥当と考えられる納得のできる値として推定され、結果的に、(1)式の形の重回帰式の係数を定めることができる。これは、部分モデル(サブシステム)を結合していくモジュラー型のアプローチそのものである。
上記重回帰モデルにより係数を定めることができる理由は、一般に、入出力関係を表すデータとして同じデータを用いる場合、モデル構造が簡単(≒モデルに含まれるパラメータ数が少ない)なモデルの方が、可同定性が高くなる事(=パラメータ値を一意に同定可能)が知られているためである。本実施形態のモジュラー型時系列データ予測装置は、上記の特性を利用している。すなわち、同じ入出力データであったとしても、その可同定性はモデルの構造やパラメータ数に依存し、パラメータ数が多いモデルでは可同定でなくなる場合であっても、パラメータ数が少なければ可同定になる場合がある、ということを利用している。
上記の例でいえば、本実施形態のモジュラー型時系列データ予測装置では、2入力の重回帰モデルでは多重共線性により可同定でなくなるが、これを1入力の単回帰の組み合わせと考えることにより、多重共線性の問題を本質的に回避して、可同定なモデルを積み上げることによって合理的なモデル構築を可能にしている。
次に、本実施形態のモジュラー型時系列データ予測装置と、いわゆる「Explainable AI(XAI)」との関係および相違について説明する。
いわゆるXAIでは、全体を一つの大きなシステムとして捉えてブラックボックス的にモデルを構築して、後から、そのモデルを適切に解釈可能な形に分解して説明を加えようとするトップダウン型のアプローチ(あるいは帰納的なアプローチ)である。これに対し、本実施形態のモジュラー型時系列データ予測装置の考え方は、個別には比較的簡単に理解できる合理的な説明が可能なモデル(サブシステム)を統合していく事により、全体の予測モデルを構築していくボトムアップ型のアプローチ(あるいは演繹的なアプローチ)である。
いわゆるXAIでは、全体を一つの大きなシステムとして捉えてブラックボックス的にモデルを構築して、後から、そのモデルを適切に解釈可能な形に分解して説明を加えようとするトップダウン型のアプローチ(あるいは帰納的なアプローチ)である。これに対し、本実施形態のモジュラー型時系列データ予測装置の考え方は、個別には比較的簡単に理解できる合理的な説明が可能なモデル(サブシステム)を統合していく事により、全体の予測モデルを構築していくボトムアップ型のアプローチ(あるいは演繹的なアプローチ)である。
上記のようなボトムアップ型のアプローチでは、サブシステム間の相互の干渉を予め考慮することができないため、予測精度面に注目した場合には、トップダウン的(XAI的)なアプローチに劣る可能性を否定できない。しかしながら、ボトムアップ型のアプローチは、合理的な説明性という意味では、先に説明したようにパラメータの可同定性を維持したモデル構築が可能になり、トップダウン的なアプローチよりも優れる可能性が高い。
また、ボトムアップ型のアプローチは、ホワイトボックス的なモデリングのアプローチとも類似している。ホワイトボックスモデリングは、基本的に演繹的、構成的なモデリング手法である。ホワイトボックスモデリングでは、物理的に意味のあるパラメータを用いて構成的にモデルを構築するのに対し、本実施形態のモジュラー型アプローチでは、物理的な意味を持たせることはできないものの、個々のパラメータを入出力と直接関係づけることにより解釈が容易なパラメータを同定し、それを用いて構成的にモデル構築している。
本実施形態のモジュラー型アプローチがXAI的なアプローチよりも合理的な説明性において優れていることは以下の事からも明らかである。例えば、先の例において、もしトップダウン的(XAI的)なアプローチを行おうとすると、正則化やPLSなどの方法を使って係数パラメータのa1やa2を求め、それに対して何等かの解釈を加える処理を施して結果を「Explainable(説明可能)」にする方法の確立を目ざす。しかしながら、上記の簡単な例からわかる様に、a1とa2が0.25という値からかけ離れた値で同定されていた場合(実際、正則化やPLSを用いると、そのようになる可能性がある)には、いかなる方法を用いて「説明」を試みても、合理的な説明が原理的に困難であることは容易にわかる。例えば、PLSや正則化を用いてパラメータ値を求めることができても、a1かa2のいずれかが負の値(例:a1=3.5、a2=-3)となっている場合には、合理的に納得のできる説明は不可能である。
つまり、本来合理的に納得できる値である0.25という値を求められていなかったとしても、可同定性が欠如している場合には、ある制約さえ満たしていればどのような内部パラメータ値であろうが予測精度だけに着目する限りは差異が見られないため、「精度の良いモデル」と見なされる。しかしながら、精度が良いモデルが、必ずしもそのモデルのパラメータ値の妥当性を補償する合理的で妥当なモデルを意味しないため、トップダウン的(XAI的)なアプローチでは原理的に合理的な説明が困難になってしまう。
なお、XAI的なアプローチでは、a1とa2とのそれぞれに対する個別の意味付けは困難であるが、a1+a2に対する意味づけは可能であるため、説明不可能であるというわけではない。ただし、XAI的なアプローチでは、内部に含まれる全てのパラメータに対して意味付けを行うことは原理的に困難になる可能性が高く、何等かの「特徴量」を抽出して、それに対して意味付けを求めることになる。
従って、「合理的な説明性の高いモデル」を構築することを主目的にするのであれば、そもそも、正則化やPLSなどの数学的に悪条件を回避する様なアルゴリズムに頼る前に、可同定性が欠如する様な状態を回避してモデルの内部構造(内部パラメータ)の一意性を保証する様にモデルを構築しておき、合理的な説明性を担保しながらモデルを統合していくボトムアップ的、構成的、演繹的なモジュラー型アプローチの方が優れていると可能性が高い。
また、このように合理的説明性を担保したモデルを構築できることに加え、最初に述べた分離可能な構造とすることで、モデルの調整も容易にすることができ、これらを合わせて合理的説明可能かつ調整可能なモデル構築を行うことが可能になる。
本実施形態のモジュラー型時系列データ予測装置は、以上の洞察に基づいて成されたものである。
本実施形態のモジュラー型時系列データ予測装置は、以上の洞察に基づいて成されたものである。
以下、実施形態のモジュラー型時系列データ予測装置について、図面を参照して詳細に説明する。
図1は、一実施形態のモジュラー型時系列データ予測装置を適用した雨量流入予測システムを概略的に示す図である。
図1は、一実施形態のモジュラー型時系列データ予測装置を適用した雨量流入予測システムを概略的に示す図である。
以下では、本実施形態のモジュラー型時系列データ予測装置を搭載したシステムの実施イメージとその効果をより明確にするために、雨水流入予測システムを対象として説明する。なお、本実施形態のモジュラー型時系列データ予測装置が適用される対象プロセスは雨水流入予測システムに限定されるものではなく、2種類以上の計測項目の時系列データが存在する任意のプロセスに対して適用することができる。モジュラー型時系列データ予測装置は、例えば、鉄鋼プラント、石油化学プラント、食品プラント、製薬プロセス、半導体製造プロセス、発電プラント、交通監視設備、空調監視設備、などのプロセスに適用することができる。
雨量流入予測システムの対象プロセス1は、流量計11と、幹線流量計12と、幹線水位計131-13Kと、地上雨量計141-14Mと、レーダ雨量計15と、ポンプ井水位計161、162と、流入渠17と、雨水ポンプ井18と、雨水ポンプ19と、流入ゲート110と、を含む都市雨水排水プロセス1である。レーダ雨量計15は、Q×Pメッシュのメッシュ毎のレーダ雨量計1511~15QPを含む。
都市雨水排水プロセス1の各種センサ11-15、161、162は、所定の周期(例えば、30秒、60秒)で、プロセスの状態を表す量や運転操作に関わる量などの計測を行う。各種センサ11-15、161、162で計測された値は、モジュラー型時系列データ予測装置で収集および保存される。
流入渠17は、下水管やポンプ場から送られてきた雨水が流れ込む渠である。
雨水ポンプ井18は、流入渠17から流入した雨水が貯められる貯水槽である。雨水ポンプ井18に流入する雨水は、沈砂池で一緒に流れてきた砂などが予め取り除かれていてもよい。
雨水ポンプ19は、雨水ポンプ井18に貯められた雨水を強制的に河川等へ送り出すポンプである。
雨水ポンプ井18は、流入渠17から流入した雨水が貯められる貯水槽である。雨水ポンプ井18に流入する雨水は、沈砂池で一緒に流れてきた砂などが予め取り除かれていてもよい。
雨水ポンプ19は、雨水ポンプ井18に貯められた雨水を強制的に河川等へ送り出すポンプである。
流入ゲート110は、流入渠17と雨水ポンプ井18との間の流路に設けられ、開閉することにより流入渠17から雨水ポンプ井18に流入する雨水の量を調整するように動作する。流入ゲート110の前後には水位計161、162が設置され、例えば、水位計161、162により測定された値に応じて流入ゲート110の動作が制御される。
図2は、一実施形態のモジュラー型時系列データ予測装置の一構成例を概略的に示す図である。
本実施形態のモジュラー型時系列データ予測装置は、データ収集保存部2と、データ抽出部3と、ペアワイズ予測モデル同定部4と、予測モデル合成法定義部5と、ペアワイズ出力変数予測部6と、合成出力変数予測部7と、予測誤差評価部8と、ペアワイズ予測モデル修正部9と、出力予測結果観測部10と、を備えている。
本実施形態のモジュラー型時系列データ予測装置の構成は、例えば、CPUやMPUなどのプロセッサを少なくとも1つと、プロセッサにより実行されるプログラムが格納されるメモリとを備えた演算装置であって、種々の機能をソフトウェアにより、若しくは、ソフトウェアとハードウエアとの組み合わせにより実現することが出来る。
本実施形態のモジュラー型時系列データ予測装置は、データ収集保存部2と、データ抽出部3と、ペアワイズ予測モデル同定部4と、予測モデル合成法定義部5と、ペアワイズ出力変数予測部6と、合成出力変数予測部7と、予測誤差評価部8と、ペアワイズ予測モデル修正部9と、出力予測結果観測部10と、を備えている。
本実施形態のモジュラー型時系列データ予測装置の構成は、例えば、CPUやMPUなどのプロセッサを少なくとも1つと、プロセッサにより実行されるプログラムが格納されるメモリとを備えた演算装置であって、種々の機能をソフトウェアにより、若しくは、ソフトウェアとハードウエアとの組み合わせにより実現することが出来る。
データ収集保存部2は、各種センサ11-15、161、162で計測されたプロセス変数の値を収集し保存する。データ収集保存部2は、出力変数データ選択部21と、入力変数データ選択部22とを備えている。
出力変数データ選択部21は、各種センサ11-15、161、162で計測された値の少なくともいずれかの量を出力変数として選択する。
入力変数データ選択部22は、各種センサ11-15、161、162で計測された値から、複数の入力変数として選択する。なお複数の入力変数の数は、複数のプロセス変数の数以下であればよい。
入力変数データ選択部22は、各種センサ11-15、161、162で計測された値から、複数の入力変数として選択する。なお複数の入力変数の数は、複数のプロセス変数の数以下であればよい。
データ抽出部3は、データ収集保存部2に保存された時系列データの中から所定のデータを抽出する。データ抽出部3は、オフライン予測モデル同定用データ抽出部31と、オンライン予測用データ抽出部32と、評価用データ抽出部33とを備えている。
オフライン予測モデル同定用データ抽出部31は、データ収集保存部2から、所定の同定周期若しくは外部からの要求により、予測モデル同定用のデータとして各種プロセスセンサ11-15、161、162で計測された(当該時刻よりも過去の)所定の期間のデータを抽出する。
オンライン予測用データ抽出部32は、データ収集保存部2から、所定の予測周期で、各種プロセスセンサ11-15、161、162の計測値からオンラインの予測に必要となる予測用データを、リアルタイムで抽出する。
評価用データ抽出部33は、データ収集保存部2から、所定の評価周期もしくは外部からの要求により、予測誤差評価に必要となる各種プロセスセンサ11-15、161、162の計測値を、リアルタイムで抽出する。
ペアワイズ予測モデル同定部4は、オフライン予測モデル同定用データ抽出部31で抽出された計測値のデータを用いて、出力変数データ選択部21で選択された出力変数および入力変数データ選択部22で選択された入力変数を変数毎にペアワイズとし、後述する方法により予測モデルを同定する。
予測モデル合成法定義部5は、ペアワイズ予測モデル同定部4で同定したペアワイズの予測モデルの出力値(ペアワイズ予測値)から、出力変数(出力変数が複数の場合は出力変数毎)を合成する方法を定義する。
ペアワイズ出力変数予測部6は、ペアワイズ予測モデル同定部4で同定したペアワイズの予測モデルに対して、時間の進行方向における所定の監視周期で、オンライン予測用データ抽出部32から抽出した入力変数各々に対応する予測用データを入力して、入力変数毎に出力変数(出力変数が複数の場合は出力変数毎)のペアワイズの予測を行う。
合成出力変数予測部7は、ペアワイズ出力変数予測部6の出力であるペアワイズの予測値に対して、予測モデル合成法定義部5で定義した方法に従って合成した出力変数の予測値を出力する。
予測誤差評価部8は、ペアワイズ出力変数予測部6の出力であるペアワイズの予測値と、合成出力変数予測部7の出力であるペアワイズの予測値とを合成した予測値を保存し、所定の周期もしくは所定のタイミングで、評価用データ抽出部33により抽出した出力変数の実績値と保存したペアワイズおよび合成予測値との誤差を、各々評価する。
ペアワイズ予測モデル修正部9は、予測誤差評価部8によって評価した予測誤差に基づいて、ペアワイズ予測モデルの調整の必要性の有無を判断し、必要なペアワイズ予測モデルの修正を行う。ペアワイズ予測モデル修正部9は、例えば、ペアワイズ予測モデルのパラメータ調整、もしくは、ペアワイズ予測モデルの削除、もしくは、ペアワイズ予測モデルの統合、を行うことができる。ペアワイズ予測モデル修正部9は、修正したペアワイズ予測モデルをペアワイズ出力変数予測部6に供給する。
出力予測結果観測部10は、例えばモニタなどの表示手段と、ユーザインタフェースとを備え、合成出力変数予測部7による合成出力変数の予測値を少なくとも観測者に提示することができる。出力予測結果観測部10は、例えば、パーソナルコンピュータ、タブレット端末や、スマートフォンなどの携帯端末であってもよい。出力予測結果観測部10は、必要に応じて、ペアワイズ出力変数予測部6によるペアワイズ予測モデルの複数の出力と、それから合成される予測分布情報の少なくとも一方を、観測者に更に提示してもよい。
次に、本実施形態のモジュラー型時系列データ予測装置の動作について説明する。
例えば都市雨水排水プロセス1では、各種プロセスセンサ11-15、161、162によって、所定の周期でプロセスの情報が計測され、データ収集保存部2に計測値(プロセス変数の値)が供給される。データ収集保存部2は、収集した計測値を、予め決められたフォーマットに従って、時系列データとして保存する。
例えば都市雨水排水プロセス1では、各種プロセスセンサ11-15、161、162によって、所定の周期でプロセスの情報が計測され、データ収集保存部2に計測値(プロセス変数の値)が供給される。データ収集保存部2は、収集した計測値を、予め決められたフォーマットに従って、時系列データとして保存する。
出力変数データ選択部21は、各種プロセスセンサ11-15、161、162で計測されたプロセス変数の中から予測の対象となる変数を選択して、選択した変数を出力変数と設定する。
例えば都市雨水排水プロセス1では、通常雨水ポンプ井18への流入量を予測の対象とする場合には、出力変数データ選択部21は、流量計(雨水ポンプ井流入量計)11で計測されている雨水ポンプ井流入量が予測対象の出力変数として選択する。
例えば都市雨水排水プロセス1では、通常雨水ポンプ井18への流入量を予測の対象とする場合には、出力変数データ選択部21は、流量計(雨水ポンプ井流入量計)11で計測されている雨水ポンプ井流入量が予測対象の出力変数として選択する。
また、流量計(雨水ポンプ井流入量計)11が設置されていない様な場合には、雨水ポンプ井水位計161、162の計測値と、雨水ポンプ19および流入ゲート110の運用状態とから、雨水ポンプ井18への流入量を計算によって求めて、これを予測対象として選択しても良い。この場合、出力変数データ選択部21は、流入量の計算値を計算して保存しておき、流入量の計算値を出力変数として選択する。
また、例えば都市雨水排水プロセス1の複数の項目を予測する場合には、出力変数データ選択部21は、複数の項目を出力変数として選択してもよい。都市雨水排水プロセス1においても、例えば、雨水ポンプ井18に複数の箇所から雨水の流入があるような場合は、出力変数データ選択部21は、複数の流入箇所各々における流入量を出力変数として選択してもよい。
入力変数データ選択部22は、出力変数データ選択部21で選択した出力変数に影響を与える可能性のある測定値を入力変数の候補として選択する。影響を与えるか否かについて事前情報が全く無い場合には、入力変数データ選択部22は、全てのプロセス変数を機械的に入力変数としておいても良い。入力変数データ選択部22は、出力変数として選択した変数を入力変数としても選択して構わない。例えば、時系列のデータ解析では、出力変数として選択した変数自身の過去の値から未来の値を予測する所謂自己回帰を検討することも多くある。このような場合には、流量計(雨水ポンプ井流入量計)11による雨水ポンプ井流入量は、出力変数と入力変数との両方として選択されても構わない。
また、出力変数データ選択部21で選択した出力変数が複数ある場合、入力変数データ選択部22は、ある出力変数に対して別の出力変数を入力変数として選択してもよい。このようにすることで、複数の出力変数同士の相関を考慮することが可能になる。一方、複数の出力変数が異なる流入箇所の流入量であるときには、プロセスの物理的な構造上、出力変数同士の関連が無い場合が多い。この場合、入力変数データ選択部22は、ある出力変数に対する入力変数として別の出力変数を除外しておくこともできる。
入力変数データ選択部22は、その他の変数に対しても、例えば、物理的な管渠の接続関係などの情報から明らかに出力変数に無関係な変数や因果関係(入力(原因)と出力(結果)との関係)が物理的に成立しない変数を予め入力変数から除外しておくこともできる。
例えば、雨水ポンプ井水位計161、162により計測される雨水ポンプ井水位は、出力変数として選択した流量計(雨水ポンプ井流入量計)11による雨水ポンプ井流入量の影響を受けて変化することは物理的に考えて明らかであるから、雨水ポンプ井水位を入力変数から除外しておく。
入力変数データ選択部22は、物理的な構造の関係が明確でない変数について予め出力変数データ選択部21で選択した出力変数との相関解析を予備的に実施し、相関係数の絶対値が所定の値(例えば0.5)未満の変数を、予め入力変数から除外しておいてもよい。さらに、入力変数データ選択部22は、相関解析を実施する際に、解析対象の変数が測定された時間をずらしながら相関解析を行い、最も相関の高くなる時のずらした時間Lを求め、時間Lの符号情報から、因果関係を満たさない符号を持つ変数を入力変数から予め除外しておいてもよい。なお、先の雨水ポンプ井水位の変数は、このような解析によっても入力変数から除外される可能性が高い。
本実施形態のモジュラー型時系列データ予測装置では、入力変数データ選択部22は、雨水水位計161、162により計測された雨水ポンプ井水位のみを除いた流量計(雨水ポンプ井流入量計)11による雨水ポンプ井流入量と、幹線流量計12により計測された幹線流入量と、幹線水位計131-13Kにより計測された幹線水位と、地上雨量計141-14Mにより計測された地上雨量と、Q×Pメッシュのメッシュ1511-15QP各々におけるレーダ雨量とを入力変数として選択するものとする。
なお、本実施形態のモジュラー型時系列データ予測装置において、例えば、K個の幹線水位計131~13Kは、最近、下水管渠内の水位情報を正確に把握しようとする行政の動向などから、Kの値が10~100程度の比較的多数の幹線水位計が設置されることも多く、また、気象レーダは国土交通省が全国にXRAINという気象レーダを設置していることから、多くの場所で取り入れることができ、対象とする雨水排水区の大きさにもよるが、Q×Pメッシュのメッシュサイズは100程度になることも珍しくないため、実際に計測される時系列データの項目数(変数の数)は、かなり多くなる。また、本実施形態のモジュラー型時系列データ予測装置が適用される都市雨水排水プロセス1に限らず、計測変数の項目数が数百~数千に及ぶことは一般的である。
本実施形態のモジュラー型時系列データ予測装置では、入力変数は、例えば、雨水ポンプ井水位計161、162による計測値のみを除いた流量計(雨水ポンプ井流入量計)11による雨水ポンプ井流入量の計測値と、幹線流量計12による幹線流入量の計測値と、幹線水位計131~13KによるK個の幹線水位の計測値と、地上雨量計141~14MによるM個の地上雨量の計測値と、Q×Pメッシュ1511~15QPの各メッシュのレーダ雨量の計測値と、である。以下では、入力変数の数をp(正の整数)として説明する。
オフライン予測モデル同定用データ抽出部31は、予測モデルを構築する人が指定した所定期間のデータを用いて要求されたタイミングで、先に出力変数データ選択部21と入力変数データ選択部22との該当する所定期間分の時系列データを取得する。オフライン予測モデル同定用データ抽出部31は、予測モデルの構築を定期的に行い、モデルを定期的に更新したい場合には、時間の進行方向において所定の周期TLで所定期間の入出力変数の時系列データを定期的に抽出しても良い。このようにして抽出されたデータセットをuk、k=1、2、…pと記載する。ukは列ベクトルであるとし、各行が時系列データの各時刻のサンプルに対応する。同様に出力変数である雨水流入量の時系列データセットyとしておく、オフライン予測モデル同定用データ抽出部31は、データセットukとデータセットyとを所定の周期で抽出する作用を持つ。
次に、ペアワイズ予測モデル同定部4では、ukとyとを用いて、(u1、y)、(u2、y)、…(up、y)のp個のペア毎に予測モデルを同定する。本実施形態のモジュラー型時系列データ予測装置では、出力変数yは、流量計(雨水ポンプ井流入量計)11による雨水ポンプ井流入量のみの一つであるが、複数(h個とする)の出力yがある場合には、h個の各々のyに対して同様の処理を行えば良く、以降の処理は出力変数の数がl個に増えた場合でも、h回繰り返して同じ処理を行えばよいだけであるので、一般性を失う事なく出力変数は1個であるとしておく。
最も典型的な予測モデルは、次式で表せる線形回帰の形をしたモデルを仮定して同定を行うものである。
y(t)=a1×uk(t-L)+a2×uk(t-L-1)+…+an×uk(t-L-n+1)+c (5)
ここで、ak、k=1、2、…、nとcは同定すべきパラメータである。cは入力ukと出力yの平均の差を表すバイアスパラメータであるが、予めukとyから平均値を除去して平均を0としておくことにより、常に0とすることができるので、以下では簡単のためにcは0としておく。
y(t)=a1×uk(t-L)+a2×uk(t-L-1)+…+an×uk(t-L-n+1)+c (5)
ここで、ak、k=1、2、…、nとcは同定すべきパラメータである。cは入力ukと出力yの平均の差を表すバイアスパラメータであるが、予めukとyから平均値を除去して平均を0としておくことにより、常に0とすることができるので、以下では簡単のためにcは0としておく。
上記(5)式において、ukがy自身ではない場合を、ディジタル信号処理やシステム同定の分野では、有限インパルス応答モデル(FIR)モデルと呼ぶ(入力変数に対する重み付き移動平均モデルと呼ぶこともできる)。一方、ukがy自身と一致する場合は自己回帰モデル(ARモデル)と呼ばれる。ここで重要な事は(5)式の線形回帰の形であっても、複数の入力変数による有限インパル応答を足し合わせ、さらに自己回帰を足し合わせた多入力1出力(複数の出力を同時に考える場合は多入力多出力)の線形回帰の形(ARX(Autoregressive eXogenous Input)などの形)ではなくて、あくまでも1入力1出力の形の線形回帰モデルになっている点であり、これが「ペアワイズ」の意味である。
ペアワイズで同定することの第一の利点は、「複数の入力」という概念自身が無くなるため、先に述べた(複数の入力変数間の関係から生じる)多重共線性の問題が決して生じる事はなく、この問題を本質的に回避できることを含め、パラメータの可同定性の観点では多入力1出力のシステムよりも有利であるからである。これを若干補足する。システム同定の分野では、入力信号ukに含まれる周波数成分の量(数)によって同定可能なパラメータの数が変化することが持続的励振条件(PE(Persistent Exciting)条件)と呼ばれるパラメータの可同定性(=パラメータを一意に決定可能)条件として知られており、同定(推定)すべきパラメータakの数が多くなるほど、入力信号ukは多数の周波数成分を含まなければならない事が知られている。
一方、本実施形態のモジュラー型時系列データ予測装置において、入力変数の値は、予測モデルの構築者が制御(調整)できるものではなく、所与の観測(計測)情報であるから、どれくらいの周波数成分を含んでいるか否かは、外部条件として与えられるものであり、パラメータの数は少なければ少ないほど可同定性が向上し、逆にパラメータの数が多ければ多いほど可同定性が劣化する可能性が高くなる。より厳密に述べると、パラメータの数が少ない場合に可同定性が悪くなる可能性は0%であり、パラメータ数が少ないものの可同定性の方が高いか、悪くてもパラメータ数が多いものの可同定性と同じである。従って、上記(5)式の1入力1出力のモデルではn個のパラメータを含むのに対し、p入力1出力の予測モデルを考えるとp×n個のパラメータを含むことになるため、入力数pが多くなればなるほどパラメータを一意に同定することが難しくなる可能性が高くなる事がわかる。
さらに、先に述べた入力変数間の従属関係(相関関係)があると、多重共線性の問題が生じ、同定すべきパラメータ数の多寡に関わらず可同定性が失われ、入力情報は所与のものであるから、いかなる方法を用いても原理的にパラメータ値を一意に決定することが不可能になる。これに対し、1入力1出力の(5)式のモデルでは多重共線性の問題が原理的に生じることはなく、(5)式の形でパラメータが同定できなくなる場合は、入力ukに十分な周波数情報が含まれていない時に限られ、この場合は、次数nを調整して同定すべきパラメータ数を調整することで、パラメータakを必ず可同定(一意に同定可能)にすることができる。
例えばukの値が一定値である場合、uk(t-L)=uk(t-L-1)=…=uk(t-L-n+1)となるが、このような場合はn=1として一つのパラメータa1だけを同定する様にすればパラメータの可同定性を維持することができる。一般的には、所与のukを用いて上記の持続的励振条件を調べることで同定可能なパラメータ数nの上限値nmaxを求めることができるので、n≦nmaxの範囲内で、例えば、AIC(赤池情報量規範)やBIC(ベイズ情報量規範)、MDL基準(最小符号化基準)などの各種の規範および基準を用いたり、交差検証(クロスバリデーション)を行ったりすることで、nの適正値を決めて、(5)式の可同定性を常に維持するこができる。
ペアワイズで同定することの第二の利点は、たとえ、出力変数を除く入力変数が1つの場合であっても、この方法では出力変数を入力とする自己回帰ARモデルと入力変数(出力変数ではない入力変数)を入力とするFIRモデルとの二つのモデルを分離して、同定することによる効果である。このような場合、ペアワイズでない通常の方法では、ARXモデルなどを用いて、自己回帰成分(AR成分)と重み付き移動平均成分(FIR成分/MA成分)とを同時に同定することになるが、多くの場合、出力変数と入力変数との関係(相関)よりも出力変数自身の相関(自己相関)の影響の方が強い場合が多く、予測に対する自己回帰成分の影響が極めて強くなる場合が多い。このような場合、入力変数の値が急変した場合にその変化に追従できなくなり、予測が遅れるという現象が生じやすい。
本実施形態のモジュラー型時系列データ予測装置の場合、例えば、出力変数はポンプ井流入量であり、入力変数として1か所の地上雨量を仮定すると、流入量に対する自己回帰成分の影響の方が地上雨量に対するFIR成分の影響よりも強くなり、降雨の変化により地上雨量の値が急変した場合でも流入量に変化が現れるまでは予測値の変化が小さく、予測値が実測値に対して遅れるという現象が生じやすい。これに対し、流入量に対する自己回帰モデル(ARモデル)と地上雨量に対する有限インパルス応答モデル(FIRモデル)とを各々個別に同定するペアワイズの方法では、監視員等のユーザは、予測に対するAR成分の影響とFIR成分の影響とを調整することができる。
ペアワイズで同定することの第三の利点は、第二の利点と密接に関係するが、各入力変数に対する(1)式のモデルの意味が明確であることである。ここで、先と同じ様に、出力変数を除く入力変数が1つの場合である下記(2)式の単純なARXモデルを考えてみる。
y(t)=a1×y(t-1)+ a2×y(t-2)+…+any×y(t-ny)
+b1×u(t-L)+ b2×u(t-L-1)+…+bnu×u(t-L-nu+1) (6)
y(t)=a1×y(t-1)+ a2×y(t-2)+…+any×y(t-ny)
+b1×u(t-L)+ b2×u(t-L-1)+…+bnu×u(t-L-nu+1) (6)
この時、上記(6)式をARXモデルとして同定すると、以下のyに関する項(Yと定義する)とuに関する項(Uと定義する)は、YとUの意味は明確ではなく、YとUとを個別に意味付けして説明することは困難である。
Y:=a1×y(t-1)+ a2×y(t-2)+…+any×y(t-ny) (7)
U:=b1×u(t-L)+ b2×u(t-L-1)+…+bnu×u(t-L-nu+1) (8)
一方、(7)式と(8)式とは、各々(5)式の形をしているので、ペアワイズで同定を行うと、YとUとは各々yの予測値を意味することになるので、YとUとは、各々別の説明変数(入力変数)に対する出力yの予測値であるという解釈をすることができる。
U:=b1×u(t-L)+ b2×u(t-L-1)+…+bnu×u(t-L-nu+1) (8)
一方、(7)式と(8)式とは、各々(5)式の形をしているので、ペアワイズで同定を行うと、YとUとは各々yの予測値を意味することになるので、YとUとは、各々別の説明変数(入力変数)に対する出力yの予測値であるという解釈をすることができる。
このようにペアワイズの同定を行うと、(7)式と(8)式とから、(6)式の形の予測モデルを構成的(Constructive)に構築することは容易であり、例えば(7)式のUと(8)式のYとの平均(=(U+Y)÷2)として構成することで、(6)式の形の予測モデルを得る事ができ、これは、再度yの予測値を意味していることは明らかである。
上記のように、予測モデルの形(構造)が同じであっても同定の手順を変える事で、その意味づけは変化する。ペアワイズで同定することの利点は、意味付けのしやすいモデルを構築することができる点である。
以上が多入力1出力ではなく1入力1出力のペアワイズで予測モデルを同定することの利点である。加えて、ペアワイズ予測モデルとして、(5)式のFIRもしくはARの形の回帰式の形で同定することで、ディジタル信号処理やシステム同定分野で従来から良く知られている周波数解析や安定解析などの様々な解析手法を、ペアワイズに直接流用できるため、その分野の知識があれば、ペアワイズでより詳細な解析や調整も可能になる。
以上が多入力1出力ではなく1入力1出力のペアワイズで予測モデルを同定することの利点である。加えて、ペアワイズ予測モデルとして、(5)式のFIRもしくはARの形の回帰式の形で同定することで、ディジタル信号処理やシステム同定分野で従来から良く知られている周波数解析や安定解析などの様々な解析手法を、ペアワイズに直接流用できるため、その分野の知識があれば、ペアワイズでより詳細な解析や調整も可能になる。
ペアワイズ予測モデルとして、(5)式の様なARモデル/FIRモデルをさらに分解し、時間遅れを考慮した変数を各々一つの説明変数と見なした単回帰モデルに分解して予測モデルを構築することもできる。
y(t)=a×uk(t-L)+c (9)
上記(9)式の様な単純な形にすると、(9)式に含まれるパラメータは、遅れ時間L、入力変数ukに対する出力変数の倍率(ゲイン、比例係数)a、入力変数と出力変数の平均値の差(入力が0の場合の出力の値、バイアス)cの三つだけとなり、全てのパラメータの意味が明確になる。
y(t)=a×uk(t-L)+c (9)
上記(9)式の様な単純な形にすると、(9)式に含まれるパラメータは、遅れ時間L、入力変数ukに対する出力変数の倍率(ゲイン、比例係数)a、入力変数と出力変数の平均値の差(入力が0の場合の出力の値、バイアス)cの三つだけとなり、全てのパラメータの意味が明確になる。
(5)式の形では、遅れ時間Lとバイアスcの解釈は(9)式と同様であるが、ak、k=1、2、…、nという回帰係数(パラメータ)の個々の意味までは明確でなく、システム同定やディジタル信号処理に関する知識を用いて、(5)式自身を解釈することはできるものの、個々のakの直感的な解釈は困難である。これに対し、(9)式の単回帰式の場合はaの解釈も容易であり、単純に、平均値を0に処理した場合の入力変数に対する出力変数の倍率を示している。これにより、全てのパラメータの意味を容易に解釈できる様になるため、(9)式を用いると、何等かのアルゴリズムを適用した場合のパラメータ同定結果が直感的な感覚とずれるような場合(例えば、異常データが混入している場合などにそのような事が起こり得る)、例えば予測モデルの構築者が直感的に手動でパラメータを調整することも容易になり、また、予測モデルの中身を他者に説明することも極めて容易になる。
本実施形態のモジュラー型時系列データ予測装置において、例えば、出力変数はポンプ井流入量であり、入力変数として1か所の地上雨量であると仮定すると、(9)式のcは地上雨量がゼロの場合の流入量を意味し、これは、汚水と雨水が同一管渠を通る合流式下水の場合は、汚水量に相当するものであるという解釈ができる。そして、(9)式のaは流入量から汚水量を引き去った場合に、地上雨量計による雨量に対する雨水流入量の倍率を示しており、地上雨量と流入量との相関が十分に高ければ、このような単純な倍率も意味を持つ。さらに、(9)式の遅れ時間Lは、地上雨量の観測点から雨水ポンプ井までの雨水の流れの遅れ時間を意味し、下水管に流出するまでの流出の遅れと下水管路を流れて雨水ポンプ井に到達するまでの流下時間の和を意味すると解釈することができる。
このように、解釈や説明が極めて容易になることが(9)式の様な単純なモデルを用いることの利点であり、実際に(9)式の形だけを用いた最も簡単なモデルは、流入量予測の簡易モデルとして、特許文献1の中でも採用されている。
上記のように直感的な解釈が可能になることにより、手動でのパラメータ調整が容易になる具体的な例について説明する。
図3および図4は、幹線水位データと流入量データとの関係の一例を概略的に示す散布図である。
図3および図4は、幹線水位データと流入量データとの関係の一例を概略的に示す散布図である。
図3および図4において、横軸は、幹線水位計131~13Kのいずれかに相当するある箇所の幹線水位データ、縦軸は流入量計11に相当する予測対象となる流入量データとした散布図の一例を示している。
図3に示す直線は、最小2乗法によって、単回帰のゲインとバイアスを求めたものであり、この場合ゲインa=9.6024、バイアスc=23.3142となっている。これは、予測値の平均2乗誤差(MSEもしくはRMSE)が最小になるという意味で、最適な回帰直線である。
図3に示す直線は、最小2乗法によって、単回帰のゲインとバイアスを求めたものであり、この場合ゲインa=9.6024、バイアスc=23.3142となっている。これは、予測値の平均2乗誤差(MSEもしくはRMSE)が最小になるという意味で、最適な回帰直線である。
図3において、例えば幹線水位の値が-1.5であるとき、回帰直線に従うと流入量の値はおおよそ9程度と予測されることになる。しかしながら、図3の散布図によれば幹線水位の値が-1.5であるときの流入量は、9よりも大きい値である場合が圧倒的に多いことがわかる。
上記のように、回帰直線による予測値と測定値との間にズレが生じる理由は、幹線水位が-2.5乃至-2付近のデータ数が圧倒的に多いため、この付近のデータにより適合する様に係数が推定されるためである。この回帰直線は、理論上は正しいものであるが、実際の流入量予測においては、必ずしも好ましくない場合が多い。なぜなら、流入量の予測値は、雨水排水ポンプの制御や運転支援に用いられるため、浸水リスクを回避するためには、流入量が多くなる場合をできる限り正確に予測したいという暗黙の要請があるためであり、流入量が少ない場合の予測精度はそれほど重要にならない場合が多いためである。
しかし、このような暗黙の要請をアルゴリズムに組みこむ場合には、流入量が多い箇所の予測誤差に重みをつけるなどの処置が必要であり、その重みの調整なども含め、一般にかなりの労力を要する。
図4では、最小2乗法と、最小2乗法の他の10種類の様々な回帰手法とを、図3と同じデータに適用した結果得られた回帰直線L1-L11の一例を示している。この中で、1番乃至10番の方法は、Pythonと呼ばれるソフトウェアのパッケージに組み込まれている様々なアドバンストな回帰手法であるが、いずれの方法による回帰直線L1-L10も、流入量が多いときに適切な予測値が得られるものではなかった。
11番目の方法は、幹線水位が-2.2乃至-1.5までの範囲のデータだけを抽出し、他のデータを使用せずに、通常の最小2乗法を適用した方法である。これは、本願発明者らが、幹線水位-2.5乃至-2付近のデータに適合するように回帰直線が求められることを理解した上で、その付近のデータを意図的に捨てることによって、流入量が多い場合のデータに適合するように意図的に調整した方法である。
このように、アルゴリズムに頼った調整を行おうとすると、例えばデータの範囲を限定する、あるいは、データを間引くなどの様々な工夫を施さないと、本来望むパラメータが得られないことがある。しかし、11番目の方法により得られる直線L11は、2点を決めれば直線が唯一に定まるという極めて単純な原理を用いれば、人間が目視で直線L11のようなラインを引くことにはほとんど労力を要しない。上記のように、単回帰という極めてわかりやすく、また図3や図4に示す散布図上に明確に可視化できるパラメータを用いることによって、手動での調整が極めて容易にできることがわかる。
すなわち、単回帰モデルを用いると、倍率(ゲイン)とバイアスという容易に理解できるパラメータを用いて解釈できて、調整が容易になる。
すなわち、単回帰モデルを用いると、倍率(ゲイン)とバイアスという容易に理解できるパラメータを用いて解釈できて、調整が容易になる。
さらに、単回帰を用いると別の観点で解釈することも可能になる。すなわち、統計分野における基本的な統計量である、平均、標準偏差、および、相関、という三つのパラメータを用いて解釈することが可能になる。
単回帰モデルの回帰係数(比例係数、傾き、倍率、ゲイン)aは、相関係数rと次式の関係にあることが広く知られている。
a=r×(σy/σu) (10)
ここで、r、σy、σuは、各々、入力変数と出力変数の相関係数、出力変数の標準偏差、入力変数の標準偏差、である。また、回帰係数aとバイアスcは、単回帰モデルの予測2乗誤差が最小となる様に最小2乗法で同定した場合、以下の関係になることも広く知られている。
c=μy-a×μu (11)
ここで、μy、μuは、各々、出力変数の平均と入力変数の平均、を表す。
a=r×(σy/σu) (10)
ここで、r、σy、σuは、各々、入力変数と出力変数の相関係数、出力変数の標準偏差、入力変数の標準偏差、である。また、回帰係数aとバイアスcは、単回帰モデルの予測2乗誤差が最小となる様に最小2乗法で同定した場合、以下の関係になることも広く知られている。
c=μy-a×μu (11)
ここで、μy、μuは、各々、出力変数の平均と入力変数の平均、を表す。
従って、(9)式において、入力変数と出力変数を各々、平均と標準偏差を用いて正規化しておけば、回帰係数aと相関係数rの値は一致し、c=0となる。
この関係が意味することは、単回帰モデルのゲインとバイアスでモデルを説明することと、各変数(入力変数と出力変数)の平均、標準偏差と、入力と出力との相関係数、のみでモデルを説明することは等価であるということである。従って、予め、変数毎に適切な平均と標準偏差の推定値を用いて各変数のデータを正規化しておけば、単回帰の入出力関係は、次式のように、相関係数のみで、全てを説明できることになる。
y´(t)=a×uk´(t-L)=r×uk´(t-L) (12)
ここで、y´=(y-μy)/σy、uk´=(uk-μuk)/σukで定義される正規化された出力変数と入力変数である。
この関係が意味することは、単回帰モデルのゲインとバイアスでモデルを説明することと、各変数(入力変数と出力変数)の平均、標準偏差と、入力と出力との相関係数、のみでモデルを説明することは等価であるということである。従って、予め、変数毎に適切な平均と標準偏差の推定値を用いて各変数のデータを正規化しておけば、単回帰の入出力関係は、次式のように、相関係数のみで、全てを説明できることになる。
y´(t)=a×uk´(t-L)=r×uk´(t-L) (12)
ここで、y´=(y-μy)/σy、uk´=(uk-μuk)/σukで定義される正規化された出力変数と入力変数である。
このように解釈すると、例えば、入出力の相関が無い場合、すなわちr=0となる様な入出力関係を持つ場合は、その出力に対して対応する入力は全く予測能力を持たない事が(11)式の関係から明確に理解することができる。また、現実のデータは、様々なノイズやアウトライア(外れ値)で汚染されている(コンタミされている)場合も多いが、このような時、平均、標準偏差、相関係数、の三つのパラメータを外れ値に対してロバストに推定するロバスト推定法を用いて推定することにより、外れ値に対してロバストな予測モデルを構築することもできる。もちろん、より複雑なモデルのパラメータ推定に対する様々なロバスト推定法も開発されているが、一般にロバスト推定した結果を解釈することは容易ではないのに対し、平均、標準偏差、相関係数という三つのパラメータに対するロバスト推定の場合は、何等かのロバスト推定を適用した推定結果の良否の判断や解釈も容易になるため、ノイズやアウトライアにコンタミされたデータに対する予測モデル構築を行いたい場合の解釈や調整も容易になる。
このように、(9)式や(10)~(12)式は、容易に理解しやすいパラメータのみで構成されているが、(5)式と比較しても、極めて単純な構造をしており、実際の応用においては、その表現能力が十分でなく、十分な予測精度が得られない可能性があった。そこで、(9)式(あるいは(10)~(12)式)の形から、(5)式の形を合成することを検討する。(9)式のモデルを用いて(5)式の形を得るためには、(9)式の遅れ時間Lを可変にして足し合わせればよい。すなわち、まず、(9)式における遅れ時間L((5)式のLと区別するため以下ではL´とする)の最小の値を(5)式のLと対応させてL´=Lとし、L以上の遅れ時間について、L´=L+1、L´=L+2、…L´=L+n-1としたn個の(9)式のモデルを加え合わせれば(5)式の形のモデルが得られる。
すなわち、以下の(13)式のn個の単回帰モデルを加え合わせれば、形式的に(5)式と等価な(14)式のFIRもしくはARモデルの形のモデルが得られる。
y(t)=a1×uk(t-L)+c (13_1)
y(t)=a2×uk(t-L-1)+c (13_2)
y(t)=an×uk(t-L-n+1)+c (13_n)
y(t)=1/n×(a1×uk(t-L)+a2×uk(t-L-1)+…+an×uk(t-L-n+1))+c (14)
y(t)=a1×uk(t-L)+c (13_1)
y(t)=a2×uk(t-L-1)+c (13_2)
y(t)=an×uk(t-L-n+1)+c (13_n)
y(t)=1/n×(a1×uk(t-L)+a2×uk(t-L-1)+…+an×uk(t-L-n+1))+c (14)
異なるn個の遅れ時間のn個の単回帰モデルを作成した上で、その平均モデルを作成すれば、形式的には(1)式と同じFIRあるいはARモデルの形のモデルが得られる。この際、上記の様に単純に平均化処理を行うこともできるが、先の述べた単回帰のゲインとバイアスと相関係数との関係を用いて、以下の様に重み付き平均を行うこともできる。
まず、予め、入出力変数を正規化して、(9)式を(12)式の形で表しておく。すると、回帰係数は相関係数と一致するので、この回帰係数(の絶対値)で重みづけした重み付き平均値として(14)式に類似した(15)式のFIR/ARモデルの形の式を得ることができる(以下では、正規化した変数と正規化していない変数との記号を区別せず、文脈に応じてu(t)、y(t)は正規化した変数を表すこととする。)
y(t)=(a1×uk(t-L)+…+an×uk(t-L-n+1))/(|a1|+|a2|+…+|an|)) (15)
y(t)=(a1×uk(t-L)+…+an×uk(t-L-n+1))/(|a1|+|a2|+…+|an|)) (15)
このようにすると、相関が強くなる遅れ時間による予測に対して重みをつけたペアワイズ予測モデルの合成が可能になる。なお、回帰係数と相関係数の関係が(12)式の様に陽に関係つけられない場合でも相関係数絶対値もしくは回帰係数の絶対値で重みづけ平均化処理を行うことは可能であるが、(12)式の関係があることで、重みづけを行う事の意味がより明確で説得性の高いものとなる。また、(15)式において、相関係数の絶対値ではなく相関係数の2乗によって重みづけを行っても良い。
ペアワイズ予測モデルとして、上記の二つの合成法は、いずれの方法でも入出力間の線形の関係しか表現することができないが、入出力間に非線形の関係がある場合も多い。このような場合は、(5)式の代わりに、適当な非線形関数φ(・)を用いて、以下の(16)式の様な非線形回帰の形として非線形変換すればよい。
y(t)=a1×φ(uk(t-L))+…+an×φ(uk(t-L-n+1))+c (16)
これは、機械学習の分野で広く知られている非線形回帰のテクニックであり、パラメータak、k=1、2、…、nに関する線形性さえ維持していれば、線形回帰の様々な手法を直接適用することができる。また、機械学習の分野でよく知られている様に、非線形関数φ(・)を直接指定しなくても、(16)式を計算する際に必要となる計画行列と呼ばれる行列に対して、φ(・)を指定することと等価に変換できるカーネル関数(類似度関数)を直接指定しても良い。ただし、その場合には、直感的な解釈性が若干低下する可能性がある事には注意する必要がある。
これは、機械学習の分野で広く知られている非線形回帰のテクニックであり、パラメータak、k=1、2、…、nに関する線形性さえ維持していれば、線形回帰の様々な手法を直接適用することができる。また、機械学習の分野でよく知られている様に、非線形関数φ(・)を直接指定しなくても、(16)式を計算する際に必要となる計画行列と呼ばれる行列に対して、φ(・)を指定することと等価に変換できるカーネル関数(類似度関数)を直接指定しても良い。ただし、その場合には、直感的な解釈性が若干低下する可能性がある事には注意する必要がある。
(16)式の様な非線形回帰を用いる事で、ある入力変数と出力変数との間に非線形関係がある場合にも対応することが可能になる。なお、ペアワイズ予測モデルの同定では、入力変数毎に出力変数との回帰モデルを構築するので、ある入力変数と出力変数の関係は非線形であるが、別の入力変数と出力変数との関係は線形であるような場合には、入力変数毎に(5)式と(16)式とを単に使い分けるだけでよく、これによって、一般的には説明性も向上することが期待される。以上の一連の作用が、本実施形態におけるペアワイズ予測モデル同定部4の作用である。
次に予測モデル合成法定義部5では、ペアワイズ予測モデル同定部4で定義した予測モデルの合成法を定義する。
最も簡単な予測モデルの合成法は、ペアワイズ予測モデルで各入力変数に対して出力変数を予測するモデルが構築されているため、その平均により予測モデルの合成を行う方法である。この場合は、予測モデル合成法定義部5では、ペアワイズ予測モデル同定部4のp個の予測出力(以下では、各入力変数u1、u2、…、upに対する予測出力をy1、y2、…、ypとする。)を平均化する次式を定義することになる。
y(t)=mean(y1(t)、y2(t)、…、yp(t))=1/p×(y1(t)+y2(t)+…+yp(t)) (17)
(17)式の定義は、最も基本的な定義方法であるが、この定義を改良することにより、合成した予測出力に対する信頼性を向上させたり、意図的に予測に傾向(バイアス)を持たせた予測を行ったりすることが可能になる。これを以下に順に説明する。
最も簡単な予測モデルの合成法は、ペアワイズ予測モデルで各入力変数に対して出力変数を予測するモデルが構築されているため、その平均により予測モデルの合成を行う方法である。この場合は、予測モデル合成法定義部5では、ペアワイズ予測モデル同定部4のp個の予測出力(以下では、各入力変数u1、u2、…、upに対する予測出力をy1、y2、…、ypとする。)を平均化する次式を定義することになる。
y(t)=mean(y1(t)、y2(t)、…、yp(t))=1/p×(y1(t)+y2(t)+…+yp(t)) (17)
(17)式の定義は、最も基本的な定義方法であるが、この定義を改良することにより、合成した予測出力に対する信頼性を向上させたり、意図的に予測に傾向(バイアス)を持たせた予測を行ったりすることが可能になる。これを以下に順に説明する。
まず、(17)式の様な単純な平均化処理(標本平均)を行うと、各ペアワイズ予測モデルの予測値が平等に扱われているため、予測精度の良いペアワイズ予測モデルの予測値と予測精度の悪いペアワイズ予測モデルの予測値が混合されて、合成した全体の予測精度が劣化してしまう可能性がある。また、実際の運用においては、ある入力変数の計測データの信頼性が低く、別の入力変数の計測データの信頼性が高いという場合も稀ではないが、信頼性の高い入力変数に対するペアワイズ予測モデルの予測値も信頼性の低い入力変数に対するペアワイズ予測モデルの予測値も平均化してしまうことで、合成した予測出力の予測精度が劣化してしまう。極端な場合、例えば、ある入力変数のセンサの故障や不具合などにより、当該入力変数の時系列データにアウトライア(外れ値、異常値)が多量に含まれる様になる場合などには、予測精度が劣化するだけでなくアウトライアに引きずられて合成した予測自身が無意味になる(破綻してしまう)可能性がある。
このような状況に対応するために予測モデル合成法定義部5で、以下の様な複数のアプローチをすることができる。
一つ目は、主に後者のアウトライアに対して合成した予測値自身が破綻しない事を重視するアプローチであり、(17)式の標本平均処理に変えて、ロバスト推定を採用する方法である。
一つ目は、主に後者のアウトライアに対して合成した予測値自身が破綻しない事を重視するアプローチであり、(17)式の標本平均処理に変えて、ロバスト推定を採用する方法である。
(17)式で用いる「標本平均」という処理は、y1、y2、…、ypのp個の予測値の中から代表値を推定する方法の一つであり、このような代表値の推定を統計分野では位置母数の推定と呼ぶ。すなわち、「標本平均」は、位置母数の推定方法の一つであり、統計的には、標本平均は推定効率の観点では良い性質を持つことが知られているが、一方でアウトライアに対するロバスト性の面では最も脆弱(非ロバスト)であることも知られている。従って、外れ値に対するロバスト性を向上させるためには、(17)式をロバストな位置母数推定で置き換える方が良い。
ロバスト推定の分野では外れ値に対するロバスト性を評価するいくつかの指標が知られているが、最も直感的に理解しやすい指標としてブレークダウンポイントという指標がある。これは、統計的推定に用いるデータの中に何パーセントアウトライアが混入することを許容するかという指標であり、用いるデータのX%を、仮想的なアウトライアを想定して∞に置き換えた場合、推定量(平均などの位置母数を推定する場合は位置母数の推定値)が∞になる(=破綻する)かならないかの境界(∞になる直前)のXの値をブレークダウンポイントと呼ぶ。
例えば、「標本平均」という位置母数の推定に対しては、ただ1点のデータを∞に置換するだけでその標本平均値も∞になるので、「標本平均」のブレークダウンポイントは0%である。ロバスト統計ではブレークダウンポイントの最大値は50%であることが知られており、最大のブレークダウンポイントを持つ位置母数推定量として「中央値(メジアン)」が知られている。従って、予測値が外れ値に影響されないようにすることを最大の目的とする場合には、(17)式の「平均」を「中央値」で置換して、予測モデル合成法定義部5の定義とすることもできる。
一方、「中央値」はロバストではあるが、推定効率が悪い事が知られており、直感的にもp個の変数で予測している中の一つの予測値しか用いないため、精度の高い予測を行うことを重視する場合に最良の方法ではないことは容易に推測できる。そのため、推定効率の向上とロバスト性の向上を両立するための各種のロバスト推定方法が知られている。最も単純な方法は、「トリム平均(刈込平均)」と呼ばれる処理であり、p個のデータの上位と下位とのα%を削除した上で平均をとる方法である。直感的に明らかな様にαを大きくするとロバスト性は向上し、その極限として中央値推定があり、αを小さくして0に近づけるとロバスト性が低下しその極限が通常の平均(標本平均)であることは明らかである。従って、このトリム平均を予測モデル合成法定義部5の定義とすることもできる。この場合、αの調整によってロバスト性を調整できるが、適切に調整する事自身が難しい場合も考えられるため、別のロバスト推定方法をとることもできる。
このような方法の代表的な例として、p個のy1、y2、…、ypの中からすべての組み合わせの(p(p-1)/2個)の、二つのyiとyj(i≠j)との平均を計算した上で、その中央値を採用するホッジスレーマン推定量(HL推定)と呼ばれる位置母数推定を行うこともできる。この方法は、ブレークダウンポイントが約30%でロバスト性が高いうえに推定効率も良いことが知られており、このHL推定を予測モデル合成法定義部5の定義とすることもできる。
他のロバスト推定法として、一般には多変量データに対して用いられる手法であるMCDと呼ばれる手法が知られているがこれを適用することもできる。MCDは、p個のデータの所定の割合(通常50%~75%)のデータを取り出し、そのすべての組み合わせの中で分散が最小になるデータを用いて推定を行う手法であり、この方法を用いて平均を推定することもできる。MCDもロバスト性が高く推定効率の良い方法として知られている。一方、HL推定やMCD推定は、p個のデータの中から取り出すデータの数(HL推定場合2個、MCDの場合はp/2~3p/4個程度)の全ての組み合わせに対する計算が必要となるため、ペアワイズモデルの入力変数の数pが増加すると処理時間が飛躍的に増加するため、予測モデル合成法定義部5の定義として用いても、実際に予測を行う際にリアルタイム性を確保できない可能性がある。
このような場合に対応するため、「全ての組み合わせ」に対して処理を行うのではなく、リアルタイムで処理が可能な回数だけ繰り返し処理を行う様にブートストラップ法を用いて代用しても良い。すなわち、p個のy1、y2、…、ypの中から、所定のk個のデータをランダムに繰り返し抽出し、所定の回数(N回)、その平均値を求める。繰り返し求めたN個の平均値の中から、HL推定と同じように、例えばN個の平均値の中央値を採用する。このようなブートストラップ法による平均値推定を、予測モデル合成法定義部5の定義とすることで、リアルタイム性を維持しながら、推定効率が良くロバスト性の高い予測値の合成を行うことができると考えられる。
上記以外にも外れ値に重みを付けて推定を行うM推定などのロバスト推定を位置母数推定に適用した方法を予測モデル合成法定義部5の定義とすることもできる。
上記以外にも外れ値に重みを付けて推定を行うM推定などのロバスト推定を位置母数推定に適用した方法を予測モデル合成法定義部5の定義とすることもできる。
二つ目は、先のアプローチが主に外れ値(外れた予測値)に対してロバストに予測値の代表値を推定することを目的としていたのに対し、実際に予測精度の良い予測値に重みを付けて合成した予測を行うことを目的としたアプローチである。基本的な考え方は、(17)式の単純な標本平均に変えて、重み付き平均を行う方法であり、次式で表される式で合成予測出力を定義する。
y(t)=w1×y1(t)+w2×y2(t)+…+wp×yp(t) (18)
ここで、wk、k=1、2、…、pは、重みであり、w1+w2+…+wp=1となる制約を満たす。予測モデル合成法定義部5の定義では、この重みの設定法を定義する必要がある。
y(t)=w1×y1(t)+w2×y2(t)+…+wp×yp(t) (18)
ここで、wk、k=1、2、…、pは、重みであり、w1+w2+…+wp=1となる制約を満たす。予測モデル合成法定義部5の定義では、この重みの設定法を定義する必要がある。
これには、ペアワイズ予測モデル同定部4で、ペアワイズ予測モデルを同定した際の予測出力と実際の出力の間の(重)相関係数、あるいは、その2乗である決定係数を用いることができる。これにより、少なくともペアワイズの予測モデル同定時のデータに対して、予測精度の高い予測能力を持つモデルの予測出力を重視した重み付き平均で予測出力を合成することができる。
また、ペアワイズの予測モデルを、例えば上述の(15)式のように構築した場合には、各遅れ時間毎の説明変数と出力変数の間の相関係数が既に算出されているので、統合したペアワイズ予測モデル((14)式に相当)毎に、各相関係数の絶対値の平均値rk、k=1、2、3、…、pを求め、rkに応じて、重みwkをwk=rk/(r1+r2+…+rp)と定義することもできる。
このように、ペアワイズ予測モデルの予測能力に応じて重み付き平均値の重みを決定して(18)式により、予測モデルの合成法を定義する動作が、本実施形態における予測モデル合成法定義部5の動作の一例である。
三つ目は、二つ目と同様に(17)式の、各ペアワイズ予測モデルの重み付き平均値により、合成した予測出力を定義するが、この重みを、データを用いて同定(推定)することで決定するものである。この際、ペアワイズ予測モデルの同定時に用いたデータと同定したパラメータ値を用いると、その同定データに対する各ペアワイズ予測モデルの予測値y1、y2、…、ypの時系列データが得られる。
そして、(17)式を、この予測値y1、y2、…、ypを入力として、同定に用いた実際の出力データを出力yと見なすと、(17)式自身が重みwk、k=1、2、…、pを回帰係数とする重回帰モデルの形式になっている。従って、これらの予測時系列データと同定に用いた出力データとを用いて、重みwk、k=1、2、…、pを重回帰の方法で同定することができる。ただし、重みは各々正でなければならないという不等式制約とw1+w2+…+wp=1という等式制約を満たす必要があるため、通常の重回帰で用いる最小2乗法は適用できないが、混合線形推定法などの制約条件を考慮できる線形回帰の推定法を用いて、重みを推定することが可能になる。
なお、この方法は、(13)式から(14)式を得る代わりに(15)式を得る箇所にも同様の考え方を適用することができるので、単回帰モデル⇒FIR/ARモデル⇒多入力の伝達関数モデル、という3段階の段階的な学習(推定)によって予測モデルを合成することができる。
本実施形態において、上記のようにして重みwk、k=1、2、…、pを推定した(17)式によって予測モデルの合成法を定義する方法が、予測モデル合成法定義部5の動作の他の例である。
四つ目は、一つ目から三つ目とは異なり、合成した予測に傾向(バイアス)を持たせたい場合のアプローチであり、意図的に過大、あるいは、過小の予測を行う事を目的とする場合の予測出力の合成法である。
上記の予測を行いたい状況は、現実の問題ではしばしば遭遇する。例えば、本実施形態で行う雨水の流入量予測の場合、この予測情報は雨水排水ポンプの起動や停止のタイミングを図るための支援情報として利用される場合が多い。このような場合、実際に流入する雨水流入量より過小な流入量の予測を行ってしまうことはリスク回避の観点から致命的になる事がある。
このような場合、もちろん、正確な流入量を予測することが好ましい事は言うまでもないが、それが現実的に困難な場合は、多少多めの量を予測しておく方が安全である。
このような場合、もちろん、正確な流入量を予測することが好ましい事は言うまでもないが、それが現実的に困難な場合は、多少多めの量を予測しておく方が安全である。
また、例えば雨水排水ポンプの制御において流入量予測を利用する際には、ポンプの起動タイミングを図るために過大側の予測をすることが好ましいが、停止タイミングを図るためには過小側の予測をすることが好ましいため、過小側の予測を行いたい場合もある。
また、例えば、水や電力の需要予測などの場合は、あまり過小に予測を行ってしまうと、水や電力の供給に支障をきたす可能性が考えられるため、若干過大側に予測を行っておく方が安全な場合もある。
また、このように、過大側あるいは過小側の予測を行える様にしておくと、純粋に過大な予測や過小の予測を行いたいという場合以外にも、実際に予測システムを運用して長期評価した場合に、予測値が過小に出る傾向がある場合には過大側に再調整したり、逆に予測値が過大に出る場合には過小側に再調整したりして微調整を行うことも可能になる。
このような動機に基づいて行う予測出力の合成法を以下に示す。これを行うためには、各ペアワイズ予測モデルの出力予測値y1、y2、…、ypに対して、先に述べた位置母数の推定方法を定義する以外に、標準偏差や分散に対応する尺度母数の推定方法も同時に定義する。例えば、位置母数を平均とした場合には、尺度母数として、標準偏差を採用する様に定義する。また、位置母数としてロバスト推定法のメジアン(中央値)を用いた場合には、尺度母数としてMAD(中央値絶対偏差:Median Absolute Deviation)を採用する。さらに、HL推定、MAD推定、ブートストラップ推定、M推定などを用いた場合には、各々の位置母数に対する尺度母数の推定方法を定義する。このようにして、y1、y2、…、ypに対して適用する位置母数と尺度母数の推定法を予め定義しておく。これを以下では各々μとσと表す(通常μは平均値、σは標準偏差の意味で使われることが多いが、ここでは、必ずしも平均と標準偏差ではなく、ここで定義した位置母数と尺度母数を表すものとする)。
次に、各ペアワイズ予測モデルの予測出力の集合、y1、y2、…、ypに対して、その集合のμ±kσの値を合成した予測出力とするように、kの値と符号とを指定する。例えば、k=1としてμ+σを指定すると、y1、y2、…、ypの中の位置母数から尺度母数の1倍だけ過大な値を合成した予測出力とすることを意味する。逆にk=-1とすると、y1、y2、…、ypの中の位置母数から尺度母数の1倍だけ過小な値を合成した予測出力とすることを意味する。もし、位置母数と尺度母数として平均と標準偏差を用いた場合は、y1、y2、…、ypが正規分布に従うと仮定すると、μ+σは約68%の値(μが50%の値)、μ-σは約32%の値を取り出すことに対応する。
なお、平均や標準偏差を用いなくても、y1、y2、…、ypの中のK%にあたるK%分位点を取り出して合成出力とするという方法も、位置母数と尺度母数を適切に定義すれば、μ±kσの形式で書けるので、このように設定した分位点と直接抽出するという方法で指定することも可能である。例えば、過大側の予測を行いたい場合に予測の異常値をある程度除外する事を想定して90%分位点にあたる予測値を採用し過大予測を行う様に合成した予測出力を定義したり、10%分位点にあたる予測値を採用し過小予測を行う様に合成した予測出力を定義したりすることもできる。
予測モデル合成法定義部5の定義方法として四つの方法を説明したが、予測モデル合成法定義部5は、これらのいずれか、あるいは、その組み合わせによって予想モデルの合成法と定義することができる。なお、以上の一連の作用は、オフラインで過去のデータを用いて定期的もしくは非定期的に実行される。
次に、オフラインで同定および定義した定義式を用いて、オンラインで予測を行う。
オンラインの予測は、時間の進行方向における所定の周期TH(<<TL(オフライン同定を定期的に行う場合の周期))で行われる。
オンラインの予測は、時間の進行方向における所定の周期TH(<<TL(オフライン同定を定期的に行う場合の周期))で行われる。
まず、オンライン予測用データ抽出部32から周期THで、入力変数のデータを抽出する。本実施形態では、入力変数は、流量計(雨水ポンプ井流入量計)11による雨水ポンプ井流入量と、幹線流量計12による幹線流入量と、幹線水位計131~13KによるK個の幹線水位計と、地上雨量計141~14MによるM個の地上雨量と、Q×Pメッシュ1511~15QPの各メッシュにおけるレーダ雨量と、を含む。
次に、ペアワイズ出力変数予測部6では、ペアワイズ予測モデル同定部4で同定した各入力変数と出力変数のペアワイズ予測モデルを用いて、p個の予測出力を計算する。これは、(5)式、(14)式、(15)式などの形で同定したパラメータ値が同定されたペアワイズ予測モデルにオンライン予測用データ抽出部32で抽出したオンラインの入力変数データを入力することで直ちに計算できる。これがペアワイズ出力変数予測部6の作用である。
次に、合成出力変数予測部7では、予測モデル合成法定義部5で定義した(17)式や(18)式に相当する合成法の定義式を用いて、ペアワイズ出力変数予測部6から出力される各入力変数に対するペアワイズの予測出力を入力することで合成した予測出力が計算される。
次に、予測誤差評価部8は、合成出力変数予測部7で計算された予測出力の予測結果とペアワイズ出力変数予測部6のペアワイズの予測出力結果とを取得し、時間の進行方向における所定の周期THで時系列データとして保存する。ここで、保存された予測結果は、誤差評価を行うタイミングまで保持される。予測誤差の評価を行うタイミングは、予測モデルのユーザなどが外部から指示しても良いし、時間の進行方向における所定の周期TM(例えばTH<TM<TL)で実施しても良い。いずれの方法であっても、予測誤差の評価を行うタイミングで、評価用データ抽出部33で、予測誤差評価部8に保存された予測出力の時系列データと対応する時間の計測された出力変数、すなわち、本実施形態では雨水流入量の時系列データを抽出する。
そして、予測誤差評価部8では、抽出された出力(雨水流入量)の時系列データとp個のペアワイズ予測モデルの予測出力との誤差、および、合成した予測出力の誤差を評価する。この際、誤差の評価は、通常時系列解析で行われるMSE(平方2乗誤差=L2誤差)やその平方根であるRMSEなどであっても良いが、目的に応じて、誤差評価の評価基準を適宜設定しておいても良い。例えば、ピークの雨水流入量の予測精度が重要である場合はピーク流入量の差や誤差が最大になる場合の差(L∞誤差=∞ノルム誤差)を誤差評価基準としても良いし、Nash-Sutcliffe係数などの評価基準を用いても良い。また、流入量が急激に増加する場合の遅れを改善する事を目的としたい場合には、予測出力と実績出力との位相差を計算し、位相差を小さくすることを誤差の評価基準としても良い。いずれにしても、予測誤差評価部8は、設定した誤差評価基準に基づいて、所定のタイミングもしくは周期TMで、時系列データとp個のペアワイズ予測モデルの予測出力との誤差、および、合成した予測出力との誤差を評価する。
次に、ペアワイズ予測モデル修正部9では、予測誤差評価部8で評価した予測誤差に基づいて、ペアワイズ予測モデルのパラメータの調整、あるいは、特定のペアワイズモデルの削除、もしくは、複数のペアワイズモデルの統合を行う。
以下に、上記のような調整機能、削除機能、および、統合機能を付加したモジュラー型予測をより具体化した方法の例について説明する。
一つ目のパラメータ調整について、ペアワイズ予測モデル修正部9は、予測誤差評価部8で評価した予測誤差に基づいて、著しく予測精度が悪い(予測誤差が所定の許容誤差より大きい)と判断される入出力ペアを抽出し、抽出された入出力ペアのペアワイズ予測モデルの再同定を行う。予測精度の悪さの判定は、(1)予測誤差に対するしきい値(許容誤差)を直接指定してしきい値を超過したペアワイズ予測モデルを予測精度の悪いモデル(あるいは入出力ペア)と判断する、(2)合成した予測出力の誤差とペアワイズ予測モデルの各予測誤差とを比較し、合成した予測誤差に対して相対的に著しく誤差が大きい(予め設定された許容範囲を超える)予測誤差を持つペアワイズ予測モデルを予測精度の悪いモデルと判断する、(3)ペアワイズ予測モデルの予測誤差(+合成した予測誤差も含んでも良い)に対する位置母数μeと尺度母数σeを計算し、μe±kσe(kは設定するパラメータで2~3ぐらいで通常は設定する)から外れた予測誤差を持つペアワイズ予測モデルを予測精度の悪いモデルと判断する、などの方法で行うことができる。(1)-(3)のような判断によって、パラメータの再調整(再同定)が必要と判断された場合、ペアワイズ予測モデル修正部9は、所定の期間の直近のデータを用いてパラメータの再同定を行う。
一つ目のパラメータ調整について、ペアワイズ予測モデル修正部9は、予測誤差評価部8で評価した予測誤差に基づいて、著しく予測精度が悪い(予測誤差が所定の許容誤差より大きい)と判断される入出力ペアを抽出し、抽出された入出力ペアのペアワイズ予測モデルの再同定を行う。予測精度の悪さの判定は、(1)予測誤差に対するしきい値(許容誤差)を直接指定してしきい値を超過したペアワイズ予測モデルを予測精度の悪いモデル(あるいは入出力ペア)と判断する、(2)合成した予測出力の誤差とペアワイズ予測モデルの各予測誤差とを比較し、合成した予測誤差に対して相対的に著しく誤差が大きい(予め設定された許容範囲を超える)予測誤差を持つペアワイズ予測モデルを予測精度の悪いモデルと判断する、(3)ペアワイズ予測モデルの予測誤差(+合成した予測誤差も含んでも良い)に対する位置母数μeと尺度母数σeを計算し、μe±kσe(kは設定するパラメータで2~3ぐらいで通常は設定する)から外れた予測誤差を持つペアワイズ予測モデルを予測精度の悪いモデルと判断する、などの方法で行うことができる。(1)-(3)のような判断によって、パラメータの再調整(再同定)が必要と判断された場合、ペアワイズ予測モデル修正部9は、所定の期間の直近のデータを用いてパラメータの再同定を行う。
本実施形態において、例えば、M個の地上雨量データの中のいずれかの地上雨量データと雨水流入量との間のペアワイズ予測モデルの精度が悪化したと判断された場合、ペアワイズ予測モデル修正部9は、地上雨量と雨水流入量との間のペアワイズ予測モデルのパラメータのみを再同定する。このように部分的な同定による再調整を行う事ができる点は、モジュラー型の構成(構成的なアプローチ)の大きな利点である。
二つ目の特定のペアワイズモデルの削除は、一つ目のパラメータの調整と類似の手続きで行うことができる。基本的な手順としては、ペアワイズ予測モデル修正部9は、一つ目のパラメータの調整で述べた方法と類似の方法で誤差評価を行い、誤差が著しく大きくなり精度が劣化したペアワイズ予測モデルを抽出し、抽出したペアワイズモデルを削除する。誤差評価の方法は一つ目の方法と同じであるが、ペアワイズ予測モデル修正部9は、(1)劣化の程度を判断するしきい値を一つ目のものよりも大きくして、一つ目の調整よりも明らかに予測精度が劣化しているものを削除対象とする、(2)まず一つ目のパラメータの再調整を行い、パラメータの再調整を行っても精度劣化の判断に用いた基準をクリアできない場合に削除対象とする、などの方法で削除すべきペアワイズ予測モデルを抽出する。
本実施形態においては、例えば、幹線水位計などは水位計が水没することなどにより、水位を正確に測ることができなくなりアウトライア(異常値)が多量に混入してしまうことなどが考えられる。この場合、水没した水位計によるペアワイズ予測モデルは分離されることになる。このように、予測精度に悪影響を与えている部分を部分的に切り離すことで予測精度の改善を図れることもモジュラー型の構成(構成的なアプローチ)の大きな利点である。
3つ目のペアワイズ予測モデルの統合は、以下の様に実施される。まず、ペアワイズ予測モデル修正部9は、ペアワイズ予測モデルの予測誤差を相互に比較評価し、その予測誤差同士の差の絶対値が所定のしきい値以下となる相互に類似の予測精度を持つペアワイズ予測モデルを、相互に類似する予測モデルとして抽出しグループ化する。
次に、ペアワイズ予測モデル修正部9は、このようにしてグループ化されたペアワイズ予測モデルの持つパラメータの値をグループ内で相互に比較する。例えば、係数パラメータak、k=1、2、…、nをベクトル化してA=[a1、a2、…、an]とし、相互に比較する予測モデルの係数パラメータベクトル同士の類似度を例えば内積などにより算出し、所定のしきい値により類似していると判断されたペアワイズ予測モデルを統合の対象候補とする。
本実施形態のモジュラー型の構成的なアプロ―チでは、パラメータの可同定性(一意決定可能性)が基本的に担保されているので、予測結果が類似しておりパラメータ値も類似していれば、その入力変数はほぼ同じ値を持つはずであるが、念のため、ペアワイズ予測モデル修正部9は、候補として抽出されたペアワイズ予測モデルの入力変数同士の相関係数を求めて、所定のしきい値(例えば、相関係数0.95)以上である場合に最終的に統合すべきペアワイズ予測モデルと判断する。そして、ペアワイズ予測モデル修正部9は、最終的に統合すべきと判断された複数のペアワイズ予測モデルの入力変数の平均をとった変数を、新たな一つの合成された入力変数として定義する。
そして、ペアワイズ予測モデル修正部9は、(1)先に比較した係数ベクトルの平均値を採用する、あるいは、(2)新たに合成した平均化された入力変数と出力変数に対して再同定を行って推定する、のいずれかの方法を用いて、合成された入力変数に対するペアワイズ予測モデルの係数を決定する。
上記のような統合は、本実施形態のモジュラー型時系列データ予測装置においては、例えば次の様なケースにおいて生じると考えられる。入力変数の中にQ×Pメッシュ1511~15QPの各々のレーダ雨量のデータが含まれているケースにおいて、レーダ雨量のメッシュのサイズは、例えばXRAINと呼ばれるXバンドの気象レーダの場合約250m×250mであるため、隣接するメッシュの雨量データはほとんど同じ値をとる場合が多い。このような場合、相互に隣接するメッシュの雨量データは統合の対象となる可能性が高く、Q×P個の入力変数が、少数の互いに独立した合成入力変数として統合されることになる可能性が高い。このような統合が起こると、レーダ雨量のメッシュ状のデータが、雨水流入量に与える影響の類似性に応じて統合されるため、入力変数の数を減らせるだけでなく、説明性が向上する事が期待できる。
また、このような場合は、先に述べた多重共線性の問題が生じているケースであるが、これに対し、例えば、Lassoと呼ばれる方法などを用いて入力変数を自動的に選択すると、いくつかのメッシュの雨量データが代表的な入力変数として選択され、選択された入力変数に対応するメッシュの雨量と極めて類似するメッシュの雨量データは無視されることになる。雨量データが極めて類似しているので、このように代表メッシュの雨量を入力変数として選択しても通常は問題がない。しかしながら、近年頻発している局所的な豪雨が見られる場合、例えば、たまたま選択したメッシュの雨量はあまり多くないのに、入力変数として選択されなかった近くのメッシュの雨量が非常に大きい様な局所的な降雨が生じた場合、あるメッシュを代表メッシュとして選択してしまうと、流入量が増加しないという誤った予測を行う可能性がある。
一方、ここで述べた入力変数の統合を行うと、類似の影響を与えるメッシュの雨量が平均化され、これは統合されたメッシュの平均雨量という意味を持つため、このようなケースでも妥当な予測結果を与える事が期待できると同時に、合理的な説明性も向上する。このように統合すべき入力変数をまとめて解釈し、説明しやすくするようにできる事もモジュラー型の構成(構成的なアプローチ)の大きな利点である。
最後に、合成出力変数予測部7で算出(合成)された予測出力は、出力予測結果観測部10に送られ、例えばいわゆるトレンドグラフと呼ばれる時系列データとしてユーザであるプラントの管理者や運転員に提示される。この際、リアルタタイム監視している監視の現在までの値は、実際の出力値である雨水流入量とその予測値を同時に表示し、将来の予測値は予測出力のみを表示する様にしておくことが好ましい。
また、合成した予測出力(予測雨水流入量)だけでなく、ペアワイズ予測モデルのp個の予測出力を用いて、その位置母数μyと尺度母数σyを計算し、μy±Kσy(Kは設定パラメータで2~3程度で設定する)の範囲を同時表示したり、あるいは、ペアワイズ予測モデルの予測出力の中から最大値(あるいはロバスト性を考えて最大値に近い95%点にもっとも近い予測値)や最小値(あるいはロバスト性を考えて最小値に近い5%点にもっとも近い予測値)を直接予測出力の範囲として同時に表示したりしても良い。
本実施形態のモジュラー型時系列データ予測装置によれば、一般にはブラックボックスモデルと呼ばれるモデルに対して、その内部パラメータに対して物理法則に矛盾しない合理的な解釈が可能になる様な仕組みを導入することで、ホワイトボックスモデルに近い、モデルの解釈と合理的な説明性を向上させて、部分的な調整、削除、統合などのモデルの維持管理持管理(メンテナンス)の煩雑さを大幅に低減できるという効果が得られる点である。
なお、このような効果を得るための根本的なアイデアは、本来、入出力データしか与えられていない内部構造が完全に未知である状況下で、出力に関連する可能性のある要素を要素毎に出力と直接関係づけ(ペアワイズ予測モデルに対応)、各要素と出力の関係を、最終的な入出力関係の内部構造として定義される様に、構成的に予測モデルを構築していることである。このようなアイデアを用いることで、本実施形態のモジュラー型時系列データ予測装置では、本来知ることのできない内部構造に立ち入ることなく、疑似的にホワイトボックスモデルと同様の説明性と調整の容易さを実現している。
また、このモジュラー型時系列予測モデルの構築方法は、統計的機械学習分野でよく知られている。先に述べたバギングやブースティングなどのアンサンブル学習と呼ばれる手法とも類似している。バギングやブースティングは、弱学習器(弱学習モデル)と呼ばれる、簡単に構成できるが、精度が必ずしも高くないモデルをたくさん集めることで、精度の高い強学習器(強学習モデル)と呼ばれる予測モデルを構築する方法であり、バギングは弱学習器に対して多数決などの処理で強学習を行う方法、ブースティングは、弱学習器の結果を用いて強学習器を学習させる方法である。ただし、バギングやブースティングでは、学習(同定)に利用データを変えることでたくさんの弱学習器を生成しているが、本実施形態では、同定データを変えるのではなく、説明変数を変えることで、弱学習器に相当するペアワイズ予測モデルを同定している点で異なっている。なお、本実施形態のモジュラー型時系列データ予測装置において、合成出力の予測方法を定義する際に、ロバスト位置母数推定や予め指定した重みを用いた重み付き平均を用いるような方法はバギングの考え方に近く、重み付き平均の重みを再同定によって求める方法はブースティングの考え方に近い。
以下、本実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの複数の実施例について説明する。
(第1実施例)
図5は、一実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの第1実施例について説明するための図である。
図5では、(1)式のFIRモデルやARモデルをペアワイズ予測モデルとした例を示している。
(第1実施例)
図5は、一実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの第1実施例について説明するための図である。
図5では、(1)式のFIRモデルやARモデルをペアワイズ予測モデルとした例を示している。
本実施例のペアワイズ予測モデルMAは、回帰ブロックA1-ANと、予測出力の合成ブロックA5と、を備えている。
例えば回帰ブロックA1-A3のそれぞれは、入力変数に対して重み付き移動平均成分(FIR成分/MA成分)を予測値として出力する。例えば回帰ブロックANは、入力変数(過去の出力変数)に対して自己回帰成分(AR成分)を予測値として出力する。
例えば回帰ブロックA1-A3のそれぞれは、入力変数に対して重み付き移動平均成分(FIR成分/MA成分)を予測値として出力する。例えば回帰ブロックANは、入力変数(過去の出力変数)に対して自己回帰成分(AR成分)を予測値として出力する。
予測出力の合成ブロックA5は、回帰ブロックA1-ANから出力された予測値を取得し、予測値を合成した合成予測値(Y(t+1)=F(Y1(t+1),Y2(t+1),…,Yn(t+1)))を演算して出力する。
例えば、通常のARXモデルなどによってパラメータを同定すると自己回帰成分(AR成分)の重みが大きくなり、FIR成分に相当する入力変数の値が急激に変化した場合に、予測出力が計測出力に追従できなくなり予測に位相遅れが生じる事がある。しかし、本実施例の構成をとることにより、AR成分の影響を調整することが可能になり、位相遅れを改善して計測出力の値が変化する前に出力の変化を予測することが可能になる。
例えば、通常のARXモデルなどによってパラメータを同定すると自己回帰成分(AR成分)の重みが大きくなり、FIR成分に相当する入力変数の値が急激に変化した場合に、予測出力が計測出力に追従できなくなり予測に位相遅れが生じる事がある。しかし、本実施例の構成をとることにより、AR成分の影響を調整することが可能になり、位相遅れを改善して計測出力の値が変化する前に出力の変化を予測することが可能になる。
(第2実施例)
図6は、一実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの第2実施例について説明するための図である。
図6は、合成した予測に傾向(バイアス)を持たせたい場合に用いられる予測モデルの一例として、意図的に過大の予測および過小の予測を行う場合の予測モデルを示している。
図6は、一実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの第2実施例について説明するための図である。
図6は、合成した予測に傾向(バイアス)を持たせたい場合に用いられる予測モデルの一例として、意図的に過大の予測および過小の予測を行う場合の予測モデルを示している。
本実施例のペアワイズ予測モデルMBは、回帰ブロックB1-BNと、予測出力の合成ブロックB5と、を備え、合成ブロックB5が、過大推定予測値(Y+(t+1)=F+(Y1(t+1),Y2(t+1),…,Yn(t+1)))と、通常推定予測値(Y(t+1)=F(Y1(t+1),Y2(t+1),…,Yn(t+1)))と、過小推定予測値(Y-(t+1)=F-(Y1(t+1),Y2(t+1),…,Yn(t+1)))と、を出力する点において上述の第1実施例のペアワイズ予測モデルMAと異なっている。
本実施例では、ペアワイズ予測モデルの位置母数と尺度母数とを用いて、通常の予測に加えて過大な予測および過小な予測を行うことにより、安全側の予測(リスクを回避する側の予測)などのニーズ(要求)に応じて、予測結果を微調整することができる。
図7および図8は、第2実施例のペアワイズ予測モデルの効果を説明するための図である。
図7は、回帰手法の一種であるL1正則化を導入したLasso回帰による雨水流入予測値と測定値との一例を示している。
図8は、第2実施例の予測モデルによる雨水流入予測値と測定値との一例を示している。
図7は、回帰手法の一種であるL1正則化を導入したLasso回帰による雨水流入予測値と測定値との一例を示している。
図8は、第2実施例の予測モデルによる雨水流入予測値と測定値との一例を示している。
図7に示す雨水流入予測値と測定値とを比較すると、予測誤差は小さいが、若干推定値が測定値よりも小さい結果となっており、過小推定になっている事がわかる。
一方、図8に示す本実施例の手法で得られた通常の推定(平均的な推定)、過大推定、および、過少推定を行った結果と測定値とを比較すると、通常の推定値(予測流入量)では予想誤差は小さいく、ピーク付近での予測精度は良いものの流入量が増加するタイミングでの予測値が測定値を下回っている。
一方、図8に示す本実施例の手法で得られた通常の推定(平均的な推定)、過大推定、および、過少推定を行った結果と測定値とを比較すると、通常の推定値(予測流入量)では予想誤差は小さいく、ピーク付近での予測精度は良いものの流入量が増加するタイミングでの予測値が測定値を下回っている。
このような場合、この予測値を制御に用いる場合には、もう少し流入量急増時の予測を過大に評価しておきたい場合がある。このとき、過大推定予測を行うと、全体的にMSE(平均2乗誤差)での予測精度は劣化するが、流入量急増時に急増する可能性を示唆する予測が可能になることがわかる。また、通常の予測値よりも過小に評価をしておきたい場合には、過小推定予測を行うと、予測値が小さくなるように調整することができる。
このような過大推定および過小推定により予測値を調整可能とすると、長期の運用の結果、通常の予測値が実測値に対して過小推定、あるいは、過大推定になっている様な傾向(バイアス)が見られる場合、通常の予測値を過大推定値、あるいは過小推定値に置き換えることで、バイアスを補正してより精度の高い予測を行うことができる。
(第3実施例)
図9は、一実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの第3実施例について説明するための図である。
本実施例では、ペアワイズ予測モデルとして、(5)式のようなARモデル/FIRモデルをさらに分解し、時間遅れを考慮した変数を各々一つの説明変数と見なした単回帰モデルに分解して予測モデルを構築している。
図9は、一実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの第3実施例について説明するための図である。
本実施例では、ペアワイズ予測モデルとして、(5)式のようなARモデル/FIRモデルをさらに分解し、時間遅れを考慮した変数を各々一つの説明変数と見なした単回帰モデルに分解して予測モデルを構築している。
本実施例のペアワイズ予測モデルMCは、複数の時間遅れブロックと、複数の単回帰ブロックと、相関重み付き平均ブロックと、を含む予測出力ブロックC1-CNを入力変数1-N毎に含み、当該複数の予測出力ブロックC1-CNから出力された予測値を合成した最終合成予測値を出力する予測出力値の合成ブロックC(N+1)を更に備える。
予測出力ブロックC1は、複数(m)の時間遅れブロック111-11mと、複数(m)の単回帰ブロック121-12mと、相関重み付き平均ブロック13と、を含む。
複数の時間遅れブロック111-11mは、入力変数1(例えば時刻tにおける値)に対してそれぞれに割り当てられた遅れ時間L11-L1m1だけ遅れた変数(例えば時刻t-L11、…、t-L1m1における値)を出力する。
複数の時間遅れブロック111-11mは、入力変数1(例えば時刻tにおける値)に対してそれぞれに割り当てられた遅れ時間L11-L1m1だけ遅れた変数(例えば時刻t-L11、…、t-L1m1における値)を出力する。
複数の単回帰ブロック121-12mは、単回帰係数を用いて、時間遅れブロック111-11mから出力された変数に基づく予測値を演算する出力する。
相関重み付き平均ブロックは、複数の単回帰ブロック121-12mから出力された複数の予測値の相関重み付き平均値を演算して出力する。
相関重み付き平均ブロックは、複数の単回帰ブロック121-12mから出力された複数の予測値の相関重み付き平均値を演算して出力する。
予測出力ブロックC2-CNは、上記予測出力ブロックC1と同様の構成であるため、個々の説明は省略する。
予測出力値の合成ブロックC(N+1)は、予測出力ブロックC1-CNから出力された予測値を合成(例えば、相関重み付き平均)して、最終推定予測値(Y(t+1)=F(Y1(t+1),Y2(t+1),…,Yn(t+1)))を出力する。
予測出力値の合成ブロックC(N+1)は、予測出力ブロックC1-CNから出力された予測値を合成(例えば、相関重み付き平均)して、最終推定予測値(Y(t+1)=F(Y1(t+1),Y2(t+1),…,Yn(t+1)))を出力する。
本実施例のペアワイズ予測モデルは、遅れ、ゲイン、バイアスの3項でモデルを説明することが可能であり、パラメータの微調整が容易である。また、単回帰係数と相関係数とには陽な関係があるため、各入力変数と各出力変数との相関係数を求めるだけで全体の予測モデルの構築ができ、予測モデルの解釈が容易である。
更に、本実施例のペアワイズ予測モデルでは、第1実施例のペアワイズ予測モデルと同様の効果を得ることができる。
更に、本実施例のペアワイズ予測モデルでは、第1実施例のペアワイズ予測モデルと同様の効果を得ることができる。
図10および図11は、第3実施例のペアワイズ予測モデルの効果の一例を説明するための図である。
図10および図11では、第3実施例のペアワイズ予測モデルに対し10個の入力変数を入力して雨水流入予測を行った例である。ここでは、各入力変数に対する遅れ時間を一つだけ考えた最も単純な場合の予測結果の例を示している。
図10および図11では、第3実施例のペアワイズ予測モデルに対し10個の入力変数を入力して雨水流入予測を行った例である。ここでは、各入力変数に対する遅れ時間を一つだけ考えた最も単純な場合の予測結果の例を示している。
入力変数の一つである降雨強度(地上雨量)の遅れ時間を推定すると33分という推定結果が得られた。図11は、上記推定結果を用いた雨水流入量の予測結果の一例を示している。図11の予測結果によれば、全体的には、ある程度の精度で予測ができているが、流入量が急増する部分での予測の遅れ(位相遅れ)が大きい事がわかる。
この時、位相遅れが大きい理由として、遅れ時間33分というものが大きすぎる可能性が推測できる。この推測に基づき、図10に、遅れ時間を15分と調整した場合の雨水流入量の予測結果の一例を示す。図10に示す予測結果によれば、遅れ時間を調整することにより、全体的な予測精度には大きな影響を与えることなく、流入量急増時の位相遅れを大幅に改善できていることがわかる。
このように、本実施例のペアワイズ予測モデルによれば、容易に解釈可能なパラメータのみで予測モデルを構築可能であり、説明性が向上するだけでなく予測モデルの微調整を容易に行うことができる。
また、単回帰係数と相関係数とは、入出力データを各々平均と標準偏差とで正規化すれば、一致することが知られており、これを逆に非正規化すれば、回帰係数とバイアスとに換算することができる。したがって、各入出力データの平均、標準偏差、および、入出力の相関のみを計算すれば、遅れ時間を指定するだけで、全ての予測モデルのパラメータを同定することができる。
さらに、合成出力の計算を相関係数の重み付き平均で計算する方法とすることにより、予測モデルを構築するために必要な計算は、遅れ時間の推定(指定)と、入出力変数の平均および標準偏差と、入出力変数の相関と、のみであり、予測モデルに含まれるパラメータの全てを、この4種類の推定値(遅れ時間、平均、標準偏差、相関係数)だけを用いて説明することが可能になる。
(効果)
上記のように、本実施形態によれば、入力変数と出力変数との相関関係が変化したり劣化したり、入力変数の時系列データに多量のアウトライアが含まれたりするなどの理由で予測精度が十分でない場合に、全体の予測モデルを再同定することなく、合成した予測出力値に影響を与えるペアワイズの予測モデルの部分的な調整や、悪影響を与える部分の削除を容易に行う事ができると同時に、これにより、予測精度劣化の理由を合理的に説明することが可能になる。
上記のように、本実施形態によれば、入力変数と出力変数との相関関係が変化したり劣化したり、入力変数の時系列データに多量のアウトライアが含まれたりするなどの理由で予測精度が十分でない場合に、全体の予測モデルを再同定することなく、合成した予測出力値に影響を与えるペアワイズの予測モデルの部分的な調整や、悪影響を与える部分の削除を容易に行う事ができると同時に、これにより、予測精度劣化の理由を合理的に説明することが可能になる。
また、入力変数間の多重共線性が極めて強くほぼ同じ説明能力を持つ入力変数が含まれる場合にも、パラメータの可同定(一意決定可能性)を維持することができ、物理法則に矛盾しない合理的な説明が可能な予測モデルを構築できると同時に、ほぼ同じ説明能力を持つ入力変数を統合することで、予測精度を維持しながら、説明性を向上させることができる。
すなわち、本実施形態によれば、予測結果を合理的に説明するとともに、容易に予測結果を調整可能とするモジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラムを提供することができる。
すなわち、本実施形態によれば、予測結果を合理的に説明するとともに、容易に予測結果を調整可能とするモジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラムを提供することができる。
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。
1…都市雨水排水プロセス、11…流入量計、12…幹線流量計、131-13K…幹線水位計、141~14M…地上雨量計、15、1511~15QP…レーダ雨量計、17…流入渠、110…流入ゲート、18…雨水ポンプ井、19…雨水ポンプ、161、162…雨水ポンプ井水位計、2…データ収集保存部、21…出力変数データ選択部、22…入力変数データ選択部、3…データ抽出部、31…オフライン予測モデル同定用データ抽出部、32…オンライン予測用データ抽出部、33…評価用データ抽出部、4…ペアワイズ予測モデル同定部、5…予測モデル合成法定義部、6…ペアワイズ出力変数予測部、7…合成出力変数予測部、8…予測誤差評価部、9…ペアワイズ予測モデル修正部、10…出力予測結果観測部、111-11m…時間遅れブロック、121-12m…単回帰ブロック、13…相関重み付き平均ブロック、A1-AN…回帰ブロック、A5…合成ブロック、B1-BN…回帰ブロック、B5…合成ブロック、C1-CN…予測出力ブロック、C(N+1)…合成ブロック。
Claims (12)
- 複数のプロセス変数を所定の周期で計測する複数のプロセスセンサを有するシステム又はプロセスに適用される装置であって、
複数の前記プロセス変数の時系列データを所定の周期で収集し保存するとともに、複数の前記プロセス変数の中から予測対象となる少なくとも一つの出力変数を選択する出力変数データ選択部と、複数の前記プロセス変数の中から複数の入力変数の候補を選択する入力変数データ選択部と、を含むデータ収集保存部と、
前記時系列データから抽出された前記出力変数と複数の前記入力変数との同定用データを用いて、1入力1出力のペア毎にペアワイズ予測モデルのパラメータを同定して複数の前記ペアワイズ予測モデルを定義するペアワイズ予測モデル同定部と、
複数の前記ペアワイズ予測モデルから出力されるペアワイズ予測値の合成法を定義する予測モデル合成法定義部と、
時間の進行方向に所定の周期又はリアルタイムで前記時系列データから抽出された複数の前記入力変数の予測用データを、複数の前記ペアワイズ予測モデルに入力して、複数の前記入力変数のそれぞれに対応する前記ペアワイズ予測値を演算するペアワイズ出力変数予測部と、
前記合成法により複数の前記ペアワイズ予測値を合成して前記出力変数の予測値を演算する合成出力変数予測部と、
を有するモジュラー型時系列データ予測装置。 - 前記ペアワイズ予測モデル同定部は、前記出力変数以外の複数の前記入力変数の同定用データに有限インパルス応答モデルを適用し、前記出力変数に自己回帰モデルを適用する、請求項1記載のモジュラー型時系列データ予測装置。
- 前記ペアワイズ予測モデル同定部は、複数の前記入力変数の各々ついて1又は複数の遅れ時間を組み込んだ単回帰モデルを用いて、1入力1出力の伝達関数モデルを構築し、前記単回帰モデルを合成したモデルを前記ペアワイズ予測モデルとする、請求項1記載のモジュラー型時系列データ予測装置。
- 前記ペアワイズ予測モデル同定部は、複数の前記入力変数各々の同定用データに対して非線形変換を行い、前記ペアワイズ予測モデルを構築する、請求項1記載のモジュラー型時系列データ予測装置。
- 前記予測モデル合成法定義部は、複数の前記ペアワイズ予測値の重み付き平均を演算することを前記合成法とし、前記合成法における重みを、前記ペアワイズ予測モデル同定部によって同定された複数の前記ペアワイズ予測モデル各々の予測精度を表す指標に基づいて決定する、請求項1記載のモジュラー型時系列データ予測装置。
- 前記予測モデル合成法定義部は、複数の前記ペアワイズ予測値の重み付き平均を演算することを前記合成法とし、前記ペアワイズ予測モデル同定部で用いた同定用データと同じデータを用いて、前記合成法における重みの値を同定して決定する、請求項1記載のモジュラー型時系列データ予測装置。
- 前記ペアワイズ出力変数予測部から出力される複数の前記ペアワイズ予測値と、前記合成出力変数予測部から出力される前記出力変数の予測値とに対して、予め指定した所定の周期毎あるいは所定期間で、前記時系列データから抽出した実績値との誤差を評価する予測誤差評価部と、
複数の前記ペアワイズ予測値の誤差と、前記出力変数の予測値の誤差とを比較することにより、調整すべき前記ペアワイズ予測モデル、統合すべき前記ペアワイズ予測モデル、および、分離すべき前記ペアワイズ予測モデルを決定して、前記ペアワイズ出力変数予測部で用いられる複数の前記ペアワイズ予測モデルを修正するペアワイズ予測モデル修正部と、
を更に有する請求項1記載のモジュラー型時系列データ予測装置。 - 前記ペアワイズ予測モデル修正部は、前記予測誤差評価部によって評価された複数の誤差の少なくとも一つが所定の許容誤差を超過したときに、当該誤差に対応する前記入力変数の候補に関する前記ペアワイズ予測モデルの再同定を行い、前記ペアワイズ出力変数予測部で用いられる複数の前記ペアワイズ予測モデルを修正する請求項7記載のモジュラー型時系列データ予測装置。
- 前記ペアワイズ予測モデル修正部は、前記予測誤差評価部によって評価された複数の前記誤差の少なくとも一つが所定の許容誤差を超過したときに、当該誤差に対応する前記入力変数に関する前記ペアワイズ予測モデルを削除して、前記ペアワイズ出力変数予測部で用いられる複数の前記ペアワイズ予測モデルを修正する請求項7記載のモジュラー型時系列データ予測装置。
- 前記ペアワイズ予測モデル修正部は、前記予測誤差評価部によって評価された複数の前記誤差を比較し、複数の前記誤差の差が所定の値以下であるときに、当該誤差に対応する前記入力変数の和あるいは平均値を一つの新たな入力変数として定義し、前記ペアワイズ出力変数予測部で当該新たな入力変数に関する前記ペアワイズ予測モデルが用いられるように、複数の前記ペアワイズ予測モデルを修正する請求項7記載のモジュラー型時系列データ予測装置。
- 複数のプロセス変数を所定の周期で計測する複数のプロセスセンサを有するシステム又はプロセスに適用される方法であって、
複数の前記プロセス変数の時系列データを所定の周期で収集して保存し、
複数の前記プロセス変数の中から予測対象となる少なくとも一つの出力変数を選択し、
複数の前記プロセス変数の中から複数の入力変数の候補を選択し、
前記時系列データから抽出された前記出力変数と複数の前記入力変数との同定用データを用いて、1入力1出力のペア毎にペアワイズ予測モデルのパラメータを同定して複数の前記ペアワイズ予測モデルを定義し、
複数の前記ペアワイズ予測モデルから出力されるペアワイズ予測値の合成法を定義し、
時間の進行方向に所定の周期又はリアルタイムで前記時系列データから抽出された複数の前記入力変数の予測用データを、複数の前記ペアワイズ予測モデルに入力して、複数の前記入力変数のそれぞれに対応する前記ペアワイズ予測値を演算し、
前記合成法により複数の前記ペアワイズ予測値を合成して前記出力変数の予測値を演算する、
モジュラー型時系列データ予測方法。 - コンピュータに、請求項11に記載の方法を実行させるモジュラー型時系列データ予測プログラム。
Priority Applications (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2021079844A JP2022173863A (ja) | 2021-05-10 | 2021-05-10 | モジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラム |
PCT/JP2022/018118 WO2022239609A1 (ja) | 2021-05-10 | 2022-04-19 | モジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラム |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2021079844A JP2022173863A (ja) | 2021-05-10 | 2021-05-10 | モジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラム |
Publications (1)
Publication Number | Publication Date |
---|---|
JP2022173863A true JP2022173863A (ja) | 2022-11-22 |
Family
ID=84029588
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2021079844A Pending JP2022173863A (ja) | 2021-05-10 | 2021-05-10 | モジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラム |
Country Status (2)
Country | Link |
---|---|
JP (1) | JP2022173863A (ja) |
WO (1) | WO2022239609A1 (ja) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR102640983B1 (ko) * | 2022-12-21 | 2024-02-23 | 재단법인차세대융합기술연구원 | 극단치와 증감 추세를 반영하여 시계열 데이터를 기호화하는 분석 서버 및 그것의 데이터 분석 방법 |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116090388B (zh) * | 2022-12-21 | 2024-05-17 | 海光信息技术股份有限公司 | 芯片内部电压预测模型生成方法、预测方法及相关装置 |
CN117036112B (zh) * | 2023-10-09 | 2023-12-22 | 石家庄坤垚科技有限公司 | 一种土地规划用的地理信息系统及方法 |
CN118014719A (zh) * | 2024-04-08 | 2024-05-10 | 南京启尚数字科技有限公司 | 一种基于线性回归模型的企业信用智能分析方法和系统 |
Family Cites Families (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP4102131B2 (ja) * | 2002-07-26 | 2008-06-18 | 株式会社東芝 | 予測モデルシステム |
JP2005222444A (ja) * | 2004-02-09 | 2005-08-18 | Toshiba Corp | 統計的予測値演算方法および装置 |
DE102009005845A1 (de) * | 2009-01-21 | 2010-07-22 | Friedrich-Schiller-Universität Jena | Verfahren zur Indentifizierung insbesondere unbekannter Substanzen durch Massenspektrometrie |
WO2014123993A1 (en) * | 2013-02-05 | 2014-08-14 | Yokogawa Corporation Of America | System, method and apparatus for determining properties of product or process streams |
JP6261960B2 (ja) * | 2013-11-13 | 2018-01-17 | 株式会社東芝 | 雨水排水ポンプ制御装置、雨水排水システム、および雨水排水ポンプ制御プログラム |
JP2018116545A (ja) * | 2017-01-19 | 2018-07-26 | オムロン株式会社 | 予測モデル作成装置、生産設備監視システム、及び生産設備監視方法 |
-
2021
- 2021-05-10 JP JP2021079844A patent/JP2022173863A/ja active Pending
-
2022
- 2022-04-19 WO PCT/JP2022/018118 patent/WO2022239609A1/ja active Application Filing
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR102640983B1 (ko) * | 2022-12-21 | 2024-02-23 | 재단법인차세대융합기술연구원 | 극단치와 증감 추세를 반영하여 시계열 데이터를 기호화하는 분석 서버 및 그것의 데이터 분석 방법 |
Also Published As
Publication number | Publication date |
---|---|
WO2022239609A1 (ja) | 2022-11-17 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
WO2022239609A1 (ja) | モジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラム | |
US11238356B2 (en) | Method of predicting streamflow data | |
CN111222698B (zh) | 面向物联网的基于长短时记忆网络的积水水位预测方法 | |
Mehdizadeh et al. | Hybrid artificial intelligence-time series models for monthly streamflow modeling | |
Han et al. | Bayesian flood forecasting methods: A review | |
Jung et al. | Current status and future advances for wind speed and power forecasting | |
Sahoo et al. | Application of support vector regression for modeling low flow time series | |
US20150206427A1 (en) | Prediction of local and network-wide impact of non-recurrent events in transportation networks | |
Kim et al. | Comparative analysis of long short-term memory and storage function model for flood water level forecasting of Bokha stream in NamHan River, Korea | |
KR102446493B1 (ko) | 확률분포에 기반한 가뭄 예측 방법 및 이를 위한 장치 | |
JP2007205001A (ja) | 流量予測装置 | |
CN112561191A (zh) | 预测模型的训练、预测方法、装置、设备、程序和介质 | |
Tsai et al. | Including spatial distribution in a data‐driven rainfall‐runoff model to improve reservoir inflow forecasting in Taiwan | |
US20210088369A1 (en) | Blockage detection using machine learning | |
JP4102131B2 (ja) | 予測モデルシステム | |
Kasiviswanathan et al. | Flood frequency analysis using multi-objective optimization based interval estimation approach | |
Okuda et al. | Non-parametric prediction interval estimate for uncertainty quantification of the prediction of road pavement deterioration | |
Wang et al. | Multi-model integrated error correction for streamflow simulation based on Bayesian model averaging and dynamic system response curve | |
Darvishi Salookolaei et al. | Application of grey system theory in rainfall estimation | |
Pham | Tracking the uncertainty in streamflow prediction through a hydrological forecasting system | |
Turki et al. | A Markova-Chain Approach to Model Vehicles Traffic Behavior | |
JP4423607B2 (ja) | 降雨予測システム及びそれを利用した排水ポンプ運転支援システム | |
Katipoğlu | Prediction of future hydrological droughts with tree-based algorithms | |
Sahoo et al. | Performance comparison of LS-SVM and ELM-based models for precipitation prediction in Barak valley: A case study | |
Mohammed Abdelkader et al. | Condition prediction of concrete bridge decks using Markov chain Monte Carlo-based method |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
RD01 | Notification of change of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7421 Effective date: 20230105 |
|
A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20240314 |