基于PSO的PID参数自整定与MATLAB/Simulink仿真实践

简介：面向需要在Simulink中快速完成PID参数整定的控制工程学习者与自动化方向学生，这份PDF文档围绕粒子群算法（PSO）整定PID参数展开，把智能优化思路与工程调参实践串成一条完整线索。正文先讲PSO的算法原理与迭代步骤，包括速度与位置更新、pBest与gBest的维护、适应度评价与参数权重c1、c2、w的取舍，再进入PID自整定部分，对比基于M文件的整定流程与Simulink仿真环境下的PSO优化实现，并说明如何通过修改评价函数来约束超调量、调节时间等指标。文件包共1个PDF文件，约529KB，篇幅紧凑，便于随查随用；已有4895人学习下载。读者可据此理解PSO与Simulink联合仿真的整体框架，掌握从参数范围设定、种群与迭代次数配置到结果图判读的完整流程，并借鉴文中的评价函数设计思路，将“所需即所得”的调参目标落到自身系统上。

1. 调参调到怀疑人生，换个思路让算法替你撞

调 PID 参数这件事干过的都清楚有多玄学：Kp 往上抬一点超调就炸，往回收一点响应又木得像个木头，三个参数六个方向，手工试凑到最后常常是「能跑就行」。这份资源包换了个思路——不再靠人盯曲线，而是让粒子群算法（PSO）自己去撞参数空间，拿 Simulink 或 M 文件搭出来的被控对象当考场，一轮一轮打分，收敛出来的 Kp/Ki/Kd 再丢回模型里验证。它适合两类人：手上有现成 Simulink 被控模型、参数怎么调都别扭的；以及想搞明白智能优化算法怎么和控制系统仿真接起来的。整套东西不依赖工具箱黑盒，PSO 是自己写的 M 文件，评价函数就是几行 log，想改优化方向随时能动手。

2. 项目骨架与运行链路：从 set_para.m 出门，绕一圈回到三张 figure

拿到压缩包先别急着点运行，把目录翻一遍比什么都省事。这套代码的组织方式很典型：一个入口脚本负责填参数，一个优化器负责迭代，一个仿真引擎负责打分，指标计算和适应度函数各自独立成文件。搞清楚谁调谁，后面改任何一处都不会踩空。

2.1 资源包里每个文件管什么

按职责拆开看，文件大致分成四层，对应「配置 → 优化 → 仿真 → 评价」四个动作：

看这张表能发现一个设计取向：把「仿真」和「评价」彻底分开了。PID_sim.m 只负责吐出一条响应曲线，超调和调节时间怎么算全丢给 parameter_cal.m。好处是换评价标准不用碰仿真代码，坏处是文件多了容易漏——这也是为什么入口被固定成 set_para.m。

作者在说明里强调了一句「请从 set_para.m 文件开始运行，其他 M 文件和 slx 放在一个文件夹」，这句话不是客套。MATLAB 的搜索路径是按当前工作目录走的，你把 set_para.m 单独拖到桌面双击，PID_access.m 就找不到了，报的错还特别含糊，只说函数未定义。

2.2 技术选型：为什么这套东西押在 MATLAB 上

同样一件事用 Python 也能做，scipy 有 c2d，numpy 有矩阵运算，为什么还是 MATLAB 更顺手？原因在并行粒度的位置。

PSO 每一代要对 m 个粒子各做一次完整仿真，这 m 次仿真之间没有任何数据依赖，天然可以并行。在 MATLAB 里这件事几乎不需要额外代码——PID_sim.m 里把 kp、ki、kd 都 reshape 成 M×1 的列向量，仿真状态 y、u、Ee 全部是 M×N 的矩阵，一行y(:,k)=...就同时推进了 M 条独立轨迹。这就是原文说的「MATLAB 的矩阵运算已经集成了并行算法」的实际含义，不是调了 parfor，而是把循环维度摊到了矩阵的列上。

第二个原因是 Simulink 联动。复杂系统（比如双闭环直流调速）用传递函数手写离散化太费劲，直接在 Simulink 里拖模型更直观，而 MATLAB 和 Simulink 共享工作区变量，M 文件里改一个 global 变量，模型里的增益模块立刻跟着变。跨语言做这件事得走进程间通信或者文件交换，工程量完全不是一个量级。

