将 CSTR 仿真逻辑解耦为传感器和 PID 控制

Question

非常感谢您提供的所有有用的教程。我正在尝试重新学习 10 年前回到 Chem E 时的基本 PID 控制。 我想解耦你的脚本。你能帮忙吗？

是否可以将此脚本分成两部分，其中一个脚本输出 Ca 和 Treactor，另一个输出 Tc？目前这只是出于我自己的好奇心。

耦合：

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

# from IMC tuning
Kc = 4.61730615181 * 2.0
tauI = 0.913444964569 / 4.0
tauD = 0.0

# define CSTR model
def cstr(x, t, u, Tf, Caf):
    Tc = u
    Ca = x[0]
    T = x[1]
    q = 100
    V = 100
    rho = 1000
    Cp = 0.239
    mdelH = 5e4
    EoverR = 8750
    k0 = 7.2e10
    UA = 5e4
    rA = k0 * np.exp(-EoverR / T) * Ca

    # Calculate concentration derivative
    dCadt = q / V * (Caf - Ca) - rA
    # Calculate temperature derivative
    dTdt = q / V * (Tf - T) \
           + mdelH / (rho * Cp) * rA \
           + UA / V / rho / Cp * (Tc - T)

    # Return xdot:
    xdot = np.zeros(2)
    xdot[0] = dCadt
    xdot[1] = dTdt
    return xdot

# Steady State Initial Conditions for the States
Ca_ss = 0.87725294608097
T_ss = 324.475443431599
x0 = np.empty(2)
x0[0] = Ca_ss
x0[1] = T_ss

# Steady State Initial Condition
u_ss = 300.0
# Feed Temperature (K)
Tf = 350
# Feed Concentration (mol/m^3)
Caf = 1

# Time Interval (min)
t = np.linspace(0, 10, 301)

# storage for recording values
op = np.ones(len(t)) * u_ss  # controller output
pv = np.zeros(len(t))  # process variable
e = np.zeros(len(t))  # error
ie = np.zeros(len(t))  # integral of the error
dpv = np.zeros(len(t))  # derivative of the pv
P = np.zeros(len(t))  # proportional
I = np.zeros(len(t))  # integral
D = np.zeros(len(t))  # derivative
sp = np.ones(len(t)) * T_ss  # set point

# define a setpoint ramp or steps
for i in range(15):
    sp[i * 20:(i + 1) * 20] = 300 + i * 7.0
sp[300] = sp[299]

# Upper and Lower limits on OP
op_hi = 350.0
op_lo = 250.0

pv[0] = T_ss

# Loop through time steps    
for i in range(len(t) - 1):
    delta_t = t[i + 1] - t[i]
    e[i] = sp[i] - pv[i]
    if i >= 1:  # calculate starting on second cycle
        dpv[i] = (pv[i] - pv[i - 1]) / delta_t
        ie[i] = ie[i - 1] + e[i] * delta_t
    P[i] = Kc * e[i]
    I[i] = Kc / tauI * ie[i]
    D[i] = -Kc * tauD * dpv[i]
    op[i] = op[0] + P[i] + I[i] + D[i]
    if op[i] > op_hi:  # check upper limit
        op[i] = op_hi
        ie[i] = ie[i] - e[i] * delta_t  # anti-reset windup
    if op[i] < op_lo:  # check lower limit
        op[i] = op_lo
        ie[i] = ie[i] - e[i] * delta_t  # anti-reset windup
    ts = [t[i], t[i + 1]]
    u[i + 1] = op[i]
    y = odeint(cstr, x0, ts, args=(u[i + 1], Tf, Caf))
    Ca[i + 1] = y[-1][0]
    T[i + 1] = y[-1][1]
    x0[0] = Ca[i + 1]
    x0[1] = T[i + 1]
    pv[i + 1] = T[i + 1]

op[len(t) - 1] = op[len(t) - 2]
ie[len(t) - 1] = ie[len(t) - 2]
P[len(t) - 1] = P[len(t) - 2]
I[len(t) - 1] = I[len(t) - 2]
D[len(t) - 1] = D[len(t) - 2]

