JP7531450B2 - 区分線形モデル作成装置、区分線形モデル作成方法およびプログラム - Google Patents

Description

本開示は、区分線形モデル作成装置、区分線形モデル作成方法およびプログラムに関する。

火力発電所などのプラントやエンジンなどでは、弁などの非線形性の強い機器が用いられる。この機器の状態推定には、カルマンフィルタなどのデータ同化手法が用いられる。一方、カルマンフィルタは状態方程式で記述可能な線形モデルしか扱うことができず、非線形性の強い機器のデータ同化を扱うことが難しい。特許文献１には、非線形性を有する機器の状態推定を、非線形カルマンフィルタによって実現する技術が開示されている。しかし、非線形カルマンフィルタによる状態推定は計算負荷が高い。非線形性を有する機器の特性を記述する他の方法として、区分線形モデルを用いる方法が提供されている（非特許文献１）。この方法は、機器の特性を示す曲線を区分ごとに一次関数で近似してモデル化する方法である。

特開２０１９－２７３２９号公報

Rafael Rui, Tohid Ardeshiri and Alexandre Bazanella, "Identification of Piecewise Affine State-Space Models via Expectation Maximization", 2016［閲覧：2021年4月5日］ＵＲＬ：https://arxiv.org/abs/1609.09789

区分線形モデルでは、各区分それぞれにおいて独立にデータ同化する。その結果、区分線形モデルの連続性が損なわれる場合がある。例えば、機器特性を示すある１つの点の値が突出して大きくなったり小さくなったりすると、その点を節点とする区分の一次関数の傾きが突出した点に引きずられて変化し、区分線形モデルの形状に尖った部分が現れ、滑らかではなくなることがある。すると、機器特性が本来滑らかな曲線で表される場合、機器特性の記述力が損なわれることになる。

本開示は、上記課題を解決することができる区分線形モデル作成装置、区分線形モデル作成方法およびプログラムを提供する。

本開示の区分線形モデル作成装置は、区分線形モデルにおける区分を規定する節点の真値をカルマンフィルタによって推定するにあたり、前記区分線形モデルの滑らかさに関するペナルティ項をカルマンゲインの算出式に追加して、カルマンゲインおよび誤差共分散行列を算出して、前記節点の真値を推定し、推定された前記節点を結んで区分線形モデルを作成する計算部、を備える。

本開示の区分線形モデル作成方法は、区分線形モデルにおける区分を規定する節点の真値をカルマンフィルタによって推定するにあたり、前記区分線形モデルの滑らかさに関するペナルティ項をカルマンゲインの算出式に追加して、カルマンゲインおよび誤差共分散行列を算出して前記節点の真値を推定し、推定された前記節点を結んで区分線形モデルを作成する。

本開示のプログラムは、コンピュータに区分線形モデルにおける区分を規定する節点の真値をカルマンフィルタによって推定するにあたり、前記区分線形モデルの滑らかさに関するペナルティ項をカルマンゲインの算出式に追加して、カルマンゲインおよび誤差共分散行列を算出して前記節点の真値を推定し、推定された前記節点を結んで区分線形モデルを作成する処理、を実行させる。

上述の区分線形モデル作成装置、区分線形モデル作成方法およびプログラムによれば、滑らかさを考慮した区分線形モデルを作成することにより、機器の非線形特性を高精度に推定・記述することができる。

実施形態の区分線形モデル作成装置の一例を示すブロック図である。 実施形態に係る区分線形モデルの一例を示す図である。 実施形態に係る局所フィッティングを説明するための図である。 実施形態に係るカーブフィッティングを説明するための図である。 カルマンフィルタに用いる式を示す図である。 実施形態に係るカルマンフィルタの拡張について説明する第１の図である。 実施形態に係るカルマンフィルタの拡張について説明する第２の図である。 実施形態に係るフィッティング処理の一例を示すフローチャートである。 実施形態に係る区分線形モデル作成処理の一例を示すフローチャートである。 実施形態に係る区分線形モデル作成装置のハードウェア構成の一例を示す図である。

＜実施形態＞
以下、本開示の実施形態に係る区分線形モデル方法について、図１～図９を参照しながら説明する。
（構成）
図１は、実施形態の区分線形モデル作成装置の一例を示すブロック図である。
図示するように、区分線形モデル作成装置１０は、データ取得部１１と、入力受付部１２と、計算部１３と、出力部１４と、記憶部１５と、を備える。
データ取得部１１は、観測データを取得する。観測データとは、例えば、プラントが備える機器の特性を示すデータである。
入力受付部１２は、閾値や係数などの設定情報やユーザによる指示操作を受け付ける
計算部１３は、区分線形モデルを作成する処理を実行する。計算部１３は、区分線形モデルの作成に際し、カーブフィッティングと局所フィッティングを使い分けて、精度良く機器特性の曲線を推定する。カーブフィッティングと局所フィッティングについては後述する。
出力部１４は、計算部１３が作成した区分線形モデルをグラフなどに表示して、表示装置や電子ファイルに出力する。
記憶部１５は、データ取得部１１が取得した観測データ、計算部１３が作成した区分線形モデル、閾値や係数など各種データを記憶する。

