Pythonで考えるDCモータの制御（２）モータのＰＩ制御について

はじめに

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

とりあえず分割することにしました

Pythonで考えるDCモータの制御が漠然とだらだらと思いつくままに、書き続けたためか、自分の文章力がないためか長くなってきたので、分割することにしました。でも、一記事にまとめるのも未だに未練があるのです。全部書き終わったら、様子を見て一記事にまとめるかもしれません。

ブログは日記だとすると、考えたり、活動したりするたびに分割して書いたほうが良いのでしょうが後で自分で見るメモと言う考えもあるので、その場合は一つにまとめておいて方が良いような気がします。

~~しばらくは、（１）の記事は、以前のままにして、その後半を本記事にコピーして加筆しています。~~ （１）の記事の本記事と重複している部分は完全に削除しました。（2020/5/11）

DCモータのＰＩ制御について

前回の記事で、ＤＣモータのモデルを伝達関数で表す際に、モータの微分方程式をそのまま伝達関数にした２次遅れモデルとインダクタンスを無視した、１次遅れモデルを作って、モータのステップ応答を各モデルで計算した場合にどうなるのか比較してみました。

kouhei-no-homepage.hatenablog.com

比較の結果はマクロ的に見れば違いはなく、ミクロ的にみるとやはり違いはあるものの、「速度制御の観点ならば簡単なモデルで設計等を進めても良いのでは？」というところまで来ました。

ただし、あくまでモデルの形による違いで、実際のモータについての議論はしていませんが、これについては、以前実施のモータの実験結果とと数学モデルから計算した結果の比較について書いて記事にありますように、おおむね問題はないかと考えています。

kouhei-no-homepage.hatenablog.com

今回は、モータの速度制御をＰＩ制御で行った場合を、モデルを用いて考えていきたいと思います。１次遅れモデルを使うとステップ上の目標ならば、閉ループの計算が解析的に可能でゲインを決める際も見通しの良い設計ができてよく、その結果を２次のモデルに適用しても、同じ結果になることを話したいと思います。

速度のPI制御システムの伝達関数

f:id:kouhei_ito:20200125190521j:plain — DCモータの速度のPI制御

DCモータの速度をフィードバックしてPI制御器で制御しようと思います。

マイクロマウス等のロボットを制御していると加速度そのものを制御したくなります。加速度は力であり、力（トルク）とは、すなわち電流なので電流を制御するのが良いのかなと妄想し始めます。既に妄想し始めて１７年以上・・・いまだに実現していません。

この記事では、まずは速度フィードバックのみの制御を考えて、余裕があったらインナーループとして電流制御を入れてみようと思っています。

閉ループの伝達関数

モータの伝達関数を $G_m$ 、PIコントローラの伝達関数を $C_{PI}$ とすると、PI制御の目標から速度までの伝達関数は

${\displaystyle \frac{C_{PI} G_m}{ 1+ C_{PI} G_m} }$

と書けます。

PIコントローラの伝達関数は比例ゲイン $K_P$ 、積分ゲイン $K_I$ とすると

${\displaystyle C_{PI}(s)=K_P + \frac{K_I}{s} = \frac{K_P s + K_I}{s} }$

これらから制御系の伝達関数を求めていきますが、前に倣ってモータモデルを２次遅れモデルの場合と１次遅れモデルの場合に分けて記述していきます。

モータモデルを２次遅れモデルとしたときの速度PI制御系の伝達関数

制御系（閉ループ）の伝達関数を $G_c(s)$ とします。

２次遅れモデルで考えた閉ループ伝達関数は以下になります。

${\displaystyle G_c(s)=\frac{KK_Ps+KK_I}{JL s^3 +(DL+JR)s^2+(DR+K^2+KK_P)s+KK_I}\\ }$

３次の伝達関数となり、もはや手計算では解くことは困難です。

モータモデルを１次遅れモデルとしたときの速度PI制御系の伝達関数

1次遅れモデルで考えた閉ループ伝達関数は以下になります。

${\displaystyle G_c(s)=\frac{KK_Ps+KK_I}{JR s^2 +(DR+K^2+KK_P)s+KK_I}\\ }$

モータモデルを１次遅れとすることで、閉ループ伝達関数は２次の伝達関数なので頑張れば手計算でも計算できそうです♪

１次遅れモータモデルによるPI速度制御系の設計について

PID制御のゲインを決めるのは限界感度法などこれまでいろいろと提案されてきました。今回のようにPI制御ならばパラメータが二つになり、どちらか一方を決めると片方を決める目安が出せそうだというのが、本記事でお話しすることです。

モータの１次モデルから導出した閉ループ伝達関数は２次の式です。その応答は伝達関数の極で決まるというのはこれまでの触れたことでした。

極の姿

２次の伝達関数の場合、極の姿には以下の三つがあります。

