Pythonで考えるDCモータの制御（１）モデルによるステップ応答の比較

モデルによるステップ応答の比較アイキャッチ画像

はじめに

こんにちは、こうへいです

初めてチャレンジしたロボットランサーから、 DCモータの制御をこれまで取り扱う事が多かったんですが同じ事を何度も繰り返しているので、一度整理して見たくなりました。これまでの記事と確実に重複しますが、それを承知でまとめてみます。

こうして整理してみると新たな何かも見えてくるかもしれません。制御のことを考えて行くと次数が高くなると、計算機の力に頼らないといけなくなります。そこで、今回はPythonで計算して確認してみようと思います。

プログラム言語はなんでもいいのです、制御を扱うのでMatlabでもいいですし Cで書いてもいいです。Pythonならば比較的どのプラットフォームでも使えて割と簡単にやりたい事がかけるので、制御工学の検証には便利かなと思っています。

DCモータのモデルを１次モデルで近似する事が多いのですが１次モデルと２次モデルの比較をしてみたいと思います

１次モデルを用いてPI制御器を設計すると、２次方程式が解けると応答の計算ができて比例ゲインと積分ゲインの関係式も手でかけたりして割と見通しが良いかなと思っています。

２次モデルだとPI制御器を設計するにも３次方程式がでてきて、解の公式があるものの覚えていないのですっきりしません。（ただし、３次までなら解の公式があるので、人間が覚えてやるにはしんどいですが、プログラムにしてしまえば何度でも大丈夫なので、この際、Pythonの力も借りながら、検討してみようかと思い始めています。記事の後半に追加するかもしれません。）

１次モデルで設計して２次モデルに適用してどうなるのかも見てみたいです。

そしてモータのステップ応答や、PI制御を施した際の応答、さらにはランプ入力に対する応答、マイクロマウスなどのロボット制御によく用いられる台形加減速制御について考えていきたいと思います。

はじめに
想定するモータ
DCモータの数学モデル
DCモータの伝達関数
- 2次遅れモデル
- 1次遅れモデル
DCモータのステップ応答
モータのステップ応答２次遅れモデルと１次遅れモデル比較
- 1724T006SRモータの機械的時定数と電気的時定数
- Pythonでステップ応答計算
DCモータの極について
- 極の判別式
おわりに

想定するモータ

自分がエンコーダ付きを持っているので。 Faulhaber社の1724T006SRモータを想定します。

f:id:kouhei_ito:20200113143437p:plain — 1724T006SRモータの仕様

各パラメータを単位をSIに揃えて、次のようにしました。

R=3.41
K=6.59e-3
L=75e-6
D=1.4e-7
J=1e-7

動粘性係数のDは仕様には載っていません。後ほど推定方法も明らかにします。

DCモータの数学モデル

以前に同じことを記事にしました。

rikei-tawamure.com

f:id:kouhei_ito:20200105205719j:plain — DCモータのモデル

DCモータの数学モデルについて再掲載します。

回路の方程式

${\displaystyle L \frac{d i}{d t} + R i +K \omega = e }$

回転に関する運動方程式

${\displaystyle J \frac{d \omega}{d t} + D \omega +T_L = K i }$

各変数は以下の通りです

${i}$ : 電流(A)
${e}$ : 入力電圧 (V)
${\omega}$ : 軸の角速度（rad/s）
${T_L}$ : 負荷トルク（Nm）

各パラメータは以下の通りです　

${L}$ : インダクタンス(H)
${R}$ : 巻線抵抗(Ω)
${K}$ : 逆起電力係数（Vs/rad）、トルク定数（N/A）（同じ値）
${D}$ : 動粘性係数（Nms/rad）
${J}$ : 慣性モーメント（Kg m²）

DCモータの伝達関数

2次遅れモデル

上記の数学モデルからDCモータの伝達関数を求めると次のようになります。

${\displaystyle \omega(s) =\frac{K}{JL s^2 + (DL+JR)s + DR+K^2} e(s)\\ }$

1次遅れモデル

一般的にL=0としてDCモータの伝達関数を1次遅れで近似することがあります。近似の伝達関数は以下のようになります。

${\displaystyle \omega(s) =\frac{K}{JRs + DR+K^2} e(s)\\ }$

DCモータのステップ応答

上記の伝達関数からDCモータのステップ応答を求めてみます。この際に入力については以下のようにします。

${\displaystyle e(s)=\frac{e}{s}\\ }$

2次遅れモデルの場合

応答の式は以下のようになります。

${\displaystyle \omega(s) =\frac{K}{JL s^2 + (DL+JR)s + DR+K^2} \frac{e}{s}\\ }$

この式を逆ラプラス変換すると時間の関数 $\omega(t)$ が求まります。逆ラプラス変換するためには制御工学ではヘビサイド展開定理と言うものを使って求められます。この展開定理をすると部分分数展開して伝達関数を有利関数の和の形にする事ができます。