本実施形態の区分線形モデル作成装置１０は、非線形な特性を有する対象を区分線形モデルによってモデル化する。非線形な特性の一例として、本実施形態ではバルブの流量特性を取り上げる。バルブの流量特性とは、例えば、Ｍｖ値（弁開度など）とＣｖ値の関係である。Ｃｖ値とは流体の流れやすさを示す値である。このＭｖ値とＣｖ値の関係が示すバルブの流量特性は非線形（例えば、バルブを２倍の開度で開いたときに流量が今の２倍となるという関係が成立しない。）であることが知られている。このような非線形な特性を記述するのに区間線形モデルが用いられることがある。図２を参照する。図２の縦軸はＭｖ値、横軸はＣｖ値である。ここで、図２に示す点Ｐ１～Ｐ９は、例えば、バルブをあるＭｖ値で制御したときに実際に観測されたＣｖ値である。観測点Ｐ１～Ｐ９を順に結んでできた線は、バルブの流量特性を示す区分線形モデルである。つまり、観測データが得られなかった範囲については点同士を結んでできた線上の点がこのバルブの流量特性を表していると考える。以降、この区分線形モデルを関数Ｆｘと呼ぶ。また、区分線形モデルの各区分を規定する点Ｐ１～Ｐ９のことを節点と呼ぶ。

一般に観測データにはノイズが重畳しているため、実際には、図２の枠２Ａに示すように、区分線形モデルで定義された機器特性に対して得られる観測データは、区分線形モデルの線を中心に値が分布する。これに対し、機器特性の真値を求める必要がある。例えば、カルマンフィルタを用いて、観測データの真値を推定する。真値とは、Ｍｖ値とＣｖ値の組合せを推定対象の状態量としたときに、推定誤差が最小となるときのＭｖ値とＣｖ値の組合せのことである。観測データが得られると、カルマンフィルタにより真値（つまり図２のグラフ上の点）を推定し、例えば、これまで真値とみなされていた点Ｐ６に代えて、新たに真値と推定された点Ｐ６´を採用し、点Ｐ５と点Ｐ６´、点Ｐ６´と点Ｐ７を結んで関数Ｆｘを更新する。これにより、実際の観測データに応じて適応的に区分線形モデル（関数Ｆｘ）を更新することができる。バルブなどは経年劣化することから、流量特性が年数を経るとともに変化するが、この方法であれば、運用中に計測された観測データから現在の流量特性を表す区分線形モデルを得ることができる。また、この方法には、特許文献１に開示されている非線形カルマンフィルタのように膨大な計算量を必要とせず、例えば、プラントの制御装置に搭載されたコンピュータでもオンラインで実行可能な計算負荷で対応することができる。

しかし、一般的なカルマンフィルタによって推定された真値に従って、区分線形モデルを更新すると、本当のバルブの流量特性から外れてしまう場合がある。ここで、図３を参照する。図３の縦軸はＭｖ値、横軸はＣｖ値である。図３に示す点Ｐ１～Ｐ７は、ある時刻までに計測されたＭｖ値とＣｖ値の観測データから推定された真値である。点Ｐ１～Ｐ７について隣り合う点同士を結ぶことで関数Ｆｘが得られる。本来、バルブの流量特性は、滑らかな曲線で表され、関数Ｆｘはこの曲線を良好に近似しているとする。バルブの運用中に、Ｍｖ値とＣｖ値の観測データが蓄積され、一般的なカルマンフィルタにより、点Ｐ４に代わる新たな点Ｐ４´が推定されたとする。すると、一般的な区分線形モデルの更新では、新たに推定された真値の周辺だけを局所的に更新する。具体的には、点Ｐ３と点Ｐ４を結ぶ線を点Ｐ３と点Ｐ４´を結ぶ線に置き換え、点Ｐ４と点Ｐ５を結ぶ線を点Ｐ４´と点Ｐ５を結ぶ線に置き換えて関数Ｆｘを更新する。更新された箇所を破線にて示す。このように更新すると、点Ｐ３～点Ｐ５の間が鋸状となる。バルブの流量特性が滑らかな曲線で表されることを考えると、更新後の関数Ｆｘによって示される点Ｐ３～点Ｐ５の間の流量特性は精度が十分でない可能性が高い。つまり、このように局所的なフィッティングによって得られる関数Ｆｘでは、点Ｐ３～点Ｐ５の間の機器特性の記述力が損なわれていることになる。このような問題に対し、本実施形態では、機器特性が示す滑らかな曲線を維持したまま区分線形モデルを作成・更新する方法を提供する。以下、本実施形態の区分線形モデルの作成・更新方法をカーブフィッティング、図３を用いて説明した区分線形モデルの作成・更新方法を局所フィッティングと呼ぶ。

