【倒立振子part5】パラメータ同定について

前回で状態空間モデルが完成したので今回は変数の値を求めていく。
重量など実測することで判明するものもあるが、粘性摩擦抵抗などのように測定できないものは近似モデルと実測が同じ結果になるようにパラメータ同定を行う。

1. 慣性ロータ系のパラメータ同定
- 1.1 DCモータの定数と内部抵抗測定
- 1.2 慣性ロータ系のパラメータ同定
2.振子のパラメータ同定
結論
次回

1. 慣性ロータ系のパラメータ同定

1.1 DCモータの定数と内部抵抗測定

慣性ロータ系のパラメータ同定をするために、使用するモータの逆起電力定数 $k_b$ 、トルク定数 $k_t$ 、端子間抵抗 $R$ を求める。
maxonなどの高級モータだとデータシートに記載されているが、
今回はAliexpressで購入した詳細が不明であるフライホイール付きのDCモータを使用する。
https://ja.aliexpress.com/item/4000314159240.html?spm=a2g0o.order_list.order_list_main.65.21ef585aU55amr&gatewayAdapt=glo2jpn

" $\tau_\omega=k_t\cdot{I_m}$ "より「 $\tau_\omega-I_m$ 特性」を得ることでトルク定数が求まる。
測定のために図のような治具を作成しトルクの測定を行う。
また無負荷状態での一定回転状態にて" $\omega=\frac{1}{k_b}V_m+{切片}$ "より「 $\omega-V_m$ 特性」を得ることで逆起電力定数が求まる。
電圧と電流値は安定化電源の値を使用し回転数については回転計、トルクはおもりの計算から導出する。

試験結果を下記に示す。

試験結果のグラフより
$\displaystyle{k_t=0.00799{N\cdot{m}/A}　,　k_b=\frac{1}{110.33273}=0.00906{V\cdot{s}/rad}}$

モータの電気的特性を考慮すると
$V_m=R\cdot{I_m}+k_b\cdot\omega$
となるので端子間抵抗は
$\displaystyle{R=\frac{V_m-k_b\cdot\omega}{I_m}=\frac{4-0.00906\times395.736}{0.396}=1.047{Ω}}$

1.2 慣性ロータ系のパラメータ同定