第三个原因是 c2d。sys=tf(0.998,[0.021,1])拿到连续传递函数，dsys=c2d(sys,Tt,'z')一行零阶保持器离散化，[num,den]=tfdata(dsys,'v')再取出系数向量。这三步用 Python 的 control 库也能写，但 MATLAB 的这套 API 更贴控制工程的教材表述，改参数的时候不容易搞混。

2.3 一次完整迭代的调用链

真正跑起来之后，数据是这样流动的：

% 简化后的调用关系，不是可以直接跑的代码 X = (xmax-xmin).*rand(m,3)+xmin; % 第一代：m 组 PID 参数，随机撒在范围内 last_f = PID_access(X); % 打分：m 个适应度值 for i = 1:n v = w*v + c1*rand*(P-X) + c2*rand*(ones(m,1)*Pg-X); % 速度更新 X = X + v*dt; % 位置更新 new_f = PID_access(X); % 新一代打分 P(update_j,:) = X(update_j,:); % 个体历史最优 Pg = (new_fmin<fmin)*new_Pg + (new_fmin>=fmin)*Pg; % 种群最优 end

PID_access(X)是整个链路的咽喉，它内部做三件事：先把 X 的三列拆成 kp、ki、kd，再调PID_sim(kp,ki,kd,false)拿到 M 条响应曲线，最后parameter_cal从曲线里算出 ts 和 sigma，代进 log 公式。也就是说，一次PID_access调用等于 m 次闭环仿真，二十代迭代就是 4000 次——这个量级在普通笔记本上跑几分钟，属于可以接受的范围。

链路里有个容易被忽略的点：PID_access的返回值越小越好（它是评价函数不是适应度），而 PSO 内部一律用new_f < last_f判优。这个符号约定在 PSO.m 和 PID_access.m 之间必须一致，一旦有人手贱把评价函数改成「越大越好」却没同步改 PSO.m，算法会立刻朝反方向狂奔，参数顶到边界上不动。我第一次改这块就翻过车，后面会细说。

2.4 跑起来的三步

环境上原作者标注的是 MATLAB R2017a，版本差异主要影响 Simulink 模型的打开和部分语法兼容，M 文件部分一般问题不大。跑之前确认三件事：文件夹里有全部 M 文件和 slx；当前工作目录是这个文件夹；命令行窗口能看到迭代打印。

% 第一步：切到资源包目录（换成你自己的路径） cd('D:\work\pso_pid'); % 第二步：从参数入口开始，它会设置种群规模、迭代次数等 set_para; % 第三步：如果只想快速验证算法是否正常，用精简版主函数 dimension = 3; % PID 三个参数，维度固定是 3 n = 20; % 迭代代数，先小一点看收敛趋势 m = 200; % 种群规模，越大越稳但越慢 [Pg, fmin] = PSO(dimension, n, m)

set_para.m和main.m做的是同一件事，区别是前者参数更全、会把绘图开关一起打开，后者只是一段最小可运行代码。新手建议先跑main.m那三行，看到命令行吐出 Pg 和 fmin 再回头跑完整版，能省掉一半排查时间。

运行过程中每迭代一次都会打印当前最优参数对应的 ts 和 sigma，这是判断程序是不是还活着的唯一信号。作者特意提醒「如果运行时间过长怀疑程序死机可以观察迭代次数是否有变化」，因为这套代码在仿真阶段没有任何进度条，界面上就是一串数字慢慢往上爬，第一次跑的人很容易以为卡死了直接 Ctrl+C。三张 figure 是在迭代全部结束后一次性画出来的：响应曲线随迭代次数的变化、调节时间曲线、超调量曲线。

3. 核心算法拆解：PSO.m 里那段向量化的速度-位置迭代

PSO 本身的数学不复杂，难的是把它写成矩阵形式还不出错。这一章把 PSO.m 拆成四块：算法步骤和代码怎么一一对应、速度更新那一行里的三个因子分别在干什么、惯性权重为什么必须随迭代衰减、以及判优时用乘法代替 if 的技巧。

3.1 算法四步和代码的对应关系

先把教科书上的四个步骤摆出来，再逐条对代码：

第一步确定参数维度和范围。这一步在代码里体现为xmax、xmin、vmax、vmin四个向量，以及X=(xmax-xmin).*rand(m,dimension)+xmin这一行初始化。

