特許 7444426 状態推定評価装置、方法、及び、プログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-02-27

(45)【発行日】2024-03-06

(54)【発明の名称】状態推定評価装置、方法、及び、プログラム

(51)【国際特許分類】

   G05B 23/02 20060101AFI20240228BHJP

   G06N 7/01 20230101ALI20240228BHJP

   G06N 20/00 20190101ALI20240228BHJP

【ＦＩ】

G05B23/02 Z

G06N7/01

G06N20/00

【請求項の数】 5

(21)【出願番号】P 2019166640

(22)【出願日】2019-09-12

(65)【公開番号】P2021043813

(43)【公開日】2021-03-18

【審査請求日】2022-07-22

(73)【特許権者】

【識別番号】506209422

【氏名又は名称】地方独立行政法人東京都立産業技術研究センター

(74)【代理人】

【識別番号】100200229

【弁理士】

【氏名又は名称】矢作 徹夫

(72)【発明者】

【氏名】金田 泰昌

(72)【発明者】

【氏名】鈴木 聡

【審査官】稲垣 浩司

(56)【参考文献】

【文献】特開２０１９－１３３４８６（ＪＰ，Ａ）

【文献】特開２０１６－１５７３２９（ＪＰ，Ａ）

【文献】特開２００１－２８２３０９（ＪＰ，Ａ）

【文献】特表２０１９－５１０２１５（ＪＰ，Ａ）

【文献】特開２０１８－０５５１６９（ＪＰ，Ａ）

【文献】特開平０７－２１８６１１（ＪＰ，Ａ）

【文献】特開２０１２－１８５１１１（ＪＰ，Ａ）

【文献】西山清，知識ベース 知識の森，１群５編６章カルマンフィルタ，日本，電子情報通信学会，2011年03月04日，https://www.ieice-hbkb.org/files/01/01gun_05hen_06m.pdf

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０５Ｂ ２３／０２

Ｇ０６Ｎ ７／０１

Ｇ０６Ｎ ２０／００

(57)【特許請求の範囲】

【請求項1】

推定対象である対象装置から観測された所定の時刻での観測値から、前記対象装置の前記所定の時刻での状態量を推定する状態推定評価装置であって、
前記所定の時刻で動作する値である状態量を動作点とし、前記動作点を条件値とする条件付き確率モデルを記憶する記憶部と、
新たな時刻での観測値を入力とし、前記条件付き確率モデルを用いて前記新たな時刻での状態量を推定し、かつ、前記新たな時刻での新たな動作点を出力する推定部と、
前記新たな動作点を入力とし、前記条件付き確率モデルの不確かさに関する確率的モデルを用いて前記新たな時刻での確信度を算出する算出部と、
を備え、
前記推定部が、相互に依存する観測更新部および時間更新部と、時間的な整合性を保つバッファを有し、
前記記憶部及び前記推定部がソフトセンサとして動作する状態推定評価装置。

【請求項2】

推定した前記所定の時刻での状態量及び前記確信度を表示する表示部をさらに備える請求項１に記載の状態推定評価装置。

【請求項3】

前記条件付き確率モデルを学習する学習部をさらに備える請求項１に記載の状態推定評価装置。

【請求項4】

推定対象である対象装置から観測された所定の時刻での観測値から、前記対象装置の前記所定の時刻での状態量を推定する状態推定評価方法であって、
新たな時刻での観測値を入力とし、前記所定の時刻で動作する値である状態量を動作点とし、前記動作点を条件値とする条件付き確率モデルを用いて前記新たな時刻での状態量を推定し、かつ、前記新たな時刻での新たな動作点を出力する推定ステップと、
前記新たな動作点を入力とし、前記条件付き確率モデルの不確かさに関する確率的モデルを用いて前記新たな時刻での確信度を算出する算出ステップと、を備え、
前記推定ステップが、相互に依存する観測更新ステップおよび時間更新ステップを有し、時間的な整合性を保つバッファを用い、
前記推定ステップがソフトセンサとして動作する状態推定評価方法。

【請求項5】

推定対象である対象装置から観測された所定の時刻での観測値から、前記対象装置の前記所定の時刻での状態量を推定する状態推定評価プログラムであって、
新たな時刻での観測値を入力とし、前記所定の時刻で動作する値である状態量を動作点とし、前記動作点を条件値とする条件付き確率モデルを用いて前記新たな時刻での状態量を推定し、かつ、前記新たな時刻での新たな動作点を出力する推定ステップと、
前記新たな動作点を入力とし、前記条件付き確率モデルの不確かさに関する確率的モデルを用いて前記所定の時刻での確信度を算出する算出ステップと、を備え、
前記推定ステップが、相互に依存する観測更新ステップおよび時間更新ステップを有し、時間的な整合性を保つバッファを用い、
前記推定ステップが、ソフトセンサとして動作するコンピュータに実行させる状態推定評価プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明の実施形態は、ソフトセンサを実現するために必要な状態推定を評価する技術に関する。