図４に、更新前の関数Ｆｘを実線で示し、カーブフィッティングによって更新した関数Ｆｘ´を破線で示す。本実施形態のカーブフィッティングでは、カルマンフィルタの誤差共分散行列の算出時に、区分線形モデルの滑らかさのペナルティ項を導入する。より具体的には、更新前の関数Ｆｘと更新後の関数Ｆｘ´の差の一階微分および二階微分の少なくとも１つについて、更新前の関数Ｆｘと更新後の関数Ｆｘ´のノルム（二乗ノルム）をペナルティ項として追加して誤差共分散行列を算出することで、観測データとの誤差を最小にし、かつ滑らかさを維持可能な区分線形モデルを算出する。この計算により、局所フィッティングでは点Ｐ４´だけが突出して移動していたのに対し、図４に示すカーブフィッティングでは、点Ｐ３と点Ｐ６が、観測データとの誤差を最小にしつつ、滑らかさを維持するようにそれぞれ点Ｐ３´と点Ｐ６´へ移動し、その結果、関数Ｆｘ´が得られる。ペナルティ項のうち一階微分の項は、元の関数Ｆｘの形状を維持したまま、関数Ｆｘが点Ｐ４から点Ｐ４´を含むような位置へ並行移動させるように働く。二階微分の項は、元の関数Ｆｘの形状（曲率）を、全ての観測データの分布に沿った形状・曲がり具合へ調整するように働く。なお、ペナルティ項には任意の大きさの重み付けを付すことができる。

図５に一般的なカルマンフィルタで用いられる公式を示す。式（１）はカルマンゲインＫの導出に用いられる式である。式（１）は、更新後の誤差共分散行列ＰをカルマンゲインＫで微分した式を表し、カルマンフィルタでは、式（１）の値を０にするようなカルマンゲインＫが導出される。具体的には、カルマンゲインＫは式（２）により計算される。また、更新後の誤差共分散行列Ｐは、式（３）で表される。各式に用いられる記号の意味は枠５Ａ内に示してある。誤差共分散行列Ｐは、関数Ｆｘ（Ｍｖ値とＣｖ値の組みのデータ）の推定値の誤差の分布を表し、例えば、この分布の中心を真値として推定する。

次に図６Ａ、図６Ｂを参照してペナルティ項について説明する。図６Ａの上側に関数Ｆｘの一例を示す。関数Ｆｘの点（Ｃｖ値、Ｍｖ値）を（Ｃｉ、Ｍｉ）、（Ｃｉ＋１、Ｍｉ＋１）、・・・、等で表すと、関数Ｆｘの一階微分、二階微分（Ｍｖ値で微分）は、それぞれ図６Ａの式（４）、式（５）で表すことができる。また、一階微分用の演算子Ｊを図６Ｂの式（６）のようにおくと、関数Ｆｘの一階微分はＪｘのように記述できる。関数Ｆｘの二階微分についても行列形式の二階微分用の演算子Ｈ（図示せず）を用いて簡単な形式で記述できる。ここで、更新前の関数Ｆｘと更新後の関数Ｆｘ´の差の一階微分を考える。具体的には、一階微分領域における、更新前の関数Ｆｘに対する更新後の関数Ｆｘのノルムの二乗Ｎ_１は、図６Ｂの式（７）で表すことができる。二階微分領域における、更新前の関数Ｆｘに対する更新後の関数Ｆｘのノルムの二乗Ｎ_２についても二階微分の演算子Ｈを用いて同様の式で表すことができる（図示せず）。従って、ペナルティ項を加えてカルマンゲインＫを導出する式は、図６Ｂの式（８）のようになる。式（８）のκ、λはそれぞれ、一階微分のペナルティに対する重み係数、二階微分のペナルティに対する重み係数であり、任意の実数を設定することができる。例えば、κ、λに０を設定すれば、局所フィッティングのときのカルマンゲインＫを導出する式となる。例えば、λに０を設定し、一階微分のペナルティだけを考慮すると、更新後の関数Ｆｘは、更新前の関数Ｆｘの全体をある方向へシフトさせたような形状となる。反対に、κに０を設定し、二階微分のペナルティだけを考慮すると、更新後の関数Ｆｘは、更新前の関数Ｆｘの全体形状（曲率）を変更するように修正される。κとλを適切に調整すると、一階微分と二階微分の両方を用いると、局所と大域の両面で精度の良い推定結果が得られる。

（動作）
［関数Ｆｘの更新処理］
次に図７を参照して、カーブフィッティング、局所フィッティングの流れについて説明する。これらの処理は、一般的にカルマンフィルタを適用して状態推定を行う処理と同様である。観測データは、Ｍｖ値とＣｖ値の組合せのデータであり、推定する状態量はＭｖ値に対するＣｖ値の真値である。データ取得部１１は、Ｍｖ値とＣｖ値を組合せた観測データを例えば所定の時間間隔で取得しているとする。

