Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationEventAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
9
Help us understand the problem. What are the problem?

More than 3 years have passed since last update.

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

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
9
Help us understand the problem. What are the problem?