(57)【要約】      （修正有）【課題】設定領域内の特性値を簡易な手法で精度よく推定可能とする。【解決手段】特性値推定装置は、第１計算モデルを入力するモデル入力部と、第２計算モデルに対してセンサデータに基づいて学習を行うモデル学習部と、第１計算モデルと第２計算モデルとのいずれか一方を選択するモデル切替部と、特性値を計算するためのシミュレーションを実行して得られた値と、センサのセンサデータと、を比較して、計算モデルの誤差を計算する予測値計算部と、誤差が所定値以下になるまで未確定のパラメータの確率分布を修正する、確率分布修正部と、センサが配置されていない場所に仮想的に配置される仮想センサのセンサデータを推定する仮想センサ値推定部と、センサが配置されていない場所での特性値の詳細分布を推定する特性値分布推定部と、特性値分布推定部で推定された特性値の詳細分布の信頼度を計算する信頼度計算部と、を備える。【選択図】図１
2018049390.tif 000002
【特許請求の範囲】
【請求項１】
特性値を監視する一つ以上のセンサで検出されたセンサデータを入力する、センサデータ入力部と、
前記特性値を計算するための第１計算モデルを入力する、モデル入力部と、
前記特性値を計算するための第２計算モデルに対して前記センサデータに基づいて学習を行う、モデル学習部と、
前記第１計算モデルと前記第２計算モデルとのいずれか一方を選択する、モデル切替部と、
前記モデル切替部が選択した計算モデルと、前記計算モデルに含まれる複数のパラメータのうち未確定のパラメータの確率分布と、前記センサのセンサデータと、に基づいて、前記特性値を計算するためのシミュレーションを実行して得られた値と、前記センサのセンサデータと、を比較して、前記計算モデルの誤差を計算する予測値計算部と、
前記誤差が所定値以下になるまで前記未確定のパラメータの確率分布を修正する、確率分布修正部と、
前記誤差に基づいて、前記センサが配置されていない場所に仮想的に配置される仮想センサのセンサデータを推定する仮想センサ値推定部と、
前記誤差と、前記仮想センサのセンサデータと、前記センサのセンサデータと、に基づいて、前記センサが配置されていない場所での前記特性値の詳細分布を推定する特性値分布推定部と、
前記特性値分布推定部で推定された前記特性値の詳細分布の信頼度を計算する信頼度計算部と、を備える、特性値推定装置。
【請求項２】
前記信頼度に基づいて、前記誤差をより小さくするための前記センサの配置プランを計算するセンサ配置計算部と、
前記計算された配置プランに関する情報を出力するセンサ配置出力部と、
前記誤差に基づいて、前記モデル切替部が選択した計算モデルを修正するモデル修正部と、
前記修正された計算モデルに関する情報を出力するモデル修正出力部と、を備える、請求項１に記載の特性値推定装置。
【請求項３】
前記特性値分布推定部は、離散的に配置された前記センサのセンサデータを補間することにより、前記センサが配置されていない場所に配置された前記仮想センサのセンサデータを推定し、
前記信頼度計算部は、前記特性値分布推定部が推定した前記仮想センサのセンサデータの前記信頼度を計算する、請求項１または２に記載の特性値推定装置。
【請求項４】
前記未確定のパラメータの確率分布を事前分布として入力する確率分布入力部と、
前記確率分布修正部が修正した前記未確定のパラメータの確率分布を事後分布として記録する確率分布記録部と、を備え、
前記予測値計算部は、前記確率分布修正部が前記未確定のパラメータの確率分布を修正するまでは、前記確率分布入力部にて入力された前記未確定のパラメータの確率分布を事前分布として用いて前記シミュレーションを実行し、前記確率分布修正部が前記未確定のパラメータの確率分布を修正した後に前記センサデータ入力部から新たなセンサデータが入力された場合には、前記確率分布記録部に記録された確率分布を事前分布として用いて前記シミュレーションを実行する、請求項１乃至３のいずれか一項に記載の特性値推定装置。
【請求項５】
前記特性値分布推定部は、ガウス過程を利用した回帰を行って、前記センサが配置されていない場所に配置された前記仮想センサのセンサデータを推定する、請求項１乃至３のいずれか一項に記載の特性値推定装置。
【請求項６】
前記モデル切替部は、空間を分割する粒度を調整して、所定の粒度の計算モデルを選択する、請求項１乃至５のいずれか一項に記載の特性値推定装置。
【請求項７】
前記特性値は温度であり、
前記未確定のパラメータは、人体発熱率および設定温度を含む、請求項１乃至６のいずれか一項に記載の特性値推定装置。
【請求項８】
前記センサは、温度センサであり、
前記特性値分布推定部は、前記温度センサの設置間隔よりも細かい間隔で温度分布を推定する、請求項７に記載の特性値推定装置。
【請求項９】
特性値を監視する一つ以上のセンサで検出されたセンサデータを入力するステップと、
前記特性値を計算するための第１計算モデルを入力するステップと、
前記特性値を計算するための第２計算モデルに対して前記センサデータに基づいて学習を行うステップと、
前記第１計算モデルと前記第２計算モデルとのいずれか一方を選択するステップと、
前記選択された計算モデルと、前記計算モデルに含まれる複数のパラメータのうち未確定のパラメータの確率分布と、前記センサのセンサデータと、に基づいて、前記特性値を計算するためのシミュレーションを実行して得られた値と、前記センサのセンサデータと、を比較して、前記計算モデルの誤差を計算するステップと、
前記誤差が所定値以下になるまで前記未確定のパラメータの確率分布を修正するステップと、
前記誤差に基づいて、前記センサが配置されていない場所に仮想的に配置される仮想センサのセンサデータを推定するステップと、
前記誤差と、前記仮想センサのセンサデータと、前記センサのセンサデータと、に基づいて、前記センサが配置されていない場所での前記特性値の詳細分布を推定するステップと、
前記推定された前記特性値の詳細分布の信頼度を計算するステップと、を備える、特性値推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、特性値推定装置および特性値推定方法に関する。
【背景技術】
【０００２】
従来の、システムモニタリングの技術は、未確定パラメータをシミュレーション結果と観測データとの比較から同定してきたが、同定したパラメータの信頼度や、モデル化そのものの誤差の情報は提供できなかった。また、統計的に少ないサンプルデータから応答局面を求める従来の技術は、対象の物理モデルを無視しているので、物理モデルとのズレも含めた推定結果の信頼度情報を提供できなかった。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０１４−２０６８７０号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
本実施形態は、設定領域内の特性値を簡易な手法で精度よく推定可能な特性値推定装置および特性値推定方法を提供するものである。
【課題を解決するための手段】
【０００５】
本実施形態によれば、特性値を監視する一つ以上のセンサで検出されたセンサデータを入力する、センサデータ入力部と、
前記特性値を計算するための第１計算モデルを入力する、モデル入力部と、
前記特性値を計算するための第２計算モデルに対して前記センサデータに基づいて学習を行う、モデル学習部と、
前記第１計算モデルと前記第２計算モデルとのいずれか一方を選択する、モデル切替部と、
前記モデル切替部が選択した計算モデルと、前記計算モデルに含まれる複数のパラメータのうち未確定のパラメータの確率分布と、前記センサのセンサデータと、に基づいて、前記特性値を計算するためのシミュレーションを実行して得られた値と、前記センサのセンサデータと、を比較して、前記計算モデルの誤差を計算する予測値計算部と、
前記誤差が所定値以下になるまで前記未確定のパラメータの確率分布を修正する、確率分布修正部と、
前記誤差に基づいて、前記センサが配置されていない場所に仮想的に配置される仮想センサのセンサデータを推定する仮想センサ値推定部と、
前記誤差と、前記仮想センサのセンサデータと、前記センサのセンサデータと、に基づいて、前記センサが配置されていない場所での前記特性値の詳細分布を推定する特性値分布推定部と、
前記特性値分布推定部で推定された前記特性値の詳細分布の信頼度を計算する信頼度計算部と、を備える、特性値推定装置が提供される。
【図面の簡単な説明】
【０００６】
【図１】第１の実施形態による特性値推定装置の概略構成を示すブロック図。
【図２】図１の特性値推定装置が使用する計算モデルの一例を説明する図。
【図３】図１の特性値推定装置の全体的な処理動作を示すフローチャート。
【図４】図１の予測値計算部の処理動作を示すフローチャート。
【図５】一日分の室温変化のシミュレーション結果を示すグラフ。
【図６】ＭＣＭＣを行った結果の一例を示す図。
【図７】事後分布の平均値を用いて計算モデルのシミュレーションを行った結果と観測データとを示す図。
【図８】ガウス過程により得られた事後分布を相関プロット形式で表示した図。
【図９】得られたガウス過程δから関数をサンプリングした結果とモデルη(ｘ，θ)とを合わせてプロットした図。
【図１０】部屋に設置した温度センサの位置を２次元に射影した地図上に表示した図。
【図１１】設置された１４個の温度センサから得られたある１日の温度変化をプロットしたグラフ。
【図１２】ガウス過程回帰の一例を示す図。
【図１３】推定された部屋単位の温度の時間推移をプロットしたグラフ。
【図１４】室内の３次元グリッド上に配置した仮想センサでのセンサ値を計算してプロットした図。
【図１５】第２の実施形態による特性値推定装置１の概略構成を示すブロック図。
【発明を実施するための形態】
【０００７】
以下、図面を参照しながら、本発明の実施形態を説明する。
【０００８】
（第１の実施形態）
図１は第１の実施形態による特性値推定装置１の概略構成を示すブロック図である。図１の特性値推定装置１は、モデル入力部２と、センサデータ入力部３と、確率分布入力部４と、モデル学習部５と、仮想センサ情報入力部６と、モデル切替部７と、確率分布記録部８と、予測値計算部９と、確率分布修正部１０と、仮想センサ値推定部１１と、再現詳細度回復部（特性値分布推定部）１２と、モデル誤差記録部１３と、信頼度計算部１４と、信頼度出力部１５と、仮想センサ値出力部１６とを備えている。
【０００９】
モデル入力部２は、特性値を計算するための第１計算モデルを入力する。より具体的には、モデル入力部２は、対象となるプラントやプロセスの振る舞いを模擬する計算モデルを入力する。計算モデルは、図１の特性値推定装置１の外部または内部で生成されて、モデル入力部２に入力される。ここで、モデルとは、対象となるプラントや機器の振る舞いを模擬するシステムの記述であり、専用に設けられたシミュレータを用いて、モデルを解釈して何らかの出力が得られるものである。
【００１０】
センサデータ入力部３は、特性値を監視する一つ以上のセンサで検出されたセンサデータを入力する。より具体的な一例では、センサデータ入力部３には、プロセスに関連するセンサデータが入力される。センサの機能や種類、数は特に問わない。センタの一例は温度センサである。
【００１１】
確率分布入力部４は、計算モデルに用いられるパラメータのうち、該当するセンサデータが得られない未確定パラメータの事前確率分布を入力する。確率分布入力部４に入力された事前確率分布は、確率分布記録部８に記録される。未確定パラメータの具体例については後述する。
【００１２】
モデル学習部５は、特性値を計算するための第２計算モデルに対して、センサデータに基づいて学習を行う。より詳細には、モデル学習部５は、対象としているプロセスの振舞をセンサデータから学習する。データからシステムの応答局面を求める手法を利用できる。例えば、ガウス過程を利用すると、データから、このデータを生成する関数の確率分布を推定できる。このように、モデル学習部５は、センサデータから応用局面として学習した計算モデルを作成する。
【００１３】
モデル切替部７は、第１計算モデルと第２計算モデルのいずれか一方を選択する。また、モデル切替部７は、空間を分割する粒度を調整可能な場合は、適切な粒度に空間を分割した計算モデルを選択する。
【００１４】
より詳細には、モデル切替部７は、モデル入力部２により入力された計算モデル（第１計算モデル）と、モデル学習部５がセンサーデータから応答局面として学習したモデル（第２計算モデル）とを切り替えて予測値計算部９に与える。例えばモデルを用いた予測が観測データとの乖離が激しく、モデル自体の妥当性が疑われる場合に、学習したモデルを用いた予測に切り替える。システム全体としてモデルが明確な部分と不明な部分が混在する場合に、計算モデルと学習モデルを切り替えて予測値計算を進めることができる。
【００１５】
また、モデル切替部７は、計算モデルの粒度が調整可能な場合には、これを切り替える機能を持つ。
【００１６】
本実施形態では、後述するように、部屋を１つの点として全体の室温を計算するが、モデル化では、空間を細かく分割して分割された空間単位で室温を推定することもできる。空間を分割すれば計算量が増大するとともに、観測データを分割空間毎に準備しなくはならなくなるので、コストがかさむ。一方分割を少なくすれば、空間方向の細かい変動に予測が追従できない。モデル切替部７は、トレードオフ関係を分析し、適切な粒度のモデルを選択する。
【００１７】
仮想センサ情報入力部６は、仮想センサの設置場所やセンシングするプロセス状態量を指定する情報を入力して記録する。
【００１８】
予測値計算部９は、モデル切替部７が選択した計算モデルと、計算モデルに含まれる複数のパラメータのうち未確定のパラメータの確率分布と、センサのセンサデータと、に基づいて、特性値を計算するためのシミュレーションを実行して得られた値とセンサのセンサデータとを比較して、計算モデルの誤差を計算する。予測値計算部９は、上述した確率分布の修正とシミュレーションを繰り返し行って、誤差ができるだけ小さくなるようにする。誤差が大きければ、未確定パラメータの確率分布を修正してシミュレーションを繰り返す。
【００１９】
より具体的には、予測値計算部９は、センサデータと未確定パラメータの確率分布とに基づいて、モデル切替部７が選択した計算モデルを用いてシミュレーションを実行し、そのシミュレーション結果を出力する。
【００２０】
確率分布修正部１０は、予測値計算部９で計算された誤差が所定値以下になるまで、未確定のパラメータの確率分布を修正する。修正された確率分布は確率分布修正部１０に記録される。
【００２１】
仮想センサ値推定部１１は、予測値計算部９で計算された誤差に基づいて、センサが配置されていない場所に仮想的に配置される仮想センサのセンサデータを推定する。仮想センサ値推定部１１は、予測値計算部９のシミュレーションデータから仮想センサのセンサデータを取得する。
【００２２】
再現詳細度回復部１２は、予測値計算部９で計算された誤差と、仮想センサのセンサデータと、センサのセンサデータとに基づいて、センサが配置されていない場所での特性値の詳細分布を推定する。よって、再現詳細度回復部１２は、特性値分布推定部に対応するものである。本実施形態では、再現詳細度回復部１２が推定する特性値の詳細分布を、再現詳細度の回復と呼ぶ。より具体的には、再現詳細度回復部１２は、離散的に配置されたセンサのセンサデータを補間することにより、センサが配置されていない場所に配置された仮想センサのセンサデータを詳細に推定する。
【００２３】
再現詳細度回復部１２は、予測値計算部９による予測値の再現詳細度が観測点（センサ位置）の再現詳細度に対して低い場合に、観測点の再現詳細度まで回復させる。本実施形態における再現詳細度とは、オフィス内の温度分布の推定である。温度計はオフィス内の限られた場所に偏りをもって設置されていることが多いが、本実施形態では、このような温度計の計測値に基づいて、オフィス内の温度分布を推定する。
【００２４】
より具体的には、再現詳細度回復部１２は、部屋単位での予測値計算を行った結果を平均として、実際の観測点から得られた温度データを補間して分布を求めることにより、空間の再現詳細度を回復する。再現詳細度の回復は、離散データに対する補間処理により実現できる。本実施形態では、ガウス過程に代表される確率過程を用いた再現詳細度回復を行う。
【００２５】
ここで、ガウス過程ＧＰ(ｍ，Σ)は、中央関数ｍを、観測点座標の共変行列Σにより定義される。予測値計算部９で得られたモデル出力をモデル化誤差で補正した予測値を中央関数であるとして、これに、離散的な観測点の空間情報を共変行列Σとして与えれば、観測点の間を補完する時空間関数を推定することができる。本実施形態では、推定された部屋全体の温度予測を中央関数として、観測された空間的に離散なセンサー配置のデータから空間分布を補った推定を行う。
【００２６】
モデル誤差記録部１３は、予測値計算部９で計算された誤差を計算モデル自体の誤差として記録する。
【００２７】
信頼度計算部１４は、再現詳細度回復部１２で推定された特性値の分布の信頼度を計算する。より具体的には、信頼度計算部１４は、仮想センサのセンサデータの信頼度を計算する。本実施形態では、ガウス過程回帰により、温度の３次元分布を推定するが、このとき推定値の信頼度を３次元分布として推定することができる。ガウス過程回帰では、分散ｃ−ｋTＣN-1ｋが信頼度に相当する。
【００２８】
仮想センサ値出力部１６は、再現詳細度回復部１２が推定した特性値の分布に基づいて、仮想センサのセンサデータを出力する。
【００２９】
図２は図１の特性値推定装置１が使用する計算モデルの一例を説明する図である。図２は、オフィス内を対象とした熱交換過程をモデル化した例を示している。ある部屋における温度変化を計算するモデルである。部屋は、壁で区切られた閉空間である。また、部屋は、壁を介して外部との熱交換を行うとともに、空調装置を通じた熱の流入・流出を行う。さらに、部屋には、室内のＯＡ（Office Automation）機器や部屋内の人間の発熱する熱量が流入する。部屋は熱容量を持っているので、流入した熱量は結果として室温に反映される。例えば、空気調整器（以下、空調機）から流入した熱量により室温が上がれば、外部との温度差が広がるので、より多くの熱量が壁を介して外に流出し、室温が下がる。このように部屋を中心にした熱の流入・流出は、熱量保存則に従って互いに連携して動作する。
【００３０】
図２は部屋の熱の収支を矢印で示している。ＨＶＡＣroomは、空調機が生成する熱量である。Ｑpersonは、部屋内の人間一人当たりの発熱量である。ＱOAは、ＯＡ機器の発熱量である。Ｑoutは、部屋の外部への流出熱量である。
【００３１】
図２では、部屋内部と外部との間の壁を通じた熱流出Ｑoutと、空調機から流入および流出する熱量ＨＶＡＣout、室内のＯＡ機器ＱOAと、人間からの発熱Ｑpersonとが合算して、部屋への熱量の収支が計算される。また、部屋の熱容量Ｃpから室温Ｔaが計算される。さらに、室温Ｔaと外気温ｔｏとの差から熱流出Ｑoutが再計算されて熱収支が変化する。部屋内のすべての人間の発熱量は、以下の（１）式で表される。
【００３２】
Ｑperson＝ｐｎｕｍ(ｉ)×PowerDenPerson    …（１）
【００３３】
（１）式では、ある時刻ｉにおける部屋内の人間の数をｐｎｕｍ(ｉ)とし、人体発熱率をPowerDenPersonとしている。部屋内のすべてのＯＡ機器の発熱量は、以下の（２）式で表される。
【００３４】
ＱOA＝ｐｎｕｍ(ｉ)×PowerDenOA    …（２）
【００３５】
（２）式では、ある時刻ｉにおける部屋内のＯＡ機器の数をｐｎｕｍ(ｉ)とし、ＯＳ機器一台当たりの熱量をPowerDenOAとしている。部屋外への流出熱量Ｑoutは、以下の（３）式で表される。
【００３６】
Ｑout＝ｑ×(ｔｏ(ｉ)−Ｔa(ｉ))    …（３）
【００３７】
（３）式において、Ｔa(ｉ)は時刻ｉにおける室温、ｔｏ(ｉ)は時刻ｉにおける外気温である。またqは壁固有の熱伝導率であり固定値とする。空調機の生成熱量ＨＶＡＣroomは、以下の（４）式で表される。
【００３８】
ＨＶＡＣroom＝−ControllerGain×(Ｔa−Ｔsc)    …（４）
【００３９】
ControllerGainは係数、Ｔaは室温、Ｔscは設定温度である。空調機の部屋への熱量ＨＶＡＣoutは、以下の（５）式で表される。
【００４０】
ＨＶＡＣout＝HVACroom÷ＣＯＰ÷Ｃｏｍｐefficiency    …（５）
【００４１】
（５）式において、HVACroomは空調機が生成する熱量、ＣＯＰは熱効率である。時刻ｉにおける空調機電力ｑ(ｉ)は、以下の（６）式で表される。
【００４２】
ｑ(ｉ)＝ＨＶＡＣout÷１０００    …（６）
【００４３】
部屋の全体の熱量Ｑは、以下の（７）式で表される。ここでdtをシミュレーション時間区切り、dTを時間区切りdt前後の室温の変化とする。
【００４４】
Ｑ＝Ｑperson＋ＱOA＋ＨＶＡＣroom＋Ｑout    …（７）
ｄＴ＝Ｑ÷Ｃp×ｄｔ    …（８）
【００４５】
時刻(ｉ＋１)における室温Ｔa(ｉ＋１)は、以下の（９）式で表される。
【００４６】
Ｔa(ｉ＋１)＝Ｔa(ｉ)＋ｄＴ    …（９）
【００４７】
上述した（１）の人体発熱率PowerDenPersonと（４）式の設定温度Ｔscは、事前に値が分かっていない未確定パラメータである。未確定パラメータに何らかの値を設定してシミュレーションを実行すると、シミュレーションにて計算された室温と室温の観測値（実測値）とを比較できる。
【００４８】
未確定パラメータを種々に変更することで、シミュレーションによる室温と観測値との誤差が変化する。そこで、未確定パラメータの値を種々に修正し、室温の計算結果と室温の観測値とが合致するように、人体発熱率pdp＝PowerDenPersonと設定温度Ｔscの値を調整することができる。
【００４９】
図３は図１の特性値推定装置１の全体的な処理動作を示すフローチャートである。図１の特性値推定装置１では、物理的または化学的な保存則などの最低限守らなければならない法則を用いてモデル化した計算モデルを用意し、この計算モデルのパラメータを、センサによるセンサデータを用いてキャリブレーションを行う（ステップＳ１）。ステップＳ１のキャリブレーションは、図１の予測値計算部９が行う。予測値計算部９と確率分布修正部１０は、計算モデルのパラメータの中で、センサデータでキャリブレーションできない未確定パラメータの確率分布の修正も行う。
【００５０】
次に、仮想センサ値推定部１１は、予測値計算部９が計算した誤差に基づいて、仮想センサのセンサデータを推定する（ステップＳ２）。次に、再現詳細度回復部１２は、予測値計算部９が計算した誤差と、仮想センサのセンサデータと、センサのセンサデータとに基づいて、再現詳細度を回復した仮想センサのセンサデータを推定する（ステップＳ３）。次に、信頼度計算部１４は、再現詳細度回復部１２にて最後に推定された仮想センサのセンサデータの信頼度を計算する（ステップＳ４）。
【００５１】
図４は図１の予測値計算部９の処理動作を示すフローチャートである。まず、センサデータ入力部３に入力されたセンサのセンサデータを取得する（ステップＳ１１）。次に、確率分布入力部４に入力された未確定パラメータの確率分布を取得する（ステップＳ１２）。次に、センサデータと、未確定パラメータの確率分布とを用いて、計算モデルに基づいてシミュレーションを実行する（ステップＳ１３）。次に、シミュレーション結果とセンサデータとを比較する（ステップＳ１４）。シミュレーション結果とセンサデータとの誤差が所定値以下であれば、図４の処理を終了し、誤差が所定値より大きければ、未確定パラメータの確率分布を修正し（ステップＳ１５）、ステップＳ１３以降の処理を繰り返す。
【００５２】
次に、本実施形態の具体例について説明する。観測値をｙ、計算モデルをη(ｘ，θ)、モデル化誤差をδ(ｘ)、観測誤差をε(ｘ)、モデルにおける未確定パラメータをθとすると、以下の（１０）式が得られる。
【００５３】
ｙ(ｘ)＝η(ｘ，θ)＋δ(ｘ)＋ε(ｘ)    …（１０）
【００５４】
モデル化誤差δとは、計算モデルηを用いて観測値ｙへのキャリブレーションを行った結果、最後に残る誤差部分に相当する。モデルη自体が実プロセスとずれている場合、あるいは予期せぬ変動要因がある場合にモデル化誤差δが生成される。
【００５５】
キャリブレーションとは、モデルηと観測データｙ＝(ｙ(ｘ1)，…，ｙ(ｘn))が与えられたときに、観測データｙを良く説明する、未確定パラメータθとモデル化誤差δ(x)を同時に推定する処理であるとする。
【００５６】
以下では、予測値計算部９が、観測データとシミュレーション結果とを比較し、観測データとシミュレーション結果とが合致するように未確定パラメータを修正する例を具体的に説明する。
【００５７】
未確定パラメータの値がある確率分布を有するものとして、この確率分布をシミュレーションしながら修正する処理を検討する。この修正処理をキャリブレーションと呼ぶこととする。キャリブレーションでは、未確定パラメータの事前の確率分布 ｐ(θ)に対して、観測データが推定値に合うように確率分布を修正しながらシミュレーションを繰り返し、最終的に事後分布 ｐ(θ|ｙ)を求めることとする。
【００５８】
本実施形態では、事後確率分布を求めるために、ＭＣＭＣ（マルコフ連鎖モンテカルロ法）を利用する。ＭＣＭＣでは、パラメータθの値を、事前確率分布に従ってランダムに選択し、シミレーション結果と観測データの比較から、次のパラメータθの選択を更新しつつ、シミュレーションを繰り返す。結果としてパラメータθの選択履歴をヒストグラム化すると、事後確率分布が得られる。
【００５９】
図５は、人体発熱率ｐｄｐ＝１００、設定温度Ｔsc＝２７を正解値として、一日分の室温変化のシミュレーション結果を示すグラフである。図５の横軸は時間、縦軸は温度である。図５には、温度変化にノイズとモデル誤差を人為的に導入して生成した、１時間毎の観測データｙをエラー幅付きで表示している。図５からわかるように、シミュレーション結果と観測データとは乖離している。
【００６０】
ここで、人体発熱率ｐｄｐと設定温度Ｔscの事前分布は、以下のガウス分布Ｎに従うものとする。
【００６１】
ｐ(pdp)＝Ｎ(μ＝80，σ2＝1/0.002)    …（１１）
ｐ(tsc)＝Ｎ(μ＝15，σ2＝1/0.01)    …（１２）
【００６２】
このガウス分布Ｎにおいて、μは平均、σ2は分散であり、τ＝１／σ2はprecisionと呼ばれる。１変数ｘのガウス分布は、以下の（１３）式で表される。
【００６３】
【数１】
2018049390.tif 000003
【００６４】
本実施形態で事後確率分布を求めるために用いたＭＣＭＣとは、多変量の確率分布からサンプルを抽出するためのアルゴリズムの一つである。具体的には、ｐ(Ａ｜Ｂ)を、条件ＢのときのＡの確率であるとすると、ベイズ則より、以下の（１４）式が成立することが知られている。
【００６５】
ｐ(Ａ｜Ｂ)〜ｐ(Ｂ｜Ａ)ｐ(Ａ)    …（１４）
【００６６】
ここで、時間ｘで観測された室温がｙであるとすると、ｙが観測されたときの人体発熱率ｐｄｐと設定温度Ｔscの事後分布は、以下の（１５）式のように記述できる。
【００６７】
ｐ(pdp，tsc｜ｙ)〜ｐ(ｙ｜ｘ，pdp，tsc)ｐ(pdp)ｐ(tsc)    …（１５）
【００６８】
ここで、人体発熱率ｐｄｐと設定温度Ｔscが与えられたときのシミュレーション結果と観測室温から、ｐ(ｙ｜ｘ，pdp，tsc)を求めることができる。
【００６９】
多変量ガウス分布Ｎ(ｙ｜ｍ，Ｋ)におけるｍは平均関数（ここではシミュレーション結果）、Ｋは観測誤差の共分散行列であり、Ｋii＝σi2δijに対して、確率分布すなわち尤度は、以下の（１６）式で表される。
【００７０】
【数２】
2018049390.tif 000004
【００７１】
尤度を最大にするｍを求めるためには、対数尤度として探索すればよい。対数尤度は、以下の（１７）式で表される。
【数３】
2018049390.tif 000005
【００７２】
ＭＣＭＣでは、パラメータｐｄｐ、ｔｓｃを設定したシミュレーション結果を中央関数ｍとし、対数尤度が高くなるようにパラメータｐｄｐ、ｔｓｃを変化させながら計算を繰り返し、結果としてパラメータｐｄｐ、ｔｓｃの試行の軌跡が事後分布ｐ(pdp，tsc｜ｙ)となるような計算を行う。
【００７３】
図６はＭＣＭＣを行った結果の一例を示す図である。人体発熱率のパラメータｐｄｐの平均値＝７７．４１３、設定温度のパラメータｔｓｃの平均値＝２７．３９の事後分布を得ることができた。
【００７４】
図７は事後分布の平均値を用いて計算モデルのシミュレーションを行った結果と観測データとを示す図である。図７に示すように、シミュレーション結果と観測データとがほぼ一致していることがわかる。
【００７５】
ＭＣＭＣを用いてパラメータｐｄｐ、ｔｓｃの事後分布を求める代わりに、ガウス過程（Gausian Process）を用いることにより連続値領域でモデル化することもできる。ガウス過程は、確率分布が多次元正規分布に従う確率過程であり、観測データの元となった関数の確率分布を推定する。空間分布を推定したり、時間変化を加味した時空間分布を推定することに用いられる。地理空間情報の分野においては、地図上の離散点で観測されたデータから、空間上の分布を推定する技術（Kriging）として利用されてきた。
【００７６】
ガウス過程は、関数の分布を表しており、中央値に相当する中央関数ｍ(ｘ)と、カーネル関数ｋ／２により特徴付けられる共変行列Ｋij=ｋ(ｘi，ｘj )によって特徴づけられる、多変量ガウス分布である。
【００７７】
ベクトルｙが中央関数ｍ、共変行列Ｋで定義される多変量ガウス分布に従うことは、以下の（１８）式で表すことができる。
【００７８】
ｙ〜Ｎ(ｍ，Ｋ)    …（１８）
Ｄ次元ベクトルｙに対する多変量ガウス分布の確率分布は、以下の（１９）式で表される。
【数４】
2018049390.tif 000006
【００７９】
ここで、観測値をｙ、計算モデルをη(ｘ，θ)、モデル化誤差をδ(ｘ)、観測誤差をε(ｘ)、モデルにおける未確定パラメータをθとすると、上述した（１０）式が成立する。
【００８０】
モデル誤差δ(ｘ)は、以下の（２０）式に示すガウス過程に従うものとする。
δ(ｘ)〜Ｎδ(０，Ｋ)    …（２０）
【００８１】
尤度ｐ(ｙ｜ｘ，θ，σ2)が最大となる未確定パラメータθとガウス過程δ(ｘ)（すなわち、カーネル関数のハイパーパラメータ）を求めるには、尤度の対数を最大にするアルゴリズムがよく用いられる。
【００８２】
ここでは、ｙ(ｘ)−η(ｘ，θ)がノイズε(ｘ)付きガウス過程に従うとすると、対数尤度は以下の（２１）式のように記述できる。
【数５】
2018049390.tif 000007
【００８３】
（２１）式のベクトルｒは、ｙ(ｘ)−η(ｘ，θ)に対応する（２２）式に示す差分ベクトルである。
【数６】
2018049390.tif 000008
【００８４】
共変行列Ｋは、もともとのモデル化誤差δ(ｘ)のカーネル関数をｋ(ｔi，ｔj)であるとすると、以下の（２３）式で表される。
Ｋij＝σi2δij＋ｋ(ｔi，ｔj)    …（２３）
【００８５】
ここで、σi2は観測誤差εの分散、δijはクロネッカーのδであるとする。カーネル関数ｋ(ｔi，ｔj)がMatern-3/2カーネルであるとすると、ｒ＝√(ｔi−ｔj)2であり、以下の（２４）式のようになる。
【数７】
2018049390.tif 000009
【００８６】
４つのパラメータａ、τ、ｐｄｐ、ｔｓｃの事前分布から、計算モデルとモデル誤差を考慮し、対数尤度を最大にするような探索を行い、これらのパラメータの事後分布を求めることができる。
【００８７】
図８はガウス過程により得られた事後分布を相関プロット形式で表示した図である。また、図９は得られたガウス過程δから関数をサンプリングした結果とモデルη(ｘ，θ)とを合わせてプロットした図である。
【００８８】
以上のようにして、予測値計算部９は、未確定パラメータθの分布を推定するとともに、モデル化誤差δを同時に推定することができる。本実施形態では、ガウス過程に用いるカーネル関数を特定の関数形式すなわち、Matern-3/2カーネルを採用したが、これ以外のカーネル関数を選択することもできる。
【００８９】
得られた事後分布は、確率分布記録部８に記録される。記録された事後分布は、新たなセンサデータが得られた場合の事前分布として活用される。また、得られたモデル化誤差δは、モデル誤差記録部１３に記録される。
【００９０】
仮想センサ値出力部１６は、予測値計算部９が予測値を計算するために計算モデルηを評価した中間結果を利用して、実際にはセンサによる観測が行われていない物理量、あるいは観測ができないプロセス特性を推定する。仮想センサ値出力部１６は、仮想センサ情報入力部６により入力された仮想センサの定義式を評価して、目的とする仮想センサ値を計算する。
【００９１】
本実施形態では、例えば部屋に溜まる熱量Ｑを直接計測することはできないが、計算過程で得られるパラメータＱの値を仮想センサ値として利用することができる。また「流入熱量/電力量」として定義できる空調効率は直接観測できないプロセス特性であるが、流入熱量は予測値計算部９において熱バランスを計算するときに求める物理量であるため、予測値計算部９の結果を利用して空調効率を推定することができる。
【００９２】
本実施形態では、計算モデルηの出力である室温を、そのまま仮想センサ値出力部１６の出力とする。
【００９３】
最も観測データを説明できるパラメータθが求まると、観測データyとη(x,θ)の差分がモデル化誤算δとなる。モデル化誤差とは、オフィスの例だと、例えば、午後からモデルには含まれない熱源が動作して室温を上げるという時系列パターンである。モデル化誤差は観測データがある箇所に離散的に記録される場合もあるし、連続的に記録される場合もある。モデル化誤差δは、モデル誤差記録部１３に記録される。
【００９４】
例えば、上述した本実施形態の予測値計算部９は、ＭＣＭＣを用いた場合には、観測を行った離散箇所(1…n)しかモデル化誤差δを計算できないが、ガウス過程を用いた場合ではδモデル化誤差(x)を連続的に関数として計算できる。
【００９５】
予測値計算部９により得られた予測値から仮想センサ値出力部１６により推定された仮想センサ値を求める場合、空間粒度が詳細であるとは限らない。再現詳細度回復部１２は、空間詳細度が離散的に配置された観測データを補間して空間詳細度を回復する。例えば本実施形態では、室温は１つの部屋を１つの点としてモデル化したものであり、部屋の温度は一様であると仮定したモデルになっている。再現詳細度回復部１２は、この推定された室温と、部屋に配置された少ない数の温度センサー値から空間温度分布を推定する。
【００９６】
本実施形態では、オフィスに設置した１４個の温度センサから、オフィス全体の３Ｄ温度分布を回復する。
【００９７】
図１０はオフィスの部屋に設置した温度センサの位置を２次元に射影した地図上に表示したものである。図１０は、オフィス内のセンサが全部で１４個ある例を示している。
【００９８】
図１１は、設置された１４個の温度センサから得られたある１日の温度変化をプロットしたグラフである。再現詳細度回復部１２は、再現詳細度回復として、ガウス過程を利用した回帰を行う。ガウス過程とは、関数ｙ(ｘ)上の確率分布として定義され、任意の点集合ｘ＿１，ｘ＿２，…ｘ＿Ｎに対するｙ(ｘ)の値の同時分布が、ガウス分布に従うものである。
【００９９】
ガウス過程による回帰の実現は、Ｎ個の観測データ(ｘ＿１，ｙ＿１)、(ｘ＿２，ｙ＿２),…, (ｘ＿Ｎ，ｙ＿Ｎ)が得られたとき、任意の点ｘ＿(Ｎ＋１)におけるｙの値の推定をｙ＿(Ｎ＋１)として、観測データのベクトルをｙ＿(Ｎ＋１)＝(ｙ＿１，…，ｙ＿ｎ，ｙ＿(ｎ＋１))とし、またΝ(ｘ｜μ，Σ) が、xが平均μ、共分散行列Σに従う、多変量ガウス分布であるとして、ｙ＿(Ｎ＋１)＝(ｙ＿１，…，ｙ＿ｎ，ｙ＿(ｎ＋１))の同時分布は、以下の（２５）式で表される。
【０１００】
ｐ(ｙN+1)＝Ｎ(ｙN+1｜０，ＣN+1)    …（２５）
【０１０１】
ここで、ＣN+1は、(Ｎ＋１)×(Ｎ＋１)の共分散行列であり、以下の（２６）式で表される。
ｃ(ｘn，ｘm)＝ｋ(ｘn，ｘm)＋β-1δnm    …（２６）
【０１０２】
ｋ（ｘn，ｘm)はカーネル関数と呼ばれ、近接するｘにおける値ｙも相関するであろうという性質に基づいて、近接相関具合を関数として表現したものである。
【０１０３】
この同時分布を分割して条件付きガウス分布を求めることができ、ＣN+1を以下の（２７）式のように分割する。
【数８】
2018049390.tif 000010
【０１０４】
ここで、ｋ＝(ｋ(ｘ1，ｘN+1)，ｋ(ｘ2，ｘN+1)，…，ｋ(ｘN，ｘN+1)であり、スカラー値ｃ＝ｋ(ｘN+1，ｘN+1)＋β-1である。ガウス分布の分割の公式より、観察データｙN＝(ｙ1，…，ｙn)が得られたときの条件付きの次のｙN+1の確率分布は、以下の（２８）式に示すガウス分布に従う確率分布として計算できる。
【数９】
2018049390.tif 000011
【０１０５】
すなわち、観測データにより、分布の中央関数と共分散行列が更新されることにより、観測データと整合性のあるガウス過程すなわち関数の確率分布が得られる。
【０１０６】
推定された中央関数ｋTＣN-1ｙは、直観的には既存の観測点からｘN+1への距離を重みとして、予測値を外挿して値を推定することに相当する。
【０１０７】
図１２はガウス過程回帰の一例を示す図である。１次元の場合の観測データは、以下の通りとしている。
【０１０８】
Ｘobs＝［-2.1, -1.5, 0.3, 1.8, 2.5]
Ｙobs＝[ 0.11579398, 0.67447562, 0.58653244, -0.93075898, -1.33088595]
【０１０９】
カーネル関数は以下の（２９）式で表される。
【数１０】
2018049390.tif 000012
【０１１０】
ここで、αとβは、推定する関数の形を決めるハイパーパラメータである。グラフでは、観測点、推定値を誤差としてプロットしたものである。
【０１１１】
２次元の場合は、観測データに対するカーネル関数は、上述した（２９）式と同様であり、２次元の分布を推定する。このとき、誤差すなわち分散も同時に推定される。
【０１１２】
以下、本実施形態によるガウス過程回帰を利用した詳細度回復の様子を示す。ガウス過程回帰では、ガウス分布における平均値として、予測値計算部９と仮想センサ値推定部１１で得られた室温の変化を用いる。
【０１１３】
図１３は、推定された部屋単位の温度の時間推移をプロットしたグラフである。室温を平均値関数とみたて、３次元空間に対応するカーネル関数を用いて、回帰を行った。予測結果は、室内の３次元グリッド上に配置した仮想センサでのセンサ値を計算し、これをプロットしたものが図１４である。
【０１１４】
図１４に示すように、本実施形態によれば、部屋単位の室温をベースとして、観測データを回帰して３次元の温度分布を回復することができる。
【０１１５】
本実施形態の予測値計算部９は、ＭＣＭＣまたはガウス過程を利用したキャリブレーションを行う例を説明したが、シミュレーション結果から未確定パラメータの確率分布を修正する他の手段を用いてもよい。例えば、カルマンフィルタを用いたパラメータの予測などが考えられる。
【０１１６】
このように、本実施形態では、センサデータと未確定パラメータの確率分布から、計算モデルによるシミュレーション結果がセンサデータと合致するまで、未確定パラメータの確率分布を修正するキャリブレーションを行う。次に、計算モデルによるシミュレーション結果とセンサデータとの誤差に基づいて、所望の仮想センサのセンサデータ（仮想センサ値）を推定する。次に、推定した仮想センサ値に対して、再現詳細度を回復した推定値の計算を行い、これに並行して推定値の信頼度を計算する。これにより、本実施形態によれば、設定領域内の特性値をセンサデータを利用して精度よく推定することができる。また、センサの数が少なくても、各センサのセンサデータを利用して、設定領域内の特性値の分布を精度よく推定できる。
【０１１７】
（第２の実施形態）
図１５は第２の実施形態による特性値推定装置１の概略構成を示すブロック図である。図１５の特性値推定装置１は、図１の構成に加えて、モデル修正部１７と、モデル修正出力部１８と、センサ配置計算部１９と、センサ配置出力部２０とを備えている。
【０１１８】
モデル修正部１７は、予測値計算部９が計算した誤差に基づいて、モデル切替部７が選択した計算モデルを修正する。モデル修正部１７の一例としては、ガウス過程でモデル誤差δが計算された場合には、予め計算モデルからモデル誤差δをガウス分布として除去した計算モデルを作成することが考えられる。
【０１１９】
モデル修正出力部１８は、モデル修正部１７が修正した計算モデルを出力して、修正した計算モデルを工場やプラント等の作業者に提示する。
【０１２０】
センサ配置計算部１９は、信頼度計算部１４が計算した信頼度に基づいて、予測値計算部９が計算した誤差をより小さくするためのセンサの配置プランを計算する。すなわち、センサ配置計算部１９は、信頼度計算部１４が計算した信頼度を利用して、信頼度を向上させるのに最適なセンサ設置箇所を計算する。センサ配置計算部１９の実現手法の一例は、最も信頼度の低い場所を特定し、センサを追加で設置する。
【０１２１】
センサ配置出力部２０は、センサ配置計算部１９が計算したセンサ配置プランを出力して、このセンサ配置プランを工場やプラント等の作業者に提示する。
【０１２２】
このように、第２の実施形態では、少ないセンサの追加で、特性値の予測精度を向上させることができる。また、予測値計算部９が計算した誤差に基づいて計算モデルを修正することで、誤差を小さくすることができる。
【０１２３】
上述した実施形態で説明した特性値推定装置の少なくとも一部は、ハードウェアで構成してもよいし、ソフトウェアで構成してもよい。ソフトウェアで構成する場合には、特性値推定装置の少なくとも一部の機能を実現するプログラムをフレキシブルディスクやＣＤ−ＲＯＭ等の記録媒体に収納し、コンピュータに読み込ませて実行させてもよい。記録媒体は、磁気ディスクや光ディスク等の着脱可能なものに限定されず、ハードディスク装置やメモリなどの固定型の記録媒体でもよい。
【０１２４】
また、特性値推定装置の少なくとも一部の機能を実現するプログラムを、インターネット等の通信回線（無線通信も含む）を介して頒布してもよい。さらに、同プログラムを暗号化したり、変調をかけたり、圧縮した状態で、インターネット等の有線回線や無線回線を介して、あるいは記録媒体に収納して頒布してもよい。
【０１２５】
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。
【符号の説明】
【０１２６】
１  特性値推定装置、２  モデル入力部、３  センサデータ入力部、４  確率分布入力部、５  モデル学習部、６  仮想センサ情報入力部、７  モデル切替部、８  確率分布記録部、９  予測値計算部、１０  確率分布修正部、１１  仮想センサ値推定部、１２  再現詳細度回復部、１３  モデル誤差記録部、１４  信頼度計算部、１５  信頼度出力部、１６  仮想センサ値出力部、１７  モデル修正部、１８  モデル修正出力部、１９  センサ配置計算部、２０  センサ配置出力部
【図１】
2018049390.tif 000013
【図２】
2018049390.tif 000014
【図３】
2018049390.tif 000015
【図４】
2018049390.tif 000016
【図５】
2018049390.tif 000017
【図６】
2018049390.tif 000018
【図７】
2018049390.tif 000019
【図８】
2018049390.tif 000020
【図９】
2018049390.tif 000021
【図１０】
2018049390.tif 000022
【図１１】
2018049390.tif 000023
【図１２】
2018049390.tif 000024
【図１３】
2018049390.tif 000025
【図１４】
2018049390.tif 000026
【図１５】
2018049390.tif 000027