【背景技術】

【0002】

ソフトセンサとは、推定対象の装置である対象装置の状態量をリアルタイムに計測できるセンサがない場合に、推定技術を用いてソフト的に状態量のセンサを実現する技術である。例えば、化学プラントでは製品の品質管理を行う場合サンプリングによる分析をする必要がある。一般的に分析は時間や費用等のコストが大きくなってしまう問題があり、またリアルタイムでの管理が困難であるという問題がある。このような場合、ソフトセンサを用いて品質（対象装置の状態量に相当）を推定し、管理を行なっている。

【0003】

ソフトセンサは、操業データ（対象装置を簡易に観測できる観測値に相当）と製品の品質との相関関係を表すモデルを作成し、そのモデルを用いて品質を推定する。モデルは第一原理等の法則からホワイトボックス的に得られることは少ない。多くの場合、製造工程から得られたデータ（観測値に相当）をデータ解析することでブラックボックス的にモデルが作成されている。そのため、ソフトセンサで用いるモデルの妥当性や、新たに得られた操業データ（新たな点ともいう）に対するソフトセンサの出力値の妥当性をどのように正当化するかが問題となる。

【0004】

このような問題に対して特許文献１では、線形回帰モデルを基礎としたソフトセンサが実現されている。線形回帰モデルが基礎となっているため、説明変数に対してどのような結果となるかが想像しやすく、ソフトセンサの挙動の説明が行いやすい。そのため、解釈可能性が高いという意味でモデルの妥当性を正当化しやすい。

【0005】

特許文献２では、ソフトセンサの入力値と保存された複数のデータ（事例ベース）との類似度を算出することで、ソフトセンサの出力値や用いるモデルを評価している。

【0006】

特許文献３では、保存した複数のソフトセンサの出力値と、後日分析で得られた複数の実測値とで回帰分析を行い、その相関から点数化することでソフトセンサの出力値を評価している。回帰分析による回帰モデルを用いるため、補間効果によりデータが存在しない領域もカバーすることができる。

【先行技術文献】

【特許文献】

【0007】

【文献】特許５７０７２３０号公報

【文献】再公表特許ＷＯ０２／００６９５３

【文献】特開２００９－２３０２０９号公報

【発明の概要】

【発明が解決しようとする課題】

【0008】

しかしながら、特許文献１は、線形回帰モデルの妥当性を定量的に評価することができないため、推定した状態量も定量的に評価できないという課題がある。特許文献２は、保存された複数のデータとの類似度だけでソフトセンサを評価しているため、データが少ない領域ではソフトセンサの評価が下がってしまうという課題がある。特許文献３は、実測値を得るために分析が必要である。通常、分析結果はリアルタイムでは出力されず、予測した後日にその実測値が得られる。そのため、リアルタイムにソフトセンサの出力を評価することができないという課題がある。

【0009】

本発明は、このような課題に着目して鋭意研究され完成されたものであり、その目的は、ソフトセンサを実現するために必要な状態推定技術であって、推定した状態量の妥当性を定量的に評価するとともに、データが少ない領域でもソフトセンサを正しく評価し、リアルタイムに評価可能な技術を提供することにある。

【課題を解決するための手段】

【0010】

上記課題を解決するために、本発明は、推定対象である対象装置から観測された所定の時刻での観測値から、前記対象装置の前記所定の時刻での状態量を推定する状態推定評価装置であって、前記所定の時刻で動作する値である動作点を条件値とする条件付き確率モデルを記憶する記憶部と、前記観測値を入力とし、前記条件付き確率モデルを用いて前記所定の時刻での状態量を推定し、かつ、前記動作点を出力する推定部と、前記動作点を入力とし、前記条件付き確率モデルを用いて前記所定の時刻での確信度を算出する算出部と、を備える状態推定評価装置である。

【0011】

他の本発明は、推定対象である対象装置から観測された所定の時刻での観測値から、前記対象装置の前記所定の時刻での状態量を推定する状態推定評価方法であって、前記観測値を入力とし、前記所定の時刻で動作する値である動作点を条件値とする条件付き確率モデルを用いて前記所定の時刻での状態量を推定し、かつ、前記動作点を出力し、前記動作点を入力とし、前記条件付き確率モデルを用いて前記所定の時刻での確信度を算出する状態推定評価方法である。

