C++によるarxモデル同定プログラムの作成

 実験データからarxモデルを得るプログラムをC++で作ったことに対する解説記事です。
理論的な話も含まれるので、手っ取り早くソースコードが見たい方はこちらから飛んでください。

 記事の構成は以下の通りです

各処理の詳細

前処理

低周波外乱の除去

  • データのオフセット
  • データのトレンド

といった低周波外乱は精度の高いモデルを得るのに好ましくないので、事前に除去しておきます。ここでは低周波外乱の中でも直流成分(=オフセット)のみの除去を目的に、全データから平均値を差し引きます。
 つまり、

\begin{eqnarray}
入力データ: u_1(k)=u_{raw}(k)-\frac{1}{N}\sum_{k=1}^N u_{raw}(k)
\tag{2.1}
\end{eqnarray}

\begin{equation}
出力データ: y_1(k)=y_{raw}(k)-\frac{1}{N}\sum_{k=1}^N y_{raw}(k)
\tag{2.2}
\end{equation}

というアルゴリズムをプログラムに実装します。

高周波外乱の除去

  • センサノイズ
  • モデリングで無視する高次モード

のような高周波成分も、事前に除去しておかなければモデルの精度を悪化させます。今回は入出力データにローパスフィルタ(LPF)を通すことにより高周波成分を除去します。
 つまり、

\begin{equation}
入力データ: u_2(k)=au_2(k)+(1-a)u_1(k)
\tag{2.3}
\end{equation}

\begin{equation}
出力データ: y_2(k)=ay_2(k)+(1-a)y_1(k)
\tag{2.4}
\end{equation}

\begin{equation}
a=\frac{T}{\Delta T-T}
\tag{2.5}
\end{equation}

\begin{equation}
T=\frac{1}{2\pi f}
\tag{2.6}
\end{equation}

というアルゴリズムをプログラムに実装します。ここで、

  • $f$:カットオフ周波数[Hz]
  • $\Delta T$:サンプリング周期[sec]
  • $T$:ローパスフィルタの時定数[sec]

です。

データ分割

 上記2つの前処理を行った後、クロスバリデーションのために入出力データを

  • モデル構築用データ:$u_{21}$、$y_{21}$
  • モデル検証用データ:$u_{22}$、$y_{22}$

に2分割します。この時、過渡状態の入出力データをシステム同定プロセスから省くために測定開始から1秒経過後のデータを2分割します。

むだ時間推定

 むだ時間を推定する方法として最も簡単な方法は、インパルス応答を見ることです。このインパルス応答は、相関解析法を用いれば一般的な入力を用いた実験データであってもインパルス応答を推定することが可能です。
 相関解析法によるインパルス応答$\hat{g}(\tau)$の推定値は、以下に示す式で処理することで得られます。

\begin{equation}
\hat{g}(\tau)=\frac
{
\sum_{k=1}^N y_{21}(k+\tau)u_{21}(k)
}
{
\sum_{k=1}^N u_{21}^2(k)
}
\tag{2.7}
\end{equation}

この処理をプログラムに実装しました。このインパルス応答の推定値を見ることにより、むだ時間推定を行います。

arxモデル同定

 arxモデルによる入出力関係を以下に示します。

\begin{equation}
\begin{split}
y(k)=&-a_1y(k-1)-a_2y(k-2)-\cdots-a_{n_a}y(k-n_a) \\
&+b_1u(k)+b_2u(k-1)+\cdots+b_{n_b}u(k-n_b+1) \\
&+w(k)
\end{split}
\tag{2.8}
\end{equation}

arxモデルは、以下のような図で表すこともあります。

ここで、

\begin{equation}
A(q)=1+a_1q^{-1}+ \cdots + a_{n_a}q^{-n_a}
\tag{2.9}
\end{equation}

\begin{equation}
B(q)=b_1+ b_2q_{-1} \cdots + b_{n_b}q^{-n_b}
\tag{2.10}
\end{equation}

です。

 プログラムの目的であるarxモデルを求めるということは、上式の

  • 入力次数$n_a$と出力次数$n_b$
  • 入力係数$b_1 \cdots b_{n_b}$と出力係数$a_1 \cdots a_{n_a}$

を求めることです。このうち、入出力係数は最小二乗法により解析的に求められますが、入出力次数を解析的に求めることはできません。よってあらかじめ次数の候補をいくつか列挙し、その候補全てに対して

  • 入出力次数組の係数同定
  • 推定値の算出
  • 損失関数の算出

を行い、その中で最も損失関数の小さいモデルを、システムの振る舞いを最も正確に記述するarxモデルとして選定することとしました。

入出力各次数の係数同定

 まず表現の簡略化のため、arxモデルによる1段先予測値$\hat{y}_{21}(k)$

\begin{equation}
\begin{split}
\hat{y}_{21}(k)=&-a_1y_{21}(k-1)-a_2y_{21}(k-2)-\cdots-a_{n_a}y_{21}(k-n_a) \\
&+b_1u_{21}(k-1)+b_2u_{21}(k-2)+\cdots+b_{n_b}u_{21}(k-n_b)
\end{split}
\tag{2.11}
\end{equation}

をベクトル

\begin{equation}
\boldsymbol{\theta}=
\begin{bmatrix}
a_1 && \cdots && a_{n_a} && b_1 && \cdots b_{n_b}
\end{bmatrix}^T
\tag{2.12}
\end{equation}

