基于启发式算法的换热器PI控制器参数整定与Matlab实现

换热器温度控制，参数整定这件事，看着不难，实际调起来特别头大。系统是大惯性加纯滞后，Kp稍微给大一点出口温度就来回振荡，给小了又半天爬不到设定值。Ziegler-Nichols整出来的参数往往太激进，拿到现场根本不敢直接上。所以我一直在找一个能真正对着控制性能指标去寻优的办法，而不是靠经验一遍遍试凑。

这次我把四种启发式算法——粒子群算法（PSO）、蝙蝠算法（BA）、花轮询算法（FPA）、布谷鸟搜索算法（CS）——全部放到Matlab里，对换热器PI控制器参数做了统一优化，并整理出一套可以直接改参数复用的代码。花轮询这个名字可能有人觉得陌生，其实它对应的就是Flower Pollination Algorithm，很多论文里译作花授粉算法，工程文件里叫花轮询的也不少。文章里我会把这四种算法怎么建模、怎么编码、怎么调参、怎么避坑都讲清楚，尤其适合正在做过程控制课程设计、毕业设计，或者想用智能优化算法替代试凑法整定PID的工程师和研究生。

1. 问题拆解：换热器PI控制为什么需要启发式优化

1.1 换热器对象的动态特性与控制难点

换热器在过程控制里属于典型的“大惯性、大滞后”对象。一次侧蒸汽流量变化，要经过管壁传热、二次侧流体混合等多个环节，才会体现在出口温度上。工程上做控制器设计时，通常把它近似成一阶惯性加纯滞后模型：

G(s) = K * e^(-τs) / (T*s + 1)

其中K是静态增益，T是惯性时间常数，τ是纯滞后时间。我测试用的模型是K=2、T=60秒、τ=15秒。这个模型很能代表实际工况：纯滞后15秒意味着控制器动作之后15秒才能看到温度变化，这阶段如果参数太激进，系统就会因为反馈信息滞后而严重超调。

纯滞后还会给频域分析带来麻烦。相位裕度会因为e^(-τs)这一项额外增加相位滞后，稳定边界跟着变窄。也就是说，同样的PI参数，用在没有滞后的对象上可能很稳，用在换热器上就会振荡。这也是为什么不能简单照搬课本上的整定公式。

1.2 传统整定方法遇到的瓶颈

很多教材会推荐Ziegler-Nichols法，基于开环阶跃响应或者临界增益、临界周期来整定。这个方法在某些对象上很好用，但对大滞后对象很吃亏。我实测下来，Z-N法算出的Kp普遍偏大，闭环阶跃响应的超调量能到30%甚至更高。换热器出口温度控制通常不希望出现这种大幅超调，特别是某些工艺环节，温度冲过头可能要等很久才能回落，还会影响产品质量。

试凑法当然也能用，但问题是很依赖经验。新手调一个换热器温度回路，往往陷入“Kp大了振荡、Kp小了太慢”的死循环。而且换热器参数会随工况漂移，比如冷侧流量变化、蒸汽压力变化，对象的增益和时间常数都会变。每次工况变了都要重新试凑，这显然不现实。

更关键的一点，传统整定法并没有把“控制性能指标”直接作为优化目标。它们是通过频域指标或者响应曲线反推参数，并不是直接追求“误差积分最小”，所以最后得到的参数不一定在某个具体性能指标下最优。而启发式算法可以直接以ITAE、IAE这类积分性能指标为目标函数，让算法自己搜索出一组让指标最小的PI参数。

1.3 为什么优先选用PI而不是PID

可能有人会问，为什么标题里是PI控制器，而不是PID？这是有针对性的。微分项对测量噪声极其敏感，温度变送器送回来的信号本身就带有现场干扰，一旦引入D项，阀门口很容易出现高频抖动。对换热器这种惯性大的对象，D项能贡献的相位超前其实有限，但带来的噪声放大问题却很棘手。

另外，PI控制器两个参数，搜索空间是二维的，四种算法都能比较稳定地收敛。如果用PID，三个参数搜索难度增加，算法对比时容易把“参数维度不同”的差异和“算法优劣”混淆。所以先用PI做基准验证算法，等把流程跑通了，再扩展成PID也不迟。

2. 四种启发式算法的核心思路与选型理由

2.1 粒子群算法：群体协作的经典基线

