Transfer function / state space model of spring / mass / damper system and simulation by Python

Overview

Control engineering study notes. The transfer function / state space model is obtained from the equation of motion of the mass point of the spring-mass-damper system, and the Python library of the control system "[Python Control Systems Library](https://python-control.readthedocs.io/en/" 0.8.3 / index.html) ”is used for simulation. 2020-04-27_23h38_44.png

--Related entry: Transfer function / state space model of RLC series circuit and simulation by Python

Equation of motion of mass point with one degree of freedom in spring-mass-damper system

m\ddot{y}(t) + c\dot{y}(t) + ky(t) = f(t)

$ \ ddot {y} (t) $ is acceleration, $ \ dot {y} (t) $ is velocity, and $ y (t) $ is displacement. $ m $ is the mass, $ c $ is the damping coefficient, and $ k $ is the spring constant (both are physical constants and therefore non-negative values).

Model the system with a transfer function

Laplace transform the force $ f (t) $ as the ** input ** to the system and the displacement $ y (t) $ as the ** output ** from the system, and the system's ** transfer function ** $ G ( Finding s) = Y (s) / F (s) $ gives: (At this time, the initial values of velocity and displacement are zero, that is, $ \ dot {y} (0) = 0 $ and $ y. (0) = 0 $).

ms^2Y(s) + c sY(s) + kY(s) = F(s)
\big(ms^2+c s + k\big) Y(s) = F(s)
G(s) = \frac{Y(s)}{F(s)} = \frac{1}{m s^2 + c s + k}

Modeling / simulation with Python

Use Python's ** Python Control Systems Library ** for simulation. Set the transfer function to model the system and simulate the ** impulse response **.

Preparation (Google Colab.environment)


!pip install --upgrade sympy
!pip install slycot
!pip install control

Impulse response simulation (transfer function)


import numpy as np
import control.matlab as ctrl
import matplotlib.pyplot as plt

m = 1    #mass[kg]Non-negative
c = 1    #Attenuation coefficient[N/m]Non-negative
k = 10   #Spring constant[Ns/m]Non-negative

sys = ctrl.tf((1),(m,c,k)) #Transfer function
print(sys)

t_range = (0,10) #Simulate the range from 0 to 10 seconds
y, t = ctrl.impulse(sys, T=np.arange(*t_range, 0.01))

plt.figure(figsize=(7,4),dpi=120,facecolor='white')
plt.hlines(0,*t_range,colors='gray',ls=':')
plt.plot(t,y)
plt.xlim(*t_range)
plt.show()

Execution result


           1
-----------------------
1e-08 s^2 + 2e-05 s + 1

f1.png

Modeled with a state space model

In the state space model, multiple observables (outputs) can be set. In addition, it is possible to calculate with any initial value.

Model the system into the following state-space model $ \ mathcal {P} $.

python