# Construct results and save data file
data = np.vstack((t, u, T, Ca, op)).T  # vertical stack and transpose data
np.savetxt('CSTR_simulator.txt', data, delimiter=',')

plt.ioff()
plt.show()
par

我尝试将其作为 cstr：

import numpy as np
from scipy.integrate import odeint

# Define CSTR model
def cstr(x, t, u, Tf, Caf):
    Ca = x[0]
    T = x[1]

    q = 100
    V = 100
    rho = 1000
    Cp = 0.239
    mdelH = 5e4
    EoverR = 8750
    k0 = 7.2e10
    UA = 5e4
    rA = k0 * np.exp(-EoverR / T) * Ca

    dCadt = q / V * (Caf - Ca) - rA
    dTdt = q / V * (Tf - T) + mdelH / (rho * Cp) * rA + UA / V / rho / Cp * (u - T)

    xdot = np.zeros(2)
    xdot[0] = dCadt
    xdot[1] = dTdt
    return xdot

# Simulation function
def simulate_cstr(u, Tf, Caf, x0, t):
    Ca = np.ones(len(t)) * x0[0]
    T = np.ones(len(t)) * x0[1]
    for i in range(len(t) - 1):
        ts = [t[i], t[i + 1]]
        y = odeint(cstr, x0, ts, args=(u[i], Tf, Caf))
        Ca[i + 1] = y[-1][0]
        T[i + 1] = y[-1][1]
        x0[0] = Ca[i + 1]
        x0[1] = T[i + 1]
    return Ca, T

这是控制器

import numpy as np
import matplotlib.pyplot as plt
from cstr_reactor import simulate_cstr  # Assuming your file is named cstr_reactor.py

# PID parameters and initial conditions
Kc = 4.61730615181 * 2.0
tauI = 0.913444964569 / 4.0
tauD = 0.0