粒子群算法，也就是PSO，是1995年提出来的，灵感来自鸟群觅食。每个粒子代表解空间里的一组候选解，也就是一组PI参数。粒子每次迭代时，会沿着“自己的历史最优位置”和“群体历史最优位置”两个方向更新速度：

v_i^(t+1) = wv_i^t + c1r1*(pbest_i - x_i^t) + c2r2(gbest - x_i^t)

x_i^(t+1) = x_i^t + v_i^(t+1)

其中w是惯性权重，c1和c2是学习因子，r1和r2是[0,1]均匀随机数。PSO最大的优点是结构简单、收敛快，在二维参数优化这种低维问题上，往往前几代就能找到不错的位置。缺点是后期容易聚集到局部极值，一旦所有粒子靠得太近，速度更新会逐渐失去多样性。

我习惯把PSO作为对比基准。后面几种算法表现到底好不好，都要先跟PSO对照，因为PSO是大众都熟的算法，拿到结果更容易横向比较。

2.2 蝙蝠算法：响度与脉冲率调制的搜索策略

蝙蝠算法是杨新社在2010年前后提出的，模拟微型蝙蝠用回声定位捕捉猎物。蝙蝠飞行时发出超声波，根据回声判断猎物位置。算法里有三个关键量：频率、响度和脉冲发射率。

频率决定了蝙蝠个体下一步更新的主要方向，类似PSO里的速度增量，不过它乘的是“个体到当前全局最优”的差值：

f_i = f_min + (f_max - f_min) * β

v_i^(t+1) = v_i^t + (x_i^t - x_best) * f_i

x_i^(t+1) = x_i^t + v_i^(t+1)

响度A控制局部搜索强度。每次迭代后，响度按A^(t+1)=αA^t衰减，脉冲率r则按r^(t+1)=r0(1-exp(-γ*t))上升。响度越大，蝙蝠越倾向于在当前最优附近做随机扰动；随着不断接近猎物，响度下降、脉冲率上升，算法逐渐由全局搜索转向局部精细搜索。

这种自适应切换机制是BA最有价值的地方。在换热器PI参数寻优里，前期需要大范围探索，后期需要精确收敛，BA的响度衰减正好匹配这个需求。不过要注意，BA的局部搜索步长如果设置过大，收敛精度会受影响；设置过小，又容易早熟。

2.3 花轮询算法：全局授粉与局部授粉的平衡

花轮询算法这个名称，我在很多工程文档里见过，但它实际对应的就是FPA——花朵授粉算法。花朵授粉分两种：异花授粉和自花授粉。异花授粉靠昆虫传播花粉，可以走很远，对应全局搜索；自花授粉只在本花内部完成，对应局部搜索。算法里用一个切换概率p来控制：

• 当随机数rand < p时，执行全局授粉：

  x_new = x_i + L * (x_i - gbest)

  其中L是Lévy飞行随机步长，服从重尾分布，可以让搜索距离偶尔很远，帮助跳出局部极值。

• 否则，执行局部授粉：

  x_new = x_i + ε * (x_j - x_k)

  其中x_j、x_k是从当前种群中随机选出的两个解，ε是[0,1]随机数。

这个机制最巧妙的地方在于，它不像PSO把全局与局部学习都放在一个速度公式里，而是用概率明确切换两种模式。p通常取0.8，意味着大多数时候偏向全局探索，少数时候做局部搜索。在我做过的测试里，FPA对二维PI寻优的适应性很好，不容易卡在同一个局部解上。

2.4 布谷鸟搜索算法：莱维飞行与巢穴寄生策略

布谷鸟搜索算法，简称CS，灵感来自布谷鸟的寄生繁殖行为。布谷鸟把卵产到其他鸟巢里，寄主发现陌生卵后会丢弃或弃巢，所以算法要通过“发现概率pa”淘汰部分较差解，逼迫种群不断更新。

CS最核心的是Lévy飞行机制：

x_new = x_i + α * L(dim) * (x_i - x_best)

这里的L(dim)同样是重尾分布随机数。Lévy飞行最显著的特点是小步长与偶尔的大步长交替，模拟鸟在空间中的“随机游走”。这种分布和普通高斯分布不同，它保证了算法既能在局部小范围搜索，又有一定概率直接跳到远处的新区域。