図７は、実施形態に係るフィッティング処理の一例を示すフローチャートである。
まず、入力受付部１２が、各種設定を受け付ける（ステップＳ１）。例えば、入力受付部１２は、式（８）のκ、λ、後述する各種閾値を取得し、記憶部１５に登録する。
次に、計算部１３が、初期化処理を行う（ステップＳ２）。初期化処理では、計算部１３が、関数Ｆｘの初期化、誤差共分散行列の初期化、状態ノイズと観測ノイズの初期化などを行う。関数Ｆｘの初期化は、これまでにデータ取得部１１が取得した観測データに基づいて、計算部１３が、区間線形モデルを作成する。例えば、図２の点Ｐ１～Ｐ９が得られていれば、各点を隣接する点を結んで関数Ｆｘを設定する（関数Ｆｘの初期化）。
次に計算部１３が、滑らかさに関するペナルティ項を追加した式を作成する（ステップＳ３）。計算部１３は、ステップＳ１で設定されたκ、λを用いて図６Ｂの式（８）を作成する。計算部１３は、局所フィッティングを行うときにはκ、λに０を設定し、カーブフィッティングを行うときにはκ、λに設定された値をセットする。

次に計算部１３が、状態方程式と観測方程式を作成し、予測処理を行う（ステップＳ４）。
状態方程式：ｘ_ｔ＝Ａｘ_ｔ－１＋Ｂｕ_ｔ＋状態遷移分のノイズ
ここで、ｘ_ｔは時刻ｔのｘ、ｘは真値を推定したい状態変数（Ｃｖ値とＭｖ値の組合せ）、ｕ_ｔは入力変数、Ａは状態遷移行列、Ｂは駆動行列である。
観測方程式：ｙ_ｔ＝Ｃｘ_ｔ＋観測ノイズ
ここで、ｙは観測される変数（Ｃｖ値）、Ｃは観測行列である。
例えば、計算部１３は、状態方程式に基づいて、誤差共分散行列Ｐに状態遷移分のノイズ（状態遷移誤差）を付加する（Ｐ_ｔ＝Ｐ_ｔ－１＋状態遷移誤差）。次に、計算部１３は、現在のＣｖ値を取得して観測Ｃｖ値として設定し、現在のＭｖ値を取得して観測行列Ｃを算出する。この観測行列は、上記の観測方程式における係数Ｃである。計算部１３は、観測行列Ｃと観測方程式に基づいてＣｖ値を予測する。

次に計算部１３が、カルマンゲインと誤差共分散行列を算出する（ステップＳ５）。例えば、計算部１３は、観測されたＣｖ値と予測されたＣｖ値の観測誤差ｅｔを計算する（ｅ_ｔ＝ｙ_ｔ－Ｃ_ｔｘ_ｔ）。また、計算部１３は、更新前のｘ_ｔと、ｘ_ｔの推定値のＣｖ値方向の差を計算する。次に計算部１３は、観測誤差の共分散行列（式（１）の右辺２項目のカッコ内）を計算する。次に計算部１３は、式（８）に基づいてカルマンゲインＫを算出する。

次に計算部１３は、関数Ｆｘを更新する（ステップＳ６）。例えば、計算部１３は、カルマンゲインＫと観測誤差ｅ_ｔを用いて、ｘ_ｔ＋１の推定値を更新する（ｘ_ｔ＋１＝ｘ_ｔ＋Ｋｅ_ｔ）。次に計算部１３は、次の時刻ｔ＋１の誤差共分散行列Ｐｔ＋１を計算する。

ここで、Ｆｘ推定値が制限を超える場合には、補正を行ってもよい。例えば、許容する更新幅に制限を設けておき、更新前の関数Ｆｘの値ｘ_ｔ＋１とｘ_ｔ＋１推定値の差が所定の閾値以上離れている場合、ｘ_ｔ＋１の推定値を、最大で変化したときの値に置き換える。例えば、（Ｃｖ値、Ｍｖ値）とすると、（１，１）から（１．１、１）までの変化を許容すると定めている場合に、ｘ_ｔ＋１推定値が（１．５、１）となるような場合、計算部１３は、推定された（１．５、１）を（１．１、１）に補正する。

また、区分線形モデルの状態誤差分散Ｐの更新時に、各点の状態誤差分散に上下限を設定してもよい。例えば、観測データがあまり得られない範囲（例えば、図３、図４の点Ｐ７）の点に対して許容する状態誤差分散の上下限値は狭く設定する。これにより、状態遷移誤差が足し込まれ蓄積されることを防止することができる。

