粒子群算法整定PID参数：MATLAB/Simulink双闭环实践

简介：这份PDF文档讲解如何用粒子群算法（PSO）在Simulink中整定PID参数，面向需要减少Kp、Ki、Kd手工试凑的高校学生和工程技术人员。内容涵盖PSO的参数范围设定、速度与位置更新、适应度评估及pBest、gBest更新，并对比基于M文件的参数整定与复杂系统仿真优化，详解以PSO自动寻优PID的过程，以及c1、c2和惯性权重w对探索与开发平衡的影响。读者可掌握按超调量、调节时间等指标设计评价函数的调参思路，多输入多输出系统同样适用。压缩包共1个PDF文件，约529KB，篇幅紧凑，适合算法入门与仿真实践对照；文档另附运行顺序、迭代观察和结果出图提示。已有4895人学习。

1. 从手工试凑到四千次仿真：我把 PID 整定交给粒子群的理由

改 PID 的痛苦做控制的人都懂：三个参数、六个调整方向，动一个另外两个就得重来，查完表还得凭手感修，换一套工况又得从头再试一遍。这套基于 Simulink 的粒子群算法（PSO）整定 PID 参数资源，核心思路是把这件事外包给算力——你写一个评价函数，明确告诉算法"超调要小、调节时间要短"，剩下由几百只"鸟"在参数框里替你搜索。资源里给的是 MATLAB M 文件加一个 slx 仿真模型，从set_para.m起步，跑完吐三张 figure。它适合手上已有可跑通的 Simulink 被控对象、指标能量化成超调和调节时间的人；不适合连模型都还没搭起来的纯新手。

2. 粒子群算法落地：速度更新式里的五个参数怎么定

2.1 为什么整定 PID 用 PSO 而不是查表加手调

传统工程整定法（比如直流调速里的典型 I 型、典型 II 型系统设计）本质上是把被控对象近似成一阶或二阶环节，再套公式算出一组 PI 参数。这条路的前提是模型足够"干净"——纯延迟小、参数不漂、结构固定。一旦对象里有饱和、有摩擦、有未建模的动态，算出来的参数就只是起点，后面还得靠人肉试凑。

PSO 和遗传算法同属群体迭代优化，区别在于它没有交叉、变异那套操作，需要整定的超参数更少：惯性权重 w、两个学习因子 c1 和 c2，就这三个。少了算子组合的玄学，收敛也更快。更实际的一点是它的"所需即所得"——你不需要知道对象内部原理，只要能把"好不好"量化成一个标量，算法就朝这个方向压。多输入多输出系统里想让 A 通道超调优先、B 通道调节时间优先，改评价函数里的权重就行，不用重新推公式。

代价当然是算力。一次 200 粒子 × 20 代的搜索就是 4000 次仿真，纯 M 文件的 PID 引擎还能扛，挂 Simulink 模型就得按分钟算了。所以常见做法是先用轻量 M 文件版本把参数框和评价函数调通，再切到 Simulink 去整复杂对象。

2.2 速度—位置更新：一份能直接跑的 MATLAB 实现

PSO 的核心就两行公式：速度由"上一步速度的惯性 + 奔向个体历史最优 + 奔向全局历史最优"三部分合成，位置是位置加速度乘步长。写成代码时最容易出问题的是限幅和向量化——限幅写错粒子会飞出参数框，向量化写错就得循环几百次调用仿真，慢一个数量级。

下面这份实现把整个种群拼成矩阵一次算完，是后面所有整定工作的骨架：