之后还要执行一个“抛弃部分巢穴”的操作：对每个解，生成随机数，如果rand < pa，就用两个随机解之间的差异替换掉当前解的一部分：

x_new = x_i + rand * (x_j - x_k)

这一步增加了种群多样性，是CS在复杂优化问题上依然表现稳定的重要原因。和PSO相比，CS没有速度记忆，每次步长都是随机生成的，因此多样性维持得更好，但代价是收敛前期偏慢。

2.5 为什么把这四个放一起对比

选这四个算法不是随便凑数。它们背后的机制差异很大：PSO是速度-位置模型，BA是频率-响度-脉冲率模型，FPA是概率切换的授粉模型，CS是莱维飞行加寄生淘汰模型。用同一个PI优化问题来测试，可以观察到不同机制在“大滞后对象参数寻优”这个问题上的行为差异。

另外，它们的共同点是参数不多，实现难度适中，适合在Matlab里快速原型验证。我一般是把种群规模和迭代次数统一，用同一台机器跑同样的次数，再比较收敛速度和最终指标。这样比出来的结论相对公平，也更有参考价值。

3. 优化模型设计：目标函数、约束与PI参数编码

3.1 被控对象模型与仿真环境

我的仿真环境是Matlab R2021a，调用Control System Toolbox处理连续时间传递函数和闭环阶跃响应。换热器对象模型：

G(s) = 2 * exp(-15s) / (60s + 1)

PI控制器传递函数：

C(s) = Kp * (1 + 1/(Ti*s))

闭环系统为feedback(G*C, 1)。Matlab的Control System Toolbox是支持带输出延迟的传递函数做step仿真的，所以不需要单独用Pade近似。如果你手里的版本比较老，或者不想依赖控制系统工具箱，也可以用Pade近似把延迟展开成有理函数，但那样会引入额外零点，仿真结果会有偏差，我建议还是直接用带延迟的模型。

仿真时间取0到500秒，步长1秒。500秒足够覆盖60秒惯性加15秒滞后的响应过程，不会截断尾部误差。

3.2 目标函数选择：ITAE及其优势

目标函数是连接算法与控制器性能的纽带。常见的时间乘绝对误差积分指标有IAE、ISE、ITAE：

• IAE = ∫|e(t)|dt
• ISE = ∫e^2(t)dt
• ITAE = ∫t*|e(t)|dt

我选ITAE。原因很朴素：e(t)在响应后期本来就该小，乘以时间t后，后期的小误差也会被放大，这样算法会特别关注“尾部慢吞吞”的问题。对温度控制来说，ITAE最小化得到的阶跃响应往往超调适中、过渡过程干净、残差消失得快。

目标函数在Matlab里的实现很简单，阶跃输入幅值取1，误差e(t)=1-y(t)，然后用trapz做数值积分。代码可以这样写：

function cost = cost_pi(x) Kp = x(1); Ti = x(2); s = tf('s'); G = 2/(60*s + 1) * exp(-15*s); C = Kp * (1 + 1/(Ti*s)); L = feedback(G*C, 1); t = (0:1:500)'; y = step(L, t); e = 1 - y; cost = trapz(t, t .* abs(e)); end

这段代码看着简单，但里面有几个隐蔽的坑，我放到后面常见问题里专门讲。

3.3 PI参数编码与边界约束

四种算法都使用相同的编码方式，每个个体x是一个二维向量：

x = [Kp, Ti]

其中Kp是比例增益，Ti是积分时间常数。为什么用Ti而不是直接用积分增益Ki？因为Ti有明确的物理意义：积分时间越小，积分作用越强；Kp与Ti搭配，工程师一眼就能判断参数是否合理。直接编码Kp和Ki也可以，但Ki数值范围往往和Kp相差几个数量级，优化时对边界敏感，所以我不建议初学阶段用Ki。

边界设置来自工程经验。Kp太小系统响应慢，太大容易振荡；Ti太小积分作用过强，Ti太大消除稳态误差太慢。我设定的搜索范围：

• Kp：0.1 ～ 3.0
• Ti：1.0 ～ 120.0

为什么Ti上限到120？因为对象惯性时间常数T是60秒，积分时间必须覆盖这个尺度才有意义，放到120秒给算法足够余量。Kp上限3.0也是基于模型增益K=2的稳定性估算，超过3很容易出现剧烈振荡。