異なる二つの実数解
共役複素解
重解

これらは式の上でどうなるのかを、これから見ていきます。今述べた、３つの場合と言うのは２次方程式の判別式の話なので、閉ループ伝達関数の特性方程式の判別式を導きましょう。まず、特性方程式は以下のようになりますね。

${\displaystyle JR s^2 +(DR+K^2+KK_P)s+KK_I=0\\ }$

すると、判別式は $D_m$ は次のようになるはずです。

${\displaystyle D_m=(DR+K^2+KK_P)^2-4JRKK_I\\ }$

$K_P$ を固定して $K_I$ を求める

判別式の正負を見ていきますが、ここでは $K_P$ を決定した場合、 $K_I$ がどうなるのかを知りたいため、最終的に $K_I$ について、式を解いていきます。

異なる二つの実数解の場合

判別式が正の場合ですので以下のようになります。

${\displaystyle (DR+K^2+KK_P)^2-4JRKK_I \gt 0 \\ }$

これより、以下のようになり

${\displaystyle K_I \lt \frac{(DR+K^2+KK_P)^2}{4JRK} \\ }$

以上の範囲で $K_I$ を決めると、極が異なる実数解になります。

極が異なる二つの実数解の場合は、応答は振動的にはなりません。目標値に近づくのにやや時間がかかるような応答になります。

共役複素解の場合

この場合は先ほどの逆になります。

${\displaystyle K_I \gt \frac{(DR+K^2+KK_P)^2}{4JRK} \\ }$

極が共役複素解の場合は、応答は振動的にはなります。立ち上がりの早い応答になりますが、オーバーシュートします。

重解の場合

以上からあと残りのパターンはイコールになる場合なので、以下のようになります。

${\displaystyle K_I = \frac{(DR+K^2+KK_P)^2}{4JRK} \\ }$

極が重解の場合は全く振動しません。後で書くかもしれませんが、この場合は減衰係数が０の場合に相当します。結局のところ、この値を境に増やすか、減らすかでゲインを決めていきます。

機械的時定数とモータゲインで記述しなおし

すみません、話の流れで後になりましたが、１次遅れのモータモデルを機械的時定数とモータゲインで記述することもできます。これらの機械的時定数とモータゲインはたいてい仕様書に出てるので、こちらを使うほうがスムーズに事が運ぶかもしれません。１次遅れモデルは次のように書きます。

${\displaystyle G_m(s) = \frac{K_m}{\tau_m s +1} \\ }$

これを使った場合、今までの $K_I$ の範囲を決める式は、次のように書き換えることができます。結果だけ書きます。

異なる二つの実数解の場合

${\displaystyle K_I \lt \frac{(1+K_mK_P)^2}{4\tau_m K_m} \\ }$

共役複素解の場合

${\displaystyle K_I \gt \frac{(1+K_mK_P)^2}{4\tau_m K_m} \\ }$

重解の場合

${\displaystyle K_I = \frac{(1+K_mK_P)^2}{4\tau_m K_m} \\ }$

比例ゲインについて

比例ゲインを決めておけば積分ゲインの目安がわかると言うのが、これまでの内容でした。

であれあば比例ゲインはどう決めれば良いのか。当記事の手法やMatlabなどでシミュレーションして応答を見て思考錯誤します。ただし、出力の応答がよければそれで良いかと言うと、モータであれば制御入力の電圧が定格を超えていないかが気にするところです。

比例ゲインを大きくすると、初動の時に大きな制御入力が入ります。実機だったら焼けることもあるかもしれません。それを考えて比例ゲインを決めないとならないので割と保守的になりまし、それが目安になると考えます。

コントローラの中のプログラムで制御入力が定格を超えないようにリミッタをかける方法もあります。この場合は、初期の頃は「最大定格で駆動！」のようにフィードバックが切られた素の特性が剥き出しになります。

本記事でも、リミッタかけた場合も比較して良し悪しを触れられればと考えています。

速度のPI制御系のステップ応答

それでは、比例ゲイン、積分ゲインが決められたあとの速度PI制御のステップ応答を調べてみます。

モータモデルが２次遅れモデルの場合

速度PI制御の閉ループ伝達関数から。ステップ応答は以下の式で表す事ができます。

${\displaystyle \omega(s)=\frac{KK_Ps+KK_I}{JL s^3 +(DL+JR)s^2+(DR+K^2+KK_P)s+KK_I}\frac{\omega_r}{s}\\ }$

ここで $\omega_r$ は目標値です。前と同様に、伝達関数の極を $\alpha$ 、 $\beta$ 、 $\gamma$ とすると上式は以下のように表す事ができます。

${\displaystyle \omega(s)=\frac{(KK_Ps+KK_I) \omega_r}{JL(s-\alpha)(s-\beta)(s-\gamma)s}\\ }$