function [Pg,fmin,hist] = PSO_pid(dim,n,m,xmin,xmax,vmax) % dim 参数维度（整定 Kp/Ki/Kd 时为 3） % n 迭代代数；m 种群规模 % xmin/xmax 1×dim 行向量，逐维限幅；vmax 标量速度上限 w = 1; c1 = 2; c2 = 2; dt = 0.3; % 惯性权重、两个学习因子、位移步长 hist = zeros(n,3); % 每代记录 [最优适应度 ts sigma] rng(1); % 固定随机种子，否则每次结果都不一样 v = zeros(m,dim); % 初速度全零 X = (xmax-xmin).*rand(m,dim) + xmin; % 初始位置在参数框内均匀撒点 P = X; % 个体历史最优位置 last_f = PID_access(X); % 整个种群一次算完，向量化关键 [fmin,idx] = min(last_f); Pg = X(idx,:); % 全局最优位置 for k = 1:n v = w*v + c1*rand*(P-X) + c2*rand*(ones(m,1)*Pg - X); v = min(max(v,-vmax),vmax); % 速度限幅 X = X + v*dt; % 位置更新，dt 决定搜索步长 X = min(max(X,xmin),xmax); % 逐维位置限幅 new_f = PID_access(X); better = new_f < last_f; P(better,:) = X(better,:); % 只覆盖变好的粒子 [new_fmin,idx] = min(new_f); if new_fmin < fmin fmin = new_fmin; Pg = X(idx,:); end last_f = new_f; [~,~,ts,sigma] = PID_sim(Pg(1),Pg(2),Pg(3),false); hist(k,:) = [fmin,ts,sigma]; fprintf('第 %d/%d 代 适应度 %.6f Kp=%.4f Ki=%.4f Kd=%.4f\n', ... k,n,fmin,Pg(1),Pg(2),Pg(3)); w = w - k/n*0.7*w; % 惯性权重随代数线性衰减 end end

逻辑上分四段看：初始化把 m 个粒子撒在[xmin,xmax]里，速度清零；速度更新那行三部分对应惯性、认知、社会；位置更新后立刻做两层限幅，先限速度再限位置，顺序反了会把粒子推出边界；最后比较适应度，只有变好的粒子才更新个体最优，全局最优则用一次比较搞定，避免用find遍历全种群。

参数怎么定，我的经验是：

rand这里取的是标量随机数，整个种群共用同一组系数，这是原实现的写法。想让粒子间独立性更强，可以换成rand(m,1)逐粒子的随机数，收敛曲线会稍平滑一点，但差距不大，不值得为这个纠结。

2.3 惯性权重衰减：收敛慢时第一个要动的开关

如果只跑不做衰减，你会看到典型的翻车现场——粒子在极值点附近来回震荡，一部分粒子甚至有明显的反向运动，几十代下来最优值几乎不动。原因在于位置更新完全依赖速度，而惯性项让速度持续累积，到了极值附近本该减速，实际却在"过冲"。

衰减式w = w - k/n*0.7*w的含义是：第 k 代时惯性权重降到初值的 0.3 倍左右。前期 w 接近 1，粒子敢往远处跑；后期 w 变小，速度被压住，搜索变平滑。实测下来同一组算力，加衰减后最终适应度往往能压下去一截，代价是前期探索略收敛。

需要留意的是w衰减和vmax会互相干扰。如果vmax定得很小，衰减的效果就被掩盖了；反之vmax很大而w不衰减，粒子会长时间在边界上蹭。判断标准很简单：看位置分布图，粒子应该在前几代散得开、后面几代越聚越紧，一路不收或者一路不散都说明这两个参数没配好。

3. 适应度函数才是方向盘：把超调和调节时间压成一个标量

3.1 为什么评价函数要取对数而不是简单加权求和

整定 PID 时通常关心四个指标：超调量 σ、上升时间、峰值时间、调节时间 ts。指标越多，加权求和越难调——超调量是百分比（0～100 量级），调节时间是秒（0.01～1 量级），直接相加等于让超调量单方面主导。

资源里给的思路是对数压缩：

$$f_{access}=\ln\left(\frac{t_s}{5\times10^{-2}}+1\right)+\ln\left(\frac{\sigma}{10^{-4}}+1\right)$$

这个形式有两个好处。一是量纲被拉平，两个指标变成可比的无量纲量；二是导数在期望值附近的斜率差异被保留下来——对 ts 求导在 ts=0.05 处约 2.5×10⁻²，对 σ 求导在 σ=10⁻⁴ 处约 0.5×10⁻⁴，也就是说 ts 每偏一点带来的适应度变化，远大于 σ 同样比例偏移带来的变化。同时两个函数都加了 1，参数趋近 0 时导数为 1，不会发散。

