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

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

1、问题定义

PID 控制器：

u(t)=Kpe(t)+Ki∫e(t)dt+Kdde(t)dtu(t) = K_p e(t) + K_i \int e(t)dt + K_d \frac{de(t)}{dt}u(t)=Kp​e(t)+Ki​∫e(t)dt+Kd​dtde(t)​

优化目标：寻找一组Kp,Ki,KdK_p, K_i, K_dKp​,Ki​,Kd​，使系统性能指标最优。

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

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

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

3、PSO 优化流程

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

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

1、直流电机传递函数

G(s)=KJs+BG(s) = \frac{K}{Js + B}G(s)=Js+BK​

其中：

• (KKK)：电机转矩系数
• (JJJ)：转动惯量
• (BBB)：阻尼系数

%% plant_model.mfunctionsys=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.mfunction[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.mfunctionfitness=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; % 加权组合% 惩罚不稳定系统ifany(isnan(y))||any(isinf(y))fitness=1e6;end% 惩罚过大控制量ifKp>100||Ki>100||Kd>100fitness=fitness+1e3;endend

3、PSO 主算法

%% pso_pid_main.mfunction[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主循环 ==========foriter=1:max_iterfprintf('迭代 %d/%d:\n',iter,max_iter);% 动态惯性权重w=w_max-(w_max-w_min)*iter/max_iter;fori=1:n_particles% 计算适应度fitness=pid_fitness(particles(i,:),plant,sim_time,dt);% 更新个体最优iffitness<pbest_fitness(i)pbest_fitness(i)=fitness;pbest(i,:)=particles(i,:);end% 更新全局最优iffitness<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);endend% 更新粒子速度和位置fori=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,:));endconvergence_curve(iter)=gbest_fitness;% 早停条件ifiter>10&&std(convergence_curve(iter-9:iter))<1e-6fprintf('收敛稳定，提前结束迭代\n');convergence_curve(iter+1:max_iter)=gbest_fitness;break;endendbest_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.mfunctionvalidate_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',[1001001200400])% 1. 阶跃响应subplot(1,3,1)y=step(sys_cl,t);plot(t,y,'b-','LineWidth',2)hold onplot(t,ones(size(t)),'r--','LineWidth',1)xlabel('Time (s)')ylabel('Amplitude')title('Step Response')grid onlegend('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 onend

五、主程序调用

%% main_pid_tuning.mclear;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 onstep(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、约束处理

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

3、自适应PSO

% 根据收敛情况调整参数ifconvergence_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;

七、与现有方法的对比