w = 1; c1 = 2; c2 = 2; % 惯性权重和两个学习因子，经典取值 dt = 0.3; % 位移仿真间隔，小于 1 相当于把步子迈小 vmax = 1; % 速度限幅，防止粒子一步跨太远 xmax = [15, 50, 2]; % Kp、Ki、Kd 的上界 xmin = [0.2, 0, 0]; % 下界，Kp 不能取 0，否则比例环节失效 v = zeros(m, dimension); % 初始速度全零 X = (xmax-xmin).*rand(m,dimension)+xmin; % 在盒约束内均匀撒点

这里的dt=0.3是个关键参数。理论上位置更新是X = X + v，但直接加会让粒子在参数空间里跳得太猛，尤其在参数范围是[0.2,15]这种量级差异很大的情况下，一次跳跃就可能从 Kp=1 撞到 Kp=14。乘上 0.3 相当于给搜索加了个阻尼，原文说的「希望仿真得慢一点，搜索平滑一点，可以适当减小 dt」就是这个意思。

第二步按速度公式更新速度并计算新位置。第三步计算每个粒子的适应度，更新个体最优 P 和种群最优 Pg。第四步记录每一代的最优解用于画图。这四步在代码里被压进了一个 for 循环：

for i = 1:n % step2：速度更新 + 限幅 + 位置更新 + 限幅 v = w*v + c1*rand*(P-X) + c2*rand*(ones(m,1)*Pg-X); v = (v>vmax).*vmax + (v>=-vmax & v<=vmax).*v + (v<-vmax).*(-vmax); X = X + v*dt; X = (X>xmax).*xmax + (X>=xmin & X<=xmax).*X + (X<xmin).*(xmin); % step3：新一代打分，更新个体最优和全局最优 new_f = PID_access(X); update_j = find(new_f < last_f); P(update_j,:) = X(update_j,:); [new_fmin, min_i] = min(new_f); new_Pg = X(min_i,:); Pg = (new_fmin<fmin)*new_Pg + (new_fmin>=fmin)*Pg; last_f = new_f; fmin = min(new_fmin, fmin); end

3.2 速度更新一行代码的三个因子

v = w*v + c1*rand*(P-X) + c2*rand*(ones(m,1)*Pg-X)这一行是整个算法的核心，拆开看是三股力量的叠加。

第一项w*v是惯性项，保留上一次的飞行方向。w 取 1 意味着完全不衰减，粒子会保持原来的速度继续冲。这一项的存在让 PSO 有「记忆」，比纯随机的搜索更高效，但也是震荡的根源——如果粒子已经冲到最优解附近，惯性还在推着它往外跑，就会出现原文动图里那种「部分粒子反向运动」的现象。

第二项c1*rand*(P-X)是个体认知项，把粒子往它自己历史上到过的最好位置拉。这里的P是 m×3 矩阵，每行是一个粒子自己的历史最优，P-X就是指向自身最优的向量。rand是 m×3 的随机矩阵，为每个维度单独生成随机数，这是 PSO 随机性的来源之一。

第三项c2*rand*(ones(m,1)*Pg-X)是社会认知项，把粒子往整个种群的历史最优拉。注意ones(m,1)*Pg这一步：Pg 是 1×3 的行向量，乘上 m×1 的全一列向量之后变成 m×3，每行都是同一个 Pg，这样才能和 X 逐元素相减。如果不做这个广播，MATLAB 会报维度不匹配。

三股力量各自的系数是 w、c1、c2，经典取值就是 1、2、2。调整这三个数会明显改变搜索行为：c1 大偏探索、每个粒子各飞各的；c2 大偏开发、整个种群很快往一个点收拢。原文在「修改参数价值权重」那一节说的就是这个——如果希望算法在前期多铺开、后期多收敛，常规做法是让 c1 随迭代递减、c2 递增，或者只调 w。

紧接着的速度限幅那一行值得单独说：

v = (v>vmax).*vmax + (v>=-vmax & v<=vmax).*v + (v<-vmax).*(-vmax);