【0012】

他の本発明は、推定対象である対象装置から観測された所定の時刻での観測値から、前記対象装置の前記所定の時刻での状態量を推定する状態推定評価プログラムであって、前記観測値を入力とし、前記所定の時刻で動作する値である動作点を条件値とする条件付き確率モデルを用いて前記所定の時刻での状態量を推定し、かつ、前記動作点を出力するステップと、前記動作点を入力とし、前記条件付き確率モデルを用いて前記所定の時刻での確信度を算出するステップと、をコンピュータに実行させる状態推定評価プログラムである。

【発明の効果】

【0013】

本発明によれば、ソフトセンサを実現するために必要な状態推定技術であって、推定した状態量の妥当性を定量的に評価するとともに、データが少ない領域でもソフトセンサを正しく評価し、リアルタイムに評価可能な技術を提供することができる。

【図面の簡単な説明】

【0014】

【図1】本発明の実施形態に係る状態推定評価装置の機能ブロック図である。

【図2】実施例１に係る状態推定器３の機能ブロック図である。

【図3】実施例１に係るモデル学習部６のフローチャートである。

【図4】実施例１に係る確信度算出部９のフローチャートである。

【図5】実施例２に係る状態推定器３０の機能ブロック図である。

【図6】実施例２に係る観測更新部１２のフローチャートである。

【図7】実施例２に係る時間更新部１３のフローチャートである。

【図8】実施例２に係るモデル学習部６のフローチャートである。

【図9】実施例２に係る確信度算出部９のフローチャートである。

【図10】実施例２に係る学習された推定対象の確率的モデルｐ（ｙ_ｋ｜ｘ_ｋ）の一例を示すグラフである。

【図11】実施例２に係るモデルの不確かさに関する確率的モデルｐ（ｈ^＊｜ｘ_ｋ）の一例を示すグラフである。

【図12】モデル化誤差が小さい場合のシミュレーションで用いる観測値の一部を示すグラフである。

【図13】モデル化誤差が小さい場合のシミュレーションによる推定結果一部を示すグラフである。

【図14】モデル化誤差が小さい場合のシミュレーションにより算出した状態推定器確信度の一部を示すグラフである。

【図15】モデル化誤差が大きい場合のシミュレーションで用いる観測値の一部を示すグラフである。

【図16】モデル化誤差が大きい場合のシミュレーションによる推定結果一部を示すグラフである。

【図17】モデル化誤差が大きい場合のシミュレーションにより算出した状態推定器確信度の一部を示すグラフである。

【発明を実施するための形態】

【0015】

図面を参照しながら本発明の実施の形態を説明する。ここで、各図において共通する部分には同一の符号を付し、重複した説明は省略する。また、図形は、長方形が処理部を表し、円柱がデータベースを表す。実線の矢印は処理の流れを表す。

【0016】

処理部及びデータベースは機能ブロック群であり、ハードウェアでの実装に限られず、ソフトウェアとしてコンピュータに実装されていてもよく、その実装形態は限定されない。例えば、パーソナルコンピュータ等のクライアント端末と有線又は無線の通信回線（インターネット回線など）に接続された専用サーバにインストールされて実装されていてもよいし、いわゆるクラウドサービスを利用して実装されていてもよい。

【0017】

（状態推定評価装置）
図１は、本発明の実施形態に係る状態推定評価装置の機能ブロック図である。状態推定評価装置１は、各データベースとして、確率的モデル記憶部５と、学習データ記憶部７を備える。また、状態推定評価装置１は、各処理部として、状態推定器３、モデル学習部６、確信度算出部９を備える。

【0018】

ここで、状態推定器３及び確率的モデル記憶部５はソフトセンサ１００（点線で表示）としても動作する。すなわち、状態推定器３は、推定対象の装置である対象装置に設けられたセンサ等によって観測された所定の時刻での観測値２を入力とし、確率的モデル記憶部５の確率的モデルを用いて対象装置の状態量４をソフトウェア的に推定し、所定の時刻での状態量４を出力する。

【0019】

状態推定評価装置１は、観測値２を入力とし、状態推定器３及び確率的モデル記憶部５を用いて状態量４を計算し、さらに確信度算出部９も用いて状態推定器確信度１０も計算する。また、状態推定評価装置１は、さらに表示部（不図示）を備え、推定した所定の時刻での状態量４と、状態推定器確信度１０を表示する表示部を備えていてもよい。状態推定評価装置１の利用者は推定した状態量の妥当性を定量的に、かつ、リアルタイムに評価することが可能になる。

