Pythonで伝達関数の解析をする話

こんばんは。

Python は、ちょっとしたことから凝ったことまで簡単にできるのが素晴らしいところです。
ただ、実行速度の遅さは難点ですね。

さて、今回は、ラプラス領域で記述された伝達関数の解析をする話です。
紹介するのは、 Python の Control ライブラリです。

次のように import しておきます。
import matplotlib.pyplot as plt
from control.matlab import tf, bode
import control as ctrl


ラプラス変換が何なのかは、ここでは省略します。

ラプラス領域での伝達関数は、多項式の分数で表すことができます。
伝達関数は、次のように定義します。
sys = tf(num, den)

num, den は、それぞれ、多項式の係数のリストです。

例えば、 G(s) = (s^2 + 2s + 3) / (4s^2 - 5s - 6) であれば、
sys = tf([1, 2, 3], [4, -5, -6])
とします。

直列（畳み込み）と並列（加算）は、次のようです。
sys = sys1 * sys2
sys = sys1 + sys2

フィードバックループは、次のようです。
sys = feedback(sys_g, sys_bf, sign=-1)
NFB なら sign=-1 、 PFB なら sign=1 ですし、他の倍数とすることもできます。
sys_bf は、 FB ループ内に入る伝達関数です。例えば、ノイズがループに乗る場合にバンドパスフィルタを入れたりできます。

Bode 線図は、次のようです。
mag, phase, omega = bode(sys, dB=True, Hz=True, deg=True, plot=True)
plt.show()

極と零点は、次のようです。（複素数のリストが返ってきます。）
poles = ctrl.poles(sys)
zeros = ctrl.zeros(sys)



実際に伝達関数を定義してみましょう。

大学で必ず習う、コイルに電圧を掛ける場合の、コイルを流れる電流を考えましょう。
入力を電圧 e 、出力を電流 i 、伝達関数を G(s) とします。
回路方程式は、抵抗を R [ohm] 、インダクタンスを L [H] としたら、
e = Ri + L di/dt
ですから、伝達関数は、出力 / 入力という定義より、
G(s) = 1 / (R + Ls)
と書けます。
従って、 Python では、

resistance:float = 10.0
inductance:float = 1 / 10**-3
sys_G = tf([1], [inductance, resistance])
※ 抵抗とインダクタンスはテキトーです。

となります。


次に、バネマスダンパ系を考えましょう。
入力を力 f 、出力を位置 x 、伝達関数を G(s) とします。
運動方程式は、質量 m [kg] 、粘性係数 η[kg/s] 、バネ定数 k [kg/s/s] としたら、
f = m d^2x/dt^2 + η dx/dt + kx
ですから、伝達関数は、
G(s) = 1 / (ms^2 + ηs + k)
と書けます。
従って、 Python では、

mass:float = 1.0
eta:float = 0.4
coe_k:float = 0.01
sys_G = tf([1], [mass, eta, coe_k])
※ 各定数はテキトーです。

となります。


とっても簡単ですね。

コンピュータとは、計算する機械ですから、ゲームばかりではなく、計算に使ってあげてもよいのですよ。


以上。