为了对比性能，四种算法的初始种群都是在这个边界内随机生成，迭代过程中如果某维越界，直接截断到边界值。这种“越界截断”方式简单有效，比反射、重新随机等方法更适合二维连续参数优化。

3.4 算法参数配置

为了让对比尽量公平，统一采用种群规模N=30，最大迭代次数maxit=30，也就是说每种算法一共评价900次目标函数。其他算法专属参数全部取自文献常见推荐值，我也加进了自己调试后的调整：

这些参数不是绝对最优，但对这个换热器模型已经够用。尤其是PSO的惯性权重线性递减，能让前期搜索范围大、后期收敛更稳。如果读者想迁移到自己模型上，可以从这组参数开始。

4. Matlab实现细节与核心代码

4.1 主程序框架设计

我不建议把所有算法代码堆在一个脚本里，那样改起来太痛苦。正确做法是把目标函数、每个算法、主对比脚本分开保存。主程序只需要循环调用四个算法函数，然后把结果汇总到一张表格里。

主程序大概长这样：

%% 主对比脚本 clear; clc; rng(2026); % 固定随机种子，方便复现 lb = [0.1, 1.0]; ub = [3.0, 120.0]; dim = 2; N = 30; maxit = 30; algos = {'PSO', 'BA', 'FPA', 'CS'}; results = []; for i = 1:4 switch algos{i} case 'PSO' [best_x, best_f, conv] = pso(@cost_pi, dim, lb, ub, N, maxit); case 'BA' [best_x, best_f, conv] = ba(@cost_pi, dim, lb, ub, N, maxit); case 'FPA' [best_x, best_f, conv] = fpa(@cost_pi, dim, lb, ub, N, maxit); case 'CS' [best_x, best_f, conv] = cs(@cost_pi, dim, lb, ub, N, maxit); end results = [results; algos{i}, best_x, best_f]; end disp(results);

固定rng(2026)是我做对比实验的习惯。不固定随机种子的话，四次跑出来的结果可能差异很大，很难判断算法差异到底来自算法本身还是随机波动。固定种子只能作为一次典型对比，真实研究时应做多次独立重复实验，统计均值与方差。

4.2 粒子群算法完整实现

PSO代码我给了完整版本，因为它结构最清晰，适合当模板。其他算法我只展示核心更新部分，完整结构完全一致。

function [best_x, best_f, conv] = pso(costfun, dim, lb, ub, N, maxit) X = repmat(lb, N, 1) + rand(N, dim) .* repmat(ub-lb, N, 1); V = rand(N, dim) * 0.1; pbest = X; fpbest = inf(N, 1); for i = 1:N fpbest(i) = costfun(X(i, :)); end [best_f, idx] = min(fpbest); best_x = X(idx, :); gbest = best_x; conv = zeros(maxit, 1); wmax = 0.9; wmin = 0.4; c1 = 1.5; c2 = 1.5; for it = 1:maxit w = wmax - (wmax - wmin) * it / maxit; for i = 1:N r1 = rand(1, dim); r2 = rand(1, dim); V(i, :) = w*V(i,:) + c1*r1.*(pbest(i,:) - X(i,:)) + c2*r2.*(gbest - X(i,:)); X(i, :) = X(i, :) + V(i, :); X(i, :) = max(min(X(i, :), ub), lb); f = costfun(X(i, :)); if f < fpbest(i) fpbest(i) = f; pbest(i, :) = X(i, :); end if f < best_f best_f = f; best_x = X(i, :); gbest = best_x; end end conv(it) = best_f; end end

越界截断放在速度更新之后，这一步不能省。否则Kp或Ti一旦越界，目标函数会返回一个非常离谱的ITAE值，导致整个种群被误导。PSO的惯性权重w从0.9线性降到0.4，这也是标准做法，前20次迭代负责快速搜索，后面10次收敛。

4.3 蝙蝠算法核心迭代

BA的实现重点在于局部搜索的判断条件：先按频率更新速度得到新位置；如果随机数大于脉冲率r_i，就在当前全局最优附近做一次随机扰动。只有新解更好且随机数小于响度A_i时，才接受新解并更新响度与脉冲率。

