企业官网建设流程全解析

简介：面向自动控制与PID整定学习者的这份PDF文档，围绕Simulink环境下用粒子群算法（PSO）自动寻优PID的Kp、Ki、Kd展开，适合希望摆脱手工查表试凑、需要处理多输入多输出系统调参的读者。文档先梳理PSO的参数范围、种群规模、迭代次数等关键设置，以及种群初始化、速度与位置更新、适应度评估和个体/全局最优更新流程，再结合M文件与Simulink模型说明PSO优化PID的实现思路，并讨论惯性权重、学习因子c1与c2对探索和开发平衡的影响。读者可据此调整评价函数，把超调、调节时间等性能诉求转化为优化目标，理解“所需即所得”的整定方式，同时获得迭代输出、结果曲线与代码运行注意事项等排错线索。资源包共1个PDF文件，约529KB，篇幅集中、便于按目录查阅；目前已有4895人学习下载。

1. 为什么把 PSO 和 PID 摆在一起：三个参数六方向的手调困局

凌晨两点，你盯着 Simulink 里那条阶跃响应，Kp 再大一点超调就压不住，积分时间再长一点调节时间又拖成一条长尾。粒子群算法整定 PID 这件事，说到底就是把「三个参数、六个方向」的手工试凑，换成一个能按指标打分的自动搜索器。这份资源给的是 MATLAB 下的一套可跑代码：PSO 主循环、PID 仿真引擎、指标计算，以及一个双闭环调速系统的 Simulink 模型。适合已经写过 PID、被查表法和人工试凑拖过工期的人；也适合想把多输入多输出系统的「A 输出超调小、B 输出调得快」写成评价函数的人。它不是一键万能钥匙，评价函数写歪，照样翻车。

2. PSO 的 MATLAB 骨架：速度更新、惯性递减与限幅三件事

2.1 从鸟群到向量化：一次迭代里到底算了什么

粒子群算法在 1995 年被提出，和遗传算法同属群体迭代，但它的参数更少：没有交叉、没有变异，只有速度、位置、个体历史最优和种群历史最优。资源里用鸟群找食物做比喻，这个比喻对写代码的帮助有限，真正要盯住的是下面三个向量：

v = w*v + c1*rand*(P - X) + c2*rand*(ones(m,1)*Pg - X); X = X + v*dt;

第一项w*v是惯性，保留上一次飞行的方向；第二项c1*rand*(P-X)是认知项，把粒子往它自己曾经到过的最好位置拉；第三项c2*rand*(Pg-X)是社会项，把粒子往整个种群当前最好位置拉。P是 m×dimension 的矩阵，每一行是一个粒子的历史最优；Pg是 1×dimension 的行向量，用ones(m,1)*Pg扩展成 m 行，才能和X逐元素相减。

这里有一个容易被忽略的细节：资源代码里的rand没有带参数，MATLAB 返回一个标量，也就是说同一代里所有维度共用同一个随机数。更常见的写法是c1*rand(m,dimension).*(P-X)，让每个维度独立随机。两者都能收敛，但标量随机的探索方向更单一，遇到维度之间数量级差得大的参数（比如 Kp 是 0.2~15，Ki 是 0~50），容易出现某个维度被拖慢的情况。我一般会在自己改代码时换成按维度随机，资源原版保留标量写法也不影响跑通，知道差别就行。

速度更新之后是位置更新，dt是位移仿真间隔，资源里取 0.3。dt小于 1 会让搜索更平滑，但代数不变的情况下总位移变小，收敛变慢；dt大于 1 则可能一步跨过最优解。这个参数没有理论最优值，通常 0.1~0.5 之间试。

2.2 把 PSO.m 逐行拆开：初值、限幅与历史最优

资源里的PSO.m是一个无约束全局粒子群函数，输入是目标函数句柄、维度、迭代代数、种群规模、位置上下限和速度上下限。下面这版是我按原逻辑整理过的，注释补全了，索引也改了一处：

