【倒立振子part10】ブラシレスDCモータのパラメータ同定part2

前回はブラシレスDCモータの逆起電力定数、トルク定数、内部抵抗の測定を行った。
今回は粘性摩擦抵抗、回転子の慣性モーメントの導出を行う。

導出方法
ステップ応答実験
パラメータ同定
- シュミレーション結果
次回

導出方法

今回は「DCモータのモデルにステップ入力したときの実験値とシュミレーションの差の二乗和を求めてパラメータ同定を行う」という方法で進めていく。
モータのモデルについては下記にて導出した。
marurunru.hatenablog.com
導出した式の中で今回は無負荷のため $\tau_L=0$ としてモデルをまとめると

$\displaystyle{J_m\ddot\theta_r(t)+\left(\mu_m+\frac{k_tk_e}{R}\right)\dot\theta_r(t)=\frac{k_t}{R}v_a(t)}$
となる。式の処理を簡単にするために変数を整理したモデルを作成すると
$\displaystyle{\dot\omega_r(t)+a\cdot\omega_r(t)=b\cdot{v_a(t)}}$
ただし
$\displaystyle{a=\frac{\left(\mu_m+\frac{k_tk_e}{R}\right)}{J_m}　,　b=\frac{\frac{k_t}{R}}{J_m}}$
となる。今回は応答値として角速度での評価をするため角速度に対してのモデルに変換する。
このモデルの伝達関数は
$\displaystyle{P(s)=\frac{b}{s+a}}$
でありこれは一時遅れ系の伝達関数と同じなので最適伝達関数の一時遅れ系にステップ入力したときの応答値と実測値が合致するような未知パラメータ $a,b$ を求める。

ステップ応答実験

パラメータ同定を行うためにモータにステップ入力した際の応答値(時間-角速度)の実験を行う。
下記で作成したプログラムを元にステップ応答値を求める。

今回は安定化電源の電圧を8V、Duty比を0%(フルDuty)とした。
本来なら実験専用のプログラムを作ったほうがいいのだろうが、データを取れたらいいという考えのもと垂れ流しのシリアル通信からデータをコピーして持ってきた。
(実験用プログラムについてはまた別途考えてまとめようと思う。)
その結果

が得られた。
パラメータ同定するために入力を単位ステップ入力にしたいので回転数も印可電圧で割った結果を使用する。

パラメータ同定

パラメータ同定を行う際に最適モデルと実験値の最小二乗法で導出するが基本的にはMATLAB等のシュミレーションソフトを使用する。
しかしMATLABを持っていないので類似した機能を持つscilabというソフトを使用する。
scilabにはSciNotesでプログラムを書いて実行することでシュミレーションできる。
実験値から最適値を導出するプログラムを下記に示す。

//実験結果の読み込み
Vm=8;                                       //印可下ステップ電圧
Num=read('Vm_8V.txt', -1,2);                //テキストファイルの2列を全部読み取る
y_step=Num(:,2);                            //ステップ入力したときの応答値[rpm]
om_unit=(y_step*2*%pi)/60/Vm;               //単位ステップ入力したときの応答値[rad/s]
preTime=Num(:,1);                           //マイコン内時間[μs]
Time=(preTime-preTime(1,:))/1000000;        //応答時間[s]

//一時遅れ系の最小二乗法
function er=error_ff(x,Time,ym);
    P=syslin('c',x(2)/(%s+x(1)));           //伝達関数の定義
    y=csim(ones(Time'),Time,P);             //シュミレーション：出力[1×要素数]=csim(制御入力[1×要素数],時間[要素数×1],伝達関数)
    er=(y-ym)*(y-ym)';                      //誤差の2乗面積
endfunction

//最適値の導出
x0=[7,70];                                  //パラメータの初期値(初めは適当)
[f,xopt,gropt]=leastsq(list(error_ff,Time,om_unit'),x0);    //最小二乗法を使用した最適探索
disp(xopt)
a=xopt(1);b=xopt(2);
Popt=syslin('c',xopt(2)/(%s+xopt(1)));      //得られた最適値での伝達関数
yopt=csim(ones(Time'),Time,Popt);           //最適値での出力
plot2d(Time',yopt,2)                        //最適値のグラフ
plot2d(Time',om_unit,5)                     //実験結果のグラフ
xgrid();
title('ステップ入力応答'); //タイトル
xlabel('時間[s]'); //X軸ラベル
ylabel('角速度[rad/s]'); //Y軸ラベル

//パラメータの導出
kt=0.0362;                                  //トルク定数[N・m/A]
ke=0.0362;                                  //逆起電力定数[V・sec/rad]
R=13.72;                                    //抵抗[Ω]
Jm=(kt/R)/b;
um=a*Jm-(kt*ke)/R;

disp("a = "+string(a))
disp("b = "+string(b))
disp("Jm = "+string(Jm)+"[kg・m^2]")
disp("um = "+string(um)+"[-]")

プログラムの解説
1～7行目：実験結果のテキストファイルを読み込んで時間と実験値を行列に格納する。
9～14行目：シュミレーションと実験値の差の二乗和を計算する関数。
18行目：最小二乗法を使用した最適探索。xoptは伝達関数内の最適なa,bのこと。leastsqについては下記を参照にして行と列を合わせることが必要。
leastsq(list(error_ff,Time,om_unit'),x0)の中のリストではlist(error_ff,Time[41×1],om_unit'[1×41])となる。
phys.sci.hokudai.ac.jp

シュミレーション結果

伝達関数の未知パラメータは
$a=114.60078　,　b=3109.0526$
となった。
シュミレーションの応答値と実験値の応答値の比較は下図の通りであり、ほぼ同等であり最適値として問題ないと判断する。

よって物性値の未知パラメータ
$J_m=\displaystyle{ \frac{\frac{k_t}{R}}{b}}=8.4865\times10^{-7}　[kg\cdot{m}^2$ ]
$\mu_m=\displaystyle{a\cdot{J_m}-\frac{k_tk_e}{R}=1.7423\times10^{-6}}　[-$ ]
となる。

次回

慣性ロータをつけた状態で同様の試験を行い慣性ロータの慣性モーメントと粘性摩擦抵抗を求める。