# Control loop
def pid_control(sp, T_ss, u_ss, t, Tf, Caf, x0):
    u = np.ones(len(t)) * u_ss
    op = np.ones(len(t)) * u_ss
    pv = np.zeros(len(t))
    e = np.zeros(len(t))
    ie = np.zeros(len(t))
    dpv = np.zeros(len(t))
    P = np.zeros(len(t))
    I = np.zeros(len(t))
    D = np.zeros(len(t))

    # Initialize Ca and T arrays
    Ca = np.ones(len(t)) * x0[0]
    T = np.ones(len(t)) * x0[1]

    # Upper and Lower limits on OP
    op_hi = 350.0
    op_lo = 250.0

    pv[0] = T_ss

    # Create plot
    plt.figure(figsize=(10, 7))
    plt.ion()
    plt.show()

    for i in range(len(t) - 1):
        delta_t = t[i + 1] - t[i]
        e[i] = sp[i] - pv[i]
        if i >= 1:
            dpv[i] = (pv[i] - pv[i - 1]) / delta_t
            ie[i] = ie[i - 1] + e[i] * delta_t
        P[i] = Kc * e[i]
        I[i] = Kc / tauI * ie[i]
        D[i] = -Kc * tauD * dpv[i]
        op[i] = op[0] + P[i] + I[i] + D[i]
        if op[i] > op_hi:
            op[i] = op_hi
            ie[i] = ie[i] - e[i] * delta_t
        if op[i] < op_lo:
            op[i] = op_lo
            ie[i] = ie[i] - e[i] * delta_t
        u[i + 1] = op[i]

        # Use the current Ca and T for the initial condition of the next step
        x0 = [Ca[i], T[i]]
        Ca, T = simulate_cstr(u, Tf, Caf, x0, t)
        pv[i + 1] = T[i + 1]

        # Debugging information
        if i % 50 == 0:
            print(f"Time: {t[i]:.2f}, Setpoint: {sp[i]:.2f}, PV: {pv[i]:.2f}, OP: {op[i]:.2f}, Ca: {Ca[i]:.2f}, T: {T[i]:.2f}")

        # Plot the results
        plt.clf()
        plt.subplot(4, 1, 1)
        plt.plot(t[0:i], u[0:i], 'b--', linewidth=3)
        plt.ylabel('Cooling T (K)')
        plt.legend(['Jacket Temperature'], loc='best')
        plt.subplot(4, 1, 2)
        plt.plot(t[0:i], Ca[0:i], 'g-', linewidth=3)
        plt.ylabel('Ca (mol/L)')
        plt.legend(['Reactor Concentration'], loc='best')

        plt.subplot(4, 1, 3)
        plt.plot(t[0:i], sp[0:i], 'r--', linewidth=2, label='Set Point')
        plt.plot(t[0:i], T[0:i], 'k:', linewidth=3, label='Reactor Temperature')
        plt.ylabel('T (K)')
        plt.xlabel('Time (min)')
        plt.legend(loc='best')
        plt.subplot(4, 1, 4)
        plt.plot(t[0:i], op[0:i], 'r--', linewidth=3, label='Controller Output (OP)')
        plt.plot(t[0:i], P[0:i], 'g:', linewidth=2, label='Proportional (Kc e(t))')
        plt.plot(t[0:i], I[0:i], 'b.-', linewidth=2, label='Integral (Kc/tauI * Int(e(t))')
        plt.plot(t[0:i], D[0:i], 'k-.', linewidth=2, label='Derivative (-Kc tauD d(PV)/dt)')
        plt.legend(loc='best')
        plt.ylabel('Output')
        plt.draw()
        plt.pause(0.01)

    op[len(t) - 1] = op[len(t) - 2]
    ie[len(t) - 1] = ie[len(t) - 2]
    P[len(t) - 1] = P[len(t) - 2]
    I[len(t) - 1] = I[len(t) - 2]
    D[len(t) - 1] = D[len(t) - 2]

    # Save data to file
    data = np.vstack((t, u, T, Ca, op)).T
    np.savetxt('data_doublet_steps.txt', data, delimiter=',')

    return u, T

# Main execution
if __name__ == "__main__":
    t = np.linspace(0, 10, 301)
    x0 = [0.87725294608097, 324.475443431599]
    sp = np.ones(len(t)) * 324.475443431599

    # Define a setpoint ramp or steps
    for i in range(15):
        sp[i * 20:(i + 1) * 20] = 300 + i * 7.0
    sp[300] = sp[299]

    u, T = pid_control(sp, 324.475443431599, 300.0, t, 350, 1, x0)

但是我从单个脚本中获得了不同的值。任何帮助将不胜感激。谢谢你。

我尝试将脚本分成两部分，如上所示。 这些值接近原始值但不准确，然后反应完全失控。所以我分割它们的方式显然是错误的，但我不知道在哪里。

Answer 1

这里是一个脚本，它根据您的要求演示了两个功能，其中包括执行器计算和测量检索（模拟阶段）两个基本阶段。 PID 控制使反应器不稳定，在不稳定区域时需要小的设定值变化。

# PID control based on new measurement
pid.setpoint=sp[i]
u[i+1] = pid(T[i]) + u_bias

# Simulate system for one time step with u[i+1] actuator value
ts = [t[i],t[i+1]]
y = odeint(cstr,x0,ts,args=(u[i+1],Tf,Caf))

它使用

simple_pid

包，但您也可以将 PID 算法编码为单独的自定义函数。输入是设定点和当前测量值。该函数的输出是新的执行器值。大多数工业应用都使用 PID 方程的速度形式，如 APMonitor 课程网站所示。

使用不稳定的系统学习 PID 控制可能会令人沮丧，因此我建议学习一些稳定的其他案例研究（混合浓度或水平），或者使用具有 Python 或 Matlab 物理硬件接口的温度控制实验室/模拟链接。一旦回到困难的 CSTR 问题，还有其他方法，例如线性和非线性模型预测控制。