一文读懂：控制界的万能公式——PID算法到底是什么？

%% PID温度控制系统模拟
% 场景：电热水壶温度控制，目标温度60℃
clear; clc; close all;
%% 1. 系统模型参数（一阶滞后系统 + 纯延迟）
% 实际物理系统: G(s) = K / (tau*s + 1) * exp(-L*s)
K = 1.2;        % 系统增益（℃/W），表示每瓦功率能产生的温升
tau = 30;       % 时间常数（秒），系统响应快慢
L = 5;          % 纯延迟时间（秒），加热到温度传感器响应的延迟
% 离散化参数
dt = 0.1;       % 采样时间（秒）
sim_time = 300; % 仿真总时长（秒）5分钟
N = sim_time / dt; % 仿真步数
% 初始化变量数组
t = (0:N-1) * dt;           % 时间向量
u = zeros(1, N);            % 控制量（加热功率，0-100%）
y = zeros(1, N);            % 实际温度输出
y(1) = 20;                  % 初始温度20℃
setpoint = 60;              % 目标温度60℃
%% 2. PID参数（经过整定的值）
Kp = 2.5;      % 比例系数 - 响应速度
Ki = 0.08;     % 积分系数 - 消除稳态误差  
Kd = 8.0;      % 微分系数 - 抑制超调
% 积分项变量
integral = 0;
prev_error = 0;
% 限幅参数
u_max = 100;    % 最大加热功率100%
u_min = 0;      % 最小加热功率0%
% 抗积分饱和标志
anti_windup = 1; % 1:启用抗积分饱和
%% 3. 仿真循环 - PID控制
for k = 2:N
    % 当前误差
    error = setpoint - y(k-1);
    % --- PID计算 ---
    % 比例项
    P = Kp * error;
    % 积分项（带抗积分饱和）
    integral = integral + Ki * error * dt;
    if anti_windup
        % 如果控制量已经饱和，停止积分累加
        u_temp = P + integral + Kd * (error - prev_error)/dt;
        if (u_temp >= u_max && error > 0) || (u_temp <= u_min && error < 0)
            integral = integral - Ki * error * dt; % 撤销本次积分
        end
    end
    % 微分项（使用测量值微分，避免设定点突变引起的冲击）
    % 实际工程中常用: D = -Kd * (y(k-1) - y(k-2))/dt
    if k > 2
        D = -Kd * (y(k-1) - y(k-2)) / dt;
    else
        D = 0;
    end
    % PID输出
    u(k) = P + integral + D;
    % 控制量限幅
    u(k) = max(u_min, min(u_max, u(k)));
    % --- 系统模型：一阶滞后 + 延迟 ---
    % 获取延迟后的控制量（延迟L秒）
    delay_steps = round(L / dt);
    if k > delay_steps
        u_delayed = u(k - delay_steps);
    else
        u_delayed = 0;
    end
    % 一阶系统差分方程（欧拉法）
    dT = (K * u_delayed - (y(k-1) - 20)) / tau;
    y(k) = y(k-1) + dT * dt;
    % 添加少量测量噪声（可选，取消注释即可启用）
    % y(k) = y(k) + randn * 0.05;
    % 保存误差用于下次微分计算
    prev_error = error;
end
%% 4. 性能指标计算（修正后的版本）
% 稳态误差（最后50秒的平均误差）
steady_start = round(N - 500/dt); % 最后50秒的起始索引
if steady_start < 1
    steady_start = 1;
end
steady_state = mean(abs(setpoint - y(steady_start:end)));
% 超调量
overshoot = (max(y) - setpoint) / setpoint * 100;
if overshoot < 0
    overshoot = 0;
end
% 上升时间（从10%到90%的上升时间，或首次到达设定值的时间）
% 这里简化为首次达到设定值的时间
rise_idx = find(y >= setpoint, 1);
if isempty(rise_idx)
    rise_time = sim_time;
else
    rise_time = t(rise_idx) - t(find(y > 20, 1));