前回導出した下記の運動方程式で慣性ロータのみを考えるため振子が安定している状態、つまり $\theta_p=0$ の時についての運動方程式は
$(J_r+J_m)\ddot{\theta_r}+(\mu_r+\mu_m+\frac{k_tk_e}{R})\dot\theta_r=\frac{k_t}{R}v_a$
慣性ロータの数式モデルを下記に変更
$\ddot{\theta_r(t)}=-a\dot\theta_r(t)+bv_a(t)$
ただし
$\displaystyle{ \begin{eqnarray} \left\{\begin{array}{l} \displaystyle{a=\frac{\mu_r+\mu_m+\frac{k_tk_e}{R}}{J_r+J_m}}\\ \displaystyle{b=\frac{\frac{k_t}{R}}{J_r+J_m}} \end{array} \right. \end{eqnarray}}$
とする。

伝達関数を求めるためにラプラス変換すると
$\displaystyle{\theta_r(s)=P(s)v_a(s)　,　P(s)=\frac{b}{s(s+a)}}$
となる。このシステムを同定するためにPコントローラを使用する。
$v(t)=K_pe(t)　\Leftrightarrow　v(s)=K_pe(s)$

このP制御したシステムの伝達関数はバネダンパー系の2時遅れシステムと同等のため固有角周波数と減衰係数が以下のように定義できる。
$\displaystyle{G(s)=\frac{bK_p}{s^2+as+bk_p}=\frac{\omega_{n1}^2}{s^2+2\zeta_1\omega_{n1}s+\omega_{n1}^2}}$

$\displaystyle{ \begin{eqnarray} \left\{\begin{array}{l} 固有角周波数：\omega_{n1}=\sqrt{bk_p}\\ 減衰係数：\displaystyle{\zeta_1=\frac{a}{2\omega_{n1}}=\frac{a}{2\sqrt{bk_p}}}\\ \end{array} \right. \end{eqnarray}}$

" $k_p\rightarrow大$ "とすると $\displaystyle{\zeta_1=\frac{a}{2\sqrt{bk_p}}\rightarrow0}$ この時出力はオーバーシュートする。
入力限界 $v_{max}$ を超えずに目標値 $r_c$ に最適な比例ゲインk_pの設計方法は下式で求められる。
$\displaystyle{|v(0)|=k_p|e(0)|=k_p|r_c|\leq{v_{max}}\Rightarrow{k_p}\leq\frac{v_{max}}{|r_c|}}$
このシステムの流れは下図の通りである。

パラメータ同定のイメージは出力結果から $A_{max},T_p　\Rightarrow　\zeta_1,\omega_{n1}　\Rightarrow　a,b$ を決定する。
ここで
$\displaystyle{ \begin{eqnarray} \left\{\begin{array}{l} \displaystyle{T_p=\frac{\pi}{\omega_{n1}\sqrt{1-\zeta_1^2}}}\\ \displaystyle{A_{max}=r_ce^{-\zeta_1\omega_{n1}T_p}} \end{array} \right. \end{eqnarray}}$ 　　　⇒　　　 $\begin{eqnarray} \left\{\begin{array}{l} \displaystyle{\xi_1=-\frac{1}{T_p}\log_e\frac{A_{max}}{r_c}}\\ \displaystyle{\omega_{n1}=\sqrt{\left(\frac{\pi}{T_p}\right)^2+\xi_1^2}}\\ \displaystyle{\zeta_1=\frac{\xi_1}{\omega_{n1}}} \end{array} \right. \end{eqnarray}$
したがってシステムの未知パラメータ $a,b$ が求まる。
$\displaystyle{ \begin{eqnarray} \left\{\begin{array}{l} \displaystyle{a=2\zeta_1\omega_{n1}=\frac{\mu_r+\mu_m+\frac{k_tk_e}{R}}{J_r+J_m}}\\ \displaystyle{b=\frac{\omega_{n1}^2}{k_p}=\frac{\frac{k_t}{R}}{J_r+J_m}} \end{array} \right. \end{eqnarray}}$
以上より、物性パラメータは
$\displaystyle{ \begin{eqnarray} \left\{\begin{array}{l} \displaystyle{\mu_r+\mu_m+\frac{k_tk_e}{R}=a({J_r+J_m})}\\ \displaystyle{Jm=\frac{\frac{k_t}{R}}{b}-J_r} \end{array} \right. \end{eqnarray}}$

2.振子のパラメータ同定

前回導出した下記の運動方程式で振子のみを考えるため慣性ロータが止まっている状態 $\ddot\theta_r=0$ で振子の基準を $\theta_p=\phi_p+\pi$ について考える(安定系システムのため同定しやすい)
その時のモデルを下記に示すが一般的な振子システムと同等なモデルとなる。

$(\frac{1}{4}m_r{lr}^2+\frac{1}{4}mp{lp}^2+J_r+J_p)\ddot{\phi_p}+\mu_p\dot{\phi_p}+(m_rl_r+m_pl_p)g\phi_p=0$
ここで
$\begin{eqnarray} \left\{\begin{array}{l} \displaystyle{\omega_{n2}^2=\frac{(m_rl_r+m_pl_p)g}{(\frac{1}{4}m_r{lr}^2+\frac{1}{4}mp{lp}^2+J_r+J_p)}}\\ \displaystyle{2\zeta_2\omega_{n2}=\frac{\mu_p}{\frac{1}{4}m_r{lr}^2+\frac{1}{4}mp{lp}^2+J_r+J_p}}\\ \end{array} \right. \end{eqnarray}$
と定義すると上記のシステムは
$\ddot\phi_p+2\zeta_2\omega_{n2}\dot\phi_p+\omega_{n2}^2\phi_p=0$
となり、これは一般的な振子のシステムである。
この微分方程式を $\phi_p(0)=\phi_{p,0}　,　\dot\phi(0)=0$ として解くと
$\displaystyle{\phi_p(t)=e^{-\zeta_2\omega_{n2}t}\left(\cos\omega_{d2}t+\frac{\zeta_p}{\sqrt{1-\zeta_p^2}}\sin\omega_{d2}t\right)\phi_{p,0}}$
となり時間変化グラフは下記の様になることが分かる。

またこのグラフで振動周期 $T\left(=T_1=T_2=･･･\right)$ と減衰率 $\lambda\left(=\frac{A_2}{A_1}=\frac{A_3}{A_2}=･･･\right)$ は一定になるという特性がある。
$\begin{eqnarray} \left\{\begin{array}{l} \displaystyle{T=\frac{2\pi}{\omega_{n2}\sqrt{1-\zeta_2^2}}}\\ \displaystyle{\lambda=e^{-\zeta_2\omega_{n2}T_p}} \end{array} \right. \end{eqnarray}$ 　　　⇒　　　 $\begin{eqnarray} \left\{\begin{array}{l} \displaystyle{\xi_2=-\frac{1}{T}\log_e\lambda}\\ \displaystyle{\omega_{n2}=\sqrt{\left(\frac{2\pi}{T}\right)^2+\xi_2^2}}\\ \displaystyle{\zeta_2=\frac{\xi_2}{\omega_{n2}}} \end{array} \right. \end{eqnarray}$
以上より未知であった物理パラメータは下記によって同定が可能になる。
$\begin{eqnarray} \left\{\begin{array}{l} \displaystyle{\left(\frac{1}{4}m_r{lr}^2+\frac{1}{4}mp{lp}^2+J_r+J_p\right)=\frac{(m_rl_r+m_pl_p)g}{\omega_{n2}^2}}\\ \displaystyle{\mu_p=2\zeta_2\omega_{n2}\left(\frac{1}{4}m_r{lr}^2+\frac{1}{4}mp{lp}^2+J_r+J_p\right)}\\ \end{array} \right. \end{eqnarray}$

結論

上記のようなDCモータ内部のパラメータ同定、慣性ロータのパラメータ同定、振子のパラメータ同定によって状態空間モデルでの不明点であった物理量を測定が可能となる。

次回

使用する部品のプログラムと本体の設計と開発を進める。

P.S
倒立振子で学ぶ制御工学の本を読みながらやっているが、1ページどころか1行のレベルで読みごたえがあって流し読みをしてると理解できない。
ただ理解できた時の知識が晴れ渡る感覚は楽しくなる。
今はまだ基本知識がある程度あるから理解できるけど、離散化とかLMI等のあんまり分かってないところはもっと時間がかかりそうな気がしてる。。。