python-controlは、MATLABのControl System Toolboxに相当するPythonライブラリです。伝達関数や状態空間モデルの定義から、ステップ応答、ボード線図、根軌跡、ナイキスト線図まで、制御工学で必要な解析をPythonで行えます。
本記事では、python-controlの基本的な使い方を実例とともに解説します。
本記事の内容
- python-controlのインストールと基本操作
- 伝達関数・状態空間モデルの定義
- 時間応答の解析（ステップ応答・インパルス応答）
- 周波数応答の解析（ボード線図・ナイキスト線図）
- 根軌跡の描画
- フィードバック制御系の構成
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
インストール
pip install control matplotlib numpy
伝達関数の定義
基本的な伝達関数
import control as ctrl
import numpy as np
# 方法1: 分子・分母の係数で定義
# G(s) = (2s + 1) / (s^2 + 3s + 2)
G1 = ctrl.tf([2, 1], [1, 3, 2])
print("G1 =", G1)
# 方法2: s変数を使って定義
s = ctrl.tf('s')
G2 = (2*s + 1) / (s**2 + 3*s + 2)
print("G2 =", G2)
# 極と零点の確認
print(f"極: {ctrl.poles(G1)}")
print(f"零点: {ctrl.zeros(G1)}")
print(f"DCゲイン: {ctrl.dcgain(G1)}")
状態空間モデル
$$ \dot{\bm{x}} = \bm{A}\bm{x} + \bm{B}\bm{u}, \quad \bm{y} = \bm{C}\bm{x} + \bm{D}\bm{u} $$
import control as ctrl
import numpy as np
# 状態空間モデルの定義
A = np.array([[0, 1], [-2, -3]])
B = np.array([[0], [1]])
C = np.array([[1, 0]])
D = np.array([[0]])
sys_ss = ctrl.ss(A, B, C, D)
print("状態空間モデル:", sys_ss)
# 伝達関数に変換
sys_tf = ctrl.tf(sys_ss)
print("伝達関数:", sys_tf)
時間応答の解析
import numpy as np
import matplotlib.pyplot as plt
import control as ctrl
# 2次振動系
wn = 2.0
zeta = 0.3
G = ctrl.tf([wn**2], [1, 2*zeta*wn, wn**2])
t = np.linspace(0, 10, 500)
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
# ステップ応答
t_step, y_step = ctrl.step_response(G, T=t)
axes[0, 0].plot(t_step, y_step, 'b', linewidth=2)
axes[0, 0].set_title('Step Response')
axes[0, 0].set_xlabel('Time [s]')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].axhline(y=1.0, color='k', linestyle='--', alpha=0.3)
# インパルス応答
t_imp, y_imp = ctrl.impulse_response(G, T=t)
axes[0, 1].plot(t_imp, y_imp, 'r', linewidth=2)
axes[0, 1].set_title('Impulse Response')
axes[0, 1].set_xlabel('Time [s]')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].grid(True, alpha=0.3)
# 任意入力への応答（正弦波）
u_sin = np.sin(2 * t)
t_forced, y_forced = ctrl.forced_response(G, T=t, U=u_sin)
axes[1, 0].plot(t_forced, u_sin, 'g--', label='Input', linewidth=1)
axes[1, 0].plot(t_forced, y_forced, 'b', label='Output', linewidth=2)
axes[1, 0].set_title('Forced Response (sinusoidal input)')
axes[1, 0].set_xlabel('Time [s]')
axes[1, 0].set_ylabel('Amplitude')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
# 初期値応答（状態空間モデルで）
sys_ss = ctrl.tf2ss(G)
t_init, y_init = ctrl.initial_response(sys_ss, T=t, X0=[1, 0])
axes[1, 1].plot(t_init, y_init, 'm', linewidth=2)
axes[1, 1].set_title('Initial Condition Response')
axes[1, 1].set_xlabel('Time [s]')
axes[1, 1].set_ylabel('Amplitude')
axes[1, 1].grid(True, alpha=0.3)
plt.suptitle(f'2nd Order System (ωn={wn}, ζ={zeta})', fontsize=14)
plt.tight_layout()
plt.show()
周波数応答の解析
ボード線図
import numpy as np
import matplotlib.pyplot as plt
import control as ctrl
# 複数のシステムのボード線図
s = ctrl.tf('s')
G1 = 1 / (s + 1)
G2 = 10 / (s**2 + 2*s + 10)
G3 = 100 / (s * (s + 1) * (s + 10))
systems = [G1, G2, G3]
labels = ['1/(s+1)', '10/(s²+2s+10)', '100/(s(s+1)(s+10))']
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))
omega = np.logspace(-2, 3, 1000)
for G, label in zip(systems, labels):
mag, phase, w = ctrl.frequency_response(G, omega)
ax1.semilogx(w, 20 * np.log10(np.abs(mag)), linewidth=2, label=label)
ax2.semilogx(w, np.degrees(np.unwrap(np.angle(mag))), linewidth=2, label=label)
ax1.set_ylabel('Gain [dB]')
ax1.set_title('Bode Plot')
ax1.legend()
ax1.grid(True, which='both', alpha=0.3)
ax1.axhline(y=0, color='k', linewidth=0.5)
ax2.set_xlabel('Frequency [rad/s]')
ax2.set_ylabel('Phase [deg]')
ax2.legend()
ax2.grid(True, which='both', alpha=0.3)
ax2.axhline(y=-180, color='r', linewidth=0.5, linestyle='--')
plt.tight_layout()
plt.show()
ナイキスト線図
import numpy as np
import matplotlib.pyplot as plt
import control as ctrl
# 開ループ伝達関数
L = ctrl.tf([10], [1, 3, 3, 1])
# ナイキスト線図（手動描画）
omega = np.logspace(-2, 2, 1000)
mag, phase, w = ctrl.frequency_response(L, omega)
plt.figure(figsize=(8, 8))
plt.plot(np.real(mag), np.imag(mag), 'b', linewidth=2, label='ω > 0')
plt.plot(np.real(mag), -np.imag(mag), 'b--', linewidth=1, alpha=0.5, label='ω < 0')
plt.plot(-1, 0, 'r+', markersize=15, markeredgewidth=2, label='(-1, 0)')
plt.xlabel('Real')
plt.ylabel('Imaginary')
plt.title('Nyquist Plot')
plt.legend()
plt.grid(True, alpha=0.3)
plt.axis('equal')
plt.tight_layout()
plt.show()
フィードバック制御系の構成
import numpy as np
import matplotlib.pyplot as plt
import control as ctrl
# プラント
P = ctrl.tf([1], [1, 2, 1, 0]) # P(s) = 1/(s^3 + 2s^2 + s)
# PIDコントローラ
Kp, Ki, Kd = 30, 10, 8
C = ctrl.tf([Kd, Kp, Ki], [1, 0])
# 開ループ伝達関数
L = P * C
# 閉ループ伝達関数
G_cl = ctrl.feedback(L, 1)
# 安定余裕
gm, pm, wpc, wgc = ctrl.margin(L)
print(f"ゲイン余裕: {20*np.log10(gm):.2f} dB")
print(f"位相余裕: {pm:.2f} deg")
# 閉ループ極
print(f"閉ループ極: {ctrl.poles(G_cl)}")
# ステップ応答
t = np.linspace(0, 10, 1000)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 目標値応答
t_out, y_out = ctrl.step_response(G_cl, T=t)
axes[0].plot(t_out, y_out, 'b', linewidth=2)
axes[0].axhline(y=1.0, color='k', linestyle='--', alpha=0.3)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Output y(t)')
axes[0].set_title('Closed-loop Step Response')
axes[0].grid(True, alpha=0.3)
# 外乱応答
G_dist = ctrl.feedback(P, C)
t_d, y_d = ctrl.step_response(G_dist, T=t)
axes[1].plot(t_d, y_d, 'r', linewidth=2)
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Output (disturbance response)')
axes[1].set_title('Disturbance Response')
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
まとめ
本記事では、python-controlライブラリの基本的な使い方を解説しました。
ctrl.tf()で伝達関数、ctrl.ss()で状態空間モデルを定義できるctrl.step_response()やctrl.impulse_response()で時間応答を解析できるctrl.frequency_response()で周波数応答（ボード線図）を解析できるctrl.root_locus()で根軌跡を描画できるctrl.feedback()でフィードバック制御系を簡単に構成できるctrl.margin()でゲイン余裕・位相余裕を計算できる
次のステップとして、以下の記事も参考にしてください。