有利関数の和にできると簡単に逆ラプラス変換することができます。その結果は、指数関数の和の形となって応答の時間関数が求められます。ここのところの流れが僕は好きで「工学楽しいよね♪」と思うところです。

伝達関数の分母＝０の式を特性方程式と呼びますが、特性方程式の解を極と呼び、制御工学ではシステムの挙動を示す重要なものです。

ところで、極は特性方程式の解なので因数定理から極の値を使うと特性方程式は因数分解できます。以下についてはそれを使って式を変形しています。

2次遅れモデルにおいてはDCモータの伝達関数の二つの極を $\alpha$ と $\beta$ としたとき上式は因数定理により以下のように書き換えられます。

${\displaystyle \omega(s) =\frac{Ke}{JL (s-\alpha)(s-\beta)s}\\ }$

ちなみに、安定なシステムは全ての極の実部が負です。もし $\alpha$ と $\beta$ が実数ならば $\alpha \lt 0$ と $\beta \lt 0$ です。

上式を部分分数展開すると

${\displaystyle \omega(s) =\frac{K_1}{s-\alpha} + \frac{K_2}{s-\beta} + \frac{K_3}{s} \\ }$

これらの係数 $K_x$ は

${\displaystyle K_1 =\lim_{s \to \alpha} \frac{Ke}{JL(s-\beta)s} = \frac{Ke}{JL(\alpha-\beta) \alpha}\\ }$

${\displaystyle K_2 =\lim_{s \to \beta} \frac{Ke}{JL(s-\alpha)s} = \frac{Ke}{JL(\beta-\alpha) \beta}\\ }$

${\displaystyle K_3 =\lim_{s \to \beta} \frac{Ke}{JL(s-\alpha)(s-\beta)} = \frac{Ke}{JL\alpha \beta}= \frac{Ke}{DR + K^2}\\ }$

最後の $K_3$ についてはチート技（と言って良いものか微妙ですが・・・）があって、この値が応答の定常値という事がわかっているので先ほど出した、もとの応答を表す式

${\displaystyle \omega(s) =\frac{K}{JL s^2 + (DL+JR)s + DR+K^2} \frac{e}{s}\\ }$

上式に次のように最終値の定理を施して

${\displaystyle \omega(\infty) = \lim_{s \to 0} s \frac{K}{JL s^2 + (DL+JR)s + DR+K^2} \frac{e}{s}= \frac{Ke}{DR + K^2}\\ }$

とすると $K_3$ が求まります。

これで、 $\omega(s)$ を逆ラプラス変換する準備ができまして、結果は次のようになります。

${\displaystyle \omega(t) = K_1 e^{\alpha t} + K_2 e^{\beta t} + K_3\\ }$

1次遅れモデルの場合

続いて、１次遅れモデルです。

${\displaystyle \omega(s) =\frac{Ke}{(JRs + DR+K^2)s} \\ }$

${\displaystyle \omega(s) =\frac{Ke}{JR\left(s + \frac{DR+K^2}{JR}\right)s} \\ }$

これを部分分数展開すると

${\displaystyle \omega(s) =\frac{K_1}{s + \frac{DR+K^2}{JR}} + \frac{K_2}{s} \\ }$

$K_x$ を求めると

${\displaystyle K_1 = \lim_{s \to -\frac{DR+K^2}{JR}} \frac{Ke}{JRs}=-\frac{Ke}{DR + K^2}\\ }$

${\displaystyle K_2= \lim_{s \to 0} \frac{Ke}{JR \left( s+ \frac{DR+K^2}{JR}\right) }=\frac{Ke}{DR + K^2}\\ }$

ということで逆ラプラス変換すると１次遅れモデルのステップ応答は以下です。

${\displaystyle \omega(t) = K_1 e^{- \frac{DR+K^2}{JR}} + K_2\\ }$

$K_1$ と $K_2$ に上で求めたものを代入してしまうと、以下のように応答式が求まります。

${\displaystyle \omega(t) = \frac{Ke}{DR+K^2} \left( 1 - e^{-\frac{DR+K^2}{JR} t} \right)\\ }$

この時間応答式からも明らかに、1次遅れモデルでもモータの定常角速度は

${\displaystyle \omega(\infty) = \frac{Ke}{DR+K^2}\\ }$

となります。これは２次遅れモデルのそれとも一致し、１次、２次のモデルとも、少なくとも定常値は一致することがわかります。

1次遅れのモデルでは

${\displaystyle K_m = \frac{Ke}{DR+K^2} \\ }$

これをモータのゲイン

${\displaystyle \tau_m = \frac{JR}{DR+K^2}\\ }$

これをモータの時定数として

