A Method to Estimate Initial States, Inputs and Parameters for Diagnosis of Equipment

Abstract

The equipment and the structures in steel works are used in a very long term. Therefore, the inspections and the maintenance are indispensable. To make this more effective, much effort has been devoted to develop methods of fault detection and diagnosis. However, the effectiveness of those methods may be limited due to the tradeoff between false positive and false negative reactions. To mitigate the tradeoff, this paper considers to use model sets involving parameters and disturbances that are expected to change. More specifically, we first propose a method to estimate the parameters and the disturbance so as to minimize the deviation between the real output and the output generated by the model sets. The resulting residual enables us to conduct the health diagnosis. The effectiveness of the proposed framework is examined by using numerical examples.

1. はじめに

製鉄所には様々な設備や建造物がある。それらは長期間に渡って運用されるため，それらの保守・点検が不可欠である。通常，設備の保守・点検は定期的あるいは問題が生じた際に行われる。しかしながら，異常は問題が周辺に波及する前に，早期に発見されるほうが望ましい。点検を頻繁に行えば早期発見は可能になるが，保守費用の増加や操業効率の低下を招く。そのため，何らかのセンサーを用いた常時監視によって異常検出の早期化や異常診断の自動化を可能にするための研究が従来から行われている^1,2)。また，最近では構造物のヘルスモニタリングなどの研究も盛んになってきている³⁾。モニタリング手法の研究に加えて，高速・高精度な画像計測^7,8)も可能になってきており，異常診断に必要な要素技術のレベルは高まってきている。

従来から知られている異常検出・診断は，カルマンフィルターの出力における事前推定値と観測値との間の差の急激な変化などを検出の拠り所としている。このような手法において変化に敏感に反応するように異常検出システムを設計すると，設備の経年劣化などの事前に想定される変化に対して頻繁にアラームを発生させてしまう。一方，変化に対して鈍感になるように設計すると，本当に検出したい異常を見逃すリスクを高めてしまう。このような状況に対し，変化が想定される要因を組み込んだモデルを用いることで，得られているデータが想定される変化に起因するものか否かを判別できれば，このようなトレードオフの緩和が可能になると考えられる。そこで本論文では，想定される変化を含むモデルを用いることで，想定される変化の推定を行うと同時に，モデルとデータの乖離に基づいて想定外の変化の存在を検出する手法を考える。

設備で想定される変化にはパラメータの変化や操業に起因する外乱によるものなどが挙げられる。このうち，パラメータの推定には拡張カルマンフィルターによる方法⁴⁾や非線形最適化に基づく手法などがある^5,6)。パラメータ推定時における外乱の存在に対しては，外乱の影響を取り除く方法が提案されている⁹⁾。一方，文献¹⁰⁾では多点同時計測を行うことで，外乱そのものの情報を用いることなくデータだけから異常の検出を行っている。このような結果を踏まえると，多点同時計測とモデルを組み合わせれば，パラメータだけでなく外乱の推定も可能になることが期待できる。

以上より，本論文では文献⁶⁾の手法を拡張することにより，パラメータ・外乱に依存する対象モデルと多チャンネルデータを用いることで，初期状態・パラメータ・入力を同時に推定する手法を提案する。診断する設備は一般に動的システムであるので，任意の時間区間のデータを用いて診断する際には初期状態の推定も必要となる。また，提案する推定手法を設備診断に適用する方法も提案し，数値例を用いてその有効性を検証する。

本論文の構成は以下の通りである。まず，2章で本論文で提案する設備診断の枠組を提案し，その枠組を実現するために初期状態・入力・パラメータを推定する手法が必要となることを示す。つぎに，3章において基本的な推定手法を提案し，その有効性を数値例を用いて検証する。さらに，4章で出力データにノイズが存在する場合を想定し，フィルターの導入と入力の周波数帯域制約を行うための手法を提案する。また，その有効性を数値例を用いて検証する。

本論文の記法は標準的なものである。Rは実数の集合をあらわす。また，ベクトル列u={u[0],...,u[N-1]}に対し，vec(u)は

  

$$\text{vec}\left(u\right):=\left[\begin{array}{c}u\left[0\right]\\ ⋮\\ u\left[N-1\right]\end{array}\right]$$

をあわらす。

2. モデル集合に基づく設備診断・異常検出

設備の異常検出・診断を行う際，カルマンフィルターの出力における事前推定値と観測値との間の差の急激な変化などを検出の拠り所とすることが多い。しかしながら，このような方法で診断システムの感度を高めると，事前に想定される変化に対して頻繁にアラームを発生させてしまう。一方，感度を下げると，本当に検出したい異常を見逃すリスクを高めてしまう。このように，設備の診断を変化だけに依拠すると偽陽性(false positive)と偽陰性(false negative)の割合を変えることは可能であっても，両者の間のトレードオフを回避することができない。

ここで，事前に想定される変化をモデル化し，そのモデルで再現できるものを正常と判断すれば，偽陽性アラームの発生を減少させることができる。また，それによって診断システムの感度を高めることも可能になるので，偽陰的な振舞いを減らすことも同時に可能となる。すなわち，事前に想定される変化のモデルを用いて事前に想定される変化とそれ以外の変化を分離して識別することができれば，偽陽性・偽陰性間のトレードオフの緩和が可能になると考えられる。そこで本論文ではこのような識別を可能にする具体的に手法を提案する。

Fig.1は以上の考え方に基づいた設備診断の枠組を図示したものである。まず，対象となる設備のモデルを構築する。ただし，ここでのモデルは事前に想定される変化も含むモデル集合である。このとき，定められたモデル集合から生成可能な出力の集合Y_mを考えることができる。もし実際の設備から得られる出力がY_mに属していれば設備は正常と判断する。この場合，設備から得られるデータと一致するモデルを特定することができ，その結果を見ることで必要な対策を検討することができる。一方，Y_mに属していなければ，事前に想定していない変化が生じていると判断することができる。もしこのような場合が生じれば，その時点で設備の保守・点検を行えばよい。このように，Fig.1の枠組を実現することができれば，より効率的かつ精密に設備の異常を診断することが可能となる。