更关键的是不对称性：指标一旦小于期望值，函数继续下降但导数越来越大，最多到 1；一旦大于期望值，函数上升但速率越来越慢。结果是算法天然倾向于"宁可多满足，不要勉强凑"，多个指标里只要有一个严重超标，代价就顶不住。

3.2 PID 引擎与指标提取：三支 M 文件的分工

评价函数不能凭空算，得先有一份能对参数向量批量仿真的引擎。这三个文件是一条链：PID_sim.m做仿真出曲线，parameter_cal.m从曲线抠指标，PID_access.m把指标合成一个数。

function [f_infty,tp,ts,sigma] = PID_sim(kp,ki,kd,debug) % kp/ki/kd 支持向量输入，一次并行算 M 组参数 Tt = 5e-4; % 仿真步长 T0 = 1e-2; % 控制周期（离散 PID 的采样时间） Tf = 1e-3; % 微分项一阶滤波时间常数 alp = Tf/(Tt+Tf); % 离散化后的滤波系数 total_t = 1; % 单次仿真时长，太短则曲线进不了稳态 N = floor(total_t/Tt); M = numel(kp); kp = reshape(kp,M,1); ki = reshape(ki,M,1); kd = reshape(kd,M,1); sys = tf(0.998,[0.021,1]); % 被控对象，一阶惯性近似直流电机转速 dsys = c2d(sys,Tt,'z'); % 离散化 [num,den] = tfdata(dsys,'v'); e_1 = zeros(M,1); Ee = zeros(M,1); u_1 = zeros(M,1); y_1 = zeros(M,1); ud_1 = zeros(M,1); r = 1500; % 阶跃给定 y = zeros(M,N); u = zeros(M,N); for k = 1:N y(:,k) = -den(2).*y_1 + num(2)*u_1 + num(1).*u(:,k); % 差分方程输出 e0 = r - y(:,k); ud = alp.*ud_1 + (1-alp)/T0.*kd.*(e0-e_1); % 微分项带滤波 u(:,k) = kp.*e0 + T0*ki.*Ee + ud; % 位置式 PID Ee = Ee + e0; u_1 = u(:,k); y_1 = y(:,k); e_1 = e0; ud_1 = ud; end if debug, plot(linspace(0,total_t,N),y), end [f_infty,tp,ts,sigma] = parameter_cal(y,Tt,2e-2,debug); end

这段代码的关键在M = numel(kp)和几个zeros(M,1)：所有状态量都堆成列向量，循环一次就把 m 组参数的响应同时算完。这就是为什么 MATLAB 版本不需要额外写并行——矩阵运算本身就吃了向量化红利。如果改成逐组循环调用仿真，200 个粒子光调度开销就够呛。

parameter_cal.m负责从响应矩阵里抠指标，思路是倒序找最后一个超出误差带的点：

function [f_infty,tp,ts,sigma] = parameter_cal(y,T0,delt_err,debug) M = size(y,1); N = size(y,2); f_infty = y(:,N); % 用末点近似稳态值 err = y - f_infty*ones(1,N); ferr = fliplr(abs(err)); % 倒序取绝对值 [~,ts_i] = max(ferr > delt_err*f_infty,[],2); ts_i = N*ones(M,1) - ts_i; ts = ts_i*T0; % 调节时间 [fp,tp] = max(y,[],2); tp = tp*T0; tp(abs(fp-f_infty) < 1e-5) = NaN; % 过阻尼无峰值，置 NaN 防污染 sigma = (fp - f_infty)./f_infty; % 超调量 end

delt_err = 2e-2就是 ±2% 的稳态误差带，这是控制教材里的常规取法，工程上按需改成 5% 也行，改完 ts 会整体变小，注意别和上一次的结果直接比。tp置 NaN 那行是必须的，过阻尼系统没有峰值，不处理的话max会返回末点，超调量算出来是负的，评价函数直接乱掉。

最后PID_access.m只有三行：

function y = PID_access(para) kp = para(:,1); ki = para(:,2); kd = para(:,3); [~,~,ts,sigma] = PID_sim(kp,ki,kd,false); y = log(ts/5e-2 + 1) + log(sigma/1e-2 + 1); % 越小的参数越优 end

