粒子群算法（PSO）整定 PID 参数MATLAB 程序

粒子群算法（PSO）整定 PID 参数MATLAB 程序

---

一、PID + PSO 整定思路（先对齐）

1、问题定义

PID 控制器：

$u(t) = K_p e(t) + K_i \int e(t)dt + K_d \frac{de(t)}{dt}$

优化目标：寻找一组 $K_p, K_i, K_d$，使系统性能指标最优。

2、性能指标（适应度函数）

常用指标（选一个或加权组合）：

• ISE（Integral Squared Error）：($\int_0^\infty e^2(t)dt$)
• IAE（Integral Absolute Error）：($\int_0^\infty |e(t)|dt$)
• ITAE（Integral Time Absolute Error）：(\int_0^\infty t|e(t)|dt)
• 超调量 + 调节时间（工程最常用）

3、PSO 优化流程

初始化粒子群（随机PID参数）
↓
计算每个粒子的适应度（仿真闭环系统）
↓
更新个体最优 & 全局最优
↓
更新粒子速度与位置
↓
判断是否收敛或达到最大迭代次数
↓
输出最优PID参数

---

二、被控对象模型（以直流电机为例）

1、直流电机传递函数

$G(s) = \frac{K}{Js + B}$

其中：

• ($K$)：电机转矩系数
• ($J$)：转动惯量
• ($B$)：阻尼系数

%% plant_model.m
function sys = plant_model()
% 直流电机模型
J = 0.01;      % 转动惯量 (kg·m²)
B = 0.1;       % 阻尼系数 (N·m·s)
K = 0.01;      % 转矩系数 (N·m/A)

% 传递函数: G(s) = K / (J*s + B)
sys = tf(K, [J B]);
end

---

三、PSO 核心算法

1、粒子结构与初始化

%% pso_pid_init.m
function [particles, velocities] = pso_pid_init(n_particles, param_bounds)
% 初始化PSO粒子群
% n_particles: 粒子数量
% param_bounds: 参数边界 [Kp_min, Ki_min, Kd_min; Kp_max, Ki_max, Kd_max]

n_params = 3;  % Kp, Ki, Kd

% 随机初始化粒子位置（PID参数）
particles = rand(n_particles, n_params) .* ...
            (param_bounds(2,:) - param_bounds(1,:)) + param_bounds(1,:);

% 初始化粒子速度
velocities = zeros(n_particles, n_params);

% 显示初始化信息
fprintf('PSO初始化完成：%d个粒子，参数边界：\n', n_particles);
fprintf('Kp: [%.2f, %.2f]\n', param_bounds(1,1), param_bounds(2,1));
fprintf('Ki: [%.2f, %.2f]\n', param_bounds(1,2), param_bounds(2,2));
fprintf('Kd: [%.2f, %.2f]\n', param_bounds(1,3), param_bounds(2,3));
end

2、适应度函数（闭环仿真）

%% pid_fitness.m
function fitness = pid_fitness(pid_params, plant, sim_time, dt)
% 计算PID控制器的适应度
% pid_params: [Kp, Ki, Kd]
% plant: 被控对象传递函数
% sim_time: 仿真时间
% dt: 仿真步长

Kp = pid_params(1);
Ki = pid_params(2);
Kd = pid_params(3);

% 创建PID控制器
C = pid(Kp, Ki, Kd);

% 闭环系统
sys_cl = feedback(C * plant, 1);

% 仿真时间向量
t = 0:dt:sim_time;

% 阶跃响应
[y, t] = step(sys_cl, t);

% 计算误差
r = ones(size(y));  % 参考输入（单位阶跃）
e = r - y;

% ---------- 性能指标选择 ----------
% 1. ISE (Integral Squared Error)
fitness = sum(e.^2) * dt;

% 2. IAE (Integral Absolute Error)
% fitness = sum(abs(e)) * dt;

% 3. ITAE (Integral Time Absolute Error)
% time_weight = t';
% fitness = sum(time_weight .* abs(e)) * dt;

% 4. 工程指标：超调量 + 调节时间
% [y_peak, peak_idx] = max(y);
% overshoot = (y_peak - 1) / 1 * 100;  % 超调量(%)
% settling_idx = find(abs(y - 1) > 0.02, 1, 'last');
% settling_time = t(settling_idx);
% fitness = overshoot + 10*settling_time;  % 加权组合