这里没有用min(max(v,-vmax),vmax)，而是用三个逻辑掩码做加权。三个条件是互斥且完备的（大于上界、区间内、小于下界），每行每个元素只会命中其中一个分支，所以结果是正确的截断。写法丑但快，避免了嵌套函数调用的开销——在 m=200、dimension=3、迭代 20 代的规模下这点开销感知不到，但如果把种群扩到 2000，这个细节就开始有意义了。

位置限幅同理，X=(X>xmax).*xmax+(X>=xmin & X<=xmax).*X+(X<xmin).*(xmin)，把跑出盒约束的粒子按回边界上。注意第三个分支写的是.*(xmin)，不是.*xmin，虽然结果一样，但说明作者写这段的时候思路是「命中下界的元素赋值为 xmin」。

3.3 惯性权重衰减：w 从定值改成 i/n 衰减之后

原文里有一个对比实验很说明问题：同一测试函数、同样的种群规模和迭代次数，唯一区别是 w 是否随迭代衰减，结果 fmin 差了一大截。

% 衰减版本，放在 for 循环末尾 w = w - i/n*0.7*w;

这行的含义是：第 i 代时，w 相比初始值衰减了i/n*0.7的比例。跑到最后一代时 w 只剩初始值的 30%，惯性项被压到很小，粒子基本只受个体最优和全局最优牵引，自然就收敛到一点上。

为什么必须衰减？回到速度公式看：如果 w 恒为 1，粒子冲过最优点之后，惯性项还在把它往外推，而认知项和社會项又把它往回拉，两股力量打架，粒子就在最优点附近来回震荡，收敛不了。原文的测试函数是f=(1-x1)^2+(x1^2-2x2)^2，理论极小值在 (1, 0.5)，不衰减的版本跑完之后粒子分布还很散，衰减版本几乎全部聚到中心点。

不过衰减曲线选 0.7 这个系数是有讲究的。衰减太快（比如 0.9）前期探索不足，容易早熟收敛到局部最优；衰减太慢（比如 0.3）后期还在乱飞，收敛精度上不去。0.7 是个经验值，在 PID 整定这种三维、范围不算极端的场景下比较稳。如果换成高维参数整定，常见做法是改用线性递减w = w_max - (w_max-w_min)*i/n，把上下界显式写出来，调起来更直观。

顺便提醒一句：PSO.m 里那行衰减被注释掉了（原文代码里前面有个%），意味着默认跑的是不衰减版本。第一次跑发现粒子散成一片、fmin 卡住不降，先把这行的注释去掉再试。

3.4 适应度比较与 Pg 更新：为什么用乘法代替 if

判优这段是整份代码里最容易看错的地方：

update_j = find(new_f < last_f); % 找出这一代比上一代更好的粒子索引 P(update_j,:) = X(update_j,:); % 只更新这些粒子的个体最优 [new_fmin, min_i] = min(new_f); % 这一代的最优值和它的粒子编号 new_Pg = X(min_i,:); % 这一代的候选全局最优 Pg = (new_fmin<fmin)*new_Pg + (new_fmin>=fmin)*Pg; % 二选一 fmin = min(new_fmin, fmin); % 保留历史最小值

find(new_f < last_f)返回的是逻辑真值的索引，只有这一代确实变好的粒子才更新 P。这比「无条件覆盖 P」正确得多——PSO 的个体最优定义是「历史最优」，不是「当前最优」，搞混了算法就退化成随机搜索。

Pg那一行用了一个布尔乘法技巧。(new_fmin<fmin)是逻辑标量，真为 1 假为 0，两路相乘相加正好构成一个二选一开关：新值更好就取 new_Pg，否则保持原 Pg。逻辑上等价于if new_fmin < fmin, Pg = new_Pg; end，但写成表达式更紧凑，也避免了在向量化代码里混进控制流。

这里有个隐含前提必须重申：new_f < last_f意味着适应度越小越好。如果哪天你把PID_access改成返回「越大越好」的评分，这两处判优全部要反过来写，否则算法会稳定地挑出最差的参数。我在改评价函数的时候就栽在这上面，后面避坑那章会展开。

4. PID自整定实现：PID_sim.m 到 PID_access.m 的指标流水线

PSO 只是个搜索框架，真正决定整定效果的是被控对象建模准不准、指标算得对不对、评价函数调得合不合理。这一章沿着「传递函数 → 离散化 → 批量仿真 → 指标提取 → 标量评分」这条流水线走一遍。

