理系的な戯れ

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

Pythonで考えるDCモータの制御(2)モータのPI制御について

はじめに

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

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

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

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

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

DCモータのPI制御について

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

kouhei-no-homepage.hatenablog.com

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

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

kouhei-no-homepage.hatenablog.com

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

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

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

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

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

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

閉ループの伝達関数

モータの伝達関数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}
}

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

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

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

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

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

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

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

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

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

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

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

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

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

極の姿

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

  1. 異なる二つの実数解
  2. 共役複素解
  3. 重解

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

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

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

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

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

{\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制御のステップ応答を調べてみます。

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

速度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\\
}

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

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

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

これで1次モデルにおける速度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

#1次遅れモデル係数
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

#1次遅れモデル極計算
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))

#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とすると 極の計算結果から、1次モデル、2次モデルともに近い値の極が二つあり、2次モデルには速いモードを示す極がもう一つある事が判ります。 お互いに近しい極は共役複素解である事がわかり、応答が振動的になる事が予想できます。

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

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

#速度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

#1次遅れモデル係数
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

#1次遅れモデル極計算
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先生 (警告が出て気になりますが、こんかいは追求しないことにします。 これらの計算は必ず虚部が打ち消されて計算結果が実数になります。この辺りの計算も楽しいのですが、今回は省略。)

PI制御の結果:1次モデルの極が異なる実数になる場合 (リニア制御)

f:id:kouhei_ito:20200126174210p:plain
極が異なる実数解になる場合のPI制御結果比較

PI制御の結果:1次モデルの極が共役複素解になる場合 (リニア制御)

f:id:kouhei_ito:20200126174342p:plain
極が共役複素解になる場合のPI制御結果比較

これらの結果から、計算するモデルがどちらであっても、速いモードの影響はマクロ的にはみられません。 また、比例ゲインと積分ゲインをどんな値にしても両者の計算モデルはマクロ的には一致した結果が出ます。 簡易的なシミュレーションや設計においては1次モデルを使っても良いのではという感じがしてきます。

当然、拡大すると、下図のように誤差はあります。

f:id:kouhei_ito:20200126175054p:plain
拡大した両者の比較

これで、PI制御についてモデルの応答を確認する事ができました。

モータを簡易的な1次遅れモデルで近似して用いたとしても速度制御の計算結果は検討に使えると思います。また、2次方程式の判別式からゲインの目安も 解析的に求まるので見通しが良いのではと思います。

3次方程式の解の公式や判別式を適用するとわざわざ近似モデルを使わなくても良いという話になります。趣味としてはそこも確かめてみたいところですが、 2次の判別式から導いた簡単な計算で、振動的なるか、ならないかの目安が出せますので、そこまで追求する必要もないのかもしれません。

ここまで来ると負荷を考慮した話をしないとならないような気がしてきますが、とりあえず、この後はランプ入力に追随する制御の話に進みたいと思います。

おわりに

今回は、PI制御について考えました。

  • PI制御についても、1次モデルと2次モデルの速度の観点からの違いはほぼない
  • 簡単な1次モデルを用いた検討で見通し良くゲインが決められ、2次モデルに検討しても良好な結果が得られる

こんな結果が得られました。

以後は

  • ランプ入力に対する2型PI制御 リニア制御
  • 台形加減速制御
  • 電流連続モードのPWM制御
  • モータドライバに関する妄想

等々を細々と続けていきたいと思います。