Fig. 1.

Framework of the diagnosis using model sets. (Online version in color.)

Fig.1の枠組を実現するためには，設備から得られた出力がY_mに属しているか否かを判定し，Y_mに属していない場合にY_mと得られた出力との乖離の程度を評価する必要がある。そのためには，モデル集合から得られる出力のうち，設備から得られたデータと最も近くなるような出力を求めればよい。これはモデル集合の中から出力誤差を最小化するモデルを推定することと等価である。一方，モデル集合は初期状態，外乱入力，パラメータなどによって記述される。そこで以降では，まず出力誤差を最小化する初期状態・入力・パラメータを推定する手法を提案する。その上で，それを設備診断に適用した場合の有効性を数値例を用いて検証する。

3. 初期状態・パラメータ・入力の同時推定

3・1　推定手法

対象システムに対して，何らかの入力u[k]∈R^m，r[k]∈R^ν(k=0,...,N-1)が印加され，その出力データy_r[k]∈R^p(k=0,...,N-1)が得られたものとする。ただし，rは目標値や制御入力などの既知の入力であるのに対し，uは外乱などの未知の入力である。

出力データに加えて，推定したいパラメータに依存する対象モデルが与えられたものとする。とくにここでは，モデルとして次の離散時間状態方程式および出力方程式が

  

$$x\left[k+1\right]=f\left(x\left[k\right],q,u\left[k\right],r\left[k\right]\right)$$(1)

  

$$y\left[k\right]=g\left(x\left[k\right],q,u\left[k\right],r\left[k\right]\right)$$(2)

で与えられているものとする。ただし，x[k]∈Rⁿ，y[k]∈R^p，q∈R^aはそれぞれ状態変数，出力，パラメータである。また，関数

  

$$f\left(\overline{x},\overline{q},\overline{u},\overline{r}\right),g\left(\overline{x},\overline{q},\overline{u},\overline{r}\right)$$(3)

はいずれもx，q，uに対して連続微分可能であるとする。式(1)，(2)のモデルを用いることにより，与えられた初期状態x[0]=x_e，パラメータq=q_e，入力データu[k]，r[k](k=0,...,N-1)に対して，出力信号y[k](k=0,...,N-1)を求めることができる。

ここでは，対象システムの出力データ

  

$$y_r=\left\{y_r\left[0\right],\dots ,y_r\left[N-1\right]\right\}$$(4)

とモデル出力

  

$$y=\left\{y\left[0\left],\dots ,y\right[N-1\right]\right\}$$(5)

との間の誤差が小さくなるようなモデルの初期状態x[0]=x_e，パラメータq=q_e，入力u={u[0],...,u[N-1]}を推定することを考える。ただし，モデル出力を求める際のr[k](k=0,...,N-1)は既知とする。上記の枠組で初期状態・パラメータ・入力を推定するためには，少くとも得られる変数の数が推定したい変数の数よりも多くなければならない。よって，pN≥n+a+mNであることが必要である。この必要条件として，さらにp>mであることが導かれる。すなわち，推定すべき入力の数よりも多くの独立な出力をとる必要がある。また，上記の枠組みではパラメータや未知入力だけでなく初期状態も推定するため，使用するデータは計測した出力データのどの部分を切り取って使用しても構わない。データ計測開始時に特別な条件を設定する必要もなく，どのタイミングでデータ計測を開始しても構わない。

本論文では初期状態・パラメータ・入力の推定問題を最適化問題に帰着させる。具体的には，y_rとyの間の誤差

  

$$e=\left\{e\left[0\left],\dots ,e\right[N-1\right]\right\},$$(6)

  

$$e\left[k\right]=y\left[k\right]-y_r\left[k\right]\left(k=0,\dots ,N-1\right)$$(7)

の大きさを以下の評価関数J

  

$$J\left(x_e,q_e,u_e\right)=\frac{1}{2}{\displaystyle \sum }_{k=0}^{N-1}{\Vert e\left[k\right]\Vert }^2$$(8)

で評価する。ただし， ‖・‖ はユークリッドノルムである。u_e∈R^mNは

  

$$u_e=\text{vec}\left(u\right)$$(9)

である。逆にu_eが与えられたとき，対応する入力は

  

$$u\left[k\right]=U_k{u}_e,$$

  

$$U_k=\left[\begin{array}{ccc}{0}_{m\times \left(mk\right)}& I_m& 0_{m\times m\left(N-1-k\right)}\end{array}\right]$$

で与えられる。vec(・)を用いると，式(8)の評価関数を以下のように書くことができる。

  

$$J\left(x_e,q_e,u_e\right)=\frac{1}{2}{\Vert \text{vec}\left(e\right)\Vert }^2$$(10)

本論文では式(10)のJの最小化問題の解

  

$$\left(x_e,q_e,u_e\right)=\arg \underset{x_e,q_e,u_e}{\min }J\left(x_e,q_e,u_e\right)$$(11)

を初期状態・パラメータ・入力の推定値とする。

本論文ではガウス・ニュートン法¹¹⁾に基づいて式(10)の最小化を行う。すなわち，第lステップの解候補

  

$$\left[\begin{array}{c}{x}_e^{\left(l\right)}\\ q_e^{\left(l\right)}\\ u_e^{\left(l\right)}\end{array}\right]$$

に対し，第l+1ステップの解候補を次式で定める。

  

$$\left[\begin{array}{c}{x}_e^{\left(l+1\right)}\\ q_e^{\left(l+1\right)}\\ u_e^{\left(l+1\right)}\end{array}\right]=\left[\begin{array}{c}{x}_e^{\left(l\right)}\\ q_e^{\left(l\right)}\\ u_e^{\left(l\right)}\end{array}\right]-{\left(M^{\left(l\right)}\right)}^{†}\text{vec}\left(e\right)$$(12)