4.1 被控对象离散化：tf、c2d、tfdata 三段式

PID_sim.m 开头几行建立被控对象：

Tt = 5e-4; % 仿真采样周期，0.5ms T0 = 1e-2; % 控制采样周期，10ms Tf = 1e-3; % 微分环节的滤波器时间常数 alp = Tf/(Tt+Tf); % 一阶低通滤波系数 total_t = 1; % 仿真总时长 1s N = floor(total_t/Tt); % 样本总数，2000 个点 sys = tf(0.998,[0.021,1]); % 一阶惯性环节，模拟直流电机转速 dsys = c2d(sys,Tt,'z'); % 零阶保持器离散化 [num,den] = tfdata(dsys,'v');% 取出分子分母系数向量

tf(0.998,[0.021,1])对应传递函数0.998/(0.021s+1)，时间常数 21ms，增益接近 1。这是一阶惯性环节的标准形式，原文说它模拟的是直流电机转速模型——电机的电气时间常数远小于机械时间常数时，从电压到转速可以近似成一阶惯性。

c2d(sys,Tt,'z')的第二个参数是采样周期，第三个参数是离散化方法。'z'是 MATLAB 的默认零阶保持器（zero-order hold），物理含义是「每个采样周期内输入保持不变」，正好对应数字控制器的实际行为。如果换成'tustin'（双线性变换），频域特性会略有不同，但对这种一阶对象影响很小。

tfdata(dsys,'v')里的'v'表示返回向量形式而不是 cell 数组。一阶离散系统的 den 是[1, -a]两个元素，num 是[0, b]或[b0, b1]两个元素，取num(1)、num(2)、den(2)就能直接写差分方程，不需要再解析结构体。

这里有个可以改的地方：如果你的被控对象是二阶或更高阶，den会变成三个以上元素，y(:,k)=-1*den(2).*y_1+...那一行就不够用了，需要改成按阶数展开的历史值数组。这也是为什么这套代码只适合拿来验证算法思路，真要上工程得自己重写仿真引擎。

4.2 PID_sim.m 逐段拆：矩阵列并行的仿真循环

核心循环如下，M 是这一批要仿真的参数组数，N 是时间点数：

M = numel(kp); % 这一批要仿真多少组参数 kp = reshape(kp,M,1); % 统一成列向量，方便逐元素运算 ki = reshape(ki,M,1); kd = reshape(kd,M,1); 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) = -1*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

逐段说。reshape(kp,M,1)是把输入统一成列向量，这样后面kp.*e0才能和e0（M×1）逐元素相乘，得到 M 组并行的控制量。

状态变量的初始化全是zeros(M,1)，每个状态都是一列 M 个值，对应 M 条独立轨迹。这是整个并行设计的根基——没有 for 循环遍历粒子，靠的是矩阵的每一行代表一个粒子。

对象差分方程y(:,k) = -den(2).*y_1 + num(2)*u_1 + num(1).*u(:,k)，是从Y(z)/U(z) = (num(1)+num(2)z^-1)/(1+den(2)z^-1)反推出来的时域形式。注意u(:,k)此刻还是预分配的 0（要等下一行才赋值），所以实际用的是「上一拍的输入」，等价于给对象加了一拍延迟。这不是 bug，离散控制本来就有这个延迟，但如果哪天你把循环顺序调整了，这个隐含时序就被破坏了。

微分项那行ud = alp.*ud_1 + (1-alp)/T0.*kd.*(e0-e_1)值得单独讲。纯微分对噪声极其敏感，实际工程里几乎没人直接用kd*(e0-e_1)/T0，都要加一阶低通滤波。这里的alp = Tf/(Tt+Tf) = 1e-3/(5e-4+1e-3) ≈ 0.667，滤波器越接近 1 平滑越强、相位滞后越大。这是标准的「不完全微分」实现，很多教材上写 PID 只写理想形式，一到工程落地就得补这一块。

积分项T0*ki.*Ee用的是控制周期 T0 而不是仿真周期 Tt，因为积分是对控制节拍累积的。这里 T0=10ms、Tt=0.5ms，比例是 20:1，意味着控制量每 20 个仿真步才更新一次——但代码里其实每一步都在算u(:,k)，没有真正做零阶保持。这是一个简化，对一阶对象影响不大，换成快速响应对象就需要补上采样保持逻辑。