部分分数展開すると、以下のようになります。

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

係数をもとめます。

${\displaystyle K_1=\lim_{s \to \alpha}\frac{(KK_Ps+KK_I) \omega_r}{JL(s-\beta)(s-\gamma)s}=\frac{(KK_P\alpha+KK_I) \omega_r}{JL(\alpha-\beta)(\alpha-\gamma)\alpha}\\ }$

${\displaystyle K_2=\lim_{s \to \beta}\frac{(KK_Ps+KK_I) \omega_r}{JL(s-\alpha)(s-\gamma)s}=\frac{(KK_P\beta+KK_I) \omega_r}{JL(\beta-\alpha)(\beta-\gamma)\beta}\\ }$

${\displaystyle K_3=\lim_{s \to \gamma}\frac{(KK_Ps+KK_I) \omega_r}{JL(s-\alpha)(s-\beta)s}=\frac{(KK_P\gamma+KK_I) \omega_r}{JL(\gamma-\alpha)(\gamma-\beta)\gamma}\\ }$

${\displaystyle K_4=\lim_{s \to 0}\frac{(KK_Ps+KK_I) \omega_r}{JL(s-\alpha)(s-\beta)(s-\gamma)}=\omega_r\\ }$

以上から、逆ラプラス変換すると、以下になります。

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

これで２次モデルにおける速度PI制御におけるステップ応答が求まりました。

モータモデルが1次遅れモデルの場合

１次モデルにおける速度PI制御の閉ループ伝達関数から。ステップ応答は以下の式で表す事ができます。

${\displaystyle \omega(s)=\frac{KK_Ps+KK_I}{JR s^2 +(DR+K^2+KK_P)s+KK_I} \frac{\omega_r}{s}\\ }$

ここでも前と同様に、伝達関数の極を $\alpha$ 、 $\beta$ とすると上式は以下のように表す事ができます。

${\displaystyle \omega(s)=\frac{(KK_Ps+KK_I) \omega_r}{JR(s-\alpha)(s-\beta)s}\\ }$

部分分数展開すると、以下のようになります。

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

係数をもとめます。

${\displaystyle K_1=\lim_{s \to \alpha}\frac{(KK_Ps+KK_I) \omega_r}{JR(s-\beta)s}=\frac{(KK_P\alpha+KK_I) \omega_r}{JR(\alpha-\beta)\alpha}\\ }$

${\displaystyle K_2=\lim_{s \to \beta}\frac{(KK_Ps+KK_I) \omega_r}{JR(s-\alpha)s}=\frac{(KK_P\beta+KK_I) \omega_r}{JR(\beta-\alpha)\beta}\\ }$

${\displaystyle K_3=\lim_{s \to 0}\frac{(KK_Ps+KK_I) \omega_r}{JR(s-\alpha)(s-\beta)}=\omega_r\\ }$

以上から、逆ラプラス変換すると、以下になります。

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

これで１次モデルにおける速度PI制御におけるステップ応答が求まりました。

余談

何度も同じことをやっています。だからこそ計算機の力を借りる方が楽だと言う思いに行き着きます。僕らはあとでPython先生の力を借りていとも簡単に極を求めて、上の方程式を計算すれば応答を求められます。しかも、今はその式を展開してこの記事のようにコーディングしなくても良い世界で素晴らしいです。

しかし、古典制御理論を作った人たちはそれができませんでした。敬服します。

速度PI制御の極を計算

さて、DCモータの素のステップ応答を求めた際と同様に、速度PI制御系の極を求めてみましょう。

#速度PI制御システムの極の計算
import numpy as np

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

#1次遅れモデルのパラメータ
K_m=K/(D*R+K**2)
tau_m=J*R/(D*R+K**2)

kp=6/500

#比例ゲイン基準の積分ゲインの境界値
print('Ki boundary={}'.format((1+K_m*kp)**2/4/tau_m/K_m))

ki=1.7

#１次遅れモデル係数
A1=J*R
B1=D*R+K**2+K*kp
C1=K*ki

#2次遅れモデル係数
A2=J*L
B2=D*L+J*R
C2=D*R+K**2+K*kp
D2=K*ki

#１次遅れモデル極計算
p1=np.roots([A1, B1, C1])

#２次遅れモデル極計算
p2=np.roots([A2, B2, C2, D2])

#極の表示
print('Closed Loop Pole of 1st Order Model:{}'.format(p1))
print('Closed Loop Pole of 2nd Order Model:{}'.format(p2))

#Ki boundary=1.6827052018576538
#Closed Loop Pole of 1st Order Model:[-180.33064516+18.28198156j -180.33064516-18.28198156j]
#Closed Loop Pole of 2nd Order Model:[-45105.24969617 +0.j     -181.40848525+14.40762598j    -181.40848525-14.40762598j]

