前阵子搭了一个优化演示工程帮一个朋友把 IEEE 30 节点系统的网损降下来。算法选的粒子群底层潮流计算交给 Matpower 引擎整套跑在 Matlab 里。借这个机会我把整个项目的建模思路、代码骨架、以及调试过程中踩过的坑完整记录下来。严格来说IEEE 30 节点系统是一个输电网测试算例但在实际工作中很多做配电网算法验证的同学也拿它当作标准网架来用所以这个标题组合非常经典。文章不空谈原理重点讲清楚我怎么建模、怎么把粒子群和 Matpower 串起来、以及最终如何分析和判断结果希望能给正在做电力系统智能优化方向的同学一个直接可参考的模板。1. 这个优化问题到底在算什么有功-无功联合优化的数学模型与工程动因1.1 为什么要把有功和无功放在同一个问题里很多刚开始接触电网优化的同学会有一个困惑有功功率和无功功率不是两回事吗有功管系统频率无功管电压水平为什么非要放在一起优化其实这是从运行机理上就耦合在一起的。电网里任何一条线路、一台变压器它的损耗都可以写成 (P_{loss} I^2 R) 的形式而电流 I 的幅值同时由有功潮流和无功潮流共同决定。也就是说线路上的无功流动同样会带来发热损耗。如果你只调整有功出力不管无功分布那么某些线路上的无功环流会一直存在网损降不下去反过来只调无功有功的不合理分配又会让重载线路继续超负荷。所以真正要降网损必须同时调整有功出力和无功分布这就是有功-无功联合优化的由来。配电网的场景里这个问题更明显。配电网的线路 R/X 比值大无功流动造成的电压跌落和线损占比比输电网高得多。分布式电源接入之后有功和无功的耦合关系又进一步复杂化。所以在配电网规划、降损、电压治理这些课题里有功-无功联合优化基本上是一个绕不开的基础问题。1.2 目标函数与决策变量怎么选这次工程里目标函数取的是系统总网损最小。原因很朴素网损是可观测、可计算的直接经济指标也是最好和基态潮流做对比的结果量。你要是想换成发电成本最小、电压偏差最小方法完全一样只是改动适应度函数里的一个表达式而已。决策变量的选取是建模里最关键的一步。IEEE 30 节点系统一共有 30 个节点、41 条支路6 台发电机。我采用的是最常见的变量组合方式除平衡节点外的 5 台发电机有功出力6 台发电机的机端电压幅值。也就是说一个粒子是一个 11 维向量里面的每个分量被提交给 Matpower 潮流计算最终得到一组潮流结果和一个网损值。为什么用机端电压代替无功出力因为这两个量在潮流计算里是强耦合的你设定了机端电压潮流解出来后自然会得到对应的无功注入。直接拿无功出力当决策变量反而麻烦因为 Matpower 在计算时会有 PV 节点和 PQ 节点之间的自动切换逻辑直接用无功出力很容易撞上这个切换边界。这一点后面踩坑部分我会单独说。1.3 约束条件从哪来潮流等式与设备不等式优化问题不能只有目标函数约束条件才是让它像电网而不是像数学游戏的关键。约束分两类等式约束就是潮流方程本身配电网或者输电网的任何运行状态都必须满足节点功率平衡。这个我不会自己手写而是委托给 Matpower 的runpf函数来完成跑通了就代表等式约束被满足了。不等式约束包括发电机有功上下限、发电机无功上下限、节点电压幅值上下限、线路传输容量上限。这些限制在 Matpower 的case30数据文件里都有初始值。bus表里有每个节点的正常电压上限Vmax和下限Vmingen表里有每台发电机的Pmax、Pmin、Qmax、Qmin。在粒子群算法框架下不等式约束大部分靠罚函数处理小部分靠边界截断处理这个设计展开讲就是第 4 章的内容。让我点一下数学上的关键逻辑目标函数要优化的是网损决策变量是发电机有功和机端电压约束是潮流方程加设备上下限。把这三句话写清楚这个优化问题才算是被正确建模了。很多同学代码写不明白往往不是粒子群没掌握而是决策变量和约束没有从 Matpower 数据结构上想透。2. 为什么用粒子群而不是直接调用内点法2.1 Matpower自带最优潮流的局限Matpower 自带的runopf函数是用内点法求解最优潮流OPF的精度高、速度快在干净的标准算例上效果相当好。那为什么还要自己做粒子群一个很现实的原因是实际工程项目里的优化模型往往和教科书上的标准 OPF 不完全一样。你可能要加一个不光滑的约束比如某个节点的电压合格率分段惩罚可能要把目标函数改成某种启发式指标可能要对某些设备动作次数做惩罚。这些非标准、非光滑、甚至不可导的改动对传统内点法是致命的因为内点法需要构造拉格朗日函数并求梯度或者海森矩阵一旦模型变得不光滑数值上非常容易发散。粒子群这类群智能算法则是另一套逻辑它不要求目标函数可导不要求约束光滑它只把目标函数当成一个黑盒来调用。潮流算得出来就给一个具体的网损值算不出来就给一个巨大的惩罚值。这种黑盒优化的特性让它在科研探索阶段显得特别灵活。2.2 粒子群适合这种问题的真正原因工程上选择粒子群主要有三个原因。第一实现成本低。一个有速度-位置更新公式加适应度计算的粒子群核心代码写得好也就一两百行。相比于内点法复杂的矩阵求导和修正方程粒子群的代码量要小得多调起来也更直观。第二天然支持并行。粒子群每一代要计算多个粒子的适应度而这些计算彼此独立完全可以用parfor并行跑。在我这个例子里每个粒子都要调用一次runpf用并行池可以把整体耗时压缩到原来的三分之一左右。第三全局搜索能力比传统梯度类算法好。电网优化问题本质上是有多个局部最优点的非线性非凸问题内点法非常依赖初值初值不好就容易收敛到局部解或者压根不收敛。粒子群通过种群协作的方式在解空间里撒点搜索虽然不保证找到全局最优但工程实践里找到的解往往是可用的而且稳定性好。2.3 它与遗传算法的简单对比也有人问我为什么不直接上遗传算法两个算法我都实现过结论是在电网潮流作为黑盒这个场景下粒子群通常更省事。遗传算法的选择、交叉、变异三个算子每个都需要调试合理的概率参数否则要么收敛太慢要么种群过早失去多样性。粒子群的参数相对集中核心就是惯性权重和学习因子尤其惯性权重可以按迭代次数线性递减整体调参路径清晰得多。另外遗传算法的二进制编码在边界连续优化问题里会有精度和编码长度之间的矛盾而粒子群直接用实数向量操作不需要解码天然贴合发电机出力用连续实数表示的需求。当然如果问题本身是离散变量主导的比如电容器分组投切、变压器分接头挡位遗传算法的离散处理反而更方便。连续变量主导的有功-无功优化粒子群是性价比最高的选择。3. 环境准备与数据摸底Matpower、case30和配电网改造3.1 十分钟搞定Matpower环境Matpower 是一个纯 Matlab 实现的开源电力系统分析工具包不需要编译下载解压就能用。从官网或者 GitHub 下载最新稳定版本之后把解压目录整个加入 Matlab 路径即可。%% 将Matpower路径加入Matlab环境 addpath(genpath(D:/matpower7.1)); savepath; % 保存路径避免每次启动重新设置 mpver(ver); % 查看并确认Matpower版本运行mpver(ver)能看到默认的case30等标准算例已经内置。我建议每次用 Matpower 之前都显式确认版本号因为不同版本对runpf返回结果的字段略有差异确认版本可以少很多莫名其妙的数组越界报错。3.2 读懂case30数据结构bus、gen、branchMatpower 的核心数据结构是mpc结构体它包含bus、gen、branch三个大矩阵以及baseMVA这个基准容量。case30.m文件本质上就是这些矩阵的文本初始化。bus矩阵的关键列第 1 列节点编号第 2 列节点类型1 是 PQ 节点2 是 PV 节点3 是平衡节点第 3、4 列是有功负荷和无功负荷第 8 列是节点电压幅值初始值第 12、13 列是正常运行时电压上限和下限。gen矩阵的关键列第 1 列是发电机所在节点号第 2 列是有功出力第 3 列是无功出力第 4、5 列是无功上限和下限第 6 列是机端电压幅值设定值第 9、10 列是有功上限和下限。branch矩阵记录线路或者变压器的参数包括首末端节点、电阻、电抗、电纳、长期允许载流量、变压器变比等。第一件事永远是用loadcase(case30)载入数据然后用runpf跑一次基态潮流看看初始网损是多少。基态潮流的结果是整个优化项目的比较基准后面所有优化效果的评估都靠它说话。mpc loadcase(case30); % 载入IEEE30节点数据 res0 runpf(mpc); % 基态潮流 loss0 sum(res0.branch(:, 14) res0.branch(:, 16)); % 基态网损单位MW fprintf(基态网损: %.4f MW\n, loss0);res.branch矩阵的第 14 列是首端有功注入第 16 列是末端有功注入两者相加就是该支路的损耗。Matpower 返回的success字段等于 1 才代表潮流收敛这个字段在后续迭代里非常重要。3.3 再说一遍这是输电网算例怎么改成配电网风格我必须把话说清楚IEEE 30 节点系统本身是输电网络线路阻抗、电压等级、负荷水平都不是典型配电网。但是在算法开发和验证阶段它的拓扑结构、发电机数量、约束类型足够丰富完全可以通过改造来模拟配电网特征。常见的改造思路有三种把部分发电机改成分布式电源模型降低其出力上限并修改无功调节范围将线路参数调大电阻分量提高 R/X 比值让电压问题更突出增加部分节点的无功负荷让无功分布更加紧张。这样改完之后还是在case30的框架下做优化但结果更接近配电网降损的真实需求。如果你已经有某个具体的配电网拓扑数据也可以用自己的潮流数据替换掉case30粒子群的代码框架完全不用动只需要改loadcase那一个地方。4. 粒子群与潮流计算耦合编码、适应度函数、主循环4.1 粒子编码设计一串向量对应一台发电机的运行方式粒子群优化的第一个核心步骤是编码。我的设计是把每台非平衡发电机的有功出力和全部发电机的机端电压幅值拼接成一维向量。以case30为例6 台发电机节点编号分别是 1、2、5、8、11、13其中节点 1 是平衡节点。那么决策向量定义为前 5 个分量节点 2、5、8、11、13 这 5 台发电机的有功出力后 6 个分量节点 1、2、5、8、11、13 这 6 台发电机的机端电压幅值。总维度是 11。这里有一个细节有些代码会把平衡节点的有功也拿来做变量这是不对的。平衡节点的有功在潮流里是补差角色系统总有功不平衡时由它自动平衡把它当决策变量会和潮流方程矛盾。为了让粒子初始化范围合理我直接从mpc.gen里读取Pmax、Pmin、Vmax、Vmin来生成上下界。这样粒子群搜索空间和物理设备允许范围自动对齐省去手动维护边界数据的麻烦。% 读取设备上下限 genBus mpc.gen(:, 1); isSlack zeros(size(genBus)); % 平衡节点由bus矩阵type列为3的节点确定 slackIdx find(mpc.bus(:, 2) 3); isSlack ismember(genBus, slackIdx); % 有功决策量范围: 非平衡发电机 PgMin mpc.gen(~isSlack, 10); PgMax mpc.gen(~isSlack, 9); % 电压决策量范围: 全部发电机 VgMin 0.94; VgMax 1.06;4.2 适应度函数网损加罚函数一个值得反复调试的核心适应度函数是粒子群和电网模型之间的唯一接口。我的设计分三层。第一层把粒子的 11 个分量写回mpc.gen调用runpf计算潮流。如果潮流不收敛直接返回一个很大的惩罚值比如1e6。这一步一定要包在try-catch里因为粒子群探索阶段经常会产生不合理的控制变量组合让潮流计算直接报错退出如果不捕获异常整个优化循环就崩了。第二层计算网损。在潮流收敛的情况下用首末端有功之和得到系统网损。第三层罚函数处理越限。我主要检查两个地方节点电压幅值是否越界以及发电机无功出力是否越限。节点电压在res.bus第 8 列发电机无功在res.gen第 3 列。越限越严重惩罚值越大这样粒子群在迭代过程中会自然避开不合理区域。function loss psoFitness(x, mpcBase) mpc mpcBase; genBus mpc.gen(:, 1); slackIdx find(mpc.bus(:, 2) 3); isSlack ismember(genBus, slackIdx); % 前5维写到非平衡机有功 mpc.gen(~isSlack, 2) x(1:5); % 后6维写到全部机端电压 mpc.gen(:, 6) x(6:11); try res runpf(mpc); if res.success ~ 1 loss 1e6; return; end % 网损 loss sum(res.branch(:, 14) res.branch(:, 16)); % 电压越限惩罚 vmin min(res.bus(:, 8)); vmax max(res.bus(:, 8)); if vmin 0.95 loss loss 2000 * abs(vmin - 0.95); end if vmax 1.05 loss loss 2000 * abs(vmax - 1.05); end % 发电机无功越限惩罚 if min(res.gen(:, 3) - res.gen(:, 4)) 0 || min(res.gen(:, 5) - res.gen(:, 3)) 0 loss loss 3000; end catch loss 1e6; end end这个适应度函数里罚系数 2000、3000、1e6 都是我反复试验后比较稳的组合。罚函数设计是这个项目里最需要耐心调试的部分权重太小约束形同虚设权重太大粒子群会过早向可行域收缩丧失搜索能力。后面我会单独讲怎么判断罚函数是否合适。4.3 速度位移更新与边界处理粒子群的标准更新公式不用担心直接用经典版本v w * v c1 * r1 * (pbest - x) c2 * r2 * (gbest - x) x x v其中w是惯性权重c1是自我认知学习因子c2是社会学习因子r1、r2是 0 到 1 的随机数。我的参数设定如下表参数取值说明惯性权重 w0.9 到 0.4 线性递减前期全局搜索后期局部精细搜索学习因子 c12.0保留自我经验学习因子 c22.0群体信息共享种群规模40网络规模不大40 个粒子足够最大迭代次数200配合收敛曲线判断是否提前结束速度上限 vmax0.3每个决策量按范围归一化后限制步长线性递减的惯性权重是工程实践里非常稳定的配置。迭代初期w大粒子跑得快可以探索大范围迭代后期w小粒子在局部精细搜索。很多失败的粒子群应用案例问题都出在全程用同一个w导致搜索行为要么太粗放要么太保守。边界处理上我的做法是粒子越界时直接把对应分量拉回到边界值并且把该维速度设为 0。这个重置边界策略比随机反弹策略更容易让粒子在可行域边界附近稳定下来因为在电力系统优化里很多最优解往往就在发电机出力上限附近。4.4 主循环骨架与计算量预估主循环的骨架大致是这样% 初始化粒子群 nPop 40; nVar 11; maxIter 200; x repmat(xMin, nPop, 1) rand(nPop, nVar) .* repmat(xMax - xMin, nPop, 1); v zeros(nPop, nVar); pbest x; pbestVal arrayfun((i) psoFitness(x(i, :), mpc), 1:nPop); gbestIdx find(pbestVal min(pbestVal), 1); gbest pbest(gbestIdx, :); gbestVal pbestVal(gbestIdx); trace zeros(maxIter, 1); % 主迭代 for iter 1:maxIter w 0.9 - 0.5 * (iter / maxIter); % 线性递减 for i 1:nPop r1 rand(nVar, 1); r2 rand(nVar, 1); v(i, :) w * v(i, :) 2.0 * r1 .* (pbest(i, :) - x(i, :)) 2.0 * r2 .* (gbest - x(i, :)); v(i, v(i, :) vmax) vmax; v(i, v(i, :) -vmax) -vmax; x(i, :) x(i, :) v(i, :); x(i, x(i, :) xMax) xMax; x(i, x(i, :) xMin) xMin; fval psoFitness(x(i, :), mpc); if fval pbestVal(i) pbest(i, :) x(i, :); pbestVal(i) fval; end end [minVal, idx] min(pbestVal); if minVal gbestVal gbestVal minVal; gbest pbest(idx, :); end trace(iter) gbestVal; end计算量方面40 个粒子乘以 200 代就是 8000 次潮流计算。在普通笔记本上一次case30潮流大约几十毫秒所以整个优化过程大概一到三分钟。如果需要跑大规模配电网建议改用parfor对粒子循环做并行化实测能缩短到原来的四分之一左右。5. 跑通之后怎么分析结果迭代曲线、电压分布、网损对比5.1 从基态潮流开始先记住优化前的网损不管优化结果多漂亮没有基准值就毫无意义。我的习惯是把基态潮流结果、基态网损、基态最低最高电压全部打印出来作为后续对比的 Baseline。基态网损由runpf直接得到然后优化结束之后把gbest重新写入mpc.gen再跑一次潮流同样统计网损和电压。两者做差除基态就是降损率。这里有一个细节优化停止时gbest对应的潮流并不一定已经用最终粒子再算过一遍因为pbestVal可能是粒子在某一代的位置计算出来的。所以最终结果必须先重算潮流再统计指标否则容易出现结果和实际不一致的尴尬。5.2 迭代曲线与收敛判断粒子群迭代过程里我把每一代的最优适应度记录在trace数组里用semilogy画收敛曲线。注意这里我通常用对数坐标因为早期适应度可能从几千快速下降到十几个 MW线性坐标会让后期曲线贴在横轴上看不出细节。一个正常的收敛曲线应该是前期快速下降中期平缓波动后期基本稳定。如果曲线后期还在明显持续下降说明迭代次数不够要加大maxIter。如果曲线很快就完全不动了说明种群可能早熟所有粒子已经挤到了同一个局部最优点。figure; plot(1:maxIter, trace, LineWidth, 1.5); xlabel(迭代次数); ylabel(最优网损/MW); title(粒子群收敛曲线); grid on;5.3 用runopf兜底验证粒子群结果是否可信做完粒子群优化我会顺手用 Matpower 自带的runopf再解一次同样的问题。这个对照有双重价值。第一runopf的内点法收敛精度高它的网损结果可以作为理论参考值。粒子群如果能跑到和内点法接近的网损说明算法实现没有问题如果差得远大概率是编码范围或罚函数出了问题。第二runopf返回的最优运行点也可以拿来和粒子群最优解做交叉验证。把两组解分别回写到mpc.gen里跑潮流对比它们的约束满足情况。如果粒子群给的解在电压约束、无功约束上全部满足而网损又接近内点法那这个优化结果就可以放心落地方案了。这里要提醒一点runopf是确定性求解粒子群是随机求解两者结果有差异很正常不必追求完全一致。粒子群的工程意义在于处理内点法搞不定的非标准模型而不是替代它。6. 实操中绕不开的坑随机种子、无功越限与潮流不收敛6.1 换一台电脑结果就不一样随机性管理粒子群是随机算法第一次跑出网损 10.5MW第二次跑出 11.8MW这是正常的。但问题是你不希望提交结果时还能复现当时的实验。我刚开始跑的时候吃过这个亏换一个随机种子结果就变评审时完全答不上来自己的数据是怎么来的。解决办法很简单优化循环开始前固定随机流rng(7);固定种子之后同一套参数每次跑的结果完全一致可以复现。我还会把最终结果连同rng种子的状态一并保存save(sprintf(pso_case30_%s.mat, datestr(now, yyyymmdd_HHMM)), gbest, gbestVal, trace, loss0);另外不要只跑一次就信数据。工程实践里我会用 3 到 5 个不同种子各跑一遍取最优结果。这样既保证可复现性又保留随机算法的全局搜索优势。6.2 发电机无功越限与PV-PQ切换问题这是我在这个项目里踩过最深的一个坑。粒子群把机端电压调得很高时发电机可能需要发出很大的无功但每台发电机都有Qmax限制。Matpower 在潮流计算过程中如果发电机无功达到上限会自动把该节点从 PV 节点切换成 PQ 节点并让无功固定在边界值上。这意味着你设定的机端电压实际上并没有达到优化结果的物理含义已经变了。最直接的判断方法是优化结束后检查res.gen里的无功出力是不是落在Qmin到Qmax区间内。一旦有越限要么罚函数起作用要么说明电压上限设得太宽松应该把电压决策范围从 0.94 到 1.06 收紧为 0.95 到 1.05。我在适应度函数里对无功越限做了 3000 的固定惩罚效果好于连续罚函数。因为无功越限是硬性不合理固定大额惩罚能让种群快速淘汰这种个体。6.3 当粒子让潮流算不出来时怎么办粒子群里总有一些个体特别激进比如让某台发电机出力特别小而某个大负荷节点的功率缺额全部压在一条重载线路上潮流可能直接发散。遇到这种情况代码会抛异常如果不处理整个for循环就断了迭代矩阵还少一行最后程序崩得非常难看。我的方案就是前面代码里展示的try-catch加success判断。返回一个超大惩罚值1e6让这个粒子在竞争中自然淘汰。这个操作可以说是粒子群和潮流计算耦合项目里最重要的工程习惯。没有它你会花大量时间在和粒子群斗智斗勇上而不是真正关注优化结果。6.4 罚函数权重与早熟问题罚函数权重过大粒子会很快放弃探索那些暂时越限但可能引导到好解的区域整个种群迅速收缩到一小块安全区域陷入局部最优。我遇到过一次把电压罚函数权重设为 5000 的情况结果迭代曲线在 20 代就平稳了最终网损反而比基态只降了几个百分点。把权重降到 2000 之后收敛曲线明显更健康迭代到 150 代才稳定网损也降得更低。罚函数调参的经验是先放宽后收紧。第一轮优化尽量让粒子群多探索把罚函数权重设小一点第二轮再把关键约束的权重加上去在可行域范围内精细搜索。这个由松到紧的策略可以显著降低早熟概率。6.5 把优化结果回写到MPC再做闭式校验最后一步我强烈建议做拿优化出的gbest重新构造一份mpc然后打印完整的潮流结果包括每个节点的电压幅值、每台发电机的实际出力和无功、每条支路的负载率。不要只看一个网损数字就结束。这份闭式校验数据才是你写报告、写论文、向别人说明优化效果时真正有价值的材料。mpc_opt mpcbase; genBus mpc_opt.gen(:, 1); slackIdx find(mpc_opt.bus(:, 2) 3); isSlack ismember(genBus, slackIdx); mpc_opt.gen(~isSlack, 2) gbest(1:5); mpc_opt.gen(:, 6) gbest(6:11); res_opt runpf(mpc_opt); loss_opt sum(res_opt.branch(:, 14) res_opt.branch(:, 16)); fprintf(优化后网损: %.4f MW, 降损率: %.2f%%\n, loss_opt, (loss0 - loss_opt) / loss0 * 100);这个习惯让我绕过了很多次看起来优化了实际上潮流根本不收敛的尴尬建议你也养成。最后再分享一点我个人的工作习惯每次调完参数我会顺手把随机种子、参数配置、最终网损和迭代曲线全部存成一个.mat文件文件名里带上日期和参数关键字。隔几周再回头看仍然能准确还原当时做了什么事情。这不是什么高深技巧但对做优化仿真的人来说绝对能节省大量返工时间。
粒子群算法结合Matpower的IEEE30节点系统网损优化实践
前阵子搭了一个优化演示工程帮一个朋友把 IEEE 30 节点系统的网损降下来。算法选的粒子群底层潮流计算交给 Matpower 引擎整套跑在 Matlab 里。借这个机会我把整个项目的建模思路、代码骨架、以及调试过程中踩过的坑完整记录下来。严格来说IEEE 30 节点系统是一个输电网测试算例但在实际工作中很多做配电网算法验证的同学也拿它当作标准网架来用所以这个标题组合非常经典。文章不空谈原理重点讲清楚我怎么建模、怎么把粒子群和 Matpower 串起来、以及最终如何分析和判断结果希望能给正在做电力系统智能优化方向的同学一个直接可参考的模板。1. 这个优化问题到底在算什么有功-无功联合优化的数学模型与工程动因1.1 为什么要把有功和无功放在同一个问题里很多刚开始接触电网优化的同学会有一个困惑有功功率和无功功率不是两回事吗有功管系统频率无功管电压水平为什么非要放在一起优化其实这是从运行机理上就耦合在一起的。电网里任何一条线路、一台变压器它的损耗都可以写成 (P_{loss} I^2 R) 的形式而电流 I 的幅值同时由有功潮流和无功潮流共同决定。也就是说线路上的无功流动同样会带来发热损耗。如果你只调整有功出力不管无功分布那么某些线路上的无功环流会一直存在网损降不下去反过来只调无功有功的不合理分配又会让重载线路继续超负荷。所以真正要降网损必须同时调整有功出力和无功分布这就是有功-无功联合优化的由来。配电网的场景里这个问题更明显。配电网的线路 R/X 比值大无功流动造成的电压跌落和线损占比比输电网高得多。分布式电源接入之后有功和无功的耦合关系又进一步复杂化。所以在配电网规划、降损、电压治理这些课题里有功-无功联合优化基本上是一个绕不开的基础问题。1.2 目标函数与决策变量怎么选这次工程里目标函数取的是系统总网损最小。原因很朴素网损是可观测、可计算的直接经济指标也是最好和基态潮流做对比的结果量。你要是想换成发电成本最小、电压偏差最小方法完全一样只是改动适应度函数里的一个表达式而已。决策变量的选取是建模里最关键的一步。IEEE 30 节点系统一共有 30 个节点、41 条支路6 台发电机。我采用的是最常见的变量组合方式除平衡节点外的 5 台发电机有功出力6 台发电机的机端电压幅值。也就是说一个粒子是一个 11 维向量里面的每个分量被提交给 Matpower 潮流计算最终得到一组潮流结果和一个网损值。为什么用机端电压代替无功出力因为这两个量在潮流计算里是强耦合的你设定了机端电压潮流解出来后自然会得到对应的无功注入。直接拿无功出力当决策变量反而麻烦因为 Matpower 在计算时会有 PV 节点和 PQ 节点之间的自动切换逻辑直接用无功出力很容易撞上这个切换边界。这一点后面踩坑部分我会单独说。1.3 约束条件从哪来潮流等式与设备不等式优化问题不能只有目标函数约束条件才是让它像电网而不是像数学游戏的关键。约束分两类等式约束就是潮流方程本身配电网或者输电网的任何运行状态都必须满足节点功率平衡。这个我不会自己手写而是委托给 Matpower 的runpf函数来完成跑通了就代表等式约束被满足了。不等式约束包括发电机有功上下限、发电机无功上下限、节点电压幅值上下限、线路传输容量上限。这些限制在 Matpower 的case30数据文件里都有初始值。bus表里有每个节点的正常电压上限Vmax和下限Vmingen表里有每台发电机的Pmax、Pmin、Qmax、Qmin。在粒子群算法框架下不等式约束大部分靠罚函数处理小部分靠边界截断处理这个设计展开讲就是第 4 章的内容。让我点一下数学上的关键逻辑目标函数要优化的是网损决策变量是发电机有功和机端电压约束是潮流方程加设备上下限。把这三句话写清楚这个优化问题才算是被正确建模了。很多同学代码写不明白往往不是粒子群没掌握而是决策变量和约束没有从 Matpower 数据结构上想透。2. 为什么用粒子群而不是直接调用内点法2.1 Matpower自带最优潮流的局限Matpower 自带的runopf函数是用内点法求解最优潮流OPF的精度高、速度快在干净的标准算例上效果相当好。那为什么还要自己做粒子群一个很现实的原因是实际工程项目里的优化模型往往和教科书上的标准 OPF 不完全一样。你可能要加一个不光滑的约束比如某个节点的电压合格率分段惩罚可能要把目标函数改成某种启发式指标可能要对某些设备动作次数做惩罚。这些非标准、非光滑、甚至不可导的改动对传统内点法是致命的因为内点法需要构造拉格朗日函数并求梯度或者海森矩阵一旦模型变得不光滑数值上非常容易发散。粒子群这类群智能算法则是另一套逻辑它不要求目标函数可导不要求约束光滑它只把目标函数当成一个黑盒来调用。潮流算得出来就给一个具体的网损值算不出来就给一个巨大的惩罚值。这种黑盒优化的特性让它在科研探索阶段显得特别灵活。2.2 粒子群适合这种问题的真正原因工程上选择粒子群主要有三个原因。第一实现成本低。一个有速度-位置更新公式加适应度计算的粒子群核心代码写得好也就一两百行。相比于内点法复杂的矩阵求导和修正方程粒子群的代码量要小得多调起来也更直观。第二天然支持并行。粒子群每一代要计算多个粒子的适应度而这些计算彼此独立完全可以用parfor并行跑。在我这个例子里每个粒子都要调用一次runpf用并行池可以把整体耗时压缩到原来的三分之一左右。第三全局搜索能力比传统梯度类算法好。电网优化问题本质上是有多个局部最优点的非线性非凸问题内点法非常依赖初值初值不好就容易收敛到局部解或者压根不收敛。粒子群通过种群协作的方式在解空间里撒点搜索虽然不保证找到全局最优但工程实践里找到的解往往是可用的而且稳定性好。2.3 它与遗传算法的简单对比也有人问我为什么不直接上遗传算法两个算法我都实现过结论是在电网潮流作为黑盒这个场景下粒子群通常更省事。遗传算法的选择、交叉、变异三个算子每个都需要调试合理的概率参数否则要么收敛太慢要么种群过早失去多样性。粒子群的参数相对集中核心就是惯性权重和学习因子尤其惯性权重可以按迭代次数线性递减整体调参路径清晰得多。另外遗传算法的二进制编码在边界连续优化问题里会有精度和编码长度之间的矛盾而粒子群直接用实数向量操作不需要解码天然贴合发电机出力用连续实数表示的需求。当然如果问题本身是离散变量主导的比如电容器分组投切、变压器分接头挡位遗传算法的离散处理反而更方便。连续变量主导的有功-无功优化粒子群是性价比最高的选择。3. 环境准备与数据摸底Matpower、case30和配电网改造3.1 十分钟搞定Matpower环境Matpower 是一个纯 Matlab 实现的开源电力系统分析工具包不需要编译下载解压就能用。从官网或者 GitHub 下载最新稳定版本之后把解压目录整个加入 Matlab 路径即可。%% 将Matpower路径加入Matlab环境 addpath(genpath(D:/matpower7.1)); savepath; % 保存路径避免每次启动重新设置 mpver(ver); % 查看并确认Matpower版本运行mpver(ver)能看到默认的case30等标准算例已经内置。我建议每次用 Matpower 之前都显式确认版本号因为不同版本对runpf返回结果的字段略有差异确认版本可以少很多莫名其妙的数组越界报错。3.2 读懂case30数据结构bus、gen、branchMatpower 的核心数据结构是mpc结构体它包含bus、gen、branch三个大矩阵以及baseMVA这个基准容量。case30.m文件本质上就是这些矩阵的文本初始化。bus矩阵的关键列第 1 列节点编号第 2 列节点类型1 是 PQ 节点2 是 PV 节点3 是平衡节点第 3、4 列是有功负荷和无功负荷第 8 列是节点电压幅值初始值第 12、13 列是正常运行时电压上限和下限。gen矩阵的关键列第 1 列是发电机所在节点号第 2 列是有功出力第 3 列是无功出力第 4、5 列是无功上限和下限第 6 列是机端电压幅值设定值第 9、10 列是有功上限和下限。branch矩阵记录线路或者变压器的参数包括首末端节点、电阻、电抗、电纳、长期允许载流量、变压器变比等。第一件事永远是用loadcase(case30)载入数据然后用runpf跑一次基态潮流看看初始网损是多少。基态潮流的结果是整个优化项目的比较基准后面所有优化效果的评估都靠它说话。mpc loadcase(case30); % 载入IEEE30节点数据 res0 runpf(mpc); % 基态潮流 loss0 sum(res0.branch(:, 14) res0.branch(:, 16)); % 基态网损单位MW fprintf(基态网损: %.4f MW\n, loss0);res.branch矩阵的第 14 列是首端有功注入第 16 列是末端有功注入两者相加就是该支路的损耗。Matpower 返回的success字段等于 1 才代表潮流收敛这个字段在后续迭代里非常重要。3.3 再说一遍这是输电网算例怎么改成配电网风格我必须把话说清楚IEEE 30 节点系统本身是输电网络线路阻抗、电压等级、负荷水平都不是典型配电网。但是在算法开发和验证阶段它的拓扑结构、发电机数量、约束类型足够丰富完全可以通过改造来模拟配电网特征。常见的改造思路有三种把部分发电机改成分布式电源模型降低其出力上限并修改无功调节范围将线路参数调大电阻分量提高 R/X 比值让电压问题更突出增加部分节点的无功负荷让无功分布更加紧张。这样改完之后还是在case30的框架下做优化但结果更接近配电网降损的真实需求。如果你已经有某个具体的配电网拓扑数据也可以用自己的潮流数据替换掉case30粒子群的代码框架完全不用动只需要改loadcase那一个地方。4. 粒子群与潮流计算耦合编码、适应度函数、主循环4.1 粒子编码设计一串向量对应一台发电机的运行方式粒子群优化的第一个核心步骤是编码。我的设计是把每台非平衡发电机的有功出力和全部发电机的机端电压幅值拼接成一维向量。以case30为例6 台发电机节点编号分别是 1、2、5、8、11、13其中节点 1 是平衡节点。那么决策向量定义为前 5 个分量节点 2、5、8、11、13 这 5 台发电机的有功出力后 6 个分量节点 1、2、5、8、11、13 这 6 台发电机的机端电压幅值。总维度是 11。这里有一个细节有些代码会把平衡节点的有功也拿来做变量这是不对的。平衡节点的有功在潮流里是补差角色系统总有功不平衡时由它自动平衡把它当决策变量会和潮流方程矛盾。为了让粒子初始化范围合理我直接从mpc.gen里读取Pmax、Pmin、Vmax、Vmin来生成上下界。这样粒子群搜索空间和物理设备允许范围自动对齐省去手动维护边界数据的麻烦。% 读取设备上下限 genBus mpc.gen(:, 1); isSlack zeros(size(genBus)); % 平衡节点由bus矩阵type列为3的节点确定 slackIdx find(mpc.bus(:, 2) 3); isSlack ismember(genBus, slackIdx); % 有功决策量范围: 非平衡发电机 PgMin mpc.gen(~isSlack, 10); PgMax mpc.gen(~isSlack, 9); % 电压决策量范围: 全部发电机 VgMin 0.94; VgMax 1.06;4.2 适应度函数网损加罚函数一个值得反复调试的核心适应度函数是粒子群和电网模型之间的唯一接口。我的设计分三层。第一层把粒子的 11 个分量写回mpc.gen调用runpf计算潮流。如果潮流不收敛直接返回一个很大的惩罚值比如1e6。这一步一定要包在try-catch里因为粒子群探索阶段经常会产生不合理的控制变量组合让潮流计算直接报错退出如果不捕获异常整个优化循环就崩了。第二层计算网损。在潮流收敛的情况下用首末端有功之和得到系统网损。第三层罚函数处理越限。我主要检查两个地方节点电压幅值是否越界以及发电机无功出力是否越限。节点电压在res.bus第 8 列发电机无功在res.gen第 3 列。越限越严重惩罚值越大这样粒子群在迭代过程中会自然避开不合理区域。function loss psoFitness(x, mpcBase) mpc mpcBase; genBus mpc.gen(:, 1); slackIdx find(mpc.bus(:, 2) 3); isSlack ismember(genBus, slackIdx); % 前5维写到非平衡机有功 mpc.gen(~isSlack, 2) x(1:5); % 后6维写到全部机端电压 mpc.gen(:, 6) x(6:11); try res runpf(mpc); if res.success ~ 1 loss 1e6; return; end % 网损 loss sum(res.branch(:, 14) res.branch(:, 16)); % 电压越限惩罚 vmin min(res.bus(:, 8)); vmax max(res.bus(:, 8)); if vmin 0.95 loss loss 2000 * abs(vmin - 0.95); end if vmax 1.05 loss loss 2000 * abs(vmax - 1.05); end % 发电机无功越限惩罚 if min(res.gen(:, 3) - res.gen(:, 4)) 0 || min(res.gen(:, 5) - res.gen(:, 3)) 0 loss loss 3000; end catch loss 1e6; end end这个适应度函数里罚系数 2000、3000、1e6 都是我反复试验后比较稳的组合。罚函数设计是这个项目里最需要耐心调试的部分权重太小约束形同虚设权重太大粒子群会过早向可行域收缩丧失搜索能力。后面我会单独讲怎么判断罚函数是否合适。4.3 速度位移更新与边界处理粒子群的标准更新公式不用担心直接用经典版本v w * v c1 * r1 * (pbest - x) c2 * r2 * (gbest - x) x x v其中w是惯性权重c1是自我认知学习因子c2是社会学习因子r1、r2是 0 到 1 的随机数。我的参数设定如下表参数取值说明惯性权重 w0.9 到 0.4 线性递减前期全局搜索后期局部精细搜索学习因子 c12.0保留自我经验学习因子 c22.0群体信息共享种群规模40网络规模不大40 个粒子足够最大迭代次数200配合收敛曲线判断是否提前结束速度上限 vmax0.3每个决策量按范围归一化后限制步长线性递减的惯性权重是工程实践里非常稳定的配置。迭代初期w大粒子跑得快可以探索大范围迭代后期w小粒子在局部精细搜索。很多失败的粒子群应用案例问题都出在全程用同一个w导致搜索行为要么太粗放要么太保守。边界处理上我的做法是粒子越界时直接把对应分量拉回到边界值并且把该维速度设为 0。这个重置边界策略比随机反弹策略更容易让粒子在可行域边界附近稳定下来因为在电力系统优化里很多最优解往往就在发电机出力上限附近。4.4 主循环骨架与计算量预估主循环的骨架大致是这样% 初始化粒子群 nPop 40; nVar 11; maxIter 200; x repmat(xMin, nPop, 1) rand(nPop, nVar) .* repmat(xMax - xMin, nPop, 1); v zeros(nPop, nVar); pbest x; pbestVal arrayfun((i) psoFitness(x(i, :), mpc), 1:nPop); gbestIdx find(pbestVal min(pbestVal), 1); gbest pbest(gbestIdx, :); gbestVal pbestVal(gbestIdx); trace zeros(maxIter, 1); % 主迭代 for iter 1:maxIter w 0.9 - 0.5 * (iter / maxIter); % 线性递减 for i 1:nPop r1 rand(nVar, 1); r2 rand(nVar, 1); v(i, :) w * v(i, :) 2.0 * r1 .* (pbest(i, :) - x(i, :)) 2.0 * r2 .* (gbest - x(i, :)); v(i, v(i, :) vmax) vmax; v(i, v(i, :) -vmax) -vmax; x(i, :) x(i, :) v(i, :); x(i, x(i, :) xMax) xMax; x(i, x(i, :) xMin) xMin; fval psoFitness(x(i, :), mpc); if fval pbestVal(i) pbest(i, :) x(i, :); pbestVal(i) fval; end end [minVal, idx] min(pbestVal); if minVal gbestVal gbestVal minVal; gbest pbest(idx, :); end trace(iter) gbestVal; end计算量方面40 个粒子乘以 200 代就是 8000 次潮流计算。在普通笔记本上一次case30潮流大约几十毫秒所以整个优化过程大概一到三分钟。如果需要跑大规模配电网建议改用parfor对粒子循环做并行化实测能缩短到原来的四分之一左右。5. 跑通之后怎么分析结果迭代曲线、电压分布、网损对比5.1 从基态潮流开始先记住优化前的网损不管优化结果多漂亮没有基准值就毫无意义。我的习惯是把基态潮流结果、基态网损、基态最低最高电压全部打印出来作为后续对比的 Baseline。基态网损由runpf直接得到然后优化结束之后把gbest重新写入mpc.gen再跑一次潮流同样统计网损和电压。两者做差除基态就是降损率。这里有一个细节优化停止时gbest对应的潮流并不一定已经用最终粒子再算过一遍因为pbestVal可能是粒子在某一代的位置计算出来的。所以最终结果必须先重算潮流再统计指标否则容易出现结果和实际不一致的尴尬。5.2 迭代曲线与收敛判断粒子群迭代过程里我把每一代的最优适应度记录在trace数组里用semilogy画收敛曲线。注意这里我通常用对数坐标因为早期适应度可能从几千快速下降到十几个 MW线性坐标会让后期曲线贴在横轴上看不出细节。一个正常的收敛曲线应该是前期快速下降中期平缓波动后期基本稳定。如果曲线后期还在明显持续下降说明迭代次数不够要加大maxIter。如果曲线很快就完全不动了说明种群可能早熟所有粒子已经挤到了同一个局部最优点。figure; plot(1:maxIter, trace, LineWidth, 1.5); xlabel(迭代次数); ylabel(最优网损/MW); title(粒子群收敛曲线); grid on;5.3 用runopf兜底验证粒子群结果是否可信做完粒子群优化我会顺手用 Matpower 自带的runopf再解一次同样的问题。这个对照有双重价值。第一runopf的内点法收敛精度高它的网损结果可以作为理论参考值。粒子群如果能跑到和内点法接近的网损说明算法实现没有问题如果差得远大概率是编码范围或罚函数出了问题。第二runopf返回的最优运行点也可以拿来和粒子群最优解做交叉验证。把两组解分别回写到mpc.gen里跑潮流对比它们的约束满足情况。如果粒子群给的解在电压约束、无功约束上全部满足而网损又接近内点法那这个优化结果就可以放心落地方案了。这里要提醒一点runopf是确定性求解粒子群是随机求解两者结果有差异很正常不必追求完全一致。粒子群的工程意义在于处理内点法搞不定的非标准模型而不是替代它。6. 实操中绕不开的坑随机种子、无功越限与潮流不收敛6.1 换一台电脑结果就不一样随机性管理粒子群是随机算法第一次跑出网损 10.5MW第二次跑出 11.8MW这是正常的。但问题是你不希望提交结果时还能复现当时的实验。我刚开始跑的时候吃过这个亏换一个随机种子结果就变评审时完全答不上来自己的数据是怎么来的。解决办法很简单优化循环开始前固定随机流rng(7);固定种子之后同一套参数每次跑的结果完全一致可以复现。我还会把最终结果连同rng种子的状态一并保存save(sprintf(pso_case30_%s.mat, datestr(now, yyyymmdd_HHMM)), gbest, gbestVal, trace, loss0);另外不要只跑一次就信数据。工程实践里我会用 3 到 5 个不同种子各跑一遍取最优结果。这样既保证可复现性又保留随机算法的全局搜索优势。6.2 发电机无功越限与PV-PQ切换问题这是我在这个项目里踩过最深的一个坑。粒子群把机端电压调得很高时发电机可能需要发出很大的无功但每台发电机都有Qmax限制。Matpower 在潮流计算过程中如果发电机无功达到上限会自动把该节点从 PV 节点切换成 PQ 节点并让无功固定在边界值上。这意味着你设定的机端电压实际上并没有达到优化结果的物理含义已经变了。最直接的判断方法是优化结束后检查res.gen里的无功出力是不是落在Qmin到Qmax区间内。一旦有越限要么罚函数起作用要么说明电压上限设得太宽松应该把电压决策范围从 0.94 到 1.06 收紧为 0.95 到 1.05。我在适应度函数里对无功越限做了 3000 的固定惩罚效果好于连续罚函数。因为无功越限是硬性不合理固定大额惩罚能让种群快速淘汰这种个体。6.3 当粒子让潮流算不出来时怎么办粒子群里总有一些个体特别激进比如让某台发电机出力特别小而某个大负荷节点的功率缺额全部压在一条重载线路上潮流可能直接发散。遇到这种情况代码会抛异常如果不处理整个for循环就断了迭代矩阵还少一行最后程序崩得非常难看。我的方案就是前面代码里展示的try-catch加success判断。返回一个超大惩罚值1e6让这个粒子在竞争中自然淘汰。这个操作可以说是粒子群和潮流计算耦合项目里最重要的工程习惯。没有它你会花大量时间在和粒子群斗智斗勇上而不是真正关注优化结果。6.4 罚函数权重与早熟问题罚函数权重过大粒子会很快放弃探索那些暂时越限但可能引导到好解的区域整个种群迅速收缩到一小块安全区域陷入局部最优。我遇到过一次把电压罚函数权重设为 5000 的情况结果迭代曲线在 20 代就平稳了最终网损反而比基态只降了几个百分点。把权重降到 2000 之后收敛曲线明显更健康迭代到 150 代才稳定网损也降得更低。罚函数调参的经验是先放宽后收紧。第一轮优化尽量让粒子群多探索把罚函数权重设小一点第二轮再把关键约束的权重加上去在可行域范围内精细搜索。这个由松到紧的策略可以显著降低早熟概率。6.5 把优化结果回写到MPC再做闭式校验最后一步我强烈建议做拿优化出的gbest重新构造一份mpc然后打印完整的潮流结果包括每个节点的电压幅值、每台发电机的实际出力和无功、每条支路的负载率。不要只看一个网损数字就结束。这份闭式校验数据才是你写报告、写论文、向别人说明优化效果时真正有价值的材料。mpc_opt mpcbase; genBus mpc_opt.gen(:, 1); slackIdx find(mpc_opt.bus(:, 2) 3); isSlack ismember(genBus, slackIdx); mpc_opt.gen(~isSlack, 2) gbest(1:5); mpc_opt.gen(:, 6) gbest(6:11); res_opt runpf(mpc_opt); loss_opt sum(res_opt.branch(:, 14) res_opt.branch(:, 16)); fprintf(优化后网损: %.4f MW, 降损率: %.2f%%\n, loss_opt, (loss0 - loss_opt) / loss0 * 100);这个习惯让我绕过了很多次看起来优化了实际上潮流根本不收敛的尴尬建议你也养成。最后再分享一点我个人的工作习惯每次调完参数我会顺手把随机种子、参数配置、最终网损和迭代曲线全部存成一个.mat文件文件名里带上日期和参数关键字。隔几周再回头看仍然能准确还原当时做了什么事情。这不是什么高深技巧但对做优化仿真的人来说绝对能节省大量返工时间。