简介围绕Simulink环境下粒子群算法PSO整定PID参数的一份PDF学习资料面向已具备MATLAB基础、希望摆脱人工试凑调参的控制工程学习者与自动化方向学生。内容从PSO的群体迭代思想、速度与位置更新步骤讲起结合M文件自整定和Simulink仿真模型说明如何把Kp、Ki、Kd作为待优化变量以超调量、调节时间等性能指标构造评价函数并分析惯性权重及c1、c2取值对探索与开发平衡的影响末尾还提示了代码运行与迭代观察的注意事项。压缩包共1个PDF文件约529KB篇幅集中适合作为算法原理与实现流程的入门导读。目前已有4895人学习读者可据此理解PSO优化PID的完整链路并对照代码开展自己的参数整定实验。1. 调参调到怀疑人生的时候可以换一群鸟去试我记得最清楚的一次调 PID是给一台带负载扰动的电机做转速环。Kp 大一点超调小一点爬得像蜗牛Ki 一加就开始抖最后三个参数在六个方向上互相拉扯改了十几个版本谁也没比第一版好多少。粒子群算法PSO在 Simulink 里整定 PID 参数这件事本质就是把这轮拉锯交给一群会互相通报位置的鸟每飞一步就是一次仿真谁的位置对应的响应曲线更好全场就朝谁靠。这份资源给的是 MATLAB 2017a 下两套能直接跑的完整闭环——一套纯 M 文件的 PID 引擎一套双闭环直流调速的 slx 模型外面套同一个 PSO 找参。适合已经被手工整定折磨过、想把参数寻优自动化的人也适合想学智能优化怎么真正接上仿真模型的读者。入口只有一个文件set_para.m从它开始跑别走岔。2. 工程骨架拆解一个 SRC 从 set_para.m 到图窗的完整链路这个资源拿过来第一眼看不懂不是代码难是文件之间的调用关系没画出来。MATLAB 工程不像 Java 有明确的 main 方法它靠文件名和工作区变量串起来一旦顺序错了报错会指向一堆莫名其妙的地方。我通常的做法是先画一张调用图再决定从哪个文件按下 F5。2.1 七个文件的分工与运行顺序把资源包解开核心是下面这几样东西它们的职责边界很清楚文件职责输入输出set_para.m参数总入口声明全局变量与搜索边界无工作区变量PSO.m粒子群主循环负责位置速度更新维度、代数、种群数最优参数 Pg、最优适应度PID_access.m适应度函数把一组参数映射成一个标量m×3 或 m×4 参数矩阵m×1 适应度列向量PID_sim.m批量 PID 仿真引擎一次跑 M 组参数kp、ki、kd仿真曲线与指标parameter_cal.m从响应曲线里量出 ts、sigma 等指标输出序列 y稳态值、峰值时间、调节时间、超调main.m简化版主函数M 文件路线下的直接入口无三张 figuredouble_circle.slx / PSO_double_circle.slx双闭环调速的 Simulink 模型工作区变量电流、转速曲线调用链是这样的set_para.m 把种群规模 m、迭代代数 n、参数上下限 xmax/xmin、速度限幅 vmax 塞进工作区PSO.m 拿到这些数后生成初始种群然后每一代调用一次 PID_access.mPID_access.m 内部把整代参数矩阵丢给 PID_sim.m 批量仿真PID_sim.m 跑完曲线后调 parameter_cal.m 量指标量出来的指标回到 PID_access.m 算成适应度适应度再回到 PSO.m 决定谁当 Pg。整个工程只有一个真正的外部输入就是被控对象的传递函数或 slx 模型。其余全是可替换的零件这也是后面二次开发好下手的原因。提示不要单独点开 PSO.m 直接运行。没有 set_para.m 铺好的全局变量PSO.m 里引用的 xmax、vmax 全是未定义报错信息会指向速度更新那一行很容易误判成矩阵维度问题。2.2 向量化的粒子群一次迭代算完 200 组 PIDPSO 里最容易写错也最值得看的是主循环那几个矩阵运算。它没有对粒子做 for 循环而是把 m 个粒子堆成一个 m×dimension 的矩阵一次运算全部推进。这样做的好处是能吃到 MATLAB 的 BLAS 加速200 个粒子的规模下比逐个循环快一个数量级。% 位置与速度按行组织每一行是一个粒子的候选 PID 参数 X (xmax - xmin) .* rand(m, dimension) xmin; % 初始位置在框内均匀撒点 v zeros(m, dimension); % 初始速度归零 P X; % 每个粒子的历史最优位置 last_f PID_access(X); % 整代一次算完适应度 for i 1:n % n 为迭代代数 % 速度更新惯性项 个体认知项 群体认知项 v w*v c1*rand(m,dimension).*(P - X) ... c2*rand(m,dimension).*(ones(m,1)*Pg - X); % 速度限幅防止粒子一步飞出有效搜索区 v (v vmax).*vmax (v -vmax v vmax).*v (v -vmax).*(-vmax); % 位置更新dt 相当于一次迭代的位移步长 X X v*dt; % 位置限幅撞到边界就贴边停住 X (X xmax).*xmax (X xmin X xmax).*X (X xmin).*xmin; 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); Pg (new_fmin fmin)*X(min_i,:) (new_fmin fmin)*Pg; last_f new_f; fmin min(new_fmin, fmin); end几个参数值得逐个抠一下。c1和c2原稿都取 2w取 1这是最经典的取值组合。c1管我历史上哪儿最好偏向个体探索c2管全场哪儿最好偏向群体收敛。两个都大粒子就来回拉扯把c2调大收敛快但容易早熟。dt原稿给 0.3这个数不是物理时间是一次迭代允许粒子移动多大比例的距离。设成 1 的时候粒子经常一步跨过头再被速度惯性拖回来表现就是曲线在最优解附近反复横跳很久不收敛。ones(m,1)*Pg是把 1×3 的全局最优复制成 m×3好让整个种群朝同一个点靠。R2016b 之后 MATLAB 支持隐式扩展其实直接写Pg - X也行但显式复制在旧版本上更保险。至于(v vmax).*vmax ...这种限幅写法是典型的逻辑索引做分段函数比min(max(v,-vmax),vmax)可读性差一点但能一次处理 NaN 以外的所有情况也避开了min/max在两个同尺寸矩阵上逐元素比较的语义陷阱。还有一个容易被忽略的细节原稿里写的是c1*rand*(P-X)无参的rand只返回一个标量等于这一代所有粒子共享同一个随机系数。这种写法收敛更整齐但也丢掉了粒子间的多样性。我一般改成rand(m,dimension)每个维度每个粒子独立随机种群散得更开代价是收敛代数要相应加一些。2.3 适应度函数把曲线好不好看压成一个数PSO 不懂什么叫超调它只认一个数字越小越好。所以真正决定整定质量的不是算法是 PID_access.m 里那个公式。function y PID_access(para) % para 的每一行是一组 [kp, ki, kd] 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为什么要套log而不是直接线性相加因为调节时间和超调量的量级差得太远。ts 是 0.1 秒量级sigma 是 0.01 量级直接相加的话 ts 一项就把 sigma 完全盖住了PSO 会疯狂压 ts 而不管超调。取对数把 0 到 1 秒的 ts 压到 0 到 3 的区间把 0 到 0.5 的 sigma 压到 0 到 4 的区间两者才有可比性。更细的一层是log(x/k 1)这个结构。它对 x 的导数在 x0 处等于 1/k在 x 很大处趋近于 0。也就是说当指标已经小到接近期望值时继续压它的收益急剧下降当指标远大于期望值时每一点改善都能换来明显适应度提升。这正好符合工程直觉调节时间从 5 秒压到 1 秒很值钱从 0.05 秒压到 0.04 秒基本没意义还可能把系统调到病态。分母里的5e-2和1e-2就是期望值等于告诉 PSO调节时间做到 50 毫秒、超调做到 1% 就算达标再往下不额外奖励。想换口味就改这两个数不用动算法一行。适应度方向也要注意PSO 默认是找最小值而这个y本身就是越小越好所以打分函数直接返回 y最优记录用min和。如果哪天你改成越大越好的评分记得整体取负别只改了比较符。2.4 惯性权重衰减为什么原稿要乘 0.7固定w1跑迭代会看到粒子在最优解附近长时间打转甚至反向飞走。原稿给的补救是让 w 随代数线性衰减% 放在 for i 1:n 循环末尾i 是当前代数 w w - i/n*0.7*w; % 每代缩水末代约剩 0.3 倍前几代 w 大粒子敢跨区域探索全局搜索能力强后几代 w 小惯性拖拽弱粒子被 pBest 和 gBest 牢牢吸住做局部精修。用原稿的测试函数f(1-x1)^2(x1^2-2x2)^2对比不衰减时 fmin 会卡在 0.1 量级加上衰减后能压到 1e-3 以下粒子分布图也从散成一圈变成聚成一个点。这条衰减曲线不是唯一解。R2017a 之后 MATLAB 自带particleswarm函数内部用的是自适应惯性权重加邻域拓扑比自己手写省事。但手写版的好处是每一步都能打印、能改、能插自己的约束做项目时反而更可控。3. M 文件版 PID 引擎在接 Simulink 之前先把闭环跑通很多人一上来就搭 slx结果模型跑不通、PSO 也不收敛两个问题搅在一起根本分不清是谁的错。我习惯的顺序是先用 M 文件写一个最小闭环确认 PSO 能在这上面找到合理参数再换成 Simulink 模型。这样出了问题至少知道是算法侧还是模型侧。3.1 一阶惯性对象的离散化与差分方程PID_sim.m 是整条链路里唯一真正跑物理的地方。它做的事情是把一个连续传递函数离散化然后用差分方程一步一步推控制量和输出M 组参数同时推进。function [f_infty, tp, ts, sigma] PID_sim(kp, ki, kd, debug) Tt 5e-4; % 仿真步长 T0 1e-2; % 控制周期PID 每隔这么久算一次 Tf 1e-3; % 微分项低通滤波的时间常数 alp Tf/(Tt Tf); % 滤波系数越接近 1 滤得越狠 total_t 1; % 仿真总时长 N floor(total_t/Tt); M numel(kp); % 这一批要同时仿多少组参数 % 被控对象一阶惯性环节模拟直流电机转速模型 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] -a1*y[k-1] b1*u[k-1] y(:,k) -den(2)*y_1 num(2)*u_1; e0 r - y(:,k); % 本拍误差 % 不完全微分的离散实现避免纯微分放大噪声 ud alp*ud_1 (1-alp)/T0*kd.*(e0 - e_1); % 位置式 PID 输出 u(:,k) kp.*e0 T0*ki.*Ee ud; 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这段代码里有三个时间常数关系不能乱。Tt是数值积分的步长越小越准但越慢T0是控制周期也就是离散 PID 的实际采样间隔Tf是微分滤波时间作用是把kd*(e0-e_1)这个差分会放大高频噪声的问题压下去alp越接近 1微分项越钝。有个非常容易踩的坑积分项写的是T0*ki.*Ee而Ee是每个仿真步都累加一次的。这意味着积分增益实际被乘了T0而累加用的是Tt。如果哪天你把Tt从 5e-4 改成 1e-4 而忘了动T0同样的 ki 值在系统里的等效强度会变PSO 前面搜出来的参数直接作废。所以这三个数一旦定下来就别动要动就整体重跑。M numel(kp)这个写法是为了兼容两种调用PSO 传一整代进来时 M200单独调试传一个标量进来时 M1。整个函数里所有的zeros(M,1)和列向量运算都是为了让这两种情况共用一套代码debug M1才打印指标也是这个原因。3.2 parameter_cal.m从一条曲线里量出调节时间和超调仿真跑出来只是 M×N 的数值矩阵PSO 看不懂曲线得有人把调得怎么样翻译成数字。这个翻译工作的质量直接决定整定结果靠不靠谱。function [f_infty, tp, ts, sigma] parameter_cal(y, T0, delt_err, debug) % y 的每一行是一条响应曲线T0 是采样时间delt_err 是误差带宽度 M size(y,1); N size(y,2); f_infty y(:,N); % 用末点近似稳态值 err y - f_infty*ones(1,N); % 各时刻与稳态的偏差 ferr fliplr(abs(err)); % 倒序ferr(:,1) 是最后一拍 % 从尾部往前找第一个超出误差带的位置即最后一次越界的时刻 [~, 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; % 超调量相对值 if debug M 1 fprintf(稳态值 %.4f超调 %.2f%%调节时间 %.4f s\n, ... f_infty, sigma*100, ts); end endfliplr加max这个组合是全篇最巧的一处。要判断最后一次进入误差带是什么时候正向扫描得写循环加状态机倒过来之后第一个超过阈值的点就是答案max对逻辑数组返回第一个 true 的索引一行搞定。这是 MATLAB 里处理最后一次满足条件的标准套路值得记住。误差带的宽度delt_err传的是2e-2也就是稳态值的 2%。工程上调速系统的调节时间一般按 ±2% 或 ±5% 算取 2% 更严格PSO 会被逼着把曲线压得更贴。如果发现怎么迭代 ts 都不降先怀疑这个数是不是定得太狠。f_infty y(:,N)这个近似在仿真总时长足够长时没问题但如果total_t1而系统时间常数偏大末点还没到稳态稳态值本身就偏后面的 ts 和 sigma 全是错的。判据很简单稳态值应该接近给定的 r1500差太多就说明仿真时间不够。峰值时间的 NaN 处理也必要。过阻尼系统根本不存在峰值max会返回末点如果不屏蔽sigma 会算成 0看着挺好实际是误判。原稿把这种情况统一置 NaN让后续统计能区分没有超调和超调恰好为 0。3.3 main.m 与三个可调旋钮M 文件路线的入口就三行但每一行都对应一个需要根据算力调整的量dimension 3; % 待整定参数个数kp、ki、kd n 20; % 迭代代数 m 200; % 种群规模 [Pg, fmin] PSO(dimension, n, m);dimension决定了搜索空间的维度也决定了xmax、xmin、vmax这几个向量的长度改的时候必须一起改否则隐式扩展会静默出错或者直接报维度不匹配。双闭环那个场景里维数是 4两个 PI 各两个参数就是因为这个。m和n的乘积就是仿真总次数。m200、n20 意味着 4000 次单组仿真每次 2000 个时间步在本机上大概跑几分钟到十几分钟取决于是否每代都画图。原稿特意提醒每次迭代都会打印出来怀疑死机就看迭代次数有没有变这不是废话是真的会遇到——把PID_sim的 debug 打开200 条曲线画在一个 figure 上渲染时间能占掉总耗时的八成。调参顺序我一般是先固定 n15 把 m 拉到 300确认搜索空间够宽、粒子能覆盖到合理区域再回头把 m 压到 100、n 提到 40 做精修。种群大解决漏掉好解代数多解决收敛不彻底两件事分开调比一起加更容易看出效果。跑完之后工作区里会留下三张图响应曲线随迭代的变化、调节时间随迭代的变化、超调量随迭代的变化。前一张用来看收敛趋势后两张用来看是不是某个指标被牺牲了。如果 ts 一路下降而 sigma 一路飙升说明适应度函数里 ts 的权重压得太重回去调log里的分母。4. Simulink 双闭环调速把 PI 参数从数字变成可注入的变量M 文件版验证完算法接下来才是真正有价值的场景模型复杂到写不出差分方程只能靠 Simulink 搭。这时候 PSO 的角色从驱动方程变成驱动模型接口从函数调用变成工作区变量这一步是整个工程最容易翻车的地方。4.1 双闭环模型的结构与四个待整定增益这个案例选的是直流调速系统的经典双闭环结构外环转速环ASR用 PI 调节器负责稳住转速、抵抗负载扰动内环电流环ACR也用 PI负责让电流快速跟上给定值启动时把电流顶在过载限幅上让电机以最大加速度冲起来。模型的物理参数在原稿里给得比较全额定 800V/2164A、1500kW额定转速 44r/min过载系数 1.5电枢电阻 0.1Ω主回路电感 3.25mH变流器放大倍数 80飞轮惯量 8.6e5。这些数决定模型的动态特性改动任何一个之前整定出来的 PI 参数都得重来。用经典工程设计方法先算一版基线得到四个增益% [Knp, KnI, Kip, KiI]分别对应转速环比例、转速环积分、电流环比例、电流环积分 PI_eng [7.3660, 1.5658, 0.8257, 30.2773];这四个数非常重要它给了 PSO 一个对照基准。整定跑完如果结果比这组还差说明要么适应度函数写歪了要么搜索边界没框住好解要么模型参数搭错了。别跳过这一步直接信 PSO 的输出。搭模型时有两个细节决定了能不能被 M 文件驱动。第一个是四个 PI 增益必须写成变量名而不是常数也就是 ASR 的比例环节填Knp、积分环节填KnIACR 同理。第二个是模型文件名要和代码里的名字对上原稿区分了double_circle.slx人工调参用和PSO_double_circle.slxPSO 驱动用这个分离很聪明避免优化过程把手工版改坏。4.2 global 变量加 sim 命令让 M 文件真正驱动 slxSimulink 模型的变量默认从 base 工作区解析M 文件函数里的局部变量它看不见。所以要么用global要么用assignin(base, ...)把值推过去。原稿选的是前者。function y PID_access_sim(para) % para 的每一行是一组 [Knp, KnI, Kip, KiI] global Knp KnI Kip KiI M size(para, 1); y zeros(M, 1); for j 1:M % 把这一组参数写进全局变量模型仿真时会读它们 Knp para(j,1); KnI para(j,2); Kip para(j,3); KiI para(j,4); % SrcWorkspace 指current让模型从当前工作区取变量 out sim(PSO_double_circle.slx, ... StopTime, 2, ... SrcWorkspace, current, ... SaveOutput, on); % 从记录信号里取转速和电流两条曲线 n_curve out.get(n_out); i_curve out.get(i_out); [n_sigma, ts_n, i_sigma, ts_i] cal_n(n_curve, i_curve); % 四项指标加权权重按工程关注度分配 y(j) 2*log(ts_n/0.5 1) 1*log(n_sigma/0.01 1) ... 1*log(ts_i/0.02 1) 1*log(i_sigma/0.2 1); end end跟 M 文件版最大的区别在性能上。PID_sim.m能一次仿 200 组因为它是纯矩阵运算而sim命令一次只能跑一个模型实例所以这里必须写 for 循环200 个粒子就是 200 次模型加载、初始化、求解、清理。单次仿真的时间从毫秒级涨到几十毫秒一代的耗时直接翻十倍以上。SrcWorkspace这个参数容易漏。默认值是base而global声明的变量在 base 工作区确实可见所以很多时候不加也能跑。但如果你的函数被parfor包起来worker 的工作区和 base 是隔离的不加这个参数或者不重新assignin模型读到的是上一次的旧值——这种 bug 最恶心结果看着合理但完全不可复现。cal_n.m负责从两条曲线里提取四个指标转速超调 σn、转速调节时间 tsn、电流超调 σi、电流调节时间 tsi。逻辑和parameter_cal.m一样只是输入从矩阵变成 timeseries取值方式要改成n_curve.Data这种。写的时候注意 Simulink 输出的时间轴不一定均匀如果是变步长求解器得先插值到均匀网格再算误差带。4.3 多指标权重改适应度函数就是改整定目标双闭环比单环有意思的地方在于指标之间是打架的。电流环压得猛启动快但电流超调大转速环压得狠稳态准但启动过程拖沓。这时候适应度函数就不只是打个分而是表个态。下面是我在这类多目标场景里常用的调整手法% 方案一偏重转速稳态质量允许电流稍微冲一点 y 3*log(ts_n/0.5 1) 2*log(n_sigma/0.01 1) ... 1*log(ts_i/0.02 1) 0.5*log(i_sigma/0.2 1); % 方案二偏重启停性能电流绝不能超限太多 y 1*log(ts_n/0.5 1) 1*log(n_sigma/0.01 1) ... 2*log(ts_i/0.02 1) 4*log(i_sigma/0.2 1);权重大小不是拍脑袋它等价于这一项偏离期望值一档我愿意用其他项偏离几档来换。系数从 1 提到 4意味着这一项在总分里的话语权翻了四倍PSO 会优先保它代价是别的指标变差。原稿里专门有一节讲在优化过程中修改参数价值权重讲的就是这件事需求变了不用重写算法改系数就行。还有一层更细的控制是给某些指标加硬约束而不是加权。比如电流超调绝对不能超过 25%超了直接判死刑if any(i_sigma 0.25) y(j) 1e6; % 直接给一个不可能被选中的大罚值 else y(j) 正常加权公式; end罚函数法的好处是简单直接缺点是如果初始种群全都违反约束PSO 就没有梯度可循等于在瞎飞。所以罚值不能设得太大以至于把差异抹平我一般用1e3 100*i_sigma这种带梯度的形式让违反程度轻的粒子还有机会活下来。4.4 和工程设计法对照确认优化没跑偏整定跑完拿出两组数放在一起看指标工程设计法基线PSO 整定结果说明转速超调 σn约 4e-4 %接近 0两者都基本无超调转速调节时间 tsn0.7875 s通常更短PSO 会优先压这一项电流超调 σi20.38 %取决于权重权重给大才能压下去电流调节时间 tsi0.0172 s通常持平电流环本身已经很快这张表的作用是验证不是炫耀。如果 PSO 的 tsn 比基线还长八成是搜索边界xmax、xmin没框住基线附近的值或者vmax太小导致粒子根本飞不到那儿。先把基线值代进适应度函数算一次分再跟 PSO 的最优分对比分都算不对就说明适应度函数本身有 bug跟算法无关。5. 避坑与排查这几个坑我逐个踩过5.1 换台电脑就报矩阵维度不匹配现象在自己机器上跑得好好的PSO.m换到另一台机器上第一代就报错位置停在X (X xmax).*xmax ...这一行。原因xmax [15,50,2]是 1×3 行向量X是 200×3 矩阵两者逐元素比较依赖 MATLAB 的隐式扩展这个特性从 R2016b 才开始有。原稿基于 2017a没问题对面机器要是 2014b 或者更老就只能报错。解决要么统一升级到 2016b 以上要么把这几个限幅改成显式扩展bsxfun(gt, X, xmax)或者干脆repmat(xmax, m, 1)一次复制好。后者多占一点内存但代码在任何版本上都能跑我一般选这个。5.2 曲线画出来了但迭代看着像卡死现象命令窗口打印了几行就停住CPU 占用很高figure 窗口一片空白等了十分钟没动静。原因PID_sim的 debug 参数传了 true每一代都把 200 条曲线画到同一个 figure 上还带hold on。200×20 等于 4000 次绘制渲染开销远大于计算开销看起来像死机实际是在画图。解决优化阶段一律debugfalse只在最后用最优参数单独跑一次仿真画图。另外原稿说的观察迭代次数是否变化确实管用在循环里加一句fprintf(iter %d, fmin %.4f\n, i, fmin)只要这行在动程序就活着。5.3 适应度出现 NaN整个种群被污染现象跑到某一代之后fmin变成 NaN之后所有代都是 NaNPg 也不再更新。原因log的输入出现零或负数。ts可能算成 0参数刚好让系统瞬间到位sigma可能算成负数欠阻尼系统的反向超调(fp-f_infty)/f_infty在 fp 小于稳态值时就是负的。log(0)是 -Inflog(负数)直接给复数再参与min比较整代就废了。解决在PID_access里加一道保护ts max(ts, 1e-6)、sigma abs(sigma) 1e-9同时对非有限值直接赋罚值。这类保护看着土但在优化循环里是必需品省得你半夜被一个 NaN 折腾两小时。5.4 调节时间恒定等于仿真总时长现象不管参数怎么变算出来的 ts 永远是 0.9995 秒接近total_t。原因有两个。一是f_infty y(:,N)取的末点还没到稳态稳态值本身偏小误差曲线全程都在误差带外max找到的越界位置永远是最靠前那一拍。二是parameter_cal的边界处理当ferr delt_err*f_infty全为 false 时max返回索引 1ts_i N-1直接退化成全程时长。解决把稳态值改成末段 10% 的均值mean(y(:, round(0.9*N):N), 2)并且在算 ts 之前加一句判断全在误差带内就直接给 0。这两处改动很小但对收敛质量影响很大因为一个永远等于总时长的指标等于没提供任何区分度。5.5 Simulink 版换成 parfor 之后结果全乱现象为了加速把PID_access_sim里的 for 改成 parfor跑完发现最优参数每次都不一样甚至同一组参数两次打分不同。原因global变量在并行池的每个 worker 里是独立的主进程改了Knpworker 里读到的还是初始值或者上一次的值。sim命令在 worker 上还需要单独加载模型重复加载会带来额外的不确定性。解决并行化时不要用 global改成把参数打包成结构体通过sim的SrcWorkspace或者set_param注入或者用Simulink.SimulationInput数组加parsim这是 R2017b 之后官方推荐的批量仿真方式每个仿真任务自带独立参数不会互相污染。改造成本比 for 循环高但省下来的时间是实打实的。6. 二次开发换对象、改约束、把 PSO 换掉三条路怎么选拿到这套代码想做自己的项目最常见的三个需求是换被控对象、加工程约束、换优化算法。三条路的改动量差别很大我按经验排个序。换对象是最简单的改动集中在一个函数里。M 文件路线下把PID_sim.m里的sys tf(0.998,[0.021,1])换成你自己的传递函数或者干脆用ss建状态空间模型其余代码一行不用动。如果对象是非线性的、带饱和或死区的那就别在 M 文件里硬写了直接上 Simulink 路线把模型搭好四个增益留成变量接PID_access_sim就行。这里有个判断标准能用tf或ss三个参数描述清楚的走 M 文件需要查表、限幅、状态机的走 slx。加约束改的是适应度函数不是算法。常见的约束有三种参数本身的上下限靠xmax、xmin卡性能指标的硬指标靠罚函数参数之间的耦合关系比如要求Kip/KiI 0.02这种得在PID_access里算完参数再判断违反就赋罚值。需要注意的是约束加得越多可行域越小初始种群的撒点范围也要跟着收否则前几代全是罚值等于白跑。换算法是最省事的因为接口已经解耦好了。整套代码只要求优化器提供输入一组参数、输出一个标量的接口PID_access就是这个接口。想换成遗传算法直接把PSO.m换成ga(PID_access, 3, [],[],[],[], xmin, xmax)连适应度函数都不用重写。想用 MATLAB 自带的particleswarm也是一行的事还能白嫖它的向量化选项和并行选项。我一般先用自带的跑一遍当对照再决定要不要手写。分享一个我自己的习惯不管换什么对象、改什么约束动手之前都先把基线参数代进适应度函数算一次分把这个数记在脚本注释里。整定跑完把最优分和基线分并排看如果最优分只比基线好一点点说明搜索空间没铺开或者代数不够如果最优分比基线还差那就是有 bug别急着相信结果。这个动作花不了一分钟但能挡掉八成看起来收敛了其实全是错的情况。另外提一个容易忽略的验证技巧把 PSO 找到的最优参数代回模型跑一次带扰动的工况比如在 0.5 秒时突加 20% 负载看转速跌落和恢复时间。优化过程只用了阶跃响应扰动响应没进适应度函数两者不一定同步变好。这一步不做交付的时候很容易被打回来。希望帮到你。本文还有配套的精品资源点击获取