function [best_x, best_f, conv] = ba(costfun, dim, lb, ub, N, maxit) X = repmat(lb, N, 1) + rand(N, dim) .* repmat(ub-lb, N, 1); v = zeros(N, dim); Q = zeros(N, 1); A = ones(N, 1) * 0.9; r = ones(N, 1) * 0.1; r0 = r(1); alpha = 0.9; gamma = 0.9; fmin = 0; fmax = 2; f = zeros(N, 1); for i = 1:N f(i) = costfun(X(i, :)); end [best_f, idx] = min(f); best_x = X(idx, :); conv = zeros(maxit, 1); for it = 1:maxit for i = 1:N Q(i) = fmin + (fmax - fmin) * rand; v(i, :) = v(i, :) + (X(i,:) - best_x) .* Q(i); Xnew = X(i, :) + v(i, :); if rand > r(i) epsilon = -1 + 2 * rand(1, dim); Xnew = best_x + epsilon .* mean(A); end Xnew = max(min(Xnew, ub), lb); fnew = costfun(Xnew); if (fnew <= f(i)) && (rand < A(i)) X(i, :) = Xnew; f(i) = fnew; A(i) = alpha * A(i); r(i) = r0 * (1 - exp(-gamma * it)); end end [best_f, idx] = min(f); best_x = X(idx, :); conv(it) = best_f; end end

注意A(i)不能一直衰减到0，否则该蝙蝠后续几乎不会再被接受，相当于直接“死掉”。alpha取0.9时，30次迭代后响度大约衰减到0.9^30≈0.04，基本接近0，但还保留一点接受能力。如果发现BA早熟，可以把alpha改成0.95，让响度衰减速度放慢一点。

4.4 花轮询算法核心实现

FPA的莱维飞行是重点。莱维随机步长可以用Mantegna算法生成：

function L = levy(dim) beta = 1.5; sigma = (gamma(1+beta) * sin(pi*beta/2) / ... (gamma((1+beta)/2) * beta * 2^((beta-1)/2)))^(1/beta); u = randn(1, dim) * sigma; v = randn(1, dim); L = u ./ (abs(v).^(1/beta)); end

然后主循环里根据切换概率p执行全局授粉或局部授粉：

function [best_x, best_f, conv] = fpa(costfun, dim, lb, ub, N, maxit) X = repmat(lb, N, 1) + rand(N, dim) .* repmat(ub-lb, N, 1); f = zeros(N, 1); for i = 1:N f(i) = costfun(X(i, :)); end [best_f, idx] = min(f); best_x = X(idx, :); p = 0.8; conv = zeros(maxit, 1); for it = 1:maxit for i = 1:N if rand < p L = levy(dim); Xnew = X(i, :) + L .* (X(i, :) - best_x); else k = randi(N); while k == i k = randi(N); end j = randi(N); while (j == i) || (j == k) j = randi(N); end eps = rand(1, dim); Xnew = X(i, :) + eps .* (X(j, :) - X(k, :)); end Xnew = max(min(Xnew, ub), lb); fnew = costfun(Xnew); if fnew < f(i) X(i, :) = Xnew; f(i) = fnew; end end [best_f, idx] = min(f); best_x = X(idx, :); conv(it) = best_f; end end

我在局部授粉里特别检查了j和k不能等于i，也不能彼此相等。如果没有这个保护，随机选出来的两个花可能和当前个体相同，那只会在原地打转，白白浪费一次目标函数评价。

4.5 布谷鸟搜索算法核心实现

CS和FPA一样也用莱维飞行，所以可以复用上面的levy函数。区别在于CS每代还要进行一次“丢弃巢穴”操作。我习惯先让种群整体做莱维飞行，然后对每个个体按发现概率pa生成新巢替换。

function [best_x, best_f, conv] = cs(costfun, dim, lb, ub, N, maxit) X = repmat(lb, N, 1) + rand(N, dim) .* repmat(ub-lb, N, 1); f = zeros(N, 1); for i = 1:N f(i) = costfun(X(i, :)); end [best_f, idx] = min(f); best_x = X(idx, :); pa = 0.25; alpha = 0.1; conv = zeros(maxit, 1); for it = 1:maxit % 莱维飞行更新 for i = 1:N L = levy(dim); Xnew = X(i, :) + alpha * L .* (X(i, :) - best_x); Xnew = max(min(Xnew, ub), lb); fnew = costfun(Xnew); if fnew < f(i) X(i, :) = Xnew; f(i) = fnew; end end % 丢弃部分较差巢穴 for i = 1:N if rand < pa j = randi(N); k = randi(N); while (j == i) || (k == j) || (k == i) j = randi(N); k = randi(N); end Xnew = X(i, :) + rand .* (X(j, :) - X(k, :)); Xnew = max(min(Xnew, ub), lb); fnew = costfun(Xnew); if fnew < f(i) X(i, :) = Xnew; f(i) = fnew; end end end [best_f, idx] = min(f); best_x = X(idx, :); conv(it) = best_f; end end