function [Pg, fmin, Savef, SaveData] = PSO(f, dimension, n, m, xmax, xmin, vmax, vmin) % f 目标函数句柄，输入 m×dimension，输出 m×1 % dimension 参数维度，PID 是 3，双闭环 PI 是 4 % n 迭代代数 % m 种群规模 % xmax/xmin 位置上下限，1×dimension % vmax/vmin 速度上下限，1×dimension Savef = zeros(n+1, 1); SaveData = zeros(m, dimension, 10); % 只留 10 代快照，用于画分布图 w = 1; c1 = 2; c2 = 2; dt = 0.3; v = zeros(m, dimension); % 初始速度为 0 X = (xmax - xmin) .* rand(m, dimension) + xmin; P = X; % 每个粒子的历史最优 last_f = f(X); [fmin, min_i] = min(last_f); Pg = X(min_i, :); % 种群历史最优 Savef(1) = fmin; N = 0; for i = 1:n v = w*v + c1*rand*(P - X) + c2*rand*(ones(m,1)*Pg - X); v = (v > vmax).*vmax + (v >= vmin & v <= vmax).*v + (v < vmin).*vmin; X = X + v*dt; X = (X > xmax).*xmax + (X >= xmin & X <= xmax).*X + (X < xmin).*xmin; new_f = f(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); Savef(i+1) = fmin; % 原代码写 i，第一代会被覆盖 if mod(i, max(1, floor(n/10))) == 0 N = N + 1; SaveData(:, :, N) = X; end % w = w - i/n*0.7*w; % 打开这行就是惯性递减 end end

几个参数的实际含义：dimension在 PID 场景是 3（Kp、Ki、Kd），双闭环是 4（转速环 Kp、Ki，电流环 Kp、Ki）；m取 50~200，太小容易早熟，太大单代仿真时间线性上升；n取 20~100，资源里 M 文件版 PID 用了 20 代、200 个粒子，双闭环模型更重，通常把m降到 30~50。

限幅那两行是逐元素写法，(v > vmax).*vmax把超上限的维度压回上限，中间项保留原值，第三项把低于下限的维度抬回下限。这种写法在vmax是标量时也能跑，但三个维度共用一个速度上限会出问题：Kp 的范围是 0.2~15，Ki 的范围是 0~50，Kd 是 0~2，统一限幅 1 对 Kp 合适，对 Ki 就太小，对 Kd 又太大。我一般会把vmax写成和维度同长的向量，比如[1, 5, 0.5]，每个维度按参数数量级单独限。

还有一个细节：Savef在原代码里从i=1开始写，而Savef(1)在循环前已经赋过初值，结果第一代的最优值被覆盖，收敛曲线少一个点。改成i+1就对了。mod(i, floor(n/10))在n<10时floor(n/10)=0，会触发除零或mod报错，加一个max(1, ...)是必要的防御。

2.3 用一个测试函数验证收敛：惯性权重到底改了什么

资源里给了一个二维测试函数：

f = @(x) (1 - x(:,1)).^2 + (x(:,1).^2 - 2*x(:,2)).^2; rng(42); % 固定随机种子，方便对比 [Pg, fmin] = PSO(f, 2, 100, 200, 40, -40, 40, -40); disp(Pg); disp(fmin);

这个函数的极小值在(1, 0.5)，理论最小值是 0。调用时位置限幅取[-40, 40]，速度限幅也取[-40, 40]，种群 200，迭代 100 代。跑完看Pg是否接近[1, 0.5]，fmin是否接近 0，这是验证 PSO 骨架有没有写错的最快办法。

真正值得对比的是惯性权重w。资源正文说「随着迭代次数减小w可以增加收敛速度」，代码里给的递减式是：

w = w - i/n*0.7*w;

这行如果直接按字面理解，会以为 100 代之后w从 1 降到 0.3。实际不是。它是乘性更新，每一代都在上一代的基础上乘一个因子(1 - 0.7*i/n)，100 代累乘下来大约是exp(-0.35) ≈ 0.70，只降到 0.7 左右，不是 0.3。想让它线性降到 0.4，应该写成：

