制御工学の計算に使える Python-Control のライブラリ基本機能を紹介します。
Python-Control は Python 環境で Matlab のように演算ができます。
ライブラリ情報
Python-Control ライブラリ内の関数は下記サイトに情報がまとめられています。
Matlab 風なので参考として Matlab の資料をみるのもよいです。
ただし Matlab と全く同じではないので注意してください。
伝達関数モデル作成
ライブラリ「control」から「matlab」をインポートします。
伝達関数（Transfer Function）のモデルは関数 tf([num], [den]) で作成できます。
num は分子、den は分母を入れます。
from control import matlab
Ga = matlab.tf([1], [1, 2])
Gb = matlab.tf([1], [1, 2, 3])
Gc = matlab.tf([1, 2], [1, -2, 3, -4])
print(Ga)
print(Gb)
print(Gc)実行結果
1
-----
s + 2
1
-------------
s^2 + 2 s + 3
s + 2
---------------------
s^3 - 2 s^2 + 3 s - 4
実用性を踏まえて２次系のモデルを作ります。
$$ G(s) = \frac{\omega_{n}^2}{s + 2 \zeta \omega_{n} s + \omega_{n}^2 } $$
パラメータを定義すれば扱いやすくなります。
zeta = 2
wn = 3
G = matlab.tf([wn*wn], [1, 2*zeta*wn, wn*wn])
print(G)実行結果
9
--------------
s^2 + 12 s + 9
伝達関数モデルの関数定義
形の決まった伝達関数モデルは関数で定義すると便利です。
以下は２次系のモデルを関数化した例です。
# 伝達関数モデル作成 (2次系)
def tf_2nd_order(k, zeta, wn):
num = [k*wn*wn]
den = [1, 2*zeta*wn, wn*wn]
return matlab.tf(num, den)
G = tf_2nd_order(1, 2, 3)
print(G)実行結果
9
--------------
s^2 + 12 s + 9
PID 補償器なども関数化すると便利です。
# 伝達関数モデル作成 (PID補償器)
def tf_controller_pid(kp, ki, kd):
# 比例補償器
num_p = [kp]
den_p = [1]
Gp = matlab.tf(num_p, den_p)
# 積分補償器
num_i = [ki]
den_i = [1, 0]
Gi = matlab.tf(num_i, den_i)
# 微分補償器
num_d = [kd, 0]
den_d = [1]
Gd = matlab.tf(num_d, den_d)
return Gp + Gi + Gd
ボード線図作成
ボード線図の作図は関数 bode( ) を使用します。
括弧の中には伝達関数を入れます。
図を出力するためにライブラリ「matplotlib」から「pyplot」をインポートします。
# Python上のコード
from control import matlab
from matplotlib import pyplot as plt
zeta = 2
wn = 3
G = matlab.tf([wn*wn], [1, 2*zeta*wn, wn*wn])
matlab.bode(G) # ボード線図のプロット
plt.show()実行結果
パラメータは適当ですが、2次系の伝達関数でゲインが -40dB となっており、位相は 0° から -180° まで変化しているので正しく出力できているのではないかと思います。
ゲイン余裕・位相余裕
計算には関数 margin(sys) または margin(mag, phase, w) を使います。
ゲイン余裕は位相が -180° のときにゲインが 0dB に対して何dBのマージンがあるのかを示したものです。
位相余裕はゲインが 0dB のときに位相が -180° に対して何度のマージンがあるのかを示したものです。
margin(sys) の引数
sys は単入力単出力（SISO）システムの伝達関数を入れます。
margin(mag, phase, w) の引数
ゲイン（magnitude）、位相（phase 単位はdeg）、周波数（frequencies 単位は rad/sec ）をボード線図の周波数応答データから入力します。
戻り値 Return
演算結果として次のパラメータが返ってきます。
gm：ゲイン余裕
pm：位相余裕 [deg]
wg：ゲイン余裕の周波数（位相 -180° の周波数）
wp：位相余裕の周波数（ゲイン 0dB の周波数）
今回は margin(sys) で使ってみます。
from control import matlab
from matplotlib import pyplot as plt
sys = matlab.tf([1], [1, 2, 1, 0])
print(sys)
gm, pm, wg, wp = matlab.margin(sys)
print(gm, pm, wg, wp)
matlab.bode(sys) # ボード線図のプロット
plt.show()結果
例を計算させるとゲイン余裕は 2.0dB、位相余裕は 21.39° です。
ボード線図を描くと合っていそうな感じです。
1
---------------
s^3 + 2 s^2 + s
2.0 21.386389751875043 1.0 0.6823278038280193
注意事項
前の例で使った2次系の伝達関数モデルはゲイン余裕と位相余裕が計算されず「inf inf nan nan」と出力されました。推測ですが、ゲインが 0dBを跨ぎ、位相が -180° を跨いでいないと上手く計算してくれないのではないかと思います。
ナイキスト線図作成
ナイキスト線図は関数 nyquist(sys) を使います。
引数sysには伝達関数を入れます。
ナイキスト線図は周波数応答 G( jω) の実部が横軸、虚部が縦軸の極座標系において、角周波数 ω を 0 から ∞ まで変化させた軌跡です。
from control import matlab
from matplotlib import pyplot as plt
sys = matlab.tf([2, 5, 1], [1, 2, 3])
print(sys)
matlab.nyquist(sys)
plt.show()結果
2 s^2 + 5 s + 1
---------------
s^2 + 2 s + 3
ステップ応答波形
ステップ応答波形は関数 step(sys, t) を使います。
2次系の伝達関数を例にグラフを描きます。
引数 sys は伝達関数、t は時間を入れます。
今回、時間は数値演算ライブラリ numpyを使っています。
# 数値演算ライブラリ
import numpy as np
# 制御工学ライブラリ
from control import matlab
# グラフ描画
from matplotlib import pyplot as plt
# 伝達関数モデル作成
zeta = 0.6
wn = 3.0
G = matlab.tf([wn*wn], [1.0, 2.0*zeta*wn, wn*wn])
print(G)
# ステップ応答
t = np.linspace(0, 5, 1000)
(yout, T) = matlab.step(G, t)
plt.plot(T, yout)
plt.grid()
plt.axhline(1, color="b", linestyle="--")
plt.xlim(0, 5)
plt.show()実行結果
パラメータは適当ですが、2次系の伝達関数で ζ＝0.6 としたのでちょとオーバーシュートした波形になります。
実際の出力結果もそのような波形になっています。
9
---------------
s^2 + 4.8 s + 9
ブロック線図（直列結合）
ブロック線図の直列結合を演算します。
関数 series( ) を使用するか、伝達関数 Ga * Gb でOKです。
# Pytho上のコード
# 制御工学ライブラリ
from control import matlab
# ブロック線図 直列
Ga = matlab.tf([1], [1, 1])
Gb = matlab.tf([2, 2], [2])
print(Ga)
print(Gb)
print(Ga * Gb)
Gab = matlab.series(Ga, Gb)
print(Gab)実行結果
1
-----
s + 1
2 s + 2
-------
2
2 s + 2
-------
2 s + 2
2 s + 2
-------
2 s + 2
ブロック線図（並列結合）
ブロック線図の並列結合を演算します。
関数 parallel( ) を使用するか、伝達関数 Ga + Gb でOKです。
# Pytho上のコード
# 制御工学ライブラリ
from control import matlab
# ブロック線図 並列
Ga = matlab.tf([1], [1, 1])
Gb = matlab.tf([2, 2], [2])
print(Ga)
print(Gb)
print(Ga + Gb)
Gab = matlab.parallel(Ga, Gb)
print(Gab)実行結果
1
-----
s + 1
2 s + 2
-------
2
2 s^2 + 4 s + 4
---------------
2 s + 2
2 s^2 + 4 s + 4
---------------
2 s + 2
ブロック線図（フィードバック）
ブロック線図のフィードバック結合を演算します。
関数 feedback( ) を使用します。
ネガティブフィードバックである点に注意してください。
# Pytho上のコード
# 制御工学ライブラリ
from control import matlab
# ブロック線図 フィードバック
Ga = matlab.tf([1], [1, 1])
Gb = matlab.tf([2, 2], [2])
print(Ga)
print(Gb)
Gab = matlab.feedback(Ga, Gb)
print(Gab)実行結果
1
-----
s + 1
2 s + 2
-------
2
2
-------
4 s + 4
まとめ
ここまでに記載した機能を使えば、伝達関数モデルを複数作成し、システムの閉ループ伝達関数を導出して、ステップ応答を導出することができると思います。また、ボード線図を作図して、ゲイン余裕や位相余裕を見ればシステムの安定性も検討できると思います。さらに、PI補償器やPID補償器を入れればシステムの制御設計も活用できるのではないかと思います。