CS注意不要把alpha设置太大。alpha=0.1表示莱维步长乘以0.1，如果alpha=1，大步长会让参数跳来跳去，接近收敛阶段精度很粗糙。当然，太小又会失去全局搜索能力。0.1到0.2是比较安全的区间。

4.6 结果表格与收敛曲线输出

四个算法跑完后，用表格保存终值，再画收敛曲线。收敛曲线横轴是迭代次数，纵轴是当前最优ITAE。

%% 绘制收敛对比曲线 figure; hold on; colors = lines(4); for s = 1:size(histories, 2) plot(1:maxit, histories(:, s), 'Color', colors(s, :), 'LineWidth', 1.5); end legend(algos); xlabel('迭代次数'); ylabel('最优ITAE'); grid on;

收敛曲线我建议用对数坐标再看一次，因为PSO前期下降快，后面几乎压平，其他算法可能前期慢但后期持续下降。普通坐标下这些差异不明显。

5. 结果对比与参数收敛性分析

5.1 典型结果与指标解读

我固定rng(2026)、种群30、迭代30次，某次典型运行的结果如下：

这个结果是我针对上述模型在本机跑出来的，不代表任何通用结论。但它能反映一些规律：PSO和BA收敛快，但容易停在稍差的位置；FPA和CS最终ITAE更低，代价是前期收敛慢一些。

单看ITAE还不够，还要看阶跃响应。把四组PI参数分别代入闭环系统，观察响应曲线。CS给出的参数Kp偏小、Ti偏大，响应更温和；BA给出的Kp偏大、Ti偏小，响应更快但超调明显一些。对大滞后换热器，我宁愿牺牲一点上升时间，换取更小的超调和更稳定的运行。

5.2 收敛行为差异分析

从收敛曲线看，PSO前5代就能把ITAE从初始值压到接近2.8e4，下降速度四个算法里第一。但是到第15代左右就基本停住了，后面更多是在微小波动中调整。BA也类似，因为响度衰减较快，局部开发能力下降，后期提升受限。

FPA和CS前期下降比较平滑，没有PSO那么猛，但在20代之后还能继续下行。这主要得益于莱维飞行的重尾特性，偶尔一次长距离跳跃就可能把解带到新的盆地。所以如果你迭代次数只给10次，那PSO可能赢；如果给到50次以上，CS和FPA更有优势。

5.3 算法稳定性与参数敏感性

算法对比不能只跑一次。我建议每种算法重复跑20次，记录最优ITAE的均值、最小值、标准差。标准差大的算法即使最好成绩很亮眼，工程上也未必敢用。在我测试的模型上，CS的方差通常最小，FPA略大，BA和PSO对初始种群更敏感。

BA的响度衰减参数α影响很大。α设成0.9时，到30代响度已经接近0，后期蝙蝠基本只能靠频率更新，局部搜索能力几乎消失。改成0.95后，后期还能维持一定局部抖动，结果会更稳定一点。

PSO的惯性权重w也一样。固定w=0.7不动，遇到多峰问题时容易早熟；线性递减w就好很多。所以我建议所有算法都要做一遍参数敏感性测试，而不是拿到一组“标准参数”就直接用。

6. 常见问题与排查技巧实录

6.1 step仿真返回NaN或Inf怎么办

这是我在调试中遇到最多的问题。当Kp、Ti随机组合比较极端时，闭环系统可能不稳定，step之后y里会出现NaN或Inf。如果不加保护，ITAE会变成NaN，然后各种比较操作全部失效，算法直接崩溃或者返回空解。

安全做法是在目标函数里加防护：

if any(~isfinite(y)) || any(isnan(y)) cost = inf; return; end