ただし，M^(l)はvec(e)の最適化変数に対するヤコビアン

  

$$\begin{array}{l}{M}^{\left(l\right)}={\left[\begin{array}{ccc}\frac{\partial vec\left(e\right)}{\partial x\left[0\right]}& \frac{\partial vec\left(e\right)}{\partial q}& \frac{\partial vec\left(e\right)}{\partial u_e}\end{array}\right]|}_{\begin{array}{c}x\left[0\right]=x_e^{\left(l\right)},\\ q=q_e^{\left(l\right)},\\ u_e=u_e^{\left(l\right)}\end{array}}\\ ={\left[\begin{array}{ccc}\frac{\partial y\left[0\right]}{\partial x\left[0\right]}& \frac{\partial y\left[0\right]}{\partial q}& \frac{\partial y\left[0\right]}{\partial u_e}\\ ⋮& ⋮& ⋮\\ \frac{\partial y\left[N-1\right]}{\partial x\left[0\right]}& \frac{\partial y\left[N-1\right]}{\partial q}& \frac{\partial y\left[N-1\right]}{\partial u_e}\end{array}\right]|}_{\begin{array}{c}x\left[0\right]=x_e^{\left(l\right)},\\ q=q_e^{\left(l\right)},\\ u_e=u_e^{\left(l\right)}\end{array}}\end{array}$$

であり，(M^(l))^†はM^(l)の擬似逆行列¹²⁾である。ここで，y_rは最適化変数とは無関係なので，eに対するヤコビアンはyに対するヤコビアンに一致することに注意されたい。

M^(l)は出力系列の初期状態，パラメータ，入力系列に対するヤコビアンであるので，M^(l)は式(1)，(2)のシステムの動特性を踏まえて求める必要がある。実際，M^(l)は以下のシステム

  

$$\xi \left[k+1\right]=A\left[k\right]\xi \left[k\right]+E\left[k\right]\left[\begin{array}{ccc}{0}_{a\times n}& I_a& 0_{a\times mN}\end{array}\right]+B\left[k\right]\left[\begin{array}{ccc}{0}_{m\times n}& 0_{m\times a}& U_k\end{array}\right]$$(13)

  

$$\eta \left[k\right]=C\left[k\right]\xi \left[k\right]+F\left[k\right]\left[\begin{array}{ccc}{0}_{a\times n}& I_a& 0_{a\times mN}\end{array}\right]+D\left[k\right]\left[\begin{array}{ccc}{0}_{m\times n}& 0_{m\times a}& U_k\end{array}\right]$$(14)

  

$$\xi \left[0\right]=\left[\begin{array}{ccc}{I}_n& 0_{n\times a}& 0_{n\times mN}\end{array}\right]$$(15)

の応答を用いて

  

$$M^{\left(l\right)}=\left[\begin{array}{c}\eta \left[0\right]\\ ⋮\\ \eta \left[N-1\right]\end{array}\right]$$(16)

で与えられる。ただし，A[k]，B[k]，C[k]，D[k]，E[k]，F[k]は式(3)の関数f，gのx，q，uに対する編微分係数

  

$$\begin{array}{l}A\left[k\right]=\frac{\partial f}{\partial \overline{x}},E\left[k\right]=\frac{\partial f}{\partial \overline{q}},B\left[k\right]=\frac{\partial f}{\partial \overline{u}},\hfill \\ C\left[k\right]=\frac{\partial g}{\partial \overline{x}},F\left[k\right]=\frac{\partial g}{\partial \overline{q}},D\left[k\right]=\frac{\partial g}{\partial \overline{u}}\hfill \end{array}$$(17)

である。式(17)の偏微分はいずれもx=x^(l)[k]，q=q_e^(l)，u=u_e^(l)[k]，r=r[k]で評価したものであり，x^(l)[k]はx[0]=x_e^(l)，q=q_e^(l)，u[i]=u_e^(l)[i](i=0,...,k-1)としたときの式(1)の応答である。よって，式(16)のM^(l)は，第ℓステップの解候補x_e^(l)，q_e^(l)，u_e^(l)を用いて式(1)，(13)，(14)，(15)の差分方程式の解を求めることによって得られる。

式(13)，(14)のξ[k]，η[k]はそれぞれ

  

$$\xi \left[k\right]=\left[\begin{array}{ccc}\frac{\partial x\left[k\right]}{\partial x\left[0\right]}& \frac{\partial x\left[k\right]}{\partial q}& \frac{\partial x\left[k\right]}{\partial u_e}\end{array}\right]$$(18)

  

$$\eta \left[k\right]=\left[\begin{array}{ccc}\frac{\partial y\left[k\right]}{\partial x\left[0\right]}& \frac{\partial y\left[k\right]}{\partial q}& \frac{\partial y\left[k\right]}{\partial u_e}\end{array}\right]$$(19)

をあらわしており，式(13)，(14)は式(18)，(19)を漸化式で表現したものである。実際，式(2)より

  

$$\begin{array}{l}\eta \left[k\right]=\left[\begin{array}{ccc}\frac{\partial y\left[k\right]}{\partial x\left[0\right]}& \frac{\partial y\left[k\right]}{\partial q}& \frac{\partial y\left[k\right]}{\partial u_e}\end{array}\right]\\ =\frac{\partial g}{\partial \overline{x}}\left[\begin{array}{ccc}\frac{\partial x\left[k\right]}{\partial x\left[0\right]}& \frac{\partial x\left[k\right]}{\partial q}& \frac{\partial x\left[k\right]}{\partial u_e}\end{array}\right]+\frac{\partial g}{\partial \overline{q}}\left[\begin{array}{ccc}\frac{\partial q}{\partial x\left[0\right]}& \frac{\partial q}{\partial q}& \frac{\partial q}{\partial u_e}\end{array}\right]+\frac{\partial g}{\partial \overline{u}}\left[\begin{array}{ccc}\frac{\partial u\left[k\right]}{\partial x\left[0\right]}& \frac{\partial u\left[k\right]}{\partial q}& \frac{\partial u\left[k\right]}{\partial u_e}\end{array}\right]\\ =\frac{\partial g}{\partial \overline{x}}\xi \left[k\right]+\frac{\partial g}{\partial \overline{x}}\left[\begin{array}{ccc}{0}_{a\times n}& I_a& 0_{a\times mN}\end{array}\right]+\frac{\partial g}{\partial \overline{u}}\left[\begin{array}{ccc}{0}_{m\times n}& 0_{m\times a}& U_k\end{array}\right]\end{array}$$(20)