また、点ごとに更新速度を変更しても構わない。例えば、観測データがあまり得られない範囲については、状態誤差分散の上下限値を狭く設定し、状態遷移誤差を大きく設定することで真値の推定を早く（更新速度を速く）することができる。反対に観測データが頻繁に得られる範囲（例えば、図３、図４の点Ｐ３～Ｐ５）では、状態遷移誤差を小さく設定することで、更新速度を遅くすることができる。また、このとき観測データが頻繁に得られる範囲では、状態誤差分散の上下限値を広く設定することで、（更新に時間がかかるが）大幅な変更を許容することになる。例えば、あまり観測データが得られない範囲（あまり使用しない弁開度）については、状態誤差分散に上下限を狭く設定し、状態遷移誤差を大きく設定することで、観測データが得られたときにすぐに更新するがあまり大きな変化は許容しないようにし、よく使われる範囲（観測データが頻繁に得られる範囲）では、その逆の設定を行うことで、少しずつであるが長期的には大きな変化を許容するといった設定にすることができる。よく使われる範囲では、経年劣化などにより長期的には大きく特性が変化することがあり得るので、このような設定は効果的である。

また、カルマンゲインに対して、点ごとに異なる重みを付けて更新量を更新してもよい。例えば、あまり大幅に更新したくない節点については小さな重みを付すことができる。また、一次微分、二次微分のペナルティ項に対し、点ごとに重み係数κ、λの値を変えてもよい。

関数Ｆｘの真値（Ｍｖ値とＣｖ値の真値、図３、図４の点Ｐ４´等の更新後の点）を更新すると、計算部１３は、更新後の隣り合う点を結んで関数Ｆｘを更新する。
計算部１３は、関数Ｆｘの更新を継続する間は、上記のステップＳ３～Ｓ６の処理を繰り返し行う。

［全体処理］
次に図８を参照して、局所フィッティング、カーブフィッティングの切り替えを含む関数Ｆｘの作成処理全体の一例について説明する。
図８は、実施形態に係る区分線形モデル作成処理の一例を示すフローチャートである。
一例として、「調整モード」、「全体Ｆｉｔｔｉｎｇ」、「局所Ｆｉｔｔｉｎｇ」の３つの処理モードを用意する。まず「調整モード」において関数Ｆｘの初期値を設定し、その後、「全体Ｆｉｔｔｉｎｇ」の処理モードにおいては、カーブフィッティングを実行し、関数Ｆｘの全体の形状を時々刻々と取得される観測データに応じて調整する。そして、必要に応じて、「局所Ｆｉｔｔｉｎｇ」に切り替えて関数Ｆｘの精度を維持する。

（調整モード）
まず、入力受付部１２は各種設定を受け付け、計算部１３はカルマンフィルタの初期化等を行う（ステップＳ１０）。この処理は、図７のステップＳ１、２と同様であるが、図８のフローチャートでは、関数Ｆｘの初期化をステップＳ１０～Ｓ１３の調整モードで行う。次にデータ取得部１１が、観測データを取得する（ステップＳ１１）。例えば、データ取得部１１は、バルブの動作範囲の全域にわたる観測データを取得する。データ取得部１１は、観測データを記憶部１５に記録する。次に、計算部１３が、記憶部１５から観測データを読み出して、カーブフィッティングを行う（ステップＳ１２）。カーブフィッティングの処理については、図７のフローチャートで説明したとおりである。出力部１４は、カーブフィッティングの結果（区分線形モデル、関数Ｆｘ）と観測データとを重畳して表示装置などに出力する。ユーザは、例えば、関数Ｆｘの形状と観測データの分布を見比べる等して、カーブフィッティングの精度の確認等を行う。カーブフィッティングが完了し、関数Ｆｘの初期値が設定できると、ユーザはその旨を区分線形モデル作成装置１０へ入力する。関数Ｆｘの設定が不十分（精度が良好ではないなど）な場合（ステップＳ１３；Ｎｏ）、ユーザは各種閾値や係数などを調整する。また、新たな観測データが得られる場合は観測データを追加して、計算部１３は、関数Ｆｘの初期値が設定できるまで、カーブフィッティングを繰り返し行う。このように初期段階でカーブフィッティングを行うことで、精度の高い機器特性の区分線形モデルを得ることができる。例えば、機器が使用される環境などの影響で、メーカが提供する機器特性と実際にプラント等に導入した機器の特性に差があるような場合、調整モードの実行によって、関数Ｆｘの初期値を設定することにより、速やかに高精度な区分線形モデルを得ることができる。

関数Ｆｘの設定が完了すると（ステップＳ１３；Ｙｅｓ）、「全体Ｆｉｔｔｉｎｇ」の処理モードに移行する。その前に、計算部１３は、運転停止するか否かを判定する（ステップＳ１４）。運転停止する場合（ステップＳ１４；Ｙｅｓ）、関数Ｆｘの更新を終了する。運転停止しない場合（ステップＳ１４；Ｎｏ）、次の全体フィッティングの処理モードへ移行する。なお、ここでは一例として「調整モード」から「全体Ｆｉｔｔｉｎｇ」へ移行する際に終了判定を行うこととしたが、例えば、ステップＳ１１～Ｓ１３のループ処理の実行中に運転停止命令が入力された場合、関数Ｆｘの更新を終了してもよい。