${\displaystyle \omega(t) = K_m \left( 1 - e^{-\frac{1}{t_m} t} \right)\\ }$

と表記することもあります。

モータの定常角速度と動粘性係数 $D$

モータの定常角速度がわかりました。両モデルとも一致します。通常モータの仕様書に動粘性係数 $D$ は記述されていませんが。今検討で使用しているモータの数学モデルには含まれているので、推定しなければなりません。

定常角速度は仕様書に無負荷回転数が載っているので、そこから未知数 $D$ を求めます。定常角速度を $\omega_{\infty}$ とすると。

${\displaystyle D = \frac{K}{\omega_{\infty}} \frac{e-K\omega_{\infty}}{R}\\ }$

上式の右辺の後半の分母が $R$ の部分は定常状態の電流を計算する式です。

以上は、モータのモデル式を以下のように微分項を定常状態ならば０であるとして、連立して解いても同じ結果が得られます。

${\displaystyle R i +K \omega_\infty = e\\ D \omega_\infty +T_L = K i\\ }$

余談かもしれませんが、モータの仕様書には、無負荷時の電流値が載っています。上の式の後半が電流値なので、仕様書の電流値を用いても $D$ は求まりますが、電流値は有効桁数が１桁しかなく、すこし違う値になります。

モータのステップ応答２次遅れモデルと１次遅れモデル比較

Pythonで多項式方程式の解を求めるためにはnumpyを使って以下のようにします。

import numpy as np

#方程式　Ax^2 + B x + C =0 の根を求める。

np.roots([A, B, C])

２次遅れモデルの極を計算してみましょう。

import numpy as np

#モータの諸元設定
R=3.41
K=6.59e-3
L=75e-6
D=1.4e-7
J=1e-7

#２次方程式の係数設定
A=J*L
B=D*L+J*R
C=D*R+K**2

#極を求める
pole=np.roots([A, B, C])

#極の表示
print(pole)

#時定数の表示
print(-1/pole)

#[-45338.94883714   -129.11782952]
#[2.20560914e-05 7.74486377e-03]

1724T006SRモータの機械的時定数と電気的時定数

上記のプログラムで極は約-45339と-129と判りました。絶対値の大きい方が収束が早く、小さい方が収束が遅いモードになります。逆数の大きさは時定数と呼ばれ単位は時間で収束の時間的目安になります。速いモードの時定数は22usで遅いモードの時定数は7.7msです。こちらは仕様書の機械的時定数8msにほぼ一致しています。 ~~早い方の時定数は電気的時定数と呼ばれます。~~ すいません、早いほうのモードの時定数は名前はまだないというのが正解かもしれません。電気的時定数は業界ではL/Rとなっていて、今回のお話の中で出てきた絶対値が大きいほうの極の逆数では無いです。（さらに確認したいと思います。）ただL/Rを計算してもその値は有効桁２桁で表示すると22usとなりほぼ同じですので、ほぼ問題なしです。

Pythonでステップ応答計算

それでは、先ほど示したステップ応答の式を１次モデルと２次モデルそれぞれに対して計算して比較してみましょう。

import numpy as np
import matplotlib.pyplot as plt

#モータの諸元設定
R=3.41
K=6.59e-3
L=75e-6
D=1.4e-7
J=1e-7

#入力電圧
e=1

#特性方程式の係数
A=J*L
B=D*L+J*R
C=D*R+K**2

#極の算出
pole=np.roots([A, B, C])

#時間応答の指数関数の係数
k1=K*e/J/L/(pole[0]-pole[1])/pole[0]
k2=K*e/J/L/(pole[1]-pole[0])/pole[1]
k3=K*e/(D*R+K**2)

#１次遅れモデルのゲインと時定数
K_m=K/(D*R+K**2)
tau_m=J*R/(D*R+K**2)

#計算結果の表示
print('A, B, C=' ,A, B, C)
print('pole=', pole)
print('K_m, tau_m, 1/tau_m=',K_m, tau_m, -1/tau_m)
print('k1,k2,k3=', k1, k2, k3)

#ステップ応答の計算 y1：１次モデル　y2:２次モデル
t=np.linspace(0.00,0.1,1000)
y1=K_m*e*(1-np.exp(-t/tau_m))
y2=k1*np.exp(pole[0]*t) + k2*np.exp(pole[1]*t) + k3

plt.plot(t, y1 , linewidth=4)
plt.plot(t, y2)

plt.grid()
plt.show()

#A, B, C= 7.499999999999998e-12 3.410105e-07 4.390550000000001e-05
#pole= [-45338.94883714   -129.11782952]
#K_m, tau_m, 1/tau_m= 150.09509059229478 0.0077666807119836905 -128.7551319648094
#k1,k2,k3= 0.4286667719668311 -150.5237573642616 150.09509059229478