w = 1 - (1 - 0.4)*i/n; % 从 1 线性降到 0.4

两种写法对收敛的影响不一样：乘性递减前期降得快、后期降得慢，线性递减全程均匀。资源里的测试结果说「减小惯性因子后粒子分布更稳定，几乎全部收敛到中心点」，这个结论方向是对的，但如果你照抄那行乘性公式，实际降幅比想象小，别指望它把w压到很低。

另外，w固定为 1 时，粒子在最优点附近来回振荡，部分粒子甚至反向运动，这是惯性过大导致的。实际调的时候，我一般先跑固定w=0.8看收敛，再开递减对比，别一上来就改公式。

3. 把 PID 响应折成一个数：评价函数与 M 文件仿真引擎

3.1 指标提取：稳态值、调节时间、超调与峰值时间

PSO 要优化 PID，必须先把一条响应曲线折成一个标量。资源里的parameter_cal.m做了这件事：

function [f_infty, tp, ts, sigma] = parameter_cal(y, T0, delt_err, debug) % y m×N 的响应矩阵，每行一个粒子 % T0 采样时间 % delt_err 调节时间的相对误差带，资源里取 2e-2 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; % 过阻尼无峰值 sigma = (fp - f_infty)./f_infty; % 超调量 end

这段代码有三个隐含假设，读的时候要心里有数。第一，稳态值取y(:,N)，也就是仿真最后一拍的值，如果仿真时间不够、系统还在上升，稳态值会偏小，调节时间和超调都会算歪。第二，delt_err是相对误差带，2e-2 表示稳态值 ±2% 以内算进入稳态，这个值越小，调节时间越长。第三，tp在过阻尼时被置为NaN，如果后面评价函数直接拿tp做运算，NaN会污染整代粒子，所以资源里的评价函数没有用峰值时间，只用ts和sigma。

我一般会在parameter_cal里加一道保护：f_infty = max(abs(f_infty), eps)，防止稳态值为 0 时除零。另外，如果仿真时间允许，稳态值用最后 5%~10% 的平均比取单点稳，尤其是响应尾部还有轻微波动的时候。

3.2 评价函数的对数设计：为什么 5e-2 和 1e-4 是分水岭

资源里把评价函数写成：

function y = PID_access(para) kp = para(:,1); ki = para(:,2); kd = para(:,3); [~, ~, ts, sigma] = PID_sim(kp, ki, kd, false); ts_ref = 5e-2; sigma_ref = 1e-2; % 正文公式写的是 1e-4，代码里是 1e-2 y = log(ts/ts_ref + 1) + log(sigma/sigma_ref + 1); end

评价函数越小越好，PSO 里直接取min，不需要再取负。这一点要先确认，因为有的实现把评价函数取负当适应度，如果两步都取负，方向就反了。

对数设计的好处是量纲归一和敏感度可控。对log(ts/5e-2 + 1)求导，在ts=5e-2处的斜率是1/(2*5e-2)=10，ts越小斜率越大，说明算法在小调节时间区间更敏感；ts越大斜率越小，说明大调节时间区间不会因为数值爆炸而主导整个评价。超调项同理。两个参数都小于期望值时，评价函数下降；都大于期望值时，评价函数上升；但上升速率越来越慢，下降速率越来越快。结果就是算法倾向于把两个指标都压到期望值附近，而不是牺牲一个换另一个。

这里有一个必须指出的不一致：资源正文的公式写的是sigma/1e-4，代码里是sigma/1e-2，两者对超调的惩罚力度差 100 倍。照代码跑，超调 1% 时log(1e-2/1e-2+1)=log(2)=0.69；照正文公式跑，log(1e-2/1e-4+1)=log(101)=4.62。后者会逼着算法把超调压得更狠，但也会让评价函数对超调过度敏感，调节时间被牺牲。我一般以代码为准，先跑通再决定要不要改成1e-4，改的时候只改一处，别两个文件各写一个值。

3.3 批量仿真引擎：把 M 个粒子塞进一次循环