（全体Ｆｉｔｔｉｎｇ）
「全体Ｆｉｔｔｉｎｇ」では、調整モードで設定した関数Ｆｘを、実際の機器特性を精度よく反映したオリジナルな関数Ｆｘとして設定する（ステップＳ１５）。オリジナルな関数Ｆｘが設定できると、この関数Ｆｘをプラントの監視や制御に利用することができる。オリジナル関数Ｆｘの設定後も、データ取得部１１は、観測データを取得する（ステップＳ１６）。計算部１３は、設定された更新速度、閾値、重み係数に従って、全ての点を対象にカーブフィッティングを行う（ステップＳ１７）。計算部１３は、カーブフィッティングによって更新される前の値と更新後の値を各点について計算し、計算結果を閾値と比較する。この閾値には、各点ごとに異なる値が設定されていてもよい。計算部１３は、どの点についても閾値を上回る変化がない場合（ステップＳ１８；Ｎｏ）、取得された観測データに基づくカーブフィッティングを継続して行う。出力部１４は、関数Ｆｘが更新される度に、更新後の関数Ｆｘを表示装置などに出力する。

計算部１３は、どれか１つの点でも閾値を上回る変化があった場合（ステップＳ１８；Ｙｅｓ）、「局所Ｆｉｔｔｉｎｇ」の処理モードに移行する。その前に、計算部１３は、運転停止するか否かを判定する（ステップＳ１９）。運転停止する場合（ステップＳ１９；Ｙｅｓ）、関数Ｆｘの更新を終了する。運転停止しない場合（ステップＳ１９；Ｎｏ）、次の「局所Ｆｉｔｔｉｎｇ」の処理モードへ移行する。なお、ここでは一例として「全体Ｆｉｔｔｉｎｇ」から「局所Ｆｉｔｔｉｎｇ」へ移行する際に終了判定を行うこととしたが、ステップＳ１６～Ｓ１８のループ処理の実行中に運転停止命令が入力された場合、即座に関数Ｆｘの更新を終了してもよい。

（局所Ｆｉｔｔｉｎｇ）
「局所Ｆｉｔｔｉｎｇ」の処理モードでは、フィッティング方式をカーブフィッティングから局所フィッティングに切り替える。例えば、バルブの一部が損傷したような場合、ある開度のときだけＣｖ値が大きく変化してしまうことがある。このような場合に、ペナルティ項を設定して、損傷の影響による特定の範囲のＭｖ値に対するＣｖ値の大幅な変化を区分線形モデル全体に及ぼすようにすると、他の範囲については、精度を低下させることになる。従って、例えば、ある特定の範囲において更新前後で閾値を上回る変化が現れたときには、カーブフィッティングよりも局所フィッティングを適用した方が全体として精度を確保できると考え、フィッティング方式をカーブフィッティングから局所フィティングに切り替えて、関数Ｆｘの形状が、図３の破線で例示するような形状となることを許容する。局所フィティングに切り替えた後は、データ取得部１１は、観測データを取得する（ステップＳ２０）。計算部１３は、各点それぞれに対し、他の点とは独立して局所フィッティングを行う（ステップＳ２１）。計算部１３は、運転停止するまでの間（ステップＳ２２；Ｎｏ）、局所フィッティングを継続して実行する。運転停止すると（ステップＳ２２；Ｙｅｓ）、関数Ｆｘの更新を終了する。なお、ステップＳ１７のフィッティング方式の切り替え判定は、自動で行ってもよいし、ユーザが行ってもよい。例えば、バルブの一部損傷が検知された場合、ユーザは、区分線形モデル作成装置１０へ、局所フィッティングへの切り替えを指示する入力を行い、この入力に基づいて、計算部１３は、カーブフィッティングから局所フィッティングへの切り替えを行ってもよい。

（効果）
以上説明したように、本実施形態によれば、カルマンゲインＫを導出する式に滑らかさに関するペナルティ項を追加する。これにより、モデル化する対象の曲線との誤差を最小にしつつ、対象曲線を近似する滑らかな区分線形モデルを作成・更新することができる。これにより、非線形な特性を有する対象の値を精度よく推定することができる。また、本実施形態のフィッティング方法は、比較的少数の行列演算で実行することができるため、簡易に実装可能であり、高速に動作する。したがって、制御装置上でリアルタイムに演算するなどの用途に利用することができる。

