python-control を使って基本的な制御系解析を行う方法について触れます。
古典制御の復習をしようと思って「演習で学ぶ基礎制御工学」という本を手にとって勉強していました。この参考書の問題を解くときに、ボード線図やステップ応答をグラフで描画してみたいなと思ったので、やり方をまとめます。
python を使うにあたって、たとえば sympy がないよといったエラーメッセージが出たときには適宜 pip3 でインストールするようにしてください。
ModuleNotFoundError: No module named 'sympy' $ pip3 install sympy
今回扱うのは、ボード線図・極配置・時間応答・根軌跡の描き方です。
以下の目次は、冒頭紹介した「演習で学ぶ基礎制御工学」内の問題番号と合わせています。
演習 7.17 2次遅れ要素の時間応答
問題
$$G(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_ns + \omega_n^2}$$
$\zeta = 0.5, \omega_n = 0.1, 0.5, 1$ としたときの極配置・ボード線図・ステップ応答を求めよ。
プロット
極・ゼロ配置
ボード線図
ステップ応答
ソースコード
from scipy import signal
import matplotlib.pyplot as plt
import control
from control.matlab import *
import numpy as np
# Definition of transfer function's numerator and denominator
zeta = 0.5
wn1 = 0.1
num1 = [wn1**2]
den1 = [1, 2*zeta*wn1, wn1**2]
wn2 = 0.5
num2 = [wn2**2]
den2 = [1, 2*zeta*wn2, wn2**2]
wn3 = 1.0
num3 = [wn3**2]
den3 = [1, 2*zeta*wn3, wn3**2]
# Plot pzmap
G1 = control.tf(num1, den1)
G2 = control.tf(num2, den2)
G3 = control.tf(num3, den3)
#control.pzmap(G1)
#control.pzmap(G2)
#control.pzmap(G3)
p1 = control.pole(G1)
z1 = control.zero(G1)
p2 = control.pole(G2)
z2 = control.zero(G2)
p3 = control.pole(G3)
z3 = control.zero(G3)
plt.scatter(p1.real, p1.imag, marker='x', label='G1 pole')
plt.scatter(p2.real, p2.imag, marker='x', label='G2 pole')
plt.scatter(p3.real, p3.imag, marker='x', label='G3 pole')
plt.scatter(z1.real, z1.imag, marker='o')
plt.scatter(z2.real, z2.imag, marker='o')
plt.scatter(z3.real, z3.imag, marker='o')
plt.axis([-1, 1, -1, 1])
plt.gca().set_aspect("equal")
plt.xlabel("Re")
plt.ylabel("Im")
plt.title("Pole Zero Map")
plt.legend()
plt.grid()
plt.show()
# Draw bode diagram
G1 = signal.lti(num1, den1)
G2 = signal.lti(num2, den2)
G3 = signal.lti(num3, den3)
warr = np.linspace(0.001, 10**3, 10**5)
w, mag1, phase1 = signal.bode(G1, warr)
w, mag2, phase2 = signal.bode(G2, warr)
w, mag3, phase3 = signal.bode(G3, warr)
plt.subplot(2, 1, 1) # Gain diagram
plt.semilogx(w, mag1, label='G1')
plt.semilogx(w, mag2, label='G2')
plt.semilogx(w, mag3, label='G3')
plt.ylabel("Magnitude [dB]")
plt.legend()
plt.grid()
plt.subplot(2, 1, 2) # Phase diagram
plt.semilogx(w, phase1, label='G1')
plt.semilogx(w, phase2, label='G2')
plt.semilogx(w, phase3, label='G3')
plt.xlabel("w [rad/sec]")
plt.ylabel("Phase [deg]")
plt.legend()
plt.grid()
plt.show()
# Plot step response
G1 = signal.TransferFunction(num1, den1)
G2 = signal.TransferFunction(num2, den2)
G3 = signal.TransferFunction(num3, den3)
t = np.linspace(0, 100, 1000)
u = np.ones_like(t) # ramp input case: u = 1*t
tout, resp1, x = signal.lsim(G1, u, t)
tout, resp2, x = signal.lsim(G2, u, t)
tout, resp3, x = signal.lsim(G3, u, t)
plt.plot(t, u, label='Input')
plt.plot(t, resp1, label='G1')
plt.plot(t, resp2, label='G2')
plt.plot(t, resp3, label='G3')
plt.legend()
plt.title('Step Response')
plt.xlabel('Time [s]')
plt.ylabel('Response')
plt.grid()
plt.show()極・ゼロ配置
はじめ以下の方法でプロットしようとしたのですが、同じグラフに複数の極・ゼロ配置をプロットできませんでした。
control.pzmap(G1)そこで、control.pole() 関数と control.zero() 関数で極・ゼロを抽出して複素平面上に散布図でプロットするようにしました。
p1 = control.pole(G1)
z1 = control.zero(G1)
plt.scatter(p1.real, p1.imag, marker='x', label='G1 pole')
plt.scatter(z1.real, z1.imag, marker='o')ボード線図
公式のマニュアルを参照して、warr でプロットの範囲を指定するようにしました。今回の場合だと、$\left[ 10^{-3}, 10^3 \right]$ です。
G1 = signal.lti(num1, den1)
warr = np.linspace(0.001, 10**3, 10**5)
w, mag1, phase1 = signal.bode(G1, warr)※上記 warr を設定しないと、$\omega_n = 0.1$ の場合のボード線図が正しく表示されませんでした。python-control のバグなのかもしれません。
ステップ応答
時間の配列 t とステップ入力の配列 u を定義して、signal.lsim() という関数で時間応答 resp を出力しています。
G1 = signal.TransferFunction(num1, den1)
t = np.linspace(0, 100, 1000)
u = np.ones_like(t)
tout, resp1, x = signal.lsim(G1, u, t)思ったこと
| 極・ゼロ配置 | G1 = control.tf(num1, den1) |
| ボード線図 | G1 = signal.lti(num1, den1) |
| ステップ応答 | G1 = signal.TransferFunction(num1, den1) |
上記のようにそれぞれ別の方法で伝達関数を定義してしまいました。もしかしたら、記述の仕方を統一できるかもしれませんね…。ご存知の方いらしたら教えてください。
演習 9.24 根軌跡の適用
問題
一巡伝達関数が以下で与えられる制御系の根軌跡を描け。
$$ G(s)H(s) = \frac{K}{s^4 + 12s^3 + 64s^2 +128s}$$
プロット
ソースコード
from control.matlab import *
from matplotlib import pyplot as plt
def main():
num = [1]
den = [1, 12, 64, 128, 0]
sys = tf(num, den)
rlocus(sys)
# plt.axis([-15, 15, -15, 15])
plt.gca().set_aspect("equal")
plt.show()
if __name__ == "__main__":
main()結構複雑な軌跡ですが、このサイトに従ってパラメータを変えてみたら案外簡単に描けました。
まとめ
学生時代によく matlab を使ってグラフをプロットしていましたが、意外と python-contol でもグラフ描画できるということを学びました。
コメント