4.3 parameter_cal.m：稳态值、超调、调节时间怎么从曲线里抠出来

仿真跑完拿到 M×N 的 y 矩阵，接下来要把它压缩成每组参数对应的几个标量指标：

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

f_infty = y(:,N)用最后一个采样点当稳态值，这个做法只在「仿真时间足够长、系统已经稳定」的前提下成立。如果 total_t 设得太短，最后一拍还在上升，稳态值就被低估，超调量算出来会偏大甚至为负。这是调参时最常见的翻车点之一。

调节时间的算法比较巧妙：把误差序列取绝对值后翻转，从末尾往前扫，max(ferr > delt_err*f_infty, [], 2)找到的是「最后一个超出阈值的点」在翻转序列里的位置，再用 N 减回去就得到正序的「第一次进入误差带的位置」。delt_err默认传的是 2e-2，即 2% 误差带，这是控制工程里的常规取法。

tp(abs(fp-f_infty)<1e-5) = NaN这行处理的是过阻尼情况——曲线单调上升，峰值就是稳态值本身，这时候峰值时间没有物理意义，标成 NaN 让后续统计忽略掉。但注意 PID_access.m 里直接用sigma计算，如果 sigma 因为 fp 和 f_infty 相等而变成 0，log 里的log(0/1e-2+1)=0，没问题；如果 sigma 变成负数（稳态值取大了），log的参数就可能是负的，直接报错。这也是为什么稳态值必须取准。

4.4 评价函数与权重：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); % debug 关掉，不画图不打印 y = log(ts/5e-2 + 1) + log(sigma/1e-2 + 1); % 评价度，越小越好 end

这是一个「越小越好」的评分，PSO 内部用new_f < last_f判优，两边口径一致。

log 形式的妙处在于对量级不敏感。ts 可能是 0.05 也可能是 0.5，sigma 可能是 0 也可能是 0.5，直接相加会让量级大的那一项完全主导。取 log 之后，两项都被压到相近的区间，谁都没法一家独大。

/5e-2和/1e-2这两个除数是归一化尺度，物理含义是「期望值」。ts 除以 0.05 意味着调节时间 50ms 算一档，sigma 除以 0.01 意味着超调 1% 算一档。原文在讨论里用的版本是log(sigma/1e-4+1)，代码里写的是1e-2，这两个值的差别很大：用 1e-4 时算法会疯狂压超调（因为超调项的梯度被放大了），用 1e-2 时更平衡。想清楚自己要什么，就调这两个分母。

log 还有个性质：log(x+1)在 x 很小时近似等于 x，在 x 很大时增长变慢。原文推导过导数，在期望值附近∂f/∂ts = 0.5/5e-2 = 10，也就是说 ts 每增加一点，评分的惩罚是线性的、不会突变。这种平滑性对 PSO 很重要——如果评分函数有阶跃，粒子在边界附近会来回跳，收敛不了。

想加第三个指标（比如上升时间、稳态误差），直接在这行里再挂一个 log 就行，但要注意三个分母的相对大小，它们决定了算法更在意哪一项。这就是原文说的「修改参数价值权重」——不用改算法，改评价函数就够了。

5. Simulink 双闭环落地与避坑排查：模型改造、版本、收敛三件事

M 文件版跑通只能算热身，真实场景多半是一个 Simulink 模型摊在那儿，参数藏在各个 Gain 模块里。这一章先说双闭环模型怎么改造成可被 M 文件驱动，再把实际跑下来最常撞的五个坑摆出来。

5.1 把双闭环模型改成可写参数的形式

双闭环直流调速系统由转速外环和电流内环构成：外环是 PI 调节器，管转速、抗负载扰动；内环也是 PI，管电流快速跟随。启动过程是这样的——电流先冲到过载值，电机以最大加速度加速，转速到给定值后电流迅速回落，维持在负载电流上，转速稳住。

原文给出的电机参数是：额定 800V/2164A/1500kW，额定转速 44r/min，过载系数 1.5，励磁 110V/191A，电枢电阻 0.1Ω，主回路电阻 0.1014Ω，电感 3.25mH，变流器放大倍数 80，通态压降 3V，飞轮惯量 8.6×10^5 N·m²。按工程设计方法算出来的 PI 参数是[Knp, KnI, Kip, KiI] = [7.3660, 1.5658, 0.8257, 30.2773]，对应的指标是转速超调 4.07×10^-4 %、转速调节时间 0.7875s、电流超调 20.38%、电流调节时间 0.0172s。