であり，同様に，式(1)より

  

$$\begin{array}{l}\xi \left[k+1\right]=\left[\begin{array}{ccc}\frac{\partial x\left[k+1\right]}{\partial x\left[0\right]}& \frac{\partial x\left[k+1\right]}{\partial q}& \frac{\partial x\left[k+1\right]}{\partial u_e}\end{array}\right]\\ =\frac{\partial f}{\partial \overline{x}}\left[\begin{array}{ccc}\frac{\partial x\left[k\right]}{\partial x\left[0\right]}& \frac{\partial x\left[k\right]}{\partial q}& \frac{\partial x\left[k\right]}{\partial u_e}\end{array}\right]+\frac{\partial f}{\partial \overline{q}}\left[\begin{array}{ccc}\frac{\partial q}{\partial x\left[0\right]}& \frac{\partial q}{\partial q}& \frac{\partial q}{\partial u_e}\end{array}\right]+\frac{\partial f}{\partial \overline{u}}\left[\begin{array}{ccc}\frac{\partial u\left[k\right]}{\partial x\left[0\right]}& \frac{\partial u\left[k\right]}{\partial q}& \frac{\partial u\left[k\right]}{\partial u_e}\end{array}\right]\\ =\frac{\partial f}{\partial \overline{x}}\xi \left[k\right]+\frac{\partial f}{\partial \overline{x}}\left[\begin{array}{ccc}{0}_{a\times n}& I_a& 0_{a\times mN}\end{array}\right]+\frac{\partial f}{\partial \overline{u}}\left[\begin{array}{ccc}{0}_{m\times n}& 0_{m\times a}& U_k\end{array}\right]\end{array}$$(21)

である。

以上が提案する推定手法である。ただし，式(10)の評価関数は変数に関して凸とは限らないので，ガウスニュートン法によって大域的な最適解が得られるとは限らないことに注意されたい。

3・2　数値例に基づく有効性の検証

ここでは，数値例を用いて3・1節の推定手法の有効性を検証する。用いる数値例はFig.2の2慣性系である。k₁，k₂，ℓ₁₀，ℓ₂₀，c₁，c₂，m₁，m₂∈Rはそれぞれバネ定数，バネの自然長，ダンパ定数，質量である。質量m₁，m₂それぞれの壁面からの変位をp₁(t)，p₂(t)∈Rとする。このシステムの状態方程式は以下のように与えられる。

  

$${\dot{x}}_c\left(t\right)=f_c\left(x_c\left(t\right),q,u_c\left(t\right)\right)$$(22)

  

$$y\left(t\right)=C_c{x}_c\left(t\right)$$(23)

  

$$f_c\left(x_c\left(t\right),q,u_c\left(t\right)\right)=\left[\begin{array}{c}{\dot{p}}_1\left(t\right)\\ {\dot{p}}_2\left(t\right)\\ -\frac{k_1}{m_1}{\delta }_1\left(t\right)-\frac{c_1}{m_1}{\dot{\delta }}_1\left(t\right)+\frac{k_2}{m_1}{\delta }_2\left(t\right)+\frac{c_2}{m_1}{\dot{\delta }}_2\left(t\right)\\ -\frac{k_2}{m_2}{\delta }_2\left(t\right)-\frac{c_2}{m_2}{\dot{\delta }}_2\left(t\right)+\frac{u}{m_2}\end{array}\right]$$

  

$${\delta }_1\left(t\right)=p_1\left(t\right)-{\mathcal{\ell }}_{10}$$

  

$${\delta }_2\left(t\right)=p_2\left(t\right)-p_1\left(t\right)-{\mathcal{\ell }}_{20}$$

  

$$C_c=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right]$$

Fig. 2.

Spring-mass-damper system used as the numerical example. (Online version in color.)

ただし，u_c(t)∈Rは外乱，x_c(t)∈R⁴，y(t)∈R²，q∈R²はそれぞれ以下のように定義される状態変数，出力および被推定パラメータである。

  

$$x_c\left(t\right)=\left[\begin{array}{c}y\left(t\right)\\ \dot{y}\left(t\right)\end{array}\right],y\left(t\right)=\left[\begin{array}{c}{p}_1\left(t\right)\\ p_2\left(t\right)\end{array}\right],q=\left[\begin{array}{c}{k}_1\\ k_2\end{array}\right]$$(24)

Fig.2のパラメータのうち，ダンパ定数と質量は変化しない定数であり，その値は以下のものを仮定する。

  

$$c_1=c_2=0.1,m_1=m_2=1,$$

とする。一方，バネ定数は変化することが想定されるパラメータ，自然長は変化は想定していないが，実際には変化してしまうパラメータとする。

上記の対象に対して初期状態・パラメータ・入力の推定を行う。推定を行うために，まず式(22)，(23)に対して前進差分に基づく離散化を行うことで式(1)，(2)を構成した。具体的には

  

$$\left[\begin{array}{c}f\left(x\left[k\right],q,u\left[k\right],r\left[k\right]\right)\\ g\left(x\left[t\right],q,u\left[k\right],r\left[k\right]\right)\end{array}\right]=\left[\begin{array}{c}x\left[k\right]+\Delta t{f}_c\left(x\left[k\right],q,u\left[k\right]\right)\\ C_{c}x\left[k\right]\end{array}\right]$$(25)