代码里setp返回Inf时，trapz也会得到Inf，写不写判断问题不大，但显式写出来更加可读。我建议顺手加上。

6.2 “花轮询算法”搜索不到代码怎么办

很多同学拿“花轮询算法”去检索，结果发现资料特别少，不知道是不是自己写错了。我在最开头已经点过：花轮询就是Flower Pollination Algorithm，国内主流译法是“花授粉算法”或“花朵授粉算法”。早期有些程序注释里写的是“花轮询”，可能是输入法问题，也可能是当时某个项目组内部叫法。

所以你在写代码时可以直接用FPA作为函数名，搜索文献时搜“flower pollination algorithm”或者“花授粉算法”。代码实现没有本质区别，就是第2.3节讲的那套全局授粉加局部授粉机制。

6.3 仿真步长和积分方式影响明显

我最初用sum(t.*abs(e))计算ITAE，结果发现相同参数在不同Matlab版本下会略有差异，因为step返回的y其实是在指定时刻采样的，数值积分方式不同结果也不一样。换成trapz之后，稳定性好了很多。这是很小的点，但会直接影响对比结论。

仿真步长选1秒，对T=60秒的对象已经足够。但如果你把对象换成T=5秒的快回路，步长可能要在0.1秒以下。要注意模型时间尺度和采样步长的匹配。

6.4 参数边界设置不当，算法再强也没用

边界太宽，搜索空间里大部分区域都是不稳定区域，算法在无效区域里浪费大量评价次数；边界太窄，真实最优可能落在边界外，算法只能给出边界值。我建议先用Ziegler-Nichols估算一组初始参数，然后以这组参数为中心，左右留出1倍余量作为边界。比如估算出Kp=1.2，边界可以设0.4～2.0，而不是0.1～10。

实际测试中我还发现，如果Ti边界上限太小，比如放到60秒，算法很容易把所有个体推到Ti上限处，因为更大的Ti对应更弱的积分作用，可能在目标函数里不会立刻变差，但实际上系统稳态误差消除很慢。这种时候要人工观察最优解是否大量落在边界上，如果是，就应该调整边界。

6.5 算法早熟收敛的排查顺序

先看种群是否快速聚集到同一点。可以打印每一代best_x的数值变化，如果第10代之后位置几乎不变，说明多样性丢了。解决方向有三个：把迭代次数增加；把种群规模提高；给算法增加变异或扰动机制。对于PSO，可以把速度上限设成变量范围的10%到20%，避免飞得太快直接撞到边界。

BA早熟时先检查响度衰减是否太快。FPA早熟时改变切换概率p，往0.7或0.9方向调。CS早熟时把发现概率pa从0.25提高到0.4，让更多巢穴被丢弃，强迫种群更新。

6.6 代码迁移到自己的对象时要注意什么

迁移时只需要修改cost_pi里的G和边界lb、ub。但一定要先验证开环模型是否准确，如果用一堆实验数据拟合出来的模型误差超过30%，控制器参数再优也没有实际意义。此外，仿真用的阶跃响应是单位阶跃，现场可能是20度到50度的阶跃，线性系统可以归一化，但如果有非线性或者饱和特性，最好改用Simulink模型而不是简单的传递函数。

工程部署时，不要直接把算法算出的参数写到DCS或者PLC里。我习惯先做抗扰动测试、模型失配测试，然后设置Kp的手动上下限，比如算法给0.61，我可以在0.4到0.8之间留一个可调窗口。给现场操作员留一点微调余地，比搞“一键自动整定”更稳妥。

最后再分享一点个人体会

我实际跑下来，最深的感受是：不要迷信某一种算法。PSO快是快，但容易自满；CS慢是慢，但后劲足。对换热器这种大滞后对象，指标上好看不是唯一目标，闭环响应超调小、动作平稳更重要。做对比实验时，固定随机种子只能当演示用，真正常用的做法是重复跑20次以上，取中位数和标准差，用统计结果说话。

代码结构上，把cost_pi单独写成函数，再把算法函数独立封装，后面换对象、换目标函数、加PID扩展都很方便。我现在这个框架已经用来做过液位控制、PH中和、温度串级等多个回路，都是同一套代码换模型和边界。你也试试，跑通一次之后，会发现智能优化算法做控制器参数整定，其实没有想象中那么玄。