要让 PSO 能驱动这个模型，改造步骤是这样的：

% 1. 在 Simulink 里把四个 PI 增益模块的值改成变量名 % ASR 的比例增益填 Knp，积分增益填 KnI % ACR 的比例增益填 Kip，积分增益填 KiI % 模块不能填数字，否则 PSO 改不动 % 2. M 文件里声明全局变量，让 Simulink 能读到 global Knp KnI Kip KiI % 3. 每一代把候选参数写进全局变量，然后启动仿真 Knp = Pg(1); KnI = Pg(2); Kip = Pg(3); KiI = Pg(4); sim('PSO_double_circle'); % 模型名要和 slx 文件名一致 % 4. 从工作区取仿真结果，交给指标计算函数 [n_sigma, ts_n, I_sigma, ts_ti, n_err] = cal_n(...);

global这一步是必须的。Simulink 模型在仿真时会去 base 工作区找变量值，M 文件函数内部的局部变量它看不见，只有声明成 global 或者用assignin('base', 'Knp', ...)显式写进 base 工作区才行。原作者在说明里写的「将变量声明为全局（global），此时变量会在工作区里面，simulink 就可以调用了」，说的就是这个机制。

另一个容易忽略的点是模型另存。原始模型double_circle.slx是用来验证工程设计参数的，改造成 PSO 驱动之后要另存为PSO_double_circle.slx，否则你调一次参数就把原始模型覆盖了，回头想对比手工调参的结果都找不到基准。

cal_n.m的作用和 M 文件版的parameter_cal.m一样，只是指标从两项变成四项（转速超调、转速调节时间、电流超调、电流调节时间）。因为双闭环的电流环响应快得多，采样精度要单独设置，不能和转速环共用一个误差带。

5.2 五条踩坑记录：现象、原因、解决

下面这些是我和身边人实际撞过的，按出现频率排。

坑一：换台机器就报「未定义函数或变量 PID_access」

现象：把文件夹拷到另一台电脑，双击set_para.m运行，命令行报找不到PID_access或PID_sim。

原因：MATLAB 的函数搜索路径只包含当前工作目录和已添加的路径。双击 M 文件时，当前目录可能被设成了桌面或者其他位置，而PID_access.m躺在资源包文件夹里。另外如果只下载了 M 文件、漏了 slx，报错会更隐蔽——代码跑到sim()那一步才炸。

解决：先cd到资源包目录，再addpath(pwd)保险一次，最后确认dir *.m能列出全部文件。养成习惯，每次开新会话先pwd看一眼。

坑二：MATLAB 版本不同，Simulink 模型直接打不开

现象：在 R2021b 里双击double_circle.slx，弹出「该模型由较新版本创建」或者打开之后模块位置错乱。

原因：Simulink 模型文件是二进制格式，高版本存的文件低版本读不了；反过来低版本存的高版本能读，但可能丢格式。原作者标的是 R2017a，用 2017a 之后的版本打开老模型一般没问题，用 2017a 打开新存的模型就悬了。

解决：在 Simulink 里用「文件 → 导出 → 早期版本」另存一份，或者干脆按 5.1 的步骤在新版本里重建模型——模块不多，重建一遍比折腾兼容性快。M 文件部分如果有语法报错，多半是~占位符或者某些函数名在旧版本不存在，逐行对着改就行。

坑三：跑了几分钟以为死机，Ctrl+C 之后前功尽弃

现象：命令行窗口停在一行不动，CPU 占用很高但没有任何输出。

原因：一次PID_access要跑 m 次完整仿真，每次 2000 个时间步，m=200 的话就是 40 万次循环，二十代下来接近 800 万次。中间没有任何进度提示，视觉上就是静止的。

解决：看命令行当前的迭代编号有没有在变。如果确实要中断，别用 Ctrl+C 硬杀，先在PSO.m的循环里加一句fprintf('iter %d, fmin %.4f\n', i, fmin)，让每代都有输出。另外第一次调试把 m 降到 20、n 降到 5，整个流程几十秒跑完，确认逻辑对了再上规模。