【0020】

以下では、状態推定器３、モデル学習部６及び確信度算出部９の実施例を２つ説明する。また、確率的モデル記憶部５が記憶する確率的モデルも実施例毎に異なるため、確率的モデルについても数式を用いて説明する。

【0021】

（実施例１）
図２は、実施例１に係る状態推定器の機能ブロック図である。状態推定器３は、観測値２から、確率的モデル記憶部５の確率的モデルを用いて、期待値を算出する期待値算出部１１を備える。

【0022】

状態推定器３は、期待値算出部１１で算出された期待値を状態量４として出力する。また、状態量４を計算する際に用いた観測値２を動作点８として出力する。出力の時間間隔は約１秒である。

【0023】

期待値算出部１１では、確率的モデル記憶部５に保存された状態推定器ｇ（ｙ_ｋ）の確率的モデルを用いて、状態量４を計算する。具体的には、確率的モデル記憶部５に格納された状態推定器ｇ（ｙ_ｋ）の確率的モデルｐ（ｘ_ｋ｜ｙ_ｋ）の平均値関数を計算し、その平均値関数を状態推定器ｇ（ｙ_ｋ）とする。そして、その平均値関数から計算された値を状態量４として出力する。つまり、式（１）を計算し、そこから得られた値を状態量４として出力する。

【数1】

【0024】

ここで、確率的モデルｐ（ｘ_ｋ｜ｙ_ｋ）は、時刻ｋでの観測値ｙ_ｋが条件として与えられている時の状態量ｘ_ｋの条件付き確率モデルである。観測値ｙ_ｋは動作点ともいう。動作点とは、各時刻で動作している値（この場合、観測値ｙ_ｋ）を中心に、物事（この場合、状態量ｘ_ｋ）を考える場合の概念である。

【0025】

図３は、実施例１に係るモデル学習部のフローチャートである。モデル学習部６は、状態推定器３と確信度算出部９とで用いる確率的モデル記憶部５の確率的モデルを学習する（Ｓ１１）。具体的な学習方法は以下の通りである。今、時刻ｋの状態量をｘ_ｋ、観測値をｙ_ｋとする。学習データ記憶部７の学習データとしてＤ＝{（ｘ_１，ｙ_１），（ｘ_２，ｙ_２），・・・，（ｘ_Ｎ，ｙ_Ｎ）}が与えられているとき、状態推定器３における状態推定器ｇ（ｙ_ｋ）の確率的モデルを学習する。事前分布として状態推定器ｇ（y_ｋ）の出力値列ｇ＝[ｇ（y_１）・・・ｇ（y_Ｎ）]^Ｔの分布ｐ（ｇ）を式（２）で与える。

【数2】

【0026】

ここで、Ｎ（μ，Σ）は平均μ、共分散Σの正規分布を表す。また、Ｋはグラム行列を表す。グラム行列Ｋは観測値列ｙ＝{ｙ_１・・・ｙ_Ｎ}から計算され、各観測値間の分散に相当する値を要素とする行列である。今、状態推定器ｇ（ｙ_ｋ）の出力に加法的にノイズ（例えば、観測ノイズ）が加わるとし、その結果得られる状態量がｘ_ｋであると考える。そして、その加法ノイズがＮ（０，σ^２Ｉ）に従うとすると、ｇに対応した状態量４の観測値列ｘ＝{ｘ_１・・・ｘ_Ｎ}の分布は式（３）で計算できる。

【数3】

【0027】

次に、新たな観測値ｙ^＊に対する状態量ｘ^＊の分布を計算する。ｘ^＊とｘの同時分布は式（４）で計算できる。

【数4】

【0028】

ここでｋ_＊およびｋ_＊＊は観測値列ｙ＝{ｙ_１・・・ｙ_Ｎ}および新たな観測値ｙ^＊から計算されるグラム行列の新たな要素である。よって、ｘ^＊の分布は式（５）として得られる。

【数5】

【0029】

したがって、新たな観測値ｙ^＊に対する状態量ｘ^＊はｘ^＊＝ｋ_＊ ^Ｔ（Ｋ＋σ^２Ｉ）^－１ｘで与えられる。以上より、状態推定器ｇ（ｙ_ｋ）の確率的モデルを、例えば式（６）で与えることができる。

【数6】

【0030】

そして、得られた確率的モデル（式（６））を確率的モデル記憶部５へ保存する（Ｓ１２）。状態推定器の確率的モデル（式（６））はガウス分布で与えられているため、式（１）は式（７）の通り容易に計算できる。

【数7】

【0031】