である。ただし，Δt=0.04である。式(25)で与えられるシステムを用いて出力データを生成する。データを生成する際のパラメータは以下の通りである。

  

$$x\left[0\right]={\left[\begin{array}{cccc}1.1& 2.2& 0& 0\end{array}\right]}^T,q={\left[\begin{array}{cc}1& 1\end{array}\right]}^T,{\mathcal{\ell }}_{10}={\mathcal{\ell }}_{20}=1$$

外乱にはFig.3の実線の信号を入力する。得られた出力をFig.4に示す。データ数Nは512であり，信号の時間区間は0から20.44秒までである。

Fig. 3.

True input for the system and the initial input for optimization. (Online version in color.)

Fig. 4.

Output data obtained from the system. (Online version in color.)

得られたFig.4の出力に対し，3・1節の提案手法を用いて推定を行う。x[0]およびqの初期解には以下のものを与える。

  

$$x\left[0\right]={\left[\begin{array}{cccc}2& 3& 1& -1\end{array}\right]}^T,$$(26)

  

$$q={\left[\begin{array}{cc}0.7& 0.8\end{array}\right]}^T$$(27)

また，外乱の初期解にはFig.3の破線の値を用いる。推定結果として得られた初期状態，パラメータは

  

$$x\left[0\right]={\left[\begin{array}{cccc}1.100& 2.200& 0.000& 0.000\end{array}\right]}^T$$(28)

  

$$k_1=1.000,k_2=1.000$$(29)

である。一方，推定された入力をFig.5に示す。以上の結果から初期状態，パラメータ，入力を高い精度で推定できることが確認できる。ここで，出力データを生成する際のバネの自然長ℓ₀₁を0.9から1.1の間で変化させ，上記の推定を繰り返すことで得られた評価関数値Jの値をFig.6に示す。Fig.6より，評価関数値がk₁の変化に対しては変化しないのに対して，ℓ₀₁の対しては|ℓ₀₁-1|の値に応じて大きくなっている。このことから，評価関数値を観察することで想定外の変化が生じていることを検出することが可能である。

Fig. 5.

Estimated input given by the optimization. (Online version in color.)

Fig. 6.

Values of the objective function for each k₁ and l₀₁. (Online version in color.)

4. 出力データにノイズがある場合における同時推定

4・1　出力のフィルタリングと入力周波数制約

3章では初期状態・パラメータ・入力を同時に推定する手法を提案し，さらに数値例を用いて提案手法の有効性を示した。ただし，3章では出力データが正確に得られる場合を想定している。一方，現実に出力データを得る場合にはデータにノイズが含まれる場合がある。そのような場合にも3・1節の手法が有効であるとは限らない。そこでまず，ノイズが存在する場合における提案手法の振舞いを数値例に基づいて検証する。

ノイズが存在する出力データを作成するために，ℓ₀₁=ℓ₀₂=1として得られた出力データに平均0，分散r=0.01の独立正規分布の雑音を加算する。このとき，Fig.4に対応する出力データはFig.7である。このデータに対して3・1節の手法を適用して推定を行い，Fig.6に対応する計算を行ったところ，Fig.8の結果が得られた。Fig.8の2つの曲面のうち，上側の曲面がここで得られた結果であり，下側の曲面は比較のためにFig.6を同じものを示したものである。Fig.8から，出力ノイズのために評価関数の値が大きくなり，自然長の変化を検出できなくなっていることが確認できる。

Fig. 7.

Output data obtained from the system with noise. (Online version in color.)

Fig. 8.

Values of the objective function for each k₁ and l₀₁ for the case that noise is added to the output. The lower surface is the one given in Fig.6. (Online version in color.)

出力にノイズが存在する場合，誤差eにローパスフィルタを適用することでノイズの影響を低減させることが考えられる。すなわち，Fig.9のフィルターFの出力zを改めてeとして前節の手法を適用する。ただし，その場合には，フィルターの動特性も含めたヤコビアンを求める必要がある。

Fig. 9.

Filter to attenuate noise.

また，出力y_rの要素ごとのオーダが異なっている場合，誤差eの各要素を重みづけして評価することが行われる。具体的には，式(10)を次式

  

$$J\left(x_e,q_e,u_e\right)=\frac{1}{2}{\displaystyle \sum }_{k=0}^{N-1}{\Vert W_{o}e\left[k\right]\Vert }^2$$(30)

のように変更する。ただし，W_o∈$R^{\tilde{p}}{}^{\times p}$は出力に対する重みである。式(30)はFig.9のフィルターFを重み行列W_oとした特別な場合に対応する。よって，Fig.9のFに一般的な線形システムを仮定してヤコビアンを導出しておけば，フィルターや重みを用いる場合に対して推定手法を適用することが可能となる。

ここで，Fig.9のフィルターFがつぎの離散時間状態方程式で与えられたものとする。

  

$$x_f\left[k+1\right]=A_f{x}_f\left[k\right]+B_{f}e\left[k\right]$$(31)

  

$$z\left[k\right]=C_f{x}_f\left[k\right]+D_{f}e\left[k\right]$$(32)

ただし，x_f[k]∈R^nf，e[k]∈R^p，z[k]∈${\text{R}}^{\tilde{p}}$はそれぞれフィルタの状態，入力，出力である。フィルタの初期状態x_f[0]は既知とする。定数の出力重みとするときにはn_f=0である。このとき，$\frac{1}{2}{\Sigma }_{k=0}^{N-1}{\Vert z\left[k\right]\Vert }^2$を最小化するためのガウス・ニュートン法は，式(12)のeをzに置き換えた次式

  

