▶ 株式会社東芝の特許一覧 ▶ 東芝電機サービス株式会社の特許一覧
特開2022-173863モジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラム
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2022173863
(43)【公開日】2022-11-22
(54)【発明の名称】モジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラム
(51)【国際特許分類】
G05B 23/02 20060101AFI20221115BHJP
【ＦＩ】
G05B23/02 P
【審査請求】未請求
【請求項の数】12
【出願形態】ＯＬ
(21)【出願番号】P 2021079844
(22)【出願日】2021-05-10
(71)【出願人】
【識別番号】000003078
【氏名又は名称】株式会社東芝
(71)【出願人】
【識別番号】598076591
【氏名又は名称】東芝インフラシステムズ株式会社
(74)【代理人】
【識別番号】100108855
【弁理士】
【氏名又は名称】蔵田 昌俊
(74)【代理人】
【識別番号】100103034
【弁理士】
【氏名又は名称】野河 信久
(74)【代理人】
【識別番号】100075672
【弁理士】
【氏名又は名称】峰 隆司
(74)【代理人】
【識別番号】100153051
【弁理士】
【氏名又は名称】河野 直樹
(74)【代理人】
【識別番号】100162570
【弁理士】
【氏名又は名称】金子 早苗
(72)【発明者】
【氏名】山中 理
(72)【発明者】
【氏名】平岡 由紀夫
【テーマコード（参考）】
3C223
【Ｆターム（参考）】
3C223AA06
3C223FF05
3C223FF22
3C223FF24
(57)【要約】
【課題】 予測結果を合理的に説明するとともに、容易に予測結果を調整可能とするモジュラー型時系列データ予測装置を提供する。
【解決手段】 実施形態による装置は、複数のプロセス変数の時系列データを所定の周期で収集保存し、複数の前記プロセス変数の中から予測対象となる出力変数を選択し、複数のプロセス変数の中から入力変数の候補を選択する機能２と、時系列データから抽出された出力変数と複数の入力変数との同定用データを用いて、１入力１出力のペアワイズ予測モデルを定義する機能４と、複数のペアワイズ予測モデルの予測値の合成法を定義する機能５と、複数のペアワイズ予測モデルに複数の入力変数の予測用データを入力して、複数の入力変数のそれぞれに対応する予測値を演算する機能６と、定義された合成法により複数の予測値を合成することにより出力変数の予測値を演算する機能７と、を有する。
【選択図】図２
【特許請求の範囲】
【請求項1】
複数のプロセス変数を所定の周期で計測する複数のプロセスセンサを有するシステム又はプロセスに適用される装置であって、
複数の前記プロセス変数の時系列データを所定の周期で収集し保存するとともに、複数の前記プロセス変数の中から予測対象となる少なくとも一つの出力変数を選択する出力変数データ選択部と、複数の前記プロセス変数の中から複数の入力変数の候補を選択する入力変数データ選択部と、を含むデータ収集保存部と、
前記時系列データから抽出された前記出力変数と複数の前記入力変数との同定用データを用いて、１入力１出力のペア毎にペアワイズ予測モデルのパラメータを同定して複数の前記ペアワイズ予測モデルを定義するペアワイズ予測モデル同定部と、
複数の前記ペアワイズ予測モデルから出力されるペアワイズ予測値の合成法を定義する予測モデル合成法定義部と、
時間の進行方向に所定の周期又はリアルタイムで前記時系列データから抽出された複数の前記入力変数の予測用データを、複数の前記ペアワイズ予測モデルに入力して、複数の前記入力変数のそれぞれに対応する前記ペアワイズ予測値を演算するペアワイズ出力変数予測部と、
前記合成法により複数の前記ペアワイズ予測値を合成して前記出力変数の予測値を演算する合成出力変数予測部と、
を有するモジュラー型時系列データ予測装置。
複数の前記プロセス変数の時系列データを所定の周期で収集し保存するとともに、複数の前記プロセス変数の中から予測対象となる少なくとも一つの出力変数を選択する出力変数データ選択部と、複数の前記プロセス変数の中から複数の入力変数の候補を選択する入力変数データ選択部と、を含むデータ収集保存部と、
前記時系列データから抽出された前記出力変数と複数の前記入力変数との同定用データを用いて、１入力１出力のペア毎にペアワイズ予測モデルのパラメータを同定して複数の前記ペアワイズ予測モデルを定義するペアワイズ予測モデル同定部と、
複数の前記ペアワイズ予測モデルから出力されるペアワイズ予測値の合成法を定義する予測モデル合成法定義部と、
時間の進行方向に所定の周期又はリアルタイムで前記時系列データから抽出された複数の前記入力変数の予測用データを、複数の前記ペアワイズ予測モデルに入力して、複数の前記入力変数のそれぞれに対応する前記ペアワイズ予測値を演算するペアワイズ出力変数予測部と、
前記合成法により複数の前記ペアワイズ予測値を合成して前記出力変数の予測値を演算する合成出力変数予測部と、
を有するモジュラー型時系列データ予測装置。
【請求項2】
前記ペアワイズ予測モデル同定部は、前記出力変数以外の複数の前記入力変数の同定用データに有限インパルス応答モデルを適用し、前記出力変数に自己回帰モデルを適用する、請求項１記載のモジュラー型時系列データ予測装置。
【請求項3】
前記ペアワイズ予測モデル同定部は、複数の前記入力変数の各々ついて１又は複数の遅れ時間を組み込んだ単回帰モデルを用いて、１入力１出力の伝達関数モデルを構築し、前記単回帰モデルを合成したモデルを前記ペアワイズ予測モデルとする、請求項１記載のモジュラー型時系列データ予測装置。
【請求項4】
前記ペアワイズ予測モデル同定部は、複数の前記入力変数各々の同定用データに対して非線形変換を行い、前記ペアワイズ予測モデルを構築する、請求項１記載のモジュラー型時系列データ予測装置。
【請求項5】
前記予測モデル合成法定義部は、複数の前記ペアワイズ予測値の重み付き平均を演算することを前記合成法とし、前記合成法における重みを、前記ペアワイズ予測モデル同定部によって同定された複数の前記ペアワイズ予測モデル各々の予測精度を表す指標に基づいて決定する、請求項１記載のモジュラー型時系列データ予測装置。
【請求項6】
前記予測モデル合成法定義部は、複数の前記ペアワイズ予測値の重み付き平均を演算することを前記合成法とし、前記ペアワイズ予測モデル同定部で用いた同定用データと同じデータを用いて、前記合成法における重みの値を同定して決定する、請求項１記載のモジュラー型時系列データ予測装置。
【請求項7】
前記ペアワイズ出力変数予測部から出力される複数の前記ペアワイズ予測値と、前記合成出力変数予測部から出力される前記出力変数の予測値とに対して、予め指定した所定の周期毎あるいは所定期間で、前記時系列データから抽出した実績値との誤差を評価する予測誤差評価部と、
複数の前記ペアワイズ予測値の誤差と、前記出力変数の予測値の誤差とを比較することにより、調整すべき前記ペアワイズ予測モデル、統合すべき前記ペアワイズ予測モデル、および、分離すべき前記ペアワイズ予測モデルを決定して、前記ペアワイズ出力変数予測部で用いられる複数の前記ペアワイズ予測モデルを修正するペアワイズ予測モデル修正部と、
を更に有する請求項１記載のモジュラー型時系列データ予測装置。
複数の前記ペアワイズ予測値の誤差と、前記出力変数の予測値の誤差とを比較することにより、調整すべき前記ペアワイズ予測モデル、統合すべき前記ペアワイズ予測モデル、および、分離すべき前記ペアワイズ予測モデルを決定して、前記ペアワイズ出力変数予測部で用いられる複数の前記ペアワイズ予測モデルを修正するペアワイズ予測モデル修正部と、
を更に有する請求項１記載のモジュラー型時系列データ予測装置。
【請求項8】
前記ペアワイズ予測モデル修正部は、前記予測誤差評価部によって評価された複数の誤差の少なくとも一つが所定の許容誤差を超過したときに、当該誤差に対応する前記入力変数の候補に関する前記ペアワイズ予測モデルの再同定を行い、前記ペアワイズ出力変数予測部で用いられる複数の前記ペアワイズ予測モデルを修正する請求項７記載のモジュラー型時系列データ予測装置。
【請求項9】
前記ペアワイズ予測モデル修正部は、前記予測誤差評価部によって評価された複数の前記誤差の少なくとも一つが所定の許容誤差を超過したときに、当該誤差に対応する前記入力変数に関する前記ペアワイズ予測モデルを削除して、前記ペアワイズ出力変数予測部で用いられる複数の前記ペアワイズ予測モデルを修正する請求項７記載のモジュラー型時系列データ予測装置。
【請求項10】
前記ペアワイズ予測モデル修正部は、前記予測誤差評価部によって評価された複数の前記誤差を比較し、複数の前記誤差の差が所定の値以下であるときに、当該誤差に対応する前記入力変数の和あるいは平均値を一つの新たな入力変数として定義し、前記ペアワイズ出力変数予測部で当該新たな入力変数に関する前記ペアワイズ予測モデルが用いられるように、複数の前記ペアワイズ予測モデルを修正する請求項７記載のモジュラー型時系列データ予測装置。
【請求項11】
複数のプロセス変数を所定の周期で計測する複数のプロセスセンサを有するシステム又はプロセスに適用される方法であって、
複数の前記プロセス変数の時系列データを所定の周期で収集して保存し、
複数の前記プロセス変数の中から予測対象となる少なくとも一つの出力変数を選択し、
複数の前記プロセス変数の中から複数の入力変数の候補を選択し、
前記時系列データから抽出された前記出力変数と複数の前記入力変数との同定用データを用いて、１入力１出力のペア毎にペアワイズ予測モデルのパラメータを同定して複数の前記ペアワイズ予測モデルを定義し、
複数の前記ペアワイズ予測モデルから出力されるペアワイズ予測値の合成法を定義し、
時間の進行方向に所定の周期又はリアルタイムで前記時系列データから抽出された複数の前記入力変数の予測用データを、複数の前記ペアワイズ予測モデルに入力して、複数の前記入力変数のそれぞれに対応する前記ペアワイズ予測値を演算し、
前記合成法により複数の前記ペアワイズ予測値を合成して前記出力変数の予測値を演算する、
モジュラー型時系列データ予測方法。
複数の前記プロセス変数の時系列データを所定の周期で収集して保存し、
複数の前記プロセス変数の中から予測対象となる少なくとも一つの出力変数を選択し、
複数の前記プロセス変数の中から複数の入力変数の候補を選択し、
前記時系列データから抽出された前記出力変数と複数の前記入力変数との同定用データを用いて、１入力１出力のペア毎にペアワイズ予測モデルのパラメータを同定して複数の前記ペアワイズ予測モデルを定義し、
複数の前記ペアワイズ予測モデルから出力されるペアワイズ予測値の合成法を定義し、
時間の進行方向に所定の周期又はリアルタイムで前記時系列データから抽出された複数の前記入力変数の予測用データを、複数の前記ペアワイズ予測モデルに入力して、複数の前記入力変数のそれぞれに対応する前記ペアワイズ予測値を演算し、
前記合成法により複数の前記ペアワイズ予測値を合成して前記出力変数の予測値を演算する、
モジュラー型時系列データ予測方法。
【請求項12】
コンピュータに、請求項１１に記載の方法を実行させるモジュラー型時系列データ予測プログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明の実施形態は、モジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラムに関する。
【背景技術】
【0002】
上下水道システム、雨水排水システム、電力システム、交通システムなどのインフラシステム、あるいは、鉄鋼プロセス、石油化学プロセス、半導体製造プロセス、などのプロセス系の産業プラントなどでは、通常、複数のプロセス状態を測定する複数のオンラインセンサが設置されている。プロセス監視制御システム（ＳＣＡＤＡ: Supervisory Control And Data Acquisition）と呼ばれるシステムは、上記のインフラシステムや産業プラントに設置されたセンサ群の計測により得られるプロセスデータ（流量、温度、水質、操作量など）を取得し、時系列データとしてサーバ上に保持している。プロセス監視制御システムは、通常、これらの時系列データをトレンドグラフとして監視員に提供することが多い。
【0003】
監視員は、上記のような直接的な時系列データのトレンド監視に加えて、予測情報を利用するプラント監視を行うことがある。予測情報を利用するプラント監視システムによれば、監視中のリアルタイムの実測値（計測値）だけでなく、監視対象の将来の値を予測して、予測結果を運転管理者に提供することにより、適切な運転計画の実現を支援することができる。予測情報を利用するプラント監視システムは、プラントの異常兆候検出および異常診断と並んで、アドバンスト監視の典型的な手段である。
【0004】
予測情報を利用するプラント監視が用いられるインフラシステムの例として、例えば、浄水システムにおける水需要の予測、下水処理場への汚水流入予測、雨天時の下水処理場および雨水排水ポンプ場への雨水流入予測、雨天時の河川水位やダム水位の予測、あるいは、気象レーダなどによる降雨の短期予測（ナウキャスト）、交通システムにおける交通量予測や渋滞予測、電力システムにおける電力需要の予測、あるいは、ビルなどの施設内の食堂や店舗などの顧客数予測、などのシステムがあり、インフラ領域の様々な分野で広く使われている。
【0005】
予測情報を利用するプラント監視システムは、予測情報そのものを付加価値として提供するだけでなく、プラントの運転計画やプラントの制御に予測情報を利用することで、より効率の良い運用やより安全性の高い運用を実現するためにも用いられ得る。例えば、水需要や電力需要等の需要予測は、水運用計画や発電および蓄電計画（ＥＭＳ：エネルギーマネージメントシステム）などの最適化に利用することで、効率の良い運用を実現するために用いることができる。また、例えば降雨予測や流入予測は、雨水排水ポンプの運転制御のための入力情報として予測値を用いることで、浸水の回避や抑制などのリスク低減に貢献する。
【0006】
インフラシステムだけでなくプロセス系のプラントである鉄鋼プロセスや石油化学プロセスなどにおいても、製品品質や歩留まりの予測などに時系列データ解析を用いることも多く、その情報は生産効率（歩留まり）向上のために利用されることも多い。
【0007】
時系列データの予測は、インフラシステムやプロセス系の産業プラントで極めて重要な役割を果たし、プラントの監視データを利用したデータドリブンの（ブラックボックス的アプローチの）予測方法として、従来から様々な方法が適用されている。
【0008】
従来からプラント監視データなどの時系列データを用いた予測技術は広く用いられているが、近年の統計的機械学習分野を中心とするＡＩ関連の理論、技術、手法の急速な進歩に伴い、時系列データ解析に対するＡＩ手法も同時並行的に進展し、分野および領域を問わず、様々なアドバンストなＡＩ手法を時系列データの予測にも適用する動きが加速している。特に、ＮＮの発展形である深層学習ネットワーク（ＤＮＮ：ディープラーニングネットワーク）などは、昨今注目を集めている代表的な方法であり、ＤＮＮはＡＩの代名詞としても使われる場合もある。
【0009】
ＡＩの代表的な手法であるＤＮＮの他にも、重回帰分析で問題になる多重共線性による不安定化（悪条件問題）を回避する方法として（先に述べたＰＬＳとは異なる観点で）導入された正則化に基づく各種の方法、特に、古くからあるＬ２正則化（リッジ回帰）に加え、多数の説明変数候補の中から必要な説明変数を自動選択することのできるＬ１正則化手法であるLassoなどの方法、Lassoとは異なる視点でベイズ推論的な手法で説明変数の自動選択を行うRVM（適合ベクトルマシン）に基づく回帰手法（予測手法）であるＲＶＲ（適合ベクトル回帰）などの方法、最適化を用いて識別問題（分類問題）をロバスト化したＳＶＭ（サポートベクトルマシン）を回帰問題（予測問題）に転用したＳＶＲ（サポートベクトル回帰）などの方法、従来の重回帰などのような線形回帰を非線形に拡張し、さらに確率分布を考えてノンパラメトリック回帰問題に帰着させ、点予測ではなく予測分布を出力するガウス過程回帰などの方法、あるいは、様々な予測サブシステムを組み合わせて予測性能を向上させるバギングやブースティングと呼ばれる各種の方法（ランダムフォレスト、アダブースト、ＸＧブーストなど）など、多様なアドバンストな時系列データの予測手法も用いられる様になってきている。
【0010】
このように、アドバンストなＡＩに基づく予測手法が次々と開発され、従来から用いられる方法と比較して各段に予測精度が向上するという報告も数多くなされる様になり、ＡＩ関連技術に関する期待が高まっている。
【0011】
一方で、ＡＩ関連技術が実際の現実的なシステム（インフラシステムやプロセス系のプラント）において、必ずしも容易に展開できない理由の一つとして、統計的およびＡＩ的な解析手法が持つ本質的な性質であるブラックボックス的なアプローチの欠点として指摘されている、「結果に対する説明性の欠如」の問題がある。そして、この「結果に対する説明性の欠如」の問題に対する意識が急速に高まっており、この問題を解決すべく、最近では「説明可能ＡＩ（ＸＡＩ：eXplainable AI)の概念が米国ＤＡＲＰＡにより提唱され、注目が集まっている。
【先行技術文献】
【特許文献】
【0012】
【特許文献1】特許第６２６１９６０号公報
【発明の概要】
【発明が解決しようとする課題】
【0013】
インフラ系のシステムでは、物理的、化学的な知見に立脚したホワイトボックス的なモデルを用いる事も多いが、これは、「Accountability（説明責任）」が求められることに関連していると考えられる。実際、ホワイトボックス的なモデルによる予測精度がブラックボックス的なモデルによる予測精度よりも必ずしも高いわけではなく、また、ホワイトボックス的なモデルが必ずしも、実際の物理現象を正確かつ忠実に再現しているとも限らない。それにも関わらず、ホワイトボックスモデルが時として好まれ、利用される理由は、ホワイトボックスモデルに含まれる「パラメータ」には物理的および化学的な意味があるため「合理的な説明がしやすく、他者に納得してもらいやすい」ことに加え、このような物理的および化学的な意味を持つことにより、精度が十分でない場合に、どのように調整すればよいかという方針や指針をたてやすい、すなわち、調整しやすい、という利点があるためであると考えられる。実際、物理的および化学的に意味を持つパラメータは、パラメータの値自体に物理的および化学的意味があるため、理解がしやすく、そのパラメータの値の取りうる範囲を想定することや、値の増減が結果に与える影響を物理的および化学的知見に基づいて考えることができる。
【0014】
上記のように、「合理的な説明のしやすさと調整のしやすさ」がホワイトボックスモデルを用いることの大きな動機付けになっている可能性がある。もしＡＩ的な手法によるブラックボックスモデルにおいても「合理的な説明のしやすさと調整のしやすさ」という要素を付加することができれば、その応用範囲は大きく広がると考えられる。
【0015】
本発明の実施形態は、上記事情を鑑みて成されたものであって、予測結果を合理的に説明するとともに、容易に予測結果を調整可能とするモジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラムを提供することを目的とする。
【課題を解決するための手段】
【0016】
実施形態によるモジュラー型時系列データ予測装置は、複数のプロセス変数を所定の周期で計測する複数のプロセスセンサを有するシステム又はプロセスに適用される装置であって、複数の前記プロセス変数の時系列データを所定の周期で収集し保存するとともに、複数の前記プロセス変数の中から予測対象となる少なくとも一つの出力変数を選択する出力変数データ選択部と、複数の前記プロセス変数の中から複数の入力変数の候補を選択する入力変数データ選択部と、を含むデータ収集保存部と、前記時系列データから抽出された前記出力変数と複数の前記入力変数との同定用データを用いて、１入力１出力のペア毎にペアワイズ予測モデルのパラメータを同定して複数の前記ペアワイズ予測モデルを定義するペアワイズ予測モデル同定部と、複数の前記ペアワイズ予測モデルから出力されるペアワイズ予測値の合成法を定義する予測モデル合成法定義部と、時間の進行方向に所定の周期又はリアルタイムで前記時系列データから抽出された複数の前記入力変数の予測用データを、複数の前記ペアワイズ予測モデルに入力して、複数の前記入力変数のそれぞれに対応する前記ペアワイズ予測値を演算するペアワイズ出力変数予測部と、前記合成法により複数の前記ペアワイズ予測値を合成して前記出力変数の予測値を演算する合成出力変数予測部と、を有する。
【図面の簡単な説明】
【0017】
【発明を実施するための形態】
【0018】
本実施形態では、調整可能な時系列デー予測を行うための基本的手段として、モジュラー型の時系列データ予測装置について説明する。本実施形態のモジュラー型の時系列予測装置は、予測モデル（システム）を部分モデル（サブシステム）に分解可能にして、各部分モデルを容易に着脱可能にするものである。
【0019】
モジュラー型の時系列データ予測装置によれば、例えば、複数の監視変数中のいずれか一つの変数のデータが欠測していたり、複数の監視変数にアウトライア（異常データ）が含まれていたり、監視変数がノイズ等を多量に含む不安定な計測データである場合、その影響が全体に波及することを回避できる。また、モジュラー型の時系列データ予測装置によれば、例えば、説明変数が必要であるか否（不要）かを適宜判断して、部分モデルを容易に脱着したり、予測モデルを部分的に調整したりすることができる。
【0020】
このような利点は、モジュラー型の構成によりモデルが「分解可能」であるという事によって得られる利点である。予測モデルにおいてモジュラー型の構成をとる利点は、単に分解可能ということによって得られる利点だけでなく、これに加えて、このモジュラー型の構成をとることで、物理法則に矛盾せずに合理的で納得できるモデルを構築することが可能になり、さらに、このことにより調整を容易にすることができる点である。
【0021】
予測モデルにおいてモジュラー型の構成をとる利点について以下に説明する。
一般に、ブラックボックス的な手法が説明可能（Explainable）でない理由は、入力および出力という外部情報だけしかわからず、内部構造が見えない不可視である事が原因であると考えられている。すなわち、ブラックボックス的な手法は、単に内部構造が見えないだけでなく、そもそも、外部の入出力関係からだけでは、内部構造を一意に決めることができず、同じ入出力関係を与える複数（無数）の内部構造がありうる。このため、内部構造の解釈が難しくなり、場合によっては物理法則に矛盾する不合理なモデルになってしまうことがある。これは、モデルの一意性の問題として、制御理論分野の中に含まれるであるシステム同定分野等で知られた事項である。このようなモデルの一意性の概念は、例えば、「ある入力データと出力データのペアに対して、その入出力関係を表すモデルを唯一に決定することができるか」という問題であり、モデルの可同定問題と呼ばれている。
一般に、ブラックボックス的な手法が説明可能（Explainable）でない理由は、入力および出力という外部情報だけしかわからず、内部構造が見えない不可視である事が原因であると考えられている。すなわち、ブラックボックス的な手法は、単に内部構造が見えないだけでなく、そもそも、外部の入出力関係からだけでは、内部構造を一意に決めることができず、同じ入出力関係を与える複数（無数）の内部構造がありうる。このため、内部構造の解釈が難しくなり、場合によっては物理法則に矛盾する不合理なモデルになってしまうことがある。これは、モデルの一意性の問題として、制御理論分野の中に含まれるであるシステム同定分野等で知られた事項である。このようなモデルの一意性の概念は、例えば、「ある入力データと出力データのペアに対して、その入出力関係を表すモデルを唯一に決定することができるか」という問題であり、モデルの可同定問題と呼ばれている。
【0022】
可同定問題は、モデルの構造（モデルを表現する数式の形）を予め限定した場合、例えば、モデルを線形伝達関数モデルや時系列モデルに限定した場合、そのモデルに含まれるパラメータ（例えば伝達関数に含まれる係数）を一意に決めることができるかという問題となる。この場合、可同定問題の中でも特にパラメータ可同定性の問題と呼ばれる。このパラメータ可同定性の問題は、現実のデータを扱う問題でも頻繁に現れる問題であり、以下の３つのケースは特に遭遇することが多い問題である。
【0023】
一つ目は、先に述べた多重共線性の問題であり、これは、入力変数（説明変数）間に強い相関を持つ場合の問題である。これは、例えば、雨水ポンプ場への雨水流入量を予測する場合に、入力変数としてレーダ雨量情報を用いて、近接する異なる２つのメッシュの降雨量を入力変数とするような場合などに典型的に現れる問題であり、現実のデータでは、説明変数間に強い相関を持つことは極めて多い。
【0024】
二つ目は、制御理論の分野でよく知られている、持続的励振条件（Persistent Exciting Condition、PE条件）と呼ばれる条件が成立しない場合の問題であり、これは、モデル化の対象プロセスに対して、十分な周波数成分を含まない入力データが適用されている場合に生じる問題である。極端ではあるが、現実にも認められる例として、入力データが一定値である場合、通常出力データも一定値になるが、このような一定値の入出力関係を表す伝達関数モデルは無数に存在するため、係数を一意に決めることができない。例えば、y(t)=b1×u(t-1)+b2×u(t-2)と表される簡単な伝達関数モデル（時系列モデル）を考える場合、u(t-1)=u(t-2)=ｋ（一定値）であるとb1とb2を一意に決めることはできない。この概念を一般化したものが持続的励振条件というパラメータ可同定条件であるが、このような条件が成立しない場合も多い。
【0025】
三つ目は、モデルの構造の複雑さと比較して、データの得られる量が少なくデータの質が良くない場合である。典型的な例として、ある対象が非線形的な挙動を示すことが物理的な観点から予めわかっているが、得られるデータが少なかったり、ある特定の運用条件の付近のデータしか得られなかったりすることにより、非線形の係数を同定できない様な場合である。例えば、y(t)=c1×u(t)+c2×u(t)^2と書かれることがわかっているが、データとして得られるu(t)の値が小さい値の場合、u(t)^2<<u(t)となるため、c2を同定することが困難になる様な場合である。
【0026】
上記三つの例は、パラメータの可同定性が欠落する典型的なケースであるが、現実の問題では、このような状況に遭遇することが極めて多い。モデル化の対象が、比較的小さな機械系システムおよび電気系システムの場合は、可同定性を向上させるための実験を行い、可同定性を確保した上で、モデルを構築することが常道である。しかし、上下水道、雨水排水、電力システム、などの巨大な社会インフラシステムの場合、モデル化の対象となるデータが自然現象から得られるデータ（降雨、日射量、風速など）であったり、対象プラントが稼働している状態のデータであったりすることが多く、パラメータの可同定性を向上させるための実験などを行うことができず、モデル化に利用できるデータは与えられたものであることが多いため、可同定性の欠落を回避することが本質的に難しい場合が多い。この結果、社会インフラシステムの場合には、可同定性が欠落した状態でモデルを構築せざるを得ない。
【0027】
以下、可同定性の欠落による問題の一例として「多重共線性」について、従来の技術と本実施形態のモジュラー型時系列データ予測装置とを対比する。
「多重共線性」とは、複数の説明変数（入力変数）間に相関関係（一次従属関係）があることにより、回帰モデルのパラメータを一意に決めることができなくなるという問題である。これを最も簡単な例を用いて説明する。
「多重共線性」とは、複数の説明変数（入力変数）間に相関関係（一次従属関係）があることにより、回帰モデルのパラメータを一意に決めることができなくなるという問題である。これを最も簡単な例を用いて説明する。
【0028】
一つの出力yが二つの入力ｕ１、ｕ２の重回帰で表される以下の重回帰モデルを考える。
y=a1×u1＋a2×u2 （１）
(1)式は最も簡単な重回帰モデルであり、ａ１とａ２とは回帰係数と呼ばれるパラメータである。ここで、物理的な具体的イメージを持つために、例えば、出力ｙを雨水ポンプ場への雨水流入量、入力ｕ１および入力ｕ２を各々異なる箇所の降雨量であることを想定して考えてみる。この場合、降雨量が増加すれば流入量が増加することは、物理的に考えて自明であるから、ａ１とａ２とは正の値を持たなければならないことは容易にわかる。ここで、降雨量ｕ１と降雨量ｕ２とに従属関係のある多重共線性の問題がある場合を想定する。
y=a1×u1＋a2×u2 （１）
(1)式は最も簡単な重回帰モデルであり、ａ１とａ２とは回帰係数と呼ばれるパラメータである。ここで、物理的な具体的イメージを持つために、例えば、出力ｙを雨水ポンプ場への雨水流入量、入力ｕ１および入力ｕ２を各々異なる箇所の降雨量であることを想定して考えてみる。この場合、降雨量が増加すれば流入量が増加することは、物理的に考えて自明であるから、ａ１とａ２とは正の値を持たなければならないことは容易にわかる。ここで、降雨量ｕ１と降雨量ｕ２とに従属関係のある多重共線性の問題がある場合を想定する。
【0029】
この多重共線性問題は、実際の現象においても頻繁に見られることであり、この場合、例えば、入力ｕ１と入力ｕ２とをレーダ雨量データの隣接メッシュの降雨量であることを想定してみれば、この問題が特殊な問題ではなく、むしろ自然な問題であることが容易にわかる。例えば、国土交通省が配信するＸＲＡＩＮと呼ばれる気象レーダでは、１メッシュのサイズが２５０ｍ×２５０ｍであるため、隣接メッシュの降雨量とは互いに数百メートル離れた箇所の降雨量を意味しており、その降雨量がほぼ同じであることは物理的に考えて極めて自然な事である。２つの降雨量データがレーダ雨量の隣接メッシュの降雨量ではなくて、２か所の地上雨量データであったとしても、通常雨水排水区の大きさが数キロ四方程度以内であることを考えると比較的近い値を持つことは容易に想像できる。従って、このような場合ｕ１≒ｕ２となっていることは、自然な事であり、これは多重共線性問題そのものである。
【0030】
そこで、ここでは簡単のため、ｕ１＝ｕ２という完全な多重共線性の状態にあった場合を想定する。そして、例えば、２０ｍｍ／ｈの強度の雨量に対して、１０ｍ３／ｓの雨水流入があるとする。今、ｕ１＝ｕ２であるので、ｕ：＝ｕ１＝ｕ２と定義すると、（１）式は、以下の（２）式の様に書き換えられる。
【0031】
y=(a1+a2)×u （２）
上記(2)式は単回帰の形であり、ｕ：＝ｕ１＝ｕ２であるから、ａ１＋ａ２＝０．５となる。しかし、ａ１＋ａ２＝０．５という制約条件を満足してさえいれば、ａ１とａ２とがいかなる値であっても、（２）式は全く同じ出力結果を出すため、ａ１とａ２の各値はａ１＋ａ２＝０．５を満たす範囲で任意の値を持つことができる。そのため、例えば、ａ１＝３．５、ａ２＝－３でも、ａ１＝１００、ａ２＝－９９．５となっても問題がなく、全く同じ入出力関係を与えるため、各々のパラメータ値の相違を入出力データから識別することができなくなる。
上記(2)式は単回帰の形であり、ｕ：＝ｕ１＝ｕ２であるから、ａ１＋ａ２＝０．５となる。しかし、ａ１＋ａ２＝０．５という制約条件を満足してさえいれば、ａ１とａ２とがいかなる値であっても、（２）式は全く同じ出力結果を出すため、ａ１とａ２の各値はａ１＋ａ２＝０．５を満たす範囲で任意の値を持つことができる。そのため、例えば、ａ１＝３．５、ａ２＝－３でも、ａ１＝１００、ａ２＝－９９．５となっても問題がなく、全く同じ入出力関係を与えるため、各々のパラメータ値の相違を入出力データから識別することができなくなる。
【0032】
一方、ａ１とａ２とは、各雨量に対する係数であるから、本来、すくなくとも正の値を持つべきであり、２地点の降雨量が同じであるなら、その値はほぼ等しくなるべきであることも物理的に考えれば妥当である。しかし、このように、入力変数間の相関関係（一次従属関係）により、パラメータ値を一意に決めることができず、本来のパラメータ値を正しく同定できなくなる事が多重共線性の問題の本質である。
【0033】
このような多重共線性などのパラメータ可同定性の欠如により、入出力関係から一意にパラメータを決めることができなくなる問題は、システム同定分野や統計分野（統計的機械学習分野≒AI分野）では、数学的な「悪条件問題」として扱われ、この悪条件を回避するための各種のアルゴリズムが開発されている。先に述べた正則化やＰＬＳなどの方法は、このような悪条件回避のための代表的な技術的手段（テクニック）であり、アドバンストなＡＩ手法の多くには、正則化や（ＰＬＳで用いられる）次元圧縮などの手法(テクニック)が採用されている。
【0034】
しかし、例えば、正則化やＰＬＳなどのテクニックを駆使したアドバンストなアルゴリズムを（１）式の問題に適用したとしても、上記の問題に関して、ａ１＞０、ａ２＞０、かつａ１≒ａ２となるような値を推定するとは限らず、むしろ、そのような値にならない可能性の方が高い。これは、正則化や次元圧縮などは、多重共線性（やこれを一般化した可同定性の欠如）の問題を、数学的な悪条件の問題として捉え、これを回避するためのある種の“テクニック”であり、悪条件を回避してパラメータを推定（同定）できるようにしているだけであるからである。すなわち、多重共線性を持つ重回帰問題に対して古典的な最小２乗法を適用すると、「答えを求めることができない＝パラメータ値を求めることができない」という問題になるため、悪条件を回避するというテクニック（正則化や次元圧縮）をアルゴリズムに導入することで、「答えを求めることができる＝パラメータ値を同定できる」問題にしているだけであり、同定した値が物理的に合理的に納得できるものであるか否かを直接考慮しているわけではないからである。
【0035】
実際、例えば、正則化のテクニックでは、悪条件を回避するために、パラメータを同定する際に、予測誤差を評価するだけでなく、「パラメータの値自身が大きくなりすぎない」というパラメータに関するノルム（大きさを図る指標）を小さくする、という評価指標を加えることで悪条件を回避しており、これはパラメータの物理的な意味づけや解釈を行うこととは直接関係していない。
【0036】
一方、このような問題に対して、合理的に納得のできるパラメータ値を得る方法について以下に説明する。はじめに、（１）式のモデルではなく、（３）式と（４）式との単回帰モデルを考える。
y=a1´×u1 （３）
y=a2´×u2 （４）
y=a1´×u1 （３）
y=a2´×u2 （４）
【0037】
この時、ｕ１＝ｕ２であり、ｕ１＝ｕ２＝２０ｍｍ／ｈに対して、ｙ＝１０ｍ３／ｓを出力するモデルは、単位系をそのままで考えると、ａ１´＝ａ２´＝０．５であることは明らかである。従って、このような関係を持つ場合に、（３）式と（４）式と単回帰を適用すると、ａ１´＝ａ２´＝０．５に近い値がパラメータとして同定（学習）されることになる。この０．５という値は、（回帰モデルなので物理的な意味は明確ではないが）合理的で納得できる値である。
【0038】
次に、（３）式と（４）式とを用いて、（１）式の重回帰モデルを構築することを考える。
今回の例では、入力変数が２か所の降雨量となり、その和をとっているため、入力となる降雨量が倍になるので、各降雨量に対する係数のａ１とａ２とは、一つの降雨量を入力とした単回帰モデルの係数の半分となり、ａ１＝ａ２＝０．２５となると考える事が物理的に妥当な納得できるパラメータ値であると考えられる。
今回の例では、入力変数が２か所の降雨量となり、その和をとっているため、入力となる降雨量が倍になるので、各降雨量に対する係数のａ１とａ２とは、一つの降雨量を入力とした単回帰モデルの係数の半分となり、ａ１＝ａ２＝０．２５となると考える事が物理的に妥当な納得できるパラメータ値であると考えられる。
【0039】
このようにして、（３）式と（４）式とでｕ１とｙとの関係、および、ｕ２とｙとの関係を単回帰で各々求めて、２つの予測モデルを構築する。そして、次に２つの予測モデルの出力ｙの平均を計算する。このように、構成的（Constructive）な手順を踏むと、ａ１とａ２との値は、各々０．２５と、物理的に妥当と考えられる納得のできる値として推定され、結果的に、（１）式の形の重回帰式の係数を定めることができる。これは、部分モデル（サブシステム）を結合していくモジュラー型のアプローチそのものである。
【0040】
上記重回帰モデルにより係数を定めることができる理由は、一般に、入出力関係を表すデータとして同じデータを用いる場合、モデル構造が簡単（≒モデルに含まれるパラメータ数が少ない）なモデルの方が、可同定性が高くなる事（＝パラメータ値を一意に同定可能）が知られているためである。本実施形態のモジュラー型時系列データ予測装置は、上記の特性を利用している。すなわち、同じ入出力データであったとしても、その可同定性はモデルの構造やパラメータ数に依存し、パラメータ数が多いモデルでは可同定でなくなる場合であっても、パラメータ数が少なければ可同定になる場合がある、ということを利用している。
【0041】
上記の例でいえば、本実施形態のモジュラー型時系列データ予測装置では、２入力の重回帰モデルでは多重共線性により可同定でなくなるが、これを１入力の単回帰の組み合わせと考えることにより、多重共線性の問題を本質的に回避して、可同定なモデルを積み上げることによって合理的なモデル構築を可能にしている。
【0042】
次に、本実施形態のモジュラー型時系列データ予測装置と、いわゆる「Explainable AI（ＸＡＩ）」との関係および相違について説明する。
いわゆるＸＡＩでは、全体を一つの大きなシステムとして捉えてブラックボックス的にモデルを構築して、後から、そのモデルを適切に解釈可能な形に分解して説明を加えようとするトップダウン型のアプローチ（あるいは帰納的なアプローチ）である。これに対し、本実施形態のモジュラー型時系列データ予測装置の考え方は、個別には比較的簡単に理解できる合理的な説明が可能なモデル（サブシステム）を統合していく事により、全体の予測モデルを構築していくボトムアップ型のアプローチ（あるいは演繹的なアプローチ）である。
いわゆるＸＡＩでは、全体を一つの大きなシステムとして捉えてブラックボックス的にモデルを構築して、後から、そのモデルを適切に解釈可能な形に分解して説明を加えようとするトップダウン型のアプローチ（あるいは帰納的なアプローチ）である。これに対し、本実施形態のモジュラー型時系列データ予測装置の考え方は、個別には比較的簡単に理解できる合理的な説明が可能なモデル（サブシステム）を統合していく事により、全体の予測モデルを構築していくボトムアップ型のアプローチ（あるいは演繹的なアプローチ）である。
【0043】
上記のようなボトムアップ型のアプローチでは、サブシステム間の相互の干渉を予め考慮することができないため、予測精度面に注目した場合には、トップダウン的（ＸＡＩ的）なアプローチに劣る可能性を否定できない。しかしながら、ボトムアップ型のアプローチは、合理的な説明性という意味では、先に説明したようにパラメータの可同定性を維持したモデル構築が可能になり、トップダウン的なアプローチよりも優れる可能性が高い。
【0044】
また、ボトムアップ型のアプローチは、ホワイトボックス的なモデリングのアプローチとも類似している。ホワイトボックスモデリングは、基本的に演繹的、構成的なモデリング手法である。ホワイトボックスモデリングでは、物理的に意味のあるパラメータを用いて構成的にモデルを構築するのに対し、本実施形態のモジュラー型アプローチでは、物理的な意味を持たせることはできないものの、個々のパラメータを入出力と直接関係づけることにより解釈が容易なパラメータを同定し、それを用いて構成的にモデル構築している。
【0045】
本実施形態のモジュラー型アプローチがＸＡＩ的なアプローチよりも合理的な説明性において優れていることは以下の事からも明らかである。例えば、先の例において、もしトップダウン的（ＸＡＩ的）なアプローチを行おうとすると、正則化やＰＬＳなどの方法を使って係数パラメータのａ１やａ２を求め、それに対して何等かの解釈を加える処理を施して結果を「Explainable（説明可能）」にする方法の確立を目ざす。しかしながら、上記の簡単な例からわかる様に、ａ１とａ２が０．２５という値からかけ離れた値で同定されていた場合（実際、正則化やＰＬＳを用いると、そのようになる可能性がある）には、いかなる方法を用いて「説明」を試みても、合理的な説明が原理的に困難であることは容易にわかる。例えば、ＰＬＳや正則化を用いてパラメータ値を求めることができても、ａ１かａ２のいずれかが負の値（例：ａ１＝３．５、ａ２＝－３）となっている場合には、合理的に納得のできる説明は不可能である。
【0046】
つまり、本来合理的に納得できる値である0.25という値を求められていなかったとしても、可同定性が欠如している場合には、ある制約さえ満たしていればどのような内部パラメータ値であろうが予測精度だけに着目する限りは差異が見られないため、「精度の良いモデル」と見なされる。しかしながら、精度が良いモデルが、必ずしもそのモデルのパラメータ値の妥当性を補償する合理的で妥当なモデルを意味しないため、トップダウン的（ＸＡＩ的）なアプローチでは原理的に合理的な説明が困難になってしまう。
【0047】
なお、ＸＡＩ的なアプローチでは、ａ１とａ２とのそれぞれに対する個別の意味付けは困難であるが、ａ１＋ａ２に対する意味づけは可能であるため、説明不可能であるというわけではない。ただし、ＸＡＩ的なアプローチでは、内部に含まれる全てのパラメータに対して意味付けを行うことは原理的に困難になる可能性が高く、何等かの「特徴量」を抽出して、それに対して意味付けを求めることになる。
【0048】
従って、「合理的な説明性の高いモデル」を構築することを主目的にするのであれば、そもそも、正則化やＰＬＳなどの数学的に悪条件を回避する様なアルゴリズムに頼る前に、可同定性が欠如する様な状態を回避してモデルの内部構造（内部パラメータ）の一意性を保証する様にモデルを構築しておき、合理的な説明性を担保しながらモデルを統合していくボトムアップ的、構成的、演繹的なモジュラー型アプローチの方が優れていると可能性が高い。
【0049】
また、このように合理的説明性を担保したモデルを構築できることに加え、最初に述べた分離可能な構造とすることで、モデルの調整も容易にすることができ、これらを合わせて合理的説明可能かつ調整可能なモデル構築を行うことが可能になる。
本実施形態のモジュラー型時系列データ予測装置は、以上の洞察に基づいて成されたものである。
本実施形態のモジュラー型時系列データ予測装置は、以上の洞察に基づいて成されたものである。
【0050】
以下、実施形態のモジュラー型時系列データ予測装置について、図面を参照して詳細に説明する。
図１は、一実施形態のモジュラー型時系列データ予測装置を適用した雨量流入予測システムを概略的に示す図である。
図１は、一実施形態のモジュラー型時系列データ予測装置を適用した雨量流入予測システムを概略的に示す図である。
【0051】
以下では、本実施形態のモジュラー型時系列データ予測装置を搭載したシステムの実施イメージとその効果をより明確にするために、雨水流入予測システムを対象として説明する。なお、本実施形態のモジュラー型時系列データ予測装置が適用される対象プロセスは雨水流入予測システムに限定されるものではなく、２種類以上の計測項目の時系列データが存在する任意のプロセスに対して適用することができる。モジュラー型時系列データ予測装置は、例えば、鉄鋼プラント、石油化学プラント、食品プラント、製薬プロセス、半導体製造プロセス、発電プラント、交通監視設備、空調監視設備、などのプロセスに適用することができる。
【0052】
雨量流入予測システムの対象プロセス１は、流量計１１と、幹線流量計１２と、幹線水位計１３１-１３Ｋと、地上雨量計１４１-１４Ｍと、レーダ雨量計１５と、ポンプ井水位計１６１、１６２と、流入渠１７と、雨水ポンプ井１８と、雨水ポンプ１９と、流入ゲート１１０と、を含む都市雨水排水プロセス１である。レーダ雨量計１５は、Ｑ×Ｐメッシュのメッシュ毎のレーダ雨量計１５１１～１５ＱＰを含む。
【0053】
都市雨水排水プロセス１の各種センサ１１－１５、１６１、１６２は、所定の周期（例えば、３０秒、６０秒）で、プロセスの状態を表す量や運転操作に関わる量などの計測を行う。各種センサ１１－１５、１６１、１６２で計測された値は、モジュラー型時系列データ予測装置で収集および保存される。
【0054】
流入渠１７は、下水管やポンプ場から送られてきた雨水が流れ込む渠である。
雨水ポンプ井１８は、流入渠１７から流入した雨水が貯められる貯水槽である。雨水ポンプ井１８に流入する雨水は、沈砂池で一緒に流れてきた砂などが予め取り除かれていてもよい。
雨水ポンプ１９は、雨水ポンプ井１８に貯められた雨水を強制的に河川等へ送り出すポンプである。
雨水ポンプ井１８は、流入渠１７から流入した雨水が貯められる貯水槽である。雨水ポンプ井１８に流入する雨水は、沈砂池で一緒に流れてきた砂などが予め取り除かれていてもよい。
雨水ポンプ１９は、雨水ポンプ井１８に貯められた雨水を強制的に河川等へ送り出すポンプである。
【0055】
流入ゲート１１０は、流入渠１７と雨水ポンプ井１８との間の流路に設けられ、開閉することにより流入渠１７から雨水ポンプ井１８に流入する雨水の量を調整するように動作する。流入ゲート１１０の前後には水位計１６１、１６２が設置され、例えば、水位計１６１、１６２により測定された値に応じて流入ゲート１１０の動作が制御される。
【0056】
図２は、一実施形態のモジュラー型時系列データ予測装置の一構成例を概略的に示す図である。
本実施形態のモジュラー型時系列データ予測装置は、データ収集保存部２と、データ抽出部３と、ペアワイズ予測モデル同定部４と、予測モデル合成法定義部５と、ペアワイズ出力変数予測部６と、合成出力変数予測部７と、予測誤差評価部８と、ペアワイズ予測モデル修正部９と、出力予測結果観測部１０と、を備えている。
本実施形態のモジュラー型時系列データ予測装置の構成は、例えば、ＣＰＵやＭＰＵなどのプロセッサを少なくとも１つと、プロセッサにより実行されるプログラムが格納されるメモリとを備えた演算装置であって、種々の機能をソフトウェアにより、若しくは、ソフトウェアとハードウエアとの組み合わせにより実現することが出来る。
本実施形態のモジュラー型時系列データ予測装置は、データ収集保存部２と、データ抽出部３と、ペアワイズ予測モデル同定部４と、予測モデル合成法定義部５と、ペアワイズ出力変数予測部６と、合成出力変数予測部７と、予測誤差評価部８と、ペアワイズ予測モデル修正部９と、出力予測結果観測部１０と、を備えている。
本実施形態のモジュラー型時系列データ予測装置の構成は、例えば、ＣＰＵやＭＰＵなどのプロセッサを少なくとも１つと、プロセッサにより実行されるプログラムが格納されるメモリとを備えた演算装置であって、種々の機能をソフトウェアにより、若しくは、ソフトウェアとハードウエアとの組み合わせにより実現することが出来る。
【0057】
データ収集保存部２は、各種センサ１１－１５、１６１、１６２で計測されたプロセス変数の値を収集し保存する。データ収集保存部２は、出力変数データ選択部２１と、入力変数データ選択部２２とを備えている。
【0058】
出力変数データ選択部２１は、各種センサ１１－１５、１６１、１６２で計測された値の少なくともいずれかの量を出力変数として選択する。
入力変数データ選択部２２は、各種センサ１１－１５、１６１、１６２で計測された値から、複数の入力変数として選択する。なお複数の入力変数の数は、複数のプロセス変数の数以下であればよい。
入力変数データ選択部２２は、各種センサ１１－１５、１６１、１６２で計測された値から、複数の入力変数として選択する。なお複数の入力変数の数は、複数のプロセス変数の数以下であればよい。
【0059】
データ抽出部３は、データ収集保存部２に保存された時系列データの中から所定のデータを抽出する。データ抽出部３は、オフライン予測モデル同定用データ抽出部３１と、オンライン予測用データ抽出部３２と、評価用データ抽出部３３とを備えている。
【0060】
オフライン予測モデル同定用データ抽出部３１は、データ収集保存部２から、所定の同定周期若しくは外部からの要求により、予測モデル同定用のデータとして各種プロセスセンサ１１－１５、１６１、１６２で計測された（当該時刻よりも過去の）所定の期間のデータを抽出する。
【0061】
オンライン予測用データ抽出部３２は、データ収集保存部２から、所定の予測周期で、各種プロセスセンサ１１－１５、１６１、１６２の計測値からオンラインの予測に必要となる予測用データを、リアルタイムで抽出する。
【0062】
評価用データ抽出部３３は、データ収集保存部２から、所定の評価周期もしくは外部からの要求により、予測誤差評価に必要となる各種プロセスセンサ１１－１５、１６１、１６２の計測値を、リアルタイムで抽出する。
【0063】
ペアワイズ予測モデル同定部４は、オフライン予測モデル同定用データ抽出部３１で抽出された計測値のデータを用いて、出力変数データ選択部２１で選択された出力変数および入力変数データ選択部２２で選択された入力変数を変数毎にペアワイズとし、後述する方法により予測モデルを同定する。
【0064】
予測モデル合成法定義部５は、ペアワイズ予測モデル同定部４で同定したペアワイズの予測モデルの出力値（ペアワイズ予測値）から、出力変数（出力変数が複数の場合は出力変数毎）を合成する方法を定義する。
【0065】
ペアワイズ出力変数予測部６は、ペアワイズ予測モデル同定部４で同定したペアワイズの予測モデルに対して、時間の進行方向における所定の監視周期で、オンライン予測用データ抽出部３２から抽出した入力変数各々に対応する予測用データを入力して、入力変数毎に出力変数（出力変数が複数の場合は出力変数毎）のペアワイズの予測を行う。
【0066】
合成出力変数予測部７は、ペアワイズ出力変数予測部６の出力であるペアワイズの予測値に対して、予測モデル合成法定義部５で定義した方法に従って合成した出力変数の予測値を出力する。
【0067】
予測誤差評価部８は、ペアワイズ出力変数予測部６の出力であるペアワイズの予測値と、合成出力変数予測部７の出力であるペアワイズの予測値とを合成した予測値を保存し、所定の周期もしくは所定のタイミングで、評価用データ抽出部３３により抽出した出力変数の実績値と保存したペアワイズおよび合成予測値との誤差を、各々評価する。
【0068】
ペアワイズ予測モデル修正部９は、予測誤差評価部８によって評価した予測誤差に基づいて、ペアワイズ予測モデルの調整の必要性の有無を判断し、必要なペアワイズ予測モデルの修正を行う。ペアワイズ予測モデル修正部９は、例えば、ペアワイズ予測モデルのパラメータ調整、もしくは、ペアワイズ予測モデルの削除、もしくは、ペアワイズ予測モデルの統合、を行うことができる。ペアワイズ予測モデル修正部９は、修正したペアワイズ予測モデルをペアワイズ出力変数予測部６に供給する。
【0069】
出力予測結果観測部１０は、例えばモニタなどの表示手段と、ユーザインタフェースとを備え、合成出力変数予測部７による合成出力変数の予測値を少なくとも観測者に提示することができる。出力予測結果観測部１０は、例えば、パーソナルコンピュータ、タブレット端末や、スマートフォンなどの携帯端末であってもよい。出力予測結果観測部１０は、必要に応じて、ペアワイズ出力変数予測部６によるペアワイズ予測モデルの複数の出力と、それから合成される予測分布情報の少なくとも一方を、観測者に更に提示してもよい。
【0070】
次に、本実施形態のモジュラー型時系列データ予測装置の動作について説明する。
例えば都市雨水排水プロセス１では、各種プロセスセンサ１１－１５、１６１、１６２によって、所定の周期でプロセスの情報が計測され、データ収集保存部２に計測値（プロセス変数の値）が供給される。データ収集保存部２は、収集した計測値を、予め決められたフォーマットに従って、時系列データとして保存する。
例えば都市雨水排水プロセス１では、各種プロセスセンサ１１－１５、１６１、１６２によって、所定の周期でプロセスの情報が計測され、データ収集保存部２に計測値（プロセス変数の値）が供給される。データ収集保存部２は、収集した計測値を、予め決められたフォーマットに従って、時系列データとして保存する。
【0071】
出力変数データ選択部２１は、各種プロセスセンサ１１－１５、１６１、１６２で計測されたプロセス変数の中から予測の対象となる変数を選択して、選択した変数を出力変数と設定する。
例えば都市雨水排水プロセス１では、通常雨水ポンプ井１８への流入量を予測の対象とする場合には、出力変数データ選択部２１は、流量計（雨水ポンプ井流入量計）１１で計測されている雨水ポンプ井流入量が予測対象の出力変数として選択する。
例えば都市雨水排水プロセス１では、通常雨水ポンプ井１８への流入量を予測の対象とする場合には、出力変数データ選択部２１は、流量計（雨水ポンプ井流入量計）１１で計測されている雨水ポンプ井流入量が予測対象の出力変数として選択する。
【0072】
また、流量計（雨水ポンプ井流入量計）１１が設置されていない様な場合には、雨水ポンプ井水位計１６１、１６２の計測値と、雨水ポンプ１９および流入ゲート１１０の運用状態とから、雨水ポンプ井１８への流入量を計算によって求めて、これを予測対象として選択しても良い。この場合、出力変数データ選択部２１は、流入量の計算値を計算して保存しておき、流入量の計算値を出力変数として選択する。
【0073】
また、例えば都市雨水排水プロセス１の複数の項目を予測する場合には、出力変数データ選択部２１は、複数の項目を出力変数として選択してもよい。都市雨水排水プロセス１においても、例えば、雨水ポンプ井１８に複数の箇所から雨水の流入があるような場合は、出力変数データ選択部２１は、複数の流入箇所各々における流入量を出力変数として選択してもよい。
【0074】
入力変数データ選択部２２は、出力変数データ選択部２１で選択した出力変数に影響を与える可能性のある測定値を入力変数の候補として選択する。影響を与えるか否かについて事前情報が全く無い場合には、入力変数データ選択部２２は、全てのプロセス変数を機械的に入力変数としておいても良い。入力変数データ選択部２２は、出力変数として選択した変数を入力変数としても選択して構わない。例えば、時系列のデータ解析では、出力変数として選択した変数自身の過去の値から未来の値を予測する所謂自己回帰を検討することも多くある。このような場合には、流量計（雨水ポンプ井流入量計）１１による雨水ポンプ井流入量は、出力変数と入力変数との両方として選択されても構わない。
【0075】
また、出力変数データ選択部２１で選択した出力変数が複数ある場合、入力変数データ選択部２２は、ある出力変数に対して別の出力変数を入力変数として選択してもよい。このようにすることで、複数の出力変数同士の相関を考慮することが可能になる。一方、複数の出力変数が異なる流入箇所の流入量であるときには、プロセスの物理的な構造上、出力変数同士の関連が無い場合が多い。この場合、入力変数データ選択部２２は、ある出力変数に対する入力変数として別の出力変数を除外しておくこともできる。
【0076】
入力変数データ選択部２２は、その他の変数に対しても、例えば、物理的な管渠の接続関係などの情報から明らかに出力変数に無関係な変数や因果関係（入力（原因）と出力（結果）との関係）が物理的に成立しない変数を予め入力変数から除外しておくこともできる。
【0077】
例えば、雨水ポンプ井水位計１６１、１６２により計測される雨水ポンプ井水位は、出力変数として選択した流量計（雨水ポンプ井流入量計）１１による雨水ポンプ井流入量の影響を受けて変化することは物理的に考えて明らかであるから、雨水ポンプ井水位を入力変数から除外しておく。
【0078】
入力変数データ選択部２２は、物理的な構造の関係が明確でない変数について予め出力変数データ選択部２１で選択した出力変数との相関解析を予備的に実施し、相関係数の絶対値が所定の値（例えば０．５）未満の変数を、予め入力変数から除外しておいてもよい。さらに、入力変数データ選択部２２は、相関解析を実施する際に、解析対象の変数が測定された時間をずらしながら相関解析を行い、最も相関の高くなる時のずらした時間Ｌを求め、時間Ｌの符号情報から、因果関係を満たさない符号を持つ変数を入力変数から予め除外しておいてもよい。なお、先の雨水ポンプ井水位の変数は、このような解析によっても入力変数から除外される可能性が高い。
【0079】
本実施形態のモジュラー型時系列データ予測装置では、入力変数データ選択部２２は、雨水水位計１６１、１６２により計測された雨水ポンプ井水位のみを除いた流量計（雨水ポンプ井流入量計）１１による雨水ポンプ井流入量と、幹線流量計１２により計測された幹線流入量と、幹線水位計１３１-１３Ｋにより計測された幹線水位と、地上雨量計１４１-１４Ｍにより計測された地上雨量と、Ｑ×Ｐメッシュのメッシュ１５１１-１５ＱＰ各々におけるレーダ雨量とを入力変数として選択するものとする。
【0080】
なお、本実施形態のモジュラー型時系列データ予測装置において、例えば、Ｋ個の幹線水位計１３１～１３Ｋは、最近、下水管渠内の水位情報を正確に把握しようとする行政の動向などから、Ｋの値が１０～１００程度の比較的多数の幹線水位計が設置されることも多く、また、気象レーダは国土交通省が全国にＸＲＡＩＮという気象レーダを設置していることから、多くの場所で取り入れることができ、対象とする雨水排水区の大きさにもよるが、Ｑ×Ｐメッシュのメッシュサイズは１００程度になることも珍しくないため、実際に計測される時系列データの項目数（変数の数）は、かなり多くなる。また、本実施形態のモジュラー型時系列データ予測装置が適用される都市雨水排水プロセス１に限らず、計測変数の項目数が数百～数千に及ぶことは一般的である。
【0081】
本実施形態のモジュラー型時系列データ予測装置では、入力変数は、例えば、雨水ポンプ井水位計１６１、１６２による計測値のみを除いた流量計（雨水ポンプ井流入量計）１１による雨水ポンプ井流入量の計測値と、幹線流量計１２による幹線流入量の計測値と、幹線水位計１３１～１３ＫによるＫ個の幹線水位の計測値と、地上雨量計１４１～１４ＭによるＭ個の地上雨量の計測値と、Ｑ×Ｐメッシュ１５１１～１５ＱＰの各メッシュのレーダ雨量の計測値と、である。以下では、入力変数の数をｐ（正の整数）として説明する。
【0082】
オフライン予測モデル同定用データ抽出部３１は、予測モデルを構築する人が指定した所定期間のデータを用いて要求されたタイミングで、先に出力変数データ選択部２１と入力変数データ選択部２２との該当する所定期間分の時系列データを取得する。オフライン予測モデル同定用データ抽出部３１は、予測モデルの構築を定期的に行い、モデルを定期的に更新したい場合には、時間の進行方向において所定の周期ＴＬで所定期間の入出力変数の時系列データを定期的に抽出しても良い。このようにして抽出されたデータセットをｕｋ、ｋ＝１、２、…ｐと記載する。ｕｋは列ベクトルであるとし、各行が時系列データの各時刻のサンプルに対応する。同様に出力変数である雨水流入量の時系列データセットｙとしておく、オフライン予測モデル同定用データ抽出部３１は、データセットｕｋとデータセットｙとを所定の周期で抽出する作用を持つ。
【0083】
次に、ペアワイズ予測モデル同定部４では、ｕｋとｙとを用いて、（ｕ１、ｙ）、（ｕ２、ｙ）、…（ｕｐ、ｙ）のｐ個のペア毎に予測モデルを同定する。本実施形態のモジュラー型時系列データ予測装置では、出力変数ｙは、流量計（雨水ポンプ井流入量計）１１による雨水ポンプ井流入量のみの一つであるが、複数（ｈ個とする）の出力ｙがある場合には、ｈ個の各々のｙに対して同様の処理を行えば良く、以降の処理は出力変数の数がｌ個に増えた場合でも、ｈ回繰り返して同じ処理を行えばよいだけであるので、一般性を失う事なく出力変数は１個であるとしておく。
【0084】
最も典型的な予測モデルは、次式で表せる線形回帰の形をしたモデルを仮定して同定を行うものである。
y(t)=a1×uk(t-L)+a2×uk(t-L-1)+…+an×uk(t-L-n+1)+c （５）
ここで、ａｋ、ｋ＝１、２、…、ｎとｃは同定すべきパラメータである。ｃは入力ｕｋと出力ｙの平均の差を表すバイアスパラメータであるが、予めｕｋとｙから平均値を除去して平均を０としておくことにより、常に０とすることができるので、以下では簡単のためにｃは０としておく。
y(t)=a1×uk(t-L)+a2×uk(t-L-1)+…+an×uk(t-L-n+1)+c （５）
ここで、ａｋ、ｋ＝１、２、…、ｎとｃは同定すべきパラメータである。ｃは入力ｕｋと出力ｙの平均の差を表すバイアスパラメータであるが、予めｕｋとｙから平均値を除去して平均を０としておくことにより、常に０とすることができるので、以下では簡単のためにｃは０としておく。
【0085】
上記（５）式において、ｕｋがｙ自身ではない場合を、ディジタル信号処理やシステム同定の分野では、有限インパルス応答モデル（ＦＩＲ）モデルと呼ぶ（入力変数に対する重み付き移動平均モデルと呼ぶこともできる）。一方、ｕｋがｙ自身と一致する場合は自己回帰モデル（ＡＲモデル）と呼ばれる。ここで重要な事は（５）式の線形回帰の形であっても、複数の入力変数による有限インパル応答を足し合わせ、さらに自己回帰を足し合わせた多入力１出力（複数の出力を同時に考える場合は多入力多出力）の線形回帰の形（ＡＲＸ（Autoregressive eXogenous Input）などの形）ではなくて、あくまでも１入力１出力の形の線形回帰モデルになっている点であり、これが「ペアワイズ」の意味である。
【0086】
ペアワイズで同定することの第一の利点は、「複数の入力」という概念自身が無くなるため、先に述べた（複数の入力変数間の関係から生じる）多重共線性の問題が決して生じる事はなく、この問題を本質的に回避できることを含め、パラメータの可同定性の観点では多入力１出力のシステムよりも有利であるからである。これを若干補足する。システム同定の分野では、入力信号ｕｋに含まれる周波数成分の量（数）によって同定可能なパラメータの数が変化することが持続的励振条件（ＰＥ（Persistent Exciting）条件）と呼ばれるパラメータの可同定性（＝パラメータを一意に決定可能）条件として知られており、同定（推定）すべきパラメータａｋの数が多くなるほど、入力信号ｕｋは多数の周波数成分を含まなければならない事が知られている。
【0087】
一方、本実施形態のモジュラー型時系列データ予測装置において、入力変数の値は、予測モデルの構築者が制御（調整）できるものではなく、所与の観測（計測）情報であるから、どれくらいの周波数成分を含んでいるか否かは、外部条件として与えられるものであり、パラメータの数は少なければ少ないほど可同定性が向上し、逆にパラメータの数が多ければ多いほど可同定性が劣化する可能性が高くなる。より厳密に述べると、パラメータの数が少ない場合に可同定性が悪くなる可能性は０％であり、パラメータ数が少ないものの可同定性の方が高いか、悪くてもパラメータ数が多いものの可同定性と同じである。従って、上記（５）式の１入力１出力のモデルではｎ個のパラメータを含むのに対し、ｐ入力１出力の予測モデルを考えるとｐ×ｎ個のパラメータを含むことになるため、入力数ｐが多くなればなるほどパラメータを一意に同定することが難しくなる可能性が高くなる事がわかる。
【0088】
さらに、先に述べた入力変数間の従属関係（相関関係）があると、多重共線性の問題が生じ、同定すべきパラメータ数の多寡に関わらず可同定性が失われ、入力情報は所与のものであるから、いかなる方法を用いても原理的にパラメータ値を一意に決定することが不可能になる。これに対し、１入力１出力の（５）式のモデルでは多重共線性の問題が原理的に生じることはなく、（５）式の形でパラメータが同定できなくなる場合は、入力ｕｋに十分な周波数情報が含まれていない時に限られ、この場合は、次数ｎを調整して同定すべきパラメータ数を調整することで、パラメータａｋを必ず可同定（一意に同定可能）にすることができる。
【0089】
例えばｕｋの値が一定値である場合、uk(t-L)=uk(t-L-1)=…=uk(t-L-n+1)となるが、このような場合はｎ＝１として一つのパラメータａ１だけを同定する様にすればパラメータの可同定性を維持することができる。一般的には、所与のｕｋを用いて上記の持続的励振条件を調べることで同定可能なパラメータ数ｎの上限値ｎｍａｘを求めることができるので、ｎ≦ｎｍａｘの範囲内で、例えば、ＡＩＣ（赤池情報量規範）やＢＩＣ（ベイズ情報量規範）、ＭＤＬ基準（最小符号化基準）などの各種の規範および基準を用いたり、交差検証（クロスバリデーション）を行ったりすることで、ｎの適正値を決めて、（５）式の可同定性を常に維持するこができる。
【0090】
ペアワイズで同定することの第二の利点は、たとえ、出力変数を除く入力変数が１つの場合であっても、この方法では出力変数を入力とする自己回帰ＡＲモデルと入力変数（出力変数ではない入力変数）を入力とするＦＩＲモデルとの二つのモデルを分離して、同定することによる効果である。このような場合、ペアワイズでない通常の方法では、ＡＲＸモデルなどを用いて、自己回帰成分（ＡＲ成分）と重み付き移動平均成分（ＦＩＲ成分／ＭＡ成分）とを同時に同定することになるが、多くの場合、出力変数と入力変数との関係（相関）よりも出力変数自身の相関（自己相関）の影響の方が強い場合が多く、予測に対する自己回帰成分の影響が極めて強くなる場合が多い。このような場合、入力変数の値が急変した場合にその変化に追従できなくなり、予測が遅れるという現象が生じやすい。
【0091】
本実施形態のモジュラー型時系列データ予測装置の場合、例えば、出力変数はポンプ井流入量であり、入力変数として１か所の地上雨量を仮定すると、流入量に対する自己回帰成分の影響の方が地上雨量に対するＦＩＲ成分の影響よりも強くなり、降雨の変化により地上雨量の値が急変した場合でも流入量に変化が現れるまでは予測値の変化が小さく、予測値が実測値に対して遅れるという現象が生じやすい。これに対し、流入量に対する自己回帰モデル（ＡＲモデル）と地上雨量に対する有限インパルス応答モデル（ＦＩＲモデル）とを各々個別に同定するペアワイズの方法では、監視員等のユーザは、予測に対するＡＲ成分の影響とＦＩＲ成分の影響とを調整することができる。
【0092】
ペアワイズで同定することの第三の利点は、第二の利点と密接に関係するが、各入力変数に対する（１）式のモデルの意味が明確であることである。ここで、先と同じ様に、出力変数を除く入力変数が１つの場合である下記（２）式の単純なＡＲＸモデルを考えてみる。
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) （６）
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) （６）
【0093】
この時、上記（６）式をＡＲＸモデルとして同定すると、以下のｙに関する項（Ｙと定義する）とｕに関する項（Ｕと定義する）は、ＹとＵの意味は明確ではなく、ＹとＵとを個別に意味付けして説明することは困難である。
【0094】
Y：=a1×y(t-1)+ a2×y(t-2)+…+any×y(t-ny) （７）
U：＝b1×u(t-L)+ b2×u(t-L-1)+…+bnu×u(t-L-nu+1) （８）
一方、（７）式と（８）式とは、各々（５）式の形をしているので、ペアワイズで同定を行うと、ＹとＵとは各々ｙの予測値を意味することになるので、ＹとＵとは、各々別の説明変数（入力変数）に対する出力ｙの予測値であるという解釈をすることができる。
U：＝b1×u(t-L)+ b2×u(t-L-1)+…+bnu×u(t-L-nu+1) （８）
一方、（７）式と（８）式とは、各々（５）式の形をしているので、ペアワイズで同定を行うと、ＹとＵとは各々ｙの予測値を意味することになるので、ＹとＵとは、各々別の説明変数（入力変数）に対する出力ｙの予測値であるという解釈をすることができる。
【0095】
このようにペアワイズの同定を行うと、（７）式と（８）式とから、（６）式の形の予測モデルを構成的（Constructive）に構築することは容易であり、例えば（７）式のＵと（８）式のＹとの平均（＝（Ｕ＋Ｙ）÷２）として構成することで、（６）式の形の予測モデルを得る事ができ、これは、再度ｙの予測値を意味していることは明らかである。
【0096】
上記のように、予測モデルの形（構造）が同じであっても同定の手順を変える事で、その意味づけは変化する。ペアワイズで同定することの利点は、意味付けのしやすいモデルを構築することができる点である。
以上が多入力１出力ではなく１入力１出力のペアワイズで予測モデルを同定することの利点である。加えて、ペアワイズ予測モデルとして、（５）式のＦＩＲもしくはＡＲの形の回帰式の形で同定することで、ディジタル信号処理やシステム同定分野で従来から良く知られている周波数解析や安定解析などの様々な解析手法を、ペアワイズに直接流用できるため、その分野の知識があれば、ペアワイズでより詳細な解析や調整も可能になる。
以上が多入力１出力ではなく１入力１出力のペアワイズで予測モデルを同定することの利点である。加えて、ペアワイズ予測モデルとして、（５）式のＦＩＲもしくはＡＲの形の回帰式の形で同定することで、ディジタル信号処理やシステム同定分野で従来から良く知られている周波数解析や安定解析などの様々な解析手法を、ペアワイズに直接流用できるため、その分野の知識があれば、ペアワイズでより詳細な解析や調整も可能になる。
【0097】
ペアワイズ予測モデルとして、（５）式の様なＡＲモデル／ＦＩＲモデルをさらに分解し、時間遅れを考慮した変数を各々一つの説明変数と見なした単回帰モデルに分解して予測モデルを構築することもできる。
y(t)=a×uk(t-L)+c （９）
上記（９）式の様な単純な形にすると、（９）式に含まれるパラメータは、遅れ時間Ｌ、入力変数ｕｋに対する出力変数の倍率（ゲイン、比例係数）ａ、入力変数と出力変数の平均値の差（入力が０の場合の出力の値、バイアス）ｃの三つだけとなり、全てのパラメータの意味が明確になる。
y(t)=a×uk(t-L)+c （９）
上記（９）式の様な単純な形にすると、（９）式に含まれるパラメータは、遅れ時間Ｌ、入力変数ｕｋに対する出力変数の倍率（ゲイン、比例係数）ａ、入力変数と出力変数の平均値の差（入力が０の場合の出力の値、バイアス）ｃの三つだけとなり、全てのパラメータの意味が明確になる。
【0098】
（５）式の形では、遅れ時間Ｌとバイアスｃの解釈は（９）式と同様であるが、ａｋ、ｋ＝１、２、…、ｎという回帰係数（パラメータ）の個々の意味までは明確でなく、システム同定やディジタル信号処理に関する知識を用いて、（５）式自身を解釈することはできるものの、個々のａｋの直感的な解釈は困難である。これに対し、（９）式の単回帰式の場合はａの解釈も容易であり、単純に、平均値を０に処理した場合の入力変数に対する出力変数の倍率を示している。これにより、全てのパラメータの意味を容易に解釈できる様になるため、（９）式を用いると、何等かのアルゴリズムを適用した場合のパラメータ同定結果が直感的な感覚とずれるような場合（例えば、異常データが混入している場合などにそのような事が起こり得る）、例えば予測モデルの構築者が直感的に手動でパラメータを調整することも容易になり、また、予測モデルの中身を他者に説明することも極めて容易になる。
【0099】
本実施形態のモジュラー型時系列データ予測装置において、例えば、出力変数はポンプ井流入量であり、入力変数として１か所の地上雨量であると仮定すると、（９）式のｃは地上雨量がゼロの場合の流入量を意味し、これは、汚水と雨水が同一管渠を通る合流式下水の場合は、汚水量に相当するものであるという解釈ができる。そして、（９）式のａは流入量から汚水量を引き去った場合に、地上雨量計による雨量に対する雨水流入量の倍率を示しており、地上雨量と流入量との相関が十分に高ければ、このような単純な倍率も意味を持つ。さらに、（９）式の遅れ時間Ｌは、地上雨量の観測点から雨水ポンプ井までの雨水の流れの遅れ時間を意味し、下水管に流出するまでの流出の遅れと下水管路を流れて雨水ポンプ井に到達するまでの流下時間の和を意味すると解釈することができる。
【0100】
このように、解釈や説明が極めて容易になることが（９）式の様な単純なモデルを用いることの利点であり、実際に（９）式の形だけを用いた最も簡単なモデルは、流入量予測の簡易モデルとして、特許文献１の中でも採用されている。
【0101】
上記のように直感的な解釈が可能になることにより、手動でのパラメータ調整が容易になる具体的な例について説明する。
図３および図４は、幹線水位データと流入量データとの関係の一例を概略的に示す散布図である。
図３および図４は、幹線水位データと流入量データとの関係の一例を概略的に示す散布図である。
【0102】
図３および図４において、横軸は、幹線水位計１３１～１３Ｋのいずれかに相当するある箇所の幹線水位データ、縦軸は流入量計１１に相当する予測対象となる流入量データとした散布図の一例を示している。
図３に示す直線は、最小２乗法によって、単回帰のゲインとバイアスを求めたものであり、この場合ゲインａ＝９．６０２４、バイアスｃ＝２３．３１４２となっている。これは、予測値の平均２乗誤差（ＭＳＥもしくはＲＭＳＥ）が最小になるという意味で、最適な回帰直線である。
図３に示す直線は、最小２乗法によって、単回帰のゲインとバイアスを求めたものであり、この場合ゲインａ＝９．６０２４、バイアスｃ＝２３．３１４２となっている。これは、予測値の平均２乗誤差（ＭＳＥもしくはＲＭＳＥ）が最小になるという意味で、最適な回帰直線である。
【0103】
図３において、例えば幹線水位の値が－１．５であるとき、回帰直線に従うと流入量の値はおおよそ９程度と予測されることになる。しかしながら、図３の散布図によれば幹線水位の値が－１．５であるときの流入量は、９よりも大きい値である場合が圧倒的に多いことがわかる。
【0104】
上記のように、回帰直線による予測値と測定値との間にズレが生じる理由は、幹線水位が－２．５乃至－２付近のデータ数が圧倒的に多いため、この付近のデータにより適合する様に係数が推定されるためである。この回帰直線は、理論上は正しいものであるが、実際の流入量予測においては、必ずしも好ましくない場合が多い。なぜなら、流入量の予測値は、雨水排水ポンプの制御や運転支援に用いられるため、浸水リスクを回避するためには、流入量が多くなる場合をできる限り正確に予測したいという暗黙の要請があるためであり、流入量が少ない場合の予測精度はそれほど重要にならない場合が多いためである。
【0105】
しかし、このような暗黙の要請をアルゴリズムに組みこむ場合には、流入量が多い箇所の予測誤差に重みをつけるなどの処置が必要であり、その重みの調整なども含め、一般にかなりの労力を要する。
【0106】
図４では、最小２乗法と、最小２乗法の他の１０種類の様々な回帰手法とを、図３と同じデータに適用した結果得られた回帰直線Ｌ１－Ｌ１１の一例を示している。この中で、１番乃至１０番の方法は、Pythonと呼ばれるソフトウェアのパッケージに組み込まれている様々なアドバンストな回帰手法であるが、いずれの方法による回帰直線Ｌ１－Ｌ１０も、流入量が多いときに適切な予測値が得られるものではなかった。
【0107】
１１番目の方法は、幹線水位が－２．２乃至－１．５までの範囲のデータだけを抽出し、他のデータを使用せずに、通常の最小２乗法を適用した方法である。これは、本願発明者らが、幹線水位－２．５乃至－２付近のデータに適合するように回帰直線が求められることを理解した上で、その付近のデータを意図的に捨てることによって、流入量が多い場合のデータに適合するように意図的に調整した方法である。
【0108】
このように、アルゴリズムに頼った調整を行おうとすると、例えばデータの範囲を限定する、あるいは、データを間引くなどの様々な工夫を施さないと、本来望むパラメータが得られないことがある。しかし、１１番目の方法により得られる直線Ｌ１１は、２点を決めれば直線が唯一に定まるという極めて単純な原理を用いれば、人間が目視で直線Ｌ１１のようなラインを引くことにはほとんど労力を要しない。上記のように、単回帰という極めてわかりやすく、また図３や図４に示す散布図上に明確に可視化できるパラメータを用いることによって、手動での調整が極めて容易にできることがわかる。
すなわち、単回帰モデルを用いると、倍率（ゲイン）とバイアスという容易に理解できるパラメータを用いて解釈できて、調整が容易になる。
すなわち、単回帰モデルを用いると、倍率（ゲイン）とバイアスという容易に理解できるパラメータを用いて解釈できて、調整が容易になる。
【0109】
さらに、単回帰を用いると別の観点で解釈することも可能になる。すなわち、統計分野における基本的な統計量である、平均、標準偏差、および、相関、という三つのパラメータを用いて解釈することが可能になる。
【0110】
単回帰モデルの回帰係数（比例係数、傾き、倍率、ゲイン）ａは、相関係数ｒと次式の関係にあることが広く知られている。
a=r×(σy/σu) （１０）
ここで、ｒ、σｙ、σｕは、各々、入力変数と出力変数の相関係数、出力変数の標準偏差、入力変数の標準偏差、である。また、回帰係数ａとバイアスｃは、単回帰モデルの予測２乗誤差が最小となる様に最小２乗法で同定した場合、以下の関係になることも広く知られている。
c=μy-a×μu （１１）
ここで、μｙ、μｕは、各々、出力変数の平均と入力変数の平均、を表す。
a=r×(σy/σu) （１０）
ここで、ｒ、σｙ、σｕは、各々、入力変数と出力変数の相関係数、出力変数の標準偏差、入力変数の標準偏差、である。また、回帰係数ａとバイアスｃは、単回帰モデルの予測２乗誤差が最小となる様に最小２乗法で同定した場合、以下の関係になることも広く知られている。
c=μy-a×μu （１１）
ここで、μｙ、μｕは、各々、出力変数の平均と入力変数の平均、を表す。
【0111】
従って、（９）式において、入力変数と出力変数を各々、平均と標準偏差を用いて正規化しておけば、回帰係数ａと相関係数ｒの値は一致し、ｃ＝０となる。
この関係が意味することは、単回帰モデルのゲインとバイアスでモデルを説明することと、各変数（入力変数と出力変数）の平均、標準偏差と、入力と出力との相関係数、のみでモデルを説明することは等価であるということである。従って、予め、変数毎に適切な平均と標準偏差の推定値を用いて各変数のデータを正規化しておけば、単回帰の入出力関係は、次式のように、相関係数のみで、全てを説明できることになる。
y´(t)=a×uk´(t-L)=r×uk´(t-L) （１２）
ここで、y´=(y-μy)/σy、uk´=(uk-μuk)/σukで定義される正規化された出力変数と入力変数である。
この関係が意味することは、単回帰モデルのゲインとバイアスでモデルを説明することと、各変数（入力変数と出力変数）の平均、標準偏差と、入力と出力との相関係数、のみでモデルを説明することは等価であるということである。従って、予め、変数毎に適切な平均と標準偏差の推定値を用いて各変数のデータを正規化しておけば、単回帰の入出力関係は、次式のように、相関係数のみで、全てを説明できることになる。
y´(t)=a×uk´(t-L)=r×uk´(t-L) （１２）
ここで、y´=(y-μy)/σy、uk´=(uk-μuk)/σukで定義される正規化された出力変数と入力変数である。
【0112】
このように解釈すると、例えば、入出力の相関が無い場合、すなわちｒ＝０となる様な入出力関係を持つ場合は、その出力に対して対応する入力は全く予測能力を持たない事が（１１）式の関係から明確に理解することができる。また、現実のデータは、様々なノイズやアウトライア（外れ値）で汚染されている（コンタミされている）場合も多いが、このような時、平均、標準偏差、相関係数、の三つのパラメータを外れ値に対してロバストに推定するロバスト推定法を用いて推定することにより、外れ値に対してロバストな予測モデルを構築することもできる。もちろん、より複雑なモデルのパラメータ推定に対する様々なロバスト推定法も開発されているが、一般にロバスト推定した結果を解釈することは容易ではないのに対し、平均、標準偏差、相関係数という三つのパラメータに対するロバスト推定の場合は、何等かのロバスト推定を適用した推定結果の良否の判断や解釈も容易になるため、ノイズやアウトライアにコンタミされたデータに対する予測モデル構築を行いたい場合の解釈や調整も容易になる。
【0113】
このように、（９）式や（１０）～(１２)式は、容易に理解しやすいパラメータのみで構成されているが、（５）式と比較しても、極めて単純な構造をしており、実際の応用においては、その表現能力が十分でなく、十分な予測精度が得られない可能性があった。そこで、（９）式（あるいは（１０）～（１２）式）の形から、(５)式の形を合成することを検討する。（９）式のモデルを用いて(５)式の形を得るためには、(９)式の遅れ時間Ｌを可変にして足し合わせればよい。すなわち、まず、（９）式における遅れ時間Ｌ（（５）式のＬと区別するため以下ではＬ´とする）の最小の値を（５）式のＬと対応させてＬ´＝Ｌとし、Ｌ以上の遅れ時間について、Ｌ´＝Ｌ＋１、Ｌ´＝Ｌ＋２、…Ｌ´＝Ｌ＋ｎ－１としたｎ個の（９）式のモデルを加え合わせれば（５）式の形のモデルが得られる。
【0114】
すなわち、以下の（１３）式のｎ個の単回帰モデルを加え合わせれば、形式的に（５）式と等価な（１４）式のＦＩＲもしくはＡＲモデルの形のモデルが得られる。
y(t)=a1×uk(t-L)+c （１３＿１）
y(t)=a2×uk(t-L-1)+c （１３＿２）
y(t)=an×uk(t-L-n+1)+c （１３＿ｎ）
y(t)=1/n×(a1×uk(t-L)+a2×uk(t-L-1)+…+an×uk(t-L-n+1))+c （１４）
y(t)=a1×uk(t-L)+c （１３＿１）
y(t)=a2×uk(t-L-1)+c （１３＿２）
y(t)=an×uk(t-L-n+1)+c （１３＿ｎ）
y(t)=1/n×(a1×uk(t-L)+a2×uk(t-L-1)+…+an×uk(t-L-n+1))+c （１４）
【0115】
異なるｎ個の遅れ時間のｎ個の単回帰モデルを作成した上で、その平均モデルを作成すれば、形式的には（１）式と同じＦＩＲあるいはＡＲモデルの形のモデルが得られる。この際、上記の様に単純に平均化処理を行うこともできるが、先の述べた単回帰のゲインとバイアスと相関係数との関係を用いて、以下の様に重み付き平均を行うこともできる。
【0116】
まず、予め、入出力変数を正規化して、（９）式を（１２）式の形で表しておく。すると、回帰係数は相関係数と一致するので、この回帰係数（の絶対値）で重みづけした重み付き平均値として（１４）式に類似した（１５）式のＦＩＲ／ＡＲモデルの形の式を得ることができる（以下では、正規化した変数と正規化していない変数との記号を区別せず、文脈に応じてｕ（ｔ）、ｙ（ｔ）は正規化した変数を表すこととする。）
y(t)=(a1×uk(t-L)+…+an×uk(t-L-n+1))/(|a1|+|a2|+…+|an|)) （１５）
y(t)=(a1×uk(t-L)+…+an×uk(t-L-n+1))/(|a1|+|a2|+…+|an|)) （１５）
【0117】
このようにすると、相関が強くなる遅れ時間による予測に対して重みをつけたペアワイズ予測モデルの合成が可能になる。なお、回帰係数と相関係数の関係が（１２）式の様に陽に関係つけられない場合でも相関係数絶対値もしくは回帰係数の絶対値で重みづけ平均化処理を行うことは可能であるが、（１２）式の関係があることで、重みづけを行う事の意味がより明確で説得性の高いものとなる。また、（１５）式において、相関係数の絶対値ではなく相関係数の２乗によって重みづけを行っても良い。
【0118】
ペアワイズ予測モデルとして、上記の二つの合成法は、いずれの方法でも入出力間の線形の関係しか表現することができないが、入出力間に非線形の関係がある場合も多い。このような場合は、（５）式の代わりに、適当な非線形関数φ（・）を用いて、以下の（１６）式の様な非線形回帰の形として非線形変換すればよい。
【0119】
y(t)=a1×φ（uk(t-L)）+…+an×φ（uk(t-L-n+1)）+c （１６）
これは、機械学習の分野で広く知られている非線形回帰のテクニックであり、パラメータａｋ、ｋ＝１、２、…、ｎに関する線形性さえ維持していれば、線形回帰の様々な手法を直接適用することができる。また、機械学習の分野でよく知られている様に、非線形関数φ（・）を直接指定しなくても、（１６）式を計算する際に必要となる計画行列と呼ばれる行列に対して、φ（・）を指定することと等価に変換できるカーネル関数（類似度関数）を直接指定しても良い。ただし、その場合には、直感的な解釈性が若干低下する可能性がある事には注意する必要がある。
これは、機械学習の分野で広く知られている非線形回帰のテクニックであり、パラメータａｋ、ｋ＝１、２、…、ｎに関する線形性さえ維持していれば、線形回帰の様々な手法を直接適用することができる。また、機械学習の分野でよく知られている様に、非線形関数φ（・）を直接指定しなくても、（１６）式を計算する際に必要となる計画行列と呼ばれる行列に対して、φ（・）を指定することと等価に変換できるカーネル関数（類似度関数）を直接指定しても良い。ただし、その場合には、直感的な解釈性が若干低下する可能性がある事には注意する必要がある。
【0120】
（１６）式の様な非線形回帰を用いる事で、ある入力変数と出力変数との間に非線形関係がある場合にも対応することが可能になる。なお、ペアワイズ予測モデルの同定では、入力変数毎に出力変数との回帰モデルを構築するので、ある入力変数と出力変数の関係は非線形であるが、別の入力変数と出力変数との関係は線形であるような場合には、入力変数毎に（５）式と（１６）式とを単に使い分けるだけでよく、これによって、一般的には説明性も向上することが期待される。以上の一連の作用が、本実施形態におけるペアワイズ予測モデル同定部４の作用である。
【0121】
次に予測モデル合成法定義部５では、ペアワイズ予測モデル同定部４で定義した予測モデルの合成法を定義する。
最も簡単な予測モデルの合成法は、ペアワイズ予測モデルで各入力変数に対して出力変数を予測するモデルが構築されているため、その平均により予測モデルの合成を行う方法である。この場合は、予測モデル合成法定義部５では、ペアワイズ予測モデル同定部４のｐ個の予測出力（以下では、各入力変数ｕ１、ｕ２、…、ｕｐに対する予測出力をｙ１、ｙ２、…、ｙｐとする。）を平均化する次式を定義することになる。
y(t)=mean(y1(t)、y2(t)、…、yp(t))=1/p×(y1(t)+y2(t)+…+yp(t)) （１７）
（１７）式の定義は、最も基本的な定義方法であるが、この定義を改良することにより、合成した予測出力に対する信頼性を向上させたり、意図的に予測に傾向（バイアス）を持たせた予測を行ったりすることが可能になる。これを以下に順に説明する。
最も簡単な予測モデルの合成法は、ペアワイズ予測モデルで各入力変数に対して出力変数を予測するモデルが構築されているため、その平均により予測モデルの合成を行う方法である。この場合は、予測モデル合成法定義部５では、ペアワイズ予測モデル同定部４のｐ個の予測出力（以下では、各入力変数ｕ１、ｕ２、…、ｕｐに対する予測出力をｙ１、ｙ２、…、ｙｐとする。）を平均化する次式を定義することになる。
y(t)=mean(y1(t)、y2(t)、…、yp(t))=1/p×(y1(t)+y2(t)+…+yp(t)) （１７）
（１７）式の定義は、最も基本的な定義方法であるが、この定義を改良することにより、合成した予測出力に対する信頼性を向上させたり、意図的に予測に傾向（バイアス）を持たせた予測を行ったりすることが可能になる。これを以下に順に説明する。
【0122】
まず、（１７）式の様な単純な平均化処理（標本平均）を行うと、各ペアワイズ予測モデルの予測値が平等に扱われているため、予測精度の良いペアワイズ予測モデルの予測値と予測精度の悪いペアワイズ予測モデルの予測値が混合されて、合成した全体の予測精度が劣化してしまう可能性がある。また、実際の運用においては、ある入力変数の計測データの信頼性が低く、別の入力変数の計測データの信頼性が高いという場合も稀ではないが、信頼性の高い入力変数に対するペアワイズ予測モデルの予測値も信頼性の低い入力変数に対するペアワイズ予測モデルの予測値も平均化してしまうことで、合成した予測出力の予測精度が劣化してしまう。極端な場合、例えば、ある入力変数のセンサの故障や不具合などにより、当該入力変数の時系列データにアウトライア（外れ値、異常値）が多量に含まれる様になる場合などには、予測精度が劣化するだけでなくアウトライアに引きずられて合成した予測自身が無意味になる（破綻してしまう）可能性がある。
【0123】
このような状況に対応するために予測モデル合成法定義部５で、以下の様な複数のアプローチをすることができる。
一つ目は、主に後者のアウトライアに対して合成した予測値自身が破綻しない事を重視するアプローチであり、（１７）式の標本平均処理に変えて、ロバスト推定を採用する方法である。
一つ目は、主に後者のアウトライアに対して合成した予測値自身が破綻しない事を重視するアプローチであり、（１７）式の標本平均処理に変えて、ロバスト推定を採用する方法である。
【0124】
（１７）式で用いる「標本平均」という処理は、ｙ１、ｙ２、…、ｙｐのｐ個の予測値の中から代表値を推定する方法の一つであり、このような代表値の推定を統計分野では位置母数の推定と呼ぶ。すなわち、「標本平均」は、位置母数の推定方法の一つであり、統計的には、標本平均は推定効率の観点では良い性質を持つことが知られているが、一方でアウトライアに対するロバスト性の面では最も脆弱（非ロバスト）であることも知られている。従って、外れ値に対するロバスト性を向上させるためには、（１７）式をロバストな位置母数推定で置き換える方が良い。
【0125】
ロバスト推定の分野では外れ値に対するロバスト性を評価するいくつかの指標が知られているが、最も直感的に理解しやすい指標としてブレークダウンポイントという指標がある。これは、統計的推定に用いるデータの中に何パーセントアウトライアが混入することを許容するかという指標であり、用いるデータのＸ％を、仮想的なアウトライアを想定して∞に置き換えた場合、推定量（平均などの位置母数を推定する場合は位置母数の推定値）が∞になる（＝破綻する）かならないかの境界（∞になる直前）のＸの値をブレークダウンポイントと呼ぶ。
【0126】
例えば、「標本平均」という位置母数の推定に対しては、ただ１点のデータを∞に置換するだけでその標本平均値も∞になるので、「標本平均」のブレークダウンポイントは０％である。ロバスト統計ではブレークダウンポイントの最大値は５０％であることが知られており、最大のブレークダウンポイントを持つ位置母数推定量として「中央値（メジアン）」が知られている。従って、予測値が外れ値に影響されないようにすることを最大の目的とする場合には、（１７）式の「平均」を「中央値」で置換して、予測モデル合成法定義部５の定義とすることもできる。
【0127】
一方、「中央値」はロバストではあるが、推定効率が悪い事が知られており、直感的にもｐ個の変数で予測している中の一つの予測値しか用いないため、精度の高い予測を行うことを重視する場合に最良の方法ではないことは容易に推測できる。そのため、推定効率の向上とロバスト性の向上を両立するための各種のロバスト推定方法が知られている。最も単純な方法は、「トリム平均（刈込平均）」と呼ばれる処理であり、ｐ個のデータの上位と下位とのα％を削除した上で平均をとる方法である。直感的に明らかな様にαを大きくするとロバスト性は向上し、その極限として中央値推定があり、αを小さくして０に近づけるとロバスト性が低下しその極限が通常の平均（標本平均）であることは明らかである。従って、このトリム平均を予測モデル合成法定義部５の定義とすることもできる。この場合、αの調整によってロバスト性を調整できるが、適切に調整する事自身が難しい場合も考えられるため、別のロバスト推定方法をとることもできる。
【0128】
このような方法の代表的な例として、ｐ個のｙ１、ｙ２、…、ｙｐの中からすべての組み合わせの（ｐ（ｐ－１）／２個）の、二つのｙｉとｙｊ（ｉ≠ｊ）との平均を計算した上で、その中央値を採用するホッジスレーマン推定量（ＨＬ推定）と呼ばれる位置母数推定を行うこともできる。この方法は、ブレークダウンポイントが約３０％でロバスト性が高いうえに推定効率も良いことが知られており、このＨＬ推定を予測モデル合成法定義部５の定義とすることもできる。
【0129】
他のロバスト推定法として、一般には多変量データに対して用いられる手法であるＭＣＤと呼ばれる手法が知られているがこれを適用することもできる。ＭＣＤは、ｐ個のデータの所定の割合（通常５０％～７５％）のデータを取り出し、そのすべての組み合わせの中で分散が最小になるデータを用いて推定を行う手法であり、この方法を用いて平均を推定することもできる。ＭＣＤもロバスト性が高く推定効率の良い方法として知られている。一方、ＨＬ推定やＭＣＤ推定は、ｐ個のデータの中から取り出すデータの数（ＨＬ推定場合２個、ＭＣＤの場合はｐ／２～３ｐ／４個程度）の全ての組み合わせに対する計算が必要となるため、ペアワイズモデルの入力変数の数ｐが増加すると処理時間が飛躍的に増加するため、予測モデル合成法定義部５の定義として用いても、実際に予測を行う際にリアルタイム性を確保できない可能性がある。
【0130】
このような場合に対応するため、「全ての組み合わせ」に対して処理を行うのではなく、リアルタイムで処理が可能な回数だけ繰り返し処理を行う様にブートストラップ法を用いて代用しても良い。すなわち、ｐ個のｙ１、ｙ２、…、ｙｐの中から、所定のｋ個のデータをランダムに繰り返し抽出し、所定の回数（Ｎ回）、その平均値を求める。繰り返し求めたＮ個の平均値の中から、ＨＬ推定と同じように、例えばＮ個の平均値の中央値を採用する。このようなブートストラップ法による平均値推定を、予測モデル合成法定義部５の定義とすることで、リアルタイム性を維持しながら、推定効率が良くロバスト性の高い予測値の合成を行うことができると考えられる。
上記以外にも外れ値に重みを付けて推定を行うM推定などのロバスト推定を位置母数推定に適用した方法を予測モデル合成法定義部５の定義とすることもできる。
上記以外にも外れ値に重みを付けて推定を行うM推定などのロバスト推定を位置母数推定に適用した方法を予測モデル合成法定義部５の定義とすることもできる。
【0131】
二つ目は、先のアプローチが主に外れ値（外れた予測値）に対してロバストに予測値の代表値を推定することを目的としていたのに対し、実際に予測精度の良い予測値に重みを付けて合成した予測を行うことを目的としたアプローチである。基本的な考え方は、（１７）式の単純な標本平均に変えて、重み付き平均を行う方法であり、次式で表される式で合成予測出力を定義する。
y(t)=w1×y1(t)+w2×y2(t)+…+wp×yp(t) （１８）
ここで、ｗｋ、ｋ＝１、２、…、ｐは、重みであり、ｗ１＋ｗ２＋…＋ｗｐ＝１となる制約を満たす。予測モデル合成法定義部５の定義では、この重みの設定法を定義する必要がある。
y(t)=w1×y1(t)+w2×y2(t)+…+wp×yp(t) （１８）
ここで、ｗｋ、ｋ＝１、２、…、ｐは、重みであり、ｗ１＋ｗ２＋…＋ｗｐ＝１となる制約を満たす。予測モデル合成法定義部５の定義では、この重みの設定法を定義する必要がある。
【0132】
これには、ペアワイズ予測モデル同定部４で、ペアワイズ予測モデルを同定した際の予測出力と実際の出力の間の（重）相関係数、あるいは、その２乗である決定係数を用いることができる。これにより、少なくともペアワイズの予測モデル同定時のデータに対して、予測精度の高い予測能力を持つモデルの予測出力を重視した重み付き平均で予測出力を合成することができる。
【0133】
また、ペアワイズの予測モデルを、例えば上述の（１５）式のように構築した場合には、各遅れ時間毎の説明変数と出力変数の間の相関係数が既に算出されているので、統合したペアワイズ予測モデル（（１４）式に相当）毎に、各相関係数の絶対値の平均値ｒｋ、ｋ＝１、２、３、…、ｐを求め、ｒｋに応じて、重みｗｋをｗｋ＝ｒｋ／（ｒ１＋ｒ２＋…＋ｒｐ）と定義することもできる。
【0134】
このように、ペアワイズ予測モデルの予測能力に応じて重み付き平均値の重みを決定して（１８）式により、予測モデルの合成法を定義する動作が、本実施形態における予測モデル合成法定義部５の動作の一例である。
【0135】
三つ目は、二つ目と同様に（１７）式の、各ペアワイズ予測モデルの重み付き平均値により、合成した予測出力を定義するが、この重みを、データを用いて同定（推定）することで決定するものである。この際、ペアワイズ予測モデルの同定時に用いたデータと同定したパラメータ値を用いると、その同定データに対する各ペアワイズ予測モデルの予測値ｙ１、ｙ２、…、ｙｐの時系列データが得られる。
【0136】
そして、（１７）式を、この予測値ｙ１、ｙ２、…、ｙｐを入力として、同定に用いた実際の出力データを出力ｙと見なすと、（１７）式自身が重みｗｋ、ｋ＝１、２、…、ｐを回帰係数とする重回帰モデルの形式になっている。従って、これらの予測時系列データと同定に用いた出力データとを用いて、重みｗｋ、ｋ＝１、２、…、ｐを重回帰の方法で同定することができる。ただし、重みは各々正でなければならないという不等式制約とｗ１＋ｗ２＋…＋ｗｐ＝１という等式制約を満たす必要があるため、通常の重回帰で用いる最小２乗法は適用できないが、混合線形推定法などの制約条件を考慮できる線形回帰の推定法を用いて、重みを推定することが可能になる。
【0137】
なお、この方法は、（１３）式から（１４）式を得る代わりに（１５）式を得る箇所にも同様の考え方を適用することができるので、単回帰モデル⇒ＦＩＲ／ＡＲモデル⇒多入力の伝達関数モデル、という３段階の段階的な学習（推定）によって予測モデルを合成することができる。
【0138】
本実施形態において、上記のようにして重みｗｋ、ｋ＝１、２、…、ｐを推定した（１７）式によって予測モデルの合成法を定義する方法が、予測モデル合成法定義部５の動作の他の例である。
【0139】
四つ目は、一つ目から三つ目とは異なり、合成した予測に傾向（バイアス）を持たせたい場合のアプローチであり、意図的に過大、あるいは、過小の予測を行う事を目的とする場合の予測出力の合成法である。
【0140】
上記の予測を行いたい状況は、現実の問題ではしばしば遭遇する。例えば、本実施形態で行う雨水の流入量予測の場合、この予測情報は雨水排水ポンプの起動や停止のタイミングを図るための支援情報として利用される場合が多い。このような場合、実際に流入する雨水流入量より過小な流入量の予測を行ってしまうことはリスク回避の観点から致命的になる事がある。
このような場合、もちろん、正確な流入量を予測することが好ましい事は言うまでもないが、それが現実的に困難な場合は、多少多めの量を予測しておく方が安全である。
このような場合、もちろん、正確な流入量を予測することが好ましい事は言うまでもないが、それが現実的に困難な場合は、多少多めの量を予測しておく方が安全である。
【0141】
また、例えば雨水排水ポンプの制御において流入量予測を利用する際には、ポンプの起動タイミングを図るために過大側の予測をすることが好ましいが、停止タイミングを図るためには過小側の予測をすることが好ましいため、過小側の予測を行いたい場合もある。
【0142】
また、例えば、水や電力の需要予測などの場合は、あまり過小に予測を行ってしまうと、水や電力の供給に支障をきたす可能性が考えられるため、若干過大側に予測を行っておく方が安全な場合もある。
【0143】
また、このように、過大側あるいは過小側の予測を行える様にしておくと、純粋に過大な予測や過小の予測を行いたいという場合以外にも、実際に予測システムを運用して長期評価した場合に、予測値が過小に出る傾向がある場合には過大側に再調整したり、逆に予測値が過大に出る場合には過小側に再調整したりして微調整を行うことも可能になる。
【0144】
このような動機に基づいて行う予測出力の合成法を以下に示す。これを行うためには、各ペアワイズ予測モデルの出力予測値ｙ１、ｙ２、…、ｙｐに対して、先に述べた位置母数の推定方法を定義する以外に、標準偏差や分散に対応する尺度母数の推定方法も同時に定義する。例えば、位置母数を平均とした場合には、尺度母数として、標準偏差を採用する様に定義する。また、位置母数としてロバスト推定法のメジアン（中央値）を用いた場合には、尺度母数としてＭＡＤ（中央値絶対偏差：Median Absolute Deviation）を採用する。さらに、ＨＬ推定、ＭＡＤ推定、ブートストラップ推定、Ｍ推定などを用いた場合には、各々の位置母数に対する尺度母数の推定方法を定義する。このようにして、ｙ１、ｙ２、…、ｙｐに対して適用する位置母数と尺度母数の推定法を予め定義しておく。これを以下では各々μとσと表す（通常μは平均値、σは標準偏差の意味で使われることが多いが、ここでは、必ずしも平均と標準偏差ではなく、ここで定義した位置母数と尺度母数を表すものとする）。
【0145】
次に、各ペアワイズ予測モデルの予測出力の集合、ｙ１、ｙ２、…、ｙｐに対して、その集合のμ±ｋσの値を合成した予測出力とするように、ｋの値と符号とを指定する。例えば、ｋ＝１としてμ＋σを指定すると、ｙ１、ｙ２、…、ｙｐの中の位置母数から尺度母数の１倍だけ過大な値を合成した予測出力とすることを意味する。逆にｋ＝－１とすると、ｙ１、ｙ２、…、ｙｐの中の位置母数から尺度母数の１倍だけ過小な値を合成した予測出力とすることを意味する。もし、位置母数と尺度母数として平均と標準偏差を用いた場合は、ｙ１、ｙ２、…、ｙｐが正規分布に従うと仮定すると、μ＋σは約６８％の値（μが５０％の値）、μ－σは約３２％の値を取り出すことに対応する。
【0146】
なお、平均や標準偏差を用いなくても、ｙ１、ｙ２、…、ｙｐの中のＫ％にあたるＫ％分位点を取り出して合成出力とするという方法も、位置母数と尺度母数を適切に定義すれば、μ±ｋσの形式で書けるので、このように設定した分位点と直接抽出するという方法で指定することも可能である。例えば、過大側の予測を行いたい場合に予測の異常値をある程度除外する事を想定して９０％分位点にあたる予測値を採用し過大予測を行う様に合成した予測出力を定義したり、１０％分位点にあたる予測値を採用し過小予測を行う様に合成した予測出力を定義したりすることもできる。
【0147】
予測モデル合成法定義部５の定義方法として四つの方法を説明したが、予測モデル合成法定義部５は、これらのいずれか、あるいは、その組み合わせによって予想モデルの合成法と定義することができる。なお、以上の一連の作用は、オフラインで過去のデータを用いて定期的もしくは非定期的に実行される。
【0148】
次に、オフラインで同定および定義した定義式を用いて、オンラインで予測を行う。
オンラインの予測は、時間の進行方向における所定の周期ＴＨ（＜＜ＴＬ（オフライン同定を定期的に行う場合の周期））で行われる。
オンラインの予測は、時間の進行方向における所定の周期ＴＨ（＜＜ＴＬ（オフライン同定を定期的に行う場合の周期））で行われる。
【0149】
まず、オンライン予測用データ抽出部３２から周期ＴＨで、入力変数のデータを抽出する。本実施形態では、入力変数は、流量計（雨水ポンプ井流入量計）１１による雨水ポンプ井流入量と、幹線流量計１２による幹線流入量と、幹線水位計１３１～１３ＫによるＫ個の幹線水位計と、地上雨量計１４１～１４ＭによるＭ個の地上雨量と、Ｑ×Ｐメッシュ１５１１～１５ＱＰの各メッシュにおけるレーダ雨量と、を含む。
【0150】
次に、ペアワイズ出力変数予測部６では、ペアワイズ予測モデル同定部４で同定した各入力変数と出力変数のペアワイズ予測モデルを用いて、ｐ個の予測出力を計算する。これは、（５）式、（１４）式、（１５）式などの形で同定したパラメータ値が同定されたペアワイズ予測モデルにオンライン予測用データ抽出部３２で抽出したオンラインの入力変数データを入力することで直ちに計算できる。これがペアワイズ出力変数予測部６の作用である。
【0151】
次に、合成出力変数予測部７では、予測モデル合成法定義部５で定義した（１７）式や（１８）式に相当する合成法の定義式を用いて、ペアワイズ出力変数予測部６から出力される各入力変数に対するペアワイズの予測出力を入力することで合成した予測出力が計算される。
【0152】
次に、予測誤差評価部８は、合成出力変数予測部７で計算された予測出力の予測結果とペアワイズ出力変数予測部６のペアワイズの予測出力結果とを取得し、時間の進行方向における所定の周期ＴＨで時系列データとして保存する。ここで、保存された予測結果は、誤差評価を行うタイミングまで保持される。予測誤差の評価を行うタイミングは、予測モデルのユーザなどが外部から指示しても良いし、時間の進行方向における所定の周期ＴＭ（例えばＴＨ＜ＴＭ＜ＴＬ）で実施しても良い。いずれの方法であっても、予測誤差の評価を行うタイミングで、評価用データ抽出部３３で、予測誤差評価部８に保存された予測出力の時系列データと対応する時間の計測された出力変数、すなわち、本実施形態では雨水流入量の時系列データを抽出する。
【0153】
そして、予測誤差評価部８では、抽出された出力（雨水流入量）の時系列データとｐ個のペアワイズ予測モデルの予測出力との誤差、および、合成した予測出力の誤差を評価する。この際、誤差の評価は、通常時系列解析で行われるＭＳＥ（平方２乗誤差＝Ｌ２誤差）やその平方根であるＲＭＳＥなどであっても良いが、目的に応じて、誤差評価の評価基準を適宜設定しておいても良い。例えば、ピークの雨水流入量の予測精度が重要である場合はピーク流入量の差や誤差が最大になる場合の差（Ｌ∞誤差＝∞ノルム誤差）を誤差評価基準としても良いし、Nash-Sutcliffe係数などの評価基準を用いても良い。また、流入量が急激に増加する場合の遅れを改善する事を目的としたい場合には、予測出力と実績出力との位相差を計算し、位相差を小さくすることを誤差の評価基準としても良い。いずれにしても、予測誤差評価部８は、設定した誤差評価基準に基づいて、所定のタイミングもしくは周期ＴＭで、時系列データとｐ個のペアワイズ予測モデルの予測出力との誤差、および、合成した予測出力との誤差を評価する。
【0154】
次に、ペアワイズ予測モデル修正部９では、予測誤差評価部８で評価した予測誤差に基づいて、ペアワイズ予測モデルのパラメータの調整、あるいは、特定のペアワイズモデルの削除、もしくは、複数のペアワイズモデルの統合を行う。
【0155】
以下に、上記のような調整機能、削除機能、および、統合機能を付加したモジュラー型予測をより具体化した方法の例について説明する。
一つ目のパラメータ調整について、ペアワイズ予測モデル修正部９は、予測誤差評価部８で評価した予測誤差に基づいて、著しく予測精度が悪い（予測誤差が所定の許容誤差より大きい）と判断される入出力ペアを抽出し、抽出された入出力ペアのペアワイズ予測モデルの再同定を行う。予測精度の悪さの判定は、（１）予測誤差に対するしきい値（許容誤差）を直接指定してしきい値を超過したペアワイズ予測モデルを予測精度の悪いモデル（あるいは入出力ペア）と判断する、（２）合成した予測出力の誤差とペアワイズ予測モデルの各予測誤差とを比較し、合成した予測誤差に対して相対的に著しく誤差が大きい（予め設定された許容範囲を超える）予測誤差を持つペアワイズ予測モデルを予測精度の悪いモデルと判断する、（３）ペアワイズ予測モデルの予測誤差（＋合成した予測誤差も含んでも良い）に対する位置母数μｅと尺度母数σｅを計算し、μｅ±ｋσｅ（ｋは設定するパラメータで２～３ぐらいで通常は設定する）から外れた予測誤差を持つペアワイズ予測モデルを予測精度の悪いモデルと判断する、などの方法で行うことができる。（１）－（３）のような判断によって、パラメータの再調整（再同定）が必要と判断された場合、ペアワイズ予測モデル修正部９は、所定の期間の直近のデータを用いてパラメータの再同定を行う。
一つ目のパラメータ調整について、ペアワイズ予測モデル修正部９は、予測誤差評価部８で評価した予測誤差に基づいて、著しく予測精度が悪い（予測誤差が所定の許容誤差より大きい）と判断される入出力ペアを抽出し、抽出された入出力ペアのペアワイズ予測モデルの再同定を行う。予測精度の悪さの判定は、（１）予測誤差に対するしきい値（許容誤差）を直接指定してしきい値を超過したペアワイズ予測モデルを予測精度の悪いモデル（あるいは入出力ペア）と判断する、（２）合成した予測出力の誤差とペアワイズ予測モデルの各予測誤差とを比較し、合成した予測誤差に対して相対的に著しく誤差が大きい（予め設定された許容範囲を超える）予測誤差を持つペアワイズ予測モデルを予測精度の悪いモデルと判断する、（３）ペアワイズ予測モデルの予測誤差（＋合成した予測誤差も含んでも良い）に対する位置母数μｅと尺度母数σｅを計算し、μｅ±ｋσｅ（ｋは設定するパラメータで２～３ぐらいで通常は設定する）から外れた予測誤差を持つペアワイズ予測モデルを予測精度の悪いモデルと判断する、などの方法で行うことができる。（１）－（３）のような判断によって、パラメータの再調整（再同定）が必要と判断された場合、ペアワイズ予測モデル修正部９は、所定の期間の直近のデータを用いてパラメータの再同定を行う。
【0156】
本実施形態において、例えば、Ｍ個の地上雨量データの中のいずれかの地上雨量データと雨水流入量との間のペアワイズ予測モデルの精度が悪化したと判断された場合、ペアワイズ予測モデル修正部９は、地上雨量と雨水流入量との間のペアワイズ予測モデルのパラメータのみを再同定する。このように部分的な同定による再調整を行う事ができる点は、モジュラー型の構成（構成的なアプローチ）の大きな利点である。
【0157】
二つ目の特定のペアワイズモデルの削除は、一つ目のパラメータの調整と類似の手続きで行うことができる。基本的な手順としては、ペアワイズ予測モデル修正部９は、一つ目のパラメータの調整で述べた方法と類似の方法で誤差評価を行い、誤差が著しく大きくなり精度が劣化したペアワイズ予測モデルを抽出し、抽出したペアワイズモデルを削除する。誤差評価の方法は一つ目の方法と同じであるが、ペアワイズ予測モデル修正部９は、（１）劣化の程度を判断するしきい値を一つ目のものよりも大きくして、一つ目の調整よりも明らかに予測精度が劣化しているものを削除対象とする、（２）まず一つ目のパラメータの再調整を行い、パラメータの再調整を行っても精度劣化の判断に用いた基準をクリアできない場合に削除対象とする、などの方法で削除すべきペアワイズ予測モデルを抽出する。
【0158】
本実施形態においては、例えば、幹線水位計などは水位計が水没することなどにより、水位を正確に測ることができなくなりアウトライア（異常値）が多量に混入してしまうことなどが考えられる。この場合、水没した水位計によるペアワイズ予測モデルは分離されることになる。このように、予測精度に悪影響を与えている部分を部分的に切り離すことで予測精度の改善を図れることもモジュラー型の構成（構成的なアプローチ）の大きな利点である。
【0159】
３つ目のペアワイズ予測モデルの統合は、以下の様に実施される。まず、ペアワイズ予測モデル修正部９は、ペアワイズ予測モデルの予測誤差を相互に比較評価し、その予測誤差同士の差の絶対値が所定のしきい値以下となる相互に類似の予測精度を持つペアワイズ予測モデルを、相互に類似する予測モデルとして抽出しグループ化する。
【0160】
次に、ペアワイズ予測モデル修正部９は、このようにしてグループ化されたペアワイズ予測モデルの持つパラメータの値をグループ内で相互に比較する。例えば、係数パラメータａｋ、ｋ＝１、２、…、ｎをベクトル化してＡ＝［ａ１、ａ２、…、ａｎ］とし、相互に比較する予測モデルの係数パラメータベクトル同士の類似度を例えば内積などにより算出し、所定のしきい値により類似していると判断されたペアワイズ予測モデルを統合の対象候補とする。
【0161】
本実施形態のモジュラー型の構成的なアプロ―チでは、パラメータの可同定性（一意決定可能性）が基本的に担保されているので、予測結果が類似しておりパラメータ値も類似していれば、その入力変数はほぼ同じ値を持つはずであるが、念のため、ペアワイズ予測モデル修正部９は、候補として抽出されたペアワイズ予測モデルの入力変数同士の相関係数を求めて、所定のしきい値（例えば、相関係数０．９５）以上である場合に最終的に統合すべきペアワイズ予測モデルと判断する。そして、ペアワイズ予測モデル修正部９は、最終的に統合すべきと判断された複数のペアワイズ予測モデルの入力変数の平均をとった変数を、新たな一つの合成された入力変数として定義する。
【0162】
そして、ペアワイズ予測モデル修正部９は、（１）先に比較した係数ベクトルの平均値を採用する、あるいは、（２）新たに合成した平均化された入力変数と出力変数に対して再同定を行って推定する、のいずれかの方法を用いて、合成された入力変数に対するペアワイズ予測モデルの係数を決定する。
【0163】
上記のような統合は、本実施形態のモジュラー型時系列データ予測装置においては、例えば次の様なケースにおいて生じると考えられる。入力変数の中にＱ×Ｐメッシュ１５１１～１５ＱＰの各々のレーダ雨量のデータが含まれているケースにおいて、レーダ雨量のメッシュのサイズは、例えばXRAINと呼ばれるＸバンドの気象レーダの場合約２５０ｍ×２５０ｍであるため、隣接するメッシュの雨量データはほとんど同じ値をとる場合が多い。このような場合、相互に隣接するメッシュの雨量データは統合の対象となる可能性が高く、Ｑ×Ｐ個の入力変数が、少数の互いに独立した合成入力変数として統合されることになる可能性が高い。このような統合が起こると、レーダ雨量のメッシュ状のデータが、雨水流入量に与える影響の類似性に応じて統合されるため、入力変数の数を減らせるだけでなく、説明性が向上する事が期待できる。
【0164】
また、このような場合は、先に述べた多重共線性の問題が生じているケースであるが、これに対し、例えば、Lassoと呼ばれる方法などを用いて入力変数を自動的に選択すると、いくつかのメッシュの雨量データが代表的な入力変数として選択され、選択された入力変数に対応するメッシュの雨量と極めて類似するメッシュの雨量データは無視されることになる。雨量データが極めて類似しているので、このように代表メッシュの雨量を入力変数として選択しても通常は問題がない。しかしながら、近年頻発している局所的な豪雨が見られる場合、例えば、たまたま選択したメッシュの雨量はあまり多くないのに、入力変数として選択されなかった近くのメッシュの雨量が非常に大きい様な局所的な降雨が生じた場合、あるメッシュを代表メッシュとして選択してしまうと、流入量が増加しないという誤った予測を行う可能性がある。
【0165】
一方、ここで述べた入力変数の統合を行うと、類似の影響を与えるメッシュの雨量が平均化され、これは統合されたメッシュの平均雨量という意味を持つため、このようなケースでも妥当な予測結果を与える事が期待できると同時に、合理的な説明性も向上する。このように統合すべき入力変数をまとめて解釈し、説明しやすくするようにできる事もモジュラー型の構成（構成的なアプローチ）の大きな利点である。
【0166】
最後に、合成出力変数予測部７で算出（合成）された予測出力は、出力予測結果観測部１０に送られ、例えばいわゆるトレンドグラフと呼ばれる時系列データとしてユーザであるプラントの管理者や運転員に提示される。この際、リアルタタイム監視している監視の現在までの値は、実際の出力値である雨水流入量とその予測値を同時に表示し、将来の予測値は予測出力のみを表示する様にしておくことが好ましい。
【0167】
また、合成した予測出力（予測雨水流入量）だけでなく、ペアワイズ予測モデルのｐ個の予測出力を用いて、その位置母数μｙと尺度母数σｙを計算し、μｙ±Ｋσｙ（Ｋは設定パラメータで２～３程度で設定する）の範囲を同時表示したり、あるいは、ペアワイズ予測モデルの予測出力の中から最大値（あるいはロバスト性を考えて最大値に近い９５％点にもっとも近い予測値）や最小値（あるいはロバスト性を考えて最小値に近い５％点にもっとも近い予測値）を直接予測出力の範囲として同時に表示したりしても良い。
【0168】
本実施形態のモジュラー型時系列データ予測装置によれば、一般にはブラックボックスモデルと呼ばれるモデルに対して、その内部パラメータに対して物理法則に矛盾しない合理的な解釈が可能になる様な仕組みを導入することで、ホワイトボックスモデルに近い、モデルの解釈と合理的な説明性を向上させて、部分的な調整、削除、統合などのモデルの維持管理持管理（メンテナンス）の煩雑さを大幅に低減できるという効果が得られる点である。
【0169】
なお、このような効果を得るための根本的なアイデアは、本来、入出力データしか与えられていない内部構造が完全に未知である状況下で、出力に関連する可能性のある要素を要素毎に出力と直接関係づけ（ペアワイズ予測モデルに対応）、各要素と出力の関係を、最終的な入出力関係の内部構造として定義される様に、構成的に予測モデルを構築していることである。このようなアイデアを用いることで、本実施形態のモジュラー型時系列データ予測装置では、本来知ることのできない内部構造に立ち入ることなく、疑似的にホワイトボックスモデルと同様の説明性と調整の容易さを実現している。
【0170】
また、このモジュラー型時系列予測モデルの構築方法は、統計的機械学習分野でよく知られている。先に述べたバギングやブースティングなどのアンサンブル学習と呼ばれる手法とも類似している。バギングやブースティングは、弱学習器（弱学習モデル）と呼ばれる、簡単に構成できるが、精度が必ずしも高くないモデルをたくさん集めることで、精度の高い強学習器（強学習モデル）と呼ばれる予測モデルを構築する方法であり、バギングは弱学習器に対して多数決などの処理で強学習を行う方法、ブースティングは、弱学習器の結果を用いて強学習器を学習させる方法である。ただし、バギングやブースティングでは、学習(同定)に利用データを変えることでたくさんの弱学習器を生成しているが、本実施形態では、同定データを変えるのではなく、説明変数を変えることで、弱学習器に相当するペアワイズ予測モデルを同定している点で異なっている。なお、本実施形態のモジュラー型時系列データ予測装置において、合成出力の予測方法を定義する際に、ロバスト位置母数推定や予め指定した重みを用いた重み付き平均を用いるような方法はバギングの考え方に近く、重み付き平均の重みを再同定によって求める方法はブースティングの考え方に近い。
【0171】
以下、本実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの複数の実施例について説明する。
（第１実施例）
図５は、一実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの第１実施例について説明するための図である。
図５では、（１）式のＦＩＲモデルやＡＲモデルをペアワイズ予測モデルとした例を示している。
（第１実施例）
図５は、一実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの第１実施例について説明するための図である。
図５では、（１）式のＦＩＲモデルやＡＲモデルをペアワイズ予測モデルとした例を示している。
【0172】
本実施例のペアワイズ予測モデルＭＡは、回帰ブロックＡ１－ＡＮと、予測出力の合成ブロックＡ５と、を備えている。
例えば回帰ブロックＡ１－Ａ３のそれぞれは、入力変数に対して重み付き移動平均成分（ＦＩＲ成分／ＭＡ成分）を予測値として出力する。例えば回帰ブロックＡＮは、入力変数（過去の出力変数）に対して自己回帰成分（ＡＲ成分）を予測値として出力する。
例えば回帰ブロックＡ１－Ａ３のそれぞれは、入力変数に対して重み付き移動平均成分（ＦＩＲ成分／ＭＡ成分）を予測値として出力する。例えば回帰ブロックＡＮは、入力変数（過去の出力変数）に対して自己回帰成分（ＡＲ成分）を予測値として出力する。
【0173】
予測出力の合成ブロックＡ５は、回帰ブロックＡ１－ＡＮから出力された予測値を取得し、予測値を合成した合成予測値（Y(t+1)=F(Y1(t+1),Y2(t+1),…,Yn(t+1))）を演算して出力する。
例えば、通常のＡＲＸモデルなどによってパラメータを同定すると自己回帰成分（ＡＲ成分）の重みが大きくなり、ＦＩＲ成分に相当する入力変数の値が急激に変化した場合に、予測出力が計測出力に追従できなくなり予測に位相遅れが生じる事がある。しかし、本実施例の構成をとることにより、ＡＲ成分の影響を調整することが可能になり、位相遅れを改善して計測出力の値が変化する前に出力の変化を予測することが可能になる。
例えば、通常のＡＲＸモデルなどによってパラメータを同定すると自己回帰成分（ＡＲ成分）の重みが大きくなり、ＦＩＲ成分に相当する入力変数の値が急激に変化した場合に、予測出力が計測出力に追従できなくなり予測に位相遅れが生じる事がある。しかし、本実施例の構成をとることにより、ＡＲ成分の影響を調整することが可能になり、位相遅れを改善して計測出力の値が変化する前に出力の変化を予測することが可能になる。
【0174】
（第２実施例）
図６は、一実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの第２実施例について説明するための図である。
図６は、合成した予測に傾向（バイアス）を持たせたい場合に用いられる予測モデルの一例として、意図的に過大の予測および過小の予測を行う場合の予測モデルを示している。
図６は、一実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの第２実施例について説明するための図である。
図６は、合成した予測に傾向（バイアス）を持たせたい場合に用いられる予測モデルの一例として、意図的に過大の予測および過小の予測を行う場合の予測モデルを示している。
【0175】
本実施例のペアワイズ予測モデルＭＢは、回帰ブロックＢ１－ＢＮと、予測出力の合成ブロックＢ５と、を備え、合成ブロックＢ５が、過大推定予測値（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))）と、を出力する点において上述の第１実施例のペアワイズ予測モデルＭＡと異なっている。
【0176】
本実施例では、ペアワイズ予測モデルの位置母数と尺度母数とを用いて、通常の予測に加えて過大な予測および過小な予測を行うことにより、安全側の予測（リスクを回避する側の予測）などのニーズ（要求）に応じて、予測結果を微調整することができる。
【0177】
図７および図８は、第２実施例のペアワイズ予測モデルの効果を説明するための図である。
図７は、回帰手法の一種であるL1正則化を導入したLasso回帰による雨水流入予測値と測定値との一例を示している。
図８は、第２実施例の予測モデルによる雨水流入予測値と測定値との一例を示している。
図７は、回帰手法の一種であるL1正則化を導入したLasso回帰による雨水流入予測値と測定値との一例を示している。
図８は、第２実施例の予測モデルによる雨水流入予測値と測定値との一例を示している。
【0178】
図７に示す雨水流入予測値と測定値とを比較すると、予測誤差は小さいが、若干推定値が測定値よりも小さい結果となっており、過小推定になっている事がわかる。
一方、図８に示す本実施例の手法で得られた通常の推定（平均的な推定）、過大推定、および、過少推定を行った結果と測定値とを比較すると、通常の推定値（予測流入量）では予想誤差は小さいく、ピーク付近での予測精度は良いものの流入量が増加するタイミングでの予測値が測定値を下回っている。
一方、図８に示す本実施例の手法で得られた通常の推定（平均的な推定）、過大推定、および、過少推定を行った結果と測定値とを比較すると、通常の推定値（予測流入量）では予想誤差は小さいく、ピーク付近での予測精度は良いものの流入量が増加するタイミングでの予測値が測定値を下回っている。
【0179】
このような場合、この予測値を制御に用いる場合には、もう少し流入量急増時の予測を過大に評価しておきたい場合がある。このとき、過大推定予測を行うと、全体的にＭＳＥ（平均２乗誤差）での予測精度は劣化するが、流入量急増時に急増する可能性を示唆する予測が可能になることがわかる。また、通常の予測値よりも過小に評価をしておきたい場合には、過小推定予測を行うと、予測値が小さくなるように調整することができる。
【0180】
このような過大推定および過小推定により予測値を調整可能とすると、長期の運用の結果、通常の予測値が実測値に対して過小推定、あるいは、過大推定になっている様な傾向（バイアス）が見られる場合、通常の予測値を過大推定値、あるいは過小推定値に置き換えることで、バイアスを補正してより精度の高い予測を行うことができる。
【0181】
（第３実施例）
図９は、一実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの第３実施例について説明するための図である。
本実施例では、ペアワイズ予測モデルとして、（５）式のようなＡＲモデル／ＦＩＲモデルをさらに分解し、時間遅れを考慮した変数を各々一つの説明変数と見なした単回帰モデルに分解して予測モデルを構築している。
図９は、一実施形態のモジュラー型の時系列予測装置におけるペアワイズ予測モデルの第３実施例について説明するための図である。
本実施例では、ペアワイズ予測モデルとして、（５）式のようなＡＲモデル／ＦＩＲモデルをさらに分解し、時間遅れを考慮した変数を各々一つの説明変数と見なした単回帰モデルに分解して予測モデルを構築している。
【0182】
本実施例のペアワイズ予測モデルＭＣは、複数の時間遅れブロックと、複数の単回帰ブロックと、相関重み付き平均ブロックと、を含む予測出力ブロックＣ１－ＣＮを入力変数１－Ｎ毎に含み、当該複数の予測出力ブロックＣ１－ＣＮから出力された予測値を合成した最終合成予測値を出力する予測出力値の合成ブロックＣ（Ｎ＋１）を更に備える。
【0183】
予測出力ブロックＣ１は、複数（ｍ）の時間遅れブロック１１１－１１ｍと、複数（ｍ）の単回帰ブロック１２１－１２ｍと、相関重み付き平均ブロック１３と、を含む。
複数の時間遅れブロック１１１－１１ｍは、入力変数１（例えば時刻ｔにおける値）に対してそれぞれに割り当てられた遅れ時間Ｌ１１－Ｌ１ｍ１だけ遅れた変数（例えば時刻ｔ－Ｌ１１、…、ｔ－Ｌ１ｍ１における値）を出力する。
複数の時間遅れブロック１１１－１１ｍは、入力変数１（例えば時刻ｔにおける値）に対してそれぞれに割り当てられた遅れ時間Ｌ１１－Ｌ１ｍ１だけ遅れた変数（例えば時刻ｔ－Ｌ１１、…、ｔ－Ｌ１ｍ１における値）を出力する。
【0184】
複数の単回帰ブロック１２１－１２ｍは、単回帰係数を用いて、時間遅れブロック１１１－１１ｍから出力された変数に基づく予測値を演算する出力する。
相関重み付き平均ブロックは、複数の単回帰ブロック１２１－１２ｍから出力された複数の予測値の相関重み付き平均値を演算して出力する。
相関重み付き平均ブロックは、複数の単回帰ブロック１２１－１２ｍから出力された複数の予測値の相関重み付き平均値を演算して出力する。
【0185】
予測出力ブロックＣ２－ＣＮは、上記予測出力ブロックＣ１と同様の構成であるため、個々の説明は省略する。
予測出力値の合成ブロックＣ（Ｎ＋１）は、予測出力ブロックＣ１－ＣＮから出力された予測値を合成（例えば、相関重み付き平均）して、最終推定予測値（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))）を出力する。
【0186】
本実施例のペアワイズ予測モデルは、遅れ、ゲイン、バイアスの３項でモデルを説明することが可能であり、パラメータの微調整が容易である。また、単回帰係数と相関係数とには陽な関係があるため、各入力変数と各出力変数との相関係数を求めるだけで全体の予測モデルの構築ができ、予測モデルの解釈が容易である。
更に、本実施例のペアワイズ予測モデルでは、第１実施例のペアワイズ予測モデルと同様の効果を得ることができる。
更に、本実施例のペアワイズ予測モデルでは、第１実施例のペアワイズ予測モデルと同様の効果を得ることができる。
【0187】
図１０および図１１は、第３実施例のペアワイズ予測モデルの効果の一例を説明するための図である。
図１０および図１１では、第３実施例のペアワイズ予測モデルに対し１０個の入力変数を入力して雨水流入予測を行った例である。ここでは、各入力変数に対する遅れ時間を一つだけ考えた最も単純な場合の予測結果の例を示している。
図１０および図１１では、第３実施例のペアワイズ予測モデルに対し１０個の入力変数を入力して雨水流入予測を行った例である。ここでは、各入力変数に対する遅れ時間を一つだけ考えた最も単純な場合の予測結果の例を示している。
【0188】
入力変数の一つである降雨強度（地上雨量）の遅れ時間を推定すると３３分という推定結果が得られた。図１１は、上記推定結果を用いた雨水流入量の予測結果の一例を示している。図１１の予測結果によれば、全体的には、ある程度の精度で予測ができているが、流入量が急増する部分での予測の遅れ（位相遅れ）が大きい事がわかる。
【0189】
この時、位相遅れが大きい理由として、遅れ時間３３分というものが大きすぎる可能性が推測できる。この推測に基づき、図１０に、遅れ時間を１５分と調整した場合の雨水流入量の予測結果の一例を示す。図１０に示す予測結果によれば、遅れ時間を調整することにより、全体的な予測精度には大きな影響を与えることなく、流入量急増時の位相遅れを大幅に改善できていることがわかる。
【0190】
このように、本実施例のペアワイズ予測モデルによれば、容易に解釈可能なパラメータのみで予測モデルを構築可能であり、説明性が向上するだけでなく予測モデルの微調整を容易に行うことができる。
【0191】
また、単回帰係数と相関係数とは、入出力データを各々平均と標準偏差とで正規化すれば、一致することが知られており、これを逆に非正規化すれば、回帰係数とバイアスとに換算することができる。したがって、各入出力データの平均、標準偏差、および、入出力の相関のみを計算すれば、遅れ時間を指定するだけで、全ての予測モデルのパラメータを同定することができる。
【0192】
さらに、合成出力の計算を相関係数の重み付き平均で計算する方法とすることにより、予測モデルを構築するために必要な計算は、遅れ時間の推定（指定）と、入出力変数の平均および標準偏差と、入出力変数の相関と、のみであり、予測モデルに含まれるパラメータの全てを、この４種類の推定値（遅れ時間、平均、標準偏差、相関係数）だけを用いて説明することが可能になる。
【0193】
（効果）
上記のように、本実施形態によれば、入力変数と出力変数との相関関係が変化したり劣化したり、入力変数の時系列データに多量のアウトライアが含まれたりするなどの理由で予測精度が十分でない場合に、全体の予測モデルを再同定することなく、合成した予測出力値に影響を与えるペアワイズの予測モデルの部分的な調整や、悪影響を与える部分の削除を容易に行う事ができると同時に、これにより、予測精度劣化の理由を合理的に説明することが可能になる。
上記のように、本実施形態によれば、入力変数と出力変数との相関関係が変化したり劣化したり、入力変数の時系列データに多量のアウトライアが含まれたりするなどの理由で予測精度が十分でない場合に、全体の予測モデルを再同定することなく、合成した予測出力値に影響を与えるペアワイズの予測モデルの部分的な調整や、悪影響を与える部分の削除を容易に行う事ができると同時に、これにより、予測精度劣化の理由を合理的に説明することが可能になる。
【0194】
また、入力変数間の多重共線性が極めて強くほぼ同じ説明能力を持つ入力変数が含まれる場合にも、パラメータの可同定(一意決定可能性)を維持することができ、物理法則に矛盾しない合理的な説明が可能な予測モデルを構築できると同時に、ほぼ同じ説明能力を持つ入力変数を統合することで、予測精度を維持しながら、説明性を向上させることができる。
すなわち、本実施形態によれば、予測結果を合理的に説明するとともに、容易に予測結果を調整可能とするモジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラムを提供することができる。
すなわち、本実施形態によれば、予測結果を合理的に説明するとともに、容易に予測結果を調整可能とするモジュラー型時系列データ予測装置、モジュラー型時系列データ予測方法、および、プログラムを提供することができる。
【0195】
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。
【符号の説明】
【0196】
１…都市雨水排水プロセス、１１…流入量計、１２…幹線流量計、１３１－１３Ｋ…幹線水位計、１４１～１４Ｍ…地上雨量計、１５、１５１１～１５ＱＰ…レーダ雨量計、１７…流入渠、１１０…流入ゲート、１８…雨水ポンプ井、１９…雨水ポンプ、１６１、１６２…雨水ポンプ井水位計、２…データ収集保存部、２１…出力変数データ選択部、２２…入力変数データ選択部、３…データ抽出部、３１…オフライン予測モデル同定用データ抽出部、３２…オンライン予測用データ抽出部、３３…評価用データ抽出部、４…ペアワイズ予測モデル同定部、５…予測モデル合成法定義部、６…ペアワイズ出力変数予測部、７…合成出力変数予測部、８…予測誤差評価部、９…ペアワイズ予測モデル修正部、１０…出力予測結果観測部、１１１－１１ｍ…時間遅れブロック、１２１－１２ｍ…単回帰ブロック、１３…相関重み付き平均ブロック、Ａ１－ＡＮ…回帰ブロック、Ａ５…合成ブロック、Ｂ１－ＢＮ…回帰ブロック、Ｂ５…合成ブロック、Ｃ１－ＣＮ…予測出力ブロック、Ｃ（Ｎ＋１）…合成ブロック。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