end
% 调节时间（进入±2%误差带并保持的时间）
settling_time = find_settling_time(t, y, setpoint, 0.02);
fprintf('========== 控制性能指标 ==========\n');
fprintf('稳态误差: %.2f ℃\n', steady_state);
fprintf('超调量: %.2f %%\n', overshoot);
fprintf('上升时间: %.1f 秒\n', rise_time);
fprintf('调节时间(±2%%): %.1f 秒\n', settling_time);
fprintf('最高温度: %.2f ℃\n', max(y));
fprintf('最终温度: %.2f ℃\n', y(end));
%% 5. 绘制结果
figure('Position', [100, 100, 1200, 800]);
% 子图1：温度响应曲线
subplot(2,2,1);
plot(t, y, 'b-', 'LineWidth', 2); hold on;
plot(t, setpoint*ones(1,N), 'r--', 'LineWidth', 1.5);
plot(t, 20*ones(1,N), 'k:', 'LineWidth', 1);
xlabel('时间 (秒)'); ylabel('温度 (℃)');
title('PID温度控制响应曲线');
legend('实际温度', '设定值 (60℃)', '初始温度', 'Location', 'southeast');
grid on;
xlim([0 sim_time]);
ylim([15 75]);
% 标注性能指标
text(10, 72, sprintf('超调量: %.1f%%', overshoot), 'FontSize', 10, 'BackgroundColor', 'w');
text(10, 69, sprintf('上升时间: %.1fs', rise_time), 'FontSize', 10, 'BackgroundColor', 'w');
text(10, 66, sprintf('稳态误差: %.2f℃', steady_state), 'FontSize', 10, 'BackgroundColor', 'w');
% 子图2：控制量输出（修改：不使用yline，改用line函数）
subplot(2,2,2);
plot(t, u, 'g-', 'LineWidth', 1.5);
xlabel('时间 (秒)'); ylabel('加热功率 (%)');
title('PID控制输出（占空比）');
grid on;
xlim([0 sim_time]);
ylim([-5 105]);
% 绘制上下限线（兼容旧版本MATLAB）
hold on;
h1 = line([0 sim_time], [100 100], 'Color', 'r', 'LineStyle', '--', 'LineWidth', 1);
h2 = line([0 sim_time], [0 0], 'Color', 'r', 'LineStyle', '--', 'LineWidth', 1);
% 添加文字标注
text(sim_time*0.8, 103, '饱和上限', 'Color', 'r', 'FontSize', 9);
text(sim_time*0.8, 3, '下限', 'Color', 'r', 'FontSize', 9);
hold off;
% 子图3：误差变化曲线
subplot(2,2,3);
error_signal = setpoint - y;
plot(t, error_signal, 'm-', 'LineWidth', 1.5);
xlabel('时间 (秒)'); ylabel('误差 (℃)');
title('控制误差随时间变化');
grid on;
xlim([0 sim_time]);
hold on;
line([0 sim_time], [0 0], 'Color', 'k', 'LineStyle', '--');
hold off;
ylim([-5 45]);
% 子图4：PID各分量贡献
subplot(2,2,4);
% 重新计算各分量用于展示
P_component = zeros(1,N);
I_component = zeros(1,N);
D_component = zeros(1,N);
integral_temp = 0;
prev_err = 0;
for k = 2:N
    error_temp = setpoint - y(k-1);
    P_component(k) = Kp * error_temp;
    integral_temp = integral_temp + Ki * error_temp * dt;
    I_component(k) = integral_temp;
    if k > 2
        D_component(k) = -Kd * (y(k-1) - y(k-2)) / dt;
    end
end
plot(t, P_component, 'r-', 'LineWidth', 1); hold on;
plot(t, I_component, 'b-', 'LineWidth', 1);
plot(t, D_component, 'g-', 'LineWidth', 1);
xlabel('时间 (秒)'); ylabel('控制分量');
title('PID各分量贡献');
legend('比例(P)', '积分(I)', '微分(D)', 'Location', 'east');
grid on;
xlim([0 sim_time]);
%% 6. 参数敏感性分析（可选）
figure('Position', [100, 100, 1000, 400]);
% 测试不同Kp的影响
Kp_test = [1.5, 2.5, 4.0];
colors = {'b', 'r', 'g'};
subplot(1,2,1);
hold on;
for i = 1:length(Kp_test)
    % 重新仿真
    [y_test, t_test] = run_pid_simulation(Kp_test(i), Ki, Kd, dt, sim_time);
    plot(t_test, y_test, colors{i}, 'LineWidth', 1.5);
