PythonControlでインパルス応答を求める。

###目的
Python Controlで、1自由度系(バネマスダンバ系)の伝達関数のインパルス応答を求める。

###事前準備
PythonControlをインストールする

###システムのモデル化
以下の状態方程式について考える。

m:質量[kg]
k:ばね定数[N/m]
c:粘性減衰係数[N・s/m]
f：力[N]

G(s)=\frac{Y(s)}{U(s)}=\frac{1}{ms^2+cs+k}

###インパルス応答
Python ControlはInitial応答は以下の関数で求める。

control.matlab.Initial(sys, T=None, X0=0.0, input=0, output=None, return_x=False)

• Parameters:
  • sys:LTIシステム(線型時不変システム)の 状態空間or伝達関数
  • T:時間ベクトル（optionで与えられていなければ、自動計算する）
  • X0: 初期値(default=0)、数値は配列に自動変換。
  • input: int型、このシミュレーションへの入力。
  • output: int型、このシミュレーションからの出力。

control.matlab.Initialは線形システムのインパルス応答を求めることができる。
システムに複数の入力または出力（MIMO）がある場合は、シミュレーション用に1つの入力を選択する必要がある。
オプションで、1つの出力を選択することができる。
出力が選択されていない場合、すべての出力する。
パラメータのinputとoutputがこれを行う。その他の入力はすべて0に設定され、他のすべての出力は無視される。

引用元：
control.matlab.Initial

###サンプルコード
responce_impulse.py

#!/usr/bin/env python
from control.matlab import *
from matplotlib import pyplot as plt
from scipy import arange 

def main():
    k=1.0
    m=0.1
    c=0.1
    num = [0, 0,1] 
    den = [m, c, k]
    sys1 = tf(num, den) 
    print sys1
    (y1a, T1a) = initial(sys1,X0 = [0, 1],T = arange(0, 10, 0.01))
    plt.axhline(0, color="b", linestyle="--")
    plt.plot(T1a, y1a)
    plt.show()
    
if __name__ == "__main__":
  main()

サンプルコードは以下に格納。
https://github.com/nnn112358/python-control_test

###参考
PythonControlをインストールする

PythonControlで1自由度系の伝達関数を求める。

PythonControlで2自由度系の伝達関数を求める。