search
LoginSignup
11
Help us understand the problem. What are the problem?

More than 5 years have passed since last update.

posted at

updated at

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

目的

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

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
What you can do with signing up
11
Help us understand the problem. What are the problem?