% 惩罚不稳定系统
if any(isnan(y)) || any(isinf(y))
    fitness = 1e6;
end

% 惩罚过大控制量
if Kp > 100 || Ki > 100 || Kd > 100
    fitness = fitness + 1e3;
end
end

3、PSO 主算法

%% pso_pid_main.m
function [best_params, best_fitness, convergence_curve] = pso_pid_main()
% PSO主程序：PID参数自整定

%% ========== 参数设置 ==========
n_particles = 30;       % 粒子数量
max_iter = 50;          % 最大迭代次数
w_max = 0.9; w_min = 0.4; % 惯性权重范围
c1 = 2.0; c2 = 2.0;     % 学习因子

% PID参数边界 [Kp, Ki, Kd]
param_bounds = [0.1, 0.01, 0.001;   % 最小值
                100, 10,   1];      % 最大值

sim_time = 5;           % 仿真时间(s)
dt = 0.01;              % 仿真步长(s)

%% ========== 获取被控对象 ==========
plant = plant_model();

%% ========== 初始化粒子群 ==========
[particles, velocities] = pso_pid_init(n_particles, param_bounds);

% 个体最优和全局最优
pbest = particles;
pbest_fitness = inf(n_particles, 1);
gbest = zeros(1, 3);
gbest_fitness = inf;

convergence_curve = zeros(max_iter, 1);

%% ========== PSO主循环 ==========
for iter = 1:max_iter
    fprintf('迭代 %d/%d:\n', iter, max_iter);
    
    % 动态惯性权重
    w = w_max - (w_max - w_min) * iter / max_iter;
    
    for i = 1:n_particles
        % 计算适应度
        fitness = pid_fitness(particles(i,:), plant, sim_time, dt);
        
        % 更新个体最优
        if fitness < pbest_fitness(i)
            pbest_fitness(i) = fitness;
            pbest(i,:) = particles(i,:);
        end
        
        % 更新全局最优
        if fitness < gbest_fitness
            gbest_fitness = fitness;
            gbest = particles(i,:);
            fprintf('  新最优: Kp=%.3f, Ki=%.3f, Kd=%.3f, 适应度=%.6f\n', ...
                    gbest(1), gbest(2), gbest(3), gbest_fitness);
        end
    end
    
    % 更新粒子速度和位置
    for i = 1:n_particles
        % 速度更新
        r1 = rand(1,3);
        r2 = rand(1,3);
        velocities(i,:) = w * velocities(i,:) + ...
                         c1 * r1 .* (pbest(i,:) - particles(i,:)) + ...
                         c2 * r2 .* (gbest - particles(i,:));
        
        % 速度限幅
        max_vel = (param_bounds(2,:) - param_bounds(1,:)) * 0.2;
        velocities(i,:) = max(min(velocities(i,:), max_vel), -max_vel);
        
        % 位置更新
        particles(i,:) = particles(i,:) + velocities(i,:);
        
        % 位置限幅
        particles(i,:) = max(min(particles(i,:), param_bounds(2,:)), param_bounds(1,:));
    end
    
    convergence_curve(iter) = gbest_fitness;
    
    % 早停条件
    if iter > 10 && std(convergence_curve(iter-9:iter)) < 1e-6
        fprintf('收敛稳定，提前结束迭代\n');
        convergence_curve(iter+1:max_iter) = gbest_fitness;
        break;
    end
end

best_params = gbest;
best_fitness = gbest_fitness;

fprintf('\n========== 优化结果 ==========\n');
fprintf('最优PID参数:\n');
fprintf('  Kp = %.4f\n', best_params(1));
fprintf('  Ki = %.4f\n', best_params(2));
fprintf('  Kd = %.4f\n', best_params(3));
fprintf('最优适应度: %.6f\n', best_fitness);
end

---

四、结果验证与可视化

%% validate_pid_performance.m
function validate_pid_performance(best_params, plant)
% 验证优化后的PID性能

Kp = best_params(1);
Ki = best_params(2);
Kd = best_params(3);

% 创建PID控制器
C = pid(Kp, Ki, Kd);

