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

##目的
1自由度系(バネマスダンバ系)の状態方程式を立式し、Python Controlで、状態方程式から伝達関数へと変換する。
比較として、ラプラス変換から、伝達関数を求める方法も記載も記載する。

###事前準備
[PythonControlをインストールする](https://qiita.com/nnn112358/items/ada16e29ef8282498bb7)


###状態方程式と伝達関数
状態方程式とは、以下の式によって表される。

```math
\dot{x}(t)=Ax(t)+Bu(t) \\
y(t)=Cx(t)+Du(t)
```

状態方程式をラプラス変換すると、

```math
sX(s)=AX(s)+BU(s) \\
Y(s)=CX(s)+DU(s)
```

式を変換すると、

```math
sX(s)-AX(s)=BU(s) \\
(sI-A)X(s)=BU(s) \\
X(s)={(sI-A)}^{-1} BU(s) \\
```

伝達関数G(s)は

```math
G(s)=\frac{Y(s)}{U(s)}=\frac{{(sI-A)}^{-1} BU(s) +DU(s)}{U(s)}=(sI-A)^{-1} B +DU
```
から求まる。

Python-Controlでは、
control.ss2tf関数にA,B,C,Dの行列を入力すると、上式に従い、伝達関数を求めることができる。
http://python-control.readthedocs.io/en/latest/generated/control.ss2tf.html



###状態方程式の立式
![Banemasu.jpg](https://qiita-image-store.s3.amazonaws.com/0/73020/02289ef7-19ae-31dc-aca0-90e5ed8dfc02.jpeg)


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

とおくと、
上記の運動方程式は、

```math
m\ddot{x}+c\dot{x}+kx=f \\
```

となる。


###状態方程式から伝達関数を求める。
ここで式を変形して、

```math
m\ddot{x}+c\dot{x}+kx=f \\
m\ddot{x}=-c\dot{x}-kx+f \\
\ddot{x}=-c/m \dot{x}-k/m x+f/m \\
```

とする。

```math
x_{1}=x  \\
x_{2}=\dot{x} \\
```

とおくと、　

```math
\\
\dot x_{1} ={x_{2}} \\
\dot x_{2}= -c/m \dot{x_{2}}-k/m x_{1}+f/m \\
```

これを行列で書くと

```math
\\
\begin{bmatrix}
\dot x_{1}  \\
\dot x_{2}  
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
-k/m & -c/m
\end{bmatrix}
\begin{bmatrix}
x_{1}  \\
x_{2}  
\end{bmatrix}
+
\begin{bmatrix}
0 \\
1/m 
\end{bmatrix}
f
```

さらに、

```math
y=x  とすると、\\

y=
\begin{bmatrix}
1 & 0
\end{bmatrix}
\begin{bmatrix}
x_{1}  \\
x_{2}  
\end{bmatrix}
```


よって、

```math

\\
\begin{bmatrix}
\dot x_{1}  \\
\dot x_{2}  
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
-k/m & -c/m
\end{bmatrix}
\begin{bmatrix}
x_{1}  \\
x_{2}  
\end{bmatrix}
+
\begin{bmatrix}
0 \\
1/m 
\end{bmatrix}
f \\
y=
\begin{bmatrix}
1 & 0
\end{bmatrix}
\begin{bmatrix}
x_{1}  \\
x_{2}  
\end{bmatrix}
```

状態方程式と比較すると、

```math
x(t)=
\begin{bmatrix}
x_{1}  \\
x_{2}  
\end{bmatrix}

,
u(t)=f
,

A=
\begin{bmatrix}
0 & 1 \\
-k/m & -c/m
\end{bmatrix}
,
B=
\begin{bmatrix}
0 \\
1/m 
\end{bmatrix}
,
C=
\begin{bmatrix}
1 & 0 
\end{bmatrix}
,
D=[0]
```
となる。

以上から、
python controlで伝達関数を求めるコードは以下の様になる。

-signal_degree_of_freedom_system_1.py

```
#!/usr/bin/env python
from control.matlab import *
from matplotlib import pyplot as plt
    
def main():
    k=3.0
    m=0.1
    c=0.01
    A = [[0., 1], [-k/m, -c/m]]
    B = [[0.], [1./m]]
    C = [[1., 0.0]]
    D = [[0.]]
    sys1 = ss2tf(A, B, C, D)
    print sys1
   
    bode(sys1)    
    plt.show()
    
if __name__ == "__main__":
  main()

```

実行結果

```
$ ./signal_degree_of_freedom_system_1.py
       10
----------------
s^2 + 0.1 s + 30
```

![Figure_1.png](https://qiita-image-store.s3.amazonaws.com/0/73020/ddb6004a-3228-3e71-ed51-0023191ee825.png)

###状態方程式をラプラス変換して伝達関数を求める。

状態方程式の立式から、u(t)=fとした場合、

```math
m\ddot{x}(t)+c\dot{x}(t)+kx(t)=u(t) \\
```

をラプラス変換すると

```math
ms^2X(s)+cX(s)+kX(s)=U(s) \\
```

さらに、y=xとすると、

```math
ms^2Y(s)+cY(s)+kY(s)=U(s) \\
\frac{Y(s)}{U(s)}=\frac{1}{ms^2+cs+k}

```



伝達関数は、

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

以上から、
python controlで伝達関数を求めるコードは以下の様になる。

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

def main():
  k=3.0
  m=0.1
  c=0.01
  num = [0, 0,1] 
  den = [m, c, k]
  sys = tf(num, den) 
  print sys
  bode(sys)    
  plt.show()

if __name__ == "__main__":
  main()
```

実行結果

```
$ ./signal_degree_of_freedom_system_2.py 
          1
--------------------
0.1 s^2 + 0.01 s + 3
```

![Figure_1-2.png](https://qiita-image-store.s3.amazonaws.com/0/73020/29808376-2038-db7a-bff0-252ef99b9e43.png)



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