よって、状態推定器（式（７））から、状態量４が計算できる。同様に、確信度算出部９で用いるための状態推定器ｇ（ｙ_ｋ）の不確かさに関する確率的モデル（条件付き確率モデル）をガウス分布として算出することができる。今、新たな観測値ｙ^＊に対する状態推定器ｇ（ｙ_ｋ）の値ｇ^＊の分布を計算する。状態推定器ｇ（ｙ_ｋ）の確率的モデル（式（６））の学習と同様の流れより、ｇ^＊とｘの同時分布は式（８）で計算できる。

【数8】

【0032】

よって、ｇ^＊の分布は式（９）として得られる。

【数9】

【0033】

すなわち、新たな観測値ｙ^＊に対する状態推定器ｇ（ｙ_ｋ）の値ｇ^＊の分散がｋ_＊＊＋ｋ_＊ ^Ｔ（Ｋ＋σ^２Ｉ）^－１ｋ_＊で与えられており、状態推定器ｇ（ｙ_ｋ）の不確かさに関する確率的モデルがガウス分布として計算できる。そして、得られた確率的モデル（式（９））を確率的モデル記憶部５へ保存する。

【0034】

なお、本実施例では状態推定器ｇ（ｙ_ｋ）の確率的モデル（式（６））と、状態推定器ｇ（ｙ_ｋ）の不確かさに関する確率的モデル（式（９））を分離して求めている。しかしながら、式（５）を状態推定器ｇ（ｙ_ｋ）の確率的モデルおよび状態推定器ｇ（ｙ_ｋ）の不確かさに関する確率的モデルとしてもよい。

【0035】

図４は、実施例１に係る確信度算出部９のフローチャートである。確信度算出部９は、動作点８と確率的モデル記憶部５に保存された状態推定器ｇ（ｙ_ｋ）の不確かさに関する確率的モデル（式（９））とを用いて、状態推定器確信度１０を計算する。具体的には、状態推定器ｇ（ｙ_ｋ）の不確かさに関する確率的モデル（式（９））の分散を計算し、その分散の逆数を状態推定器確信度１０とする。状態推定器ｇ（ｙ_ｋ）の不確かさに関する確率的モデル（式（９））はガウス分布で与えられているため、その分散は式（１０）の通り容易に計算できる（Ｓ２１）。

【数10】

【0036】

よって、状態推定器確信度１０は式（１１）で計算できる（Ｓ２２）。

【数11】

【0037】

（実施例２）
図５は、実施例２に係る状態推定器の機能ブロック図である。状態推定器３０は、観測更新部１２、時間更新部１３、及び、バッファ１４を備える。ここで、実施例１は推定対象である対象装置のデータを使って状態推定器３を学習（すなわち、観測値yから状態量xの関数を学習）する。これに対し、実施例２は、対象装置を学習（状態量xから観測値yの関数を学習）し、それを用いて観測更新部１２および時間更新部１３を計算することで、状態推定器３０を実現する。

【0038】

状態推定器３０は、観測更新部１２で計算された値を状態量４および動作点８として出力する。観測更新部１２では時間更新部１３で計算された結果を用い、また、時間更新部１３では観測更新部１２で計算された結果を用いる。このように、観測更新部１２と時間更新部１３は相互に依存しており、バッファ１４を用いることで代数ループを防ぎ、時間的な整合性を保っている。

【0039】

図６は、観測更新部１２のフローチャートである。観測更新部１２は、状態量４の複数の候補値（以下、粒子という）に関して、時間更新部１３で計算された粒子を、センサ等によって観測された観測値２と確率的モデル記憶部５に保存された確率的モデルに基づき更新する。そして、更新された粒子から状態量４を推定する。ここでは、粒子は１００個とし、その時間間隔は約１秒とし、粒子１００個の平均が状態量になる。

【0040】

具体的な計算方法は以下の通りである。今、時刻ｋまでに得られた観測値列をＹ_ｋ＝{ｙ_１・・・ｙ_ｋ}、時刻ｋの状態をｘ_ｋ、時間更新部１３で計算された時刻ｋの状態の粒子をｘ_{ｋ｜ｋ－１} ^（ｉ）、時間更新部１３で計算された時刻ｋの状態の分布をｐ（ｘ_ｋ｜Ｙ_ｋ－１）、確率的モデル記憶部５に保存されている確率的モデルをｐ（ｙ_ｋ｜ｘ_ｋ）とする。このとき、観測更新部１２では、式（１２）のフィルタ分布を計算する。

【数12】

【0041】