$$\left[\begin{array}{c}{x}_e^{\left(l+1\right)}\\ q_e^{\left(l+1\right)}\\ u_e^{\left(l+1\right)}\end{array}\right]=\left[\begin{array}{c}{x}_e^{\left(l\right)}\\ q_e^{\left(l\right)}\\ u_e^{\left(l\right)}\end{array}\right]-{\left(M_z^{\left(l\right)}\right)}^{†}\text{vec}\left(z\right)$$(33)

で与えられる。ただし，このときのヤコビアンM^(l)は

  

$$M_z^{\left(l\right)}={\left[\begin{array}{ccc}\frac{\partial vec\left(z\right)}{\partial x\left[0\right]}& \frac{\partial vec\left(z\right)}{\partial q}& \frac{\partial vec\left(z\right)}{\partial u_e}\end{array}\right]|}_{\begin{array}{c}x\left[0\right]=x_e^{\left(l\right)},\\ q=q_e^{\left(l\right)},\\ u_e=u_e^{\left(l\right)}\end{array}}$$(34)

である。

さらに，式(34)のヤコビアンは以下のシステム

  

$${\xi }_f\left[k+1\right]=A_f{\xi }_f\left[k\right]+B_f\eta \left[k\right]$$(35)

  

$${\eta }_f\left[k\right]=C_f{\xi }_f\left[k\right]+D_f\eta \left[k\right]$$(36)

と式(13)，(14)の応答を用いて

  

$$M_z^{\left(l\right)}=\left[\begin{array}{c}{\eta }_f\left[0\right]\\ ⋮\\ {\eta }_f\left[N-1\right]\end{array}\right]$$(37)

で与えられる。ただし，式(35)の初期値は

  

$${\xi }_f\left[0\right]=\left[\begin{array}{ccc}{0}_{n_f}& 0_{n_f\times a}& 0_{n_f\times mN}\end{array}\right]$$(38)

である。実際，式(31)，(32)より

  

$$\begin{array}{l}\left[\begin{array}{ccc}\frac{\partial x_f\left[k+1\right]}{\partial x\left[0\right]}& \frac{\partial x_f\left[k+1\right]}{\partial q}& \frac{\partial x_f\left[k+1\right]}{\partial u_e}\end{array}\right]\\ =A_f\left[\begin{array}{ccc}\frac{\partial x_f\left[k\right]}{\partial x\left[0\right]}& \frac{\partial x_f\left[k\right]}{\partial q}& \frac{\partial x_f\left[k\right]}{\partial u_e}\end{array}\right]+B_f\left[\begin{array}{ccc}\frac{\partial e\left[k\right]}{\partial x\left[0\right]}& \frac{\partial e\left[k\right]}{\partial q}& \frac{\partial e\left[k\right]}{\partial u_e}\end{array}\right]\end{array}$$(39)

  

$$\begin{array}{l}\left[\begin{array}{ccc}\frac{\partial z\left[k\right]}{\partial x\left[0\right]}& \frac{\partial z\left[k\right]}{\partial q}& \frac{\partial z\left[k\right]}{\partial u_e}\end{array}\right]\\ =C_f\left[\begin{array}{ccc}\frac{\partial x_f\left[k\right]}{\partial x\left[0\right]}& \frac{\partial x_f\left[k\right]}{\partial q}& \frac{\partial x_f\left[k\right]}{\partial u_e}\end{array}\right]+D_f\left[\begin{array}{ccc}\frac{\partial e\left[k\right]}{\partial x\left[0\right]}& \frac{\partial e\left[k\right]}{\partial q}& \frac{\partial e\left[k\right]}{\partial u_e}\end{array}\right]\end{array}$$(40)

が成り立ち，さらに

  

$$\begin{array}{l}\left[\begin{array}{ccc}\frac{\partial e\left[k\right]}{\partial x\left[0\right]}& \frac{\partial e\left[k\right]}{\partial q}& \frac{\partial e\left[k\right]}{\partial u_e}\end{array}\right]\\ =\left[\begin{array}{ccc}\frac{\partial y\left[k\right]}{\partial x\left[0\right]}& \frac{\partial y\left[k\right]}{\partial q}& \frac{\partial y\left[k\right]}{\partial u_e}\end{array}\right]=\eta \left[k\right]\end{array}$$(41)

が成り立つ。よって，

  

$${\xi }_f\left[k\right]=\left[\begin{array}{ccc}\frac{\partial x_f\left[k\right]}{\partial x\left[0\right]}& \frac{\partial x_f\left[k\right]}{\partial q}& \frac{\partial x_f\left[k\right]}{\partial u_e}\end{array}\right]$$(42)

  

$${\eta }_f\left[k\right]=\left[\begin{array}{ccc}\frac{\partial z\left[k\right]}{\partial x\left[0\right]}& \frac{\partial z\left[k\right]}{\partial q}& \frac{\partial z\left[k\right]}{\partial u_e}\end{array}\right]$$(43)

とすれば式(35)，(36)が得られる。

式(37)のヤコビアンを用いれば，フィルターを用いた場合においても初期状態・パラメータ・入力の推定を行うことができる。ここで，伝達関数$\frac{1}{3.2\text{s}+1}$を状態空間実現し，サンプル周期Δtで双一次変換に基づく離散化を行って得られたシステムをFとする。このFを用いて式(33)に基づく最適化を行うと，入力の推定値はFig.10に示されるものが得られる。Fig.10より，入力に高周波数帯域の大きな成分があらわれることが確認できる。これは，フィルターによって入力から偏差までの高周波数帯域のゲインが減少したことにより，大きな入力が許容されてしまうためである。このように，出力にノイズが存在する場合，フィルターを設定するだけでは不十分である。

Fig. 10.

Estimated input given in the case that the noise exists. (Online version in color.)

