上个月有个做电加热设备的朋友找我诉苦说厂里那台热处理炉的PID参数一直是老师傅靠手感摸出来的一批料换一种工况就得重新试试一次大半天。我听完就乐了这都什么年代了参数整定完全可以交给算法自己去搜。我给他搭了一套蜣螂算法DBO优化PID参数的仿真环境m代码驱动Simulink自动跑把Kp、Ki、Kd三个值当成优化变量几分钟就能出一组看着还不错的参数。这篇文章把整套思路和可运行的代码分享出来适合正在做控制仿真、又对参数整定头疼的朋友。先说清楚一件事这不是什么高深的前沿课题而是把一类已经被验证过的智能优化算法套到控制工程最经典的参数整定问题上。思路本身不复杂真正麻烦的是m代码和Simulink之间的数据交互、多次仿真的稳定性、以及适应度函数怎么设计才合理。这些恰恰是论文里不会细写、但实际动手一定会遇到的东西。1. PID参数为什么难整定试凑、Z-N法与智能搜索的对比1.1 三个参数互相牵扯人工试凑的效率太低PID控制器之所以能统治工业控制这么多年是因为它足够简单比例项管响应速度积分项管稳态误差微分项管超调抑制。问题是这三个参数不是独立工作的。你把Kp调大响应快了但超调跟着上来为了压超调加Kd结果系统对噪声敏感了为了消稳态误差加Ki又把之前的稳定性节奏全打乱。这种强耦合关系让试凑法变成一门手艺活老师傅能调好是因为他脑子里积累了大量的调参手感而这种手感很难复制。尤其遇到大纯滞后、强非线性、或者工况经常变化的对象人工试凑基本就是在碰运气。比如我朋友那台热处理炉加热惯性大温度反馈滞后明显用试凑法一个参数组合下去等温度曲线跑完要等十几分钟试几组就一上午过去了。1.2 Z-N法虽然快但在非理想对象上表现偏保守很多人说整定用Ziegler-Nichols不就行了这话对一半。Z-N法分为反应曲线法和临界比例度法本质是通过对象的阶跃响应或临界振荡数据用经验公式直接算出三组参数。好处是快几分钟能出结果坏处是它针对的是理想线性小滞后对象一旦对象带纯延迟、非线性、执行机构饱和Z-N法给出的参数往往超调很大甚至系统根本稳不住。以我后面测试用的对象为例一阶惯性加1.5秒纯延迟用Z-N反应曲线法算出来的参数做阶跃响应超调能到20%以上。这种参数丢到现场设备来回震荡谁看了都头大。1.3 智能算法整定的本质把调参问题翻译成优化问题智能算法整定的思路和前面完全不一样。它不跟你讲先调P再调I最后调D的手感而是把PID参数直接定义成优化变量把控制效果定义成目标函数然后让算法在这个三维空间里自动搜索。你只需要告诉它三个约束参数范围、仿真时长、性能指标怎么算。剩下的事就是让算法反复调用仿真模型不断试新的参数组合最后返回一组让目标函数最小的Kp、Ki、Kd。这个过程的价值在于不需要人工干预能处理非线性、大延迟、带饱和的对象而且每次工况变化后只要重新跑一遍优化就行。缺点也很明显——计算量大因为每一组参数都要跑一次完整仿真。但配合Simulink加速模式几分钟完成几十轮迭代是完全可以接受的。选型上为什么用DBO而不是更常见的GA或者PSO两个原因。第一GA涉及到编码解码和一堆超参数交叉率、变异率、选择策略你用它对PID参数这种连续实值问题反而绕远路PSO在低维问题上容易早熟飞着飞着就堆到局部最优出不来。第二DBO是2022年提出的新算法它的四种行为机制天然把全局探索和局部开发分开设计滚球负责跑远找新区域产卵和偷窃负责在好解周边精细搜索这种分工在PID参数这种多峰、非线性的目标函数上比较稳而且实现起来并不比PSO复杂。2. 蜣螂算法在模仿什么四种行为与数学表达2.1 滚球与跳舞全局探索的两种随机策略蜣螂俗称屎壳郎有个很有意思的行为它会把粪球滚成球体然后推着走并且依靠月光或者偏振光导航尽量让粪球沿直线滚动。如果遇到障碍物它没法直线前进就会爬到粪球上跳舞——也就是绕着一个轴转圈重新确定方向。算法把这套行为抽象成了两个更新公式。正常滚球时个体位置更新为x_i(t1) x_i(t) alpha * k * x_i(t-1) b * |x_i(t) - X_w|其中t是当前迭代次数alpha是方向系数随机取1或-1表示蜣螂偏离原来方向的角度k是偏转系数用来控制历史位置对当前位置的影响b是常数X_w是当前种群里的全局最差位置。公式里的|x_i(t) - X_w|可以理解为光照强度的变化量蜣螂倾向于远离光源弱也就是适应度差的区域这样种群整体会朝着更好的方向推进。遇到障碍时跳舞更新为x_i(t1) x_i(t) tan(theta) * |x_i(t) - x_i(t-1)|其中theta是[0, pi]范围内的随机角度。注意当theta等于0、pi/2或者pi时tan(theta)分别为0、无穷大、0这时候位置不更新模拟蜣螂跳舞后没有找到新方向的随机性。这个机制的作用是给种群注入随机扰动避免所有个体都往同一个方向集中从而减少陷入局部最优的概率。2.2 产卵、觅食、偷窃局部开发的三重保障比滚球更妙的是后面三种行为。雌蜣螂会把卵产在粪球里然后埋在土壤中一个安全区域。算法把这个区域建模成动态收缩的边界Lb* max(X* * (1 - R), Lb) Ub* min(X* * (1 R), Ub) R 1 - t / T_maxX*是当前局部最优个体Lb和Ub是问题的变量上下界。R随迭代次数从1线性降到0意味着产卵区域一开始很大后期逐渐收缩到局部最优附近。产卵更新为x_i(t1) X* b1 * (x_i(t) - Lb*) b2 * (x_i(t) - Ub*)b1、b2是两个随机向量。这个公式保证新个体在局部最优周边产生实现精细搜索。觅食行为描述的是小蜣螂从巢穴出来找食物的过程它也有一个动态收缩的最佳觅食区域边界定义和产卵类似只是中心换成了全局最优X^bLbb max(Xb * (1 - R), Lb) Ubb min(Xb * (1 R), Ub) x_i(t1) x_i(t) C1 * (x_i(t) - Lbb) C2 * (x_i(t) - Ubb)注意这里C1服从柯西分布C2服从标准正态分布目的是让搜索步伐既有随机跳跃又有局部微调。偷窃行为则模拟某些蜣螂不自己推粪球而是直接去偷别的蜣螂的粪球。这类个体围绕全局最优Xb和局部最优X*做随机游走x_i(t1) Xb S * g * (|x_i(t) - X*| |x_i(t) - Xb|)S是常数g是标准正态分布随机向量。偷窃个体的优势是它始终记住全局最优在哪里同时利用当前个体与两个最优的差异来产生新位置搜索步长能自适应变化。2.3 角色分配与迭代边界实现时最容易糊涂的地方实际编码的时候一个绕不开的问题是种群里的个体怎么分配行为角色我查阅了多份公开的复现代码发现不同实现差别还挺大。有的按20%/40%/30%/10%分有的干脆四等分。我自己测试下来一套比较稳的划分方式是行为角色占比主要作用滚球跳舞30%全局探索防止早熟产卵20%围绕局部最优精细搜索觅食20%围绕全局最优搜索偷窃30%强化全局最优附近的挖掘这个比例不是定死的。如果发现算法收敛太慢可以增加滚球比例如果发现陷入局部最优就增加偷窃和觅食的比例。重要的是理解每种角色在干什么而不是死记一个数字。每次迭代都要做一步把种群按适应度从小到大排序排在前面的个体作为产卵和觅食角色它们本身质量较高负责开发排在后面的作为滚球和偷窃角色它们需要更多探索。这就是为什么代码里要先排序再更新顺序不能乱。3. 联合仿真的数据链路m代码指挥Simulink跑仿真3.1 为什么选工作区变量传参数而不是其他方式DBO优化PID最核心的工程问题是算法在m脚本里跑Simulink模型怎么快速拿到当前的Kp、Ki、Kd并把仿真结果返回给m脚本计算适应度常用的方式有三种。第一种是把PID参数直接写成工作区变量Simulink模型引用变量名m脚本用assignin把新参数塞进base工作区再调用sim()仿真。第二种是用set_param去改PID模块里P、I、D三个参数的值。第三种是把优化算法整个放进Simulink的MATLAB Function模块里在模型内部完成循环。我强烈建议用第一种。set_param的方式每次改参数都要触发模型编译几百次迭代跑下来速度感人MATLAB Function的方式调试不方便中间变量看不到而且一旦算法里出个数组维度错误排查起来很痛苦。工作区变量方式最大的优势是PID模块里填的表达式比如Kp在每次仿真开始前自动从base工作区取值m脚本只需要在sim()之前写一行assignin(base,Kp,newKp)干净利落。这里有个容易踩的坑如果你把DBO主程序写成函数函数内部定义的局部变量Kp不会自动传给SimulinkSimulink只认base工作区。所以必须显式assignin到base否则sim()会报未定义变量Kp。我在第六章会详细展开。3.2 Simulink模型的搭建细节模型结构不复杂下面这条链就是全部Step - Sum(-) - PID Controller - Saturation - Transfer Fcn(5/(2s1)) - Transport Delay(1.5s) - To Workspace(Scope) ^ | |___________________________ 负反馈 _______________________________________|模块选型上注意几点。PID Controller模块的P、I、D参数直接填Kp、Ki、Kd不要加引号不要内联数值Simulink会在仿真启动时从base工作区解析这三个名字。Saturation限幅模块模拟执行机构的物理限制我习惯把上下限设成[-3, 3]这样优化出来的参数更贴近工程实际不会出现控制器输出无限大的理想情况。被控对象用Transfer Fcn加Transport Delay组合比用LTI System模块更直观。为了把数据导回工作区在对象输出端接一个To Workspace模块变量名设成youtSave format选Array再单独用一个To Workspace记录时间tout或者在模型配置参数里勾选输出时间。我偏好前者因为显式、可控。模型搭好后一定要用CtrlD编译一次确认没有错误再把模型文件保存为dbopid_sim.slx。主程序里通过load_system(dbopid_sim)预加载模型这样可以避免每次sim()调用都重新解析模型文件能省下大量时间。3.3 适应度函数设计ITAE加超调惩罚适应度函数是整定效果的决定性因素。我用的是经典ITAE指标也就是时间加权绝对误差积分ITAE integral( t * |e(t)| , dt )其中e(t)是设定值与被控输出的偏差。ITAE相比IAE和ISE的好处是它对稳态阶段的小误差更敏感能有效抑制长时间的低频振荡整定出来的系统超调小、阻尼特性好这是控制领域几十年实践下来的共识。但只优化ITAE还不够在带纯延迟的对象上可能收敛出的参数让系统超调10%以上而ITAE值依然不高。所以我在适应度里加了一项超调惩罚J ITAE * (1 1.5 * Overshoot/100)超调量Overshoot越大惩罚越重。这个1.5的系数可以根据你的工程偏好调整如果系统绝对不允许超调就把系数加大到3甚至5如果允许10%以内超调换更快响应就调小一点。另外还要处理一种特殊情况DBO随机生成的参数组合可能让系统发散仿真输出变成NaN或者Inf。这种个体不能直接参与排序否则会污染整个种群的比较。我会在适应度函数里显式判断一旦发现yout里有NaN或者Inf直接返回一个很大的数比如1e10。4. 可运行的DBO-PID整套m代码4.1 主程序框架与初始化下面这套代码我按MATLAB R2023a的语法写的大部分版本都能兼容。先用rng固定随机种子保证每次跑出来的结果可复现方便对比调参效果。%% DBO优化PID参数主程序 clear; clc; close all; rng(20240401); % 固定随机种子方便复现 % 加载Simulink模型 model dbopid_sim; load_system(model); % 参数搜索范围根据对象特性设定 lb [0, 0, 0]; % Kp, Ki, Kd下限 ub [3, 1.5, 2]; % Kp, Ki, Kd上限 % DBO参数 N 30; % 种群规模 T 50; % 迭代次数 dim 3; % 优化变量维度 % 角色数量划分产卵20%、觅食20%、偷窃30%、滚球30% N_brood floor(0.2 * N); N_forage floor(0.2 * N); N_steal floor(0.3 * N); N_roll N - N_brood - N_forage - N_steal; % 种群初始化 X repmat(lb, N, 1) rand(N, dim) .* repmat(ub - lb, N, 1); fitness zeros(N, 1); % 先给base工作区赋默认PID参数防止sim()首次调用报错 assignin(base, Kp, 1); assignin(base, Ki, 0.1); assignin(base, Kd, 0.1); % 初始适应度评估 for i 1:N fitness(i) dbopid_fitness(X(i, :), model); end % 按适应度排序适应度小的排前面 [fitness, order] sort(fitness); X X(order, :); % 记录全局最优与局部最优 Xbest X(1, :); fbest fitness(1); Xlocal X(1, :); Xworst X(end, :); convergence zeros(T, 1); % 保存每轮最优适应度用于画收敛曲线4.2 四种行为更新代码主循环里每次迭代分两步先按角色更新位置再重新评估和排序。我这里把四种行为都完整写了滚球行为需要用到上一代的位置prevX所以在每次迭代开始时保存一份。%% DBO主循环 for t 1:T R 1 - t / T; % 动态收缩系数随迭代从1降到0 prevX X; % 保存上一代位置供滚球/跳舞公式使用 for i 1:N if i N_brood % ----- 产卵行为在局部最优Xlocal附近动态收缩区域内生成新个体 ----- LbStar max(Xlocal .* (1 - R), lb); UbStar min(Xlocal .* (1 R), ub); b1 rand(1, dim); b2 rand(1, dim); X(i, :) Xlocal b1 .* (X(i, :) - LbStar) b2 .* (X(i, :) - UbStar); elseif i N_brood N_forage % ----- 觅食行为围绕全局最优Xbest收缩搜索 ----- Lbb max(Xbest .* (1 - R), lb); Ubb min(Xbest .* (1 R), ub); C1 tan(pi * (rand(1, dim) - 0.5)); % 柯西分布抽样 C2 randn(1, dim); % 标准正态分布 X(i, :) X(i, :) C1 .* (X(i, :) - Lbb) C2 .* (X(i, :) - Ubb); elseif i N_brood N_forage N_steal % ----- 偷窃行为围绕全局最优与局部最优的差做随机游走 ----- S 2; % 偷窃步长常数 g randn(1, dim); X(i, :) Xbest S * g .* (abs(X(i, :) - Xbest) abs(X(i, :) - Xlocal)); else % ----- 滚球跳舞行为 ----- if rand 0.9 % 正常滚球按概率取方向系数alpha为1或-1 if rand 0.5 alpha 1; else alpha -1; end k_coef 0.1; % 偏转系数 b_coef 0.3; % 光照强度影响系数 deltaX abs(X(i, :) - Xworst); X(i, :) X(i, :) alpha * k_coef * prevX(i, :) b_coef * deltaX; else % 跳舞随机角度重新定向theta为0或pi/2或pi时不更新 theta rand * pi; if abs(theta) 1e-10 || abs(theta - pi/2) 1e-10 || abs(theta - pi) 1e-10 X(i, :) X(i, :); else X(i, :) X(i, :) tan(theta) .* abs(X(i, :) - prevX(i, :)); end end end % 边界钳制确保参数不超出搜索范围 X(i, :) max(X(i, :), lb); X(i, :) min(X(i, :), ub); end % 重新评估所有个体 for i 1:N fitness(i) dbopid_fitness(X(i, :), model); end % 排序更新全局最优与局部最优 [fitness, order] sort(fitness); X X(order, :); if fitness(1) fbest fbest fitness(1); Xbest X(1, :); end Xlocal X(1, :); Xworst X(end, :); convergence(t) fbest; fprintf(iter %d/%d, fbest %.4f, [Kp Ki Kd] [%.3f %.3f %.3f]\n, ... t, T, fbest, Xbest(1), Xbest(2), Xbest(3)); end % 输出最终结果 fprintf(\nDBO优化完成\n); fprintf(最优参数: Kp %.4f, Ki %.4f, Kd %.4f\n, Xbest(1), Xbest(2), Xbest(3)); fprintf(最优适应度: %.4f\n, fbest);4.3 仿真评估函数适应度函数单独放到一个dbopid_fitness.m文件里主程序用循环调用它。function f dbopid_fitness(x, model) % 将当前PID参数写入base工作区供Simulink模型引用 assignin(base, Kp, x(1)); assignin(base, Ki, x(2)); assignin(base, Kd, x(3)); % 运行仿真固定仿真时长 try simOut sim(model, StopTime, 20, ReturnWorkspaceOutputs, on); catch f 1e10; % 仿真报错直接给最大值 return; end % 从仿真输出中提取时间和被控量 tout simOut.get(tout); yout simOut.get(yout); t tout(:); y yout(:); % 防发散处理 if any(isnan(y)) || any(isinf(y)) f 1e10; return; end % 计算ITAE e 1 - y; % 单位阶跃设定 ITAE trapz(t, t .* abs(e)); % 计算超调量 OS (max(y) - 1) / 1 * 100; if OS 0 OS 0; end % 适应度 ITAE * (1 1.5 * 超调百分比) f ITAE * (1 1.5 * OS / 100); end这里有个细节simOut.get(yout)返回的是模型里To Workspace模块导出的数据。如果你的To Workspace模块Save format设置了Timeseries而不是Array那读取方式要改成simOut.yout.getElement(yout).Values.Data。建议在模型里统一把Save format设成Array代码最省事。5. 实测对比一个纯延迟被控对象上的整定结果5.1 测试对象与仿真条件为了验证这套东西的效果我用一个典型的过程控制对象做测试一阶惯性加纯延迟传递函数为G(s) 5 / (2s 1) * e^(-1.5s)。这个对象模拟的是带传输延迟的加热腔变量之间有明显的滞后对PID参数很敏感适合用来对比整定效果。仿真条件设置如下仿真时长20秒单位阶跃给定执行器饱和限幅[-3, 3]求解器用变步长ode45最大步长限制0.05秒。DBO种群规模30迭代次数50搜索范围Kp∈[0,3]Ki∈[0,1.5]Kd∈[0,2]。为了对比我同时跑了Z-N整定法和标准PSO算法。Z-N法直接按反应曲线计算PSO的种群和迭代次数与DBO保持一致。5.2 DBO、PSO与Z-N法的阶跃响应对比在我机器上的典型结果如下表整定方法KpKiKd超调量调节时间ITAEZ-N法0.320.1070.2424.6%9.1s7.83PSO0.870.160.337.8%5.6s2.67DBO1.080.190.552.1%3.8s1.32DBO优化出来的参数把超调压到了2%以内调节时间也只有Z-N法的一半左右。这个结果说明两点一是这套联合仿真链路确实能跑通二是DBO在低维连续优化问题上表现不差至少在这个被测对象上明显优于Z-N和PSO。收敛曲线方面DBO大致在10代左右就有明显的适应度下降到20代之后基本趋于平稳而PSO到20代还在缓慢下降最终收敛值比DBO高出一截。不过我不建议把这个结果过度泛化换个被控对象、换个边界范围相对优劣可能会有变化。但至少对这类带纯延迟的过程对象DBO的表现是有竞争力的。5.3 参数配置的几条经验从多次测试来看有三条经验比较值得分享。第一搜索边界不要拍脑袋定。先用Z-N法粗算一组参数然后以这组参数为中心向上下各扩展2到3倍作为搜索范围这样算法一开始就在合理的区域内探索收敛快而且不容易跑偏。比如Z-N算出Kp0.32我就把Kp上限定到3留了接近10倍的余地既不放过可能的好解也不会因为界太宽导致算法前期疯狂试探无效区域。第二种群规模和迭代次数不是越大越好。我试过N50、T100结果比N30、T50好不到哪里去但耗时翻了三倍也试过N10、T20结果不太稳定有时会收敛到明显较差的局部最优。对于三维参数问题N20到30、T40到60是一个性价比很高的区间。第三固定随机种子极其重要。DBO本身有随机性如果不固定种子每次跑出来的最优参数会有差别你很难判断某个参数改动到底是算法改进了还是随机波动导致的。固定种子之后每次跑的结果完全一致调试和参数对比都方便得多。6. 联合仿真翻车记录我把最常见的坑一次说清6.1 sim()总是报参数未定义工作区变量的时机问题第一次尝试联合仿真的人几乎都会遇到这个报错调用sim()时提示模型里的Kp未定义。原因前面提过——sim()在base工作区解析模型变量而你如果把DBO主程序写成了函数你在函数里定义的局部变量KpSimulink根本看不到。解决办法分两步。第一步在主程序初始化时先手动assignin一次默认值确保模型第一次仿真前工作区里有Kp、Ki、Kd这三个变量。第二步在每次个体评估之前用assignin更新这三个变量。注意是在sim()之前更新不是在sim()之后。我见过有人把顺序写反结果模型永远用的是上一代参数适应度全算错了收敛曲线乱成一团。另外还要注意变量名冲突。如果你的base工作区已经有一个叫K的变量或者另一个脚本里残留了Ki的历史值有可能导致模型引用了错误变量。建议在每次优化的最开头加一行evalin(base, clear Kp Ki Kd)清掉历史残留。6.2 适应度函数里的发散处理没有兜底的算法迟早翻车DBO在搜索过程中会随机生成一些极端参数比如Kp特别大、Kd特别小这时候闭环系统可能直接发散。仿真输出yout里出现NaN或者Inf如果不管直接拿去算trapz轻则适应度算成NaN导致排序全乱重则整个循环崩溃。所以适应度函数里一定要做两件事。第一检测yout中是否有NaN/Inf有就直接返回1e10把这个个体判死刑。第二用try-catch把sim()包起来因为有些参数组合可能让求解器直接报错不捕获的话整个优化过程会中断。这两个兜底逻辑是联合仿真稳定运行的前提我一开始没加跑了几次之后程序动不动就中断血泪教训。还有个容易被忽略的点发散个体返回的1e10虽然保证了排序不会错但如果所有个体都发散了那种群就全报废了。出现这种情况优先检查lb/ub是不是设得太大或者被控对象模型本身就有问题。6.3 反复仿真越跑越慢求解器与模型配置建议DBO要跑N×T次仿真N30、T50就是1500次。如果你发现程序越跑越慢大概率不是算法问题而是仿真配置问题。最影响速度的几件事第一模型里挂了Scope或者仪表模块每次仿真都要刷新波形窗口一多直接卡成幻灯片。优化的模型里一律不要放Scope或者放完先删掉用To Workspace把数据导出来再统一画图。第二变步长求解器在某些参数下步长被压得很小求解时间暴涨。可以试试把最大步长限制设成0.05或0.1或者干脆用固定步长离散求解器。第三每次sim()调用都在反复编译模型用load_system预加载模型可以减少这部分开销。如果模型复杂度上来了建议把模型切成加速模式在Simulink界面把仿真模式从Normal改成Accelerator或者在m代码里用set_param(model, SimulationMode, accelerator)。加速模式下第一次仿真会编译后面每次仿真都复用编译产物1500次仿真跑下来能快好几倍。6.4 把优化结果落回工程现场的最后一公里DBO给出一组最优参数不代表可以直接扔到设备上用。仿真模型和真实对象之间永远有偏差主要是未建模动态、噪声特性、执行机构磨损这些因素。我的习惯是把DBO的搜索结果当成高质量初值再在真实设备上做一次局部的梯度式细调比如用单纯的爬山法或者Nelder-Mead只在小范围微调Kp、Ki、Kd很快就能收敛到现场工况下的最优。这种先用智能算法粗搜再用传统方法精调的组合方案比纯靠人在现场试要稳得多。DBO负责在全局范围内找到一片好的参数区域局部精调负责克服模型失配带来的偏差。两段式流程下来一组实际可用的PID参数通常半小时内就能搞定而不是像以前那样调一整天。另外整定结束后记得把优化过程画出来看一遍。把convergence数组画成曲线如果你发现适应度在最后十几代还在明显下降说明迭代次数不够下次加大T如果前几代就快速收敛然后一路平坦说明问题维度低、算法余量充足可以大胆减小N节省下次运行时间。这些针对性的调整比自己瞎猜参数靠谱得多。就拿我朋友那台热处理炉来说这套流程跑完出来的参数放到真实设备上只做了两次微调就交付了。他在厂里感慨说原来调参这活儿也能工业化。我倒觉得这事儿本来就不难只是大多数人被手动试凑这四个字给框住了。换个思路把参数整定交给算法去搜剩下的事就简单了。