ここで、ディラックのデルタ関数を用いて時間更新部１３で計算された状態の分布ｐ（ｘ_ｋ｜Ｙ_ｋ－１）を近似し、粒子ｘ_{ｋ｜ｋ－１} ^（ｉ）に対して式（１３）を用いて尤度を計算する（Ｓ３１）。

【数13】

【0042】

式（１３）は、フィルタ分布（式（１２））がｘ_{ｋ｜ｋ－１} ^（ｉ）と尤度ω_ｉで近似されていることを意味する。このとき、状態量４をｘ＾_ｋとすると、状態量４は式（１４）の重み付き平均で計算できる。また、観測ノイズに対する状態量４の推定誤差は式（１５）の推定誤差共分散行列で与えられる（Ｓ３２）。なお、本来は、ハット記号「^」は「ｘ」の上に記載されるものであるが、明細書で使用可能な文字コードの都合上「ｘ」と「^」を並べて記載する。

【数14】

【数15】

【0043】

図７は、実施例２に係る時間更新部１３のフローチャートである。時間更新部１３は、現時刻の粒子から、確率的モデル記憶部５に保存された確率的モデルに基づいて時刻１ステップ先の粒子を計算する。具体的な計算方法は以下の通りである。確率的モデル記憶部５に保存されている確率的モデルをｐ（ｘ_ｋ＋１｜ｘ_ｋ）とする。このとき、時間更新部１３では、式（１６）の予測分布を計算する（Ｓ４１）。

【数16】

【0044】

つまり、予測分布が、観測更新部１２で計算された尤度ω_ｉを重みとした混合分布で近似されていることを意味する。この混合分布（式（１６））からサンプリングすることで（Ｓ４２）、時刻ｋの１ステップ先の粒子ｘ_{ｋ＋１｜ｋ} ^（ｉ）が計算できる。

【0045】

図８は、実施例２に係るモデル学習部６のフローチャートである。モデル学習部６は、観測更新部１２および時間更新部１３で用いる確率的モデル記憶部５に保存された確率的モデルｐ（ｙ_ｋ｜ｘ_ｋ）およびｐ（ｘ_ｋ＋１｜ｘ_ｋ）を学習する。具体的な学習方法は以下の通りである。今、学習データ７に保存された学習データＤ＝{（ｘ_１，ｙ_１），（ｘ_２，ｙ_２），・・・，（ｘ_Ｎ，ｙ_Ｎ）}を用いて、ｙ_ｋ＝ｈ（ｘ_ｋ）＋ｖ_ｋで表される対象装置（システムともいう）の確率的モデルｐ（ｙ_ｋ｜ｘ_ｋ）を学習する。ここでｖ_ｋは観測ノイズである。事前分布として関数ｈ（ｘ_ｋ）の出力値列ｈ＝{ｈ（ｘ_１）・・・ｈ（ｘ_Ｎ）}の分布ｐ（ｈ）を式（１７）で与える。

【数17】

【0046】

ここで、Ｎ（μ，Σ）は平均μ、共分散Σの正規分布を表す。また、Ｋはグラム行列を表す。グラム行列Ｋは状態量４の観測値列Ｘ_Ｎ＝{ｘ_１・・・ｘ_Ｎ}から計算され、各状態量間の分散に相当する値を要素とする行列である。観測ノイズｖ_ｋがＮ（０，σ^２Ｉ）に従うとすると、ｈに対応した観測値列ｙ＝{ｙ_１・・・ｙ_Ｎ}の分布は式（１８）で計算できる。

【数18】

【0047】

次に、新たな点ｘ^＊に対する観測値ｙ^＊の分布を計算する。ｙ^＊とｙの同時分布は式（１９）で計算できる。

【数19】

【0048】

よって、新たな観測値ｙ^＊の分布は式（２０）として得られる。

【数20】

【0049】

したがって、点ｘ^＊に対する新たな観測値ｙ^＊はｙ^＊＝ｋ_＊ ^Ｔ（Ｋ＋σ^２Ｉ）^－１ｙで与えられる。以上より、ｙ_ｋ＝ｈ（ｘ_ｋ）＋ｖ_ｋで表されるシステムの観測ノイズによる確率的モデルｐ（ｙ_ｋ｜ｘ_ｋ）を、例えば式（２１）で与えることができる。

【数21】

そして、得られた確率的モデル（２１）式を確率的モデル記憶部５に保存する（Ｓ５１）。

【0050】