上記の結果からわかるように，フィルターを用いる際には，同時に入力uに高周波数帯域の成分があらわれることを抑制する必要がある。このためには，uのパラメトリゼーションを事前に定められた周波数帯域の信号のみがあらわれるように変更すればよい。周波数帯域を制限した入力のパラメトリゼーションは離散フーリエ逆変換(IDFT)あるいはアダマール逆変換を用いて行うことができる。

IDFTに基づく方法では入力を離散正弦波信号の線形結合でパラメトライズする。正の偶数のNに対してIDFT行列W(N)∈R^N×Nを以下のように定義する。

  

$$\begin{array}{l}{\left(W\left(N\right)\right)}_{\kappa ,\mu }=\\ \{\begin{array}{ll}1\hfill & \left(\mu =1\right)\hfill \\ \cos \left(\left(\kappa -1\right)\frac{\mu }{2}{\omega }_N\right)\hfill & \left(\mu =2,4,\cdots ,N\right)\hfill \\ \sin \left(\left(\kappa -1\right)\frac{\mu -1}{2}{\omega }_N\right)\hfill & \left(\mu =3,5,\cdots ,N-1\right)\hfill \end{array}\end{array}$$(44)

ただし，ω_N=$\frac{2\pi }{N}$である。さらに，各正弦波に対する係数v_μ∈R^m(μ=1,...,N)を用いて入力を

  

$$u\left[k\right]={\displaystyle \sum }_{\mu =1}^N{\left(W\left(N\right)\right)}_{\left(k+1\right),\mu }{v}_{\mu }$$(45)

とパラメトライズする。式(45)は次式のようにコンパクトに表現することができる。

  

$$\begin{array}{l}\text{vec}\left(u\right)=\left[\begin{array}{c}{\displaystyle \sum _{\mu =1}^N{\left(W\left(N\right)\right)}_{1,\mu )}{v}_{\mu }}\\ ⋮\\ {\displaystyle \sum _{\mu =1}^N{\left(W\left(N\right)\right)}_{1,\mu )}{v}_{\mu }}\end{array}\right]\\ =\left[\begin{array}{ccc}{\left(W\left(N\right)\right)}_{1,1}{I}_m& \cdots & {\left(W\left(N\right)\right)}_{1,N}{I}_m\\ ⋮& \ddots & ⋮\\ {\left(W\left(N\right)\right)}_{N,1}{I}_m& \cdots & {\left(W\left(N\right)\right)}_{N,N}{I}_m\end{array}\right]\left[\begin{array}{c}{v}_1\\ ⋮\\ v_N\end{array}\right]\\ =\left(W\left(N\right)\otimes I_m\right)\text{vec}\left(v\right)\end{array}$$(46)

ただし，A⊗BはAとBのクロネッカー積，すなわち

  

$$A\otimes B=\left[\begin{array}{ccc}{a}_{11}B& a_{12}B& \cdots \\ a_{21}B& a_{22}B& \cdots \\ ⋮& ⋮& \ddots \end{array}\right]$$(47)

である。IDFT行列W(N)は正則行列であるので，クロネッカー積の性質[13]からW(N)⊗I_mも正則行列である。よって，信号u[0],...,u[N-1]を定めることとフーリエ係数v₁，t_s，v_Nを定めることは等価である。

W(N)の隣接する偶数・奇数列，すなわち各λ=1,...,$\frac{N}{2}$-1に対してW(N)の第2λ，2λ+1列は周波数λ_ωNの正弦波・余弦波をあらわしている。よって，W(N)の左側から，ある奇数列分だけを用いることにより，低い周波数帯域のみに制限した入力をパラメトライズすることができる。すなわち，あるλ∈{1,...,$\frac{N}{2}$-1}を定め，フーリエ係数v₁,...,v_2λ+1を用いてvec(u)を

  

$$\text{vec}\left(u\right)={\widehat{W}}_{\lambda }\left(N\right)\left[\begin{array}{c}{v}_1\\ v_2\\ ⋮\\ v_{2\lambda +1}\end{array}\right]$$(48)

とパラメトライズすればよい。ただし，

  

$${\widehat{W}}_{\lambda }\left(N\right)=\left(W\left(N\right)\left[\begin{array}{c}{I}_{2\lambda +1}\\ 0\end{array}\right]\right)\otimes I_m$$

である。

uを式(48)でパラメトライズすると，最適化変数はuでなくv₁,...,v_2λ+1となる。この変更に応じて，ヤコビアンを計算する際の式(13)，(14)，(15)も以下のように変更する必要がある。

  

$$\xi \left[k+1\right]=A\left[k\right]\xi \left[k\right]+E\left[k\right]\left[\begin{array}{ccc}{0}_{a\times n}& I_a& 0_{a\times m\tilde{N}}\end{array}\right]+B\left[k\right]\left[\begin{array}{ccc}{0}_{m\times n}& 0_{m\times a}& U_k{\widehat{W}}_{\lambda }\left(N\right)\end{array}\right]$$(49)

  

$$\eta \left[k\right]=C\left[k\right]\xi \left[k\right]+F\left[k\right]\left[\begin{array}{ccc}{0}_{a\times n}& I_a& 0_{a\times m\tilde{N}}\end{array}\right]+D\left[k\right]\left[\begin{array}{ccc}{0}_{m\times n}& 0_{m\times a}& U_k{\widehat{W}}_{\lambda }\left(N\right)\end{array}\right]$$(50)

  

$$\xi \left[0\right]=\left[\begin{array}{ccc}{I}_n& 0_{n\times a}& 0_{n\times m\tilde{N}}\end{array}\right]$$(51)

ただし，$\tilde{N}$=2λ+1である。一方，式(35)，(36)についてはξ_f[k]，η_f[k]のサイズのNを$\tilde{N}$に変更するだけでよい。

4・2　数値例

Fig.7のデータに4・1節の手法を適用する。周波数帯制限を行うために，IDFT行列W(512)のうち，周波数の低いものから10種類の周波数成分のみを用いた。すなわち，λ=10である。これにより得られた初期状態およびパラメータは

  

