Pythonによる制御工学入門第1回

参考書の第4章「制御対象の振る舞い」を自作の数値シミュレーションで再現する

Pythonによる制御工学入門第1回

controlモジュールのインポート

参考書だと、from control.matlab import * でインポートしているが、
個人的に、この方法だと、どの関数がcontrolモジュールのものかわからないので、下記でインポートするようにしている

import control.matlab as pyctrl

時間応答

1次遅れ系(台車のモデル)

質量 $M$ , 粘性摩擦力 $\mu$ とすると、運動方程式は下記になる

$\displaystyle M \ddot{x}(t) + \mu \dot{x}(t) = u(t)$

速度 $\dot{x}(t)$ に関する伝達関数は、

$\displaystyle \begin{align} ( M s + \mu ) X(s) &= U(s) \\ P(s) &= \frac{X(s)}{U(s)} = \frac{1}{M s + \mu} \\ &= \frac{\frac{1}{\mu}}{1 + \frac{M}{\mu} s} \\ &= \frac{K}{1 + T s} \end{align}$

$K = \frac{1}{\mu}$ , $T = \frac{M}{\mu}$ としている

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import control.matlab as pyctrl

(T, K) = (0.5, 1)
P = pyctrl.tf([0, K], [T, 1])
y, t = pyctrl.step(P, np.arange(0, 5, 0.01))

fig = plt.figure(figsize=(6, 5))
ax = fig.add_subplot(111)
ax.plot(t, y)

ax.axhline(0.632, color='k',lw=0.5)
ax.axvline(0.5, color='k',lw=0.5)
ax.scatter(0.5, 0.632, color='k',lw=0.5)
ax.annotate('$(0.5, 0.632)$', xy=(0.7, 0.5))

ax.set_xlabel('t')
ax.set_ylabel('y')
ax.grid()
fig.tight_layout()

1次遅れ系(数値シミュレーション)

運動方程式から、数値計算する

$\displaystyle \begin{align} \ddot{x}(t) &= - \frac{\mu}{M} \dot{x}(t) + \frac{1}{M} u(t) \\ \ddot{x}(t) &= - \frac{1}{T} \dot{x}(t) + \frac{T}{K} u(t) \end{align}$

(T, K) = (0.5, 1)
dt = 0.001
t = [0]
vel_s = [0.0]
acc_s = [0.0]

for i in range(5000):
    acc_s += [- 1.0 /  T * vel_s[-1] + T / K * 4.0]
    vel_s += [vel_s[-1] + acc_s[-1] * dt]
    t += [t[-1] + dt]

fig = plt.figure(figsize=(6, 5))
ax = fig.add_subplot(111)

ax.plot(t, vel_s)
ax.axhline(0.632, color='k',lw=0.5)
ax.axvline(0.5, color='k',lw=0.5)
ax.scatter(0.5, 0.632, color='k',lw=0.5)
ax.annotate('$(0.5, 0.632)$', xy=(0.7, 0.5))

ax.set_xlabel('t')
ax.set_ylabel('y')
ax.grid()
fig.tight_layout()

2次遅れ系

バネマスダンパ系で考える
質量 $M$ , 粘性摩擦力 $\mu$ , ばね定数 $k$ とすると、運動方程式は下記になる

$\displaystyle M \ddot{x}(t) + \mu \dot{x}(t) + k x(t) = u(t)$

速度 $x(t)$ に関する伝達関数は、

$\displaystyle \begin{align} ( M s^2 + \mu s + k ) X(s) &= U(s) \\ P(s) &= \frac{X(s)}{U(s)} = \frac{1}{M s^2 + \mu s + k} \\ &= \frac{\frac{1}{M}}{s^2 + \frac{\mu}{M} s + \frac{k}{M}} \\ &= \frac{K \omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2} \end{align}$

$\frac{1}{M} = K \omega_n ^2$ , $\frac{\mu}{M} = 2 \zeta \omega_n$ , $\frac{k}{M} = \omega_n ^2$ としている

(K, zeta, omega_n) = (1.0, 0.4, 5)
P = pyctrl.tf([0, K * omega_n**2], [1, 2 * zeta * omega_n, omega_n**2])
y, t = pyctrl.step(P, np.arange(0, 5, 0.01))

fig = plt.figure(figsize=(6, 5))
ax = fig.add_subplot(111)
ax.plot(t, y)

ax.set_xlabel('t')
ax.set_ylabel('y')
ax.grid()
fig.tight_layout()

2次遅れ系(数値シミュレーション)

$\displaystyle \begin{align} \ddot{x}(t) &= - \frac{\mu}{M} \dot{x}(t) - \frac{k}{M} x(t) + \frac{1}{M} u(t) \\ \ddot{x}(t) &= - 2 \zeta \omega_n \dot{x}(t) - \omega_n^2 x(t) + K \omega_n^2 u(t) \end{align}$

(K, zeta, omega_n) = (1.0, 0.4, 5)
dt = 0.001
t = [0]
pos_s = [0.0]
vel_s = [0.0]
acc_s = [0.0]

for i in range(5000):
    acc_s += [- 2.0 * zeta * omega_n * vel_s[-1] - omega_n**2 * pos_s[-1] + K * omega_n**2 * 1.0]
    vel_s += [vel_s[-1] + acc_s[-1] * dt]
    pos_s += [pos_s[-1] + vel_s[-1] * dt]
    t += [t[-1] + dt]

fig = plt.figure(figsize=(6, 5))
ax = fig.add_subplot(111)

ax.plot(t, pos_s)

ax.set_xlabel('t')
ax.set_ylabel('y')
ax.grid()
fig.tight_layout()

状態空間モデルの時間応答

$\displaystyle \frac{d}{dt} \begin{bmatrix} x(t) \\ \dot{x}(t) \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -4 & -5 \end{bmatrix} \begin{bmatrix} x(t) \\ \dot{x}(t) \end{bmatrix} + \begin{bmatrix} 0 \\ 1 \end{bmatrix} u(t)$

A = [[0, 1],[-4, -5]]
B = [[0], [1]]
C = np.eye(2)
D = np.zeros([2, 1])
P = pyctrl.ss(A, B, C, D)

零入力応答

$\displaystyle x(0) = \begin{bmatrix} -0.3 \\ 0.4 \end{bmatrix}$

Td = np.arange(0, 5, 0.01)
X0 = [-0.3, 0.4]
x, t = pyctrl.initial(P, Td, X0) #ゼロ入力応答

fig = plt.figure(figsize=(6, 5))
ax = fig.add_subplot(111)

ax.plot(t, x[:,0], label='$x_1$')
ax.plot(t, x[:,1], '-.', label='$x_2$')

ax.set_xlabel('t')
ax.set_ylabel('$x_1$,$x_2$')
ax.legend()
ax.grid()
fig.tight_layout()

零入力応答(数値シミュレーション)

dt = 0.001
t = [0]
pos_s = [-0.3]
vel_s = [0.4]
acc_s = [0.0]

for i in range(5000):
    acc_s += [- 4 * pos_s[-1] - 5 * vel_s[-1]]
    vel_s += [vel_s[-1] + acc_s[-1] * dt]
    pos_s += [pos_s[-1] + vel_s[-1] * dt]
    t += [t[-1] + dt]

fig = plt.figure(figsize=(6, 5))
ax = fig.add_subplot(111)

ax.plot(t, pos_s, label='$x_1$')
ax.plot(t, vel_s, '-.', label='$x_2$')

ax.set_xlabel('t')
ax.set_ylabel('$x_1$,$x_2$')
ax.legend()
ax.grid()
fig.tight_layout()