9
10

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

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

Last updated at Posted at 2017-07-17

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

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

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

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

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

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

式を変換すると、

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

伝達関数G(s)は

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

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

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

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

となる。

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

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 \\

とする。

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

とおくと、 

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

これを行列で書くと

\\
\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=x  とすると、\\

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

よって、


\\
\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}

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

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

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

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

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

をラプラス変換すると

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

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

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

伝達関数は、

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://github.com/nnn112358/python-control_test

9
10
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
9
10

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?