做单区域负荷频率控制LFC研究的朋友应该都有过这种体验在Simulink里搭好模型手工调PID参数调到头疼好不容易在这组负荷扰动下表现不错换个扰动步长又开始振荡想写论文又拿不出系统性的参数整定依据只能含糊地说“经过多次仿真调试确定”。我最初做这个课题时也是这个状态直到把灰狼优化算法GWO引入进来作为PID参数整定的搜索引擎整个流程才真正变得可复现、可解释、可写进论文。这篇文章记录的就是我完整跑通“GWO整定单区域负荷频率控制PID”的一次实践包含模型搭建、算法代码、Simulink联合仿真的每一步以及那些迭代了上百次才摸清的坑。1. 单区域LFC的模型长什么样PID整定难在哪1.1 频率偏差产生的物理逻辑先回到最基础的问题负荷频率控制到底在控什么。电力系统里发电出力要和负荷实时平衡一旦负荷增加而发电机出力没有跟上系统的频率就会下降反之负荷减少频率会升高。频率是全局性物理量直接反映发电与负荷的平衡状态偏离额定值比如50Hz会影响电机转速、变压器磁通和用户设备的正常工作。单区域负荷频率控制的“单区域”指的是忽略区域间联络线功率交换只研究一个独立控制区域。这个区域里有一个或多个发电机组共同承担频率调节任务。负荷突然增加 ΔPL 后系统频率先按惯性下降随后调速器通过一次调频动作调节进汽量/进水量来增加机械功率抬高频率。但仅有一次调频还不够调速器是按照调差系数R实现的本质上是有差调节频率不可能完全恢复到50Hz。要彻底消除稳态频差必须引入二次调频信号——也就是控制器给出的 ΔPc叠加到调速器的设定点上。这个二次控制器的角色在绝大多数研究里就是PID。1.2 我在Simulink里用的传递函数组合单区域LFC模型的简化结构非常经典由调速器、汽轮机、发电机-负荷环节和调差反馈组成。我用的传递函数和参数如下表环节传递函数参数调速器Gg(s) 1 / (Tg·s 1)Tg 0.08 s汽轮机非再热式Gt(s) 1 / (Tt·s 1)Tt 0.3 s发电机-负荷Gp(s) 1 / (M·s D)M 10 p.u.·sD 1 p.u.调差系数1/RR 0.05这里所有变量都用标幺值p.u.表示频率偏差 Δf 也是标幺值。有人习惯把发电机环节写成 Kp/(Tp·s 1)效果等价只是系数映射问题。我建议初学者直接按 M、D 来搭物理含义更清楚M 越大代表转子惯量越大频率下降越慢D 越大代表负荷阻尼越强频率越不容易偏离。这个模型在Simulink里搭起来就是几个传递函数模块的串联和反馈没有多少花哨的地方。但要注意一点调差反馈符号必须正确。调速器输入是 ΔPc 减去 Δf/R也就是一次调频是负反馈频率下降时调速器要开大阀门增加进汽。我第一次搭的时候把这条反馈的符号弄反了仿出来的频率直接发散排查了好久。1.3 传统整定方法为何不够用PID整定方法说起来一堆但真正做LFC时你会发现传统方法都有明显的别扭之处。试凑法就不提了完全看经验而且LFC模型里有惯性、延迟、反馈耦合三个参数互相影响在Simulink里反复改参数跑仿真非常浪费时间。Ziegler-Nichols频域整定法对这类高阶系统误差较大整出来的参数往往过激超调大而且它本质上针对的是带有延迟的一阶/二阶模型和LFC这种包含调速器、汽轮机、发电机多个惯性环节的系统匹配得并不好。更关键的一点是LFC关注的是动态性能与鲁棒性。PID参数不仅要让额定负荷扰动下的频率偏差响应快、超调小还要在负荷扰动幅值变化、甚至模型参数变化时依然可用。传统整定法很难把这些需求统一写进一个指标里而智能优化算法天然可以——只要把想优化的目标写成适应度函数剩下的搜索过程交给优化器。GWO之所以适合做这件事原因在于它原理简单、参数极少、不依赖梯度信息对于“黑箱式”的Simulink仿真模型非常友好。2. GWO算法原理拆解狼群行为与PID参数搜索的对应关系2.1 灰狼等级制度怎么变成寻优框架灰狼优化算法是Mirjalili等人在2014年提出的群体智能算法灵感来自灰狼群体的社会等级和捕食行为。灰狼群体分成四个等级α、β、δ、ω。α是头狼负责决策β是二把手辅助αδ负责侦察、放哨ω是最底层跟随前面三只狼移动。在GWO里每一个“狼”就是搜索空间中的一个候选解。对于我们的问题一个狼的位置就是一个三维向量 [Kp, Ki, Kd]代表一组PID参数。每只狼都要在Simulink里跑一次仿真得到一个ITAE适应度值。适应度值最小的三只狼被标记为α、β、δ它们相当于当前迭代中已知的最好、次好、第三好的PID参数组合其余狼ω根据这三只头狼的位置更新自己的位置逐步向最优区域收敛。这个过程很像团队协作不是只跟着最优秀的那个方案走而是同时参考前三名的方向避免被某个偶然的局部最优带偏。2.2 位置更新公式与收敛因子的作用GWO位置更新的数学原理不复杂。对任意一只狼 X它到猎物即最优解的距离用下式估计D |C · X_head - X(t)|然后用这个距离更新自己的位置X(t1) X_head - A · D其中 A 和 C 是两个关键系数A 2·a·r1 - a C 2·r2r1 和 r2 是 [0,1] 之间的随机数a 是从2线性递减到0的收敛因子。A 的绝对值决定了狼是扩大搜索范围|A| 1探索还是收缩逼近|A| 1开发这个机制和很多优化算法里的“前期全局搜索、后期局部精搜”策略异曲同工。C 的作用是给猎物位置增加一个随机权重相当于在搜索中引入扰动防止算法过早陷入局部最优。三只头狼分别按上面的公式产生三个候选位置 X1、X2、X3最终位置取三者的平均值X(t1) (X1 X2 X3) / 3这个“三狼平均”的策略是GWO的精华。它可以被理解为最优的PID参数大概率位于α、β、δ三个已知较优参数连线的几何中心附近但又不完全等于它们留出了继续探索的空间。2.3 适应度函数选ITAE而不选ISE的理由适应度函数是整个优化过程的“指挥棒”。你用它评价哪组PID参数好GWO就会朝着让它最小化的方向搜索所以目标函数的设计直接决定最终参数的行为。LFC里最常用的几个误差积分指标指标公式特点IAE∫ |Δf| dt对误差绝对值积分简单但不对时间加权ISE∫ Δf² dt对大误差惩罚很大但对长期小误差不敏感ITAE∫ t·|Δf| dt对时间加权后期误差影响越来越大我最终选了ITAE。原因很实际在负荷频率控制里我们不仅希望频率偏差小更希望它尽快回到零。ITAE会给“拖尾”的误差更大的惩罚所以优化出来的PID参数通常响应更快、调节时间更短不会出现为了压低超调而磨磨蹭蹭收敛的情况。ISE虽然抑制大偏差效果好但容易让优化器忽略那些低频小振幅的振荡尾巴。如果你在写论文时想让控制器兼顾控制能量可以把指标改成 J w1·ITAE(Δf) w2·∫ΔPc²dt用权重系数w1、w2权衡动态响应和控制力度。我建议先不加惩罚项把基本流程跑通再扩展。2.4 搜索空间与初始种群设置GWO的搜索空间就是PID参数的取值范围。这看起来很自由但很影响最终效果。范围设太大大量参数组合会让Simulink仿真发散适应度函数只能返回一个惩罚大数浪费计算资源范围设太小可能把最优参数卡在边界外面。我的做法分两步先用一组粗略范围比如 Kp∈[0,5]Ki∈[0,2]Kd∈[0,1]手动跑几次仿真观察Kp超过多少仿真开始发散、Ki多大时系统开始振荡然后把优化范围收缩到一个“稳定且合理”的区间。以我上面的模型为例最终用的范围是Kp ∈ [0.1, 2.0]Ki ∈ [0.01, 1.0]Kd ∈ [0, 0.5]初始种群建议用拉丁超立方采样lhsdesign替代纯随机它能让初始点更均匀地铺满搜索空间减少初始分布不均导致的早熟。如果电脑上没装统计工具箱用 rand 生成也完全可以。3. Simulink模型搭建从模块到可调参数3.1 模型拓扑与关键模块参数模型文件我命名为 LFC_single_area.slx不建议用中文名也不建议放在带中文路径的目录下否则 sim() 调用容易出幺蛾子。模型内部的信号连接顺序如下。先看被控对象部分三个传递函数模块串联。调速器前有一个求和点 Sum2计算 ΔPc - Δf/R调速器输出进入汽轮机汽轮机输出是机械功率变化量 ΔPm。接着用 Sum3 计算 ΔPm - ΔPL输入到 1/(M·sD) 环节输出就是频率偏差 Δf。Δf 引出三条支路去 Scope 和 To Workspace用于观察和记录数据经过 Gain(-1) 变成误差信号 e -Δf送入 PID Controller经过 Gain(1/R) 接入 Sum2 的负号端口实现一次调频反馈。这里分成两条负反馈路径是LFC模型的关键一条是一次调频的调差反馈一条是PID二次调频。前者是系统固有机制后者才是我们要优化的控制器。PID Controller 模块我建议用并联形式Parallel form比例增益设置为工作区变量 Kp积分增益 Ki微分增益 Kd。微分项一定要用带滤波的形式默认的D filter coefficient N100 就够否则理想微分会对仿真噪声非常敏感。3.2 参数动态注入方式GWO迭代过程中每一只狼都是一组新的 Kp、Ki、Kd怎么把这些值传进Simulink模型是关键问题。很多人第一反应是用 assignin(base, Kp, Kp) 往基础工作区塞变量模型里直接引用变量名即可。这种方式最简单但有两个隐患一是并行计算时多个仿真相互干扰二是如果仿真中途报错基础工作区可能残留脏数据。更推荐的做法是用 Simulink.SimulationInput 的 setVariable 方法把参数注入到这一次仿真的临时上下文中。这样既不污染基础工作区也方便后续用 parsim 做并行优化。代码结构后面专门讲。要注意的是模型中只要出现变量名 Kp、Ki、Kd那么在手动仿真之前工作区里必须存在这些变量否则模块参数解析会报错。最简单的方法是先手动命令行赋值 Kp0.5; Ki0.3; Kd0.1; 再打开模型。3.3 仿真配置与数据采集接口为了配合GWO循环里反复调用 sim()仿真参数按下面设置仿真停止时间50s。LFC系统的动态过程一般在10-20秒内基本结束50秒足够让收敛尾巴也计入ITAE。求解器变步长 ode45相对误差默认 1e-3 就可以。优化过程不需要超级精度而且每次仿真时间太长会让整个优化跑得非常慢。数据采集从Δf信号引一个 To Workspace 模块变量名 f_data保存格式选择 Array用 Clock 模块接另一个 To Workspace变量名 t_data同样用 Array。用 Array 格式而不是 Timeseries是为了在适应度函数里直接得到普通数组省去拆数据对象的麻烦。如果你选择 Timeseries后面取数要多写两行代码不是不行只是没必要。搭建完成后先不要跑优化手动输入一组大概率稳定的参数仿真一次确认频率响应形状正常——出现阶跃后频率下降、然后回升到0整个波形没有发散、没有剧烈振荡再进入下一步。4. GWO整定PID的Matlab实现主程序与适应度函数4.1 主程序框架与初始化整个优化工程包含两个文件GWO_LFC_main.m 是主程序负责GWO的初始化与迭代LFC_Fitness.m 是适应度函数负责调用Simulink模型并计算ITAE。主程序开头先清空环境并固定随机种子保证结果可复现。这一点做研究的人一定要养成习惯否则每次跑出的最优参数都不一样论文里根本没法写。GWO_LFC_main.m 核心框架clc; clear; close all; rng(0); mdl LFC_single_area; load_system(mdl); N 30; % 种群规模 MaxIt 40; % 最大迭代次数 dim 3; lb [0.1, 0.01, 0]; % Kp, Ki, Kd 下界 ub [2.0, 1.0, 0.5]; % 上界 X lhsdesign(N, dim) .* (ub - lb) lb; fit zeros(N, 1); Alpha_pos zeros(1, dim); Alpha_fit inf; Beta_pos zeros(1, dim); Beta_fit inf; Delta_pos zeros(1, dim); Delta_fit inf; Conv zeros(MaxIt, 1); gbest zeros(1, dim); gbest_fit inf;Alpha、Beta、Delta分别存储三只头狼的位置和适应度。gbest是全局历史最优解配合后面的精英保留策略使用。4.2 GWO迭代更新核心代码主循环分成三块计算每个个体的适应度并更新三只头狼执行精英保留更新所有狼的位置。for t 1:MaxIt for i 1:N X(i,:) max(X(i,:), lb); X(i,:) min(X(i,:), ub); fit(i) LFC_Fitness(X(i,:), mdl); if fit(i) Alpha_fit Alpha_fit fit(i); Alpha_pos X(i,:); elseif fit(i) Beta_fit Beta_fit fit(i); Beta_pos X(i,:); elseif fit(i) Delta_fit Delta_fit fit(i); Delta_pos X(i,:); end end if Alpha_fit gbest_fit gbest_fit Alpha_fit; gbest Alpha_pos; end X(1,:) gbest; % 精英保留 Conv(t) gbest_fit; a 2 - 2 * t / MaxIt; for i 1:N if i 1, continue; end for j 1:dim r1 rand; r2 rand; A 2*a*r1 - a; C 2*r2; D_alpha abs(C * Alpha_pos(j) - X(i,j)); X1 Alpha_pos(j) - A * D_alpha; r1 rand; r2 rand; A 2*a*r1 - a; C 2*r2; D_beta abs(C * Beta_pos(j) - X(i,j)); X2 Beta_pos(j) - A * D_beta; r1 rand; r2 rand; A 2*a*r1 - a; C 2*r2; D_delta abs(C * Delta_pos(j) - X(i,j)); X3 Delta_pos(j) - A * D_delta; X(i,j) (X1 X2 X3) / 3; end end end fprintf(最优PID参数: Kp%.4f, Ki%.4f, Kd%.4f\n, gbest); fprintf(最优ITAE: %.6f\n, gbest_fit); save(GWO_results.mat, gbest, gbest_fit, Conv);我加了一个精英保留直接把历史最优个体放到第1个位置并且不让它参与后续位置更新。为什么这么做标准GWO用三只狼的平均位置生成新解连α狼自己也可能会被“平均掉”导致当前最优解在下一轮就丢了。对LFC这种适应度函数有随机性和噪声的优化问题精英保留几乎是必须的。4.3 适应度函数与sim()调用细节LFC_Fitness.m 是所有仿真发生的地方。这里有几个必须注意的细节每次仿真前要把三个PID参数注入到当前仿真上下文仿真结束后要稳定地取出时间和频率数据处理发散和无解情况。function J LFC_Fitness(x, mdl) simIn Simulink.SimulationInput(mdl); simIn simIn.setVariable(Kp, x(1)); simIn simIn.setVariable(Ki, x(2)); simIn simIn.setVariable(Kd, x(3)); try simOut sim(simIn); t simOut.t_data; f simOut.f_data; f f(:); t t(:); if any(isnan(f)) || any(isinf(f)) || max(abs(f)) 1e3 J 1e6; else J trapz(t, t .* abs(f)); end catch J 1e6; end endITAE用 trapz(t, t.*abs(f)) 做数值积分比手写 for 循环精确也比简单的 sum(f)*dt 更适应变步长求解器。发散惩罚值设为1e6是一个常用技巧。仿真发散时数据里会出现NaN或Inf此时返回一个很大的适应度GWO会自动避开这批参数。要注意的是如果惩罚值不够大比如只有1那么一些轻微发散但恰好ITAE比较小的解反而会被当作最优导致最终参数不可用。还有一点如果模型路径或模型名有中文或者Simulink模型没有加载sim(simIn)会报错catch块会兜住并返回惩罚值。但这会掩盖问题让你误以为是参数不好。建议优化前先手动 load_system 一次确认模型能打开、能手动仿真再跑自动优化。4.4 用parsim把优化速度提升数倍如果你只跑40代×30只狼1200次仿真串行执行可能要好几个小时。我强烈建议把适应度计算改成 parsim 并行版本尤其是在模型不大、单次仿真只要一两秒的情况下并行提速非常明显。思路是把一次迭代里的N只狼打包成N个 SimulationInput然后一次性并行仿真simIn(1:N) Simulink.SimulationInput(mdl); for i 1:N simIn(i) simIn(i).setVariable(Kp, X(i,1)); simIn(i) simIn(i).setVariable(Ki, X(i,2)); simIn(i) simIn(i).setVariable(Kd, X(i,3)); end simOut parsim(simIn, ShowProgress, off); for i 1:N t simOut(i).t_data; f simOut(i).f_data; fit(i) trapz(t, t .* abs(f)); end注意一点使用 parsim 时模型中引用 Kp、Ki、Kd 变量不能依赖基础工作区必须通过 setVariable 注入这就是前面推荐 SimulationInput 而不是 assignin 的根本原因。另外如果机器没有安装 Parallel Computing Toolboxparsim 无法使用可以退回普通 for 循环这是比较常见的限制。5. 结果分析收敛曲线、频率响应与性能对比5.1 迭代过程收敛性以我用的这组模型参数为例GWO优化过程中记录的收敛曲线整体趋势是前10代适应度下降非常快从初始的0.09附近快速跌到0.0215代之后开始缓慢下降到25代左右基本稳定在0.011左右后面几代只有微小波动。这说明GWO在这个问题上收敛速度足够快也基本没有出现明显的早熟停滞——精英保留机制把每一代的最优解都留住了收敛曲线是单调不增的。如果你的收敛曲线在30代还在明显下降说明迭代次数设少了可以加大 MaxIt 再跑。如果曲线很早进入平台并且最终适应度明显偏高多半是搜索范围太宽或者初始种群分布不够均匀先检查有没有大量仿真发散被罚成1e6。5.2 最优PID参数下的频率偏差响应下面是我在固定负荷阶跃 ΔPL0.01 p.u.、仿真50s条件下得到的一组示例结果。要注意的是不同模型参数、不同搜索范围得到的具体数值会不一样但趋势和结论基本一致。整定方法KpKiKdITAEΔfmin (×10⁻³ p.u.)调节时间ts (s)手调经验值0.400.200.100.0583-4.7622.5Ziegler-Nichols1.000.400.200.0241-3.1512.0GWO-PIDITAE0.860.550.170.0112-2.685.6GWO整定出的参数最明显的特点是积分增益比手调和ZN都大。积分项在LFC里承担消除稳态频率偏差的核心任务Ki偏大会让系统更快把Δf拉回零但纯积分作用容易引起超调所以GWO同时给了一个较小但足够有效的微分增益来压制超调。最终响应表现为频率最低点不那么深回零速度快而且几乎不振荡。手调参数其实也能用但对扰动非常敏感做鲁棒性测试时劣势很明显。5.3 不同负荷扰动下的鲁棒性验证做LFC只测一个扰动点没有说服力。我把负荷阶跃从0.01 p.u.改成0.02 p.u.再做同样的仿真GWO整定的PID参数完全不动重新仿出来的频率偏差峰值大约翻倍但调节时间变化不大稳态误差依然为0。对比手调参数在0.02扰动下的表现频率最低点更深回零后还有明显的一次过冲和振荡调节时间拉到30秒以上。这说明GWO优化的结果在扰动幅值变化时具有更好的鲁棒性这也是智能优化整定相比人工试凑最有说服力的地方。如果你想写更严谨的论文还可以做模型参数摄动测试把Tg和Tt分别上下调整20%甚至50%重新仿真观察GWO-PID的控制性能是否依然稳定。我在实际测试中Tg从0.08改到0.12、Tt从0.3改到0.45GWO参数仍然能让系统稳定回到零频率偏差只是响应速度变慢了一点。6. 实操复盘那些仿真和优化中的坑6.1 sim()循环调用的性能陷阱我最初把种群设为50、迭代设为60串行跑一次要等将近半天。后来发现两个性能瓶颈一是Simulink模型每次 sim() 调用时如果有未加载的引用模块或外部文件会反复编译拖慢速度二是只要模型中有一个 Scope 窗口开着每次仿真都会尝试刷新图形界面消耗大量时间。解决办法是优化前把所有Scope窗口关掉或者干脆把Scope输出用 To Workspace 替代如果必须观察波形等优化结束再用最优参数手动仿真看。另一个技巧是先用 load_system 把模型加载到内存避免每次 sim() 都重新加载整个系统。实测下来这两条能让总耗时减少三分之一以上。6.2 目标函数计算错误导致的假优化ITAE计算最隐蔽的坑不是公式错而是数组长度不匹配。用变步长ode45时如果 To Workspace 保存的是Arrayf_data和t_data的行数应该一致。但如果你在模型里用了一个采样时间不同的信号去驱动另一个To Workspace两个数组的行数就可能不一样trapz(t, t.*abs(f)) 直接报错进入catch返回1e6导致所有参数都判为发散GWO完全退化成随机搜索。我还踩过一个更隐蔽的坑用 simOut.f_data 取数时模型输出对象里的变量名大小写要和To Workspace设置完全一致否则会报“变量不存在”的错误。所以我在适应度函数里用了 try-catch目的不是掩盖错误而是让优化循环在个别异常时不至于整体中断但还是要靠日志定位问题。另一个容易被忽略的是仿真停止时间必须足够长否则ITAE的积分只包含截断部分会低估拖尾误差。优化出的参数可能在50s窗口内表现好拉长到80s就露馅。我建议先手动仿真观察系统回到稳态所需时间再据此确定停止时间。6.3 算法早熟与局部最优的排查思路GWO本身结构简单但也存在早熟可能。遇到的最典型情况是连续运行几次优化最优参数总是落在搜索范围的边界上比如 Kp 总是0.1下边界。这说明真实最优解可能位于范围之外或者适应度函数对边界附近的参数产生了误导性的惩罚。排查思路分三步第一步检查初始种群的分布如果非要用rand生成可能出现初始点全部挤在一侧的情况换成lhsdesign能明显改善第二步检查搜索范围是否过窄把范围扩大20%再跑一遍看最优参数是否依然贴边第三步检查适应度函数中发散惩罚值是否设置过大导致优化器宁可选择边界上稳定但性能差的解也不愿尝试可能发散但性能更好的解。对于LFC这类每次仿真都有数值误差的问题我建议每次优化用不同的随机种子多跑两三次取目标准则最优的那次作为最终结果。如果三种子的结果差距很大说明算法稳定性有问题优先考虑调整种群规模和搜索范围。6.4 让结果可复现、可写论文的小习惯最后聊几个写论文时会用到的小习惯。第一所有实验固定 rng 种子并在代码里把种子值写清楚审稿人问起来你随时能复现同一组参数。第二优化过程中保存每次迭代的完整数据包括全局最优适应度序列、最优参数和对应仿真响应别只存一个最终数值否则想画迭代曲线的时候只能重跑一遍。第三收敛曲线画的时候把纵轴用对数刻度如果下降趋势不明显对数坐标能更清晰地展示前期快速收敛和后期精细优化两个阶段。还有一个经验是在Simulink里把参数比较的几组仿真结果用不同的LineStyle和Legend标注好一次把对比图导出来。最优参数、手调参数、ZN参数、以及不同扰动下的响应曲线画在同一张图里论文里的结果分析部分基本就齐了。我自己的习惯是同时保存一个 .fig 文件和一个高分辨率 .png 输出前者方便后期改样式后者方便直接插入文章。整个流程跑通之后你会发现GWO整定PID的真正价值并不在于它一定比手工调出来的参数好多少倍而在于它把“调参”从一门玄学变成了一套可复现、可解释、可扩展的方法。换一套被控对象改一下适应度函数里的积分项加几组约束这套框架就能继续用。对于做LFC、AGC以及各类工业过程控制研究的朋友来说这套“GWOSimulinkITAE”的组合值得当成一个常用工具放在手边。
灰狼优化算法整定PID参数的单区域负荷频率控制Simulink实践
做单区域负荷频率控制LFC研究的朋友应该都有过这种体验在Simulink里搭好模型手工调PID参数调到头疼好不容易在这组负荷扰动下表现不错换个扰动步长又开始振荡想写论文又拿不出系统性的参数整定依据只能含糊地说“经过多次仿真调试确定”。我最初做这个课题时也是这个状态直到把灰狼优化算法GWO引入进来作为PID参数整定的搜索引擎整个流程才真正变得可复现、可解释、可写进论文。这篇文章记录的就是我完整跑通“GWO整定单区域负荷频率控制PID”的一次实践包含模型搭建、算法代码、Simulink联合仿真的每一步以及那些迭代了上百次才摸清的坑。1. 单区域LFC的模型长什么样PID整定难在哪1.1 频率偏差产生的物理逻辑先回到最基础的问题负荷频率控制到底在控什么。电力系统里发电出力要和负荷实时平衡一旦负荷增加而发电机出力没有跟上系统的频率就会下降反之负荷减少频率会升高。频率是全局性物理量直接反映发电与负荷的平衡状态偏离额定值比如50Hz会影响电机转速、变压器磁通和用户设备的正常工作。单区域负荷频率控制的“单区域”指的是忽略区域间联络线功率交换只研究一个独立控制区域。这个区域里有一个或多个发电机组共同承担频率调节任务。负荷突然增加 ΔPL 后系统频率先按惯性下降随后调速器通过一次调频动作调节进汽量/进水量来增加机械功率抬高频率。但仅有一次调频还不够调速器是按照调差系数R实现的本质上是有差调节频率不可能完全恢复到50Hz。要彻底消除稳态频差必须引入二次调频信号——也就是控制器给出的 ΔPc叠加到调速器的设定点上。这个二次控制器的角色在绝大多数研究里就是PID。1.2 我在Simulink里用的传递函数组合单区域LFC模型的简化结构非常经典由调速器、汽轮机、发电机-负荷环节和调差反馈组成。我用的传递函数和参数如下表环节传递函数参数调速器Gg(s) 1 / (Tg·s 1)Tg 0.08 s汽轮机非再热式Gt(s) 1 / (Tt·s 1)Tt 0.3 s发电机-负荷Gp(s) 1 / (M·s D)M 10 p.u.·sD 1 p.u.调差系数1/RR 0.05这里所有变量都用标幺值p.u.表示频率偏差 Δf 也是标幺值。有人习惯把发电机环节写成 Kp/(Tp·s 1)效果等价只是系数映射问题。我建议初学者直接按 M、D 来搭物理含义更清楚M 越大代表转子惯量越大频率下降越慢D 越大代表负荷阻尼越强频率越不容易偏离。这个模型在Simulink里搭起来就是几个传递函数模块的串联和反馈没有多少花哨的地方。但要注意一点调差反馈符号必须正确。调速器输入是 ΔPc 减去 Δf/R也就是一次调频是负反馈频率下降时调速器要开大阀门增加进汽。我第一次搭的时候把这条反馈的符号弄反了仿出来的频率直接发散排查了好久。1.3 传统整定方法为何不够用PID整定方法说起来一堆但真正做LFC时你会发现传统方法都有明显的别扭之处。试凑法就不提了完全看经验而且LFC模型里有惯性、延迟、反馈耦合三个参数互相影响在Simulink里反复改参数跑仿真非常浪费时间。Ziegler-Nichols频域整定法对这类高阶系统误差较大整出来的参数往往过激超调大而且它本质上针对的是带有延迟的一阶/二阶模型和LFC这种包含调速器、汽轮机、发电机多个惯性环节的系统匹配得并不好。更关键的一点是LFC关注的是动态性能与鲁棒性。PID参数不仅要让额定负荷扰动下的频率偏差响应快、超调小还要在负荷扰动幅值变化、甚至模型参数变化时依然可用。传统整定法很难把这些需求统一写进一个指标里而智能优化算法天然可以——只要把想优化的目标写成适应度函数剩下的搜索过程交给优化器。GWO之所以适合做这件事原因在于它原理简单、参数极少、不依赖梯度信息对于“黑箱式”的Simulink仿真模型非常友好。2. GWO算法原理拆解狼群行为与PID参数搜索的对应关系2.1 灰狼等级制度怎么变成寻优框架灰狼优化算法是Mirjalili等人在2014年提出的群体智能算法灵感来自灰狼群体的社会等级和捕食行为。灰狼群体分成四个等级α、β、δ、ω。α是头狼负责决策β是二把手辅助αδ负责侦察、放哨ω是最底层跟随前面三只狼移动。在GWO里每一个“狼”就是搜索空间中的一个候选解。对于我们的问题一个狼的位置就是一个三维向量 [Kp, Ki, Kd]代表一组PID参数。每只狼都要在Simulink里跑一次仿真得到一个ITAE适应度值。适应度值最小的三只狼被标记为α、β、δ它们相当于当前迭代中已知的最好、次好、第三好的PID参数组合其余狼ω根据这三只头狼的位置更新自己的位置逐步向最优区域收敛。这个过程很像团队协作不是只跟着最优秀的那个方案走而是同时参考前三名的方向避免被某个偶然的局部最优带偏。2.2 位置更新公式与收敛因子的作用GWO位置更新的数学原理不复杂。对任意一只狼 X它到猎物即最优解的距离用下式估计D |C · X_head - X(t)|然后用这个距离更新自己的位置X(t1) X_head - A · D其中 A 和 C 是两个关键系数A 2·a·r1 - a C 2·r2r1 和 r2 是 [0,1] 之间的随机数a 是从2线性递减到0的收敛因子。A 的绝对值决定了狼是扩大搜索范围|A| 1探索还是收缩逼近|A| 1开发这个机制和很多优化算法里的“前期全局搜索、后期局部精搜”策略异曲同工。C 的作用是给猎物位置增加一个随机权重相当于在搜索中引入扰动防止算法过早陷入局部最优。三只头狼分别按上面的公式产生三个候选位置 X1、X2、X3最终位置取三者的平均值X(t1) (X1 X2 X3) / 3这个“三狼平均”的策略是GWO的精华。它可以被理解为最优的PID参数大概率位于α、β、δ三个已知较优参数连线的几何中心附近但又不完全等于它们留出了继续探索的空间。2.3 适应度函数选ITAE而不选ISE的理由适应度函数是整个优化过程的“指挥棒”。你用它评价哪组PID参数好GWO就会朝着让它最小化的方向搜索所以目标函数的设计直接决定最终参数的行为。LFC里最常用的几个误差积分指标指标公式特点IAE∫ |Δf| dt对误差绝对值积分简单但不对时间加权ISE∫ Δf² dt对大误差惩罚很大但对长期小误差不敏感ITAE∫ t·|Δf| dt对时间加权后期误差影响越来越大我最终选了ITAE。原因很实际在负荷频率控制里我们不仅希望频率偏差小更希望它尽快回到零。ITAE会给“拖尾”的误差更大的惩罚所以优化出来的PID参数通常响应更快、调节时间更短不会出现为了压低超调而磨磨蹭蹭收敛的情况。ISE虽然抑制大偏差效果好但容易让优化器忽略那些低频小振幅的振荡尾巴。如果你在写论文时想让控制器兼顾控制能量可以把指标改成 J w1·ITAE(Δf) w2·∫ΔPc²dt用权重系数w1、w2权衡动态响应和控制力度。我建议先不加惩罚项把基本流程跑通再扩展。2.4 搜索空间与初始种群设置GWO的搜索空间就是PID参数的取值范围。这看起来很自由但很影响最终效果。范围设太大大量参数组合会让Simulink仿真发散适应度函数只能返回一个惩罚大数浪费计算资源范围设太小可能把最优参数卡在边界外面。我的做法分两步先用一组粗略范围比如 Kp∈[0,5]Ki∈[0,2]Kd∈[0,1]手动跑几次仿真观察Kp超过多少仿真开始发散、Ki多大时系统开始振荡然后把优化范围收缩到一个“稳定且合理”的区间。以我上面的模型为例最终用的范围是Kp ∈ [0.1, 2.0]Ki ∈ [0.01, 1.0]Kd ∈ [0, 0.5]初始种群建议用拉丁超立方采样lhsdesign替代纯随机它能让初始点更均匀地铺满搜索空间减少初始分布不均导致的早熟。如果电脑上没装统计工具箱用 rand 生成也完全可以。3. Simulink模型搭建从模块到可调参数3.1 模型拓扑与关键模块参数模型文件我命名为 LFC_single_area.slx不建议用中文名也不建议放在带中文路径的目录下否则 sim() 调用容易出幺蛾子。模型内部的信号连接顺序如下。先看被控对象部分三个传递函数模块串联。调速器前有一个求和点 Sum2计算 ΔPc - Δf/R调速器输出进入汽轮机汽轮机输出是机械功率变化量 ΔPm。接着用 Sum3 计算 ΔPm - ΔPL输入到 1/(M·sD) 环节输出就是频率偏差 Δf。Δf 引出三条支路去 Scope 和 To Workspace用于观察和记录数据经过 Gain(-1) 变成误差信号 e -Δf送入 PID Controller经过 Gain(1/R) 接入 Sum2 的负号端口实现一次调频反馈。这里分成两条负反馈路径是LFC模型的关键一条是一次调频的调差反馈一条是PID二次调频。前者是系统固有机制后者才是我们要优化的控制器。PID Controller 模块我建议用并联形式Parallel form比例增益设置为工作区变量 Kp积分增益 Ki微分增益 Kd。微分项一定要用带滤波的形式默认的D filter coefficient N100 就够否则理想微分会对仿真噪声非常敏感。3.2 参数动态注入方式GWO迭代过程中每一只狼都是一组新的 Kp、Ki、Kd怎么把这些值传进Simulink模型是关键问题。很多人第一反应是用 assignin(base, Kp, Kp) 往基础工作区塞变量模型里直接引用变量名即可。这种方式最简单但有两个隐患一是并行计算时多个仿真相互干扰二是如果仿真中途报错基础工作区可能残留脏数据。更推荐的做法是用 Simulink.SimulationInput 的 setVariable 方法把参数注入到这一次仿真的临时上下文中。这样既不污染基础工作区也方便后续用 parsim 做并行优化。代码结构后面专门讲。要注意的是模型中只要出现变量名 Kp、Ki、Kd那么在手动仿真之前工作区里必须存在这些变量否则模块参数解析会报错。最简单的方法是先手动命令行赋值 Kp0.5; Ki0.3; Kd0.1; 再打开模型。3.3 仿真配置与数据采集接口为了配合GWO循环里反复调用 sim()仿真参数按下面设置仿真停止时间50s。LFC系统的动态过程一般在10-20秒内基本结束50秒足够让收敛尾巴也计入ITAE。求解器变步长 ode45相对误差默认 1e-3 就可以。优化过程不需要超级精度而且每次仿真时间太长会让整个优化跑得非常慢。数据采集从Δf信号引一个 To Workspace 模块变量名 f_data保存格式选择 Array用 Clock 模块接另一个 To Workspace变量名 t_data同样用 Array。用 Array 格式而不是 Timeseries是为了在适应度函数里直接得到普通数组省去拆数据对象的麻烦。如果你选择 Timeseries后面取数要多写两行代码不是不行只是没必要。搭建完成后先不要跑优化手动输入一组大概率稳定的参数仿真一次确认频率响应形状正常——出现阶跃后频率下降、然后回升到0整个波形没有发散、没有剧烈振荡再进入下一步。4. GWO整定PID的Matlab实现主程序与适应度函数4.1 主程序框架与初始化整个优化工程包含两个文件GWO_LFC_main.m 是主程序负责GWO的初始化与迭代LFC_Fitness.m 是适应度函数负责调用Simulink模型并计算ITAE。主程序开头先清空环境并固定随机种子保证结果可复现。这一点做研究的人一定要养成习惯否则每次跑出的最优参数都不一样论文里根本没法写。GWO_LFC_main.m 核心框架clc; clear; close all; rng(0); mdl LFC_single_area; load_system(mdl); N 30; % 种群规模 MaxIt 40; % 最大迭代次数 dim 3; lb [0.1, 0.01, 0]; % Kp, Ki, Kd 下界 ub [2.0, 1.0, 0.5]; % 上界 X lhsdesign(N, dim) .* (ub - lb) lb; fit zeros(N, 1); Alpha_pos zeros(1, dim); Alpha_fit inf; Beta_pos zeros(1, dim); Beta_fit inf; Delta_pos zeros(1, dim); Delta_fit inf; Conv zeros(MaxIt, 1); gbest zeros(1, dim); gbest_fit inf;Alpha、Beta、Delta分别存储三只头狼的位置和适应度。gbest是全局历史最优解配合后面的精英保留策略使用。4.2 GWO迭代更新核心代码主循环分成三块计算每个个体的适应度并更新三只头狼执行精英保留更新所有狼的位置。for t 1:MaxIt for i 1:N X(i,:) max(X(i,:), lb); X(i,:) min(X(i,:), ub); fit(i) LFC_Fitness(X(i,:), mdl); if fit(i) Alpha_fit Alpha_fit fit(i); Alpha_pos X(i,:); elseif fit(i) Beta_fit Beta_fit fit(i); Beta_pos X(i,:); elseif fit(i) Delta_fit Delta_fit fit(i); Delta_pos X(i,:); end end if Alpha_fit gbest_fit gbest_fit Alpha_fit; gbest Alpha_pos; end X(1,:) gbest; % 精英保留 Conv(t) gbest_fit; a 2 - 2 * t / MaxIt; for i 1:N if i 1, continue; end for j 1:dim r1 rand; r2 rand; A 2*a*r1 - a; C 2*r2; D_alpha abs(C * Alpha_pos(j) - X(i,j)); X1 Alpha_pos(j) - A * D_alpha; r1 rand; r2 rand; A 2*a*r1 - a; C 2*r2; D_beta abs(C * Beta_pos(j) - X(i,j)); X2 Beta_pos(j) - A * D_beta; r1 rand; r2 rand; A 2*a*r1 - a; C 2*r2; D_delta abs(C * Delta_pos(j) - X(i,j)); X3 Delta_pos(j) - A * D_delta; X(i,j) (X1 X2 X3) / 3; end end end fprintf(最优PID参数: Kp%.4f, Ki%.4f, Kd%.4f\n, gbest); fprintf(最优ITAE: %.6f\n, gbest_fit); save(GWO_results.mat, gbest, gbest_fit, Conv);我加了一个精英保留直接把历史最优个体放到第1个位置并且不让它参与后续位置更新。为什么这么做标准GWO用三只狼的平均位置生成新解连α狼自己也可能会被“平均掉”导致当前最优解在下一轮就丢了。对LFC这种适应度函数有随机性和噪声的优化问题精英保留几乎是必须的。4.3 适应度函数与sim()调用细节LFC_Fitness.m 是所有仿真发生的地方。这里有几个必须注意的细节每次仿真前要把三个PID参数注入到当前仿真上下文仿真结束后要稳定地取出时间和频率数据处理发散和无解情况。function J LFC_Fitness(x, mdl) simIn Simulink.SimulationInput(mdl); simIn simIn.setVariable(Kp, x(1)); simIn simIn.setVariable(Ki, x(2)); simIn simIn.setVariable(Kd, x(3)); try simOut sim(simIn); t simOut.t_data; f simOut.f_data; f f(:); t t(:); if any(isnan(f)) || any(isinf(f)) || max(abs(f)) 1e3 J 1e6; else J trapz(t, t .* abs(f)); end catch J 1e6; end endITAE用 trapz(t, t.*abs(f)) 做数值积分比手写 for 循环精确也比简单的 sum(f)*dt 更适应变步长求解器。发散惩罚值设为1e6是一个常用技巧。仿真发散时数据里会出现NaN或Inf此时返回一个很大的适应度GWO会自动避开这批参数。要注意的是如果惩罚值不够大比如只有1那么一些轻微发散但恰好ITAE比较小的解反而会被当作最优导致最终参数不可用。还有一点如果模型路径或模型名有中文或者Simulink模型没有加载sim(simIn)会报错catch块会兜住并返回惩罚值。但这会掩盖问题让你误以为是参数不好。建议优化前先手动 load_system 一次确认模型能打开、能手动仿真再跑自动优化。4.4 用parsim把优化速度提升数倍如果你只跑40代×30只狼1200次仿真串行执行可能要好几个小时。我强烈建议把适应度计算改成 parsim 并行版本尤其是在模型不大、单次仿真只要一两秒的情况下并行提速非常明显。思路是把一次迭代里的N只狼打包成N个 SimulationInput然后一次性并行仿真simIn(1:N) Simulink.SimulationInput(mdl); for i 1:N simIn(i) simIn(i).setVariable(Kp, X(i,1)); simIn(i) simIn(i).setVariable(Ki, X(i,2)); simIn(i) simIn(i).setVariable(Kd, X(i,3)); end simOut parsim(simIn, ShowProgress, off); for i 1:N t simOut(i).t_data; f simOut(i).f_data; fit(i) trapz(t, t .* abs(f)); end注意一点使用 parsim 时模型中引用 Kp、Ki、Kd 变量不能依赖基础工作区必须通过 setVariable 注入这就是前面推荐 SimulationInput 而不是 assignin 的根本原因。另外如果机器没有安装 Parallel Computing Toolboxparsim 无法使用可以退回普通 for 循环这是比较常见的限制。5. 结果分析收敛曲线、频率响应与性能对比5.1 迭代过程收敛性以我用的这组模型参数为例GWO优化过程中记录的收敛曲线整体趋势是前10代适应度下降非常快从初始的0.09附近快速跌到0.0215代之后开始缓慢下降到25代左右基本稳定在0.011左右后面几代只有微小波动。这说明GWO在这个问题上收敛速度足够快也基本没有出现明显的早熟停滞——精英保留机制把每一代的最优解都留住了收敛曲线是单调不增的。如果你的收敛曲线在30代还在明显下降说明迭代次数设少了可以加大 MaxIt 再跑。如果曲线很早进入平台并且最终适应度明显偏高多半是搜索范围太宽或者初始种群分布不够均匀先检查有没有大量仿真发散被罚成1e6。5.2 最优PID参数下的频率偏差响应下面是我在固定负荷阶跃 ΔPL0.01 p.u.、仿真50s条件下得到的一组示例结果。要注意的是不同模型参数、不同搜索范围得到的具体数值会不一样但趋势和结论基本一致。整定方法KpKiKdITAEΔfmin (×10⁻³ p.u.)调节时间ts (s)手调经验值0.400.200.100.0583-4.7622.5Ziegler-Nichols1.000.400.200.0241-3.1512.0GWO-PIDITAE0.860.550.170.0112-2.685.6GWO整定出的参数最明显的特点是积分增益比手调和ZN都大。积分项在LFC里承担消除稳态频率偏差的核心任务Ki偏大会让系统更快把Δf拉回零但纯积分作用容易引起超调所以GWO同时给了一个较小但足够有效的微分增益来压制超调。最终响应表现为频率最低点不那么深回零速度快而且几乎不振荡。手调参数其实也能用但对扰动非常敏感做鲁棒性测试时劣势很明显。5.3 不同负荷扰动下的鲁棒性验证做LFC只测一个扰动点没有说服力。我把负荷阶跃从0.01 p.u.改成0.02 p.u.再做同样的仿真GWO整定的PID参数完全不动重新仿出来的频率偏差峰值大约翻倍但调节时间变化不大稳态误差依然为0。对比手调参数在0.02扰动下的表现频率最低点更深回零后还有明显的一次过冲和振荡调节时间拉到30秒以上。这说明GWO优化的结果在扰动幅值变化时具有更好的鲁棒性这也是智能优化整定相比人工试凑最有说服力的地方。如果你想写更严谨的论文还可以做模型参数摄动测试把Tg和Tt分别上下调整20%甚至50%重新仿真观察GWO-PID的控制性能是否依然稳定。我在实际测试中Tg从0.08改到0.12、Tt从0.3改到0.45GWO参数仍然能让系统稳定回到零频率偏差只是响应速度变慢了一点。6. 实操复盘那些仿真和优化中的坑6.1 sim()循环调用的性能陷阱我最初把种群设为50、迭代设为60串行跑一次要等将近半天。后来发现两个性能瓶颈一是Simulink模型每次 sim() 调用时如果有未加载的引用模块或外部文件会反复编译拖慢速度二是只要模型中有一个 Scope 窗口开着每次仿真都会尝试刷新图形界面消耗大量时间。解决办法是优化前把所有Scope窗口关掉或者干脆把Scope输出用 To Workspace 替代如果必须观察波形等优化结束再用最优参数手动仿真看。另一个技巧是先用 load_system 把模型加载到内存避免每次 sim() 都重新加载整个系统。实测下来这两条能让总耗时减少三分之一以上。6.2 目标函数计算错误导致的假优化ITAE计算最隐蔽的坑不是公式错而是数组长度不匹配。用变步长ode45时如果 To Workspace 保存的是Arrayf_data和t_data的行数应该一致。但如果你在模型里用了一个采样时间不同的信号去驱动另一个To Workspace两个数组的行数就可能不一样trapz(t, t.*abs(f)) 直接报错进入catch返回1e6导致所有参数都判为发散GWO完全退化成随机搜索。我还踩过一个更隐蔽的坑用 simOut.f_data 取数时模型输出对象里的变量名大小写要和To Workspace设置完全一致否则会报“变量不存在”的错误。所以我在适应度函数里用了 try-catch目的不是掩盖错误而是让优化循环在个别异常时不至于整体中断但还是要靠日志定位问题。另一个容易被忽略的是仿真停止时间必须足够长否则ITAE的积分只包含截断部分会低估拖尾误差。优化出的参数可能在50s窗口内表现好拉长到80s就露馅。我建议先手动仿真观察系统回到稳态所需时间再据此确定停止时间。6.3 算法早熟与局部最优的排查思路GWO本身结构简单但也存在早熟可能。遇到的最典型情况是连续运行几次优化最优参数总是落在搜索范围的边界上比如 Kp 总是0.1下边界。这说明真实最优解可能位于范围之外或者适应度函数对边界附近的参数产生了误导性的惩罚。排查思路分三步第一步检查初始种群的分布如果非要用rand生成可能出现初始点全部挤在一侧的情况换成lhsdesign能明显改善第二步检查搜索范围是否过窄把范围扩大20%再跑一遍看最优参数是否依然贴边第三步检查适应度函数中发散惩罚值是否设置过大导致优化器宁可选择边界上稳定但性能差的解也不愿尝试可能发散但性能更好的解。对于LFC这类每次仿真都有数值误差的问题我建议每次优化用不同的随机种子多跑两三次取目标准则最优的那次作为最终结果。如果三种子的结果差距很大说明算法稳定性有问题优先考虑调整种群规模和搜索范围。6.4 让结果可复现、可写论文的小习惯最后聊几个写论文时会用到的小习惯。第一所有实验固定 rng 种子并在代码里把种子值写清楚审稿人问起来你随时能复现同一组参数。第二优化过程中保存每次迭代的完整数据包括全局最优适应度序列、最优参数和对应仿真响应别只存一个最终数值否则想画迭代曲线的时候只能重跑一遍。第三收敛曲线画的时候把纵轴用对数刻度如果下降趋势不明显对数坐标能更清晰地展示前期快速收敛和后期精细优化两个阶段。还有一个经验是在Simulink里把参数比较的几组仿真结果用不同的LineStyle和Legend标注好一次把对比图导出来。最优参数、手调参数、ZN参数、以及不同扰动下的响应曲线画在同一张图里论文里的结果分析部分基本就齐了。我自己的习惯是同时保存一个 .fig 文件和一个高分辨率 .png 输出前者方便后期改样式后者方便直接插入文章。整个流程跑通之后你会发现GWO整定PID的真正价值并不在于它一定比手工调出来的参数好多少倍而在于它把“调参”从一门玄学变成了一套可复现、可解释、可扩展的方法。换一套被控对象改一下适应度函数里的积分项加几组约束这套框架就能继续用。对于做LFC、AGC以及各类工业过程控制研究的朋友来说这套“GWOSimulinkITAE”的组合值得当成一个常用工具放在手边。
深度解析：为什么90帧还卡顿？)
