Pythonで始めるPID制御シミュレーション：2槽タンクモデルで直感的理解

前のページで、ExcelとVBAを使ってPID制御を理解したな。
今回は、それをPythonでやってみよう。Pythonなら柔軟にパラメータをいじれて、
グラフ化も簡単だから、もっと直感的にPID制御を理解できるはずだ。

Pythonならscipyやmatplotlibを使って、さまざまな条件を即座に試せるわね。
オーバーシュートが出るか出ないか、チューニング次第でどう変わるか、
いろいろ実験できそう！

なぜPythonを使うのか

VBAやExcelでの実行も有効ですが、Pythonは以下のような利点があります：
・豊富なライブラリ：scipy.integrate.solve_ivpを用いて高度な数値解析が可能。
・手軽なチューニング：パラメータを簡単に変更し、結果を即時チェック。
・可視化が容易：matplotlibで応答特性を瞬時にプロット。

せやな。Pythonやと、パラメータをさくっと変えて、
どんな応答になるかすぐに確かめられるで。

ふむふむ。では、さっそく2槽タンクモデルとPID制御則を
Pythonで実装してみましょう！

モデルとPID制御則の再確認

2槽タンクモデルは以下の連立微分方程式で表されます：

$$\frac{dy_1}{dt} = \frac{1}{A_1}q_{in1} – \frac{1}{A_1 R_1}y_1$$
$$\frac{dy_2}{dt} = \frac{1}{A_2 R_1}y_1 – \frac{1}{A_2 R_2}y_2$$

PID制御則は、$q_{in1}(t)$を以下のように決定します。

$$q_{in1}(t) = q_{ss} + K_c \left((y_{2set}-y_2(t)) + \frac{1}{T_I}\int_0^t (y_{2set}-y_2(\tau))\,d\tau – T_D\frac{dy_2}{dt}\right)$$

偏差の積分項や微分項を簡単に扱うため、状態変数に偏差の積分値を持たせるで。

なるほど。状態変数は [y1, y2, e_int] として、
solve_ivp で解くわけね。

Pythonでの実装例

以下は実際のPythonコード例です。通常は隠しておき、必要に応じて表示できるようにしています。

Pythonコードを見る


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# パラメータ設定
A1 = 30.0
R1 = 0.08
A2 = 30.0
R2 = 0.08

q_ss = 50.0
y2set = 8.0
y1_init = 4.0
y2_init = 4.0

# PIDパラメータ例
Kc = 20.0
TI = 2.0
TD = 1.2

t_span = (0, 50)
t_eval = np.linspace(t_span[0], t_span[1], 501)

def tank_system(t, y):
    # y = [y1, y2, e_int]
    y1, y2, e_int = y
    e = y2set - y2
    dy2_dt = (1/(A2*R1))*y1 - (1/(A2*R2))*y2
    q_in1 = q_ss + Kc*(e + (e_int/TI) - TD*dy2_dt)
    dy1_dt = (1/A1)*q_in1 - (1/(A1*R1))*y1
    de_int_dt = e
    return [dy1_dt, dy2_dt, de_int_dt]

y0 = [y1_init, y2_init, 0.0]
sol = solve_ivp(tank_system, t_span, y0, t_eval=t_eval, method='RK45')

t = sol.t
y1_sol = sol.y[0]
y2_sol = sol.y[1]
e_int_sol = sol.y[2]

def q_in1_fun(y1, y2, e_int):
    e = y2set - y2
    dy2_dt = (1/(A2*R1))*y1 - (1/(A2*R2))*y2
    return q_ss + Kc*(e + (e_int/TI) - TD*dy2_dt)

q_in1_sol = q_in1_fun(y1_sol, y2_sol, e_int_sol)

plt.figure(figsize=(10,5))

plt.subplot(1,2,1)
plt.plot(t, y2_sol, label="y2(t)")
plt.axhline(y2set, color='r', linestyle='--', label='setpoint')
plt.xlabel("Time [h]")
plt.ylabel("Liquid Level y2 [m]")
plt.title("Level Response in Tank 2")
plt.grid(True)
plt.legend()

plt.subplot(1,2,2)
plt.plot(t, q_in1_sol, label="q_in1(t)")
plt.xlabel("Time [h]")
plt.ylabel("Flow Rate q_in1 [m^3/h]")
plt.title("Inflow Rate Controlled by PID")
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

どうや、Pythonでこの通りシミュレーションできたやろ？
このコードを元に、$K_c$や$T_I$、$T_D$を変えてみるといい。