\begin{equation}
\boldsymbol{\varphi}=
\begin{bmatrix}
-y_{21}(k-1) && \cdots && -y_{21}(k-n_a) && u_{21}(k-n_k) && \cdots && u_{21}(k-n_k-n_b+1)
\end{bmatrix}^T
\tag{2.13}
\end{equation}

を使い、

\begin{equation}
\hat{y}_{21}(k,\boldsymbol{\theta})=\boldsymbol{\theta}^T\boldsymbol{\varphi}
\tag{2.14}
\end{equation}

と表現しなおします。$n_k$はむだ時間です。これにより、モデル構築用データと一段先予測値との誤差は、

\begin{equation}
\begin{split}
\epsilon(k,\boldsymbol{\theta})&=y_{21}(k)-\hat{y}_{21}(k,\boldsymbol{\theta})
&=y_{21}(k)-\boldsymbol{\theta}^T\varphi
\end{split}
\tag{2.15}
\end{equation}

となり、この誤差の二乗和は

\begin{equation}
k_{start}=max(na,n_b+n_k)
\end{equation}

\begin{equation}
J_N(\boldsymbol{\theta})=
\frac{1}{N-k_{start}}\sum_{k=k_{start}}^N
\left\{
y_{21}(k)-\boldsymbol{\theta}^T\boldsymbol{\varphi}
\right\}
\tag{2.16}
\end{equation}

となります。ここでさらに

\begin{equation}
\boldsymbol{R}(N)
=
\frac{1}{N-k_{start}}\sum_{k=k_{start}}^N \boldsymbol{\varphi}(k)\boldsymbol{\varphi}^T(k)
\tag{2.17}
\end{equation}

\begin{equation}
\boldsymbol{f}(N)
=
\frac{1}{N-k_{start}}\sum_{k=k_{start}}^N \boldsymbol{\varphi}(k)y_{21}(k)
\tag{2.18}
\end{equation}

\begin{equation}
c(N)
=
\frac{1}{N-k_{start}}\sum_{k=k_{start}}^N y_{21}(k)^2
\tag{2.19}
\end{equation}

を使って書き換えると、

\begin{equation}
J_N(\boldsymbol{\theta})
=
\boldsymbol{\theta}^T\boldsymbol{R}(N)\boldsymbol{\theta}

2\boldsymbol{\theta}^T\boldsymbol{f}(N)
+
c(N)
\tag{2.20}
\end{equation}

となります。これを$\boldsymbol{\theta}$について微分すると、

\begin{equation}
\nabla J_N(\boldsymbol{\theta})
=
2\boldsymbol{R}(N)\boldsymbol{\theta}-2\boldsymbol{f}(N)
\tag{2.21}
\end{equation}

となります。これが$0$となるときの$\boldsymbol{\theta}$が、その次数組に最も適した係数です。

 以上のように$R(N)$、$f(N)$からarxモデルの係数ベクトルを出力するプログラムを作成しました。

推定値の算出

 次数組に対するarxモデルができたので、そのすべてのarxモデルとモデル検証用データ$y_{22}(k)$を用いて推定値$\hat{y}_{22}(k)$を求めるプログラムを作成します。

損失関数の算出

 算出した推定値$\hat{y}_{22}(k)$とモデル検証用データ$y_{22}(k)$で、

\begin{equation}
損失関数:
V
=
\frac{1}{N}\sum_{k=1}^N
\left\{
\hat{y}_{22}(k)-y_{22}(k)
\right\}^2
\tag{2.22}
\end{equation}

を計算します。その後は全候補の中で最も損失関数の小さいモデルを、同定対象に夫も適合するarxモデルとして選定します。

 作成するプログラムには、

  • 全候補の損失関数
  • 損失関数の最も小さいarxモデルの入出力次数・各項の係数

を出力させることとしました。

 このように実験データを2分割し、一方のデータでモデルを構築、他方のデータでモデルの評価を行うことをクロスバリデーションと呼びます。

モデルの低次元化

 上記の計算プロセスでは傾向的に高次モデルが選ばれてしまいます。そこで、得られたモデルの極・零点を求め、極零相殺の有無を示してモデルの低次元化を試みます。

極・零点の数値計算

 極・零点はDKA法を用います。DKA法とは、$n$次多項式の持つ$n$この会の内$j$番目の解を、初期値

\begin{equation}
z_j^0
=
-\frac{1}{n}\frac{c_1}{c_0}
+
r\exp\
\left\{
\sqrt{-1}
\left(
\frac{2\pi}{n}(j-1)+\frac{\pi}{2n}
\right)
\right\}
\tag{2.23}
\end{equation}

と、ニュートン法の式

\begin{equation}
z_j^{k+1}
=
z_j^k-\frac{f(z_j^k)}{f'(z_j^k)}
\tag{2.24}
\end{equation}

を用いて根を求める方法です。ここで、$c_0$、$c_1$はそれぞれ$n$乗項、$n-1$乗項の係数、$r$は任意の値となります。

 この方法を分子・分母に対して行い、極・零点を出力するプログラムを作成しました。

低次元モデルの出力

 求めた極・零点の出力結果を複素平明王にプロットすれば、いくつの極・零点の組が相殺されうるかがわかります。その結果から低次元化モデルの入力・出力の次数を指定すると、その次数についての

  • arxモデル(入出力次数と各項の係数)
  • 推定値
  • 損失関数

を出力するプログラムを作成しました。