ベクトル図とボード線図

はじめ

制御工学では、系に対して周波数ごとの安定性を議論するためにベクトル線図やボード線図といったグラフが活用される。今回は、それらをPythonの複素数計算機能を用いて周波数伝達関数から導出することを目的とする。具体的には、一次遅れ系と二次遅れ系についてグラフを描写してみる。

一次遅れ系

一次遅れ系とは、以下に示すような伝達関数で示されるような系のことをいう。

G(s)=\frac{1}{Ts+1}

今回はRL直列回路に入力電圧を与えたときの電流出力について考察する。

ベクトル図

python itiji_vector.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
omega = np.arange(0.0, 100.0, 0.001)

R=1.0
L=1.0

s=1j*omega

#伝達関数G（s）→G(jω)
G=R/(L*s+R)
x=G.real
y=G.imag

plt.title("一次遅れ系")
plt.xlabel('実軸')
plt.ylabel('虚軸')

plt.plot(x, y, 'b-')
plt.show()

これを実行すると以下のようなグラフが出力される。

ボード線図

ボード線図は、横軸を周波数（対数表示）縦軸を振幅（対数表示）、位相としたものであり、幅広い周波数範囲において特徴を大域的に捉えることに適している。

振幅について

それでは、振幅についてのプログラムを以下に示す。

python iciji_bode_gain.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
omega = np.arange(0.0, 100.0, 0.001)


s=1j*omega
T=0.1
#伝達関数G（s）→G(jω)
G=1/(s*T+1)
#ゲインの定義
g=20*np.log10(abs(G))
#横軸対数表示
plt.gca().set_xscale('log')
plt.plot(omega, g, 'b-')
plt.title('ボード線図（ゲイン）')
plt.xlabel('ω[rad/s]')
plt.ylabel('ゲイン[dB]')

plt.show()

位相について

次に、位相についてのプログラムとグラフを以下に示す。

python ichiji_bode_theta.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

omega = np.arange(0.0, 10000.0, 0.001)



s=1j*omega

#伝達関数G（s）→G(jω)
s=1j*omega
T=0.1
#伝達関数G（s）→G(jω)
G=1/(s*T+1)


#位相角の定義

#複素数の偏角（弧度法）を度数法に直す

theta = np.degrees(np.angle(G))


#横軸対数表示
plt.gca().set_xscale('log')
plt.plot(omega, theta, 'b-')
plt.title('ボード線図（位相角）')
plt.xlabel('ω[rad/s]')
plt.ylabel('位相[度]')

plt.show()

これにより、生成されるグラフを以下に示す。

二次遅れ系

二次遅れ系は、RLC直列回路のような振動する系のモデルであり以下の伝達関数で定義される。

G(s) =\frac{\omega_n^2}{s^2+2\zeta \omega_n +\omega^2_n}

ただし、$\zeta$は減衰係数で$\omega_n$は固有角速度である。

ベクトル図

以下のようなプログラムを書く

python niji_vector.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
omega = np.arange(0.0, 100.0, 0.001)

jita=0.5
omega_n=10

# k=omega_n**2/((omega_n**2-omega**2)**2+(2*omega_n*jita*omega)**2)

# x= k*(omega_n**2-omega**2)
# y= k*(-2*omega_n*jita*omega)

s=1j*omega

#伝達関数G（s）→G(jω)
G=omega_n**2/(s**2+2*jita*omega_n*s+omega_n**2)
x=G.real
y=G.imag

plt.title("二次遅れ系%f"%jita)
plt.xlabel('実軸')
plt.ylabel('虚軸')

plt.plot(x, y, 'b-')
plt.show()

これを実行すると以下のような画像が出力される。

ボード線図

一次遅れ系と同様にプログラムを実行することでグラフを描写することを試みる。

振幅について

それでは、振幅についてプログラムを示す。

python niji_bode_gain.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
omega = np.arange(0.0, 100.0, 0.001)

jita=0.5
omega_n=10

s=1j*omega

#伝達関数G（s）→G(jω)
G=omega_n**2/(s**2+2*jita*omega_n*s+omega_n**2)

#ゲインの定義
g=20*np.log10(abs(G))
#横軸対数表示
plt.gca().set_xscale('log')
plt.plot(omega, g, 'b-')
plt.title('ボード線図（ゲイン）')
plt.xlabel('ω[rad/s]')
plt.ylabel('ゲイン[dB]')

plt.show()

位相について

プログラムを以下に示す。

python niji_bode_theta.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

omega = np.arange(0.0, 10000.0, 0.001)

jita=0.5
omega_n=10

s=1j*omega

#伝達関数G（s）→G(jω)
G=omega_n**2/(s**2+2*jita*omega_n*s+omega_n**2)



#位相角の定義

#複素数の偏角（弧度法）を度数法に直す

theta = np.degrees(np.angle(G))


#横軸対数表示
plt.gca().set_xscale('log')
plt.plot(omega, theta, 'b-')
plt.title('ボード線図（位相角）')
plt.xlabel('ω[rad/s]')
plt.ylabel('位相[度]')

plt.show()

このプログラムを実行した場合以下のような画像が出力される。

まとめ

今回は、制御工学におけるベクトル線図やボード線図についてライブラリをあまり使用せずに、複素数の機能のみを用いて描写することができるかを試みた。結果、伝達関数から周波数伝達関数に変換させて、グラフを描写することができるプログラムを描写することができた。複素数は、電気工学や制御工学、流体力学などで使用されるため、それらを扱うことのできるPythonとは相性がいいと考えられる。