end
plot(t, setpoint*ones(1,N), 'k--', 'LineWidth', 1.5);
xlabel('时间 (秒)'); ylabel('温度 (℃)');
title('不同Kp参数的影响');
legend('Kp=1.5', 'Kp=2.5', 'Kp=4.0', '设定值', 'Location', 'southeast');
grid on;
xlim([0 sim_time]);
% 测试不同Ki的影响
Ki_test = [0.02, 0.08, 0.15];
subplot(1,2,2);
hold on;
for i = 1:length(Ki_test)
    [y_test, t_test] = run_pid_simulation(Kp, Ki_test(i), Kd, dt, sim_time);
    plot(t_test, y_test, colors{i}, 'LineWidth', 1.5);
end
plot(t, setpoint*ones(1,N), 'k--', 'LineWidth', 1.5);
xlabel('时间 (秒)'); ylabel('温度 (℃)');
title('不同Ki参数的影响');
legend('Ki=0.02', 'Ki=0.08', 'Ki=0.15', '设定值', 'Location', 'southeast');
grid on;
xlim([0 sim_time]);
%% 7. 3D参数敏感性分析（可选）
figure('Position', [100, 100, 800, 600]);
% 在Kp和Ki参数空间进行扫描
Kp_range = 1:0.3:4;
Ki_range = 0.02:0.01:0.16;
ISE = zeros(length(Kp_range), length(Ki_range)); % 积分平方误差
for i = 1:length(Kp_range)
    for j = 1:length(Ki_range)
        [y_test, ~] = run_pid_simulation(Kp_range(i), Ki_range(j), Kd, dt, sim_time);
        error_test = setpoint - y_test;
        ISE(i,j) = sum(error_test.^2) * dt;
    end
end
surf(Ki_range, Kp_range, ISE);
xlabel('Ki'); ylabel('Kp'); zlabel('ISE (积分平方误差)');
title('PID参数对性能的影响 (Kd=8.0固定)');
colorbar;
shading interp;
view(45, 30);
fprintf('\n最佳参数点 (最小ISE):\n');
[min_ISE, idx] = min(ISE(:));
[row, col] = ind2sub(size(ISE), idx);
fprintf('Kp = %.1f, Ki = %.3f, ISE = %.2f\n', Kp_range(row), Ki_range(col), min_ISE);

%% 辅助函数：计算调节时间
function ts = find_settling_time(t, y, setpoint, tolerance)
    error_band = setpoint * tolerance;
    settled = abs(y - setpoint) <= error_band;
    % 找到最后离开误差带的时间
    last_exit = 0;
    for i = 2:length(y)
        if settled(i-1) && ~settled(i)
            last_exit = t(i);
        end
    end
    if last_exit == 0
        ts = t(end);
    else
        % 找到最后进入误差带并保持的时间
        idx = find(settled & (t > last_exit), 1);
        if isempty(idx)
            ts = t(end);
        else
            ts = t(idx);
        end
    end
end

%% 辅助函数：运行PID仿真（修正版）
function [y, t] = run_pid_simulation_fixed(Kp, Ki, Kd, dt, sim_time)
    K = 1.2; tau = 30; L = 5;
    N = round(sim_time / dt);
    t = (0:N-1) * dt;
    y = zeros(1, N);
    y(1) = 20;
    y(2) = 20;  % 初始化第二个点，避免索引错误
    setpoint = 60;
    integral = 0;
    prev_error = setpoint - y(1);
    for k = 2:N
        error = setpoint - y(k-1);
        P = Kp * error;
        integral = integral + Ki * error * dt;
        % 微分项（使用测量值微分，处理边界条件）
        if k >= 3
            % 使用前两个测量值计算微分
            D = -Kd * (y(k-1) - y(k-2)) / dt;
        else
            % 第一个时间步，微分项为0
            D = 0;
        end
        u = P + integral + D;
        u = max(0, min(100, u));
        delay_steps = round(L / dt);
        if k > delay_steps
            u_delayed = u;
        else
            u_delayed = 0;
        end
        dT = (K * u_delayed - (y(k-1) - 20)) / tau;
        y(k) = y(k-1) + dT * dt;
        % 防止温度异常（可选）
        if y(k) < 19 || y(k) > 85
            y(k) = y(k-1);  % 如果温度异常，保持原值
        end
        prev_error = error;
    end
end