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

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

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

プログラムの構造

 今回は生データからarxモデルを導出し、それを低次元化するまでのプロセスを一貫して実行するexeファイルを作るのではなく、そのプロセスを以下のような5つのサブプロセスに分割し、各々のサブプロセスに対して1つのexeファイルができるようなプログラム構造としました。

  • 前処理
  • むだ時間推定
  • arxモデル同定
  • 極・零点の数値計算
  • 低次元化モデルの出力

 また、各処理には以下に述べるように関数のプログラムを作り、紐づけをしました。

前処理に用いる関数

 直流成分除去のため入力・出力の生データから平均値を差し引くわけなので、平均値を計算する機能と、全データから指定の値を差し引く機能を関数にします。直流成分を除去した後はLPFで高周波成分を除去するので、LPFも関数にします。

 以上をまとめると、前処理のために関数化する機能は以下のようになります。

実装する関数

  • 平均値の算出
  • データから指定の値を差し引く処理
  • LPF

むだ時間推定に用いる関数

 むだ時間推定に相関解析法を用いるので、この相関解析法を関数にします。よって、むだ時間推定のために関数化する機能は以下のようになります。

実装する関数

  • 相関解析法

arxモデル同定に用いる関数

 まず計算プロセスの中で逆行列を求める場面と、出力次数と入力次数の最大値・最小値を選ぶ場面があるので、これらを関数にします。

 続いて、式(2.17)の$\boldsymbol{R}(\boldsymbol{\theta})$と式(2.18)の$\boldsymbol{f}(\boldsymbol{\theta})$の元となる$\boldsymbol{\varphi}(k)$の計算を関数にします。さらに関数で求めた$\boldsymbol{\varphi}(k)$を用いて、$\boldsymbol{R}(\boldsymbol{\theta})$と$\boldsymbol{f}(\boldsymbol{\theta})$を求める処理も関数にします。

 上記の関数を作れば、$\boldsymbol{R}(N)$と$\boldsymbol{f}(N)$を以下の式にあてはめれば求めることでarxモデルの係数ベクトル$\boldsymbol{\theta}$を求めることができます。

\begin{equation}
\boldsymbol{\theta}
=
\boldsymbol{R}(N)^{-1}\boldsymbol{f}(N)
\tag{2.25}
\end{equation}

式(2.25)の計算処理を各次数組に対して行う処理を関数にします。また、すべての次数組に対して求めたarxモデルの係数ベクトル$\boldsymbol{\theta}$を行列にまとめる処理を関数にします。

 上記の関数で候補の次数組全ての係数ベクトル$\boldsymbol{\theta}$を求められたので、これらの係数とモデル検証用データを用いて$k$番目の推定値を計算する処理を関数にします。さらに、この関数を参照して入出力次数候補全ての推定値を行列にまとめる処理を関数にします。

 上記の関数で求めた推定値とモデル検証用データを用いて損失関数を求める処理を関数にします。またそれらの損失関数の中で、最小の損失関数値を示す入力・出力の次数と、そのときのarxモデルの係数を出力する機能を関数にします。

 以上をまとめると、arxモデル同定のための関数は以下のようになります

実装する関数

  • 行列の逆行列計算
  • 出力次数・入力次数のどちらか大きい方を選ぶ処理
  • 出力次数・入力次数のどちらか小さい方を選ぶ処理
  • $\boldsymbol{\varphi}(k)$を求める処理
  • $\boldsymbol{f}(N)$を求める処理
  • $\boldsymbol{R}(N)$を求める処理
  • 1つの入力・出力次数組に対してarxモデルの係数を求める処理
  • 全ての入力・出力次数組に対するarxモデルの係数を行列にまとめる処理
  • 1つの入力・出力次数組のarxモデルから$k$番目の推定値を求める処理
  • 全ての入力・出力次数組の全時間にわたる推定値を行列にまとめる処理
  • 推定値とモデル検証用データから損失関数を求める処理
  • 損失関数が最小となる入力・出力次数とarxモデルの係数を出力する処理

極・零点の数値計算に用いる関数

 極・零点を求める方法であるDKA法を関数にします。そのDKA法には、該当する多項式と導関数の値それぞれを求める機能を関数にします。

 以上をまとめると、極・零点の数値計算を求めるための関数は以下のようになります

実装する関数

  • DKA法で多項式の根を求める処理
  • 多項式の値を求める処理
  • 多項式の導関数を求める処理

低次元化モデル出力に用いる関数

 ここでやりたい処理は、「低次元化モデルの入力・出力次数からarxモデルの係数・推定値・損失関数を出力する」というシンプルなものなので、この処理のために関数を作るということはしませんでした。

実装する関数

無し