なお、学習データ７に保存されたデータをＤ＝{（ｘ_１，ｘ_２），（ｘ_２，ｘ_３），・・・，（ｘ_Ｎ－１，ｘ_Ｎ）}とすることで、ｘ_ｋ＋１＝ｆ（ｘ_ｋ）＋ｗ_ｋで表されるシステムの確率的モデルｐ（ｘ_ｋ＋１｜ｘ_ｋ）の学習が、モデル学習部６において同様に行える。ここで、ｗ_ｋは状態ノイズである。そして、学習された確率的モデルを確率的モデル記憶部５に保存する。

【0051】

同様に、モデルの不確かさに関する確率的モデルをガウス分布として算出することができる。今、推定された状態量ｘ^＊に対する関数値ｈ^＊の分布を計算する。推定対象の確率的モデルの学習と同様の流れより、ｈ^＊とｙの同時分布は式（２２）で計算できる。

【数22】

【0052】

ここでｋ_＊およびｋ_＊＊は状態量４の観測値列Ｘ_Ｎ＝{ｘ_１・・・ｘ_Ｎ}および推定された状態量ｘ^＊から計算されるグラム行列の新たな要素である。よって、ベイズの定理より、ｈ^＊の分布は式（２３）として得られる。

【数23】

【0053】

すなわち、推定された状態量ｘ^＊に対するｈ^＊の分散がｋ_＊＊＋ｋ_＊ ^Ｔ（Ｋ＋σ^２Ｉ）^－１ｋ_＊で与えられており、モデルの不確かさに関する確率的モデルがガウス分布として計算できる（S５２）。そして、得られたモデルの不確かさに関する確率的モデル（式（２３））を確率的モデル記憶部５に保存する（Ｓ５３）。

【0054】

なお、推定された状態量ｘ^＊に対する関数値ｆ^＊の分布も同様に計算できる。そして、モデルfの不確かさに関する確率的モデルがガウス分布として計算できる。

【0055】

本実施例では、推定対象の確率的モデルとモデルの不確かさに関する確率的モデルを別に保存した。しかしながら、モデルの不確かさに関する確率的モデルを推定対象の確率的モデルとしてもよい。

【0056】

図９は、実施例２に係る確信度算出部９のフローチャートである。確信度算出部１０では、動作点８と確率的モデル記憶部５に保存された確率的モデルの不確かさに関する確率的モデル（式（２３））とを用いて、状態推定器確信度１０を計算する。具体的には、モデルの不確かさに関する確率的モデル（式（２３））の分散を計算し、その分散の逆数を状態推定器確信度１０とする。モデルの不確かさに関する確率的モデル（式（２３））はガウス分布で与えられているため、その分散は式（２４）の通り容易に計算できる（Ｓ６１）。

【数24】

【0057】

よって、状態推定器確信度１０は次式で計算できる（Ｓ６２）。

【数25】

【0058】

（検証）
本実施形態に係る状態推定評価方法の効果について検証するため、式（２６）、式（２７）で表される非線形システムの状態推定を行う。状態推定器の実施形態として実施例２（図５の状態推定器３０の機能ブロック図）を用いる。

【数26】

【数27】

【0059】

ここで、状態ノイズｗ_ｋがＮ（０，０．０４）に、観測ノイズｖ_ｋがＮ（０，０．０４）にそれぞれ従うとする。

【0060】

－３０≦ｘ_ｋ≦０の内、２０００個のｘ_ｋに対応したｘ_ｋ＋１およびｘ_ｋを式（２６）および式（２７）から計算し、確率的モデルを学習するためのデータとして利用する。

【0061】

図１０は、このデータを用いて学習された推定対象の確率的モデルｐ（ｙ_ｋ｜ｘ_ｋ）を示す。図１１は、モデルの不確かさに関する確率的モデルｐ（ｈ^＊｜ｘ_ｋ）を示す。

【0062】

図１０および図１１平面に描かれた実線の曲線（ｔｒｕｅ）はｖ_ｋ＝０としたときの式（２６）のグラフを、破線の曲線は学習した確率的モデルより計算されたグラフを表す。

【0063】

図１０の縦方向に描かれた破線（ｄｅｎｓｉｔｙ ｏｆ ｎｏｉｓｅ）は、代表的な３点（ｘ_ｋ＝－４．０，０．０，４．０）において、学習した推定対象の確率的モデルより計算された分布を表す。すなわち、観測値に含まれるノイズに関する分布を表現しており、観測誤差に対する影響を考慮したモデルとなっている。また、図１１縦方向に描かれた破線（ｄｅｎｓｉｔｙ ｏｆ ｍｏｄｅｌｌｉｎｇ ｅｒｒｏｒ）は、代表的な３点（ｘ_ｋ＝－４．０，０．０，４．０）において、学習したモデルの不確かさに関する確率的モデルより計算された分布を表す。