図９は、実施形態の区分線形モデル作成装置のハードウェア構成の一例を示す図である。
コンピュータ９００は、ＣＰＵ９０１、主記憶装置９０２、補助記憶装置９０３、入出力インタフェース９０４、通信インタフェース９０５を備える。
上述の区分線形モデル作成装置１０は、コンピュータ９００に実装される。そして、上述した各機能は、プログラムの形式で補助記憶装置９０３に記憶されている。ＣＰＵ９０１は、プログラムを補助記憶装置９０３から読み出して主記憶装置９０２に展開し、当該プログラムに従って上記処理を実行する。また、ＣＰＵ９０１は、プログラムに従って、記憶領域を主記憶装置９０２に確保する。また、ＣＰＵ９０１は、プログラムに従って、処理中のデータを記憶する記憶領域を補助記憶装置９０３に確保する。

なお、区分線形モデル作成装置１０の全部または一部の機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより各機能部による処理を行ってもよい。ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、ＣＤ、ＤＶＤ、ＵＳＢ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。また、このプログラムが通信回線によってコンピュータ９００に配信される場合、配信を受けたコンピュータ９００が当該プログラムを主記憶装置９０２に展開し、上記処理を実行しても良い。また、上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであってもよい。なお、区分線形モデル作成装置１０は、複数のコンピュータ９００によって構成されていても良い。

以上のとおり、本開示に係るいくつかの実施形態を説明したが、これら全ての実施形態は、例として提示したものであり、発明の範囲を限定することを意図していない。これらの実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で種々の省略、置き換え、変更を行うことができる。これらの実施形態及びその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれる。

＜付記＞
各実施形態に記載の区分線形モデル作成装置、区分線形モデル作成方法およびプログラムは、例えば以下のように把握される。

（１）第１の態様に係る区分線形モデル作成装置１０は、区分線形モデルにおける区分を規定する節点の真値をカルマンフィルタによって推定するにあたり、前記区分線形モデルの滑らかさに関するペナルティ項をカルマンゲインの算出式に追加して、カルマンゲインおよび誤差共分散行列を算出して、前記節点の真値を推定し、推定された前記節点を結んで区分線形モデルを作成する計算部１３、を備える。
これにより、区分線形モデル全体の形状の滑らかさを維持したまま、区分線形モデルを作成、更新することができる。

（２）第２の態様に係る区分線形モデル作成装置１０は、（１）の区分線形モデル作成装置１０であって、前記ペナルティ項は、前記区分線形モデルの更新前と更新後の差の一階微分のノルムおよび／または二階微分のノルムを含む。
これにより、一階微分のノルムおよび／または二階微分のノルムを含めることで、現在の滑らかな形状を維持したままスライドさせたり（一階微分）、観測データの分布に沿った滑らかな形状に修正したり（二階微分）、それらをバランスさせて滑らかさを保ったまま精度の良い区分線形モデルを作成することができる。

（３）第３の態様に係る区分線形モデル作成装置１０は、（１）～（２）の区分線形モデル作成装置１０であって、前記計算部は、Ｐを更新後の誤差共分散行列、ｔｒをトレース、Ｎ_１を前記区分線形モデルの更新前と更新後の差の一階微分のノルム、Ｎ_２を前記区分線形モデルの更新前と更新後の差の二階微分のノルム、Ｋをカルマンゲイン、κとλを重み係数としたときに以下の式を満たすＫを算出する。
∂ｔｒ（Ｐ）／∂Ｋ＋κ・∂Ｎ_１／∂Ｋ＋λ・∂Ｎ_２／∂Ｋ＝０
これにより、一般的なカルマンフィルタと同様の手法でデータ同化を行うことができる。

（４）第４の態様に係る区分線形モデル作成装置１０は、（１）～（３）の区分線形モデル作成装置１０であって、前記計算部は、前記ペナルティ項を除外してカルマンゲインＫを算出する局所フィッティングと、前記ペナルティ項を追加してカルマンゲインＫを算出するカーブフィッティングと、を切り替えて実行することができる。
対象の特性に応じて局所フィッティングとカーブフィッティングを使い分けることで、区間線形モデルの精度（特性の記述力）を確保することができる。

（５）第５の態様に係る区分線形モデル作成装置１０は、（４）の区分線形モデル作成装置１０であって、前記計算部は、前記カーブフィッティングの実行中に、更新前の前記節点と更新後の前記節点の差が所定の閾値を上回る場合、前記カーブフィッティングから前記局所フィッティングに切り替える。
これにより、対象の特性の変化に速やかに対応することができる。

（６）第６の態様に係る区分線形モデル作成装置１０は、（１）～（５）の区分線形モデル作成装置１０であって、前記計算部は、更新前の前記節点と更新後の前記節点の差が所定の閾値以下となるよう、更新後の前記節点の値を上限値または下限値で上書きする。
これにより、区間線形モデルの急激かつ大幅な更新を防止することができる。例えば、外乱等の影響により一時的に観測データが変動したような場合でも、その影響を一定の範囲に抑えることができる。