注意这里 σ 的分母写的是1e-2，而正文公式里是1e-4。这不是笔误，而是权重：分母越小，同样的超调量算出来的对数项越大，惩罚越重。示例代码把超调的权重放宽了两个数量级，整出来的参数会更偏向压调节时间。想让超调优先，把分母改回1e-4甚至更小即可。这一处差异在很多复现里被当成 bug 改掉，其实改的是优化目标本身。

3.3 权重怎么改：让不同通道的指标各让一步

改成 Simulink 版本后指标变成四个：转速超调 σn、转速调节时间 tsn、电流超调 σi、电流调节时间 tsi。四项直接相加，谁的量级大谁说话，所以要先把每项除以自己的期望值：

J = log(ts_n/5e-2 + 1) + log(sig_n/1e-4 + 1) ... + log(ts_i/1e-2 + 1) + log(sig_i/5e-2 + 1);

这里的四个分母就是四个"价值权重"。分母调小 = 该指标被重视；分母调大 = 该指标可以妥协。比如启动过程更在意电流不冲过载，就把sig_i的分母压到1e-3；如果负载扰动恢复速度更重要，就把ts_n的分母压到1e-2。

我一般会先跑一版四个分母都取工程经验值的结果，看曲线哪一项明显不满意，再单独压那一项的分母重跑，而不是一次性把所有分母都调小——目标太苛刻会让整个种群挤在边界上，适应度反而降不下去。

4. 接进 Simulink 双闭环：用 global 把 M 文件和 slx 焊在一起

4.1 双闭环模型里哪些环节必须变成变量

直流调速系统的结构是外环转速、内环电流，两个调节器都是 PI。转速环负责抗负载扰动、抑制参数漂移，电流环负责让电枢电流快速跟随给定。启动时电流先顶到过载值形成平台，电机以最大加速度升速，接近给定转速时电流迅速回落，最后稳定在负载电流上。

模型参数按下面这组搭（资源里给的示例模型double_circle.slx）：

按工程设计方法能算出一组基准参数：[Knp, KnI, Kip, KiI] = [7.3660, 1.5658, 0.8257, 30.2773]，跑出来的指标是 σn ≈ 4.07×10⁻⁴%、tsn ≈ 0.7875 s、σi ≈ 20.38%、tsi ≈ 0.0172 s。这组数就是 PSO 的对照组——如果整定结果连这组都打不过，说明评价函数或者参数框出了问题，不用往下看了。

要让 PSO 动态改这四个参数，模型里的四个比例/积分增益模块不能填常数，必须填变量名（Kp_n、Ki_n、Kp_i、Ki_i），并且在 M 文件里把它们声明成global。变量进了工作区，Simulink 仿真时才能读到新值。改完另存为PSO_double_circle.slx，原模型留着做对照。

4.2 global 与 assignin：M 文件驱动 slx 的两种接法

这里有个容易翻车的细节：global声明的是"当前函数工作区"的全局变量，而 Simulink 模型默认从 base 工作区取变量。如果 PSO 是在函数里跑的，只在函数内global Kp_n是不够的，模型根本看不到。两种稳妥接法：

一是全程global，并且保证主脚本也在 base 工作区跑：

global Kp_n Ki_n Kp_i Ki_i Kp_n = 7.3660; Ki_n = 1.5658; Kp_i = 0.8257; Ki_i = 30.2773;

二是用assignin显式往 base 工作区写，函数里也能用：

assignin('base','Kp_n',X(j,1)); assignin('base','Ki_n',X(j,2)); assignin('base','Kp_i',X(j,3)); assignin('base','Ki_i',X(j,4));

我一般选第二种，原因很直接：调sim的时候如果把仿真模式设成'accelerator'或'rapid accelerator'，模型编译后会缓存参数，改工作区变量不一定立刻生效；用assignin至少能保证每次写进去的值是明确的，出问题时好排查。

驱动一次仿真就是一行：

out = sim('PSO_double_circle.slx', ... 'StopTime','1.5', ... 'SimulationMode','normal', ... 'ReturnWorkspaceOutputs','on');

