【倒立振子part28】モータのシステム同定

part26にて求めたモデルを使用してシステム同定を行ってみる。
手探り状態なので間違っているかも。。。
今回システム同定をするにあたりMATLABが必要となりそうなのでお試しで評価版を使用していく。
SITB含めても2万ちょいなので買ってもいいかもしれない。

プリ同定
- 無負荷状態
- ロータ装着状態
同定実験の設計
- シフトレジスタの個数n
- 信号の振幅
システム同定の進め方
システム同定(無負荷状態)
- システム同定の比較
- 整合性
システム同定(ロータ装着状態)
- システム同定の比較
- 整合性
パラメータ計算
状態空間モデルの導出
まとめ
参考にさせていただいた資料

プリ同定

今回使用したモータはID-549XWである。
まずシステム(今回は回路内蔵のモータ)の支配的な立ち上がりや定常のゲインといった大まかなダイナミクスをつかむことが目的となる。
プリ同定のやり方としてシステムに、ステップ応答試験やランダム加振試験をして応答値を確認する方法があるが今回はステップ応答試験での結果を参照する。
またシステム同定にて求めたモデルに対して結果と実験値を比較して妥当性の検証を行う。

試験方法としてはduty100%(12V印可計算)になるようにモータにステップ入力して応答値(角速度)を確認する。

無負荷状態

ロータ装着状態

システム結果については単位ステップになるように印可電圧で出力を割っている。
予想した通りの応答値になっているのでこれで行くとする。

同定実験の設計

伝達関数を同定するためには各周波数における正弦波入力の応答値からボード線図を使用してカーブフィッティングを行い周波数伝達関数を導出することができるが試験回数が多くなり大変である。
そこで制御対象が持つ全てのモードを励起させるために多数の周波数をもつ信号を印可することで一括で求められる。
一番望ましい入力信号は白色信号であるが物理的には実現できないので疑似的に白色信号を作成する。また線形システム同定には2値信号のみで問題ないので、疑似白色2値信号であるM系列信号を使用してシステム同定を行う。

M系列信号はMATLABでコマンド一つで作ることができる。
その時に使用する引数について検討する。
・シフトレジスタの個数
・信号の振幅

シフトレジスタの個数n

ここでM系列信号の周期 $N(=2^n-1)$ に対して+1(立ち上がり状態)をn回繰り返すのは1回だけなのでその際が最大持続パルスとなり、これが対象の立ち上がり時間 $t_r$ より大きくすることで定常ゲインを正確に出力できる。
つまり、1周期が完了するまでにn個のシフトが行われるため各シフトにサンプリング周期Tの時間がかかるため
$nT>t_r$ ⇒ $n>\frac{t_r}{T}$
ただしM系列信号のサンプリング周期とモータへの入力のサンプリング周期(マイコン側)が同じが一番望ましいがモータへの入力のサンプリング周期が短いためM系列信号のシフトレジスタが大きくなりデータ数が莫大になってしまう。そこでM系列信号のために下記のクロック周波数 $T_m$ を使用することでnをある程度小さくする。
$T_m=pT$
pを増加させることでnは小さくなるが、増加させすぎるとパワースペクトルが一定である周波数帯域が減少するため、一般的には $p\leq4$ ぐらいがちょうど良いらしい。
マイコン側での入力はサンプリング周期Tで行い、データ取りはクロック周期Tmで行うことである程度のデータ数の取得ができるようになる。
今回は立ち上がり時間Trとサンプリング周期Tより下記のようなM系列信号を使用して試験を行うことにする。
またシステム同定する際にデータを半分にして、同定用と検証用として使用するため元の信号に周期を反転させた信号を合わせたものを入力信号とする。

信号の振幅

システムへの入力は小さすぎると出力のSN比(出力とノイズの比率)が劣化される。大きすぎると制御対象によっては出力が飽和してしまい非線形となってしまう。
今回の制御対象はモータであるため、入力である印可電圧と出力である角速度は飽和しないため大きく設定して問題ないと考える。
またモータ単体での評価なので試験場所的にも大きくして問題ないためとりあえずduty100%(12V印可計算)で評価する。

システム同定の進め方

「システム同定の基礎」を読んでいると、データの読み込みや外れ値を削除したりの前処理が必要であるが、MATLABのSystem Identification Toolboxを使えばGUIで簡単に行える。
GUIの使い方
下記にSystem Identification Toolboxについてまとめられているのでこれを参考にする。
jp.mathworks.com
今回はSISOなので下記を参考にGUIの使い方を学んでいく
jp.mathworks.com

システム同定(無負荷状態)

入力データのM系列の条件

サンプリング周期	シフトレジスタ係数	クロック周期	シフトレジスタ個数	データ数
$T(s)$	$p$	$T_{m}(s)$	$n$	$N=2^n-1$
$0.04$	$1$	$0.04$	$7$	$127$

入出力データをまとめると

システム同定の比較

System Identification Toolboxを使用して1次遅れ系(sys1)と2次遅れ系(sys2)でシステム同定を行う。

1次遅れ系の伝達関数
$\displaystyle{G(s)=\frac{1031}{s+37.39}}$

2次遅れ系の伝達関数
$G(s)=\displaystyle{\frac{11550}{s^2+114s+4341}}$

整合性

1次遅れ系の推定データへの適合は94.47%
2次遅れ系の推定データへの適合は99.62%

システム同定(ロータ装着状態)

入力データのM系列の条件

サンプリング周期	シフトレジスタ係数	クロック周期	シフトレジスタ個数	データ数
$T(s)$	$p$	$T_{m}(s)$	$n$	$N=2^n-1$
$0.4$	$1$	$0.4$	$7$	$127$

