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.

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

$ \ 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 $).

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

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) = ... $.

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>

\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

\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

y = \left[
    \begin{array}{cc}
      1 & 0 \\
      0 & 1 \\
    \end{array}
  \right]   \left[
    \begin{array}{c}
      x_1 \\
      x_2
    \end{array}
  \right]

!pip install japanize-matplotlib

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()