【0064】

今、学習データは－３０≦ｘ_ｋ≦０の範囲から取得しているため、０≦ｘ_ｋの範囲ではモデル化誤差が大きくなるはずである。実際、図１１では、０≦ｘ_ｋにおけるモデル化誤差が大きくなっている。すなわち、実線（ｔｒｕｅ）と太い点線（ｍｏｄｅｌ）が乖離している。

【0065】

また、そのときの分布の幅が非常に大きくなっている。すなわち、細い点線（ｄｅｎｓｉｔｙ ｏｆ ｍｏｄｅｌｌｉｎｇ ｅｒｒｏｒ）はピークが低く、ばらつきが大きい。つまり、モデル化誤差の大きさを分布の幅によって表現している。

【0066】

まず、モデル化誤差が小さい場合として、初期値をｘ_０＝－４．０とし、粒子の数を１００としてシミュレーションを行う。図１２は、シミュレーションで用いる観測値（ｍｅａｓｕｒｅｍｅｎｔ）の一部を示すグラフである。図１３は、シミュレーションによる推定結果の一部を示すグラフである。図１４は、シミュレーションにより算出した状態推定器確信度（ｃｏｎｆｉｄｅｎｃｅ）の一部を示すグラフである。図１１より、ｘ_０＝－４．０の近傍はモデル化誤差が少ないため、精度のよい推定結果が得られている（図１３参照）。また、それに対応して状態推定器確信度も高い数字となっている（図１４参照）。このように、データが少ない領域でも新たな点に対するソフトセンサの出力値を正しく評価することが可能である。

【0067】

次に、モデル化誤差が大きい場合として、初期値をｘ_０＝４．０とし、粒子の数を１００としてシミュレーションを行う。図１５は、シミュレーションで用いる観測値（ｍｅａｓｕｒｅｍｅｎｔ）の一部を示すグラフである。図１６は、シミュレーションによる推定結果の一部を示すグラフである。図１７は、シミュレーションにより算出した状態推定器確信度（ｃｏｎｆｉｄｅｎｃｅ）の一部を示すグラフである。図１１より、ｘ_０＝４．０の近傍はモデル化誤差が大きいため、推定精度が劣化していることがわかる（図１６参照）。また、それに対応して状態推定器確信度も低い数字となっている（図１７参照）。

【0068】

（作用効果）
本実施形態の状態推定評価装置によれば、対象装置（本検証では、非線形システム）を実際に計測した観測値と推定対象である状態量との相互相関式を確率分布（例えばガウス分布）でモデル化している。そして、推定量の確率分布からばらつきを求め、信頼度を分散の逆数として計算している。このため、ばらつきが小さい場合、状態推定器の確信度が高くなる。

【0069】

シミュレーションによる検証によれば、確率分布によるモデル化誤差が小さいところは分布の幅が狭く、対象装置の観測値から、対象装置の状態量の真値を推定できており、また確信度の数値が相対的に大きいことがわかる。また、モデル化誤差が大きいところは分布の幅が広く、真値が推定できておらず、また信頼度の数値が相対的に小さいことがわかる。

【0070】

したがって、本実施形態の状態推定評価装置によれば、ソフトセンサを実現するために必要な状態推定技術であって、推定した状態量の妥当性を定量的に評価するとともに、データが少ない領域でもソフトセンサを正しく評価し、リアルタイムに評価することが可能である。

【0071】

以上、本発明の実施例（変形例を含む）について説明してきたが、これらのうち、２つ以上の実施例を組み合わせて実施しても構わない。あるいは、これらのうち、１つの実施例を部分的に実施しても構わない。さらには、これらのうち、２つ以上の実施例を部分的に組み合わせて実施しても構わない。

【0072】

また、本発明は、上記発明の実施例の説明に何ら限定されるものではない。特許請求の範囲の記載を逸脱せず、当業者が容易に想到できる範囲で種々の変形態様もこの発明に含まれる。例えば、非線形システムの一例として、リチウムイオン電池の健全性推定を行ってもよい。この場合、簡易計測可能な電流及び電圧を実際に計測し、電流で電圧を正規化した「電圧／電流」を観測値ｙ_ｋとし、電池の劣化状態（Ｓｔａｔｅ Ｏｆ Ｈｅａｌｔｈ、ＳＯＨ）を状態量ｘ_ｋとすればよい。また、他の電池の劣化診断のみならず、印刷のインク配合推定や、めっき廃液管理に用いてもよい。

【符号の説明】

【0073】

１ 状態推定評価装置
３、３０ 状態推定器
１００ ソフトセンサ

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】