StopTime要覆盖到转速完全稳定的时刻，1.5 s 对这套模型够用；SimulationMode建议先'normal'跑通，确认参数真的在变之后再换加速模式。加速模式第一次编译慢，但连续几千次仿真能省下大量时间，这是后面能不能把规模跑起来的关键。

4.3 cal_n.m：从 logsout 里抠出四个指标

模型输出用To Workspace模块记成 Dataset 格式，取出来是 timeseries。指标计算和 M 文件版本逻辑一致，只是数据源换了：

function [sig_n,ts_n,sig_i,ts_i] = cal_n(logsout) n_ts = logsout.get('n').Values; % 转速 timeseries i_ts = logsout.get('i').Values; % 电枢电流 timeseries n_end = n_ts.Data(end); sig_n = (max(n_ts.Data) - n_end)/abs(n_end); % 转速超调 idx = find(abs(n_ts.Data - n_end) > 0.02*abs(n_end),1,'last'); ts_n = n_ts.Time(idx); % 转速调节时间 i_plat = mean(i_ts.Data(round(0.6*end):round(0.8*end))); % 电流平台均值 sig_i = (max(i_ts.Data) - i_plat)/i_plat; % 电流超调 idx = find(abs(i_ts.Data - i_plat) > 0.05*i_plat,1,'last'); ts_i = i_ts.Time(idx); % 电流调节时间 end

电流环和转速环的稳态值取法不一样，这是有意的：转速稳态是末点值，电流稳态是平台段均值。原因在于电流响应在启动瞬间就到顶了，后半段虽然有微小纹波，但均值才是真正的稳态；如果也取末点，一点纹波就能让sig_i上下跳几个百分点。

0.6*end到0.8*end这个窗口是经验值，取的是电流已经回落稳定的区间。如果仿真时长改了，这个窗口要跟着调，否则会把进入稳态前的过渡段算进平台均值里。

4.4 一次完整整定的调用链和耗时预期

把四维参数的整定串起来，主入口仍然是一个脚本：

dimension = 4; % [Kp_n Ki_n Kp_i Ki_i] n = 15; m = 30; % 先小规模试跑 xmin = [1, 0.1, 0.1, 5]; xmax = [20, 10, 5, 80]; vmax = [2, 1, 0.5, 5]; % 四维量级差很大，速度上限要逐维给 [Pg,fmin,hist] = PSO_double_circle(dimension,n,m,xmin,xmax,vmax);

注意vmax在这里从标量变成了逐维向量。四个参数的量级差两个数量级（Kp 在 1～20，KiI 在 5～80），用统一速度上限会让小量级维度几乎不动。这是切到多参数系统后最容易忽略的一处改动。

耗时上心里要有数：单次 Simulink 仿真按 0.05～0.2 s 算，m=30、n=15就是 450 次，五到十秒能出结果；真要按m=200、n=20跑满 4000 次，来回就是十几分钟到半小时。所以我的习惯永远是先 30×15 跑通链路、看适应度曲线有没有下降趋势，再决定要不要加规模。

5. 整定过程避坑：五类翻车现场的现象、原因与处置

5.1 环境与版本类

现象：脚本一跑就报Undefined function or variable，或者tfdata返回的格式和预期不一样，取den(2)时报索引越界。

原因：这套代码写在 R2017a 上，之后的版本里tfdata默认返回元胞数组的行为、sim的调用语法、隐式扩展的支持范围都有变化。另外从别处拷过来的文件如果和 slx 不在同一目录，也没addpath，同样报找不到函数。

解决：把全部 m 文件和 slx 放进同一个文件夹，用set_para.m作为唯一入口启动，运行前在命令行cd到该目录；tfdata(dsys,'v')一定要带'v'参数强制返回行向量；sim用 Name-Value 形式传参，别用老的simset配置集。真跑不通就明确降到 R2017a 验证一次，确认是版本问题再逐条改语法。

现象：仿真结果里转速曲线是一条直线，或者电流一直是 0。

原因：模型里的增益模块还填着常数，PSO 改的工作区变量压根没被读到。