% 闭环系统
sys_cl = feedback(C * plant, 1);

% 仿真时间
t = 0:0.01:5;

% 阶跃响应
figure('Color','white','Position',[100 100 1200 400])

% 1. 阶跃响应
subplot(1,3,1)
y = step(sys_cl, t);
plot(t, y, 'b-', 'LineWidth', 2)
hold on
plot(t, ones(size(t)), 'r--', 'LineWidth', 1)
xlabel('Time (s)')
ylabel('Amplitude')
title('Step Response')
grid on
legend('Closed Loop', 'Reference')

% 计算性能指标
[y_peak, peak_idx] = max(y);
overshoot = (y_peak - 1) / 1 * 100;
settling_idx = find(abs(y - 1) > 0.02, 1, 'last');
settling_time = t(settling_idx);
steady_state_error = abs(1 - y(end));

fprintf('性能指标:\n');
fprintf('  超调量: %.2f%%\n', overshoot);
fprintf('  调节时间: %.3f s\n', settling_time);
fprintf('  稳态误差: %.4f\n', steady_state_error);

% 2. 控制量
subplot(1,3,2)
u = lsim(C, 1-y, t);
plot(t, u, 'g-', 'LineWidth', 2)
xlabel('Time (s)')
ylabel('Control Effort')
title('Control Signal')
grid on

% 3. 误差响应
subplot(1,3,3)
e = 1 - y;
plot(t, e, 'r-', 'LineWidth', 2)
xlabel('Time (s)')
ylabel('Error')
title('Error Response')
grid on
end

---

五、主程序调用

%% main_pid_tuning.m
clear; clc; close all;

% 运行PSO优化
[best_params, best_fitness, convergence_curve] = pso_pid_main();

% 获取被控对象
plant = plant_model();

% 验证性能
validate_pid_performance(best_params, plant);

% 绘制收敛曲线
figure('Color','white')
plot(convergence_curve, 'b-', 'LineWidth', 2)
xlabel('Iteration')
ylabel('Best Fitness')
title('PSO Convergence Curve')
grid on

% 与传统方法对比（Ziegler-Nichols）
fprintf('\n========== 与传统方法对比 ==========\n');
[Ku, Tu] = get_ultimate_gain(plant);
Kp_zn = 0.6 * Ku;
Ki_zn = 1.2 * Ku / Tu;
Kd_zn = 0.075 * Ku * Tu;
fprintf('Ziegler-Nichols PID: Kp=%.3f, Ki=%.3f, Kd=%.3f\n', Kp_zn, Ki_zn, Kd_zn);

% 验证ZN方法性能
C_zn = pid(Kp_zn, Ki_zn, Kd_zn);
sys_zn = feedback(C_zn * plant, 1);
step(sys_zn)
hold on
step(feedback(pid(best_params(1), best_params(2), best_params(3)) * plant, 1))
legend('Ziegler-Nichols', 'PSO Optimized')
title('Performance Comparison')

参考代码 通过粒子群算法实现pid参数的自整定 www.youwenfan.com/contentcsw/82832.html

六、工程级改进方案

1、多目标优化（NSGA-II）

% 同时优化多个指标
fitness = w1*overshoot + w2*settling_time + w3*control_effort;

2、约束处理

% 添加约束：控制量饱和、超调限制
if overshoot > 20%
    fitness = fitness + 1e4;  % 惩罚
end

3、自适应PSO

% 根据收敛情况调整参数
if convergence_rate < threshold
    c1 = c1 * 1.1;  % 增强局部搜索
    c2 = c2 * 0.9;  % 减弱全局搜索
end

4、离散PSO（数字控制器）

% 考虑采样时间和量化效应
dt = 0.001;  % 更小的采样时间
% 离散PID实现
uk = Kp*ek + Ki*sum(ek)*dt + Kd*(ek - ek_1)/dt;

---

七、与现有方法的对比

方法	优点	缺点	适用场景
Ziegler-Nichols	简单快速	保守，超调大	初步调试
Cohen-Coon	比ZN稍好	仍较保守	过程控制
PSO优化	全局最优，性能优异	计算量大	离线整定
GA优化	全局搜索能力强	收敛慢	复杂系统
强化学习	在线自适应	实现复杂	智能控制