入出力データをまとめると

システム同定の比較

System Identification Toolboxを使用して1次遅れ系(sys1)と2次遅れ系(sys2)でシステム同定を行う。

1次遅れ系の伝達関数
$\displaystyle{G(s)=\frac{98}{s+3.563}}$

2次遅れ系の伝達関数
$\displaystyle{G(s)=\frac{1084}{s^2+11.06s+40.98}}$

整合性

1次遅れ系の推定データへの適合は93.74%
2次遅れ系の推定データへの適合は99.65%

パラメータ計算

上記の結果から1次遅れ系でも2次遅れ系でもほぼ同等の整合性が得られることからインダクタンス成分を無視して考える。
ここでpart27で求めた無負荷状態でのモータの運動方程式では粘性項を無視していたが必要そうなので追加すると、角速度 $\Omega_m(s)$ と入力電圧 $V_a(s)$ の伝達関数は
$\displaystyle{ G_m(s)=\frac{\Omega_m(s)}{V_a(s)}=\frac{k_t}{sJ_mR+(D+k_bk_t)} }$
となる。
従ってシステム同定の結果より
$\begin{cases} \displaystyle{\frac{k_t}{J_mR}=1031}\\ \displaystyle{\frac{RD+k_tk_b}{J_mR}=37.39}\\ \end{cases}$
となる。
モータの逆起電力定数とトルク定数は同値として $k_b=k_t=0.0385$ であり、回転子の慣性モーメントは $J_m=6.592\times10^-6$ であるため、上式を計算することで未知パラメータが判明する。
$\begin{cases} \displaystyle{R=5.6645}\\ \displaystyle{D=-1.517\times10^{-5}}\\ \end{cases}$

次にロータを装着状態での粘性項をD'とすると同様に計算できる。
$\begin{cases} \displaystyle{J_r=6.2765\times10^{-5}}\\ \displaystyle{D'=-1.4577\times10^{-5}}\\ \end{cases}$

またpart27からpart25で使用する外力の $\tau_omega$ は1次遅れ系のため下記の様に求めることができる。
$\begin{eqnarray} \tau_\omega= \tau_L&=& \tau_m-J_m\dot\omega_m-D\omega\\ &=& k_ti_a(t)-J_m\dot\omega_m-D\omega_m\\ &=& \frac{k_t}{R}v_a-J_m\dot\omega_m-\left(\frac{k_tk_b}{R}+D\right)\omega_m\\ &=& \frac{k_t}{R}v_a-J_m\ddot\theta_r-\left(\frac{k_tk_b}{R}+D\right)\dot\theta_r　　(状態方程式用に変数変換)\\ \end{eqnarray}$

状態空間モデルの導出

part25から運動モデルは
$\begin{cases} J_r\ddot{\theta_r}+(m_r{l_r}^2+mp{l_p}^2+J_r+J_p)\ddot{\theta_p}+\mu_p\dot{\theta_p}-(m_rl_r+m_pl_p)g\theta_p=0\\ J_r\ddot{\theta_r}+\mu_r\dot\theta_r+J_r\ddot{\theta_p}=\tau_\omega \end{cases}$
となるため、上式と組み合わせると
$\begin{cases} J_r\ddot{\theta_r}+(m_r{l_r}^2+mp{l_p}^2+J_r+J_p)\ddot{\theta_p}+\mu_p\dot{\theta_p}-(m_rl_r+m_pl_p)g\theta_p=0\\ J_r\ddot{\theta_r}+D'\dot\theta_r+J_r\ddot{\theta_p}=\frac{k_t}{R}v_a-J_m\ddot\theta_r-(\frac{k_tk_b}{R}+D)\dot\theta_r \end{cases}$
これはpart4と同じ形のため
$\left(\begin{array}{c}\dot\theta_p\\ \ddot\theta_p\\\dot\theta_r\\\ddot\theta_r\end{array}\right)= \begin{pmatrix}0 & 1 & 0 & 0\\\frac{a_{22}(m_rl_r+m_pl_p)g}{\Delta} & \frac{-a_{22}\mu_p}{\Delta} & 0 & \frac{a_{12}(D'+D+\frac{k_tk_e}{R})}{\Delta}\\0 & 0 &0 & 1\\\frac{-a_{21}(m_rl_r+m_pl_p)g}{\Delta} & \frac{a_{21}\mu_p}{\Delta} & 0 &\frac{-a_{11}(D'+D+\frac{k_tk_e}{R})}{\Delta}\end{pmatrix} \left(\begin{array}{c}\theta_p\\ \dot\theta_p\\\theta_r\\\dot\theta_r\end{array}\right)+\left(\begin{array}{c}0\\\frac{-a_{12}\frac{k_t}{R}}{\Delta}\\0\\\frac{a_{11}\frac{k_t}{R}}{\Delta}\end{array}\right)v_a$

ただし
$\begin{eqnarray} \left\{\begin{array}{l} a_{11}=m_r{lr}^2+mp{lp}^2+J_r+J_p\\ a_{12}=J_r\\ a_{21}=J_r\\ a_{12}=J_r+J_m\\ \Delta=a_{11}a_{22}-a_{12}a_{21}\\ \end{array} \right. \end{eqnarray}$
とする。

求めたパラメータを使って計算すると
$\displaystyle{ A=\begin{bmatrix}0 & 1 & 0 & 0 \\43.3736 & -0.2079 & 0 & 0.0352 \\ 0 & 0 & 0 & 1 \\-39.2521 & 0.1881 &0 & -3.3758 \end{bmatrix} }$
$\displaystyle{ B=\begin{bmatrix}0\\-1.0312\\0\\98.9293 \end{bmatrix} }$