PSO 每一代要评估 m 个粒子，如果每个粒子单独跑一次仿真，MATLAB 的循环开销很大。资源里的PID_sim.m把 m 组 PID 参数拼成列向量，一次循环推进所有粒子：

function [f_infty, tp, ts, sigma] = PID_sim(kp, ki, kd, debug) Tt = 5e-4; % 仿真步长 0.5 ms T0 = 1e-2; % 控制周期 10 ms 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; 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=5e-4是仿真步长，T0=1e-2是控制周期，两者差 20 倍。代码里积分项用T0*ki，微分项用(1-alp)/T0*kd，但控制量u(:,k)每个仿真步都刷新，没有做零阶保持。这意味着它更接近连续控制器的离散化，而不是真实的数字控制器。如果你要严格模拟 10 ms 控制周期，得在u后面加一个零阶保持，或者把控制量更新放到mod(k, T0/Tt)==0的条件里。资源里的目的是整定参数趋势，不是做控制器代码生成，所以这个近似可以接受，但心里要清楚。

还有一个顺序问题：y(:,k)计算时用了u(:,k)，而u(:,k)在下一行才被赋值。因为u预分配为全零，第一次进入循环时u(:,k)=0；之后每次循环，当前拍的u(:,k)都还没算出来，读到的是预分配的 0。也就是说num(1)*u(:,k)这一项在整个仿真中恒为 0。对于tf(0.998,[0.021,1])离散化后num(1)=0的一阶惯性对象，这没有影响；但如果换成num(1)≠0的对象，响应会偏慢。我一般会把控制量先算出来，再推进对象，或者至少确认num(1)是不是 0。

r=1500是阶跃幅值，对应电机转速给定。sys=tf(0.998,[0.021,1])是一阶惯性对象，时间常数 0.021 s，稳态增益 0.998。这些数值来自资源里的直流电机简化模型，换对象时整段都要改，不能只改r。

3.4 入口脚本与第一次跑通

资源说从set_para.m开始运行，所有 M 文件和 slx 放在同一个文件夹。不同版本的资源入口脚本名字可能不一样，有的叫main.m，有的叫set_para.m，以压缩包里的为准。入口里通常设三个量：

dimension = 3; n = 20; m = 200; [Pg, fmin] = PSO(dimension, n, m);

第一次跑不要直接上m=200, n=20。先把m改成 20、n改成 5，确认命令行每代都有打印、三幅图能出来，再放大规模。资源里说「每次迭代都会打印出来，如果运行时间过长怀疑程序死机可以观察迭代次数是否有变化」，这句话是实话：M 文件版一次迭代要算 200 条 1 秒、步长 0.5 ms 的响应，单代耗时可能到几十秒，看起来像卡死，其实在算。

跑完会输出三幅图：响应曲线随迭代次数的变化、调节时间随迭代次数的变化、超调量随迭代次数的变化。第一幅图如果线条叠在一起看不清，在图上局部放大，或者把debug打开只画最终一代。资源里提到「响应曲线可能要在左上角局部放大才能获得该图像」，说的就是这个。

4. Simulink 双闭环调速系统联调：从 slx 到 PSO 的调用链

4.1 模型改造：PI 参数变成可写变量

M 文件版的 PID 引擎适合验证算法，真正复杂的系统还是要回到 Simulink。资源里用双闭环直流调速系统做例子：外环转速环 ASR 是 PI 调节器，内环电流环 ACR 也是 PI 调节器。启动时电流先到过载限幅，电机以最大加速度加速，转速到给定后电流迅速降下来，维持负载转矩所需的电流。

模型里要整定四个参数：Knp、KnI、Kip、KiI。改造的关键不是搭模型，而是把四个 Gain 模块里的数字换成变量名。双击 Gain 模块，把Gain一栏写成Knp，不要写7.366。四个模块都改完，另存为PSO_double_circle.slx，和 PSO 的 M 文件放在同一目录。

资源里给的电机参数如下，搭模型时按这个填：

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