개요

공학을 제어하는 학습 노트.스프링·송·감진기계의 질점의 운동 방정식에서 전달 함수/상태 공간 모델을 구하고 제어계의 파이톤 프로그램 라이브러리'Python Control Systems Library를 사용하여 시뮬레이션한다.

관련 항목: RLC 직렬 회로의 전달 함수 상태 공간 모델과 Pythhon의 시뮬레이션을 이용하여

스프링, 송어, 감진기계 중 1자유도 질점의 운동 방정식

$$ m\ddot{y}(t) + c\dot{y}(t) + ky(t) = f(t)$$
$\ddot{y}(t)는 가속도, $\dot{y}(t)는 속도, $y(t)는 변위질량은 m$c이고 쇠퇴계수는 $k$k는 용수철 상수입니다. (모두 물리 상수이기 때문에 마이너스입니다.)

전달 함수를 통해 시스템을 모델링하다

힘$f(t)$를 시스템에 입력하고 $y(t)$s를 시스템에서 출력하는 출력으로 바꾸면 라프라스는 시스템의 전달 함수를 얻기 위해 $G(s)=Y(s)/F(s)$s로 바뀐다. 그러면 다음과 같다. (이때 속도와 위치의 초기값은 0이다. 즉, $\dot{y}(0)=0달러와 $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}$$

파이썬 모델링 시뮬레이션

Python의 Python Control Systems Library를 사용하여 시뮬레이션합니다.전달 함수를 설정하여 시스템을 모델링하고 펄스 응답을 시뮬레이션한다.
준비(GoogleColab. 환경)

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

펄스 응답의 시뮬레이션 (전달 함수)

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

m = 1    # 質量 [kg] 非負
c = 1    # 減衰係数 [N/m] 非負
k = 10   # バネ係数 [Ns/m] 非負

sys = ctrl.tf((1),(m,c,k)) # 伝達関数
print(sys)

t_range = (0,10) # 0～10秒の範囲をシミュレーション
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()

실행 결과

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

상태 공간 모형화

상태 공간 모델에서 여러 관측량(출력)을 설정할 수 있습니다.또한 임의의 초기값으로 계산할 수 있다.
시스템을 다음 상태 공간 모델로 모델링$\mathcal{P} 달러.

\mathcal{P}\,:\,\left\{
\begin{array}{l}
\dot{\boldsymbol{x}} = \mathbf{A}\boldsymbol{x} + \mathbf{B}\boldsymbol{u}& 状態方程式\\
\boldsymbol{y} = \mathbf{C}\boldsymbol{x} + \mathbf{D}\boldsymbol{u}& 出力方程式（観測方程式）
\end{array}
\right.

준비로 동작 방정식 $\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)$$
$x_1(t)=y(t)$、$x_2(t)=\dot{y}(t)달러, 상태확정$\boldsymbol{x}=[x1],\,x2]^T달러.변위 $y (t) $상태 $x1달러, 속도 $\dot{y}(t)$2달러입니다.
이것보다 $\dot{x}1(t)=\dot{y}(t)=x_2달러입니다.
또한 $\dot{x}2(t)=\ddot{y}(t)=-\frac{k}{m}(t)-\fracc{c}{m}\dot{y}(t)+\frac{1}{m}(t)$.
$u(t)$를 입력하면 다음과 같은 상태 방정식이 나타납니다.

\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

$$\dot{\boldsymbol{x}}
=\mathbf{A}\boldsymbol{x} +\mathbf{B} u$$
또한 다음 출력 방정식(관측 방정식)(위치 이동에 해당하는 상태 $x1달러와 속도에 해당하는 상태 $x2달러의 형식을 관찰)을 획득한다.

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$$
그러나 $\mathbf{D}=0$입니다.

파이썬 모델링 시뮬레이션

상태 공간 모델에서 펄스 응답을 시뮬레이션합니다.전달 함수와 달리 위치 이동의 초기값으로 $0.1달러를 설정하고 속도도 출력(관측)합니다.
다이어그램에서 일본어로 사용할 준비(Google Colab. 환경)

!pip install japanize-matplotlib

펄스 응답 시뮬레이션 (상태 공간 모델과 초기값 포함)

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

m = 1    # 質量 [kg] 非負
c = 1    # 減衰係数 [N/m] 非負
k = 10   # バネ係数 [Ns/m] 非負

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（=変位）の初期値を 0.1 に設定。x2（=速度）の初期値は 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(['変位','速度']) :
  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()

실행 결과는 다음과 같습니다.