$$x\left[0\right]={\left[\begin{array}{cccc}1.098& 2.320& 1.301\times {10}^{-2}& -1.661\end{array}\right]}^T$$

  

$$k_1=1.000,k_2=1.000$$

である。上記の結果から，初期状態の推定値にはやや誤差が残るがパラメータは正確な値が得られていることがわかる。一方，得られた入力をFig.11に示す。Fig.11ではFig.10のように大きな値は表れていない。また，周波数帯域を制限しているため真の外乱入力を再現できてはいないが，おおよそ妥当な結果が得られている。

Fig. 11.

Estimated input when the input is parameterized based on the IDFT. (Online version in color.)

この場合にもバネの自然長ℓ₀₁を変化させて出力データを生成し，それらに対して4・1節の手法をそれぞれ適用し，評価関数の値を計算した。得られた評価関数値をFig.12に示す。Fig.12より，ノイズが無い場合よりも曲面は鈍っているが，この場合にもℓ₀₁の変化に対して評価関数値が増大していることが確認できる。

Fig. 12.

Values of the objective function for each k₁ and l₀₁ in the case that the IDFT parameterization is employed. (Online version in color.)

なお，IDFTではなくアダマール逆変換(付録参照)を行った場合に得られる入力をFig.13に示す。IDFTとは異なりアダマール変換では矩形的な入力を生成することが可能である。

Fig. 13.

Estimated input when the input is parameterized based on the inverse Hadamard transform. (Online version in color.)

5. おわりに

設備診断のために初期状態・パラメータ・入力を同時に推定する手法を提案した。任意のタイミングで診断を行うためには，パラメータだけでなく初期状態も推定する必要がある。また，設備は外乱の影響を受けるが，出力の次元を増やすことで外乱の推定は可能である。本論文では，上記の推定問題を最適化問題に帰着させ，ガウス・ニュートン法を実行するために必要なヤコビアンの計算アルゴリズムも提示した。とくに，出力にノイズが含まれる場合，誤差にフィルターを用いることだけでなく，入力の帯域制約も必要であることを指摘し，IDFT あるいはアダマール逆変換を利用することで帯域制約が可能であることを示した。提案した手法の有効性を数値例を用いて検証し，想定したパラメータの推定と，想定外の変化の検出が可能であることを確認した。

付録

アダマール逆変換

アダマール逆変換は，離散的な周期信号を同じ周期の矩形波の合成で表すことのできる変換である¹⁴⁾。基底関数の線形結合によって信号を表現する点はIDFTと同様である。ただし，IDFTの基底関数が三角関数であるのに対し，アダマール逆変換は{-1, 1}の値のみをとる基底関数(ウォルシュ関数)を用いる。

IDFTではなくアダマール逆変換を用いる場合には，交番数順(sequency order)のアダマール行列をH(N)を式(44)のW(N)の代わりに用いる。例としてN=4に対するアダマール行列H(4)を以下に示す。

  

$$H\left(4\right)=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& 1& -1& -1\\ 1& -1& -1& 1\\ 1& -1& 1& -1\end{array}\right]$$(52)

アダマール行列の詳細については文献¹⁴⁾などを参照されたい。

文献

• 1)  R.Isermann: Fault-Diagnosis Systems - An Introduction from Fault Detection to Fault Tolerance, Springer-Verlag, Berlin, (2006). https://doi.org/10.1007/3-540-30368-5
• 2)  R.Isermann: Fault-Diagnosis Applications - Model-Based Condition Monitoring: Actuators, Drives, Machinery, Plants, Sensors, and Fault-tolerant Systems, Springer-Verlag, Berlin, (2011). https://doi.org/10.1007/978-3-642-12767-0
• 3)  巨大構造物ヘルスモニタリング－劣化のメカニズムから監視技術とその実際まで，NTS，東京，(2015).
• 4)  足立修一，丸田一郎:カルマンフィルタの基礎，東京電機大学出版局，東京，(2012).
• 5)  L. Lee and K. Poolla: J. Dyn. Syst. Meas. Control, 121(1999), 71. https://doi.org/10.1115/1.2802444
• 6)  P.M.Quang，浅井徹，有泉亮:非線形システムの初期状態・パラメータ推定のためのヤコビアンダイナミクスの導出，第4回制御部門マルチシンポジウム，計測自動制御学会，東京，(2017), 2E3-5.
• 7)  M.Fujigaki, T.Hara, N.Ikoma and Y.Murata: J. Jpn. Soc. Exp. Mech., 12(2012), 179 (in Japanese). https://doi.org/10.11395/jjsem.12.179
• 8)  T. Aoyama, L. Li, M. Jiang, K. Inoue, T. Takaki, I. Ishii, H. Yang, C. Umemoto, H. Matsuda, M. Chikaraishi and A. Fujiwara: IEEE/ASME Trans. Mechatron., 23(2018), 179. https://doi.org/10.1109/TMECH.2017.2764504
• 9)  T.Asai, Y.Asano, Y.Nishida and S.Nishino: 34th Chinese Control Conference and SICE Annual Conference 2015, SICE, Tokyo, (2015), PWeC-11.
• 10)  S.Nakamura, T.Sakagami, S.Harada and S.Kubo: Concr. Res. Technol., 19(2008), 51 (in Japanese). https://doi.org/10.3151/crt1990.19.3_51
• 11)  J.Nocedal and S.J.Wright: Numerical Optimization, Springer-Verlag, New York, (1999). https://doi.org/10.1007/978-0-387-40065-5
• 12)  池辺八洲彦，池辺淑子，浅井信吉，宮崎佳典:現代線形代数－分解定理を中心として－，共立出版，東京，(2009).
• 13)  K.Zhou: Robust and Optimal Control, Prentice Hall, New Jersey, (1996).
• 14)  喜安善市:アダマール行列とその応用－構成理論からウォルシュ関数まで－，電子通信学会，東京，(1980).