以上は、ある一つの比例ゲインと積分ゲインの組み合わせに行ったものです。比例ゲインが6/500に対して、積分ゲインの境界は1.68程度という事がわかりました。

積分ゲインの大きさと極

積分ゲインの大きさを境界値より少し大きい1.7とすると極の計算結果から、１次モデル、２次モデルともに近い値の極が二つあり、２次モデルには速いモードを示す極がもう一つある事が判ります。お互いに近しい極は共役複素解である事がわかり、応答が振動的になる事が予想できます。

速度PI制御のステップ応答

続いて、PI制御のステップ応答を計算してみます。１次モデルで決定した比例ゲイン、積分ゲインに対しての２次モデルの応答も合わせてみます。

#速度PI制御システムのステップ応答

%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt

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

#1次遅れモデルのパラメータ
K_m=K/(D*R+K**2)
tau_m=J*R/(D*R+K**2)

K_m=K/(D*R+K**2)
tau_m=J*R/(D*R+K**2)

#比例ゲイン
kp=6/500

#比例ゲイン基準の積分ゲインの境界値
print('Ki boundary={}'.format( (D*R+K**2+K*kp)**2/4/J/R/K ))
#Ki boundary=1.6827052018576538

#積分ゲイン
ki=1.5

#目標値
ref=150

#１次遅れモデル係数
A1=J*R
B1=D*R+K**2+K*kp
C1=K*ki

#2次遅れモデル係数
A2=J*L
B2=D*L+J*R
C2=D*R+K**2+K*kp
D2=K*ki

#１次遅れモデル極計算
p1=np.roots([A1, B1, C1])

#2次遅れモデル極計算
p2=np.roots([A2, B2, C2, D2])

#極の表示
print('Closed Loop Pole of 1st Order Model:{}'.format(p1))
print('Closed Loop Pole of 2nd Order Model:{}'.format(p2))

k11=ref*(K*kp*p1[0]+K*ki)/J/R/(p1[0]-p1[1])/p1[0]
k12=ref*(K*kp*p1[1]+K*ki)/J/R/(p1[1]-p1[0])/p1[1]
k13=ref

k21=ref*(K*kp*p2[0]+K*ki)/J/L/(p2[0]-p2[1])/(p2[0]-p2[2])/p2[0]
k22=ref*(K*kp*p2[1]+K*ki)/J/L/(p2[1]-p2[0])/(p2[1]-p2[2])/p2[1]
k23=ref*(K*kp*p2[2]+K*ki)/J/L/(p2[2]-p2[0])/(p2[2]-p2[1])/p2[2]
k24=ref

#係数の表示
print('Closed Loop Coef of 1st Order Model:\n{}'.format([k11, k12, k13]))
print('Closed Loop Coef of 2nd Order Model:\n{}'.format([k21, k22, k23, k24]))

#応答計算
t=np.linspace(0,0.05,100000)
y1=k11*np.exp(p1[0]*t)+k12*np.exp(p1[1]*t)+k13
y2=k21*np.exp(p2[0]*t)+k22*np.exp(p2[1]*t)+k23*np.exp(p2[2]*t)+k24

plt.figure(figsize=(8,4))
plt.plot(t, y1, linewidth=4, label='1stOrderModel')
plt.plot(t, y2, label='2ndOrderModel')
plt.legend(loc='lower right')
plt.xlabel('Time(s)')
plt.ylabel('Anglurr Velocity(rad/s)')
plt.grid()
plt.show()

#計算結果
#kp=6/500
#ki=1.5
#Ki boundary=1.6827052018576538
#Closed Loop Pole of 1st Order Model:[-239.75178441 -120.90950591]
#Closed Loop Pole of 2nd Order Model:[-45105.16261945   -242.3144503    -120.58959692]
#Closed Loop Coef of 1st Order Model:
#[-140.09743062266793, -9.9025693773321, 150]
#Closed Loop Coef of 2nd Order Model:
#[0.7815215899368377, -140.2176714096278, -10.56385018030904, 150]

比例ゲインを6/500にすると、積分ゲインの境界値は1.68程度です。今回のプログラムでは重解の場合を計算しようとすると係数が発散します。重解の際の計算式をちゃんと用意すれば良いのですが、手抜きです。この記事がこんなに長くなるとは思わず、最初は書くつもりでしたがまたの宿題にさせてください。

Pythonでの複素数計算

極が共役複素根の時は、ステップ応答式はオイラーの公式を適用して三角関数が出てくる式を用意しないとダメなんですが、 Pythonは警告は出るのですが、複素数の計算をちゃんとやってくれてそれらしい値を出してくれます。ビバ！Python先生（警告が出て気になりますが、こんかいは追求しないことにします。これらの計算は必ず虚部が打ち消されて計算結果が実数になります。この辺りの計算も楽しいのですが、今回は省略。）