坑四：适应度不降反升，参数全部顶在上界

现象：迭代打印出来的 fmin 越来越大，最后 Pg 稳定在[15, 50, 2]这个上界上不动。

原因：判优方向搞反了。要么是PID_access被改成了「越大越好」，要么是new_f < last_f被改成了>。还有一种可能是 xmax 设得太小，真正的优解在范围之外，粒子全体撞墙。

解决：先打印一代的new_f看数值分布，正常的评价函数值应该在个位数到几十之间；如果全是几百上千，说明仿真本身发散，PID 参数给了个会让系统震荡的值。再把xmax放宽一档试，比如 Kp 上限从 15 提到 30。最后检查PID_access里两个 log 的分母是不是正数，分母为负会导致 log 输出负数，整个评分体系翻转。

坑五：超调量算出 NaN 或负值，调节时间等于仿真总长

现象：sigma出现负值，ts恒等于 total_t。

原因：f_infty = y(:,N)取的是最后一个采样点，如果仿真时长不够、系统还在上升段，这个「稳态值」就是个中间值，后面算出来的误差、峰值全是错的。另一种情况是 PID 参数给了纯积分，系统一直震荡不收敛。

解决：把total_t从 1s 延长到 2s 或 3s，再看最后一小段曲线是不是平的。parameter_cal里可以加一句plot(y(end,:))抽查，肉眼判断稳态。实在拿不准稳态值，改用最小二乘拟合后 10% 的数据段，比直接取末点稳。

坑六（Simulink 专属）：改了 global 变量，模型里没反应

现象：Knp在命令行里明明改了，sim()跑出来的响应还是老样子。

原因：模型里的 Gain 模块填的还是具体数字，不是变量名；或者变量的作用域不对，M 文件函数里global了但 base 工作区里没同步。

解决：点开 Gain 模块确认 Value 栏写的是Knp而不是7.366；在sim()之前加一句assignin('base','Knp',Knp)强制同步；跑完用whos global看一眼变量在不在。

6. 二次开发：换对象、改评价、接自己的系统

把 PSO.m 当模板用，改造点其实只有三处。

换被控对象，动PID_sim.m里的sys=tf(...)。一阶惯性改成二阶振荡环节时，den会从两个元素变成三个，对象差分方程那行要扩成y(:,k) = -den(2)*y_1 - den(3)*y_2 + num(2)*u_1 + num(3)*u_2，同时补y_2、u_2两个历史状态。参数范围xmax/xmin也要跟着重新给，二阶对象的 Kp 通常比一阶小一个量级，沿用[15,50,2]会让粒子一开始就发散，sigma全变 NaN，评分失去意义。

改优化目标，动PID_access.m里那行 log。加稳态误差就再挂一项log(ess/1e-3+1)，想压上升时间就把ts换成tr。改完必须验证方向——打印第一代和第五代的fmin，正常应该单调下降；如果上升，先查是不是哪个分母取了负数。这一步我建议写成一个小脚本，跑三代就停，别等到二十代跑完才发现方向反了。

接自己的 Simulink 模型，关键是参数注入方式的选择。global最简单但容易污染工作区，项目大了会和其他脚本打架；assignin('base',...)更可控，每次仿真前同步一次；sim的'SrcWorkspace'参数可以把变量源指到函数工作区，最干净但要求 MATLAB 版本较新。我一般用assignin，改起来直观。

还有一个值得加的功能：把每一代的 Pg 和对应指标存进数组，迭代结束后画一张收敛曲线。原代码只存了sigma_data和ts_data，没存 Pg 本身，回头想复现某一代的最好参数就找不到了。加两行Pg_history(i,:) = Pg;的事，省掉后面无数麻烦。

验证方法上，我的习惯是拿工程设计的经典参数当基准。比如双闭环那组[7.3660, 1.5658, 0.8257, 30.2773]，先算一遍它的评分，再让 PSO 跑出来的结果和它比——如果 PSO 的评分比基准还差，说明参数范围或者评价函数出了问题，别急着怀疑算法。从那以后我每次改评价函数，都强制先跑一遍基准参数做对照，确认新评分下基准依然排在前列，才敢让 PSO 全量开跑。希望帮到你。

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