Python NumPy SciPy : 状態方程式によるボード線図

前回 ではラプラス変換から伝達関数を求め、ボード線図を書く方法を紹介しました。 その他にも lti, bode 関数では状態方程式の係数行列を与えるとボード線図が書けてしまいます。 制御の分野では簡潔で見通しが良くなる理由から、連立微分方程式を行列形式の一階微分方程式である状態方程式に書き変えて、系の応答特性を解析します。

前回の一自由度粘性減衰系を例とし、微分方程式を一階微分方程式に書き変え、

\begin{eqnarray} u(t) &=& \frac{1}{k}f(t) \nonumber \\ x(t) &=& x_1(t) \nonumber \\ \dot{x}_1(t) &=& x_2(t) \nonumber \\ \dot{x}_2(t) &=& -\frac{c}{m}x_2(t)-\frac{k}{m}x_1(t)+\frac{k}{m}u(t) \nonumber \end{eqnarray}

これを行列形式で表すと

\begin{eqnarray} \frac{d}{dt}\begin{bmatrix} {x_1(t)}\\ {x_2(t)} \end{bmatrix} &=& \begin{bmatrix} {0} & {1}\\ {-\frac{k}{m}} & {-\frac{c}{m}} \end{bmatrix} \begin{bmatrix} {x_1(t)}\\ {x_2(t)} \end{bmatrix} + \begin{bmatrix} {0} \\ {\frac{k}{m}} \end{bmatrix} u(t) \nonumber \\ u(t) &=& \begin{bmatrix} {1} & {0} \end{bmatrix} \begin{bmatrix} {x_1(t)}\\ {x_2(t)} \end{bmatrix} \nonumber \end{eqnarray}

となります。各項の係数行列を定義し、lti, bode 関数を用いると、簡単にボード線図が書けます。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# パラメータ設定
m = 1
c = 1
k = 400

A = np.array([[0, 1], [-k/m, -c/m]])
B = np.array([[0], [k/m]])
C = np.array([1, 0])
D = np.array([0])

s1 = signal.lti(A, B, C, D)
w, mag, phase = signal.bode(s1, np.arange(1, 500, 1)) 

# プロット
plt.figure(1)
plt.subplot(211)
plt.loglog(w, 10**(mag/20))
plt.ylabel("Amplitude")
plt.axis("tight")
plt.subplot(212)
plt.semilogx(w, phase)
plt.xlabel("Frequency[Hz]")
plt.ylabel("Phase[deg]")
plt.axis("tight")
plt.ylim(-180, 180)
# plt.show()

コメント

Comments powered by Disqus
書籍更新情報
2017-02-18
Pythonによる科学技術計算 基礎編
1.4版への更新が可能になりました。
サポートページはこちら
電子書籍