（７）第７の態様に係る区分線形モデル作成装置１０は、（１）～（６）の区分線形モデル作成装置１０であって、前記節点ごとに状態遷移誤差、更新速度、前記カルマンゲインへの重み付け、前記ペナルティ項への重み付けを設定して、前記節点の真値を推定する。
これにより、各節点の付近で得られる観測データの特性に応じて真意を推定し、区間線形モデルの推定精度を向上することができる。

（８）第８の態様に係る区分線形モデル作成方法は、区分線形モデルにおける区分を規定する節点の真値をカルマンフィルタによって推定するにあたり、前記区分線形モデルの滑らかさに関するペナルティ項をカルマンゲインの算出式に追加して、カルマンゲインおよび誤差共分散行列を算出して前記節点の真値を推定し、推定された前記節点を結んで区分線形モデルを作成する。

（９）第９の態様に係るプログラムは、コンピュータ９００に、区分線形モデルにおける区分を規定する節点の真値をカルマンフィルタによって推定するにあたり、前記区分線形モデルの滑らかさに関するペナルティ項をカルマンゲインの算出式に追加して、カルマンゲインおよび誤差共分散行列を算出して前記節点の真値を推定し、推定された前記節点を結んで区分線形モデルを作成する処理を実行させる。

１０・・・区分線形モデル作成装置
１１・・・データ取得部
１２・・・入力受付部
１３・・・計算部
１４・・・出力部
１５・・・記憶部
９００・・・コンピュータ
９０１・・・ＣＰＵ
９０２・・・主記憶装置
９０３・・・補助記憶装置
９０４・・・入出力インタフェース
９０５・・・通信インタフェース

Claims (9)

1. 区分線形モデルにおける区分を規定する節点の真値をカルマンフィルタによって推定するにあたり、前記区分線形モデルの滑らかさに関するペナルティ項をカルマンゲインの算出式に追加して、カルマンゲインおよび誤差共分散行列を算出して、前記節点の真値を推定し、推定された前記節点を結んで区分線形モデルを作成する計算部、
   を備える区分線形モデル作成装置。
2. 前記ペナルティ項は、前記区分線形モデルの更新前と更新後の差の一階微分のノルムおよび／または二階微分のノルムを含む、
   請求項１に記載の区分線形モデル作成装置。
3. 前記計算部は、Ｐを更新後の誤差共分散行列、ｔｒをトレース、Ｎ_１を前記区分線形モデルの更新前と更新後の差の一階微分のノルム、Ｎ_２を前記区分線形モデルの更新前と更新後の差の二階微分のノルム、Ｋをカルマンゲイン、κとλを重み係数としたときに以下の式を満たすＫを算出する、
   ∂ｔｒ（Ｐ）／∂Ｋ＋κ・∂Ｎ_１／∂Ｋ＋λ・∂Ｎ_２／∂Ｋ＝０
   請求項１または請求項２に記載の区分線形モデル作成装置。
4. 前記計算部は、前記ペナルティ項を除外してカルマンゲインＫを算出する局所フィッティングと、前記ペナルティ項を追加してカルマンゲインＫを算出するカーブフィッティングと、を切り替えて実行することができる、
   請求項１から請求項３の何れか１項に記載の区分線形モデル作成装置。
5. 前記計算部は、前記カーブフィッティングの実行中に、更新前の前記節点と更新後の前記節点の差が所定の閾値を上回る場合、前記カーブフィッティングから前記局所フィッティングに切り替える、
   請求項４に記載の区分線形モデル作成装置。
6. 前記計算部は、更新前の前記節点と更新後の前記節点の差が所定の閾値以下となるよう、更新後の前記節点の値を上限値または下限値で上書きする、
   請求項１から請求項５の何れか１項に記載の区分線形モデル作成装置。
7. 前記計算部は、前記節点ごとに状態遷移誤差、更新速度、前記カルマンゲインへの重み付け、前記ペナルティ項への重み付けを設定して、前記節点の真値を推定する、
   請求項１から請求項６の何れか１項に記載の区分線形モデル作成装置。
8. 区分線形モデルにおける区分を規定する節点の真値をカルマンフィルタによって推定するにあたり、前記区分線形モデルの滑らかさに関するペナルティ項をカルマンゲインの算出式に追加して、カルマンゲインおよび誤差共分散行列を算出して前記節点の真値を推定し、推定された前記節点を結んで区分線形モデルを作成する、
   区分線形モデル作成方法。
9. コンピュータに、
   区分線形モデルにおける区分を規定する節点の真値をカルマンフィルタによって推定するにあたり、前記区分線形モデルの滑らかさに関するペナルティ項をカルマンゲインの算出式に追加して、カルマンゲインおよび誤差共分散行列を算出して前記節点の真値を推定し、推定された前記節点を結んで区分線形モデルを作成する処理、
   を実行させるプログラム。

Images