\mathcal{P}\,:\,\left\{
\begin{array}{l}
\dot{\boldsymbol{x}} = \mathbf{A}\boldsymbol{x} + \mathbf{B}\boldsymbol{u}&Equation of state\\
\boldsymbol{y} = \mathbf{C}\boldsymbol{x} + \mathbf{D}\boldsymbol{u}&Output equation (observation equation)
\end{array}
\right.

As a preparation, transform the equation of motion into $ \ ddot {y} (t) = ... $.

m\ddot{y}(t) + c\dot{y}(t) + ky(t) = f(t)
m\ddot{y}(t) = - ky(t) - c\dot{y}(t) + f(t)
\ddot{y}(t) = - \frac{k}{m}y(t) - \frac{c}{m} \dot{y}(t) +\frac{1}{m} f(t)

As $ x \ _1 (t) = y (t) $, $ x \ _2 (t) = \ dot {y} (t) $, ** state ** $ \ boldsymbol {x} = \ [x \ _1 , , x \ _2] ^ T $ is defined. ** Displacement ** $ y (t) $ is the state $ x_1 $, ** Velocity ** $ \ dot {y} (t) $ is the state $ x_2 $. </ font>

From this, $ \ dot {x} \ _1 (t) = \ dot {y} (t) = x \ _2 $.

Also, $ \ dot {x} \ _2 (t) = \ ddot {y} (t) =-\ frac {k} {m} y (t)-\ frac {c} {m} \ dot {y} (t) + \ frac {1} {m} f (t) $.

If the force $ f (t) $ is the input $ u (t) $, the following ** equation of state ** is obtained from the above.

python


\left[
    \begin{array}{c}
      \dot{x}_1 \\
      \dot{x}_2
    \end{array}
  \right]
=\left[
    \begin{array}{cc}
      0 & 1  \\
      -\frac{k}{m} & -\frac{c}{m} 
    \end{array}
  \right]
  \left[
    \begin{array}{c}
      x_1 \\
      x_2
    \end{array}
  \right]
+ \left[
    \begin{array}{c}
      0 \\
      \frac{1}{m}
    \end{array}
  \right] u

python


\dot{\boldsymbol{x}} 
=\left[
    \begin{array}{cc}
      0 & 1  \\
      -\frac{k}{m} & -\frac{c}{m} 
    \end{array}
  \right]
\boldsymbol{x} + \left[
    \begin{array}{c}
      0 \\
      \frac{1}{m}
    \end{array}
  \right] u
\dot{\boldsymbol{x}} =\mathbf{A}\boldsymbol{x} + \mathbf{B} u

In addition, the following ** output equation ** (observation equation) is obtained (state $ x_1 $ corresponding to ** displacement ** and state $ x_2 $ corresponding to ** velocity ** are observed).

python


y = \left[
    \begin{array}{cc}
      1 & 0 \\
      0 & 1 \\
    \end{array}
  \right]   \left[
    \begin{array}{c}
      x_1 \\
      x_2
    \end{array}
  \right]
y =\mathbf{C}\boldsymbol{x} + \mathbf{D} u

However, $ \ mathbf {D} = 0 $.

Modeling / simulation with Python

Simulate the ** impulse response ** in the state-space model. Unlike the case of the transfer function, $ 0.1 $ is set as the initial value of the displacement, and the velocity is also output (observed).

Preparation for using Japanese in graphs (Google Colab.environment)


!pip install japanize-matplotlib

Impulse response simulation (state space model with initial values)


import numpy as np
import control.matlab as ctrl
import matplotlib.pyplot as plt
import japanize_matplotlib

m = 1    #mass[kg]Non-negative
c = 1    #Attenuation coefficient[N/m]Non-negative
k = 10   #Spring constant[Ns/m]Non-negative

A = [[0,1],[-k/m, -c/m]]
B = [[0],[1/m]]
C = [[1,0],[0,1]]
D = 0

sys = ctrl.ss(A,B,C,D)
#print(sys)

t_range = (0,10)

# x1(=Initial value of displacement) is 0.Set to 1. x2 (=Initial value of speed) is set to 0
y, t = ctrl.impulse(sys, T=np.arange(*t_range, 0.01), X0=[0.1, 0.0])

fig, ax = plt.subplots(nrows=2, figsize=(7,5),dpi=120,sharex='col')
fig.patch.set_facecolor('white') 
fig.subplots_adjust(hspace=0.1)  
for i,s in enumerate(['Displacement','speed']) :
  ax[i].hlines(0,*t_range,colors='gray',ls=':')
  ax[i].plot(t,y[:,i])
  ax[i].set_ylabel(s)
  ax[i].set_xlim(*t_range)
plt.show()

The execution result is as follows.
f2.png

Recommended Posts

Transfer function / state space model of spring / mass / damper system and simulation by Python
Transfer function / state space model of RLC series circuit and simulation by Python
[Control engineering] Calculation of transfer functions and state space models by Python
Custom state space model in Python
[Scientific / technical calculation by Python] Fitting by nonlinear function, equation of state, scipy
Implementation of particle filters in Python and application to state space models
Implement a model with state and behavior (3) --Example of implementation by decorator
Explanation of production optimization model by Python
I implemented "Basics of Time Series Analysis and State Space Model" (Hayamoto) with pystan
[Python of Hikari-] Chapter 06-02 Function (argument and return value 1)