解决：双击每个比例/积分增益模块看参数框，必须是Kp_n这样的变量名而不是数字；把模型里Model Properties的初始化回调清空，避免回调在每次仿真前又把变量覆盖回默认值。

5.2 算法与仿真类

现象：程序跑了几分钟没有任何输出，怀疑死机。

原因：大概率不是死机，是 figure 刷屏或者单次仿真时间太长。原实现每次迭代都会plot并且hold on，200 条曲线叠上去，绘图开销能占掉大半时间。

解决：先看命令行打印的迭代计数有没有在变——只要在变就是活的。调试阶段把debug参数置false，关掉绘图只留文字输出；确认参数正确后再打开画图。另外把figure提到循环外，或者在循环里只存数据、结束后统一绘制。

现象：最优参数一直贴在xmin或xmax边界上，曲线反而越来越差。

原因：速度限幅和位置限幅不匹配。vmax定得相对dt太大，粒子一步就跨出参数框，被裁到边界后速度没清零，下一代继续往外冲，实际搜索空间只剩边界那一圈。

解决：按vmax * dt ≈ (xmax - xmin)/4来配，或者干脆把dt从 0.3 降到 0.1，让步长细一点。改完看位置分布图，粒子应该在前几代覆盖大半个框，后面才收缩；一直贴着边就是还没配好。

现象：适应度一直在降，但导出的响应曲线看不出改善，甚至超调变大了。

原因：指标提取算错了。最常见的是稳态值用末点近似，而total_t太短曲线根本没进稳态，算出来的ts是 0 或者极大值；其次是过阻尼系统tp没置 NaN，超调量算成负数。

解决：先把total_t加长到曲线明显走平，再检查parameter_cal的分支；单组参数跑一次debug=true，把四个指标打印出来和肉眼看图的结果对一遍，对不上就别急着跑优化。这一步花五分钟，能省掉后面半小时的无效搜索。

现象：同样的脚本跑两次，得到的最优参数完全不同。

原因：rand没设种子，初始种群每次都不一样。PSO 本身是随机算法，小规模种群下这种差异会很明显。

解决：跑之前rng(1)固定种子；同时把xmin/xmax收紧到工程上有意义的范围，别图省事给个[0, 100]——搜索空间越大，随机性带来的方差越大。

6. 从能跑到敢用：结果复核与三处可以继续改进的地方

跑出Pg只是开始，我判断一组参数能不能用，固定走三步：一是把Pg代回模型跑一次纯 Simulink 仿真，不看 PSO 的内部图，只看阶跃响应；二是和工程设计法的基准参数做同条件对比，四项指标列成表；三是把负载突加、给定阶跃两种工况都跑一遍，只看一种工况很容易过拟合。

对比表长这样，PSO 列的数值仅作量级参考，你换模型重跑肯定不一样：

三处值得动手改的地方。第一是给 σ 加保护：对数项在 σ→0 时导数为 1，如果某一代恰好整出零超调，这一项的梯度就消失了，算法会失去方向。改成log(max(sigma,1e-6)/w + 1)或者直接对超调加一个分段罚函数，能避免这种退化。

第二是把惯性权重换成非线性衰减，比如w = w_min + (w_max-w_min)*(1 - k/n)^2，前期降得慢、后期降得快，比线性衰减更适合这种"前松后紧"的参数搜索。

第三是别把 PSO 的结果当终值。它擅长在大范围里找盆地，不擅长在盆地底部精修。我的做法是 PSO 跑完后，把Pg作为初值丢给fminsearch之类的小范围搜索再走一轮，通常还能再压一点适应度。

还有一个血泪经验：先确认 Simulink 模型本身在基准参数下能跑出合理曲线，再交给 PSO。模型有问题的时候，PSO 会用一组离谱参数"适应"出一个好看的指标——它只是在拟合你的评价函数，不是在理解你的电机。从那以后我每次开整定前，都强制先用基准参数跑一遍单步仿真，把转速和电流曲线打印出来看一眼，确认曲线形状对得上物理直觉再放行。希望帮到你。

本文还有配套的精品资源，点击获取