理系的な戯れ

理工学系とくにロボットやドローンに関する計算・プログラミング等の話題を扱って、そのようなことに興味がある人たちのお役に立てればと思っております。

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

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

はじめに

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

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

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

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

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

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

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

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

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

想定するモータ

自分がエンコーダ付きを持っているので。 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モータの数学モデル

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

kouhei-no-homepage.hatenablog.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 m2

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)が求まります。 逆ラプラス変換するためには制御工学では ヘビサイド展開定理と言うものを使って求められます。この展開定理をすると部分分数展開して伝達関数を有利関数の和の形にする事ができます。

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

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

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

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次遅れモデルの場合

続いて、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}\\
}

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

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

K_1K_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}\\
}

となります。これは2次遅れモデルのそれとも一致し、1次、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の部分は定常状態の電流を計算する式です。

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

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

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

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

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

import numpy as np

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

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

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

import numpy as np

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

#2次方程式の係数設定
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を計算してもその値は有効桁2桁で表示すると22usとなりほぼ同じですので、 ほぼ問題なしです。

Pythonでステップ応答計算

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

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)

#1次遅れモデルのゲインと時定数
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:1次モデル y2:2次モデル
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

計算してみると1次遅れモデルの時定数も7.7msとなることが確認できました。 最終値は約150rad/sです。

ステップ応答のグラフを書いて比較すると

f:id:kouhei_ito:20200124231400p:plain
ステップ応答比較
両モデルとも完全に一致しているように見えます。初期時点の部分を拡大すると以下になります。

f:id:kouhei_ito:20200124231433p:plain
ステップ応答比較 初期部分拡大
初期の頃は電気的時定数が表す速いモードの影響がみられます。

こうしてみると、マクロ的には2次遅れモデルと1次遅れモデルと一致しています。ミクロ的には違いますが、 速度制御の観点からすると1次モデルで設計しても良いと感じます。 PWM制御などするときは電気的時定数を気にしなければなりませんので、次は制御を施した場合の両モデルの比較をしたいと思います。

DCモータの極について

DCモータの極を求めるのにPythonを使いましたが、所詮は2次方程式の解ですので、解析的な答えを書いておこうと思います。 次のようになると思います。

{\displaystyle
-\frac{DL+JR}{2JL} \pm \frac{1}{2JL}\sqrt{(DL+JR)^2 - 4JL(DR+K^2)}\\
}

上記の答えを見ますと、DCモータの各パラメータはすべて正の定数のはずですから、 仮に上記の式の値が共役複素数となったとしてもその実部はすべて負となって DCモータは安定です。電気繋げたら暴走するDCモータは基本的には無いということになりますかね。

極の判別式

極の式を見ると、後ろのルート内が気になります。判別式です。

{\displaystyle
D_m=(DL+JR)^2 - 4JL(DR+K^2)\\
}

これがマイナスになると極が共役複素数になり、もしかしたら勝手に振動する モータが作れるかもしれないんですが、どうなんだろう?

上式を展開して整理しなおすと

{\displaystyle
D_m=(DL-JR)^2 - 4JLK^2\\
}

前のカッコの中が0になれば確実に複素数になりますね。 でも、例えば今回例に取り上げている1724モータの場合で考えるとJRDLの107倍ぐらい大きく 難しいのかもしれません。

おわりに

モーターモデルを1次遅れモデルと2次遅れモデルで比較してみました。

マクロ的に見れば、違いは無く見ることができ、収束の早い極の無視できる 世界ならば、十分に1次遅れモデルでの制御系の設計をするのも ありだと考えます。

次回はDCモータのPI制御について考えていきます。