MATLAB中带CSTR的PID控制器x3输出异常问题求助

我在MATLAB中完成含CSTR的PID控制器作业，已写出初步可运行代码，但绘制x3输出时图形混乱，迭代到26次后因数值溢出无法继续运行。系统目标是通过参考信号r(n)、控制信号u(n)、误差信号e(n)将x3控制在指定范围，目前未实现该功能，求解决思路与帮助。

相关参考示意图：

• 伪代码示意图
• 需求目标示意图

现有MATLAB代码 #

clear; clc; close all;

% Initialization
Kp0 = 1.5;
Ki0 = 0.1;
Kd0 = 0.25;
u=0:1;

% CSTR system parameters
Ts = 0.1;
Da1 = 3;
Da2 = 0.5;
Da3 = 1;

x1(1:5) = 0.5;
x2(1:5) = 0.5;
x3(1:5) = 1.5;

d2 = 1;

t = 1:20;
r = 0.45 + 0.15*sin(2*pi*(1/200)*t);

for n = 3:length(r)
    % Compute error
    e(n) = r(n) - x3(n);
     
    % Update PID gains
    Kp(n+1) = Kp0 + e(n);
    Ki(n+1) = Ki0 + e(n);
    Kd(n+1) = Kd0 + e(n);

    % PID Control Signal
    u(n) = u(n-1) + Kp(n) * (e(n) - e(n-1)) + Ki(n) * e(n) + Kd(n) * (e(n) - 2*e(n-1) + e(n-2));
    
    % Saturate Control Signal
    if u(n) > max(u)
        u(n) = max(u);
    elseif u(n) < min(u)
        u(n) = min(u);
    end
    
    % Apply u to CSTR system
    k1x1 = Ts*f1_func(x1(n), x2(n));
    k1x2 = Ts*f2_func(x1(n), x2(n), d2, u(n));
    k1x3 = Ts*f3_func(x2(n), x3(n), d2);

    k2x1 = Ts*f1_func(x1(n) + k1x1/2, x2(n) + k1x2/2);
    k2x2 = Ts*f2_func(x1(n) + k1x1/2, x2(n) + k1x2/2, d2, u(n));
    k2x3 = Ts*f3_func(x2(n) + k1x2/2, x3(n) + k1x3/2, d2);

    k3x1 = Ts*f1_func(x1(n) + k2x1/2, x2(n) + k2x2/2);
    k3x2 = Ts*f2_func(x1(n) + k2x1/2, x2(n) + k2x2/2, d2, u(n));
    k3x3 = Ts*f3_func(x2(n) + k2x2/2, x3(n) + k2x3/2, d2);

    k4x1 = Ts*f1_func(x1(n) + k3x1, x2(n) + k3x2);
    k4x2 = Ts*f2_func(x1(n) + k3x1, x2(n) + k3x2, d2, u(n));
    k4x3 = Ts*f3_func(x2(n) + k3x2, x3(n) + k3x3, d2);

    x1(n + 1) = x1(n) + 1/6 * (k1x1 + 2*k2x1 + 2*k3x1 + k4x1);
    x2(n + 1) = x2(n) + 1/6 * (k1x2 + 2*k2x2 + 2*k3x2 + k4x2);
    x3(n + 1) = x3(n) + 1/6 * (k1x3 + 2*k2x3 + 2*k3x3 + k4x3);
   

    % Plotting
    figure(1);
    subplot(3,1,1);
    plot(x3, 'r');
    title('x3 - r');
    subplot(3,1,2);
    plot(r);
end

问题排查与修改思路 #

• PID增益更新逻辑错误：当前直接用Kp(n+1) = Kp0 + e(n)更新增益，会让增益随误差无限制波动，极易导致系统发散。若用固定增益PID，直接移除增益更新代码，全程使用Kp0/Ki0/Kd0；若需自适应PID，需采用合理规则（如基于误差变化率调整，而非简单叠加误差）。
• 控制信号饱和逻辑无效：当前用max(u)/min(u)作为饱和阈值，会随u数组的变化不断改变，无法起到稳定的限幅作用。应设定固定的上下限，例如：

  u_max = 5; % 根据系统特性设定合理值
  u_min = -5;
  u(n) = max(min(u(n), u_max), u_min);
  

• 变量初始化不规范：
  • e、Kp、Ki、Kd未提前初始化，会导致n=3时e(n-1)等变量为未定义状态，MATLAB默认赋值0会干扰PID计算，需提前用zeros初始化数组；
  • x1/x2/x3无需初始化前5个元素，仅需初始化第一个元素，后续迭代自动扩展即可。
• 绘图逻辑导致图形混乱：每次循环都调用plot会重复绘制、刷新图形，应将所有绘图代码移到循环结束后统一执行，同时添加参考信号r的对比曲线，便于观察跟踪效果。
• 状态方程函数需验证：f1_func/f2_func/f3_func是CSTR系统的核心动力学模型，若函数实现不符合实际物理规律（如反应速率、换热关系写错），会直接导致x3等状态量发散，需严格对照CSTR的状态方程确认函数正确性。

内容的提问来源于stack exchange，提问作者Halil İbrahim Özen