システムの同定

ARMAモデル

システムの挙動を過去の入力と出力の重ね合わせで表現したモデルは、ARMAモデルとよびます。
$$y(n)=\sum_{k=1}^{p} a(k) y(n-k)+\sum_{k=1}^{q} b(k) u(n-k)$$
このブログでは、逐次式との比較が容易なARMAモデルによるシステム同定を考えていきます。ARMAモデルのシステム同定とは、具体的にはシステムの入力に対する出力の挙動データから、数式モデルの 𝑎(𝑘),  𝑏(𝑘)  を決定することになります。

EXCELファイルのダウンロード

こちらからシステム同定のEXCELファイルをダウンロードできます。
システム同定のEXCELファイルには、システムの入力と出力の応答データから、ARMAモデルのパラメータ 𝑎, 𝑏 を最小二乗法で逐次的に推定する処理を組み込んでいます。
車両走行モデルやバネマスモデルの入力と出力の応答データから、システム同定する例を各シートに示しています。

入力と出力の応答データ U(N), Y(N) に別のシステムの応答データを設定すれば、そのシステムのARMAモデルのパラメータ 𝑎, 𝑏 を求めることができます。

システム同定Excelシートの説明

こんな感じで動きます

畳み込み演算（コンボリューション）

システム応答の畳み込み演算（convolution）をEXCELシート上で確認します。
３次のモデルの場合、過去３回の入出力データと a, b のパラメータの畳み込み演算で次の出力を計算します。EXCELシート上では下の図のような計算になります。

システム同定の確認

車両走行モデル　～加減速の挙動

システム同定（走行 加速度）のシートでは、状態方程式でモデル化した車両の加減速の挙動データを用いて車両の加減速のARMAモデルを同定します。

車両走行モデルの加速度は１次遅れのモデルにしたので、次数は１（パラメータ数２）でARMAのパラメータ同定を実行しています。1次遅れのモデルについてはこちらを参照して下さい。1次遅れモデルで $$e^{-\frac{\Delta t}{T}}=\alpha$$ とすると逐次式は

　  𝑥(𝑘+1) = α 𝑥(𝑘)+(1−α)𝑢(𝑘)

時定数1sec,  離散時間100ms より
$$\mathrm{A} 1=\alpha=e^{-\frac{\Delta t}{T}}=0.9048 \ldots$$ $$\mathrm{B} 1=1-\alpha=0.0951 \ldots$$ であり、１次遅れの真値と比べると若干の誤差はあるものの、グラフからも分かるように、ほぼ同定できています。

車両走行モデル　～速度の挙動

システム同定（走行 速度）のシートでは、状態方程式でモデル化した車両の加減速の挙動データを用いて、車両の速度のARMAモデルを同定します。特性が分かっていないシステムを同定する場合は、次数を増やしたり、演算回数を増やしてみて、同定パラメータの誤差の収束状況を確認します。不必要にパラメータの次数は上げず、なるべく少ない次数でのシステム同定を目指します。

バネマスモデル　～変位速度の挙動

システム同定（バネマス 変位速度）のシートでは、状態方程式でモデル化したバネマスモデルの荷重変化時の変位速度の挙動データを用いてARMAモデルを同定します。

バネマスモデル　～変位量の挙動

システム同定（バネマス ストローク）のシートでは、状態方程式でモデル化したバネマスモデルの荷重変化時の変位量(ストローク量)の挙動データを用いてARMAモデルを同定します。
挙動データの変化が少ない場合、逐次推定で必要となる演算回数が増えます。繰り返し演算回数を増やすと演算時間も長くなるので注意が必要です。

ARMAモデルパラメータの逐次推定方法

ARMAモデル
$$y(n)=\sum_{k=1}^{p} a(k) y(n-k)+\sum_{k=1}^{q} b(k) u(n-k)$$ において
$$x=\left[a_{1}, a_{2}, \cdots \cdot, a_{p}, b_{1}, b_{2}, \cdots \cdot, b_{q}\right]^{T}$$ $$z(n)=[y(n-1), \cdots, y(n-p), u(n-1), \cdots, u(n-q)]^{T}$$ $$Z=\begin{bmatrix} z^{T}(n+1) \\ \vdots \\ z^{T}(N) \\  \end{bmatrix}$$ $$P_N=({Z_N}^{T} Z_N)^{-1}$$ としたときに、最小二乗法で逐次的に 𝑥 を推定するアルゴリズムは $$\hat{x}_{n+1}=\hat{x}_{n}+g(n+1)\left\{\mathrm{y}(n+1)-z^{T}(n+1) \hat{x}_{n}\right\}$$ $$g(n+1)=\frac{P_{n} z(n+1)}{1+z^{T}(n+1) P_{n} z(n+1)}$$ $$P_{n}=\left[I-g(n) z^{T}(n)\right] P_{n-1}$$ です。

EXCELには、このアルゴリズムをVBAのmodule1に実装しています。EXCELの「開発」タグを選択し、「Visual Basic」からプログラムを確認できます。

詳細は、参考図書などで確認お願いします。